Successive Over Relaxation Method

THE SUCCESSIVE OVER RELAXATION METHOD Math. 609, Fall 2008 (updated: September 16, 2008) In this class, we consider the “Successive Over Relaxation Method” (SOR). This is another example of a splitting method. In this case we split A = (ω −1 D + L) + ((1 − ω −1 )D + U ), i.e., (ω −1 D + L)x = ((ω −1 − 1)D − U )x + b and obtain the iterative method (5.1)

(ω −1 D + L)xi+1 = ((ω −1 − 1)D − U )xi + b.

As in the Gauss-Seidel and Jacobi iterations, we can only apply this method to matrices with non-vanishing diagonal entries. In the above method, ω is an iteration parameter which we shall have to choose. When ω = 1, the method reduces to Gauss-Seidel. Like GaussSeidel, this method can be implemented as a sweep. We first develop a necessary condition on ω required for convergence of the SOR method. A simple computation (do it!) shows that the reduction matrix for the SOR method is GSOR = (ω −1 D + L)−1 ((ω −1 − 1)D − U ). From the corollaries of the previous class, a necessary condition for SOR to converge for all starting iterates and right hand sides is that ρ(GSOR ) < 1. The eigenvalues of GSOR are the roots of the characteristic polynomial, P (λ) = det(λI − GSOR ) = λn + Cn−1 λn−1 + · · · + C0 . Note that C0 = P (0) = (−1)n det(GSOR ) and that n Y P (λ) = (λ − λi ) i=1

where λ1 , λ2 , . . . , λi are the roots of P . Expanding the above expression shows that n Y C0 = (−1)n λi = (−1)n det(GSOR ). i=1

Note that if |C0 | ≥ 1, there has to be at least one root whose absolute value is at least one, i.e., ρ(GSOR ) ≥ 1. Consequently, a necessary condition for the convergence of the SOR iteration is that |C0 | = |det(GSOR )| < 1. 1

2

Now, the determinant of a triangular matrix is the product of the entries on the diagonal. Using this and other simple properties of the determinant gives (why?) Qn (ω −1 − 1)Di,i det((ω −1 − 1)D − U ) i=1 Q det(GSOR ) = = = (1 − ω)n . n −1 D det(ω −1 D + L) ω i,i i=1 For this product to be less that one in absolute value, it is necessary that |1 − ω| < 1, i.e., 0 < ω < 2. We restate this as a proposition. Proposition 1. A necessary condition for the SOR method to converge for any starting vector and right hand side is that 0 < ω < 2. We shall provide a convergence theorem for the SOR method. Note that the reduction matrix GSOR is generally not symmetric even when A is symmetric. Accordingly, even if A has real entries, any analysis of GSOR will have to involve complex numbers as, in general, GSOR will have complex eigenvalues. Accordingly, we shall provide an analysis for complex Hermitian matrices A. Definition 1. The conjugate transpose N ∗ of a general n × m matrix N with complex entries is the m × n matrix with entries ¯j,i . (N ∗ )i,j = N Here the bar denotes complex conjugate. An n × n matrix A with complex entries is called Hermitian (or conjugate symmetric) if A∗ = A. We shall use the Hermitian inner product (·, ·) on Cn × Cn defined by (x, y) =

n X

xi y¯i .

i=1

We shall discuss more general Hermitian inner products in later classes. We note that (x, y) = (y, x),

(αx, y) = α(x, y), and (x, αy) = α(x, ¯ y)

for all x, y ∈ Cn and α ∈ C. In addition, it is easy to see that when N is an n × m matrix, N ∗ is the unique matrix satisfying (check it by applying it to the standard basis vectors!) (N x, y) = (x, N ∗ y) for all x ∈ Cm , y ∈ Cn . In the above equality, (·, ·) denotes the Hermitian inner product on Cn on the left and the Hermitian inner product on Cm on the right. We also note the following properties of a Hermitian n × n matrix A.

3

• (Ax, x) is real since (Ax, x) = (x, Ax) = (Ax, x). • If A is positive definite1 then Ai,i > 0 since Ai,i = (Aei , ei ) > 0 where ei denotes the i’th standard basis vector for Cn . The first task of this class is to prove a convergence theorem for successive over relaxation (SOR). This proof is not very intuitive and was probably discovered by playing with the equations until something worked. Theorem 1. Let A be an n × n Hermitian positive definite matrix and ω be in (0, 2). Then the SOR method for iteratively solving Ax = b converges for any starting vector and right hand side. Proof. The reduction matrix for SOR is GSOR = (ω −1 D+L)−1 ((ω −1 −1)D− U ) and we set Q = (ω −1 D + L). We shall show that ρ(GSOR ) < 1. We note that (5.2) (I −GSOR ) = (ω −1 D+L)−1 ((ω −1 D+L)−((ω −1 −1)D−U )) = Q−1 A. Let x ∈ Cn be an eigenvector of GSOR with eigenvalue λ and set y = (I − GSOR )x = (1 − λ)x. Then, by (??), (5.3)

