Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration.
An important variation of preconditioned conjugate gradient algorithms is inexact preconditioner implemented with inner-outer iterations [G. H. Golub and M. L. Overton,  Numerical Analysis, Lecture Notes in Math. 912, Springer, Berlin, New York, 1982], where the preconditioner is solved by an inner iteration to a prescribed precision. In this paper, we formulate an inexact preconditioned conjugate gradient algorithm for a symmetric positive definite system and analyze its convergence property. We establish a linear convergence result using a local relation of residual norms. We also analyze the algorithm using a global equation and show that the algorithm may have the superlinear convergence property when the inner iteration is solved to high accuracy. The analysis is in agreement with observed numerical behavior of the algorithm. In particular, it suggests a heuristic choice of the stopping threshold for the inner iteration. Numerical examples are given to show the effectiveness of this choice and to compare the convergence bound.
Introduction
Iterative methods for solving linear systems are usually combined with a preconditioner that can be
easily solved. For some practical problems, however, a natural and efficient choice of preconditioner
may be one that can not be solved easily by a direct method and thus may require an iterative
method (called inner iteration) itself to solve the preconditioned equations. There also exist cases
where the matrix operator contains inverses of some other matrices, an explicit form of which is
not available. Then the matrix-vector products can only be obtained approximately through an
inner iteration. The linear systems arising from the saddle point problems [3] is one such example.
For these types of problems, the original iterative method will be called the outer iteration and
the iterative method used for solving the preconditioner or forming the matrix-vector products is
called the inner iteration.
A critical question in the use of the inner-outer iterations is to what precision the preconditioner
should be solved, i.e., what stopping threshold should be used in the inner iteration. Clearly, a
very high precision will render the outer iteration close to the exact case and a very low one on the
other hand could make the outer iteration irrelevant. An optimal one will be to allow the stopping
Address Computing and Computational Mathematics Program, Department of Computer Science,
Stanford University, Stanford, CA 94305, USA E-mail : golub@sccm.stanford.edu. Research supported in part by
National Science Foundation Grant DMS-9403899.
Department of Applied Mathematics, University of Manitoba, Winnipeg, Manitoba, Canada R3T 2N2.
E-mail: ye@gauss.amath.umanitoba.ca. Research supported by Natural Sciences and Engineering Research Council
of Canada.
threshold as large as possible in order to reduce the cost of inner iteration, while maintaining the
convergence characteristic of the outer iteration. In other words, we wish to combine the inner
and outer iterations so that the total number of operations is minimized. An answer to the above
question requires understanding of how the accuracy in the inner iteration affects the convergence of
the outer iteration. This has been studied by Golub and Overton [5, 6] for the Chebyshev iteration
and the Richardson iteration, by Munthe-Kaas [8] for preconditioned steepest descent algorithms,
by Elman and Golub [3] for the Uzawa algorithm for the saddle point problems, and by Giladi,
Golub and Keller [4] for the Chebyshev iteration with a varying threshold. Golub and Overton also
observed the interesting phenomenon for the preconditioned conjugate gradient algorithm (as the
outer iteration) that the convergence of CG could be maintained for very large stopping threshold
in the inner iteration; yet the convergence rate may be extremely sensitive to the change of the
threshold at certain point.
Given that the known convergence properties of CG depends strongly on a global minimization
property, the phenomenon found in [5, 6] seems very surprising and makes the inexact preconditioned
conjugate gradient an attractive option for implementing preconditioners. However, there
has been no theoretical analysis to explain these interesting phenomena, and its extreme sensitivity
to the threshold makes it hard to implement it in practice. The present paper is an effort in this
direction. We shall formulate and analyze an inexact preconditioned conjugate gradient method for
a symmetric positive definite system. By establishing a local relation between consecutive residual
norms, we prove a linear convergence property with a bound on the rate and illustrate that the
convergence rate is relatively insensitive to the change of threshold up to certain point. In par-
ticular, the result is used to arrive at a heuristic choice of the stopping threshold. We also show,
using a global relation as in [9], that the algorithm may have the superlinear convergence property
when the global orthogonality is nearly preserved, which usually occurs with smaller thresholds
and shorter iterations.
The paper is organized as follows. In section 2, we present the inexact preconditioned conjugate
gradient algorithm, and some of its properties together with two numerical examples illustrating
its numerical behaviour. We then give in section 3 a local analysis, showing the linear convergence
property, and in section 4, a global analysis, showing the superlinear convergence property. Finally,
we give some numerical examples in section 5 to illustrate our results.
We shall use the standard notation in numerical analysis. The M-norm k \Delta k M and the M -
inner product are defined by kvk
respectively. cond(A) denotes the spectral condition number of a matrix A. A + denotes the
Moore-Penrose generalized inverse of A.
2 PCG with inner-outer Iterations
We consider the preconditioned conjugate gradient (PCG) algorithm for solving with a
preconditioner M , where both A and M are symmetric positive definite. Then at each step of PCG
iteration, a search direction p n is found by first solving the preconditioned system Mz . If a
direct method is not available for solving M , an iterative method, possibly (though not necessarily)
CG itself, can be used to solve it. In this case, we find z n by the inner iteration such that
with the stopping criterion
is the stopping threshhold in the inner iteration. Here we have used the M \Gamma1 -norm
for the theoretical convenience. In practical computations, one can replace (RES) by the following
criterion in terms of the 2-norm
Because z n is only used to define the search direction p n in PCG, only its direction has any
significance. Therefore, we also propose the following direction based stopping criterion
i.e., the acute angle between Mz n and r n in the M \Gamma1 -inner product (or the acute angle between
z n and M \Gamma1 r n in the M-inner product) is at most '.
Remark 1: It can be easily proved (using a triangular relation in the M \Gamma1 -inner product) that if
(RES) is satisfied, then
i.e., (ANG) is satisfied with
sin
However, the converse implication is not necessarily true. Therefore the criterion (ANG) is less
restrictive than (RES). Namely, (ANG) may be satisfied before (RES) does in the inner itera-
tion. However, our numerical tests suggest that there is usually little difference between the two.
Our introduction of (ANG) and the use of M \Gamma1 -inner product are primarily for the theoretical
convenience.
Remark 2: If the inner iteration is carried out by CG, the residual e n is orthogonal to Mz n in
the M \Gamma1 -inner product. Then again by a triangular relation, we have
s
So in this case, (RES) and (ANG) are equivalent. It can also be proved that the CG iteration
minimizes the residual as well as the angle in (ANG) over the Krylov subspace concerned.
Now, we formulate the inexact preconditioned conjugate gradient algorithm (IPCG) as follows.
Inexact Preconditioned Conjugate Gradient Algorithm (IPCG):
For
z T
equivalently ff
end for
Clearly, if j or the above algorithm is the usual PCG. With j (or ') ? 0, we allow the
preconditioner Mz to be solved only approximately.
Remark 3. As is well-known, there are several different formulations for ff n and fi n in the algo-
rithm, which are all equivalent in the exact PCG. For the inexact case, this may no longer be true.
Our particular formulation here implies that some local orthogonality properties are maintained
even with inexact z n , as given in the next lemma. Our numerical tests also show that it indeed
leads to a more stable algorithm. For example, if fi n is computed by the old form z T
the algorithm may not converge for larger j.
Lemma 1 The sequences generated by Algorithm IPCG satisfies the following local orthogonality
Proof First
z T
Then
Thus z T
z T
Now supposing p T
and
Thus the lemma follows from p T
by an induction argument.
By eliminating p n in the recurrence, IPCG can be written with two consecutive steps as a second
order recurrence
a unified form that includes Chebyshev and second order Richardson iterations (cf [2]). From the
local orthogonality, we obtain the following local minimization property.
Proposition 1
and
Proof By the orthogonality, kr
A \Gamma1 . It
is easy to check that z T n r n 6= 0 for either (RES) or (ANG). So ff n 6= 0 and thus we have the strict
inequality.
For the second inequality, we just need to show r n+1 ? A \Gamma1 Az n and r This follows from
r T
Recall that the steepest descent method constructs r sd
Ar n from r n such that
kr sd
. A bound on kr sd
can be obtained from the Kantrovich
inequality. An inexact preconditioned version of the steepest descent method [8] is to construct
r sd
n Az n such that
kr sd
A
z
Then the second inequality of the above proposition shows that kr n+1 k A
(by choosing
2.1 A Numerical Example
To motivate our discussion in the next two sections, we now present two numerical examples to
illustrate the typical convergence behaviour of IPCG. We simulate IPCG by artificially perturbing
applying PCG with I and z is a pseudo random
vector with entries uniformly distributed in (\Gamma0:5; 0:5), and j is a perturbation parameter. The following
are the convergence curves for two diagonal matrices
(10000 \Theta 10000) and
1000] (1000 \Theta 1000) with a random constant term.
We first note that the A \Gamma1 -norm of the residuals decreases monotonically (cf Prop. 2.1) and
convergence occurs for quite large j. We further observe from the second example that for smaller
the superlinear convergence property of exact CG is recovered. In fact, they follow the
unperturbed case very closely. When j is larger, the convergence tends to be linear, with a rate
depending on j. Interestingly, the convergence rate is insensitive to the change of magnitude of j
until a certain point (0:1 in the first and 0:4 in the second), around which it becomes extremely
sensitive to a relatively small change of j.
Our explanation and analysis to the above observation will be given in two categories, although
there is no clear boundary between the two. For smaller j, some global properties (e.g. the
orthogonality among r n ) is expected to be preserved, which leads to a global near minimization
property and thus the superlinear convergence. For larger j, the global minimization property will
be lost. However, we will show that some local relation from r n to r n+1 is not destroyed, which turns
out to preserve a linear convergence property. We consider each of these two categories separately
in the next two sections.
3 Linear Convergence of IPCG
In this section we present a local relation between consecutive residual norms, which lead to a linear
convergence bound for IPCG. Our basic idea is to relate the reduction factor of IPCG

Figure

1: Convergence curves for various perturbations j (solid (bottom):
number of iterations
A-inverse-norm
of
residuals
number of iterations
A-inverse-norm
of
residuals
to the reduction factor of the steepest descent method (along the inexact preconditioned gradient
direction z n , see (3) of section 2)
r=rn \GammatAz n
From Proposition 2.1, we have oe n -
Theorem 1 Let oe n and fl n be defined as in (4) and (5) for IPCG. Then for
(a) oe
with
z
(b)
r T
r T
\Gammag n
\Gammag
where g k is defined by
and
Proof From r T
r T
Substituting ff
Apn in, we obtain
r T
Thus
r T
using z
Noting
that
from (2), we have
ff \Gamma2
z T
z T
oe
n r n , and we have used Ap
substituting
the above and (5) into (8), we obtain part (a).
For part (b), let
. From (a),
oe
oe
first step of IPCG amounts to one step of the steepest descend method),
by the definition. Now (b) follows from r T
and oe
In the exact case, j by the orthogonality and the above equations are simplified. In the
inexact case, using the local orthogonality (Lemma 2.1), j n can be bounded. We shall consider
the stopping criterion (ANG) only in this section. Since (RES) implies (ANG) with sin
results here apply to the (RES) case as well by simply replacing sin ' by j.
First we need the following lemma that is geometrically clear. (An acute angle ' 2 [0; -=2]
between two vectors u; v in an inner product ! \Delta; \Delta ? is defined by cos
resp.) be the acute angle between u and u
then the acute angle between u 1 and v 1 is at least -=2
Proof We can assume that u; u are all unit vectors. Then we write
vector orthogonal to u ( v resp.
Let
where we note that u is orthogonal to v and v orthogonal to u.
is the norm associated with
- a cos ' 1 sin
- a cos ' 1 sin
where we note that the second last expression is an increasing function of a 2 [0; 1] for
and therefore is bounded by its value at a = 1.
Lemma 3 If z n satisfies (ANG) with ' -=4, then
z T
z
Proof First (Mz
By applying the above Lemma with the inner
product defined by M \Gamma1 to the pairs Mz which satisfy (ANG), we obtain
i.e. jz T
On the other hand, jz T
Furthermore, it is easy to check that, for any vector v
A
where - min and - max denotes respectively the minimal and the maximal eigenvalues. Thus,
which completes the proof of the lemma.
We now present a bound for fl n that has been derived by Munthe-Kaas [8]. We shall use one of
the variations of the bounds in [8, Theorem 5] and we repeat some arguments of [8] below. First
the following lemma is a generalization of [1, Corollary 4]
Lemma 4 [8, Lemma 2] Suppose p and q 2 R n are such that
kpk kqk
(p
(q
where W is a symmetric positive definite matrix and -
kpk kqk
So applying the above Lemma to p; q and W , we obtain
r T
(p
(q
1\Gammasin ' . Then we have
r T
We now present a linear convergence bound for IPCG.
Theorem 2 If
converges and for even n,
where
s
Proof We consider two consecutive oe n of IPCG. By
apply. So from Theorem 1 we have
oe
where fl is as defined in Lemma 5. Then
\Gammag n+1
and
. So
So, oe n oe n+1 - (oeK) 4 . Therefore, for even n, the bound follows from kr n+1 k A
shows that IPCG converges. Finally, by expanding oeK in terms
of
using
In terms of the stopping criterion (RES), the same result holds with sin ' replaced by j (see
Remark 1). Therefore, if
converge with a rate depending
on oe (the standard PCG convergence rate). In particular, the bound indicates that the IPCG
convergence rate is relatively insensitive to the magnitude of j for smaller j but increases sharply
at certain point (see the rate curves in Fig. 2). However, the bound on the convergence rate here
tends to be pessimistic and it does not recover the classical bound for the case
it does seem to reflect the trend of how the rate changes as j changes.
To compare the bound with the actual numerical results, we consider an example similar to the
one in section 2. Namely, we consider a diagonal matrix whose eigenvalues are linearly distributed
on [1; -], and apply PCG with the same kind of random perturbation. We carry out IPCG and
compute the actual convergence rate by (kr ranging from 10 \Gamma6 to 1. In

Figure

2 the graphs of the bound and the actual computed rate are plotted (for the case
and 100).
By comparing the actual convergence rate and its bound, we observe that the bound, as a worst
case bound, follows the trend of the actual convergence rate curve quite closely. The bound reaches
1 at
In particular, j 0 seems to be a good estimate of the point at which the actual rate starts to increase
significantly (i.e., the slope is greater than 1). Note that when the slope is less than 1, any increase

Figure

2: Actual Convergence Rate and its Bound verses j
actual rate
bound
00.50.70.9stopping threshold eta=sin(theta)
convergence
rate
actual rate
bound
00.750.850.95stopping threshold eta=sin(theta)
convergence
rate
in the rate may be compensated by a comparable or more increase in j. We therefore advocate a
value around j 0 as a heuristic choice of the stopping threshold for inner iterations. Our numerical
examples in Section 5 confirm that this is indeed a reasonable strategy in balancing the numbers
of inner and outer iterations.
4 Superlinear Convergence of IPCG
The bound in the previous section demonstrates linear convergence of IPCG. We observed in section
2 that for smaller j, IPCG may actually enjoy the superlinear convergence property of exact CG.
Here we explain this phenomenon by the method of [9], i.e. by considering a global equation that
is approximately satisfied by IPCG. We remark that a global property is necessary in examining
superlinear convergence.
Let
From r respectively, we obtain
the following matrix equations for IPCG
r
where
and U n =B
Combining the equations in (11), we obtain the following equation
AM
r
Note that -
U n is a tridiagonal matrix such that e T
Therefore, the inexact case satisfies an equation similar to the exact
case with the error term AM
We rewrite (13) in a scaled form as in [9, Eq. (8)]. Let D
Then from (13)
r n+1
where
By the
stopping criteria,
Now applying the
same argument in the proof of [9, Theorem 3.5] to (14), we obtain the following theorem. (The
details are omitted here).
Theorem 3 Assume r are linearly independent and let V T
(i.e. the matrix consisting of the first n rows of -
where
r T
To interpret the above result, we note that kr T
R n k=kr n+1 k A \Gamma1 is a measure of M
orthogonality among the residual vectors. If then the residuals of PCG are orthogonal
with respect to M \Gamma1 and hence K but is not too large, the loss of
orthogonality among the residuals may be gradual and it will take modest length of run before
n are of magnitude j. Therefore, in this regime, kr n+1 k A
. Note that ffl
decreases superlinearly because of annihilation of the extreme spectrum (see [10]).
See [9] for some artificially perturbed numerical examples.
In summary, if the global M \Gamma1 -orthogonality among the residual vectors are nearly maintained
to certain step, the residual of IPCG is very close to that of exact PCG up to that point and thus
may display the superlinear convergence property.
5 Numerical Examples
In this section, we present numerical examples of inner-outer iterations, testing various choices of
the stopping threshold as compared with
- of (10), where A). For this
purpose, we shall consider
- for d ranging 0:01 to 5 as well as
We consider
with the homogeneous Dirichlet boundary condition Using uniform five-point finite
difference with the step size
N+1 , we obtain n \Theta n (with
A is the discretization of \Gammar(a(x; y)r) and is a (N \Theta N) block tridiagonal matrix.
We consider solving this equation using the block Jacobi preconditioner M (i.e. M is the block
diagonal part of A) or using the discrete Laplacian L (i.e. L is the discretization of \Gamma\Delta) as the
preconditioner. For the purpose of testing inner-outer iterations, an iterative method (i.e. SOR,
CG or preconditioned CG) is used for the preconditioners. We denote them by CG-SOR, CG-CG
and CG-PCG respectively.
We compare the number of outer (n outer ) and total inner (n inner ) iterations required to reduce
the A \Gamma1 -norm of the residual by 10 \Gamma8 . will be used in our tests.
In the first test, a(x; and the discrete Laplacian
preconditioner L is used. SOR (with the optimal parameter [11]), CG and PCG with the modified
incomplete Cholesky factorization are used in the inner iteration to solve r. The results are
listed in Tables 1, 2 and 3.
used.

Table

1: Iteration Counts for CG-SOR with discrete Laplacian preconditioner L
In the first rows of the tables, we also list the results for 1, in which case one step of inner
iteration is carried out for each outer iteration. In this way, it is closely related to those cases with
very close to 1. Interestingly, however, this extreme case is usually equivalent to applying CG
directly to the original matrix A or with a preconditioner. For example, if the inner iteration is
CG itself, one step of inner CG produces inner solution z n in the same direction of r n and thus the
outer iteration is exactly CG applied to A (see Appendix for a detailed discussion for other inner
solvers). This is why that convergence still occurs in these extreme cases and our analysis, which are
based on Mz only, would not include this case. However, if z n is chosen only to satisfy
1.0 136 136 1.0 232 232 1.0 478 478

Table

2: Iteration Counts for CG-CG with discrete Laplacian preconditioner L

Table

3: Iteration Counts for CG-PCG with discrete Laplacian preconditioner L
but not in the particular way here, the convergence is not expected. The relatively
small iteration counts are due to the fact that the original system is not too ill-conditioned.
In comparing the performance of different j with respect to the outer iteration counts, it appears
that lies right around the point at which the outer iteration count starts to increase significantly.
This confirms our convergence analysis for the outer iteration. For the total inner iteration counts,
the performance for larger j seems to be irregular among different inner solvers, and this can
be attributed to the different convergence characteristic at the extreme case different
inner solvers (see Appendix). In particular, we observe that for larger j, CG-CG (or CG-PCG)
performs better than CG-SOR. This phenomenon was also observed in [3]. Overall, j 0 seems to be
a reasonable choice in balancing the numbers of inner and outer iterations.
In the above example, - and thus j 0 remains nearly constant for different N . In the second
test, we use the block Jacobi preconditioner M and a(x;
y. Then
and 100 respectively. Both SOR and CG are used in the inner iteration. The results are listed in

Tables

4 and 5 in an Appendix. Similar behaviour was observed.
6 Conclusion
We have formulated and analyzed an inexact preconditioned conjugate gradient method. The
method is proved to be convergent for fairly large thresholds in the inner iterations. A linear
convergence bound, though pessimistic, is obtained, which leads to a heuristic choice of the stopping
threshold in the inner iteration. Numerical tests demonstrate the efficiency of the choice.
It still remains an unsolved problem to choose an optimal j that minimizes the total amount of
work (see [4]), although j 0 here provides a first approximation. Solving such a problem demands a
sharper bound in the outer iteration and analysis of the near extreme threshold cases. It is not clear
whether a better bound could be obtained from the approach of the present paper. It seems there
are more properties of IPCG awaiting for discovery. For example, better bounds for the steepest
descent reduction factor fl n may exist for IPCG, which in turn would lead to improvement to the
results here.



--R

Some inequalities involving the euclidean condition of a matrix
A Generalized Conjugate Gradient Method for the Numerical Solution of Elliptic Partial Differential Equations
Inexact and preconditioned Uzawa algorithms for saddle point problems
Inner and outer iterations for the Chebyshev algorithm Stanford SCCM Technical Report 95-12
Convergence of a two-stage Richardson iterative procedure for solving systems of linear equations
The convergence of inexact Chebyshev and Richardson iterative methods for solving linear systems.
Matrix Computations
The convergence rate of inexact preconditioned steepest descent algorithm for solving linear systems
Analysis of the Finite Precision Bi-Conjugate Gradient algorithm for Nonsymmetric Linear Systems
The rate of convergence of conjugate gradients
Iterative Solution of Large Linear Systems Academic Press
--TR

--CTR
Carsten Burstedde , Angela Kunoth, Fast iterative solution of elliptic control problems in wavelet discretization, Journal of Computational and Applied Mathematics, v.196 n.1, p.299-319, 1 November 2006
Angela Kunoth, Fast Iterative Solution of Saddle Point Problems in Optimal Control Based on Wavelets, Computational Optimization and Applications, v.22 n.2, p.225-259, July 2002
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
conjugate gradient method;inexact preconditioner;inner-outer iterations