university-logo

Iterative techniques in matrix algebra

Tsung-Min Hwang

Department of Mathematics National Taiwan Normal University, Taiwan

February 16, 2011

university-logo

Outline

university-logo

Definition

∥ · ∥ : ℝ^{n} → ℝ is a vector norm if (i) ∥x∥ ≥ 0, ∀ x ∈ ℝ^{n}, (ii) ∥x∥ = 0 if and only if x = 0, (iii) ∥αx| = |α|∥x∥ ∀ α ∈ ℝ and x ∈ ℝ^{n}, (iv) ∥x + y∥ ≤ ∥x∥ + ∥y∥ ∀ x, y ∈ ℝ^{n}.

Definition

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

∥x∥_{2} = (x^{T}x)^{1/2} = {∑_{i=1}^{n}x_{i}^{2}}^{1/2} and ∥x∥_{∞} = max_{1≤i≤n}|x_{i}|.

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

university-logo

Theorem (Cauchy-Bunyakovsky-Schwarz inequality)

For each x = [x_{1}, x_{2}, ⋯, x_{n}]^{T} and y = [y_{1}, y_{2}, ⋯, y_{n}]^{T} in ℝ^{n},

x^{T}y = ∑_{i=1}^{n}x_{i}y_{i} ≤ {∑_{i=1}^{n}x_{i}^{2}}^{1/2} {∑_{i=1}^{n}y_{i}^{2}}^{1/2} = ∥x∥_{2} · ∥y∥_{2}.

Proof: If x = 0 or y = 0, the result is immediate. Suppose x ≠ 0 or y ≠ 0. For each α ∈ ℝ,

0 ≤ ∥x − αy∥_{2}^{2} = ∑_{i=1}^{n}(x_{i} − αy_{i})^{2} = ∑_{i=1}^{n}x_{i}^{2} − 2α ∑_{i=1}^{n}x_{i}y_{i} + α^{2} ∑_{i=1}^{n}y_{i}^{2},

and

2α ∑_{i=1}^{n}x_{i}y_{i} ≤ ∑_{i=1}^{n}x_{i}^{2} + α^{2} ∑_{i=1}^{n}y_{i}^{2} = ∥x∥_{2}^{2} + α^{2}∥y∥_{2}^{2}.

university-logo

Since ∥x∥_{2} > 0 and ∥y∥_{2} > 0, we can let

α = \frac{∥x∥_{2}}{∥y∥_{2}}

to give

(2\frac{∥x∥_{2}}{∥y∥_{2}}) (∑_{i=1}^{n}x_{i}y_{i}) ≤ ∥x∥_{2}^{2} + \frac{∥x∥_{2}^{2}}{∥y∥_{2}^{2}}∥y∥_{2}^{2} = 2∥x∥_{2}^{2}.

Thus

x^{T}y = ∑_{i=1}^{n}x_{i}y_{i} ≤ ∥x∥_{2}∥y∥_{2}.

university-logo

For each x, y ∈ ℝ^{n},

∥x + y∥_{∞} = 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}| = ∥x∥_{∞} + ∥y∥_{∞}

and

∥x + y∥_{2}^{2} = ∑_{i=1}^{n}(x_{i} + y_{i})^{2} = ∑_{i=1}^{2}x_{i}^{2} + 2 ∑_{i=1}^{n}x_{i}y_{i} + ∑_{i=1}^{n}y_{i}^{2}
≤ ∥x∥_{2}^{2} + 2∥x∥_{2}∥y∥_{2} + ∥y∥_{2}^{2} = (∥x∥_{2} + ∥y∥_{2})^{2},

which gives

∥x + y∥_{2} ≤ ∥x∥_{2} + ∥y∥_{2}.

Definition

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

∥x^{(k)} − x∥ < ε, ∀ k ≥ N(ε).
T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  6 / 1

university-logo

Theorem

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

lim_{k→∞}x_{i}^{(k)} = x_{i}, ∀ i = 1, 2, …, n.

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

max_{1≤i≤n}|x_{i}^{(k)} − x_{i}| = ∥x^{(k)} − x∥_{∞} < ε, ∀ k ≥ N(ε).

This result implies that

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

Hence

lim_{k→∞}x_{i}^{(k)} = x_{i}, ∀ i.

university-logo

“⇐” For a given ε > 0, let N_{i}(ε) represent an integer with

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

Define

N(ε) = max_{1≤i≤n}N_{i}(ε).

If k ≥ N(ε), then

max_{1≤i≤n}|x_{i}^{(k)} − x_{i}| = ∥x^{(k)} − x∥_{∞} < ε.

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

university-logo

Theorem

For each x ∈ ℝ^{n},

∥x∥_{∞} ≤ ∥x∥_{2} ≤ \sqrt{n}∥x∥_{∞}.

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

∥x∥_{∞}^{2} = |x_{j}|^{2} ≤ ∑_{i=1}^{n}x_{i}^{2} = ∥x∥_{2}^{2},

so ∥x∥_{∞} ≤ ∥x∥_{2} and

∥x∥_{2}^{2} = ∑_{i=1}^{n}x_{i}^{2} ≤ ∑_{i=1}^{n}x_{j}^{2} = nx_{j}^{2} = n∥x∥_{∞}^{2},

so ∥x∥_{2} ≤ \sqrt{n}∥x∥_{∞}.

university-logo

Definition

A matrix norm ∥ · ∥ 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) ∥A∥ ≥ 0; (ii) ∥A∥ = 0 if and only if A = 0; (iii) ∥αA∥ = |α|∥A∥; (iv) ∥A + B∥ ≤ ∥A∥ + ∥B∥;

(v) ∥AB∥ ≤ ∥A∥∥B∥;

Theorem

If ∥ · ∥ is a vector norm on ℝ^{n}, then

∥A∥ = max_{∥x∥=1}∥Ax∥

is a matrix norm.

university-logo

For any z ≠ 0, we have x = z/∥z∥ as a unit vector. Hence

∥A∥ = max_{∥x∥=1}∥Ax∥ = max_{z≠0} ∥A (\frac{z}{∥z∥})∥ = max_{z≠0} \frac{∥Az∥}{∥z∥}.

Corollary

∥Az∥ ≤ ∥A∥ · ∥z∥.

Theorem

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

∥A∥_{∞} = max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}|.

university-logo

Proof: Let x be an n-dimension vector with

1 = ∥x∥_{∞} = max_{1≤i≤n}|x_{i}|.

Then

∥Ax∥_{∞} = max_{1≤i≤n} |∑_{j=1}^{n}a_{ij}x_{j} |

≤ max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}| max_{1≤j≤n}|x_{j}| = max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}|.

Consequently,

∥A∥_{∞} = max_{∥x∥_{∞}=1}∥Ax∥_{∞} ≤ max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}|.

