師大

### Iterative Methods for Solving Large Linear Systems (I)

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University

October 25, 2011

師大

### Outline

1 Jacobi and Gauss-Seidel methods

2 Successive Over-Relaxation (SOR) Method

3 Symmetric Successive Over Relaxation

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 2 / 45

師大

General procedures for the construction of iterative methods

Given a linear system of nonsingular A

Ax = b, (1)

we consider the splitting of A

A = M − N (2)

with M nonsingular. Then (1) is equivalent to M x = N x + b, or
x = M^{−1}N x + M^{−1}b ≡ T x + f.

This suggests an iterative process

x_{k+1}= T x_{k}+ f = M^{−1}N x_{k}+ M^{−1}b, (3)
where x_{0} is given. Then the solution x of (1) is determined by iteration.

師大

Remark 1

(a) Define ε_{k}= x_{k}− x. Then

ε_{k+1}= x_{k+1}− x = M^{−1}N x_{k}+ M^{−1}b − M^{−1}N x − M^{−1}b

= (M^{−1}N )ε_{k} = (M^{−1}N )^{k}ε_{0}

which implies that ρ(M^{−1}N ) < 1 if and only if {ε_{k}} → 0.

(b) Let rk= b − Axk. Then,

xk+1 = M^{−1}N xk+ M^{−1}b

= M^{−1}(M − A)x_{k}+ M^{−1}b

= x_{k}+ M^{−1}(b − Ax_{k})

= x_{k}+ z_{k}
where M zk= rk.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 4 / 45

師大

Example 1

We consider the standard splitting of A

A = D − L − R, (4)

where A = [a_{ij}]^{n}_{i,j=1}, D = diag(a_{11}, a_{22}, · · · , a_{nn}),

−L =

0 0

a21 0

... . .. . ..

an1 · · · an,n−1 0

,

−R =

0 a_{12} · · · a_{1n}
0 . .. ...

. .. a_{n−1,n}

0 0

.

師大

For a_{i,i}6= 0, i = 1, . . . , n, D is nonsingular. If we choose
M = D and N = L + R

in (2), we then obtain the Jacobi Method (Total-step Method):

xk+1= D^{−1}(L + R)xk+ D^{−1}b
or in formula

x_{k+1,j} = 1

a_{jj}(−X

i6=j

a_{ji}x_{k,i}+ b_{j}), j = 1, . . . , n, k = 0, 1, . . . .

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 6 / 45

師大

Example 2

If D − L is nonsingular in (4), then we choose M = D − L, N = R

as in (2) are possible and yields the so-called Gauss-Seidel Method (Single-Step Method):

xk+1= (D − L)^{−1}Rxk+ (D − L)^{−1}b
or in formula

x_{k+1,j} = 1

a_{jj}(−X

i<j

ajix_{k+1,i}−X

i>j

ajix_{k,i}+bj), j = 1, . . . , n, k = 1, 2, . . . .

- Total-Step Method = TSM = Jacobi method.

- Single-Step Method = SSM = Gauss-Seidel method.

師大

We consider the following points on Examples 1 and 2:

(i) flops counts per iteration step.

(ii) Convergence speed.

Let k · k be a vector norm, and kT k be the corresponding operator norm.

Then kε_{m}k

kε_{0}k = kT^{m}ε_{0}k

kε_{0}k ≤ kT^{m}k. (5)

Here kT^{m}k^{m}^{1} is a measure for the average of reduction of error εm per
iteration step. We call

Rm(T ) = − ln(kT^{m}k^{m}^{1}) = −1

mln(kT^{m}k) (6)

the average of convergence rate for m iterations.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 8 / 45

師大

The larger is R_{m}(T ), so the better is convergence rate. Let
σ = (kε_{m}k/kε_{0}k)^{m}^{1}. From (5) and (6) we get

σ ≤ kT^{m}k^{m}^{1} ≤ e^{−R}^{m}^{(T )},
or

σ^{1/R}^{m}^{(T )}≤ 1
e.

That is, after 1/R_{m}(T ) steps in average the error is reduced by a factor of
1/e. Since Rm(T ) is not easy to determine, we consider m → ∞. Since

m→∞lim kT^{m}k^{m}^{1} = ρ(T ),
it follows

R∞(T ) = lim

m→∞R_{m}(T ) = − ln ρ(T ).

