**logo**

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

September 12, 2015

**logo**

**Outline**

**1** **Norms of vectors and matrices**

**2** **Eigenvalues and eigenvectors**

**3** **The Jacobi and Gauss-Siedel Iterative Techniques**

**4** **Relaxation Techniques for Solving Linear Systems**

**5** **Error bounds and iterative refinement**

**6** **The conjugate gradient method**

**2 / 99**

**logo**

**Definition 1**

k · k : R^{n}→ R is a vector norm if
**(i)** kxk ≥ 0, ∀ x ∈ R^{n},

**(ii)** kxk = 0 if and only if x = 0,

**(iii)** kαx| = |α|kxk ∀ α ∈ R and x ∈ R^{n},
**(iv)** kx + yk ≤ kxk + kyk ∀ x, y ∈ R^{n}.
**Definition 2**

The `_{2} and `∞norms for x = [x_{1}, x_{2}, · · · , x_{n}]^{T} are defined by

kxk_{2} = (x^{T}x)^{1/2}=
( _{n}

X

i=1

x^{2}_{i}
)1/2

and kxk_{∞}= max

1≤i≤n|x_{i}|.

The `_{2} norm is also called the Euclidean norm.

**logo**

**Theorem 3 (Cauchy-Bunyakovsky-Schwarz inequality)**
For each x = [x1, x2, · · · , xn]^{T} and y = [y1, y2, · · · , yn]^{T} in R^{n},

x^{T}y =

n

X

i=1

xiyi ≤
( _{n}

X

i=1

x^{2}_{i}

)1/2( _{n}
X

i=1

y^{2}_{i}
)1/2

= kxk2· kyk_{2}.

Proof: If x = 0 or y = 0, the result is immediate.

Suppose x 6= 0 and y 6= 0. For each α ∈ R,
0 ≤ kx − αyk^{2}_{2} =

n

X

i=1

(x_{i}− αy_{i})^{2}=

n

X

i=1

x^{2}_{i} − 2α

n

X

i=1

x_{i}y_{i}+ α^{2}

n

X

i=1

y_{i}^{2},

and 2α

n

X

i=1

x_{i}y_{i}≤

n

X

i=1

x^{2}_{i} + α^{2}

n

X

i=1

y_{i}^{2}= kxk^{2}_{2}+ α^{2}kyk^{2}_{2}.

**4 / 99**

**logo**

Since kxk_{2}> 0and kyk_{2} > 0, we can let
α = kxk_{2}

kyk_{2}
to give

2kxk_{2}

kyk_{2}

n

X

i=1

x_{i}y_{i}

!

≤ kxk^{2}_{2}+kxk^{2}_{2}

kyk^{2}_{2}kyk^{2}_{2}= 2kxk^{2}_{2}.
Thus

x^{T}y =

n

X

i=1

xiyi ≤ kxk_{2}kyk_{2}.

**logo**

For each x, y ∈ R^{n},
kx + yk_{∞} = max

1≤i≤n|x_{i}+ y_{i}| ≤ max

1≤i≤n(|x_{i}| + |y_{i}|)

≤ max

1≤i≤n|x_{i}| + max

1≤i≤n|y_{i}| = kxk_{∞}+ kyk∞

and

kx + yk^{2}_{2} =

n

X

i=1

(xi+ yi)^{2}=

2

X

i=1

x^{2}_{i} + 2

n

X

i=1

xiyi+

n

X

i=1

y_{i}^{2}

≤ kxk^{2}_{2}+ 2kxk2kyk_{2}+ kyk^{2}_{2} = (kxk2+ kyk2)^{2},
which gives

kx + yk_{2} ≤ kxk_{2}+ kyk_{2}.

**6 / 99**

**logo**

**Definition 4**

A sequence {x^{(k)}∈ R^{n}}^{∞}_{k=1}is convergent to x with respect to
the norm k · k if ∀ ε > 0, ∃ an integer N (ε) such that

kx^{(k)}− xk < ε, ∀ k ≥ N (ε).

**Theorem 5**

{x^{(k)}∈ R^{n}}^{∞}_{k=1} converges to x with respect to k · k∞if and only
if

k→∞lim x^{(k)}_{i} = x_{i}, ∀ i = 1, 2, . . . , n.

Proof: “⇒” Given any ε > 0, ∃ an integer N (ε) such that

1≤i≤nmax |x^{(k)}_{i} − x_{i}| = kx^{(k)}− xk_{∞}< ε, ∀ k ≥ N (ε).

**logo**

This result implies that

|x^{(k)}_{i} − x_{i}| < ε, ∀ i = 1, 2, . . . , n.

Hence

k→∞lim x^{(k)}_{i} = xi, ∀ i.

“⇐” For a given ε > 0, let Ni(ε)represent an integer with

|x^{(k)}_{i} − x_{i}| < ε, whenever k ≥ Ni(ε).

Define

N (ε) = max

1≤i≤nN_{i}(ε).

If k ≥ N (ε), then

1≤i≤nmax |x^{(k)}_{i} − x_{i}| = kx^{(k)}− xk_{∞}< ε.

This implies that {x^{(k)}} converges to x with respect to k · k∞.

**8 / 99**

**logo**

**Theorem 6**
For each x ∈ R^{n},

kxk_{∞}≤ kxk_{2} ≤√

nkxk∞.
Proof: Let x_{j} be a coordinate of x such that

kxk^{2}_{∞}= |xj|^{2}≤

n

X

i=1

x^{2}_{i} = kxk^{2}_{2},

so kxk∞≤ kxk_{2} and
kxk^{2}_{2} =

n

X

i=1

x^{2}_{i} ≤

n

X

i=1

x^{2}_{j} = nx^{2}_{j} = nkxk^{2}_{∞},

so kxk_{2} ≤√

nkxk∞.

**logo**

**Definition 7**

A matrix norm k · k on the set of all n × n matrices is a real-valued function satisfying for all n × n matrices A and B and all real number α:

**(i)** kAk ≥ 0;

**(ii)** kAk = 0 if and only if A = 0;

**(iii)** kαAk = |α|kAk;

**(iv)** kA + Bk ≤ kAk + kBk;

**(v)** kABk ≤ kAkkBk;

**Theorem 8**

If k · k is a vector norm on R^{n}, then
kAk = max

kxk=1kAxk

is a matrix norm. _{10 / 99}

**logo**

For any z 6= 0, we have x = z/kzk as a unit vector. Hence kAk = max

kxk=1kAxk = max

z6=0

A

z kzk

= max

z6=0

kAzk kzk .

**Corollary 9**

kAzk ≤ kAk · kzk.

**Theorem 10**

If A = [a_{ij}]is an n × n matrix, then
kAk_{∞}= max

1≤i≤n n

X

j=1

|a_{ij}|.

**logo**

Proof: Let x be an n-dimension vector with 1 = kxk∞= max

1≤i≤n|x_{i}|.

Then

kAxk_{∞} = max

1≤i≤n

n

X

j=1

aijxj

≤ max

1≤i≤n n

X

j=1

|a_{ij}| max

1≤j≤n|x_{j}| = max

1≤i≤n n

X

j=1

|a_{ij}|.

Consequently,

kAk_{∞}= max

kxk∞=1kAxk_{∞}≤ max

1≤i≤n n

X

j=1

|a_{ij}|.

On the other hand, let p be an integer with

n

X

j=1

|a_{pj}| = max

1≤i≤n n

X

j=1

|a_{ij}|,

**12 / 99**

**logo**

and x be the vector with xj =

1, if a_{pj}≥ 0,

−1, if apj< 0.

Then

kxk∞= 1 and a_{pj}xj = |apj|, ∀ j = 1, 2, . . . , n,
so

kAxk_{∞}= max

1≤i≤n

n

X

j=1

a_{ij}x_{j}

≥

n

X

j=1

a_{pj}x_{j}

=

n

X

j=1

|a_{pj}|

= max

1≤i≤n n

X

j=1

|a_{ij}|.

This result implies that
kAk_{∞}= max

kxk∞=1

kAxk_{∞}≥ max

1≤i≤n n

X

j=1

|a_{ij}|.

which gives

kAk_{∞}= max

1≤i≤n n

X|a_{ij}|. _{13 / 99}

**logo**

**Exercise**

Page 441: 5, 9, 10, 11

**14 / 99**

**logo**

**Eigenvalues and eigenvectors**

**Definition 11 (Characteristic polynomial)**

If A is a square matrix, the characteristic polynomial of A is defined by p(λ) =det(A − λI).

**Definition 12 (Eigenvalue and eigenvector)**

If p is the characteristic polynomial of the matrix A, the zeros of p are eigenvalues of the matrix A. If λ is an eigenvalue of A and x 6= 0 satisfies (A − λI)x = 0, then x is an eigenvector of A corresponding to the eigenvalue λ.

**Definition 13 (Spectrum and Spectral Radius)**

The set of all eigenvalues of a matrix A is called the spectrum of A.

The spectral radius of A is

ρ(A) = max{|λ|; λis an eigenvalue ofA}.

**logo**

**Theorem 14**

If A is an n × n matrix, then
**(i)** kAk_{2} =p

ρ(A^{T}A);

**(ii)** ρ(A) ≤ kAkfor any matrix norm.

Proof: Proof for the second part. Suppose λ is an eigenvalue of Aand x 6= 0 is a corresponding eigenvector such that Ax = λx and kxk = 1. Then

|λ| = |λ|kxk = kλxk = kAxk ≤ kAkkxk = kAk, that is, |λ| ≤ kAk. Since λ is arbitrary, this implies that ρ(A) = max |λ| ≤ kAk.

**Theorem 15**

For any A and anyε > 0, there exists a matrix norm k · k such that

ρ(A) < kAk < ρ(A) + ε.

**16 / 99**

**logo**

**Definition 16**

We call an n × n matrix A convergent if

k→∞lim(A^{k})_{ij} = 0 ∀ i = 1, 2, . . . , n and j = 1, 2, . . . , n.

**Theorem 17**

The following statements are equivalent.

**1** Ais aconvergent matrix;

**2** lim

k→∞kA^{k}k = 0forsomematrix norm;

**3** lim

k→∞kA^{k}k = 0forallmatrix norm;

**4** ρ(A) < 1;

**5** lim

k→∞A^{k}x = 0foranyx.

**logo**

**Exercise**

Page 449: 11, 12, 18, 19

**18 / 99**

**logo**

**Jacobi and Gauss-Siedel Iterative Techniques**

For small dimension of linear systems, it requires for direct techniques.

For large systems, iterative techniques are efficient in terms of both computer storage and computation.

The basic idea of iterative techniques is to split the coefficient matrix A into

A = M − (M − A),

for some matrix M , which is called thesplitting matrix. Here we assume that A and M are bothnonsingular. Then the original problem is rewritten in the equivalent form

M x = (M − A)x + b.

**logo**

This suggests an iterative process

x^{(k)}= (I − M^{−1}A)x^{(k−1)}+ M^{−1}b ≡ T x^{(k−1)}+ c,

whereT is usually called theiteration matrix. The initial vector
x^{(0)}can be arbitrary or be chosen according to certain

conditions.

Two criteria for choosing the splitting matrix M are
x^{(k)}is easily computed. More precisely, the system
M x^{(k)}= y is easy to solve;

the sequence {x^{(k)}} converges rapidly to the exact
solution.

Note that one way to achieve the second goal is to choose M
so that M^{−1}approximate A^{−1}.

In the following subsections, we will introduce some of the mostly commonly used classic iterative methods.

**20 / 99**

**logo**

**Jacobi Method**

If we decompose the coefficient matrix A as A = L + D + U,

whereDis thediagonal part,Lis thestrictly lower triangular part, andU is thestrictly upper triangular part, ofA, and chooseM = D, then we derive the iterative formulation for Jacobi method:

x^{(k)}= −D^{−1}(L + U )x^{(k−1)}+ D^{−1}b.

With this method, the iteration matrixT_{J} = −D^{−1}(L + U )and
c = D^{−1}b. Each component x^{(k)}_{i} can be computed by

x^{(k)}_{i} =

b_{i}−

i−1

X

j=1

a_{ij}x^{(k−1)}_{j} −

n

X

j=i+1

a_{ij}x^{(k−1)}_{j}

,

a_{ii}.

**logo**

a11x^{(k)}_{1} + a12x^{(k−1)}_{2} + a13x^{(k−1)}_{3} + · · · + a1nx^{(k−1)}n = b1

a_{21}x^{(k−1)}_{1} +a_{22}x^{(k)}_{2} + a_{23}x^{(k−1)}_{3} + · · · + a_{2n}x^{(k−1)}_{n} = b_{2}
...
an1x^{(k−1)}_{1} + an2x^{(k−1)}_{2} + an3x^{(k−1)}_{3} + · · · +annx^{(k)}n = bn.
**Algorithm 1 (Jacobi Method)**

Given x^{(0)}, tolerance T OL, maximum number of iteration M .
Set k = 1.

While k ≤ M and kx − x^{(0)}k_{2} ≥ T OL
Set k = k + 1, x^{(0)} = x.

For i = 1, 2, . . . , n xi=

bi−Pi−1

j=1aijx^{(0)}_{j} −Pn

j=i+1aijx^{(0)}_{j}

,

aii

End For

End While ^{22 / 99}

**logo**

**Example 18**

Consider the linear system Ax = b given by

E1 : 10x1 − x2 + 2x3 = 6,

E_{2} : −x_{1} + 11x_{2} − x_{3} + 3x_{4} = 25,
E3 : 2x1 − x2 + 10x3 − x4 = −11,

E4 : 3x2 − x3 + 8x4 = 15

which has the unique solution x = [1, 2, −1, 1]^{T}.
Solving equation E_{i} for x_{i}, for i = 1, 2, 3, 4, we obtain

x_{1} = 1/10x_{2} − 1/5x_{3} + 3/5,

x2 = 1/11x1 + 1/11x3 − 3/11x_{4} + 25/11,

x_{3} = −1/5x_{1} + 1/10x_{2} + 1/10x_{4} − 11/10,

x_{4} = − 3/8x_{2} + 1/8x_{3} + 15/8.

**logo**

Then Ax = b can be rewritten in the form x = T x + c with

T =

0 1/10 −1/5 0

1/11 0 1/11 −3/11

−1/5 1/10 0 1/10

0 −3/8 1/8 0

and c =

3/5 25/11

−11/10 15/8

and the iterative formulation for Jacobi method is
x^{(k)}= T x^{(k−1)}+ c for k = 1, 2, . . . .

The numerical results of such iteration is list as follows:

**24 / 99**

**logo**

k x1 x2 x3 x4

0 0.0000 0.0000 0.0000 0.0000 1 0.6000 2.2727 -1.1000 1.8750 2 1.0473 1.7159 -0.8052 0.8852 3 0.9326 2.0533 -1.0493 1.1309 4 1.0152 1.9537 -0.9681 0.9738 5 0.9890 2.0114 -1.0103 1.0214 6 1.0032 1.9922 -0.9945 0.9944 7 0.9981 2.0023 -1.0020 1.0036 8 1.0006 1.9987 -0.9990 0.9989 9 0.9997 2.0004 -1.0004 1.0006 10 1.0001 1.9998 -0.9998 0.9998

**logo**

**Matlab code of Example**

clear all; delete rslt.dat; diary rslt.dat; diary on;

n = 4; xold = zeros(n,1); xnew = zeros(n,1); T = zeros(n,n);

T(1,2) = 1/10; T(1,3) = -1/5; T(2,1) = 1/11;

T(2,3) = 1/11; T(2,4) = -3/11; T(3,1) = -1/5;

T(3,2) = 1/10; T(3,4) = 1/10; T(4,2) = -3/8; T(4,3) = 1/8;

c(1,1) = 3/5; c(2,1) = 25/11; c(3,1) = -11/10; c(4,1) = 15/8;

xnew = T * xold + c; k = 0;

fprintf(’ k x1 x2 x3 x4 \n’);

while ( k <= 100 & norm(xnew-xold) > 1.0d-14 ) xold = xnew; xnew = T * xold + c; k = k + 1;

fprintf(’%3.0f ’,k);

for jj = 1:n

fprintf(’%5.4f ’,xold(jj));

end

fprintf(’\n’);

end **26 / 99**

**logo**

**Gauss-Seidel Method**

When computing x^{(k)}_{i} for i > 1, x^{(k)}_{1} , . . . , x^{(k)}_{i−1}have already been
computed and are likely to be better approximations to the exact
x_{1}, . . . , x_{i−1}than x^{(k−1)}_{1} , . . . , x^{(k−1)}_{i−1} . It seems reasonable to compute
x^{(k)}_{i} using these most recently computed values. That is

a11x^{(k)}_{1} + a12x^{(k−1)}_{2} + a13x^{(k−1)}_{3} + · · · + a1nx^{(k−1)}n = b1

a21x^{(k)}_{1} +a22x^{(k)}_{2} + a23x^{(k−1)}_{3} + · · · + a2nx^{(k−1)}n = b2

a31x^{(k)}_{1} +a32x^{(k)}_{2} +a33x^{(k)}_{3} + · · · + a3nx^{(k−1)}n = b3

...
an1x^{(k)}_{1} +an2x^{(k)}_{2} +an3x^{(k)}_{3} + · · · +annx^{(k)}n = bn.
This improvement induce the Gauss-Seidel method.

The Gauss-Seidel method setsM = D + Land defines the iteration as

x^{(k)}= −(D + L)^{−1}U x^{(k−1)}+ (D + L)^{−1}b.

**logo**

That is, Gauss-Seidel method usesT_{G}= −(D + L)^{−1}U as the
iteration matrix. The formulation above can be rewritten as

x^{(k)}= −D^{−1}

Lx^{(k)}+ U x^{(k−1)}− b
.

Hence each component x^{(k)}_{i} can be computed by

x^{(k)}_{i} =

b_{i}−

i−1

X

j=1

a_{ij}x^{(k)}_{j} −

n

X

j=i+1

a_{ij}x^{(k−1)}_{j}

,

a_{ii}.

For Jacobi method, only the components of x^{(k−1)}are
used to compute x^{(k)}. Hence x^{(k)}_{i} , i = 1, . . . , n,can be
computed in parallel at each iteration k.

At each iteration of Gauss-Seidel method, since x^{(k)}_{i} can
not be computed until x^{(k)}_{1} , . . . , x^{(k)}_{i−1}are available, the
method is not a parallel algorithm in nature.

**28 / 99**

**logo**

**Algorithm 2 (Gauss-Seidel Method)**

Given x^{(0)}, tolerance T OL, maximum number of iteration M .
Set k = 1.

For i = 1, 2, . . . , n xi=

bi−Pi−1

j=1aijxj−Pn

j=i+1aijx^{(0)}_{j}

,

aii

End For

While k ≤ M and kx − x^{(0)}k_{2} ≥ T OL
Set k = k + 1, x^{(0)} = x.

For i = 1, 2, . . . , n
x_{i}=

b_{i}−Pi−1

j=1a_{ij}x_{j}−Pn

j=i+1a_{ij}x^{(0)}_{j}
,

a_{ii}
End For

End While

**logo**

**Example 19**

Consider the linear system Ax = b given by

E1 : 10x1 − x2 + 2x3 = 6,

E2 : −x_{1} + 11x2 − x3 + 3x4 = 25,
E_{3} : 2x_{1} − x_{2} + 10x_{3} − x_{4} = −11,

E_{4} : 3x_{2} − x_{3} + 8x_{4} = 15

which has the unique solution x = [1, 2, −1, 1]^{T}.
Gauss-Seidel method gives the equation

x^{(k)}_{1} = _{10}^{1}x^{(k−1)}_{2} − ^{1}_{5}x^{(k−1)}_{3} + ^{3}_{5},
x^{(k)}_{2} = _{11}^{1}x^{(k)}_{1} + _{11}^{1}x^{(k−1)}_{3} − _{11}^{3}x^{(k−1)}_{4} + ^{25}_{11},
x^{(k)}_{3} = −^{1}_{5}x^{(k)}_{1} + _{10}^{1}x^{(k)}_{2} + _{10}^{1}x^{(k−1)}_{4} − ^{11}_{10},
x^{(k)}_{4} = − ^{3}_{8}x^{(k)}_{2} + ^{1}_{8}x^{(k)}_{3} + ^{15}_{8}.

**30 / 99**

**logo**

The numerical results of such iteration is list as follows:

k x1 x2 x3 x4

0 0.0000 0.0000 0.0000 0.0000 1 0.6000 2.3273 -0.9873 0.8789 2 1.0302 2.0369 -1.0145 0.9843 3 1.0066 2.0036 -1.0025 0.9984 4 1.0009 2.0003 -1.0003 0.9998 5 1.0001 2.0000 -1.0000 1.0000 The results of Example appear to imply that the

Gauss-Seidel method is superior to the Jacobi method.

This is almost always true, but there are linear systems for which the Jacobi method converges and the Gauss-Seidel method does not.

See Exercises 17 and 18 (8th edition).

**logo**
**Matlab code of Example**

clear all; delete rslt.dat; diary rslt.dat; diary on;

n = 4; xold = zeros(n,1); xnew = zeros(n,1); A = zeros(n,n);

A(1,1)=10; A(1,2)=-1; A(1,3)=2; A(2,1)=-1; A(2,2)=11; A(2,3)=-1; A(2,4)=3; A(3,1)=2; A(3,2)=-1;

A(3,3)=10; A(3,4)=-1; A(4,2)=3; A(4,3)=-1; A(4,4)=8; b(1)=6; b(2)=25; b(3)=-11; b(4)=15;

for ii = 1:n xnew(ii) = b(ii);

for jj = 1:ii-1

xnew(ii) = xnew(ii) - A(ii,jj) * xnew(jj);

end for jj = ii+1:n

xnew(ii) = xnew(ii) - A(ii,jj) * xold(jj);

end

xnew(ii) = xnew(ii) / A(ii,ii);

end

k = 0; fprintf(’ k x1 x2 x3 x4 \n’);

while ( k <= 100 & norm(xnew-xold) > 1.0d-14 ) xold = xnew; k = k + 1;

for ii = 1:n xnew(ii) = b(ii);

for jj = 1:ii-1

xnew(ii) = xnew(ii) - A(ii,jj) * xnew(jj);

end for jj = ii+1:n

xnew(ii) = xnew(ii) - A(ii,jj) * xold(jj);

end

xnew(ii) = xnew(ii) / A(ii,ii);

end

fprintf(’%3.0f ’,k);

for jj = 1:n

fprintf(’%5.4f ’,xold(jj));

end fprintf(’\n’);

end

diary off **32 / 99**

**logo**

**Lemma 20**

Ifρ(T ) < 1, then(I − T )^{−1}exists and
(I − T )^{−1}=

∞

X

i=0

T^{i}= I + T + T^{2}+ · · · .

Proof: Let λ be an eigenvalue of T , then 1 − λ is an eigenvalue of I − T . But |λ| ≤ ρ(A) < 1, so 1 − λ 6= 0 and 0 is not an eigenvalue of I − T , which means (I − T ) is nonsingular.

Next we show that (I − T )^{−1}= I + T + T^{2}+ · · ·. Since
(I − T )

m

X

i=0

T^{i}

!

= I − T^{m+1},
and ρ(T ) < 1 implies kT^{m}k → 0 as m → ∞, we have

(I − T ) lim

m→∞

m

XT^{i}

!

= (I − T )

∞

XT^{i}

!

= I.

**logo**

**Theorem 21**

Foranyx^{(0)}∈ R^{n} , the sequence produced by
x^{(k)}= T x^{(k−1)}+ c, k = 1, 2, . . . ,

convergesto theuniquesolution of x = T x + c if and only if ρ(T ) < 1.

Proof: Suppose ρ(T ) < 1. The sequence of vectors x^{(k)}produced by
the iterative formulation are

x^{(1)}= T x^{(0)}+ c

x^{(2)}= T x^{(1)}+ c = T^{2}x^{(0)}+ (T + I)c
x^{(3)}= T x^{(2)}+ c = T^{3}x^{(0)}+ (T^{2}+ T + I)c

... In general

x^{(k)}= T^{k}x^{(0)}+ (T^{k−1}+ T^{k−2}+ · · · + T + I)c.

**34 / 99**

**logo**

Since ρ(T ) < 1, lim_{k→∞}T^{k}x^{(0)}= 0for any x^{(0)} ∈ R^{n}. By
Lemma 20,

(T^{k−1}+ T^{k−2}+ · · · + T + I)c → (I − T )^{−1}c, as k → ∞.

Therefore

k→∞lim x^{(k)}= lim

k→∞T^{k}x^{(0)}+

∞

X

j=0

T^{j}

c = (I − T )^{−1}c.

Conversely, suppose {x^{(k)}} → x = (I − T )^{−1}c. Since

x − x^{(k)}= T x + c − T x^{(k−1)}− c = T (x − x^{(k−1)}) = T^{2}(x − x^{(k−2)})

= · · · = T^{k}(x − x^{(0)}).

Let z = x − x^{(0)}. Then

k→∞lim T^{k}z = lim

k→∞(x − x^{(k)}) = 0.

It follows from theorem ρ(T ) < 1.

**logo**

**Theorem 22**

IfkT k < 1, then the sequencex^{(k)}convergestoxforanyinitial
x^{(0)}and

**1** kx − x^{(k)}k ≤ kT k^{k}kx − x^{(0)}k

**2** kx − x^{(k)}k ≤ _{1−kT k}^{kT k}^{k} kx^{(1)}− x^{(0)}k.

Proof: Since x = T x + c and x^{(k)} = T x^{(k−1)}+ c,
x − x^{(k)}= T x + c − T x^{(k−1)}− c

= T (x − x^{(k−1)})

= T^{2}(x − x^{(k−2)}) = · · · = T^{k}(x − x^{(0)}).

The first statement can then be derived

kx − x^{(k)}k = kT^{k}(x − x^{(0)})k ≤ kT k^{k}kx − x^{(0)}k.

For the second result, we first show that

kx^{(n)}− x^{(n−1)}k ≤ kT k^{n−1}kx^{(1)}− x^{(0)}k for any n ≥ 1. _{36 / 99}

**logo**

Since

x^{(n)}− x^{(n−1)}= T x^{(n−1)}+ c − T x^{(n−2)}− c

= T (x^{(n−1)}− x^{(n−2)})

= T^{2}(x^{(n−2)}− x^{(n−3)}) = · · · = T^{n−1}(x^{(1)}− x^{(0)}),
we have

kx^{(n)}− x^{(n−1)}k ≤ kT k^{n−1}kx^{(1)}− x^{(0)}k.

Let m ≥ k,
x^{(m)}− x^{(k)}

=

x^{(m)}− x^{(m−1)}
+

x^{(m−1)}− x^{(m−2)}

+ · · · +

x^{(k+1)}− x^{(k)}

=T^{m−1}

x^{(1)}− x^{(0)}

+ T^{m−2}

x^{(1)}− x^{(0)}

+ · · · + T^{k}

x^{(1)}− x^{(0)}

=

T^{m−1}+ T^{m−2}+ · · · + T^{k}

x^{(1)}− x^{(0)}
,

**logo**

hence

kx^{(m)}− x^{(k)}k

≤

kT k^{m−1}+ kT k^{m−2}+ · · · + kT k^{k}

kx^{(1)}− x^{(0)}k

=kT k^{k}

kT k^{m−k−1}+ kT k^{m−k−2}+ · · · + 1

kx^{(1)}− x^{(0)}k.

Since limm→∞x^{(m)} = x,
kx − x^{(k)}k

= lim

m→∞kx^{(m)}− x^{(k)}k

≤ lim

m→∞kT k^{k}

kT k^{m−k−1}+ kT k^{m−k−2}+ · · · + 1

kx^{(1)}− x^{(0)}k

=kT k^{k}kx^{(1)}− x^{(0)}k lim

m→∞

kT k^{m−k−1}+ kT k^{m−k−2}+ · · · + 1

=kT k^{k} 1

1 − kT kkx^{(1)}− x^{(0)}k.

This proves the second result.

**38 / 99**

**logo**

**Theorem 23**

If A isstrictly diagonal dominant, then both theJacobiand
Gauss-Seidelmethodsconvergesforanyinitial vector x^{(0)}.

Proof: By assumption, A is strictly diagonal dominant, hence aii6= 0 (otherwise A is singular) and

|aii| >

n

X

j=1,j6=i

|aij|, i = 1, 2, . . . , n.

For Jacobi method, the iteration matrix TJ= −D^{−1}(L + U )has
entries

[T_{J}]_{ij} =
(−^{a}_{a}^{ij}

ii, i 6= j,

0, i = j.

Hence

kT_{J}k_{∞}= max

1≤i≤n n

X

j=1,j6=i

a_{ij}
aii

= max

1≤i≤n

1

|aii|

n

X

j=1,j6=i

|a_{ij}| < 1,

and this implies that the Jacobi method converges.

**logo**

For Gauss-Seidel method, the iteration matrix

T_{G}= −(D + L)^{−1}U. Let λ be any eigenvalue of T_{G}and y,
kyk_{∞}= 1, is a corresponding eigenvector. Thus

T_{G}y = λy =⇒ −U y = λ(D + L)y.

Hence for i = 1, . . . , n,

−

n

X

j=i+1

aijyj = λaiiyi+ λ

i−1

X

j=1

aijyj. This gives

λaiiyi = −λ

i−1

X

j=1

aijyj−

n

X

j=i+1

aijyj

and

|λ||a_{ii}||y_{i}| ≤ |λ|

i−1

X

j=1

|a_{ij}||y_{j}| +

n

X

j=i+1

|a_{ij}||y_{j}|.

**40 / 99**

**logo**

Choose the index k such that |y_{k}| = 1 ≥ |y_{j}| (this index can
always be found since kyk∞= 1). Then

|λ||a_{kk}| ≤ |λ|

k−1

X

j=1

|a_{kj}| +

n

X

j=k+1

|a_{kj}|

which gives

|λ| ≤ Pn

j=k+1|a_{kj}|

|a_{kk}| −Pk−1

j=1|a_{kj}| <

Pn

j=k+1|a_{kj}|
Pn

j=k+1|a_{kj}| = 1

Since λ is arbitrary, ρ(T_{G}) < 1. This means the Gauss-Seidel
method converges.

The rate of convergence depends on the spectral radius of the matrix associated with the method.

One way to select a procedure to accelerate convergence is to choose a method whose associated matrix has minimal spectral radius.

**logo**

**Exercise**

Page 459: 9, 10, 11

**42 / 99**

**logo**

**Relaxation Techniques for Solving Linear Systems**

**Definition 24**

Suppose ˜x ∈ R^{n} is an approximated solution of Ax = b. Theresidual
vectorrfor ˜xisr = b − A˜x.

Let the approximate solution x^{(k,i)}produced by Gauss-Seidel method
be defined by

x^{(k,i)}=h

x^{(k)}_{1} , . . . , x^{(k)}_{i−1}, x^{(k−1)}_{i} , . . . , x^{(k−1)}_{n} iT

and

r^{(k)}_{i} =h

r^{(k)}_{1i} , r_{2i}^{(k)}, . . . , r_{ni}^{(k)}i^{T}

= b − Ax^{(k,i)}

be the corresponding residual vector. Then the mth component of
r_{i}^{(k)}is

r^{(k)}_{mi} = bm−

i−1

X

j=1

amjx^{(k)}_{j} −

n

X

j=i

amjx^{(k−1)}_{j} ,

**logo**

or, equivalently,

r^{(k)}_{mi} = bm−

i−1

X

j=1

amjx^{(k)}_{j} −

n

X

j=i+1

amjx^{(k−1)}_{j} − a_{mi}x^{(k−1)}_{i} ,

for each m = 1, 2, . . . , n.

In particular, the ith component of r_{i}^{(k)}is

r^{(k)}_{ii} = b_{i}−

i−1

X

j=1

a_{ij}x^{(k)}_{j} −

n

X

j=i+1

a_{ij}x^{(k−1)}_{j} − a_{ii}x^{(k−1)}_{i} ,
so

aiix^{(k−1)}_{i} + r^{(k)}_{ii} = bi−

i−1

X

j=1

aijx^{(k)}_{j} −

n

X

j=i+1

aijx^{(k−1)}_{j}

= a_{ii}x^{(k)}_{i} .

**44 / 99**

**logo**

Consequently, the Gauss-Seidel method can be characterized
as choosing x^{(k)}_{i} to satisfy

x^{(k)}_{i} = x^{(k−1)}_{i} +r^{(k)}_{ii}
aii

.

Relaxation method is modified the Gauss-Seidel procedure to
x^{(k)}_{i} = x^{(k−1)}_{i} + ωr_{ii}^{(k)}

aii

= x^{(k−1)}_{i} + ω
aii

bi−

i−1

X

j=1

aijx^{(k)}_{j} −

n

X

j=i+1

aijx^{(k−1)}_{j} − a_{ii}x^{(k−1)}_{i}

= (1 − ω)x^{(k−1)}_{i} + ω
aii

b_{i}−

i−1

X

j=1

a_{ij}x^{(k)}_{j} −

n

X

j=i+1

a_{ij}x^{(k−1)}_{j}

(1) for certain choices of positive ω such that the norm of the

residual vector is reduced and the convergence is significantly

**logo**

These methods are called for

ω < 1: under relaxation, ω = 1: Gauss-Seidel method, ω > 1: over relaxation.

Over-relaxation methods are calledSOR (Successive

over-relaxation). To determine the matrix of the SOR method, we rewrite (1) as

a_{ii}x^{(k)}_{i} + ω

i−1

X

j=1

a_{ij}x^{(k)}_{j} = (1 − ω)a_{ii}x^{(k−1)}_{i} − ω

n

X

j=i+1

a_{ij}x^{(k−1)}_{j} + ωb_{i},
so that if A = L + D + U , then we have

(D + ωL)x^{(k)} = [(1 − ω)D − ωU ] x^{(k−1)}+ ωb
or

x^{(k)} = (D + ωL)^{−1}[(1 − ω)D − ωU ] x^{(k−1)}+ ω(D + ωL)^{−1}b

≡ T_{ω}x^{(k−1)}+ c_{ω}.

**46 / 99**

**logo**

**Example 25**

The linear system Ax = b given by

4x1 + 3x2 = 24,

3x1 + 4x2 − x3 = 30,

− x2 + 4x3 = −24,
has the solution [3, 4, −5]^{T}.

Numerical results of Gauss-Seidel method with x^{(0)}= [1, 1, 1]^{T}:

k x_{1} x_{2} x_{3}

0 1.0000000 1.0000000 1.0000000 1 5.2500000 3.8125000 -5.0468750 2 3.1406250 3.8828125 -5.0292969 3 3.0878906 3.9267578 -5.0183105 4 3.0549316 3.9542236 -5.0114441 5 3.0343323 3.9713898 -5.0071526 6 3.0214577 3.9821186 -5.0044703

**logo**

Numerical results of SOR method withω = 1.25and
x^{(0)}= [1, 1, 1]^{T}:

k x_{1} x_{2} x_{3}

0 1.0000000 1.0000000 1.0000000 1 6.3125000 3.5195313 -6.6501465 2 2.6223145 3.9585266 -4.6004238 3 3.1333027 4.0102646 -5.0966863 4 2.9570512 4.0074838 -4.9734897 5 3.0037211 4.0029250 -5.0057135 6 2.9963276 4.0009262 -4.9982822 7 3.0000498 4.0002586 -5.0003486

**48 / 99**

**logo**

Numerical results of SOR method withω = 1.6and
x^{(0)}= [1, 1, 1]^{T}:

k x_{1} x_{2} x_{3}

0 1.0000000 1.0000000 1.0000000 1 7.8000000 2.4400000 -9.2240000 2 1.9920000 4.4560000 -2.2832000 3 3.0576000 4.7440000 -6.3324800 4 2.0726400 4.1334400 -4.1471360 5 3.3962880 3.7855360 -5.5975040 6 3.0195840 3.8661760 -4.6950272 7 3.1488384 4.0236774 -5.1735127

**logo**
**Matlab code of SOR**

clear all; delete rslt.dat; diary rslt.dat; diary on;

n = 3; xold = zeros(n,1); xnew = zeros(n,1); A = zeros(n,n); DL = zeros(n,n); DU = zeros(n,n);

A(1,1)=4; A(1,2)=3; A(2,1)=3; A(2,2)=4; A(2,3)=-1; A(3,2)=-1; A(3,3)=4;

b(1,1)=24; b(2,1)=30; b(3,1)=-24; omega=1.25;

for ii = 1:n DL(ii,ii) = A(ii,ii);

for jj = 1:ii-1

DL(ii,jj) = omega * A(ii,jj);

end

DU(ii,ii) = (1-omega)*A(ii,ii);

for jj = ii+1:n

DU(ii,jj) = - omega * A(ii,jj);

end end

c = omega * (DL \ b); xnew = DL \ ( DU * xold ) + c;

k = 0; fprintf(’ k x1 x2 x3 \n’);

while ( k <= 100 & norm(xnew-xold) > 1.0d-14 ) xold = xnew; k = k + 1; xnew = DL \ ( DU * xold ) + c;

fprintf(’%3.0f ’,k);

for jj = 1:n

fprintf(’%5.4f ’,xold(jj));

end fprintf(’\n’);

end diary off

**50 / 99**

**logo**

**Theorem 26 (Kahan)**

If aii6= 0, for each i = 1, 2, . . . , n, then ρ(Tω) ≥ |ω − 1|. This implies that the SOR method can converge only if 0 < ω < 2.

**Theorem 27 (Ostrowski-Reich)**

If A ispositive definiteand the relaxation parameter ω satisfying
0 < ω < 2, then theSORiterationconvergesforanyinitial vector x^{(0)}.

**Theorem 28**

If A ispositive definiteandtridiagonal, thenρ(T_{G}) = [ρ(T_{J})]^{2}< 1and
theoptimalchoice of ω for the SOR iteration is

ω = 2

1 + q

1 − [ρ(TJ)]^{2}
.

With this choice of ω,ρ(Tω) = ω − 1.

**logo**

**Example 29**
The matrix

A =

4 3 0

3 4 −1

0 −1 4

,

given in previous example, is positive definite and tridiagonal.

Since

TJ = −D^{−1}(L + U ) =

1

4 0 0

0 ^{1}_{4} 0
0 0 ^{1}_{4}

0 −3 0

−3 0 1

0 1 0

=

0 −0.75 0

−0.75 0 0.25

0 0.25 0

,

**52 / 99**

**logo**

we have

T_{J} − λI =

−λ −0.75 0

−0.75 −λ 0.25

0 0.25 −λ

, so

det(T_{J} − λI) = −λ(λ^{2}− 0.625).

Thus,

ρ(T_{J}) =√
0.625
and

ω = 2

1 +p1 − [ρ(T_{J})]^{2} = 2
1 +√

1 − 0.625 ≈ 1.24.

This explains the rapid convergence obtained in previous example when using ω = 0.125

**logo**

**Symmetric Successive Over Relaxation (SSOR)** **Method**

Let A be symmetric and A = D + L + L^{T}. The idea is in fact to
implement the SOR formulationtwice,one forwardandone
backward, at each iteration. That is, SSOR method defines

(D + ωL)x^{(k−}^{1}^{2}^{)} = (1 − ω)D − ωL^{T} x^{(k−1)}+ ωb, (2)
(D + ωL^{T})x^{(k)} = [(1 − ω)D − ωL] x^{(k−}^{1}^{2}^{)}+ ωb. (3)
Define

Mω: = D + ωL,

N_{ω}: = (1 − ω)D − ωL^{T}.
Then from the iterations (2) and (3), it follows that

x^{(k)} = M_{ω}^{−T}N_{ω}^{T}M_{ω}^{−1}N_{ω} x^{(k−1)}+ ω M_{ω}^{−T}N_{ω}^{T}M_{ω}^{−1}+ M_{ω}^{−T} b

≡ T (ω)x^{(k−1)}+ M (ω)^{−1}b.

**54 / 99**

**logo**

But

((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 the splitting matrix is

M (ω) = 1

ω(2 − ω)(D + ωL)D^{−1} D + ωL^{T} .
The iteration matrix is

T (ω) = (D + ωL^{T})^{−1}[(1 − ω)D − ωL] (D + ωL)^{−1}(1 − ω)D − ωL^{T} .

**logo**

**Exercise**
Page 467: 2, 8

**56 / 99**

**logo**

**Error bounds and iterative refinement**

**Example 30**

The linear system Ax = b given by

1 2

1.0001 2

x1

x_{2}

=

3

3.0001

has the unique solution x = [1, 1]^{T}.

The poor approximation ˜x = [3, 0]^{T} has the residual vector
r = b − A˜x =

3

3.0001

−

1 2

1.0001 2

3 0

=

0

−0.0002

,
so krk∞= 0.0002. Although the norm of the residual vector is
small, the approximation ˜x = [3, 0]^{T} is obviously quite poor; in
fact, kx − ˜xk∞= 2.

**logo**

The solution of above example represents the intersection of the lines

`1: x1+ 2x2 = 3 and `2 : 1.0001x1+ 2x2 = 3.0001.

`1 and `_{2} are nearly parallel. The point (3, 0) lies on `_{1} which
implies that (3, 0) also lies close to `_{2}, even though it differs
significantly from the intersection point (1, 1).

**58 / 99**

**logo**

**Theorem 31**

Suppose that ˜xis an approximate solution of Ax = b, A is nonsingular matrix and r = b − A˜x. Then

kx − ˜xk ≤ krk · kA^{−1}k
and if x 6= 0 and b 6= 0,

kx − ˜xk

kxk ≤ kAk · kA^{−1}kkrk
kbk.
Proof: Since

r = b − A˜x = Ax − A˜x = A(x − ˜x) and A is nonsingular, we have

kx − ˜xk = kA^{−1}rk ≤ kA^{−1}k · krk. (4)
Moreover, since b = Ax, we have

kbk ≤ kAk · kxk.

**logo**

It implies that

1

kxk ≤ kAk

kbk. (5)

Combining Equations (4) and (5), we have kx − ˜xk

kxk ≤ kAk · kA^{−1}k
kbk krk.

**Definition 32 (Condition number)**

The condition number of nonsingular matrix A is
κ(A) = kAk · kA^{−1}k.

For any nonsingular matrix A,

1 = kIk = kA · A^{−1}k ≤ kAk · kA^{−1}k = κ(A).

**60 / 99**

**logo**

**Definition 33**

A matrix A iswell-conditionedifκ(A)is close to1, and is ill-conditionedwhen κ(A) is significantly greater than 1.

In previous example, A =

1 2

1.0001 2

. Since

A^{−1} =

−10000 10000 5000.5 −5000

, we have

κ(A) = kAk∞· kA^{−1}k_{∞}= 3.0001 × 20000 = 60002 1.

**logo**

How to estimate the effective condition number in t-digit arithmetic without having to invert the matrix A?

If the approximate solution ˜xof Ax = b is being determined using t-digit arithmetic and Gaussian elimination, then

krk = kb − A˜xk ≈ 10^{−t}kAk · k˜xk.

All the arithmetic operations in Gaussian elimination technique are performed using t-digit arithmetic, but the residual vector r are done in double-precision (i.e., 2t-digit) arithmetic.

Use the Gaussian elimination method which has already been calculated to solve

Ay = r.

Let ˜ybe the approximate solution.

**62 / 99**

**logo**

Then

˜

y ≈ A^{−1}r = A^{−1}(b − A˜x) = x − ˜x
and

x ≈ ˜x + ˜y.

Moreover,

k˜yk ≈ kx − ˜xk = kA^{−1}rk

≤ kA^{−1}k · krk ≈ kA^{−1}k(10^{−t}kAk · k˜xk) = 10^{−t}k˜xkκ(A).

It implies that

κ(A) ≈ k˜yk
k˜xk10^{t}.
**Iterative refinement**

In general, ˜x + ˜yis a more accurate approximation to the solution of Ax = b than ˜x.

**logo**

**Algorithm 3 (Iterative refinement)**

Given tolerance T OL, maximum number of iteration M , number of digits of precision t.

Solve Ax = b by using Gaussian elimination in t-digit arithmetic.

Set k = 1 while ( k ≤ M )

Compute r = b − Ax in 2t-digit arithmetic.

Solve Ay = r by using Gaussian elimination in t-digit arithmetic.

If kyk∞< T OL, then stop.

Set k = k + 1 and x = x + y.

End while

**64 / 99**