On the other hand, let p be an integer with

∑_{j=1}^{n}|a_{pj}| = max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}|,

university-logo

and x be the vector with

x_{j} = \{ \begin{matrix}1, & if a_{pj} ≥ 0, \\ −1, & if a_{pj} < 0. \end{matrix}

Then

∥x∥_{∞} = 1 and a_{pj}x_{j} = |a_{pj}|, ∀ j = 1, 2, …, n,

so ∥Ax∥_{∞} = max_{1≤i≤n} |∑_{j=1}^{n}a_{ij}x_{j} | ≥ |∑_{j=1}^{n}a_{pj}x_{j} | = |∑_{j=1}^{n}|a_{pj}|| = max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}|. This result implies that

∥A∥_{∞} = max_{∥x∥_{∞}=1}∥Ax∥_{∞} ≥ max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}|.

which gives

∥A∥_{∞} = max_{1≤i≤n} ∑_{j=1}^{n}|a_{ij}|.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  13 / 1

university-logo

Eigenvalues and eigenvectors

Definition (Characteristic polynomial)

If A is a square matrix, the characteristic polynomial of A is defined by

p(λ) = det(A − λI ).

Definition (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 ≠ 0 satisfies (A − λI )x = 0, then x is an eigenvector of A corresponding to the eigenvalue λ.

Definition (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}.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  14 / 1

university-logo

Theorem

If A is an n × n matrix, then (i) ∥A∥_{2} = \sqrt{ρ(A^{T}A)}; (ii) ρ(A) ≤ ∥A∥ for any matrix norm.

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

|λ| = |λ|∥x∥ = ∥λx∥ = ∥Ax∥ ≤ ∥A∥∥x∥ = ∥A∥,

that is, |λ| ≤ ∥A∥. Since λ is arbitrary, this implies that ρ(A) = max |λ| ≤ ∥A∥.

Theorem

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

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

university-logo

Definition

We call an n × n matrix A convergent if

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

Theorem

The following statements are equivalent. 1 A is a convergent matrix; 2 lim_{k→∞}∥A^{k}∥ = 0 for some matrix norm;

3 lim_{k→∞}∥A^{k}∥ = 0 for all matrix norm; 4 ρ(A) < 1; 5 lim_{k→∞}A^{k}x = 0 for any x.

university-logo

Iterative techniques for solving linear systems

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 the splitting matrix. Here we assume that A and M are both nonsingular. Then the original problem is rewritten in the equivalent form

Mx = (M − A)x + b.

This suggests an iterative process

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

where T is usually called the iteration matrix. The initial vector x^{(0)} can be arbitrary or be chosen according to certain conditions.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  17 / 1

university-logo

Two criteria for choosing the splitting matrix M are x^{(k)} is easily computed. More precisely, the system Mx^{(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.

university-logo

Jacobi Method

If we decompose the coefficient matrix A as

A = L + D + U,

where D is the diagonal part, L is the strictly lower triangular part, and U is the strictly upper triangular part, of A, and choose M = 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 matrix T_{J} = −D^{−1}(L + U) and c = D^{−1}b. Each component x_{i}^{(k)} can be computed by

x_{i}^{(k)} = (b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k−1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k−1)})/a_{ii}.

university-logo

a_{11}x_{1}^{(k)} + a_{12}x_{2}^{(k−1)} + a_{13}x_{3}^{(k−1)} + ⋯ + a_{1n}x_{n}^{(k−1)} = b_{1}
a_{21}x_{1}^{(k−1)} + a_{22}x_{2}^{(k)} + a_{23}x_{3}^{(k−1)} + ⋯ + a_{2n}x_{n}^{(k−1)} = b_{2}
⋮
a_{n1}x_{1}^{(k−1)} + a_{n2}x_{2}^{(k−1)} + a_{n3}x_{3}^{(k−1)} + ⋯ + a_{nn}x_{n}^{(k)} = b_{n}.

Algorithm (Jacobi Method)

Given x^{(0)}, tolerance TOL, maximum number of iteration M. Set k = 1. While k ≤ M and ∥x − x^{(0)}∥_{2} ≥ TOL

Set k = k + 1, x^{(0)} = x.

For i = 1, 2, …, n

x_{i} = (b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j}^{(0)} − ∑_{j=i+1}^{n} a_{ij}x_{j}^{(0)})/a_{ii}

End For End While

university-logo

Example

Consider the linear system Ax = b given by

E_{1} : 10x_{1} − x_{2} + 2x_{3}  =  6,
E_{2} : −x_{1} + 11x_{2} − x_{3} + 3x_{4} = 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}.

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,
x_{2} = 1/11x_{1}  + 1/11x_{3} − 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.

university-logo

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

T = [ \begin{matrix}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 \end{matrix} ] and c = [ \begin{matrix}3/5 \\ 25/11 \\ −11/10 \\ 15/8 \end{matrix} ]

and the iterative formulation for Jacobi method is

x^{(k)} = Tx^{(k−1)} + c for k = 1, 2, ….

The numerical results of such iteration is list as follows:

university-logo

k  x_{1}  x_{2}  x_{3}  x_{4}
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

university-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

university-logo

Gauss-Seidel Method

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

a_{11}x_{1}^{(k)} + a_{12}x_{2}^{(k−1)} + a_{13}x_{3}^{(k−1)} + ⋯ + a_{1n}x_{n}^{(k−1)} = b_{1}
a_{21}x_{1}^{(k)} + a_{22}x_{2}^{(k)} + a_{23}x_{3}^{(k−1)} + ⋯ + a_{2n}x_{n}^{(k−1)}  = b_{2}
a_{31}x_{1}^{(k)} + a_{32}x_{2}^{(k)} + a_{33}x_{3}^{(k)} + ⋯ + a_{3n}x_{n}^{(k−1)}  = b_{3}
⋮
a_{n1}x_{1}^{(k−1)} + a_{n2}x_{2}^{(k−1)} + a_{n3}x_{3}^{(k−1)} + ⋯ + a_{nn}x_{n}^{(k)} = b_{n}.

This improvement induce the Gauss-Seidel method. The Gauss-Seidel method sets M = D + L and defines the iteration as

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

university-logo

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

x^{(k)} = −D^{−1} (Lx^{(k)} + Ux^{(k−1)} − b).

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

x_{i}^{(k)} = (b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k−1)})/a_{ii}.

For Jacobi method, only the components of x^{(k−1)} are used to compute x^{(k)}. Hence x_{i}^{(k)}, i = 1, …, n, can be computed in parallel at each iteration k. At each iteration of Gauss-Seidel method, since x_{i}^{(k)} can not be computed until x_{1}^{(k)}, …, x_{i−1}^{(k)} are available, the method is not a parallel algorithm in nature.