R∞ is called the asymptotic convergence rate. It holds always
R_{m}(T ) ≤ R∞(T ).

師大

Example 3

Consider the Dirichlet boundary-value problem (Model problem):

−∆u ≡ −u_{xx}− u_{yy}= f (x, y), 0 < x, y < 1, (7)
u(x, y) = 0 (x, y) ∈ ∂Ω,

for the unit square Ω := {x, y|0 < x, y < 1} ⊆ R^{2} with boundary ∂Ω.

To solve (7) by means of a difference methods, one replaces the differential operator by a difference operator. Let

Ω_{h} := {(x_{i}, yi)|i, j = 1, . . . , N + 1},

∂Ωh := {(x_{i}, 0), (xi, 1), (0, yj), (1, yj)|i, j = 0, 1, . . . , N + 1},
where

xi= ih, yj = jh, i, j = 0, 1, . . . , N + 1, h := _{N +1}^{1} , N ≥ 1, an integer.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 10 / 45

師大

The differential operator −u_{xx}− u_{yy} can be replaced for all (x_{i}, y_{i}) ∈ Ω_{h}
by the difference operator:

4ui,j− u_{i−1,j} − u_{i+1,j} − u_{i,j−1}− u_{i,j+1}

h^{2} (8)

up to an error τi,j. For small h one can expect that the solution zi,j, for i, j = 1, . . . , N , of the linear system

4zi,j− z_{i−1,j}− z_{i+1,j} − z_{i,j−1}− z_{i,j+1}= h^{2}fi,j, i, j = 1, . . . , N, (9)
z0,j = zN +1,j = zi,0= zi,N +1 = 0, i, j = 0, 1, . . . , N + 1,

obtained from (8) by omitting the error τi,j, agrees approximately with the
u_{i,j}. Let

z = [z_{1,1}, z_{2,1}, · · · , z_{N,1}, z_{1,2}, · · · , z_{N,2}, · · · , z_{1,N}, · · · , z_{N,N}]^{T}
and

b = h^{2}[f1,1, · · · , f_{N,1}, f1,2, · · · , f_{N,2}, · · · , f_{1,N}, · · · , f_{N,N}]^{T}.

師大

Then (9) is equivalent to a linear system Az = b with the N^{2}× N^{2}matrix.

4 −1 −1

−1 ...

...

...

... −1

...

−1 4 −1

−1 4 −1

...

... −1

... ...

...

.. .

.. . −1

.. .

−1 −1 4

.. . ...

... −1

...

...

...

...

...

...

.. .

..

. −1

−1 4 −1

..

. −1

.. .

.. . ..

.

.. . −1

−1 −1 4

(11)

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 12 / 45

師大

Let A = D − L − R. The matrix J = D^{−1}(L + R) belongs to the Jacobi
method (TSM). The N^{2} eigenvalues and eigenvectors of J can be
determined explicitly. We can verify at once, by substitution, that N^{2}
vectors z^{(k,l)}, k, l = 1, . . . , N with components

z_{i,j}^{(k,l)} := sin kπi

N + 1sin lπj

N + 1, 1 ≤ i, j ≤ N, satisfy

J z^{(k,l)} = λ^{(k,l)}z^{(k,l)}
with

λ^{(k,l)} := 1

2(cos kπ

N + 1 + cos lπ

N + 1), 1 ≤ k, l ≤ N.

師大

J thus has eigenvalues λ^{(k,l)}, 1 ≤ k, l ≤ N . Then we have
ρ(J ) = λ_{1,1} = cos π

N + 1 = 1 −π^{2}h^{2}

2 + O(h^{4}) (12)

and

R∞(J ) = − ln(1 − π^{2}h^{2}

2 + O(h^{4})) = π^{2}h^{2}

2 + O(h^{4}).

These show that (i) TSM converges;

(ii) Diminution of h will not only enlarge the flop counts per step, but also the convergence speed will drastically make smaller.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 14 / 45

師大

### Some theorems and definitions

ρ(T ): A measure of quality for convergence.

Definition 4

A real m × n-matrix A = (a_{ik}) is called nonnegative (positive), denoted
by A ≥ 0 (A > 0), if a_{ik} ≥ 0 (> 0), i = 1, . . . , m, k = 1, . . . , n.

Definition 5

An m × n-matrix A is called reducible, if there is a subset

