Example slides from a few Prosper-based lectures

Matthias Heil & Andrew L. Hazel

                                .
  M.Heil@maths.man.ac.uk     &     ahazel@maths.man.ac.uk
              http://www.maths.man.ac.uk/˜ mheil

Department of Mathematics University of Manchester

What are Krylov subspaces and why do they matter?

We want to solve Ax = b

What are Krylov subspaces and why do they matter?

We want to solve

Ax = b

• Consider the elementary splitting: Ix = (I − A)x + b

What are Krylov subspaces and why do they matter?

We want to solve

Ax = b

• Consider the elementary splitting: Ix = (I − A)x + b

• Write this as a fixpoint iteration: x_{n+1} = (I − A)x_{n} + b

What are Krylov subspaces and why do they matter?

We want to solve

Ax = b

• Consider the elementary splitting: Ix = (I − A)x + b

• Write this as a fixpoint iteration: x_{n+1} = (I − A)x_{n} + b

• Rewrite as

x_{n+1} = x_{n} + (b − Ax_{n})

What are Krylov subspaces and why do they matter?

We want to solve

Ax = b

• Consider the elementary splitting: Ix = (I − A)x + b

• Write this as a fixpoint iteration: x_{n+1} = (I − A)x_{n} + b

• Rewrite as

x_{n+1} = x_{n} + (b − Ax_{n})

• or:

x_{n+1} = x_{n} + r_{n}

where r_{n} = b − Ax_{n} is the residual.

What are Krylov subspaces and why do they matter?

We want to solve

Ax = b

• Consider the elementary splitting: Ix = (I − A)x + b

• Write this as a fixpoint iteration: x_{n+1} = (I − A)x_{n} + b

• Rewrite as

x_{n+1} = x_{n} + (b − Ax_{n})

• or:

x_{n+1} = x_{n} + r_{n}

where r_{n} = b − Ax_{n} is the residual.

• For simplicity assume x_{0} = 0, then:

x_{n+1} = ∑_{i=0}^{n} r_{i} where r_{0} = b.

Matthias Heil & Andrew L. Hazel, Department of Mathematics, University of Manchester, UK http://www.maths.man.ac.uk/˜ mheil – p.2/10

What are Krylov subspaces and why do they matter? (2)

x_{n+1} = (I−A)x_{n} +b (1)  x_{n+1} = x_{n} +r_{n} (2)  x_{n+1} = ∑_{i=0}^{n} r_{i} (3)

x_{n+1} = (I−A)x_{n} +b (1)  x_{n+1} = x_{n} +r_{n} (2)  x_{n+1} = ∑_{i=0}^{n} r_{i} (3)

• Multiply (2) by −A and add b:

b − Ax_{n+1} = b − Ax_{n} − Ar_{n}

x_{n+1} = (I−A)x_{n} +b (1)  x_{n+1} = x_{n} +r_{n} (2)  x_{n+1} = ∑_{i=0}^{n} r_{i} (3)

• Multiply (2) by −A and add b:

b − Ax_{n+1} = b − Ax_{n} − Ar_{n}

r_{n+1} = r_{n} − Ar_{n}

x_{n+1} = (I−A)x_{n} +b (1)  x_{n+1} = x_{n} +r_{n} (2)  x_{n+1} = ∑_{i=0}^{n} r_{i} (3)

• Multiply (2) by −A and add b:

b − Ax_{n+1} = b − Ax_{n} − Ar_{n}

r_{n+1} = r_{n} − Ar_{n}

r_{n+1} = (I − A)r_{n}

What are Krylov subspaces and why do they matter? (2)

x_{n+1} = (I−A)x_{n} +b (1)  x_{n+1} = x_{n} +r_{n} (2)  x_{n+1} = ∑_{i=0}^{n} r_{i} (3)

• Multiply (2) by −A and add b:

b − Ax_{n+1} = b − Ax_{n} − Ar_{n}

r_{n+1} = r_{n} − Ar_{n}
r_{n+1} = (I − A)r_{n}
r_{n+1} = (I − A)^{n+1}r_{0}

Note that this shows that the iteration only converges if ρ(I − A) < 1.