university-logo

Algorithm (Gauss-Seidel Method)

Given x^{(0)}, tolerance TOL, maximum number of iteration M. Set k = 1. For i = 1, 2, …, n

x_{i} = (b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j} − ∑_{j=i+1}^{n} a_{ij}x_{j}^{(0)})/a_{ii}

End For While k ≤ M and ∥x − x^{(0)}∥_{2} ≥ TOL

Set k = k + 1, x^{(0)} = x.

For i = 1, 2, …, n

x_{i} = (b_{i} − ∑_{j=1}^{i−1} a_{ij}x_{j} − ∑_{j=i+1}^{n} a_{ij}x_{j}^{(0)})/a_{ii}

End For End While

university-logo

Example

Consider the linear system Ax = b given by

E_{1} : 10x_{1} − x_{2} + 2x_{3}  =  6,
E_{2} : −x_{1} + 11x_{2} − x_{3} + 3x_{4} = 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_{1}^{(k)} =  \frac{1}{10}x_{2}^{(k−1)} − \frac{1}{5}x_{3}^{(k−1)}  + \frac{3}{5},
x_{2}^{(k)} = \frac{1}{11}x_{1}^{(k)}  + \frac{1}{11}x_{3}^{(k−1)} − \frac{3}{11}x_{4}^{(k−1)} + \frac{25}{11},
x_{3}^{(k)} = −\frac{1}{5}x_{1}^{(k)} +  \frac{1}{10}x_{2}^{(k)}  + \frac{1}{10}x_{4}^{(k−1)} − \frac{11}{10},
x_{4}^{(k)} =  −  \frac{3}{8}x_{2}^{(k)} +  \frac{1}{8}x_{3}^{(k)}  + \frac{15}{8}.

university-logo

The numerical results of such iteration is list as follows:

k  x_{1}  x_{2}  x_{3}  x_{4}
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.

university-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

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  30 / 1

university-logo

Lemma

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

(I − T )^{−1} = ∑_{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 − λ ≠ 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 ) (∑_{i=0}^{m}T^{i}) = I − T^{m+1},

and ρ(T ) < 1 implies ∥T^{m}∥ → 0 as m → ∞, we have

(I − T ) ( lim_{m→∞} ∑_{i=0}^{m}T^{i}) = (I − T ) (∑_{i=0}^{∞}T^{i}) = I .

university-logo

Theorem

For any x^{(0)} ∈ ℝ^{n} , the sequence produced by

x^{(k)} = Tx^{(k−1)} + c, k = 1, 2, …,

converges to the unique solution of x = Tx + 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)} = Tx^{(0)} + c
x^{(2)} = Tx^{(1)} + c = T^{2}x^{(0)} + (T + I )c
x^{(3)} = Tx^{(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.
T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  32 / 1

university-logo

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

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

Therefore

lim_{k→∞}x^{(k)} = lim_{k→∞}T^{k}x^{(0)} + (∑_{j=0}^{∞}T^{j} )c = (I − T )^{−1}c.

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

x − x^{(k)} = Tx + c − Tx^{(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

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

It follows from theorem ρ(T ) < 1.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  33 / 1

university-logo

Theorem

If ∥T ∥ < 1, then the sequence x^{(k)} converges to x for any initial x^{(0)} and 1 ∥x − x^{(k)}∥ ≤ ∥T ∥^{k}∥x − x^{(0)}∥ 2 ∥x − x^{(k)}∥ ≤ \frac{∥T ∥^{k}}{1−∥T ∥}∥x^{(1)} − x^{(0)}∥.

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

x − x^{(k)} = Tx + c − Tx^{(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

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

For the second result, we first show that

∥x^{(n)} − x^{(n−1)}∥ ≤ ∥T ∥^{n−1}∥x^{(1)} − x^{(0)}∥ for any n ≥ 1.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  34 / 1

university-logo

Since

x^{(n)} − x^{(n−1)} = Tx^{(n−1)} + c − Tx^{(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

∥x^{(n)} − x^{(n−1)}∥ ≤ ∥T ∥^{n−1}∥x^{(1)} − x^{(0)}∥.

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)}),

university-logo

hence

∥x^{(m)} − x^{(k)}∥
≤ (∥T ∥^{m−1} + ∥T ∥^{m−2} + ⋯ + ∥T ∥^{k})∥x^{(1)} − x^{(0)}∥
= ∥T ∥^{k} (∥T ∥^{m−k−1} + ∥T ∥^{m−k−2} + ⋯ + 1)∥x^{(1)} − x^{(0)}∥.

Since lim_{m→∞} x^{(m)} = x,

∥x − x^{(k)}∥
= lim_{m→∞}∥x^{(m)} − x^{(k)}∥
≤ lim_{m→∞}∥T ∥^{k} (∥T ∥^{m−k−1} + ∥T ∥^{m−k−2} + ⋯ + 1)∥x^{(1)} − x^{(0)}∥
= ∥T ∥^{k}∥x^{(1)} − x^{(0)}∥ lim_{m→∞}(∥T ∥^{m−k−1} + ∥T ∥^{m−k−2} + ⋯ + 1)
= ∥T ∥^{k} \frac{1}{1 − ∥T ∥}∥x^{(1)} − x^{(0)}∥.

This proves the second result.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  36 / 1

university-logo

Theorem

If A is strictly diagonal dominant, then both the Jacobi and Gauss-Seidel methods converges for any initial vector x^{(0)}.

Proof: By assumption, A is strictly diagonal dominant, hence a_{ii} ≠ 0 (otherwise A is singular) and

|a_{ii}| > ∑_{j=1,j≠i}^{n} |a_{ij}|, i = 1, 2, …, n.

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

[T_{J}]_{ij} = \begin{cases}−\frac{a_{ij}}{a_{ii}} & i ≠ j \\ 0 & i = j \end{cases}

Hence

∥T_{J}∥_{∞} = max_{1≤i≤n} ∑_{j=1,j≠i}^{n} |\frac{a_{ij}}{a_{ii}} | = max_{1≤i≤n} \frac{1}{|a_{ii}|} ∑_{j=1,j≠i}^{n} |a_{ij}| < 1,

and this implies that the Jacobi method converges.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  37 / 1

university-logo

For Gauss-Seidel method, the iteration matrix T_{G} = −(D + L)^{−1}U. Let λ be any eigenvalue of T_{G} and y, ∥y∥_{∞} = 1, is a corresponding eigenvector. Thus

T_{G}y = λy ⟹ −Uy = λ(D + L)y.

Hence for i = 1, …, n,

− ∑_{j=i+1}^{n}a_{ij}y_{j} = λa_{ii}y_{i} + λ ∑_{j=1}^{i−1}a_{ij}y_{j}.

This gives

λa_{ii}y_{i} = −λ ∑_{j=1}^{i−1}a_{ij}y_{j} − ∑_{j=i+1}^{n}a_{ij}y_{j}

and

|λ||a_{ii}||y_{i}| ≤ |λ| ∑_{j=1}^{i−1}|a_{ij}||y_{j}| + ∑_{j=i+1}^{n}|a_{ij}||y_{j}|.

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

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  38 / 1

university-logo

|λ||a_{kk}| ≤ |λ| ∑_{j=1}^{k−1}|a_{kj}| + ∑_{j=k+1}^{n}|a_{kj}|

which gives

|λ| ≤ \frac{∑_{j=k+1}^{n} |a_{kj}|}{|a_{kk}| − ∑_{j=1}^{k−1} |a_{kj}|} < \frac{∑_{j=k+1}^{n} |a_{kj}|}{∑_{j=k+1}^{n} |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.

university-logo

Successive over-relaxation (SOR) method

Definition

Suppose x̃ ∈ ℝ^{n} is an approximated solution of Ax = b. The residual vector r for x̃ is r = b − Ax̃.

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

x^{(k,i)} = [x_{1}^{(k)}, …, x_{i−1}^{(k)}, x_{i}^{(k−1)}, …, x_{n}^{(k−1)}]^{T} .

and

r_{i}^{(k)} = [r_{1i}^{(k)}, r_{2i}^{(k)}, …, r_{ni}^{(k)}]^{T} = b − Ax^{(k,i)}

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

r_{mi}^{(k)} = b_{m} − ∑_{j=1}^{i−1}a_{mj}x_{j}^{(k)} − ∑_{j=i}^{n}a_{mj}x_{j}^{(k−1)},

university-logo

or, equivalently,

r_{mi}^{(k)} = b_{m} − ∑_{j=1}^{i−1}a_{mj}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{mj}x_{j}^{(k−1)} − a_{mi}x_{i}^{(k−1)},

for each m = 1, 2, …, n. In particular, the ith component of r_{i}^{(k)} is

r_{ii}^{(k)} = b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k−1)} − a_{ii}x_{i}^{(k−1)},

so

a_{ii}x_{i}^{(k−1)} + r_{ii}^{(k)} = b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k−1)}
= a_{ii}x_{i}^{(k)}.

university-logo

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

x_{i}^{(k)} = x_{i}^{(k−1)} + \frac{r_{ii}^{(k)}}{a_{ii}}.

Relaxation method is modified the Gauss-Seidel procedure to

x_{i}^{(k)} = x_{i}^{(k−1)} + ω \frac{r_{ii}^{(k)}}{a_{ii}}
= x_{i}^{(k−1)} + \frac{ω}{a_{ii}} [b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k−1)} − a_{ii}x_{i}^{(k−1)}]