I ⊂ ˜N ≡ {1, 2, . . . , n}, I 6= φ, I 6= ˜N such that i ∈ I, j 6∈ I ⇒ a_{ij} = 0.

A is not reducible ⇔ A is irreducible.

師大

Remark 2

G(A) is the directed graph associated with the matrix A. If A is an
n × n-matrix, then G(A) consists of n vertices P1, · · · , Pn and there is an
(oriented) arc Pi → P_{j} in G(A) precisely if aij 6= 0.

It is easily shown that A is irreducible if and only if the graph G(A) is
connected in the sense that for each pair of vertices (Pi, Pj) in G(A) there
is an oriented path from Pi to Pj. i.e., if i 6= j, there is a sequence of
indices i = i_{1}, i_{2}, · · · , i_{s}= j such that (a_{i}_{1}_{,i}_{2} · · · a_{i}_{s−1}_{,i}_{s}) 6= 0.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 16 / 45

師大

Theorem 6 (Perron-Frobenius) Let A ≥ 0 irreducible. Then

(i) ρ = ρ(A) is a simple eigenvalue;

(ii) There is a positive eigenvector z associated to ρ, i.e., Az = ρz, z > 0;

(iii) If Ax = λx, x ≥ 0, then λ = ρ, x = αz, α > 0;

(iv) A ≤ B, A 6= B =⇒ ρ(A) < ρ(B).

師大

Theorem 7

Let A ≥ 0, x > 0. Define the quotients:

qi(x) ≡ (Ax)i

xi

= 1 xi

n

X

k=1

aikxk, for i = 1, . . . , n.

Then

1≤i≤nmin qi(x) ≤ ρ(A) ≤ max

1≤i≤nqi(x). (13)

If A is irreducible, then it holds additionally, either

q1 = q2= · · · = qn (then x = µz, qi= ρ(A)) or

1≤i≤nmin qi(x) < ρ(A) < max

1≤i≤nqi(x). (14)

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 18 / 45

師大

Theorem 8

The statements in Theorem 7 can be formulated as: Let A ≥ 0, x > 0.

(13) corresponds:

Ax ≤ µx ⇒ ρ ≤ µ, Ax ≥ νx ⇒ ν ≤ ρ.

(15)

Let A ≥ 0, irreducible, x > 0. (14) corresponds :

Ax ≤ µx, Ax 6= µx ⇒ ρ < µ,

Ax ≥ νx, Ax 6= νx ⇒ ν < ρ. (16)

Definition 9

A real matrix B is called an M -matrix if b_{ij} ≤ 0, i 6= j and B^{−1} exists
with B^{−1}≥ 0.

師大

### Sufficient conditions for convergence of TSM and SSM

Theorem 10

Let B be a real matrix with b_{ij} ≤ 0 for i 6= j. Then the following statements are
equivalent.

(i) B is an M −matrix.

(ii) There exists a vector v > 0 so that Bv > 0.

(iii) B has a decomposition B = sI − C with C ≥ 0 and ρ(C) < s.

(iv) For each decomposition B = D − C with D = diag (di) and C ≥ 0, it
holds: di> 0, i = 1, 2, . . . , n, and ρ(D^{−1}C) < 1.

(v) There is a decomposition B = D − C, with D = diag(di) and C ≥ 0 it
holds: di> 0, i = 1, 2, . . . , n and ρ(D^{−1}C) < 1.

Further, if B is irreducible, then (vi) is equivalent to (i)-(v).

(vi) There exists a vector v > 0 so that Bv ≥ 0, 6= 0.

Proof

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 20 / 45

師大

Lemma 11

Let A be an arbitrary complex matrix and define |A| = [|aij|]. If |A| ≤ C, then ρ(A) ≤ ρ(C). Especially ρ(A) ≤ ρ(|A|).

Proof

Theorem 12

Let A be an arbitrary complex matrix. It satisfies either (Strong Row Sum Criterion):

X

j6=i

|a_{ij}| < |a_{ii}|, i = 1, . . . , n. (17)

or (Weak Row Sum Criterion):

X

j6=i

|aij| ≤ |aii|, i = 1, . . . , n,

< |ai_{0}i_{0}|, at least one i0, (18)

for A irreducible. Then TSM(Jacobi) and SSM(GS) both are convergent.

Proof

師大

Relaxation Methods (Successive Over-Relaxation (SOR) Method) Consider the parametrized linear system ωAx = ωb and consider the splitting

