Page 177 - Matrices theory and applications

9. Iterative Methods for Linear Problems
160
(at most n). However, the round-off errors pollute the final result, and we
would prefer to consider the conjugate gradient as an iterative method in
which the number N of iterations, much less than n, gives a rather good
approximation of ¯x. We shall see that the choice of N is linked to the
condition number of the matrix A.
We denote by ·, ·
the canonical scalar product on IR .When A ∈ SPD n
n
and b ∈ IR , the function n
1
x → J(x):= Ax, x
− b, x
2
is strictly convex and tends to infinity as x → +∞. It thus reaches
its infimum at a unique point ¯x, which is the unique vector where the
gradient of J vanishes. We shall denote by r (for residue) the gradient of
J: r(x)= Ax − b. Hence ¯x is the solution of the linear system Ax = b.
n
If A¯x = b and x ∈ IR , x =¯x,then
1
J(x)= J(¯x)+ A(x − ¯x),x − ¯x
>J(¯x). (9.4)
2
The conjugate gradient is thus a descent method.
We shall denote by E the quadratic form associated to A: E(x):=
n
Ax, x
. It is the square of a norm of IR . The character ⊥ A indicates
the orthogonality with respect to the scalar product defined by A.
9.5.1 A Theoretical Analysis
n
Let x 0 ∈ IR be given. We define e 0 = x 0 − ¯x, r 0 = r(x 0 )= Ae 0 .Wemay
assume that e 0 = 0; otherwise, we would already have the solution. For
k ≥ 1, let us define the vector space
H k := {P(A)r 0 | P ∈ IR[X], deg P ≤ k − 1}, H 0 = {0}.
In H k+1 , the linear subspace H k is of codimension 0 or 1. In the first case,
H k+1 = H k , and it follows that H k+2 = AH k+1 + H k+1 = AH k + H k =
H k+1 = H k and thus by induction, H k = H m for every m>k.Let us
denote by l the smallest index such that H l = H l+1 .For k< l, H k is thus
of codimension one in H k+1 , while if k ≥ l,then H k = H k+1 . It follows
that dim H k = k if k ≤ l.In particular, l ≤ n.
One can always find, by Gram–Schmidt orthonormalization, an A-
1
orthogonal basis (that is, such that Ap j ,p i
=0 if i = j) {p 0 ,... ,p l−1 }
of H l such that {p 0 ,... ,p k−1 } is a basis of H k when k ≤ l. The vectors p j ,
which are not necessarily unit vectors, are defined, up to a scalar multiple,
by
p k ∈H k+1 , p k ⊥ A H k .

1
One must distinguish in this section between the two scalar products, namely ·, ·
and A·, · .