= (1 − ω)x_{i}^{(k−1)} + \frac{ω}{a_{ii}} [b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k−1)}] (1)

for certain choices of positive ω such that the norm of the residual vector is reduced and the convergence is significantly faster.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  42 / 1

university-logo

These methods are called for ω < 1: under relaxation, ω = 1: Gauss-Seidel method, ω > 1: over relaxation.

Over-relaxation methods are called SOR (Successive over-relaxation). To determine the matrix of the SOR method, we rewrite (??) as

a_{ii}x_{i}^{(k)} + ω ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} = (1 − ω)a_{ii}x_{i}^{(k−1)} − ω ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k−1)} + ω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_{ω}.

university-logo

Example

The linear system Ax = b given by

4x_{1} + 3x_{2}  = 24,
3x_{1} + 4x_{2} − x_{3} = 30,
− x_{2} + 4x_{3} = −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
7 3.0134110 3.9888241 -5.0027940
T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  44 / 1

university-logo

Numerical results of SOR method with ω = 1.25 and 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

university-logo

Numerical results of SOR method with ω = 1.6 and 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

university-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 university-logo

Theorem (Kahan)

If a_{ii} ≠ 0, for each i = 1, 2, …, n, then ρ(T_{ω}) ≥ |ω − 1|. This implies that the SOR method can converge only if 0 < ω < 2.

Theorem (Ostrowski-Reich)

If A is positive definite and the relaxation parameter ω satisfying 0 < ω < 2, then the SOR iteration converges for any initial vector x^{(0)}.

Theorem

If A is positive definite and tridiagonal, then ρ(T_{G}) = [ρ(T_{J})]^{2} < 1 and the optimal choice of ω for the SOR iteration is

ω = \frac{2}{1 + \sqrt{1 − [ρ(T_{J})]^{2}}} .

With this choice of ω, ρ(T_{ω}) = ω − 1.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  48 / 1

university-logo

Example

The matrix

A = [ \begin{matrix}4 & 3 & 0 \\ 3 & 4 & −1 \\ 0 & −1 & 4 \end{matrix} ] ,

given in previous example, is positive definite and tridiagonal.

Since

T_{J} = −D^{−1}(L + U) = [ \begin{matrix}\frac{1}{4} & 0 & 0 \\ 0 & \frac{1}{4} & 0 \\ 0 & 0 & \frac{1}{4} \end{matrix} ] [ \begin{matrix}0 & −3 & 0 \\ −3 & 0 & 1 \\ 0 & 1 & 0 \end{matrix} ]

= [ \begin{matrix}0 & −0.75 & 0 \\ −0.75 & 0 & 0.25 \\ 0 & 0.25 & 0 \end{matrix} ] ,

university-logo

we have

T_{J} − λI = [ \begin{matrix}−λ & −0.75 & 0 \\ −0.75 & −λ & 0.25 \\ 0 & 0.25 & −λ \end{matrix} ] ,

so

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

Thus,

ρ(T_{J}) = \sqrt{0.625}

and

ω = \frac{2}{1 + \sqrt{1 − [ρ(T_{J})]^{2}}} = \frac{2}{1 + \sqrt{1 − 0.625}} ≈ 1.24.

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

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  50 / 1

university-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 formulation twice, one forward and one backward, at each iteration. That is, SSOR method defines

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

Define

\{ \begin{matrix}M_{ω}: = D + ωL, \\ N_{ω}: = (1 − ω)D − ωL^{T}. \end{matrix}

Then from the iterations (??) and (??), it follows that

x_{i+1} = (M_{ω}^{−T} N_{ω}^{T} M_{ω}^{−1}N_{ω})x_{i} + ω (M_{ω}^{−T} N_{ω}^{T} M_{ω}^{−1} + M_{ω}^{−T} )b
≡ T (ω)x_{i} + M(ω)^{−1}b.

university-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(ω) = \frac{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}].

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  52 / 1

university-logo