ωA = ωD − ωL − ωR + D − D

= (D − ωL) − ((1 − ω)D + ωR) ≡ M − N.

From (3) we have the iteration

x_{k+1} = (D − ωL)^{−1}((1 − ω)D + ωR) x_{k}+ ω(D − ωL)^{−1}b. (19)
From Remark 1 (b) the iteration (19) is equivalent to

x_{k+1} = x_{k}+ ωz_{k}
where

(D − ωL)z_{k}= r_{k}≡ b − Ax_{k}.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 22 / 45

師大

Define

Lω := (D − ωL)^{−1}((1 − ω)D + ωR) .
We may assume D = I, i.e.,

Lω:= (I − ωL)^{−1}((1 − ω)I + ωR) .

Otherwise, we can let ˜A = D^{−1}A, ˜L = D^{−1}L, ˜R = D^{−1}R. Then it holds
that

A = I − ˜˜ L − ˜R.

ω < 1: under relaxation

ω = 1: single-step method (GS) ω > 1: over relaxation.

師大

We now try to choose an ω such that ρ(L_{ω}) is small as possible. But this
is only under some special assumptions possible. we first list a few

qualitative results about ρ(L_{ω}).

Theorem 13

Let A = D − L − L^{∗} be Hermitian and positive definite. Then the
relaxation method is convergent for 0 < ω < 2.

Theorem 14

Let A be Hermitian and nonsingular with positive diagonal. If SSM converges, then A is positive definite.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 24 / 45

師大

### Determination of the Optimal Parameter ω for 2-consistly Ordered Matrices

For an important class of matrices the more qualitative assertions of
Theorems 13 and 14 can be considerably sharpened. This is the class of
consistly ordered matrices. The optimal parameter ω_{b} with

ρ(L_{ω}_{b}) = min

ω ρ(L_{ω})
can be determined. We consider A = I − L − R.

Definition 15

A is called 2-consistly ordered, if the eigenvalues of αL + α^{−1}R are
independent of α.

師大

If A is 2-consistly ordered, then L + R and −(L + R) (α = −1) has the same eigenvalues. The nonzero eigenvalues of L + R appear in pairs.

Hence

det(λI − L − R) = λ^{m}

r

Y

i=1

(λ^{2}− µ^{2}_{i}), n = 2r + m (m = 0, possible). (20)

Theorem 16

Let A be 2-consistly ordered, a_{ii}= 1, ω 6= 0. Then hold:

(i) If λ 6= 0 is an eigenvalue of L_{ω} and µ satisfies the equation

(λ + ω − 1)^{2}= λµ^{2}ω^{2}, (21)
then µ is an eigenvalue of L + R (so is −µ).

(ii) If µ is an eigenvalue of L + R and λ satisfies the equation (21), then λ is an eigenvalue of Lω.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 26 / 45

師大

Remark 3

If ω = 1, then λ = µ^{2}, and ρ((I − L)^{−1}R) = (ρ(L + R))^{2}.
Proof: We first prove the identity

det(λI − sL − rR) = det(λI −√

sr(L + R)). (22)

Since both sides are polynomials of the form λ^{n}+ · · · and
sL + rR =√

sr(r s

rL +r r

sR) =√

sr(αL + α^{−1}R),
if sr 6= 0, then sL + rR and √

sr(L + R) have the same eigenvalues. It is obviously also for the case sr = 0. The both polynomials in (22) have the same roots, so they are identical.

師大

For

det(I − ωL) det(λI − L_{ω}) = det(λ(I − ωL) − (1 − ω)I − ωR)

= det((λ + ω − 1)I − ωλL − ωR) = Φ(λ) and det(I − ωL) 6= 0, λ is an eigenvalue of Lω if and only if Φ(λ) = 0.

From (22) follows

Φ(λ) = det((λ + ω − 1)I − ω

√

λ(L + R)) and that is (from (20))

Φ(λ) = (λ + ω − 1)^{m}

r

Y

i=1

((λ + ω − 1)^{2}− λµ^{2}_{i}ω^{2}), (23)

where µ_{i} is an eigenvalue of L + R. Therefore, if µ is an eigenvalue of
(L + R) and λ satisfies (21), so is Φ(λ) = 0, then λ is eigenvalue of Lω.
This shows (ii).

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 28 / 45

師大

Now if λ 6= 0 an eigenvalue of L_{ω}, then one factor in (23) must be zero.