Matthias Heil & Andrew L. Hazel, Department of Mathematics, University of Manchester, UK http://www.maths.man.ac.uk/˜ mheil – p.3/10

x_{n+1} = (I−A)x_{n} +b (1)  x_{n+1} = x_{n} +r_{n} (2)  x_{n+1} = ∑_{i=0}^{n} r_{i} (3)

• Multiply (2) by −A and add b:

b − Ax_{n+1} = b − Ax_{n} − Ar_{n}
r_{n+1} = r_{n} − Ar_{n}

r_{n+1} = (I − A)r_{n}

r_{i} = (I − A)^{i}r_{0} (4)

x_{n+1} = (I−A)x_{n} +b (1)  x_{n+1} = x_{n} +r_{n} (2)  x_{n+1} = ∑_{i=0}^{n} r_{i} (3)

• Multiply (2) by −A and add b:

b − Ax_{n+1} = b − Ax_{n} − Ar_{n}
r_{n+1} = r_{n} − Ar_{n}

r_{n+1} = (I − A)r_{n}

r_{i} = (I − A)^{i}r_{0} (4)

• Insert (4) into (3):

x_{n+1} = ∑_{i=0}^{n}(I − A)^{i}r_{0} = ∑_{i=0}^{n}(I − A)^{i}b

What are Krylov subspaces and why do they matter? (3)

• So while we’re carrying out the iteration

x_{n} = (I − A)x_{n−1} + b

the method generates the iterates

x_{n} = ∑_{i=0}^{n−1}(I − A)^{i}b

What are Krylov subspaces and why do they matter? (3)

• So while we’re carrying out the iteration

x_{n} = (I − A)x_{n−1} + b

the method generates the iterates

x_{n} = ∑_{i=0}^{n−1}(I − A)^{i}b

• Expanding this gives

x_{n} = ∑_{i=0}^{n−1} α_{i}A^{i}b

for some coefficients α_{i}.

What are Krylov subspaces and why do they matter? (3)

• So while we’re carrying out the iteration

x_{n} = (I − A)x_{n−1} + b

the method generates the iterates

x_{n} = ∑_{i=0}^{n−1}(I − A)^{i}b

• Expanding this gives

x_{n} = ∑_{i=0}^{n−1} α_{i}A^{i}b

for some coefficients α_{i}.

• In other words:

x_{n} ∈ span(b, Ab, A^{2}b, ..., A^{n−1}b) =: K_{n}(A, b)

where K_{n}(A, b) is the Krylov subspace.

What are Krylov subspaces and why do they matter? (4)

• Our (well, Richardson’s...) method produces:

x_{n} = ∑_{i=0}^{n−1} α_{i}A^{i}b ∈ K_{n}(A, b)

for some coefficients α_{i}.

x_{n} = ∑_{i=0}^{n−1} α_{i}A^{i}b ∈ K_{n}(A, b)

for some coefficients α_{i}.

• Are there any ‘better’ choices for the coefficients α_{i}?

x_{n} = ∑_{i=0}^{n−1} α_{i}A^{i}b ∈ K_{n}(A, b)

for some coefficients α_{i}.

• Are there any ‘better’ choices for the coefficients α_{i}?

• Yes, but it depends how you define ‘better’, or indeed ‘best’ as we’re trying to find the ‘best’ approximation in a given Krylov subspace.

x_{n} = ∑_{i=0}^{n−1} α_{i}A^{i}b ∈ K_{n}(A, b)

for some coefficients α_{i}.

• Are there any ‘better’ choices for the coefficients α_{i}?

• Yes, but it depends how you define ‘better’, or indeed ‘best’ as we’re trying to find the ‘best’ approximation in a given Krylov subspace.

• Various approaches and the resulting schemes:

• Minimise the residual (MINRES, GMRES)

• Minimise the error (GMERR)

• Make the residual orthogonal to the Krylov subspace (Ritz/Galerkin: CG, CGS)

• ...

x_{n} = ∑_{i=0}^{n−1} α_{i}A^{i}b ∈ K_{n}(A, b)

for some coefficients α_{i}.

• Are there any ‘better’ choices for the coefficients α_{i}?

