GMRES in FVLib

S. Clain ^{a}

January 13, 2012

1 The theory

1.1 Projection method and Krylov Spaces

Let consider a linear problem Ax = b in ℝ^{n} with A an invertible matrix. For k ≤ n, we consider two sub vectorial spaces K_{k}, L_{k} of ℝ^{n} of dimension k. Let x_{0} ∈ ℝ^{n}, the projection method consists in providing x_{k} ∈ ℝ^{n} such that we both have

x_{k} ∈ x_{0} + K_{k}, Ax_{k} − b ⊥ L_{k}.

In other words, we seek u_{k} ∈ K_{k} such that

< A(x_{0} + u_{k}) − b, v >=< Au_{k} − r_{0}, v >= 0,  ∀v ∈ L_{k}  (1)

with r_{0} = Ax_{0} − b. If A is invertible, we have existence and uniqueness of the solution.

We introduce the Krylov space K_{k}(A, r_{0},) given by

K_{k}(A, r_{0},) = Span(r_{0}, Ar_{0}, ..., A^{k−1}r_{0})

and set K_{k} = K_{k}(A, r_{0}), L_{k} = AK_{k}(A, r_{0},) = Span(Ar_{0}, ..., A^{k−1}r_{0}, A^{k}r_{0}).

proposition x_{k} is the solution of (1) if and only if x_{k} minimizes the functional

u ∈ K_{k} → E(x_{0} + u) =< Au − r_{0}, Au − r_{0} >  (2)

with corresponds to the minimization of the residu on K_{k}. Moreover we have E(x_{0} + u_{k+1}) ≤ E(x_{0} + u_{k}), thus the residual always decreases. ■

If u_{k} ∈ K_{k} is the minimizer, then we have using definition of space L_{k}

< Au_{k} − r_{0}, Av >= 0 ∀v ∈ K_{k}.  (3)

We have the following result.

proposition If the Krylov space is invariant by the composition with A

AK_{k}(A, r_{0}) = K_{k}(A, r_{0})  (4)

then x_{k} = x_{0} + u_{k} is the solution, namely Ax_{k} − b = 0. ■

One can prove that there exists at least a k ≤ n such that (4) is achieved but k coul be very large. In the iterative method, we only require that the residual part E(x_{0} + u_{k}) is small enough.

1.2 Orthogonal basis and the Arnoldi Algorithm

When property (4) is not achieved, vectors r_{0}, Ar_{0}, ..., A^{k−1}r_{0} define a basis but the major inconvenient is that the vector are merely parallel (remember that A^{k}r_{0} converges to the first eigenvector). In order to provide a more suitable basis we build an orthogonal basis using the Arnoldi argorithm.

We start setting q_{1} = r_{0}/|r_{0}|. Assume now that q_{1}, q_{2}, .., q_{k} is a base of K_{k}(A, r_{0}). If property (4) is not achieved then we can prove that Aq_{k} ∉ K_{k}(A, r_{0}) and we write

Aq_{k} = ∑_{i=1}^{k+1} h_{jk}q_{j}

where q_{k+1} is a normalized vector, we have to determine. Using the orthogonal condition, we compute h_{jk} =< Aq_{k}, q_{j} >, j = 1, ..., k. Then we set

p_{k+1} = Aq_{k} − ∑_{i=1}^{k} h_{jk}q_{j},  h_{k+1,k} = |p_{k+1}|,  q_{k+1} = \frac{p_{k+1}}{h_{k+1,k}} .

From a matricial point of view, we have AQ_{k} = Q_{k+1}H(k + 1, k) where Q_{k} is a n × k Stiefel matrix and H(k + 1, k) is a k + 1 × k Hessenberg matrix

Q = [q_{1}, ..., q_{k}], H(k + 1, k) = [ \begin{matrix}h_{11} & … & … & … & h_{1k} \\ h_{21} & h_{22} & … & … & ⋮ \\ 0 & h_{32} & ... & … & ⋮ \\ ⋮ & … & ... & ... & ⋮ \\ 0 & … & 0 & h_{k,k.1} & h_{k,k} \\ 0 & … & … & 0 & h_{k+1,k} \end{matrix} ] .

Note that h_{k+1,k} = 0 is equivalent to property (4). Hence we can proceed as long as h_{k+1,k} > 0.

1.3 An equivalent minimizing condition

Vector u_{k} is the unique element of K_{k} wich minimizes

E(x_{0} + u) =< Au − r_{0}, Au − r_{0} >= ∥Au − r_{0}∥^{2}.

Using the orthogonal basis Q of K_{k}, we can write any u ∈ K_{k} as

u = ∑_{i=1}^{k} y_{i}^{k}q_{i}, y^{k} ∈ ℝ^{k}