Error bounds and iterative refinement

Example

The linear system Ax = b given by

[ \begin{matrix}1 & 2 \\ 1.0001 & 2 \end{matrix} ] [ \begin{matrix}x_{1} \\ x_{2} \end{matrix} ] = [ \begin{matrix}3 \\ 3.0001 \end{matrix} ]

has the unique solution x = [1, 1]^{T}. The poor approximation x̃ = [3, 0]^{T} has the residual vector

r = b − Ax̃ = [ \begin{matrix}3 \\ 3.0001 \end{matrix} ] − [ \begin{matrix}1 & 2 \\ 1.0001 & 2 \end{matrix} ] [ \begin{matrix}3 \\ 0 \end{matrix} ] = [ \begin{matrix}0 \\ −0.0002 \end{matrix} ] ,

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

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  53 / 1

university-logo

The solution of above example represents the intersection of the lines

ℓ_{1} : x_{1} + 2x_{2} = 3 and ℓ_{2} : 1.0001x_{1} + 2x_{2} = 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).

university-logo

Theorem

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

∥x − x̃∥ ≤ ∥r∥ · ∥A^{−1}∥

and if x ≠ 0 and b ≠ 0,

\frac{∥x − x̃∥}{∥x∥} ≤ ∥A∥ · ∥A^{−1}∥\frac{∥r∥}{∥b∥}.

Proof: Since

r = b − Ax̃ = Ax − Ax̃ = A(x − x̃)

and A is nonsingular , we have

∥x − x̃∥ = ∥A^{−1}r∥ ≤ ∥A^{−1}∥ · ∥r∥.  (4)

Moreover, since b = Ax, we have

∥b∥ ≤ ∥A∥ · ∥x∥.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  55 / 1

university-logo

It implies that

\frac{1}{∥x∥} ≤ \frac{∥A∥}{∥b∥}.  (5)

Combining Equations (??) and (??), we have

\frac{∥x − x̃∥}{∥x∥} ≤ \frac{∥A∥ · ∥A^{−1}∥}{∥b∥}∥r∥.

Definition (Condition number)

The condition number of nonsingular matrix A is

κ(A) = ∥A∥ · ∥A^{−1}∥.

For any nonsingular matrix A,

1 = ∥I ∥ = ∥A · A^{−1}∥ ≤ ∥A∥ · ∥A^{−1}∥ = κ(A).

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  56 / 1

university-logo

Definition

A matrix A is well-conditioned if κ(A) is close to 1, and is ill-conditioned when κ(A) is significantly greater than 1.

In previous example,

A = [ \begin{matrix}1 & 2 \\ 1.0001 & 2 \end{matrix} ] .

Since

A^{−1} = [ \begin{matrix}−10000 & 10000 \\ 5000.5 & −5000 \end{matrix} ] ,

we have

κ(A) = ∥A∥_{∞} · ∥A^{−1}∥_{∞} = 3.0001 × 20000 = 60002 ≫ 1.

university-logo

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

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

∥r∥ = ∥b − Ax̃∥ ≈ 10^{−t}∥A∥ · ∥x̃∥.

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 ỹ be the approximate solution. Then

ỹ ≈ A^{−1}r = A^{−1}(b − Ax̃) = x − x̃

and

x ≈ x̃ + ỹ.
T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  58 / 1

university-logo

Moreover,

∥ỹ∥ ≈ ∥x − x̃∥ = ∥A^{−1}r∥
≤ ∥A^{−1}∥ · ∥r∥ ≈ ∥A^{−1}∥(10^{−t}∥A∥ · ∥x̃∥) = 10^{−t}∥x̃∥κ(A).

It implies that

κ(A) ≈ \frac{∥ỹ∥}{∥x̃∥}10^{t}.

Iterative refinement

In general, x̃ + ỹ is a more accurate approximation to the solution of Ax = b than x̃.

university-logo

Algorithm (Iterative refinement)

Given tolerance TOL, 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 ∥y∥_{∞} < TOL, then stop.

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

university-logo

Example

The linear system given by

[ \begin{matrix}3.3330 & 15920 & −10.333 \\ 2.2220 & 16.710 & 9.6120 \\ 1.5611 & 5.1791 & 1.6852 \end{matrix} ] [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ \begin{matrix}15913 \\ 28.544 \\ 8.4254 \end{matrix} ]

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

Using Gaussian elimination and five-digit rounding arithmetic leads successively to the augmented matrices

[ \begin{matrix}3.3330 & 15920 & −10.333 & 15913 \\ 0 & −10596 & 16.501 & −10580 \\ 0 & −7451.4 & 6.5250 & −7444.9 \end{matrix} ]

and

[ \begin{matrix}3.3330 & 15920 & −10.333 & 15913 \\ 0 & −10596 & 16.501 & −10580 \\ 0 & 0 & −5.0790 & −4.7000 \end{matrix} ] .

university-logo

The approximate solution is

x̃^{(1)} = [1.2001, 0.99991, 0.92538]^{T}.

The residual vector corresponding to x̃ is computed in double precision to be

r^{(1)} = b − Ax̃^{(1)}
= [ \begin{matrix}15913 \\ 28.544 \\ 8.4254 \end{matrix} ] − [ \begin{matrix}3.3330 & 15920 & −10.333 \\ 2.2220 & 16.710 & 9.6120 \\ 1.5611 & 5.1791 & 1.6852 \end{matrix} ] [ \begin{matrix}1.2001 \\ 0.99991 \\ 0.92538 \end{matrix} ]

= [ \begin{matrix}15913 \\ 28.544 \\ 8.4254 \end{matrix} ] − [ \begin{matrix}15913.00518 \\ 28.26987086 \\ 8.611560367 \end{matrix} ] = [ \begin{matrix}−0.00518 \\ 0.27412914 \\ −0.186160367 \end{matrix} ] .

Hence the solution of Ay = r^{(1)} to be

ỹ^{(1)} = [−0.20008, 8.9987 × 10^{−5}, 0.074607]^{T}

and the new approximate solution x^{(2)} is

x^{(2)} = x^{(1)} + ỹ^{(1)} = [1.0000, 1.0000, 0.99999]^{T}.
T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  62 / 1

university-logo

Using the suggested stopping technique for the algorithm, we compute r^{(2)} = b − Ax̃^{(2)} and solve the system Ay^{(2)} = r^{(2)}, which gives

ỹ^{(2)} = [1.5002 × 10^{−9}, 2.0951 × 10^{−10}, 1.0000 × 10^{−5}]^{T}.

Since

∥ỹ^{(2)}∥_{∞} ≤ 10^{−5},

we conclude that

x̃^{(3)} = x̃^{(2)} + ỹ^{(2)} = [1.0000, 1.0000, 1.0000]^{T}