• Yes, but it depends how you define ‘better’, or indeed ‘best’ as we’re trying to find the ‘best’ approximation in a given Krylov subspace.

• Various approaches and the resulting schemes:

• Minimise the residual (MINRES, GMRES)

• Minimise the error (GMERR)

• Make the residual orthogonal to the Krylov subspace (Ritz/Galerkin: CG, CGS)

• ...

• An important technical question is how to construct the basis for the Krylov subspace if we abandon our simple iterative scheme.

How quickly do the methods converge?

• Most iterative methods perform poorly if applied directly to the original problem Ax = b.

How quickly do the methods converge?

• Most iterative methods perform poorly if applied directly to the original problem Ax = b.

• Idea: Solve a ‘preconditioned’ problem

PAx = Pb

where P is chosen such that the iterative method applied to PA converges more quickly.

How quickly do the methods converge? (2)

Preconditioned problem

PAx = Pb

How quickly do the methods converge? (2)

Preconditioned problem

PAx = Pb

• How do we choose/construct P?

How quickly do the methods converge? (2)

Preconditioned problem

PAx = Pb

• How do we choose/construct P?

• Simple idea: Use a (cheap!) approximation to the inverse. P ≈ A^{−1}.

How quickly do the methods converge? (2)

Preconditioned problem

PAx = Pb

• How do we choose/construct P?

• Simple idea: Use a (cheap!) approximation to the inverse. P ≈ A^{−1}.

• Pointless in its exact form!

How quickly do the methods converge? (2)

Preconditioned problem

PAx = Pb

• How do we choose/construct P?

• Simple idea: Use a (cheap!) approximation to the inverse.

P ≈ A^{−1}.

• Pointless in its exact form!

• Some popular choices:

• Diagonal preconditioner (also in block form):

P = (diagA)^{−1}

• Incomplete LU factorisation:

P = (L_{inc}U_{inc})^{−1}

• In general, the choice of preconditioner is very problem dependent.

Assessing the convergence: The minimal polynomial

• We have already seen that the convergence is related to the eigenvalues of the iteration matrix; ρ(I − A) < 1 for the Richardson iteration.

Assessing the convergence: The minimal polynomial

• How exactly are the eigenvalues of the matrix related to the convergence of Krylov subspace methods?

Assessing the convergence: The minimal polynomial

• How exactly are the eigenvalues of the matrix related to the convergence of Krylov subspace methods?

• Definition: The minimal polynomial of a N × N nonsingular (diagonalisable) matrix

B with m ≤ N distinct eigenvalues λ_{j} is given by q(t) = Π_{j=1}^{m}(t − λ_{j}).

Assessing the convergence: The minimal polynomial

• How exactly are the eigenvalues of the matrix related to the convergence of Krylov subspace methods?

• Definition: The minimal polynomial of a N × N nonsingular (diagonalisable) matrix

B with m ≤ N distinct eigenvalues λ_{j} is given by

q(t) = Π_{j=1}^{m}(t − λ_{j}).

• The minimal polynomial has the (obvious) property that q(B) = Π_{j=1}^{m}(B − Iλ_{j}) = 0.

Assessing the convergence: The minimal polynomial

• How exactly are the eigenvalues of the matrix related to the convergence of Krylov subspace methods?

• Definition: The minimal polynomial of a N × N nonsingular (diagonalisable) matrix

B with m ≤ N distinct eigenvalues λ_{j} is given by

q(t) = Π_{j=1}^{m}(t − λ_{j}).

• The minimal polynomial has the (obvious) property that q(B) = Π_{j=1}^{m}(B − Iλ_{j}) = 0.

• Expanding out the polynomial gives

q(t) = Π_{j=1}^{m}(t − λ_{j}) = ∑_{j=0}^{m} β_{j}t^{j}

for some coefficients β_{j}.

Assessing the convergence: The minimal polynomial

• How exactly are the eigenvalues of the matrix related to the convergence of Krylov subspace methods?

• Definition: The minimal polynomial of a N × N nonsingular (diagonalisable) matrix

B with m ≤ N distinct eigenvalues λ_{j} is given by

q(t) = Π_{j=1}^{m}(t − λ_{j}).