Hence Au = AQ_{k}y^{k}. Note that u has a fix dimension while dimension of vector y^{k} changes with k. Moreover let us remark that r_{0} = |r_{0}|Q_{k+1}e^{k+1} = (1, 0, ..., 0) the first canonical vector of ℝ^{k+1}. We then have

E(x_{0} + u) = ∥Au − r_{0}∥^{2}
= ∥AQ_{k}y^{k} − |r_{0}|Q_{k+1}e^{k+1}∥^{2}
= ∥Q_{k+1}H(k + 1, k)y^{k} − |r_{0}|Q_{k+1}e^{k+1}∥^{2}
= ∥Q_{k+1}{H(k + 1, k)y^{k} − |r_{0}|e^{k+1}}∥^{2}

Since Q_{k+1} is unitary, we deduce

E(x_{0} + u) = Ẽ_{k}(y^{k}) = ∥H(k + 1, k)y^{k} − |r_{0}|e^{k+1}∥^{2}.

So we reduce the problem of minimization in ℝ^{k} of the overdetermined system H(k + 1, k)y^{k} = |r_{0}|e^{k+1}.

1.4 Orthogonalization of the Hessenberg matrix

To solve the minimizination problem, we use a QR transformation of the H(k + 1, k) matrix using Givens rotations.

Step 1 One have H(2, 1) = [h_{11}, h_{21}]^{T}. We apply a Givens matrix (unitary matrix) of the form

G_{1} = [ \begin{matrix}cos(θ) & sin(θ) \\ − sin(θ) & cos(θ) \end{matrix} ] .

If we choose

cos(θ) = \frac{h_{11}}{\sqrt{h_{11}^{2} + h_{21}^{2}}} , sin(θ) = \frac{h_{21}}{\sqrt{h_{11}^{2} + h_{21}^{2}}}

We find

U(2, 1) = G_{1}H(2, 1) = [ \begin{matrix}√ & \overline{h_{11}^{2} + h_{21}^{2}} \\ 0 \end{matrix} ] , c^{1} = |r_{0}|G_{1}e_{2}^{1} = |r_{0}| [ \frac{h_{11}}{\sqrt{h_{11}^{2} + h_{21}^{2}}} ; \frac{−h_{21}}{\sqrt{h_{11}^{2} + h_{21}^{2}}} ] = [ \begin{matrix}c_{1}^{1} \\ c_{2}^{1} \end{matrix} ] .

Step k We first note that

H(k + 1, k) = [ H(k, k − 1) h_{1k} ; ⋮ ; 0 h_{k+1,k} ] .

Now assume that we have determined the Givens matrizes G_{1},...,G_{k−1} such that we have a k × k − 1 upper triangular matrix U(k, k − 1) and a vector c^{k−1} ∈ ℝ^{k} and

∥H(k, k − 1)y^{k−1} − |r_{0}|e^{k−1}∥ = ∥U(k − 1, k)y^{k−1} − |r_{0}|c^{k−1}∥

Problem of minimizing ∥H(k + 1, k)y^{k} − |r_{0}|e^{k}∥ is equivalent to minimize ∥Ũ(k + 1, k)y^{k} − |r_{0}|c̃^{k}∥ with

Ũ(k + 1, k) = [ U(k, k − 1) h_{1k} ; ⋮ ; 0 h_{k+1,k} ] , c̃^{k} = [ c_{1}^{k−1} ; ⋮ ; c_{k}^{k−1} ; 0 ] .

We consider the k + 1 × k + 1 Givens matrix

G_{k} = [ \begin{matrix}Id & 0 & 0 \\ 0 & cs_{k} & sn_{k} \\ 0 & −sn_{k} & cs_{k} \end{matrix} ] ,

if one chooses

cs_{k} = \frac{h_{kk}}{\sqrt{h_{kk}^{2} + h_{k+1,k}^{2}}} , sn_{k} = \frac{h_{k+1,k}}{\sqrt{h_{kk}^{2} + h_{k+1,k}^{2}}}

then we have

U(k + 1, k) = G_{k} [ h̃_{1k} ; U(k, k − 1) ⋮ ; h̃_{k−1,k} ; 0  h̃_{k,k} ; 0  h_{k+1,k} ] , c_{k} = G_{k} [ c_{1}^{k−1} ; ⋮ ; c_{k−1}^{k−1} ; c_{k}^{k−1} ; 0 ] ,

where we have compute the successive products h̃^{k} = G_{k−1}G_{k−2}...G_{2}G_{1}h^{k}. Note that we

We get the solution y^{k} using a backward substitution of system

U(k + 1, k)y^{k} = c^{k}

starting at line k since the last line (k+1) is degenerated. We evaluate the minimum with

x_{k} = x_{0} + ∑_{i=1}^{y} y_{i}^{k}q_{i}.