is sufficiently accurate. In the linear system

Ax = b,

A and b can be represented exactly. Realistically, the matrix A and vector b will be perturbed by δA and δb, respectively, causing the linear system

(A + δA)x = b + δb

to be solved in place of Ax = b.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  63 / 1

university-logo

Theorem

Suppose A is nonsingular and

∥δA∥ < \frac{1}{∥A^{−1}∥}.

Then the solution x̃ of (A + δA)x̃ = b + δb approximates the solution x of Ax = b with the error estimate

\frac{∥x − x̃∥}{∥x∥} ≤ \frac{κ(A)}{1 − κ(A)(∥δA∥/∥A∥)} (\frac{∥δb∥}{∥b∥} + \frac{∥δA∥}{∥A∥}) .

If A is well-conditioned, then small changes in A and b produce correspondingly small changes in the solution x. If A is ill-conditioned, then small changes in A and b may produce large changes in x.

university-logo

The conjugate gradient method

Consider the linear systems

Ax = b

where A is large sparse and symmetric positive definite. Define the inner product notation

< x, y >= x^{T}y for any x, y ∈ ℝ^{n}.

Theorem

Let A be symmetric positive definite. Then x^{∗} is the solution of Ax = b if and only if x^{∗} minimizes

g(x) =< x, Ax > −2 < x, b > .

university-logo

Proof: (“ ⇒”) Rewrite g(x) as

g(x) = < x − x^{∗}, A(x − x^{∗}) > + < x, Ax^{∗} > + < x^{∗}, Ax >
− < x^{∗}, Ax^{∗} > −2 < x, b >
= < x − x^{∗}, A(x − x^{∗}) > − < x^{∗}, Ax^{∗} >
+2 < x, Ax^{∗} > −2 < x, b >
= < x − x^{∗}, A(x − x^{∗}) > − < x^{∗}, Ax^{∗} > +2 < x, Ax^{∗} − b > .

Suppose that x^{∗} is the solution of Ax = b, i.e., Ax^{∗} = b. Then

g(x) =< x − x^{∗}, A(x − x^{∗}) > − < x^{∗}, Ax^{∗} >

which minimum occurs at x = x^{∗}.

university-logo

(“⇐”) Fixed vectors x and v, for any α ∈ ℝ,

f (α) ≡ g(x + αv)
= < x + αv, Ax + αAv > −2 < x + αv, b >
= < x, Ax > +α < v, Ax > +α < x, Av > +α^{2} < v, Av >
−2 < x, b > −2α < v, b >
= < x, Ax > −2 < x, b > +2α < v, Ax > −2α < v, b > +α^{2} < v, Av >
= g(x) + 2α < v, Ax − b > +α^{2} < v, Av > .

Because f is a quadratic function of α and < v, Av > is positive, f has a minimal value when f′(α) = 0. Since

f′(α) = 2 < v, Ax − b > +2α < v, Av >,

the minimum occurs at

α̂ = −\frac{< v, Ax − b >}{< v, Av >} = \frac{< v, b − Ax >}{< v, Av >}

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  67 / 1

university-logo

and

g(x + α̂v) = f (α̂) = g(x) − 2\frac{< v, b − Ax >}{< v, Av >} < v, b − Ax >

+ (\frac{< v, b − Ax >}{< v, Av >})^{2} < v, Av >

= g(x) − \frac{< v, b − Ax >^{2}}{< v, Av >}.

So, for any nonzero vector v, we have

g(x + α̂v) < g(x) if < v, b − Ax >≠ 0  (6)

and

g(x + α̂v) = g(x) if < v, b − Ax >= 0.  (7)

Suppose that x^{∗} is a vector that minimizes g. Then

g(x^{∗} + α̂v) ≥ g(x^{∗}) for any v.  (8)

university-logo

From (??), (??) and (??), we have

< v, b − Ax^{∗} >= 0 for any v,

which implies that Ax^{∗} = b. Let

r = b − Ax.

Then

α = \frac{< v, b − Ax >}{< v, Av >} = \frac{< v, r >}{< v, Av >}.

If r ≠ 0 and if v and r are not orthogonal, then

g(x + αv) < g(x)

which implies that x + αv is closer to x^{∗} than is x.

university-logo

Let x^{(0)} be an initial approximation to x^{∗} and v^{(1)} ≠ 0 be an initial search direction. For k = 1, 2, 3, …, we compute

α_{k} = \frac{< v^{(k)}, b − Ax^{(k−1)} >}{< v^{(k)}, Av^{(k)} >},
x^{(k)} = x^{(k−1)} + α_{k}v^{(k)}

and choose a new search direction v^{(k+1)}. Question: How to choose {v^{(k)}} such that {x^{(k)}} converges rapidly to x^{∗}? Let Φ : ℝ^{n} → ℝ be a differential function on x. Then it holds

\frac{Φ(x + εp) − Φ(x)}{ε} = ∇Φ(x)^{T}p + O(ε).

The right hand side takes minimum at

p = −\frac{∇Φ(x)}{∥∇Φ(x)∥} (i.e., the largest descent)

for all p with ∥p∥ = 1 (neglect O(ε)).

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  70 / 1

university-logo

Denote x = [x_{1}, x_{2}, …, x_{n}]^{T}. Then

g(x) =< x, Ax > −2 < x, b >= ∑_{i=1}^{n} ∑_{j=1}^{n}a_{ij}x_{i}x_{j} − 2 ∑_{i=1}^{n}x_{i}b_{i}.

It follows that

\frac{∂g}{∂x_{k}} (x) = 2 ∑_{i=1}^{n}a_{ki}x_{i} − 2b_{k}, for k = 1, 2, …, n.

Therefore, the gradient of g is

∇g(x) = [\frac{∂g}{∂x_{1}}(x), \frac{∂g}{∂x_{2}}, ⋯, \frac{∂g}{∂x_{n}}(x)]^{T} = 2(Ax − b) = −2r.

university-logo

Steepest descent method (gradient method)

Given an initial x_{0} ≠ 0.

For k = 1, 2, … r_{k−1} = b − Ax_{k−1} If r_{k−1} = 0, then stop;

else α_{k} = \frac{r_{k−1}^{T}r_{k−1}}{r_{k−1}^{T}Ar_{k−1}}, x_{k} = x_{k−1} + α_{k}r_{k−1}.

End for

Theorem

If x_{k}, x_{k−1} are two approximations of the steepest descent method for solving Ax = b and λ_{1} ≥ λ_{2} ≥ ⋯ ≥ λ_{n} > 0 are the eigenvalues of A, then it holds:

∥x_{k} − x^{∗}∥_{A} ≤ (\frac{λ_{1} − λ_{n}}{λ_{1} + λ_{n}}) ∥x_{k−1} − x^{∗}∥_{A},