Qy = Ax

and (5.4)

(Q − A)y = (Q − A)Q−1 Ax = (A − AQ−1 A)x = A(I − Q−1 A)x = AGSOR x = λAx.

Taking the inner product of (??) with y in the second place and (??) with y in the first place gives (Qy, y) = (Ax, y), and (y, (Q − A)y) = (y, λAx). These equations can be rewritten (using the fact that D is real and symmetric) ¯ ω −1 (Dy, y) + (Ly, y) = (1 − λ)(Ax, x), and −1 ¯ (ω − 1)(Dy, y) − (y, U y) = (1 − λ)λ(x, Ax). Now since A is Hermitian, (Ly, y) = (y, U y) so adding the two equations above gives (5.5)

(2ω −1 − 1)(Dy, y) = (1 − |λ|2 )(Ax, x).

1A complex n × n matrix A is positive definite if (Ax, x) > 0 for all non zero x ∈ Cn .

4

Now, x 6= 0 implies (Ax, x) > 0 so Ax 6= 0. As Q is nonsingular, (??) implies y 6= 0. In addition, the assumption on ω implies that 2ω −1 − 1 > 0 so that the left hand side of (??) is positive (from above we know that D is diagonal with positive numbers on the diagonal). For the right hand side of (??) to be positive, it is necessary that |λ| < 1, i.e., ρ(GSOR ) < 1.  Remark 1. The above proof works for other splittings of A, e.g. A = D + L + U where D is Hermitian positive definite and L∗ = U (D, L and U need not be diagonal, strictly lower triangular, and strictly upper triangular, respectively). Remark 2. Unlike the analysis for Gauss-Seidel and Jacobi given earlier, the above proof for SOR does not lead to any explicit bound for the spectral radius. All it shows is that the spectral radius has to be less than one. One step of an iterative method for solving Ax = b can be used to provide an “approximate” inverse. Specifically, we define B : Rn → Rn by defining By to be the result of one iterative step applied to solving Ax = y with initial iterate x0 = 0. All of the iterative methods that we have considered up to now have been of the form Qxi+1 = (Q − A)xi + b. With this notation, the reduction matrix G is given by G = I − Q−1 A and the corresponding approximate inverse is B = Q−1 . When A is symmetric positive definite (SPD) real (or Hermitian positive definite) it is advantageous to have approximate inverses which are of the same type. We shall consider the case when A is SPD real. The approximate inverse corresponding to the Jacobi method is B = D−1 and is SPD when A is SPD. The approximate inverse in the case of Gauss-Seidel is B = (D+L)−1 and is generally not symmetric. We can develop a symmetric approximate inverse from Gauss-Seidel by introducing the “transpose” iteration, (D + U )xi+1 = −Lxi + b. Note that the implementation of this method as a sweep is similar to the original Gauss-Seidel method except that one goes through the vector in reverse order. The pseudo-code is as follows: F U N CT ION gst(X, B, A, n) , 1 DO { F OR j = n, n − 1, . . .P Xj = (Bj − k6=j, Aj,k 6=0 Aj,k Xk )A−1 j,j .} RET U RN EN D

5

We get a symmetric approximate inverse by defining By by first applying one step of Gauss-Seidel with zero initial iterate and right hand side y and using the result x1 as an initial iterate for one step of the transpose iteration, i.e., By = x2 where (5.6)

(D + L)x1 = y, and (D + U )x2 = −Lx1 + y.

We can explicitly compute B = Q−1 by computing the reduction matrix G = I − Q−1 A for the two step procedure and identifying Q−1 . Let x be the solution of Ax = y and, as usual, set ei = x − xi with x0 = 0. The errors e1 and e0 are related by e1 = −(D + L)−1 U e0 = (I − (D + L)−1 A)e0 (the usual Gauss-Seidel reduction matrix) while the errors e2 and e1 are related by e2 = −(D + U )−1 Le1 = (I − (D + U )−1 A)e1 . Thus the error reduction matrix for the two steps is given by G = (I − (D + U )−1 A)(I − (D + L)−1 A) = I − [(D + U )−1 + (D + L)−1 − (D + U )−1 A(D + L)−1 ]A = I − (D + U )−1 [(D + L) + (D + U ) − A](D + L)−1 A = I − (D + U )−1 D(D + L)−1 A. Thus, the approximate inverse associated with the two steps is given by (5.7)

B = (D + U )−1 D(D + L)−1

and is obviously SPD when A is SPD (why?). Exercise 1. The symmetric successive over relaxation method (SSOR) is defined in an analogous fashion, i.e., by taking one step of (??) with zero initial iterate followed by onestep of the transpose iteration (ω −1 D + U )xi+1 = ((ω −1 − 1)D − L)xi + b. Compute the approximate inverse associated with SSOR. Your result should reduce to (??) when ω = 1.