Lecture AIMS Notes 2006
Peter J. Olver
7. Iterative Methods for Linear Systems
Linear iteration coincides with multiplication by successive powers of a matrix; convergence of the iterates depends on the magnitude of its eigenvalues. We discuss in some detail a variety of convergence criteria based on the spectral radius, on matrix norms, and on eigenvalue estimates provided by the Gerschgorin Circle Theorem. We will then turn our attention to the three most important iterative schemes used to accurately approximate the solutions to linear algebraic systems. The classical Jacobi method is the simplest, while an evident serialization leads to the popular GaussSeidel method. Completely general convergence criteria are hard to formulate, although convergence is assured for the important class of diagonally dominant matrices that arise in many applications. A simple modification of the GaussSeidel scheme, known as Successive Over-Relaxation (SOR), can dramatically speed up the convergence rate, and is the method of choice in many modern applications. Finally, we introduce the method of conjugate gradients, a powerful "semi-direct" iterative scheme that, in contrast to the classical iterative schemes, is guaranteed to eventually produce the exact solution.
7.1. Linear Iterative Systems.
We begin with the basic definition of an iterative system of linear equations. Definition 7.1. A linear iterative system takes the form u(k+1) = T u(k) , u(0) = a. (7.1) The coefficient matrix T has size n n. We will consider both real and complex systems, and so the iterates u(k) are vectors either in R n (which assumes that the coefficient matrix T is also real) or in C n . For k = 1, 2, 3, . . ., the solution u(k) is uniquely determined by the initial conditions u(0) = a. Powers of Matrices The solution to the general linear iterative system (7.1) is, at least at first glance, immediate. Clearly, u(1) = T u(0) = T a,
u(2) = T u(1) = T 2 a,
u(3) = T u(2) = T 3 a,
Warning: The superscripts on u(k) refer to the iterate number, and should not be mistaken for derivatives.
10/18/06
103
c 2006
Peter J. Olver
and, in general, u(k) = T k a. (7.2) Thus, the iterates are simply determined by multiplying the initial vector a by the successive powers of the coefficient matrix T . And so, unlike differential equations, proving the existence and uniqueness of solutions to an iterative system is completely trivial. However, unlike real or complex scalars, the general formulae and qualitative behavior of the powers of a square matrix are not nearly so immediately apparent. (Before continuing, the reader is urged to experiment with simple 2 2 matrices, trying to detect patterns.) To make progress, recall how we managed to solve linear systems of differential equations by suitably adapting the known exponential solution from the scalar version. In the iterative case, the scalar solution formula (2.8) is written in terms of powers, not exponentials. This motivates us to try the power ansatz u(k) = k v, (7.3)
in which is a scalar and v is a fixed vector, as a possible solution to the system. We find u(k+1) = k+1 v, while T u(k) = T (k v) = k T v.
These two expressions will be equal if and only if T v = v. Therefore, (7.3) is a nontrivial solution to (7.1) if and only if is an eigenvalue of the coefficient matrix T and v = 0 an associated eigenvector . Thus, to each eigenvector and eigenvalue of the coefficient matrix, we can construct a solution to the iterative system. We can then appeal to linear superposition to combine the basic power solutions to form more general solutions. In particular, if the coefficient matrix is complete, then this method will, as in the case of linear ordinary differential equations, produce the general solution. Theorem 7.2. If the coefficient matrix T is complete, then the general solution to the linear iterative system u(k+1) = T u(k) is given by u(k) = c1 k v1 + c2 k v2 + + cn k vn , 1 2 n (7.4)
where v1 , . . . , vn are the linearly independent eigenvectors and 1 , . . . , n the corresponding eigenvalues of T . The coefficients c1 , . . . , cn are arbitrary scalars and are uniquely prescribed by the initial conditions u(0) = a. Proof : Since we already know that (7.4) is a solution to the system for arbitrary c1 , . . . , cn , it suffices to show that we can match any prescribed initial conditions. To this end, we need to solve the linear system u(0) = c1 v1 + + cn vn = a. (7.5)
Completeness of T implies that its eigenvectors form a basis of C n , and hence (7.5) always admits a solution. In matrix form, we can rewrite (7.5) as S c = a, 10/18/06 so that 104 c = S -1 a,
c 2006 Peter J. Olver
where S = ( v1 v2 . . . vn ) is the (nonsingular) matrix whose columns are the eigenvectors. Q.E.D. Remark : Solutions in the incomplete cases rely on the Jordan canonical form. As with systems of differential equations, the formulas are more complicated, and will not be written out. Example 7.3. Consider the iterative system x(k+1) = with initial conditions x(0) = a, y (0) = b. x(k) , y (k) (7.7) The system can be rewritten in our matrix form (7.1), with T = .6 .2 .2 .6 , u(k) = a= a . b
3 5
x(k) +
1 5
y (k) ,
y (k+1) =
1 5
x(k) +
3 5
y (k) ,
(7.6)
Solving the characteristic equation det(T - I ) = 2 - 1.2 - .32 = 0 produces the eigenvalues 1 = .8, 2 = .4. We then solve the associated linear systems (T - j I )vj = 0 for the corresponding eigenvectors: 1 = .8 , v1 = 1 , 1 2 = .4 , v2 = -1 . 1
Therefore, the basic power solutions are u1 = ( .8)k
(k)
1 , 1
u2 = ( .4)k
(k)
-1 . 1 c1 ( .8)k - c2 ( .4)k , c1 ( .8)k + c2 ( .4)k b-a . 2
Theorem 7.2 tells us that the general solution is given as a linear combination, u(k) = c1 u1 + c2 u2 = c1 ( .8)k
(k) (k)
1 1
+ c2 ( .4)k
-1 1
=
where c1 , c2 are determined by the initial conditions: u(0) = c1 - c2 c1 + c2 = a , b and hence c1 = a+b , 2 c2 =
Therefore, the explicit formula for the solution to the initial value problem (7.67) is x(k) = ( .8)k a-b a+b + ( .4)k , 2 2 y (k) = ( .8)k a+b b-a + ( .4)k . 2 2
In particular, as k , the iterates u(k) 0 converge to zero at a rate governed by the dominant eigenvalue 1 = .8. Thus, (7.6) defines a stable iterative system. Figure 7.1 illustrates the cumulative effect of the iteration. The initial conditions consist of a large number of points on the unit circle x2 + y 2 = 1, which are successively mapped to points on progressively smaller and flatter ellipses, all converging towards the origin. 10/18/06 105
c 2006 Peter J. Olver
Figure 7.1.
Stable Iterative System.
Example 7.4. The Fibonacci numbers are defined by the second order iterative scheme u(k+2) = u(k+1) + u(k) , (7.8) with initial conditions u(0) = a, u(1) = b. (7.9) In short, to obtain the next Fibonacci number, add the previous two. The classical Fibonacci integers start with a = 0, b = 1; the next few are u(0) = 0, u(1) = 1, u(2) = 1, u(3) = 2, u(4) = 3, u(5) = 5, u(6) = 8, u(7) = 13, . . . . The Fibonacci integers occur in a surprising variety of natural objects, including leaves, flowers, and fruit, [46]. They were originally introduced by the Italian Renaissance mathematician Fibonacci (Leonardo of Pisa) as a crude model of the growth of a population of rabbits. In Fibonacci's model, the k th Fibonacci number u(k) measures the total number of pairs of rabbits at year k. We start the process with a single juvenile pair at year 0. Once a year, each pair of rabbits produces a new pair of offspring, but it takes a full year for a rabbit pair to mature enough to produce offspring of their own. Just as every higher order ordinary differential equation can be replaced by an equivalent first order system, so every higher order iterative equation can be replaced by a first
In general, an iterative system u(k+j) = T1 u(k+j-1) + + Tj u(k) in which the new iterate depends upon the preceding j values is said to have order j.
We ignore important details like the sex of the offspring.
10/18/06
106
c 2006
Peter J. Olver
order iterative system. In this particular case, we define the vector u(k) = u(k) u(k+1) R2,
and note that (7.8) is equivalent to the matrix system u(k+1) u(k+2) = 0 1 1 1 u(k) , u(k+1) or u(k+1) = T u(k) , where T = 0 1 1 1 .
To find the explicit formula for the Fibonacci numbers, we must determine the eigenvalues and eigenvectors of the coefficient matrix T . A straightforward computation produces 1+ 5 1- 5 1 = = 1.618034 . . . , 2 = = - .618034 . . . , 2 2 v1 =
-1+ 5 2
1
,
v2 =
-1- 5 2
1
.
Therefore, according to (7.4), the general solution to the Fibonacci system is u(k) = u u(k+1)
(k)
= c1
1+ 5 2
k
-1+ 5 2
1 + c2
+ c2
-1- 5 2
1- 5 2 a b
k
-1- 5 2
1
.
(7.10)
The initial data u(0) = c1
-1+ 5 2
1
1
=
uniquely specifies the coefficients 2 a + (1 + 5) b , c1 = 2 5 5) a + 2 b 2 5 1+ 5 2
k
2 a + (1 - 5) b c2 = - . 2 5 1- 5 2
k
The first entry of the solution vector (7.10) produces the explicit formula u
(k)
=
(-1 +
+
(1 +
5) a - 2 b 2 5
(7.11)
for the k th Fibonacci integer. It is a remarkable fact that, for every value of k, all the 5's cancel out, and the Binet formula does indeed produce the Fibonacci integers listed above. Another useful observation is that, since 5-1 1+ 5 0 < | 2 | = < 1 < 1 = , 2 2 10/18/06 107
c 2006 Peter J. Olver
for the k th Fibonacci number. For the particular initial conditions a = 0, b = 1, (7.11) reduces to the classical Binet formula k k 1- 5 1 1+ 5 (7.12) u(k) = - 2 2 5
Figure 7.2.
Fibonacci Iteration.
the terms involving k go to (and so the zero solution to this iterative system is unstable) 1 while the terms involving k go to zero. Therefore, even for k moderately large, the first 2 term in (7.11) is an excellent approximation (and one that gets more and more accurate with increasing k) to the k th Fibonacci number. A plot of the first 4 iterates, starting with the initial data consisting of equally spaced points on the unit circle, can be seen in Figure 7.2. As in the previous example, the circle is mapped to a sequence of progressively more eccentric ellipses; however, their major semi-axes become more and more stretched out, and almost all points end up going off to . 1 The dominant eigenvalue 1 = 2 1 + 5 = 1.618034 . . . is known as the golden ratio and plays an important role in spiral growth in nature, as well as in art, architecture and design, [46]. It describes the overall growth rate of the Fibonacci integers, and, in 1 fact, any sequence of Fibonacci numbers with initial conditions b = 2 1 - 5 a. -3 1 6 Example 7.5. Let T = 1 -1 -2 be the coefficient matrix for a three-1 -1 0 (k+1) (k) dimensional iterative system u = T u . Its eigenvalues and corresponding eigenvectors are 1 = - 2, 4 v1 = -2 , 1 10/18/06 2 = - 1 + i , 2- i v2 = -1 , 1 108 3 = - 1 - i , 2+ i v3 = -1 . 1
c 2006 Peter J. Olver
where b1 , b2 , b3 are arbitrary complex scalars. If we are only interested in real solutions, we can, as in the case of systems of differential equations, break up any complex solution into its real and imaginary parts, each of which constitutes a real solution. We begin by writing 2 = - 1 + i = 2 e3 i /4 , and hence
3 (- 1 + i )k = 2k/2 e3 k i /4 = 2k/2 cos 3 k + i sin 4 k . 4
Therefore, according to (7.4), the general complex solution to the iterative system is 2+ i 2- i 4 u(k) = b1 (- 2)k -2 + b2 (- 1 + i )k -1 + b3 (- 1 - i )k -1 , 1 1 1
Therefore, the complex solution 3 2 cos 3 k + sin 3 k 2 sin 4 k - cos 3 k 2- i 4 4 4 k/2 3 (- 1 + i )k -1 = 2k/2 + i 2 - cos 4 k - sin 3 k 4 1 3 3 cos k sin k
4 4
(7.13) where c1 , c2 , c3 are arbitrary real scalars, uniquely prescribed by the initial conditions.
is a combination of two independent real solutions. The complex conjugate eigenvalue 3 = - 1 - i leads, as before, to the complex conjugate solution -- and the same two real solutions. The general real solution u(k) to the system can be written as a linear combination of the three independent real solutions: 3 3 2 cos 3 k + sin 4 k 2 sin 3 k - cos 4 k 4 4 4 k/2 c1 (- 2)k -2 + c2 2k/2 + c3 2 , - cos 3 k - sin 3 k 4 4 1 cos 3 k sin 3 k
4 4
7.2. Stability.
With the solution formula (7.4) in hand, we are now in a position to understand the qualitative behavior of solutions to (complete) linear iterative systems. The most important case for applications is when all the iterates converge to 0. Definition 7.6. The equilibrium solution u = 0 to a linear iterative system (7.1) is called asymptotically stable if and only if all solutions u(k) 0 as k . Asymptotic stability relies on the following property of the coefficient matrix. Definition 7.7. A matrix T is called convergent if its powers converge to the zero matrix, T k O, meaning that the individual entries of T k all go to 0 as k . The equivalence of the convergence condition and stability of the iterative system follows immediately from the solution formula (7.2). Proposition 7.8. The linear iterative system u(k+1) = T u(k) has asymptotically stable zero solution if and only if T is a convergent matrix. 10/18/06 109
c 2006 Peter J. Olver
Proof : If T k O, and u(k) = T k a is any solution, then clearly u(k) 0 as k , (k) proving stability. Conversely, the solution uj = T k ej is the same as the j th column of T k . If the origin is asymptotically stable, then uj T k all tend to 0, proving that T k O.
(k)
0. Thus, the individual columns of Q.E.D.
To facilitate the analysis of convergence, we shall adopt a norm on our underlying vector space, R n or C n . The reader may be inclined to choose the Euclidean (or Hermitian) norm, but, in practice, the norm u
= max | u1 |, . . . , | un |
(7.14)
prescribed by the vector's maximal entry (in modulus) is usually much easier to work with. Convergence of the iterates is equivalent to convergence of their norms: u(k) 0 if and only if u(k) 0 as k .
The fundamental stability criterion for linear iterative systems relies on the size of the eigenvalues of the coefficient matrix. Theorem 7.9. A linear iterative system (7.1) is asymptotically stable if and only if all its (complex) eigenvalues have modulus strictly less than one: | j | < 1. Proof : Let us prove this result assuming that the coefficient matrix T is complete. (The proof in the incomplete case relies on the Jordan canonical form, and is outlined in the exercises.) If j is an eigenvalue such that | j | < 1, then the corresponding basis solution uj
(k)
= k vj tends to zero as k ; indeed, j uj
(k)
= k vj j
= | j |k vj
- 0
since
| j | < 1.
Therefore, if all eigenvalues are less than 1 in modulus, all terms in the solution formula (7.4) tend to zero, which proves asymptotic stability: u(k) 0. Conversely, if any eigenvalue satisfies | j | 1, then the solution u(k) = k vj does not tend to 0 as k , and j hence 0 is not asymptotically stable. Q.E.D. Consequently, the necessary and sufficient condition for asymptotic stability of a linear iterative system is that all the eigenvalues of the coefficient matrix lie strictly inside the unit circle in the complex plane: | j | < 1. Definition 7.10. The spectral radius of a matrix T is defined as the maximal modulus of all of its real and complex eigenvalues: (T ) = max { | 1 |, . . . , | k | }. We can then restate the Stability Theorem 7.9 as follows: Theorem 7.11. The matrix T is convergent if and only if its spectral radius is strictly less than one: (T ) < 1. 10/18/06 110
c 2006 Peter J. Olver
If T is complete, then we can apply the triangle inequality to (7.4) to estimate u(k) = c1 k v1 + + cn k vn 1 n
| 1 |k c1 v1 + + | n |k cn vn (T )k | c1 | v1 + + | cn | vn
(7.15) = C (T )k ,
for some constant C > 0 that depends only upon the initial conditions. In particular, if (T ) < 1, then u(k) C (T )k - 0 as k , (7.16) in accordance with Theorem 7.11. Thus, the spectral radius prescribes the rate of convergence of the solutions to equilibrium. The smaller the spectral radius, the faster the solutions go to 0. If T has only one largest (simple) eigenvalue, so | 1 | > | j | for all j > 1, then the first term in the solution formula (7.4) will eventually dominate all the others: k v1 1 k vj for j > 1 and k 0. Therefore, provided that c1 = 0, the solution (7.4) has the j asymptotic formula (7.17) u(k) c1 k v1 , 1 and so most solutions end up parallel to v1 . In particular, if | 1 | = (T ) < 1, such a solution approaches 0 along the direction of the dominant eigenvector v1 at a rate governed by the modulus of the dominant eigenvalue. The exceptional solutions, with c1 = 0, tend to 0 at a faster rate, along one of the other eigendirections. In practical computations, one rarely observes the exceptional solutions. Indeed, even if the initial condition does not involve the dominant eigenvector, round-off error during the iteration will almost inevitably introduce a small component in the direction of v1 , which will, if you wait long enough, eventually dominate the computation. Warning: The inequality (7.15) only applies to complete matrices. In the general case, one can prove that the solution satisfies the slightly weaker inequality u(k) C k for all k 0, where > (T ) (7.18)
is any number larger than the spectral radius, while C > 0 is a positive constant (whose value may depend on how close is to ). Example 7.12. According to Example 7.5, the matrix -3 T = 1 -1 1 6 -1 -2 -1 0 1 = - 2, 2 = - 1 + i , 3 = - 1 - i .
has eigenvalues
Since | 1 | = 2 > | 2 | = | 3 | = 2 , the spectral radius is (T ) = | 1 | = 2. We conclude that T is not a convergent matrix. As the reader can check, either directly, or from the solution formula (7.13), the vectors u(k) = T k u(0) obtained by repeatedly multiplying any nonzero initial vector u(0) by T rapidly go off to , at a rate roughly equal to (T )k = 2k . 10/18/06 111
c 2006 Peter J. Olver
On the other hand, the matrix T =-1T = - 3 1
1 3 1 3 1 -3
1 3 1 3
-2
2 3
with eigenvalues
0
1 = 2 , 3 1 2 = 3 - 3 =
1 3
+
1 3 1 3
i, i,
2 has spectral radius (T ) = 3 , and hence is a convergent matrix. According to (7.17), if we write the initial data u(0) = c1 v1 +c2 v2 +c3 v3 as a linear combination of the eigenvectors, 2 k then, provided c1 = 0, the iterates have the asymptotic form u(k) c1 - 3 v1 , where T v1 = ( 4, -2, 1 ) is the eigenvector corresponding to the dominant eigenvalue 1 = - 2 . 3 Thus, for most initial vectors, the iterates end up decreasing in length by a factor of almost 2 exactly 3 , eventually becoming parallel to the dominant eigenvector v1 . This is borne out T by a sample computation: starting with u(0) = ( 1, 1, 1 ) , the first ten iterates are
-.0121 .0061 , -.0030
-.0936 .0462 , -.0231
-.0036 -.0054 -.0081 .0040 , .0027 , .0018 , -.0009 -.0013 -.0020
-.0275 -.0416 -.0627 .0312 , .0208 , .0138 , -.0069 -.0105 -.0158
-.0024 .0012 . -.0006
-.0182 .0091 , -.0046
7.3. Matrix Norms.
The convergence of a linear iterative system is governed by the spectral radius or largest eigenvalue (in modulus) of the coefficient matrix. Unfortunately, finding accurate approximations to the eigenvalues of most matrices is a nontrivial computational task. Indeed, all practical numerical algorithms rely on some form of iteration. But using iteration to determine the spectral radius defeats the purpose, which is to predict the behavior of the iterative system in advance! In this section, we present two alternative approaches for directly investigating convergence and stability issues. Matrix norms form a natural class of norms on the vector space of n n matrices and can, in many instances, be used to establish convergence with a minimal effort. Matrix Norms We work exclusively with real n n matrices in this section, although the results straightforwardly extend to complex matrices. We begin by fixing a norm on R n . The norm may or may not come from an inner product -- this is irrelevant as far as the construction goes. Each norm on R n will naturally induce a norm on the vector space Mnn of all n n matrices. Roughly speaking, the matrix norm tells us how much a linear transformation stretches vectors relative to the given norm. 10/18/06 112
c 2006 Peter J. Olver
Theorem 7.13. If
is any norm on R n , then the quantity A = max { A u | u =1} (7.19)
defines a norm on Mnn , known as the natural matrix norm. Proof : First note that A < , since the maximum is taken on a closed and bounded subset, namely the unit sphere S1 = { u = 1 } for the given norm. To show that (7.19) defines a norm, we need to verify the three basic axioms of Definition 5.8. Non-negativity, A 0, is immediate. Suppose A = 0. This means that, for every unit vector, A u = 0, and hence A u = 0 whenever u = 1. If 0 = v R n is any nonzero vector, then u = v/r, where r = v , is a unit vector, so A v = A(r u) = r A u = 0. (7.20)
Therefore, A v = 0 for every v R n , which implies A = O is the zero matrix. This serves to prove the positivity property. As for homogeneity, if c R is any scalar, c A = max { c A u } = max { | c | A u } = | c | max { A u } = | c | A . Finally, to prove the triangle inequality, we use the fact that the maximum of the sum of quantities is bounded by the sum of their individual maxima. Therefore, since the norm on R n satisfies the triangle inequality, A + B = max { A u + B u } max { A u + B u } max { A u } + max { B u } = A + B . Q.E .D.
The property that distinguishes a matrix norm from a generic norm on the space of matrices is the fact that it also obeys a very useful product inequality. Theorem 7.14. A natural matrix norm satisfies Av A Furthermore, AB A B , for all A, B Mnn . u (7.22) = 1. v , for all A Mnn , v Rn. (7.21)
Proof : Note first that, by definition A u A for all unit vectors Then, letting v = r u where u is a unit vector and r = v , we have A v = A(r u) = r A u r A = v A B = max { A B u } = max { A (B u) } max { A Bu }= A A ,
proving the first inequality. To prove the second, we apply the first to compute
max { B u } = A
B .
Q.E .D.
Remark : In general, a norm on the vector space of n n matrices is called a matrix norm if it also satisfies the multiplicative inequality (7.22). Most, but not all, matrix norms used in applications come from norms on the underlying vector space. 10/18/06 113
c 2006 Peter J. Olver
The multiplicative inequality (7.22) implies, in particular, that A2 A ity is not necessarily valid. More generally:
2
; equal-
Proposition 7.15. If A is a square matrix, then Ak A k . In particular, if A < 1, then Ak 0 as k , and hence A is a convergent matrix: Ak O. The converse is not quite true; a convergent matrix does not necessarily have matrix norm less than 1, or even 1 -- see Example 7.20 below. An alternative proof of Proposition 7.15 can be based on the following useful estimate: Theorem 7.16. The spectral radius of a matrix is bounded by its matrix norm: (A) A . (7.23)
Proof : If is a real eigenvalue, and u a corresponding unit eigenvector, so that A u = u with u = 1, then A u = u = | | u = | |. Since A || A . (7.24)
is the maximum of A u over all possible unit vectors, this implies that (7.25)
If all the eigenvalues of A are real, then the spectral radius is the maximum of their absolute values, and so it too is bounded by A , proving (7.23). If A has complex eigenvalues, then we need to work a little harder to establish (7.25). (This is because the matrix norm is defined by the effect of A on real vectors, and so we cannot directly use the complex eigenvectors to establish the required bound.) Let = r e i be a complex eigenvalue with complex eigenvector z = x + i y. Define m = min Re (e i z) = (cos ) x - (sin ) y 0 2 . (7.26)
Since the indicated subset is a closed curve (in fact, an ellipse) that does not go through the origin, m > 0. Let 0 denote the value of the angle that produces the minimum, so m = (cos 0 ) x - (sin 0 ) y = Re e i 0 z Define the real unit vector u= Then Au = Re e i 0 z m = (cos 0 ) x - (sin 0 ) y , m so that u = 1. .
1 1 r Re e i 0 A z = Re e i 0 r e i z = Re e i (0 +) z . m m m r m
Therefore, keeping in mind that m is the minimal value in (7.26), A Au = Re e i (0 +) z r = | |, (7.27) Q.E.D.
c 2006 Peter J. Olver
and so (7.25) also holds for complex eigenvalues. 10/18/06 114
Explicit Formulae Let us now determine the explicit formulae for the matrix norms induced by our most important vector norms on R n . The simplest to handle is the norm v
= max { | v1 |, . . . , | vn | }.
Definition 7.17. The ith absolute row sum of a matrix A is the sum of the absolute values of the entries in the ith row:
n
si = | ai1 | + + | ain | = Proposition 7.18. absolute row sum: A
j =1
| aij |.
(7.28)
Proof : Let s = max{s1 , . . . , sn } denote the right hand side of (7.29). Given any v R n , we compute n n max | aij vj | aij vj A v = max j =1 j =1 n max | aij | max | vj | = s v .
j =1
The matrix norm of a matrix A is equal to its maximal n = max{s1 , . . . , sn } = max | aij | 1 i n . (7.29)
j =1
In particular, by specializing to v = 1, we deduce that A s. On the other hand, suppose the maximal absolute row sum occurs at row i, so
n
si =
j =1
| aij | = s.
(7.30)
Let u R n be the specific vector that has the following entries: uj = + 1 if aij > 0, while uj = - 1 if aij < 0. Then u = 1. Moreover, since aij uj = | aij |, the ith entry of A u is equal to the ith absolute row sum (7.30). This implies that A
Au
s.
Q.E .D.
Combining Propositions 7.15 and 7.18, we have established the following convergence criterion. Corollary 7.19. If all the absolute row sums of A are strictly less than 1, then A < 1 and hence A is a convergent matrix. 10/18/06 115
c 2006 Peter J. Olver
Example 7.20. Consider the symmetric matrix A = lute row sums are
1 2
+ -1 = 3 A
5 6
1 , -3 + 5 7 6 , 12
1 4
=
5 6
7 12
, so
-
1 2 1 3
-1 3
1 4
. Its two abso-
= max
=
Since the norm is less than 1, A is a convergent matrix. Indeed, its eigenvalues are 9 - 73 9 + 73 .7310 . . . , 2 = .0190 . . . , 1 = 24 24 and hence the spectral radius is 9 + 73 .7310 . . . , (A) = 24 which is slightly smaller than its norm. The row sum test for convergence is not always conclusive. For example, the matrix A=
1 2 3 5
.83333 . . . .
On the other hand, its eigenvalues are (15 601 )/40, and hence its spectral radius is 15 + 601 (A) = .98788 . . . , 40 which implies that A is (just barely) convergent, even though its maximal row sum is larger than 1. The matrix norm associated with the Euclidean norm given by largest singular value. v
2
-
-3 5
1 4
has matrix norm
A
=
11 10
> 1.
(7.31)
=
2 2 v1 + + vn is
Proposition 7.21. The matrix norm corresponding to the Euclidean norm equals the maximal singular value: A
2
= 1 = max {1 , . . . , r },
r = rank A > 0,
while
O
2
= 0.
(7.32)
Unfortunately, as we discovered in Example 7.20, matrix norms are not a foolproof test of convergence. There exist convergent matrices such that (A) < 1 and yet have matrix norm A 1. In such cases, the matrix norm is not able to predict convergence of the iterative system, although one should expect the convergence to be quite slow. Although such pathology might show up in the chosen matrix norm, it turns out that one can always rig up some matrix norm for which A < 1. This follows from a more general result, whose proof can be found in [40]. Theorem 7.22. Let A have spectral radius (A). If > 0 is any positive number, then there exists a matrix norm such that (A) A < (A) + . (7.33) Corollary 7.23. If A is a convergent matrix, then there exists a matrix norm such that A < 1. 10/18/06 116
c 2006 Peter J. Olver
Proof : By definition, A is convergent if and only if (A) < 1. Choose > 0 such that (A) + < 1. Any norm that satisfies (7.33) has the desired property. Q.E.D. Remark : Based on the accumulated evidence, one might be tempted to speculate that the spectral itself radius defines a matrix norm. Unfortunately, this is not the case. For 0 1 has zero spectral radius, (A) = 0, in violation example, the nonzero matrix A = 0 0 of a basic norm axiom.
7.4. Iterative Solution of Linear Algebraic Systems.
In this section, we return to the most basic problem in linear algebra: solving the linear algebraic system A u = b, (7.34) consisting of n equations in n unknowns. We assume that the coefficient matrix A is nonsingular, and so the solution u = A-1 b is unique. We will introduce several popular iterative methods that can be used to approximate the solution for certain classes of coefficient matrices. The resulting algorithms will provide an attractive alternative to Gaussian Elimination, particularly when dealing with the large, sparse systems that arise in the numerical solution to differential equations. One major advantage of an iterative technique is that it (typically) produces progressively more and more accurate approximations to the solution, and hence, by prolonging the iterations, can, at least in principle, compute the solution to any desired order of accuracy. Moreover, even performing just a few iterations may produce a reasonable approximation to the true solution -- in stark contrast to Gaussian Elimination, where one must continue the algorithm through to the bitter end before any useful information can be extracted. A partially completed Gaussian Elimination is of scant use! A significant weakness is that iterative schemes are not universally applicable, and their design relies upon the detailed structure of the coefficient matrix. We shall be attempting to solve the linear system (7.34) by replacing it with an iterative system of the form u(k+1) = T u(k) + c, u(0) = u0 , (7.35)
in which T is an n n matrix and c a vector. This represents a slight generalization of our earlier iterative system (7.1), in that the right hand side is now an affine function of u(k) . Suppose that the solutions to the affine iterative system converge: u(k) u as k . Then, by taking the limit of both sides of (7.35), we discover that the limit point u solves the fixed-point equation u = T u + c. (7.36) Thus, we need to design our iterative system so that (a) te solution to the fixed-point system u = T u + c coincides with the solution to the original system A u = b, and (b) the iterates defined by (7.35) are known to converge to the fixed point. Before exploring these issues in depth, let us look at a simple example. 10/18/06 117
c 2006 Peter J. Olver
Example 7.24. Consider the linear system 3 x + y - z = 3, x - 4 y + 2 z = -1, - 2 x - y + 5 z = 2, 3 b = -1 . 2 (7.37)
One easy way to convert a linear system into a fixed-point form is to rewrite it as u = I u - A u + A u = ( I - A)u + b = T u + c, In the present case, -2 -1 1 5 -2 , T = I - A = -1 2 1 -4 3 c = b = -1 . 2 where T = I - A, c = b.
which has the vectorial form A u = b, with x 3 1 -1 y , , 1 -4 u= 2 A= z -2 -1 5
The resulting iterative system u(k+1) = T u(k) + c has the explicit form x(k+1) = - 2 x(k) - y (k) + z (k) + 3, z
(k+1)
y (k+1) = - x(k) + 5 y (k) - 2 z (k) - 1, = 2x
(k)
(7.38)
+y
(k)
- 4z
(k)
+ 2.
Another possibility is to solve the first equation in (7.37) for x, the second for y, and the third for z, so that x = -1 y + 3
1 3
z + 1,
y=
1 4
x+
1 2
z+ 1, 4
z=
2 5
x+
1 5
2 y+ 5.
The resulting equations have the form of a fixed-point system 1 T =4
2 5 1 3 1 2 1 5
u = T u + c,
in which
0 -1 3
1 5
0
1 3 1 2
0
The corresponding iteration u(k+1) = T u(k) + c takes the explicit form x(k+1) = - 1 y (k) + 3 y (k+1) = z (k+1) =
1 4 2 5
,
1 1 c = 4 .
2 5
z (k) + 1, z (k) + 1 , 4
2 y (k) + 5 .
x(k) + x(k) +
(7.39)
Do the resulting iterative schemes converge to the solution x = y = z = 1? The results, starting with initial guess u(0) = (0, 0, 0), are tabulated as follows. 10/18/06 118
c 2006 Peter J. Olver
k 0 1 2 3 4 5 6 7 8 9 10 11
u(k+1) = T u(k) + b 0 3 0 15 30 261 870 6069 22500 145743 571980 3522555 0 0 -1 2 - 13 -1 - 64 -7 - 322 -4 - 1633 - 244 - 7939 - 133 - 40300 - 5665 - 196240 - 5500 - 992701 - 129238 - 4850773 - 184261 - 24457324 - 2969767
u(k+1) = T u(k) + c 0 1 1.05 1.05 1.0075 1.005 .9986 1.0004 .9995 1.0001 .9999 1.0000 0 .25 .7 .9375 .9925 1.00562 1.002 1.0012 1.0000 1.0001 .9999 1.0000 0 .4 .85 .96 1.0075 1.0015 1.0031 .9999 1.0004 .9998 1.0001 1.0000
For the first scheme, the answer is clearly no -- the iterates become wilder and wilder. Indeed, this occurs no matter how close the initial guess u(0) is to the actual solution -- unless u(0) = u happens to be exactly equal. In the second case, the iterates do converge to the solution, and it does not take too long, even starting from a poor initial guess, to obtain a reasonably accurate approximation. Of course, in such a simple example, it would be silly to use iteration, when Gaussian Elimination can be done by hand and produces the solution almost immediately. However, we use the small examples for illustrative purposes, bringing the full power of iterative schemes to bear on the large linear systems arising in applications. The convergence of solutions to (7.35) to the fixed point u is based on the behavior of the error vectors e(k) = u(k) - u , (7.40)
which measure how close the iterates are to the true solution. Let us find out how the successive error vectors are related. We compute e(k+1) = u(k+1) - u = (T u(k) + a) - (T u + a) = T (u(k) - u ) = T e(k) , e(k+1) = T e(k) , e(k) = T k e(0) . Now, the solutions to (7.35) converge to the fixed point, u(k) u , if and only if the error vectors converge to zero: e(k) 0 as k . Our analysis of linear iterative systems, as summarized in Proposition 7.8, establishes the following basic convergence result. 10/18/06 119
c 2006 Peter J. Olver
showing that the error vectors satisfy a linear iterative system
(7.41)
with the same coefficient matrix T . Therefore, they are given by the explicit formula
Proposition 7.25. The affine iterative system (7.35) will converge to the solution to the fixed point equation (7.36) if and only if T is a convergent matrix: (T ) < 1. The spectral radius (T ) of the coefficient matrix will govern the speed of convergence. Therefore, our main goal is to construct an iterative scheme whose coefficient matrix has as small a spectral radius as possible. At the very least, the spectral radius must be less than 1. For the two iterative schemes presented in Example 7.24, the spectral radii of the coefficient matrices are found to be (T ) 4.9675, ( T ) = .5.
Therefore, T is not a convergent matrix, which explains the wild behavior of its iterates, 1 whereas T is convergent, and one expects the error to roughly decrease by a factor of 2 at each step. The Jacobi Method The first general iterative scheme for solving linear systems is based on the same simple idea used in our illustrative Example 7.24. Namely, we solve the ith equation in the system A u = b, which is
n
aij uj = bi ,
j =1
for the ith variable ui . To do this, we need to assume that all the diagonal entries of A are nonzero: aii = 0. The result is 1 ui = - aii where a - ij , aii tij = 0,
n j=1 j=i
b aij uj + i = aii
n
tij uj + ci ,
j =1
(7.42)
i = j, i = j,
and
ci =
bi . aii
(7.43)
The result has the form of a fixed-point system u = T u + c, and forms the basis of the Jacobi method (7.44) u(k+1) = T u(k) + c, u(0) = u0 , named after the influential nineteenth century German analyst Carl Jacobi. The explicit form of the Jacobi iterative scheme is n b 1 (k) (k+1) aij uj + i . ui = - (7.45) aii j=1 aii
j=i
It is instructive to rederive the Jacobi method in a direct matrix form. We begin by decomposing the coefficient matrix A=L+D+U 10/18/06 120
c 2006
(7.46)
Peter J. Olver
into the sum of a strictly lower triangular matrix L, a diagonal matrix D, and a strictly upper triangular matrix U , each of which is uniquely specified. For example, when 3 1 -1 (7.47) 2, A = 1 -4 -2 -1 5 the decomposition (7.46) yields 0 0 0 0 0, L= 1 -2 -1 0 3 D = 0 0 0 0 -4 0 , 0 5 0 U = 0 0 1 -1 0 2. 0 0
Warning: The L, D, U in the elementary additive decomposition (7.46) have nothing to do with the L, D, U appearing in factorizations arising from Gaussian Elimination. The latter play no role in the iterative solution methods considered here. We then rewrite the system A u = (L + D + U ) u = b in the alternative form D u = - (L + U ) u + b.
The Jacobi fixed point equations (7.42) amounts to solving for u = T u + c, where T = - D-1 (L + U ), c = D-1 b. (7.48)
For the example (7.47), we recover the Jacobi iteration matrix as 1 T = - D-1 (L + U ) = 4
2 5
0 -1 3
1 5
0
1 3 1 2
0
Deciding in advance whether or not the Jacobi method will converge is not easy. However, it can be shown that Jacobi iteration is guaranteed to converge when the original coefficient matrix has large diagonal entries, in accordance with Definition 6.25. Theorem 7.26. If A is strictly diagonally dominant, then the associated Jacobi iteration scheme converges. Proof : We shall prove that T < 1, and so Corollary 7.19 implies that T is a convergent matrix. The absolute row sums of the Jacobi matrix T = - D-1 (L + U ) are, according to (7.43), n n 1 | a | < 1, si = | tij | = (7.49) | aii | j=1 ij j =1
j=i
.
because A is strictly diagonally dominant. Thus, result follows. 10/18/06 121
T
= max{s1 , . . . , sn } < 1, and the Q.E.D.
c 2006 Peter J. Olver
Example 7.27. Consider the linear system 4 x + y + w = 1, x + 4 y + z + v = 2, y + 4 z + w = - 1, x + z + 4 w + v = 2, y + w + 4 v = 1. The Jacobi method solves the respective equations for x, y, z, w, v, leading to the iterative scheme 1 x(k+1) = - 1 y (k) - 1 w(k) + 4 , 4 4
1 1 y (k+1) = - 1 x(k) - 4 z (k) - 1 v (k) + 2 , 4 4 1 z (k+1) = - 1 y (k) - 1 w(k) - 4 , 4 4
1 1 w(k+1) = - 1 x(k) - 4 z (k) - 1 v (k) + 2 , 4 4 1 v (k+1) = - 1 y (k) - 1 w(k) + 4 . 4 4
The coefficient matrix of the original system, 4 1 A = 0 1 0 1 4 1 0 1 0 1 4 1 0 1 0 1 4 1 0 1 0 , 1 4
is diagonally dominant, and so we are guaranteed that the Jacobi iterations will eventually converge to the solution. Indeed, the Jacobi scheme takes the iterative form (7.48), with -1 0 4 1 T = 0 -4 -1 0 4 0 -1 4 0 -1 4 0 -1 4 -1 4 0 0
1 -4 1 -4 1 -4
0 0
Note that T = 3 < 1, validating convergence of the scheme. Thus, to obtain, say, 4 four decimal place accuracy in the solution, we estimate that it would take less than log(.5 10-4 )/ log .75 34 iterates, assuming a moderate initial error. But the matrix norm always underestimates the true rate of convergence, as prescribed by the spectral radius (T ) = .6124, which would imply about log(.5 10-4 )/ log .6124 20 iterations to obtain the desired accuracy. Indeed, starting with the initial guess x(0) = y (0) = z (0) = w(0) = v (0) = 0, the Jacobi iterates converge to the exact solution x = - .1, y = .7, z = - .6, w = .7, v = - .1,
-1 4 0 , 1 -4 0
0
c= -
1 4 1 2 1 4 1 2 1 4
.
to within four decimal places in exactly 20 iterations. 10/18/06 122
c 2006 Peter J. Olver
The GaussSeidel Method The GaussSeidel method relies on a slightly more refined implementation of the Jacobi process. To understand how it works, it will help to write out the Jacobi iteration scheme (7.44) in full detail: u1 u2 u3
(k+1)
= = t21 u1 . . .
(k) (k) (k)
(k+1)
t12 u2 + t13 u3 + + t1,n-1 un-1 + t1n u(k) + c1 , n
(k)
(k)
(k)
(k)
(k+1)
+ t23 u3 + + t2,n-1 un-1 + t2n u(k) + c2 , n
(k)
(k)
(k)
= t31 u1 + t32 u2 . . .
. . .
..
(k)
+ t3,n-1 un-1 + t3n u(k) + c3 , n .. . . . .
(7.50)
.
(k) (k)
u(k+1) = tn1 u1 + tn2 u2 + tn3 u3 + + tn,n-1 un-1 n
+ cn ,
where we are explicitly noting the fact that the diagonal entries of T vanish. Observe that we are using the entries of u(k) to compute all of the updated values of u(k+1) . Presumably, if the iterates u(k) are converging to the solution u , then their individual (k+1) should be a better approximation to u entries are also converging, and so each uj j than uj
(k)
is. Therefore, if we begin the k th Jacobi iteration by computing u1
(k+1)
using the
(k)
first equation, then we are tempted to use this new and improved value to replace u1 each of the subsequent equations. In particular, we employ the modified equation u2
(k+1)
in
= t21 u1
(k+1)
+ t23 u3 + + t1n u(k) + c2 n
(k)
to update the second component of our iterate. This more accurate value should then be (k+1) , and so on. used to update u3 The upshot of these considerations is the GaussSeidel method i = 1, . . . , n, (7.51) named after Gauss (as usual!) and the German astronomer/mathematician Philipp von Seidel. At the k th stage of the iteration, we use (7.51) to compute the revised entries (k+1) (k+1) , . . . , u(k+1) in their numerical order. Once an entry has been updated, the , u2 u1 n new value is immediately used in all subsequent computations. ui = ti1 u1 + + ti,i-1 ui-1 Example 7.28. For the linear system 3 x + y - z = 3, x - 4 y + 2 z = -1, - 2 x - y + 5 z = 2,
(k+1) (k+1) (k+1)
+ ti,i+1 ui+1 + + tin u(k) + ci , n
(k)
the Jacobi iteration method was given in (7.39). To construct the corresponding Gauss Seidel scheme we use updated values of x, y and z as they become available. Explicitly, x(k+1) = - 1 y (k) + 3 y (k+1) = z (k+1) = 10/18/06
1 4 2 5 1 (k) + 1, 3 z x(k+1) + 1 z (k) + 1 , 2 4 (k+1) 1 (k+1) 2 x +5y + 5.
(7.52)
123
c 2006
Peter J. Olver
The resulting iterates starting with u(0) = 0 are 1.0222 1.1333 1.0000 u(3) = 1.0306 , u(2) = .9833 , u(1) = .5000 , 1.0150 1.0500 .9000 1.0001 1.0000 .9977 u(7) = 1.0000 , u(6) = .9994 , u(5) = .9990 , 1.0001 .9999 .9989
u(4)
u(8)
and have converged to the solution, to 4 decimal place accuracy, after only 8 iterations -- as opposed to the 11 iterations required by the Jacobi method.
.9948 = 1.0062 , .9992 1.0000 = 1.0000 , 1.0000
The GaussSeidel iteration scheme is particularly suited to implementation on a serial (k) computer, since one can immediately replace each component ui by its updated value (k+1) , thereby also saving on storage in the computer's memory. In contrast, the Jacobi ui scheme requires us to retain all the old values u(k) until the new approximation u(k+1) has been computed. Moreover, GaussSeidel typically (although not always) converges faster than Jacobi, making it the iterative algorithm of choice for serial processors. On the other hand, with the advent of parallel processing machines, variants of the parallelizable Jacobi scheme have recently been making a comeback. What is GaussSeidel really up to? Let us rewrite the basic iterative equation (7.51) by multiplying by aii and moving the terms involving u(k+1) to the left hand side. In view of the formula (7.43) for the entries of T , the resulting equation is ai1 u1
(k+1)
+ + ai,i-1 ui-1
(k+1)
+ aii ui
(k+1)
= - ai,i+1 ui+1 - - ain u(k) + bi . n (7.53)
(k)
In matrix form, taking (7.46) into account, this reads (L + D)u(k+1) = - U u(k) + b,
and so can be viewed as a linear system of equations for u(k+1) with lower triangular coefficient matrix L + D. Note that the fixed point of (7.53), namely the solution to (L + D) u = - U u + b, coincides with the solution to the original system A u = (L + D + U ) u = b. In other words, the GaussSeidel procedure is merely implementing Forward Substitution to solve the lower triangular system (7.53) for the next iterate: u(k+1) = - (L + D)-1 U u(k) + (L + D)-1 b. The latter is in our more usual iterative form u(k+1) = T u(k) + c, where T = - (L + D)-1 U, c = (L + D)-1 b. (7.54)
Consequently, the convergence of the GaussSeidel iterates is governed by the spectral radius of the coefficient matrix T . 10/18/06 124
c 2006 Peter J. Olver
Therefore, the GaussSeidel matrix is
Returning to Example 7.28, we have 3 3 1 -1 1 , 1 -4 L+D = 2 A= -2 -2 -1 5 0 -1 T = - (L + D) U = 0 0
0 0 -4 0 , -1 5 -.3333 -.0833
0 0 U= 0
1 -1 0 2. 0 0
-.1500
.2500
Its eigenvalues are 0 and .0833.2444 i , and hence its spectral radius is ( T ) .2582. This is roughly the square of the Jacobi spectral radius of .5, which tell us that the GaussSeidel iterations will converge about twice as fast to the solution. This can be verified by more extensive computations. Although examples can be constructed where the Jacobi method converges faster, in many practical situations GaussSeidel tends to converge roughly twice as fast as Jacobi. Completely general conditions guaranteeing convergence of the GaussSeidel method are hard to establish. But, like the Jacobi scheme, it is guaranteed to converge when the original coefficient matrix is strictly diagonally dominant. Theorem 7.29. If A is strictly diagonally dominant, then the GaussSeidel iteration scheme for solving A u = b converges. Proof : Let e(k) = u(k) - u denote the k th GaussSeidel error vector. As in (7.41), the error vectors satisfy the linear iterative system e(k+1) = T e(k) , but a direct estimate of T is not so easy. Instead, let us write out the linear iterative system in components: ei Let m(k) = e(k)
(k+1)
.3333 .5833 .
= ti1 e1
(k+1)
+ + ti,i-1 ei-1
(k+1)
+ ti,i+1 ei+1 + + tin e(k) . n
(k)
(k)
(7.55) (7.56)
= max | e1 |, . . . , | e(k) | n
denote the norm of the k th error vector. To prove convergence, e(k) 0, it suffices to show that m(k) 0 as k . We claim that diagonal dominance of A implies that m(k+1) s m(k) , where s= T
<1
(7.57)
denotes the matrix norm of the Jacobi matrix (not the GaussSeidel matrix), which, by (7.49), is less than 1. We infer that m(k) sk m(0) 0 as k , demonstrating the theorem. To prove (7.57), we use induction on i = 1, . . . , n. Our induction hypothesis is | ej
(k+1)
| s m(k) < m(k) | ej | m(k)
(k)
for
j = 1, . . . , i - 1. j = 1, . . . , n.
c 2006 Peter J. Olver
(When i = 1, there is no assumption.) Moreover, by (7.56), for all 125
10/18/06
We use these two inequalities to estimate | ei | ei
(k+1)
(k+1)
| from (7.55):
(k)
| | ti1 | | e1
(k+1)
| + + | ti,i-1 | | ei-1 | + | ti,i+1 | | ei+1 | + + | tin | | e(k) | n m(k) s m(k) ,
(k+1)
(k+1)
| ti1 | + + | tin |
which completes the induction step. As a result, the maximum m(k+1) = max | e1 |, . . . , | e(k+1) | n s m(k) Q.E.D.
also satisfies the same bound, and hence (7.57) follows.
Example 7.30. For the linear system considered in Example 7.27, the GaussSeidel iterations take the form
1 x(k+1) = - 1 y (k) - 1 w(k) + 4 , 4 4
1 1 1 y (k+1) = - 4 x(k+1) - 4 z (k) - 1 v (k) + 2 , 4 1 1 z (k+1) = - 1 y (k+1) - 4 w(k) - 4 , 4
1 w(k+1) = - 1 x(k+1) - 4 z (k+1) - 1 v (k) + 1 , 4 4 2 1 1 1 v (k+1) = - 4 y (k+1) - 4 w(k+1) + 4 .
Starting with x(0) = y (0) = z (0) = w(0) = v (0) = 0, the GaussSeidel iterates converge to the solution x = - .1, y = .7, z = - .6, w = .7, v = - .1, to four decimal places in 11 iterations, again roughly twice as fast as the Jacobi scheme. Indeed, the convergence rate is governed by the corresponding GaussSeidel matrix T , which is 4 1 0 1 0 0 4 1 0 1 0 0 4 1 0 0 0 0 4 1 -1 0 0 0 0 0 0 0 0 4 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 = 0 1 0 0 0 -.2500 .0625 -.0156 .0664 -.0322 0 -.2500 .0625 -.0156 .0664 -.2500 .0625 -.2656 .1289 -.0479 0 -.2500 .0625 . -.2656 .1289
Its spectral radius is ( T ) = .3936, which is, as in the previous example, approximately the square of the spectral radius of the Jacobi coefficient matrix, which explains the speed up in convergence. Successive OverRelaxation (SOR) As we know, the smaller the spectral radius (or matrix norm) of the coefficient matrix, the faster the convergence of the iterative scheme. One of the goals of researchers in numerical linear algebra is to design new methods for accelerating the convergence. In his 1950 thesis, the American mathematician David Young discovered a simple modification of the Jacobi and GaussSeidel methods that can, in favorable situations, lead to a dramatic speed up in the rate of convergence. The method, known as successive over-relaxation, and often abbreviated as SOR, has become the iterative method of choice in many modern applications, [13, 47]. In this subsection, we provide a brief overview. 10/18/06 126
c 2006 Peter J. Olver
In practice, finding the optimal iterative algorithm to solve a given linear system is as hard as solving the system itself. Therefore, researchers have relied on a few tried and true techniques for designing iterative schemes that can be used in the more common applications. Consider a linear algebraic system A u = b. Every decomposition of the coefficient matrix into the difference of two matrices leads to an equivalent system of the form M u = N u + b. (7.59) Provided that M is nonsingular, we can rewrite the system in the fixed point form u = M -1 N u + M -1 b = T u + c, where T = M -1 N, c = M -1 b. A=M -N (7.58)
Now, we are free to choose any such M , which then specifies N = A-M uniquely. However, for the resulting iterative scheme u(k+1) = T u(k) + c to be practical we must arrange that (a) T = M -1 N is a convergent matrix, and (b) M can be easily inverted. The second requirement ensures that the iterative equations M u(k+1) = N u(k) + b (7.60)
can be solved for u(k+1) with minimal computational effort. Typically, this requires that M be either a diagonal matrix, in which case the inversion is immediate, or upper or lower triangular, in which case one employs Back or Forward Substitution to solve for u(k+1) . With this in mind, we now introduce the SOR method. It relies on a slight generalization of the GaussSeidel decomposition (7.53) of the matrix into lower plus diagonal and upper triangular parts. The starting point is to write A = L + D + U = L + D - ( - 1) D - U , (L + D)u = ( - 1) D - U u + b. (7.61)
where 0 = is an adjustable scalar parameter. We decompose the system A u = b as (7.62)
It turns out to be slightly more convenient to divide (7.62) through by and write the resulting iterative system in the form ( L + D)u(k+1) = (1 - ) D - U u(k) + b, (7.63)
where = 1/ is called the relaxation parameter . Assuming, as usual, that all diagonal entries of A are nonzero, the matrix L + D is an invertible lower triangular matrix, and so we can use Forward Substitution to solve the iterative system (7.63) to recover u(k+1) . The explicit formula for its ith entry is ui
(k+1)
= ti1 u1
(k+1)
+ + ti,i-1 ui-1
(k)
(k+1)
+ ti,i+1 ui+1 + + tin u(k) + ci , n 127
c 2006
+ (1 - ) ui
(k)
+
(7.64)
10/18/06
Peter J. Olver
where tij and ci denote the original Jacobi values (7.43). As in the GaussSeidel approach, (k+1) we update the entries ui in numerical order i = 1, . . . , n. Thus, to obtain the SOR scheme (7.64), we merely multiply the right hand side of the GaussSeidel scheme (7.51) (k) by the adjustable relaxation parameter and append the diagonal term (1 - ) ui . In particular, if we set = 1, then the SOR method reduces to the GaussSeidel method. Choosing < 1 leads to an under-relaxed method, while > 1, known as over-relaxation, is the choice that works in most practical instances. To analyze the SOR scheme in detail, we rewrite (7.63) in the fixed point form u(k+1) = T u(k) + c , where T = ( L + D)-1 (1 - ) D - U , The rate of is to choose as possible. convergence c = ( L + D)-1 b. (7.66) (7.65)
convergence is governed by the spectral radius of the matrix T . The goal the relaxation parameter so as to make the spectral radius of T as small As we will see, a clever choice of can result in a dramatic speed up in the rate. Let us look at an elementary example. 2 -1 , which we decompose as A = Example 7.31. Consider the matrix A = -1 2 L + D + U , where L= 0 -1 0 0 , D= 2 0 0 2 , U= 0 -1 0 0 . 0
1 2 1 2
Jacobi iteration is based on the coefficient matrix T = - D-1 (L + U ) =
. Its 0 spectral radius is (T ) = .5, and hence the Jacobi scheme takes, on average, roughly 3.3 - 1/ log10 .5 iterations to produce each new decimal place in the solution. The SOR scheme (7.63) takes the explicit form 2 0 - 2 2 - 0 2
-1
u(k+1) =
2(1 - ) 0
u(k) + b, 2 (1 - ) 1- 1 (1 - ) 2 . (2 - )2
1 2
where GaussSeidel is the particular case = 1. The SOR coefficient matrix is T = 2 (1 - ) 0 2 (1 - ) =
1 4
To compute the eigenvalues of T , we form its characteristic equation
1 1 0 = det(T - I ) = 2 - 2 - 2 + 4 2 + (1 - )2 = ( + - 1)2 - 4 2 . (7.67)
Our goal is to choose so that (a) both eigenvalues are less than 1 in modulus, so | 1 |, | 2 | < 1. This is the minimal requirement for convergence of the method. (b) the largest eigenvalue (in modulus) is as small as possible. This will give the smallest spectral radius for T and hence the fastest convergence rate. 10/18/06 128
c 2006 Peter J. Olver
The product of the two eigenvalues is the determinant, 1 2 = det T = (1 - )2 . If 0 or 2, then det T 1, and hence at least one of the eigenvalues would have modulus larger than 1. Thus, in order to ensure convergence, we must require 0 < < 2. For GaussSeidel, at = 1, the eigenvalues are 1 = 1 , 2 = 0, and the spectral radius is 4 (T1 ) = .25. This is exactly the square of the Jacobi spectral radius, and hence the Gauss Seidel iterates converge twice as fast; so it only takes, on average, about -1/ log10 .25 = 1.66 GaussSeidel iterations to produce each new decimal place of accuracy. It can be shown that as increases above 1, the two eigenvalues move along the real axis towards each other. They coincide when = = 8 - 4 3 1.07, at which point 1 = 2 = - 1 = .07 = (T ), which is the convergence rate of the optimal SOR scheme. Each iteration produces slightly more than one new decimal place in the solution, which represents a significant improvement over the GaussSeidel convergence rate. It takes about twice as many GaussSeidel iterations (and four times as many Jacobi iterations) to produce the same accuracy as this optimal SOR method. Of course, in such a simple 2 2 example, it is not so surprising that we can construct the best value for the relaxation parameter by hand. Young was able to find the optimal value of the relaxation parameter for a broad class of matrices that includes most of those arising in the finite difference and finite element numerical solutions to ordinary and partial differential equations. For the matrices in Young's class, the Jacobi eigenvalues occur in signed pairs. If are a pair of eigenvalues for the Jacobi method, then the corresponding eigenvalues of the SOR iteration matrix satisfy the quadratic equation ( + - 1)2 = 2 2 . (7.68)
If = 1, so we have standard GaussSeidel, then 2 = 2 , and so the eigenvalues are = 0, = 2 . The GaussSeidel spectral radius is therefore the square of the Jacobi spectral radius, and so (at least for matrices in the Young class) its iterates converge twice as fast. The quadratic equation (7.68) has the same properties as in the 2 2 version (7.67) (which corresponds to the case = 1 ), and hence the optimal value of will be the 2 one at which the two roots are equal: 1 = 2 = - 1, which occurs when = 2 2 - 2 1 - 2 = . 2 1 + 1 - 2
For example, if J = .99, which is rather slow convergence (but common for iterative numerical solution schemes for partial differential equations), then GS = .9801, which is 10/18/06 129
c 2006 Peter J. Olver
Therefore, if J = max | | denotes the spectral radius of the Jacobi method, then the GaussSeidel has spectral radius GS = 2 , while the SOR method with optimal relaxation J parameter 2 , has spectral radius = - 1. = (7.69) 1 + 1 - 2 J
twice as fast, but still quite slow, while SOR with = 1.7527 has = .7527, which is dramatically faster . Indeed, since (GS )14 (J )28 , it takes about 14 GaussSeidel (and 28 Jacobi) iterations to produce the same accuracy as one SOR step. It is amazing that such a simple idea can have such a dramatic effect.
More precisely, since the SOR matrix is not diagonalizable, the overall convergence rate is slightly slower than the spectral radius. However, this technical detail does not affect the overall conclusion.
10/18/06
130
c 2006
Peter J. Olver