where ∥x∥_{A} = \sqrt{x^{T}Ax}. Thus the gradient method is convergent.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  72 / 1

university-logo

If the condition number of A (= λ_{1}/λ_{n}) is large, then \frac{λ_{1}−λ_{n}}{λ_{1}+λ_{n}} ≈ 1. The gradient method converges very slowly. Hence this method is not recommendable. It is favorable to choose that the search directions {v^{(i)}} as mutually A-conjugate, where A is symmetric positive definite.

Definition

Two vectors p and q are called A-conjugate (A-orthogonal), if p^{T}Aq = 0.

university-logo

Lemma

Let v_{1}, …, v_{n} ≠ 0 be pairwisely A-conjugate. Then they are linearly independent. Proof: From

0 = ∑_{j=1}^{n}c_{j}v_{j}

follows that

0 = (v_{k})^{T}A (∑_{j=1}^{n}c_{j}v_{j} ) = ∑_{j=1}^{n}c_{j}(v_{k})^{T}Av_{j} = c_{k}(v_{k})^{T}Av_{k},

so c_{k} = 0, for k = 1, …, n.

university-logo

Theorem

Let A be s.p.d. and v_{1}, …, v_{n} ∈ ℝ^{n}\{0} be pairwisely A-orthogonal. Give x_{0} and let r_{0} = b − Ax_{0}. For k = 1, …, n, let

α_{k} = \frac{< v_{k}, b − Ax_{k−1} >}{< v_{k}, Av_{k} >} and x_{k} = x_{k−1} + α_{k}v_{k}.

Then Ax_{n} = b and

< b − Ax_{k}, v_{j} >= 0, for each j = 1, 2, …, k − 1.

Proof: Since, for each k = 1, 2, …, n,

x_{k} = x_{k−1} + α_{k}v_{k},

we have

Ax_{n} = Ax_{n−1} + α_{n}Av_{n} = (Ax_{n−2} + α_{n−1}Av_{n−1}) + α_{n}Av_{n}
⋮
= Ax_{0} + α_{1}Av_{1} + α_{2}Av_{2} + ⋯ + α_{n}Av_{n}.
T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  75 / 1

university-logo

< Ax_{n} − b, v_{k} >
= < Ax_{0} − b, v_{k} > +α_{1} < Av_{1}, v_{k} > + ⋯ + α_{n} < Av_{n}, v_{k} >
= < Ax_{0} − b, v_{k} > +α_{1} < v_{1}, Av_{k} > + ⋯ + α_{n} < v_{n}, Av_{k} >
= < Ax_{0} − b, v_{k} > +α_{k} < v_{k}, Av_{k} >
= < Ax_{0} − b, v_{k} > +\frac{< v_{k}, b − Ax_{k−1} >}{< v_{k}, Av_{k} >} < v_{k}, Av_{k} >
= < Ax_{0} − b, v_{k} > + < v_{k}, b − Ax_{k−1} >
= < Ax_{0} − b, v_{k} >

+ < v_{k}, b − Ax_{0} + Ax_{0} − Ax_{1} + ⋯ − Ax_{k−2} + Ax_{k−2} − Ax_{k−1} >

= < Ax_{0} − b, v_{k} > + < v_{k}, b − Ax_{0} > + < v_{k}, Ax_{0} − Ax_{1} >
+ ⋯ + < v_{k}, Ax_{k−2} − Ax_{k−1} >
= < v_{k}, Ax_{0} − Ax_{1} > + ⋯ + < v_{k}, Ax_{k−2} − Ax_{k−1} > .

university-logo

For any i

x_{i} = x_{i−1} + α_{i}v_{i} and Ax_{i} = Ax_{i−1} + α_{i}Av_{i},

we have

Ax_{i−1} − Ax_{i} = −α_{i}Av_{i}.

Thus, for k = 1, …, n,

< Ax_{n} − b, v_{k} >

= −α_{1} < v_{k}, Av_{1} > − ⋯ − α_{k−1} < v_{k}, Av_{k−1} >= 0

which implies that Ax_{n} = b. Suppose that

< r_{k−1}, v_{j} >= 0 for j = 1, 2, …, k − 1.  (9)

By the result

r_{k} = b − Ax_{k} = b − A(x_{k−1} + α_{k}v_{k}) = r_{k−1} − α_{k}Av_{k}

university-logo

it follows that

< r_{k}, v_{k} > = < r_{k−1}, v_{k} > −α_{k} < Av_{k}, v_{k} >
= < r_{k−1}, v_{k} > −\frac{< v_{k}, b − Ax_{k−1} >}{< v_{k}, Av_{k} >} < Av_{k}, v_{k} >
= 0.

From assumption (??)and A-orthogonality, for j = 1, …, k − 1

< r_{k}, v_{j} >=< r_{k−1}, v_{j} > −α_{k} < Av_{k}, v_{j} >= 0

which is completed the proof by the mathematic induction. Method of conjugate directions: Let A be s.p.d., b, x_{0} ∈ ℝ^{n}. Given v_{1}, …, v_{n} ∈ ℝ^{n}\{0} pairwisely A-orthogonal.

r_{0} = b − Ax_{0}, For k = 1, …, n,

α_{k} = \frac{<v_{k},r_{k−1}>}{<v_{k},Av_{k}>}, x_{k} = x_{k−1} + α_{k}v_{k},
r_{k} = r_{k−1} − α_{k}Av_{k} = b − Ax_{k}.

End For

university-logo

Practical Implementation In k-th step a direction v_{k} which is A-orthogonal to v_{1}, …, v_{k−1} must be determined. It allows for orthogonalization of r_{k} against v_{1}, …, v_{k}. Let r_{k} ≠ 0, g(x) decreases strictly in the direction −r_{k}. For ε > 0 small, we have g(x_{k} − εr_{k}) < g(x_{k}).

If r_{k−1} = b − Ax_{k−1} ≠ 0, then we use r_{k−1} to generate v_{k} by

v_{k} = r_{k−1} + β_{k−1}v_{k−1}.  (10)

Choose β_{k−1} such that

0 = < v_{k−1}, Av_{k} >=< v_{k−1}, Ar_{k−1} + β_{k−1}Av_{k−1} >
= < v_{k−1}, Ar_{k−1} > +β_{k−1} < v_{k−1}, Av_{k−1} > .

university-logo

That is

β_{k−1} = −\frac{< v_{k−1}, Ar_{k−1} >}{< v_{k−1}, Av_{k−1} >}.  (11)

Theorem

Let v_{k} and β_{k−1} be defined in (??) and (??), respectively. Then r_{0}, …, r_{k−1} are mutually orthogonal and