Let µ satisfy (21). Then

(i) µ 6= 0: From (21) follows λ + ω − 1 6= 0, so

(λ + ω − 1)^{2} = λω^{2}µ^{2}_{i}, for one i (from (23)),

= λω^{2}µ^{2}, (from (21)).

This shows that µ = ±µi, so µ is an eigenvalue of L + R.

(ii) µ = 0: We have λ + ω − 1 = 0 and 0 = Φ(λ) = det((λ + ω − 1)I − ω√

λ(L + R)) = det(−ω√

λ(L + R)), i.e., L + R is singular, so µ = 0 is eigenvalue of L + R.

師大

Theorem 17

Let A = I − L − R be 2-consistly ordered. If L + R has only real eigenvalues and satisfies ρ(L + R) < 1, then it holds

ρ(L_{ω}_{b}) = ω_{b}− 1 < ρ(L_{ω}), for ω 6= ω_{b},
where

ω_{b} = 2

1 +p1 − ρ^{2}(L + R) (solve ω_{b} in (21)).

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

0.9975 0.998 0.9985 0.999 0.9995 1

ω

** spectral radius**

Figure: figure of ρ(Lω_{b})

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 30 / 45

師大

Theorem 18 One has in general,

ρ(Lω) =

( ω − 1, for ωb≤ ω ≤ 2
1 − ω + ^{1}_{2}ω^{2}µ^{2}+ ωµ

q

1 − ω + ^{1}_{4}ω^{2}µ^{2}, for 0 < ω ≤ ω_{b}
(24)

師大

Remark: We first prove the following Theorem proposed by Kahan: For arbitrary matrices A it holds

ρ(Lω) ≥ |ω − 1|, for all ω. (25)
Since det(I − ωL) = 1 for all ω, the characteristic polynomial Φ(λ) of L_{ω}
is

Φ(λ) = det(λI − L_{ω}) = det((I − ωL)(λI − L_{ω}))

= det((λ + ω − 1)I − ωλL − ωR).

For

n

Q

i=1

λi(Lω) = Φ(0) = det((ω − 1)I − ωR) = (ω − 1)^{n}, it follows
immediately that

ρ(Lω) = max

i |λ_{i}(Lω)| ≥ |ω − 1|.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 32 / 45

師大

Proof of Theorem: By assumption the eigenvalues µ_{i} of L + R are real
and −ρ(L + R) ≤ µi≤ ρ(L + R) < 1. For a fixed ω ∈ (0, 2) (by (25) in
the Remark it suffices to consider the interval (0,2)) and for each µ_{i} there
are two eigenvalues λ^{(1)}_{i} (ω, µ_{i}) and λ^{(2)}_{i} (ω, µ_{i}) of L_{ω}, which are obtained
by solving the quadratic equation (21) in λ.

Geometrically, λ^{(1)}_{i} (ω) and λ^{(2)}_{i} (ω) are obtained as abscissae of the points
of intersection of

the straight line g_{ω}(λ) = λ + ω − 1
ω
and

the parabola m_{i}(λ) := ±√
λµ_{i}

(see Figure 2). The line g_{ω}(λ) has the slope 1/ω and passes through the
point (1,1). If gω(λ) ∩ mi(λ) = φ, then λ^{(1)}_{i} (ω) and λ^{(2)}_{i} (ω) are conjugate
complex with modulus |ω − 1| (from (21)).

師大

Evidently

ρ(L_{ω}) = max

i (|λ^{(1)}_{i} (ω)|, |λ^{(2)}_{i} (ω)|) = max(|λ^{(1)}(ω)|, |λ^{(2)}(ω)|),
where λ^{(1)}(ω), λ^{(2)}(ω) being obtained by intersecting g_{ω}(λ) with
m(λ) := ±√

λµ, with µ = ρ(L + R) = maxi|µ_{i}|. By solving (21) with
µ = ρ(L + R) for λ, one verifies (24) immediately, and thus also the
remaining assertions of the theorem.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 34 / 45

師大

1 1

( )

gω λ ( )

b gω λ

( )

m λ

( ) i

m λ

λ

(1)

i

λ

( 2 )

i

λ θ

Figure: Geometrical view of λ^{(1)}_{i} (ω) and λ^{(2)}_{i} (ω).

師大

### Application to Finite Difference Methods: Model Problem

We consider the Dirichlet boundary-value problem (Model problem) as
in Example 3. We shall solve a linear system Az = b of the N^{2}× N^{2}
matrix A as in (11).

To Jacobi method: The iterative matrix is J = L + R = 1

4(4I − A).

It is easily seen that A is 2-consistly ordered (Exercise!).

To Gauss-Seidel method: The iterative matrix is
H = (I − L)^{−1}R.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 36 / 45

師大

From the Remark of Theorem 16 and (12) follows that
ρ(H) = ρ(J )^{2} = cos^{2} π

N + 1.

According to Theorem 17 the optimal relaxation parameter ω_{b} and ρ(L_{ω}_{b})
are given by

ω_{b} = 2

1 +q

1 − cos^{2} _{N +1}^{π}

= 2

1 + sin_{N +1}^{π}
and

ρ(Lω_{b}) = cos^{2}_{N +1}^{π}
(1 + sin_{N +1}^{π} )^{2}.

The number k = k(N ) with ρ(J )^{k}= ρ(L_{ω}_{b}) indicates that the k steps of
Jacobi method produce the same reduction as one step of the optimal
relaxation method. Clearly,

k = ln ρ(Lωb)/ ln ρ(J ).

師大

Now for small z one has ln(1 + z) = z − z^{2}/2 + O(z^{3}) and for large N we
have

cos

π

N + 1

= 1 − π^{2}

2(N + 1)^{2} + O( 1
N^{4}).

Thus that

ln ρ(J ) = π^{2}

2(N + 1)^{2} + O( 1
N^{4}).

Similarly,

ln ρ(L_{ω}_{b}) = 2[ln ρ(J ) − ln(1 + sin π
N + 1)]

= 2[− π^{2}

2(N + 1)^{2} − π

N + 1 + π^{2}

2(N + 1)^{2} + O( 1
N^{3})]

= − 2π

N + 1+ O( 1

N^{3}) (for large N ).

and

k = k(N ) ≈ 4(N + 1)

π .

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 38 / 45

師大

The optimal relaxation method is more than N times as fast as the Jacobi method. The quantities

R_{J} := − ln 10

ln ρ(J ) ≈ 0.467(N + 1)^{2}. (26)
R_{H} := 1

2R_{J} ≈ 0.234(N + 1)^{2} (27)

RL_{ωb} := − ln 10

ln ρ(L_{ω}_{b}) ≈ 0.367(N + 1) (28)
indicate the number of iterations required in the Jacobi, the Gauss-Seidel
method, and the optimal relaxation method, respectively, in order to
reduce the error by a factor of 1/10.

師大

SSOR (Symmetric Successive Over Relaxation):

A is symmetric and A = D − L − L^{T}. Let

M_{ω}: = D − ωL,

N_{ω}: = (1 − ω)D + ωL^{T}, and

M_{ω}^{T} = D − ωL^{T},
N_{ω}^{T} = (1 − ω)D + ωL.

Then from the iterations

Mωx_{i+1/2} = Nωxi+ ωb,
M_{ω}^{T}xi+1 = N_{ω}^{T}x_{i+1/2}+ ωb,
follows that

x_{i+1} = M_{ω}^{−T}N_{ω}^{T}M_{ω}^{−1}N_{ω} x_{i}+ ˜b

≡ Gx_{i}+ ω M_{ω}^{−T}N_{ω}^{T}M_{ω}^{−1}+ M_{ω}^{−T} b

≡ Gx_{i}+ M (ω)^{−1}b.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 40 / 45

師大

It holds that

((1 − ω)D + ωL) (D − ωL)^{−1}+ I

= (ωL − D − ωD + 2D)(D − ωL)^{−1}+ I

= −I + (2 − ω)D(D − ωL)^{−1}+ I

= (2 − ω)D(D − ωL)^{−1},
Thus

M (ω)^{−1} = ω D − ωL^{T}−1

(2 − ω)D(D − ωL)^{−1},
then

M (ω) = 1

ω(2 − ω)(D − ωL)D^{−1} D − ωL^{T}

(29)

≈ (D − L)D^{−1} D − L^{T} , (ω = 1).

師大

### Appendix

Proof of Theorem

(1) =⇒ (2): Let e = (1, · · · , 1)^{T}. Since B^{−1}≥ 0 is nonsingular it follows
ν = B^{−1}e > 0 and Bν = B(B^{−1}e) = e > 0.

(2) =⇒ (3): Let s > max(bii). It follows B = sI − C with C ≥ 0. There exists a ν > 0 with Bν = sν − Cν (via (2)), also sν > Cν. From the statement (15) in Theorem 8 follows ρ(C) < s.

(3) =⇒ (1): B = sI − C = s(I −^{1}_{s}C). For ρ(^{1}_{s}C) < 1 and from Theorem
2.6 (I −^{1}_{s}C)^{−1} follows that there exists a series expansion

∞

P

ν=0

(^{1}_{s}C)^{k}.
Since the terms in sum are nonnegtive, we get B^{−1} = ^{1}_{s}(I − ^{1}_{s}C)^{−1} ≥ 0.

(2) =⇒ (4): From Bν = Dν − Cν > 0 follows Dν > Cν ≥ 0 and d_{i} > 0,
for i = 1, 2, · · · , n. Hence D^{−1} ≥ 0 and ν > D^{−1}Cν ≥ 0. From (15)
follows that ρ(D^{−1}C) < 1.

(4) =⇒ (5): Trivial.

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 42 / 45

師大

(5) =⇒ (1): Since ρ(D^{−1}C) < 1, it follows from Theorem 2.6 that
(I − D^{−1}C)^{−1} exists and equals to

∞

P

k=0

(D^{−1}C)^{k}. Since the terms in sum
are nonnegative, we have (I − D^{−1}C)^{−1} is nonnegative and

B^{−1}= (I − D^{−1}C)^{−1}D^{−1} ≥ 0.

(2) =⇒ (6): Trivial.

(6) =⇒ (5): Consider the decomposition B = D − C, with d_{i} = b_{ii}. Let
{I = i | d_{i}≤ 0}. From d_{i}νi−P

k6=icikνk≥ 0 follows c_{ik} = 0 for i ∈ I,
and k 6= i. Since Bν ≥ 0, 6= 0 =⇒ I 6= {1, · · · , n}. But B is irreducible

=⇒ I = ∅ and d_{i} > 0. Hence for Dν >, 6= Cν also ν >, 6= D^{−1}Cν and
(16) show that ρ(D^{−1}C) < 1.

return

師大

Proof of Lemma 11 There is a x 6= 0 with Ax = λx and |λ| = ρ(A).

Hence

ρ(A)|xi| = |

n

X

k=1

aikxk| ≤

n

X

k=1

|a_{ik}||x_{k}| ≤

n

X

k=1

cik|x_{k}|.

Thus,

ρ(A)|x| ≤ C|x|.

If |x| > 0, then from (15) we have ρ(A) ≤ ρ(C). Otherwise, let

I = {i | x_{i} 6= 0} and C_{I} be the matrix, which consists of the ith row and
ith column of C with i ∈ I. Then we have ρ(A)|xI| ≤ C_{I}|x_{I}|. Here |x_{I}|
consists of ith component of |x| with i ∈ I. Then from |x_{I}| > 0 and (15)
follows ρ(A) ≤ ρ(C_{I}). We now fill C_{I} with zero up to an n × n matrix ˜C_{I}.
Then ˜CI ≤ C. Thus, ρ(C_{I}) = ρ( ˜CI) ≤ ρ(C) (by Theorem ??(3)).

return

T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 44 / 45

師大

Proof of Theorem ?? Let A = D − L − R. From (17) and (18) D must be nonnsingular and then as in Remark 9.3 we can w.l.o.g. assume that D = I. Now let B = I − |L| − |R|. Then (17) can be written as Be > 0.

From Theorem 10(2) and (1) follows that B is an M -matrix.

(18) can be written as Be ≥ 0, Be 6= 0. Since A is irreducible, also B, from Theorem 10 (6) and (1) follows that B is an M -matrix.

Especially, from theorem 10(1), (4) and Theorem ?? follows that
ρ(|L| + |R|) < 1 and ρ((I − |L|)^{−1}|R|) < 1.

Now Lemma 11 shows that

ρ(L + R) ≤ ρ(|L| + |R|) < 1.

So TSM is convergent. Similarly,

ρ((I − L)^{−1}R) = ρ(R + LR + · · · + L^{n−1}R)

≤ ρ(|R| + |L||R| + · · · + |L|^{n−1}|R|)

= ρ((I − |L|)^{−1}|R|) < 1.

So SSM is convergent.

return