Download Numerical Linear Algebra, Lecture Notes - Mathematics - 7 and more Study notes Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
12 GMRES
We are now going to discuss how to solve the linear system Ax = b for arbitrary matrices
A ∈ R
n×n using similar ideas as we did in the CG method.
One possibility would be to look at the normal equations A
T Ax = A
T b, which lead to the
same solution, provided the matrix A is invertible. Since A
T A is symmetric and positive
definite, we can use the CG method to solve this problem. If x
∗ is the solution and
x j
∈ x 1
+ K
j− 1
(A
T r 1
, A
T A) is a typical iteration then we have the identity
‖x
∗ − x j
2
A
T A
= (x
∗ − x j
T A
T A(x
∗ − x j
) = (Ax
∗ − Ax j
T (Ax
∗ − Ax j
= ‖b − Ax j
2
2
This explains why this method is called CGNR-method, which stands for “Conjugate
Gradient on the Normal equation to minimize the Residual”.
However, the method is only of limited use since we know that we have to expect a condition
number which is the square of the original condition number.
Nonetheless, we can generalise this idea in the following way.
Definition 12.1 Assume A ∈ R
n×n is invertible and b ∈ R
n
. Let x 1
∈ R
n and r 1
b − Ax 1. The jth iteration of the GMRES method (Generalized Minimum Residual) is
defined to be the solution of
min{‖b − Ax‖ 2 : x ∈ x 1 + Kj− 1 (r 1 , A)}.
The analysis of the GMRES method is similar to the analysis of the CG method. Since
every x ∈ x 1
+ K
j− 1
(r 1
, A) can be written as x j
= x 1
j− 1
k=
a k
A
k− 1 r 1
, we can conclude
that
b − Ax j
= b − Ax 1
j− 1 ∑
k=
a k
A
k r 1
=: Q(A)r 1
with a polynomial Q ∈ π j− 1
satisfying Q(0) = 1. Hence, minimising over all those x is
equivalent to minimising over all these Q.
Theorem 12.2 If A ∈ R
n×n is invertible and if {x j
} denote the iterations of the GMRES
method then
‖b − Axj ‖ 2 = min
Q∈π j− 1
Q(0)=
‖Q(A)r 1 ‖ 2 ≤ ‖P (A)r 1 ‖ 2 ≤ ‖P (A)‖ 2 ‖r 1 ‖ 2
for every P ∈ πj− 1 with P (0) = 1. Moreover, the method needs at most n steps to terminate
with the solution of Ax = b.
Proof: It remains to show that the method stops with the solution x
∗
. To see this, let
ϕ(t) = det(A − tI) be the characteristic polynomial of A. Then, ϕ ∈ π n
and, since A is
invertible, ϕ(0) = det A %= 0. Hence, we can choose the polynomial P := ϕ/ϕ(0) in the
last step. However, the theorem of Cayley-Hamilton gives ϕ(A) = 0 and hence P (A) = 0
so that must have b = Ax n+
, unless the method terminates before.!
In the case of the CG method, we then used the fact that there is a basis of R
n consisting
of orthonormal eigenvectors of A. This will not be true in the general case that we are
dealing with now.
Nonetheless, we will assume that the matrix A is diagonalisable, which means that there
is a matrix V ∈ C
n×n and a diagonal matrix D = diag(λ 1
,... , λ n
) ∈ C
n×n , such that
A = V DV
− 1
.
The diagonal elements of D are the eigenvalues of A. From now on, we have to work with
complex matrices.
If P ∈ πj− 1 is of the form P (t) = 1 +
j− 1
k=
akt
k then we can conclude that
P (A) = P (V DV
− 1 ) = 1 +
j− 1 ∑
k=
a k
(V DV
− 1 )
k = 1 +
j− 1 ∑
k=
a k
V D
k V
− 1
= V P (D)V
− 1 ,
and we derive the following error estimate, which is, if compared to the corresponding
result of the CG method, worse by a factor κ 2 (V ).
Corollary 12.3 Assume the matrix A ∈ R
n×n is similar to a diagonal matrix with eigen-
values λ j
∈ C. The jth iteration of the GMRES method satisfies the error bound
‖b − Axj ‖ 2 ≤ κ 2 (V ) min
P ∈π j− 1
P (0)=
max
1 ≤k≤n
|P (λk)|.
Note that, in contrast to the corresponding minimisation problem of the CG method, here,
the eigenvalues λ j
can be complex, which complicates matters. We will not further pursue
this here, instead, we will have a look at a possible implementation.
The implementation is based upon the idea to find an orthonormal basis for Kj (r 1 , A) and
then to solve the minimisation problem within this basis.
The specific structure of K j
(r 1
, A) allows us to use the orthogonalisation method of Gram-
Schmidt. To this end, we start with setting r 1 = b − Ax 1 and v 1 = r 1 /‖r 1 ‖ 2. Then, for
1 ≤ k ≤ j − 1 we compute
- wk = Avk −
k
!=
(Avk, v!) 2 v!,
- vk+1 = wk/‖wk‖ 2.
This method is also called Arnoldi method, and, for the time being, we assume that it is
well defined in every step, which means that we assume that wk %= 0 always.
Lemma 12.4 If the Arnoldi method is well defined up to step j then {v 1
,... , v j
} is an
orthonormal basis of K j
(r 1
, A).
Algorithm 6: GM RES
Input : A ∈ R
n×n invertible, b ∈ R
n , x 1
∈ R
n
r 1
:= b − Ax 1
, β = ‖r 1
2
, v 1
= r 1
/β
for j = 1, 2 ,... do
wj := Avj
for k = 1 to j do
h kj
= (w j
, v k
2
wj = wj − hkj vk
hj+1,j = ‖wj ‖ 2
if h j+1,j
= 0 then break
else v j+
= w j
/h j+1,j
Compute yj as the solution of ‖βe 1 − Hj y‖ 2
xj+1 = x 1 + Vj yj
Output : Solution x
∗ of Ax = b.