< v_{k}, Av_{i} >= 0, for i = 1, 2, …, k − 1.

That is {v_{1}, …, v_{k}} is an A-orthogonal set.

Having chosen v_{k}, we compute

α_{k} = \frac{< v_{k}, r_{k−1} >}{< v_{k}, Av_{k} >} = \frac{< r_{k−1} + β_{k−1}v_{k−1}, r_{k−1} >}{< v_{k}, Av_{k} >}
= \frac{< r_{k−1}, r_{k−1} >}{< v_{k}, Av_{k} >} + β_{k−1} \frac{< v_{k−1}, r_{k−1} >}{< v_{k}, Av_{k} >}
= \frac{< r_{k−1}, r_{k−1} >}{< v_{k}, Av_{k} >}.  (12)

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  80 / 1

university-logo

Since

r_{k} = r_{k−1} − α_{k}Av_{k},

we have

< r_{k}, r_{k} >=< r_{k−1}, r_{k} > −α_{k} < Av_{k}, r_{k} >= −α_{k} < r_{k}, Av_{k} > . Further, from (??),

< r_{k−1}, r_{k−1} >= α_{k} < v_{k}, Av_{k} >,

so

β_{k} = −\frac{< v_{k}, Ar_{k} >}{< v_{k}, Av_{k} >} = −\frac{< r_{k}, Av_{k} >}{< v_{k}, Av_{k} >}
= \frac{(1/α_{k}) < r_{k}, r_{k} >}{(1/α_{k}) < r_{k−1}, r_{k−1} >} = \frac{< r_{k}, r_{k} >}{< r_{k−1}, r_{k−1} >}.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  81 / 1

university-logo

Algorithm (Conjugate Gradient method (CG-method))

Let A be s.p.d., b ∈ ℝ^{n}, choose x_{0} ∈ ℝ^{n}, r_{0} = b − Ax_{0} = v_{0}. If r_{0} = 0, then N = 0 stop, otherwise for k = 0, 1, …

(a). α_{k} = \frac{<r_{k},r_{k}>}{<v_{k},Av_{k}>}, (b). x_{k+1} = x_{k} + α_{k}v_{k}, (c). r_{k+1} = r_{k} − α_{k}Av_{k}, (d). If r_{k+1} = 0, let N = k + 1, stop. (e). β_{k} = \frac{<r_{k+1},r_{k+1}>}{<r_{k},r_{k}>}, (f). v_{k+1} = r_{k+1} + β_{k}v_{k}.

Theoretically, the exact solution is obtained in n steps. If A is well-conditioned, then approximate solution is obtained in about \sqrt{n} steps. If A is ill-conditioned, then the number of iterations may be greater than n.

university-logo

Select a nonsingular matrix C so that

Ã = C^{−1}AC^{−T}

is better conditioned. Consider the linear system

Ãx̃ = b̃,

where

x̃ = C^{T}x and b̃ = C^{−1}b.

Then

Ãx̃ = (C^{−1}AC^{−T})(C^{T}x) = C^{−1}Ax.

Thus,

Ax = b ⇔ Ãx̃ = b̃ and x = C^{−T} x̃.

university-logo

x̃_{k} = C^{T}x_{k},

we have

r̃_{k} = b̃ − Ãx̃_{k} = C^{−1}b − (C^{−1}AC^{−T})C^{T}x_{k}
= C^{−1}(b − Ax_{k}) = C^{−1}r_{k}.

Let

ṽ_{k} = C^{T}v_{k} and w_{k} = C^{−1}r_{k}.

Then

β̃_{k} = \frac{< r̃_{k}, r̃_{k} >}{< r̃_{k−1}, r̃_{k−1} >} = \frac{< C^{−1}r_{k}, C^{−1}r_{k} >}{< C^{−1}r_{k−1}, C^{−1}r_{k−1} >}
= \frac{< w_{k}, w_{k} >}{< w_{k−1}, w_{k−1} >}.

T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  84 / 1

university-logo

Thus,

α̃_{k} = \frac{< r̃_{k−1}, r̃_{k−1} >}{< ṽ_{k}, Ãṽ_{k} >} = \frac{< C^{−1}r_{k−1}, C^{−1}r_{k−1} >}{< C^{T}v_{k}, C^{−1}AC^{−T}C^{T}v_{k} >}
= \frac{< w_{k−1}, w_{k−1} >}{< C^{T}v_{k}, C^{−1}Av_{k} >}

and, since

< C^{T}v_{k}, C^{−1}Av_{k} > = [v_{k}]^{T} CC^{−1}Av_{k} = [v_{k}]^{T} Av_{k}
= < v_{k}, Av_{k} >,

we have

α̃_{k} = \frac{< w_{k−1}, w_{k−1} >}{< v_{k}, Av_{k} >}.

Further,

x̃_{k} = x̃_{k−1} + α̃_{k} ṽ_{k}, so C^{T}x_{k} = C^{T}x_{k−1} + α̃_{k}C^{T}v_{k}

and

x_{k} = x_{k−1} + α̃_{k}v_{k}.
T.M. Hwang (Taiwan Normal Univ. )  Iterative Tech. in Matrix Algebra  February 16, 2011  85 / 1

university-logo

Continuing,

r̃_{k} = r̃_{k−1} − α̃_{k} Ãṽ_{k},

so

C^{−1}r_{k} = C^{−1}r_{k−1} − α̃_{k}C^{−1}AC^{−T}C^{T}v_{k}

and

r_{k} = r_{k−1} − α̃_{k}Av_{k}.

Finally,

ṽ_{k+1} = r̃_{k} + β̃_{k} ṽ_{k} and C^{T}v_{k+1} = C^{−1}r_{k} + β̃_{k}C^{T}v_{k},

so

v_{k+1} = C^{−T}C^{−1}r_{k} + β̃_{k}v_{k} = C^{−T}w_{k} + β̃_{k}v_{k}.

university-logo

Algorithm (Preconditioned CG-method (PCG-method))

Choose C and x_{0}. Set r_{0} = b − Ax_{0}, solve Cw_{0} = r_{0} and C^{T}v_{1} = w_{0}. If r_{0} = 0, then N = 0 stop, otherwise for k = 1, 2, …

(a). α_{k} =< w_{k−1}, w_{k−1} > / < v_{k}, Av_{k} >, (b). x_{k} = x_{k−1} + α_{k}v_{k}, (c). r_{k} = r_{k−1} − α_{k}Av_{k}, (d). If r_{k} = 0, let N = k + 1, stop.

Otherwise, solve Cw_{k} = r_{k} and C^{T}z_{k} = w_{k}, (e). β_{k} =< w_{k}, w_{k} > / < w_{k−1}, w_{k−1} >, (f). v_{k+1} = z_{k} + β_{k}v_{k}.