• The minimal polynomial has the (obvious) property that q(B) = Π_{j=1}^{m}(B − Iλ_{j}) = 0.

• Expanding out the polynomial gives

q(t) = Π_{j=1}^{m}(t − λ_{j}) = ∑_{j=0}^{m} β_{j}t^{j}

for some coefficients β_{j}.

• Applying this to q(B) = 0 gives:

0 = β_{0}I + β_{1}B + ... + β_{N} B^{m}

Matthias Heil & Andrew L. Hazel, Department of Mathematics, University of Manchester, UK http://www.maths.man.ac.uk/˜ mheil – p.8/10

Assessing the convergence: The minimal polynomial (2)

0 = β_{0}I + β_{1}B + ... + β_{m}B^{m}

0 = β_{0}I + β_{1}B + ... + β_{m}B^{m}

• Multiply through by B^{−1}:

B^{−1} = − \frac{1}{β_{0}} ∑_{j=0}^{m−1} β_{j}B^{j}.

0 = β_{0}I + β_{1}B + ... + β_{m}B^{m}

• Multiply through by B^{−1}:

B^{−1} = − \frac{1}{β_{0}} ∑_{j=0}^{m−1} β_{j}B^{j}.

• So the solution of Bx = b can be represented as

x = B^{−1}b = − \frac{1}{β_{0}} ∑_{j=0}^{m−1} β_{j}B^{j}b.

0 = β_{0}I + β_{1}B + ... + β_{m}B^{m}

• Multiply through by B^{−1}:

B^{−1} = − \frac{1}{β_{0}} ∑_{j=0}^{m−1} β_{j}B^{j}.

• So the solution of Bx = b can be represented as

x = B^{−1}b = − \frac{1}{β_{0}} ∑_{j=0}^{m−1} β_{j}B^{j}b.

• This shows that

x ∈ K_{m}(B, b).

Assessing the convergence: The minimal polynomial (3)

• This characterises the solution to Bx = b as a member of a Krylov subspace

x ∈ K_{m}(B, b),

where m is the number of distinct eigenvalues of the matrix B.

• This characterises the solution to Bx = b as a member of a Krylov subspace

x ∈ K_{m}(B, b),

where m is the number of distinct eigenvalues of the matrix B.

• The fewer distinct eigenvalues the (preconditioned) matrix B = PA has, the lower the degree of the minimal polynomial and the smaller the Krylov subspace, the solution x lives in.

• This characterises the solution to Bx = b as a member of a Krylov subspace

x ∈ K_{m}(B, b),

where m is the number of distinct eigenvalues of the matrix B.

• The fewer distinct eigenvalues the (preconditioned) matrix B = PA has, the lower the degree of the minimal polynomial and the smaller the Krylov subspace, the solution x lives in.

• Recall that during the n-th iteration, a Krylov subspace method will determine the ‘best’ approximate solution in the Krylov subspace K_{n}(B, b).

• This characterises the solution to Bx = b as a member of a Krylov subspace

x ∈ K_{m}(B, b),

where m is the number of distinct eigenvalues of the matrix B.

• The fewer distinct eigenvalues the (preconditioned) matrix B = PA has, the lower the degree of the minimal polynomial and the smaller the Krylov subspace, the solution x lives in.

• Recall that during the n-th iteration, a Krylov subspace method will determine the ‘best’ approximate solution in the Krylov subspace K_{n}(B, b).

• ⟹ The solution will be found in at most m iterations!

• This characterises the solution to Bx = b as a member of a Krylov subspace

x ∈ K_{m}(B, b),

where m is the number of distinct eigenvalues of the matrix B.

• The fewer distinct eigenvalues the (preconditioned) matrix B = PA has, the lower the degree of the minimal polynomial and the smaller the Krylov subspace, the solution x lives in.

• Recall that during the n-th iteration, a Krylov subspace method will determine the ‘best’ approximate solution in the Krylov subspace K_{n}(B, b).

• ⟹ The solution will be found in at most m iterations!

• This provides an alternative characterisation of a good preconditioner:

‘A sufficient condition for a good preconditioner is that the preconditioned matrix B = PA has a low degree minimal polynomial.’ (Murphy et al. 2000)