Krylov subspaces - Fundamentals of Numerical Computation

The power and inverse iterations have a flaw that seems obvious once it is pointed out. Given a seed vector $\mathbf{u}$ , they produce a sequence of vectors $\mathbf{u}_1,\mathbf{u}_2,\ldots$ that are scalar multiples of $\mathbf{u},\mathbf{A}\mathbf{u},\mathbf{A}^{2}\mathbf{u},\ldots$ , but only the most recent vector is used to produce an eigenvector estimate.

It stands to reason that we could do no worse, and perhaps much better, if we searched among all linear combinations of the vectors seen in the past. In other words, we seek a solution in the range (column space) of the matrix

\mathbf{K}_m = \begin{bmatrix} \mathbf{u} & \mathbf{A}\mathbf{u} & \mathbf{A}^{2} \mathbf{u} & \cdots & \mathbf{A}^{m-1} \mathbf{u} \end{bmatrix}.

(8.4.1)

In general, we expect that the dimension of the Krylov^[1] subspace $\mathcal{K}_m$ , which is the rank of $\mathbf{K}_m$ , equals $m$ , though it may be smaller.

8.4.1Properties¶

As we have seen with the power iteration, part of the appeal of the Krylov matrix is that it can be generated in a way that fully exploits the sparsity of $\mathbf{A}$ , simply through repeated matrix-vector multiplication. Furthermore, we have some important mathematical properties.

8.4.2Dimension reduction¶

The problems $\mathbf{A}\mathbf{x}=\mathbf{b}$ and $\mathbf{A}\mathbf{x}=\lambda\mathbf{x}$ are posed in a very high-dimensional space $\mathbb{R}^n$ or $\mathbb{C}^n$ . One way to approximate them is to replace the full $n$ -dimensional space with a much lower-dimensional $\mathcal{K}_m$ for $m\ll n$ . This is the essence of the Krylov subspace approach.

For instance, we can interpret $\mathbf{A}\mathbf{x}_m\approx \mathbf{b}$ in the sense of linear least-squares—that is, using Theorem 8.4.1 to let $\mathbf{x}=\mathbf{K}_m\mathbf{z}$ ,

\min_{\mathbf{x}\in\mathcal{K}_m} \| \mathbf{A}\mathbf{x}-\mathbf{b} \| = \min_{\mathbf{z}\in\mathbb{C}^m} \| \mathbf{A}(\mathbf{K}_m\mathbf{z})-\mathbf{b} \| = \min_{\mathbf{z}\in\mathbb{C}^m} \| (\mathbf{A}\mathbf{K}_m)\mathbf{z}-\mathbf{b} \|.

(8.4.4)

The natural seed vector for $\mathcal{K}_m$ in this case is the vector $\mathbf{b}$ . In the next example we try to implement (8.4.4). We do take one precaution: because the vectors $\mathbf{A}^{k}\mathbf{b}$ may become very large or small in norm, we normalize after each multiplication by $\mathbf{A}$ , just as we did in the power iteration.

Example 8.4.1 (Conditioning of the Krylov matrix)

Julia

MATLAB

Python

Example 8.4.1

First we define a triangular matrix with known eigenvalues, and a random vector $b$ .

λ = @. 10 + (1:100)
A = triu(rand(100, 100), 1) + diagm(λ)
b = rand(100);

Next we build up the first ten Krylov matrices iteratively, using renormalization after each matrix-vector multiplication.

Km = [b zeros(100, 29)]
for m in 1:29
    v = A * Km[:, m]
    Km[:, m+1] = v / norm(v)
end

Now we solve least-squares problems for Krylov matrices of increasing dimension, recording the residual in each case.

resid = zeros(30)
for m in 1:30
    z = (A * Km[:, 1:m]) \ b
    x = Km[:, 1:m] * z
    resid[m] = norm(b - A * x)
end

The linear system approximations show smooth linear convergence at first, but the convergence stagnates after only a few digits have been found.

using Plots
plot(0:29, resid; m=:o,
    xaxis=(L"m"), yaxis=(:log10, L"\| b-Ax_m \|"),
    title="Residual for linear systems", legend=:none)

Example 8.4.1

First we define a triangular matrix $\mathbf{A}$ with known eigenvalues and a random vector $\mathbf{b}$ .

lambda = 10 + (1:100);
A = diag(lambda) + triu(rand(100), 1); 
b = rand(100, 1);

Next, we build up the first ten Krylov matrices iteratively. In order to keep the columns from growing exponentially in norm, we normalize them as we go. This doesn’t affect the column space of the Krylov matrix, in principle at least.

Km = b;
for m = 1:29      
    v = A * Km(:, m);
    Km(:, m+1) = v / norm(v);
end

Now we solve least-squares problems for Krylov matrices of increasing dimension, recording the residual in each case.

warning off  
resid = zeros(30, 1);
for m = 1:30  
    z = (A * Km(:, 1:m)) \ b;
    x = Km(:, 1:m) * z;
    resid(m) = norm(b - A * x);
end

The linear system approximations show smooth linear convergence at first, but the convergence stagnates after only a few digits have been found.

clf
semilogy(resid, '.-')
xlabel('m'),  ylabel('|| b-Ax_m ||')
set(gca,'ytick',10.^(-6:2:0))
axis tight, title('Residual for linear systems')

Example 8.4.1

First we define a triangular matrix with known eigenvalues, and a random vector $b$ .

ev = 10 + arange(1, 101)
A = triu(random.rand(100, 100), 1) + diag(ev)
b = random.rand(100)

Next we build up the first ten Krylov matrices iteratively, using renormalization after each matrix-vector multiplication.

Km = zeros([100, 30])
Km[:, 0] = b
for m in range(29):
    v = A @ Km[:, m]
    Km[:, m + 1] = v / norm(v)

Now we solve least-squares problems for Krylov matrices of increasing dimension, recording the residual in each case.

from numpy.linalg import lstsq
resid = zeros(30)
resid[0] = norm(b)
for m in range(1, 30):
    z = lstsq(A @ Km[:, :m], b, rcond=None)[0]
    x = Km[:, :m] @ z
    resid[m] = norm(b - A @ x)

The linear system approximations show smooth linear convergence at first, but the convergence stagnates after only a few digits have been found.

semilogy(range(30), resid, "-o")
xlabel("$m$"),  ylabel("$\\| b-Ax_m \\|$")
title(("Residual for linear systems"));

8.4.3The Arnoldi iteration¶

The breakdown of convergence in Demo 8.4.1 is due to a critical numerical defect in our approach: the columns of the Krylov matrix (8.4.1) increasingly become parallel to the dominant eigenvector, as (8.2.4) predicts, and therefore to one another. As we saw in The QR factorization, near-parallel vectors create the potential for numerical cancellation. This manifests as a large condition number for $\mathbf{K}_m$ as $m$ grows, eventually creating excessive error when solving the least-squares system.

The polar opposite of an ill-conditioned basis for $\mathcal{K}_m$ is an orthonormal one. Suppose we had a thin QR factorization of $\mathbf{K}_m$ :

\begin{align*} \mathbf{K}_m = \mathbf{Q}_m \mathbf{R}_m & = \begin{bmatrix} \mathbf{q}_1& \mathbf{q}_2 & \cdots & \mathbf{q}_m \end{bmatrix} \begin{bmatrix} R_{11} & R_{12} & \cdots & R_{1m} \\ 0 & R_{22} & \cdots & R_{2m} \\ \vdots & & \ddots & \\ 0 & 0 & \cdots & R_{mm} \end{bmatrix}. \end{align*}

(8.4.5)

Then the vectors $\mathbf{q}_1,\ldots,\mathbf{q}_m$ are the orthonormal basis we seek for $\mathcal{K}_m$ . By Theorem 8.4.1, we know that $\mathbf{A}\mathbf{q}_m \in \mathcal{K}_{m+1}$ , and therefore

\mathbf{A} \mathbf{q}_m = H_{1m} \, \mathbf{q}_1 + H_{2m} \, \mathbf{q}_2 + \cdots + H_{m+1,m}\, \mathbf{q}_{m+1}

(8.4.6)

for some choice of the $H_{ij}$ . Note that by using orthonormality, we have

\mathbf{q}_i^* (\mathbf{A}\mathbf{q}_m) = H_{im}, \qquad i=1,\ldots, m.

(8.4.7)

Since we started by assuming that we know $\mathbf{q}_1,\ldots,\mathbf{q}_m$ , the only unknowns in (8.4.6) are $H_{m+1,m}$ and $\mathbf{q}_{m+1}$ . But they appear only as a product, and we know that $\mathbf{q}_{m+1}$ is a unit vector, so they are uniquely defined (up to sign) by the other terms in the equation.

We can now proceed iteratively.

Algorithm 8.4.1 (Arnoldi iteration)

Given matrix $\mathbf{A}$ and vector $\mathbf{u}$ :

Let $\mathbf{q}_1= \mathbf{u} \,/\, \| \mathbf{u}\|$ .
For $m=1,2,\ldots$
a. Use (8.4.7) to find $H_{im}$ for $i=1,\ldots,m$ .
b. Let
$\mathbf{v} = (\mathbf{A} \mathbf{q}_m) - H_{1m} \, \mathbf{q}_1 - H_{2m}\, \mathbf{q}_2 - \cdots - H_{mm}\, \mathbf{q}_m.$
(8.4.8)
c. Let $H_{m+1,m}=\|\mathbf{v}\|$ .
d. Let $\mathbf{q}_{m+1}=\mathbf{v}\,/\,H_{m+1,m}$ .

Return the $n\times (m+1)$ matrix $\mathbf{Q}_{m+1}$ whose columns are $\mathbf{q}_1,\dots, \mathbf{q}_{m+1}$ and the $(m+1)\times m$ matrix $\mathbf{H}_m$ .

The vectors $\mathbf{q}_1,\dots, \mathbf{q}_m$ are an orthonormal basis for the space $\mathcal{K}_m$ , which makes them ideal for numerical computations in that space. The matrix $\mathbf{H}_m$ plays an important role, too, and has a particular “triangular plus” structure,

\mathbf{H}_m = \begin{bmatrix} H_{11} & H_{12} & \cdots & H_{1m} \\ H_{21} & H_{22} & \cdots & H_{2m} \\ & H_{32} & \ddots & \vdots \\ & & \ddots & H_{mm} \\ & & & H_{m+1,m} \end{bmatrix}.

(8.4.9)

The identity (8.4.6) used over all the iterations can be collected into a single matrix equation,

\mathbf{A}\mathbf{Q}_m = \mathbf{Q}_{m+1} \mathbf{H}_m,

(8.4.10)

which is a fundamental identity of Krylov subspace methods.

8.4.4Implementation¶

An implementation of the Arnoldi iteration is given in Function 8.4.2. A careful inspection shows that inner nested loop does not exactly implement (8.4.7) and (8.4.8). The reason is numerical stability. Though the described and implemented versions are mathematically equivalent in exact arithmetic (see Exercise 6), the approach in Function 8.4.2 is more stable.

Algorithm 8.4.2 (arnoldi)

Julia

MATLAB

Python

Arnoldi iteration

arnoldi.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
"""
    arnoldi(A, u, m)

Perform the Arnoldi iteration for `A` starting with vector `u`, out
to the Krylov subspace of degree `m`. Returns the orthonormal basis
(`m`+1 columns) and the upper Hessenberg `H` of size `m`+1 by `m`.
"""
function arnoldi(A, u, m)
    n = length(u)
    Q = zeros(n, m+1)
    H = zeros(m+1, m)
    Q[:, 1] = u / norm(u)
    for j in 1:m
        # Find the new direction that extends the Krylov subspace.
        v = A * Q[:, j]
        # Remove the projections onto the previous vectors.
        for i in 1:j
            H[i, j] = dot(Q[:, i], v)
            v -= H[i, j] * Q[:, i]
        end
        # Normalize and store the new basis vector.
        H[j+1, j] = norm(v)
        Q[:, j+1] = v / H[j+1, j]
    end
    return Q, H
end

Arnoldi iteration

arnoldi.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
function [Q, H] = arnoldi(A, u, m)
% ARNOLDI   Arnoldi iteration for Krylov subspaces.
% Input:
%   A    square matrix (n by n)
%   u    initial vector
%   m    number of iterations
% Output: 
%   Q    orthonormal basis of Krylov space (n by m+1)
%   H    upper Hessenberg matrix, A*Q(:,1:m)=Q*H (m+1 by m)

n = length(A);
Q = zeros(n, m+1);  
H = zeros(m+1, m);
Q(:, 1) = u / norm(u);
for j = 1:m
  % Find the new direction that extends the Krylov subspace.
  v = A * Q(:, j);
  % Remove the projections onto the previous vectors.
  for i = 1:j
    H(i, j) = Q(:, i)' * v;
    v = v - H(i,j) * Q(:,i);
  end
  % Normalize and store the new basis vector.
  H(j+1, j) = norm(v);
  Q(:, j+1) = v / H(j+1, j);
end

Arnoldi iteration

arnoldi.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def arnoldi(A, u, m):
    """
    arnoldi(A, u, m)

    Perform the Arnoldi iteration for A starting with vector u, out to the Krylov
    subspace of degree m. Return the orthonormal basis (m+1 columns) and the upper
    Hessenberg H of size m+1 by m.
    """
    n = u.size
    Q = np.zeros([n, m + 1])
    H = np.zeros([m + 1, m])
    Q[:, 0] = u / np.linalg.norm(u)
    for j in range(m):
        # Find the new direction that extends the Krylov subspace.
        v = A @ Q[:, j]
        # Remove the projections onto the previous vectors.
        for i in range(j + 1):
            H[i, j] = Q[:, i] @ v
            v -= H[i, j] * Q[:, i]
        # Normalize and store the new basis vector.
        H[j + 1, j] = np.linalg.norm(v)
        Q[:, j + 1] = v / H[j + 1, j]

    return Q, H

Example 8.4.2 (Arnoldi iteration)

Julia

MATLAB

Python

Example 8.4.2

Here again is the linear system from Example 8.4.1.

λ = @. 10 + (1:100)
A = triu(rand(100, 100), 1) + diagm(λ)
b = rand(100);

We can use $\mathbf{b}$ as the seed vector for the Arnoldi iteration.

Q, H = FNC.arnoldi(A, b, 30)
println("Q has size $(size(Q))")
println("H has size $(size(H))")

Q has size (100, 31)
H has size (31, 30)

Here’s one validation of the key identity (8.4.10).

should_be_near_zero = opnorm(A * Q[:, 1:20] - Q[:, 1:21] * H[1:21, 1:20])

3.0873555092912766e-14

Using the Krylov matrix to project the linear system into a Kyrlov subspace in Example 8.4.1 was unable to get the residual much smaller than about 10^-4. But the Arnoldi basis gives us a stable way to work in that subspace and get better results.

z = (A * Q) \ b
x = Q * z
@show resid_norm = norm(b - A * x);

resid_norm = norm(b - A * x) = 1.5639837383349104e-9

Example 8.4.2

Here again is the linear system from Example 8.4.1.

lambda = 10 + (1:100);
A = diag(lambda) + triu(rand(100), 1); 
b = rand(100, 1);

We can use $\mathbf{b}$ as the seed vector for the Arnoldi iteration.

[Q, H] = arnoldi(A, b, 30);
disp(sprintf("Q is %d by %d", size(Q)))
disp(sprintf("H is %d by %d", size(H)))

Q is 100 by 31

H is 31 by 30

Here’s one validation of the key identity (8.4.10).

should_be_near_zero = norm(A * Q(:, 1:20) - Q(:, 1:21) * H(1:21, 1:20))

z = (A * Q) \ b;
x = Q * z;
resid_norm = norm(b - A * x)

Error using \
Matrix dimensions must agree.

Example 8.4.2

Here again is the linear system from Example 8.4.1.

ev = 10 + arange(1, 101)
A = triu(random.rand(100, 100), 1) + diag(ev)
b = random.rand(100);

We can use $\mathbf{b}$ as the seed vector for the Arnoldi iteration.

Q, H = FNC.arnoldi(A, b, 30)
print("Q has size", Q.shape)
print("H has size", H.shape)

Q has size (100, 31)
H has size (31, 30)

Here’s one validation of the key identity (8.4.10).

from numpy.linalg import norm
should_be_near_zero = norm(A @ Q[:, :20] - Q[:, :21] @ H[:21, :20])
print(should_be_near_zero)

7.366224069215788e-14

z, _, _, _ = linalg.lstsq(A @ Q, b)
x = Q @ z
resid_norm = norm(b - A @ x)
print(f"residual norm: {resid_norm:.2e}")

residual norm: 3.23e-09

In the next section, we revisit the idea of approximately solving $\mathbf{A}\mathbf{x}=\mathbf{b}$ over a Krylov subspace as suggested in Example 8.4.2. A related idea explored in Exercise 7 is used to approximate the eigenvalue problem for $\mathbf{A}$ , which is the approach that underlies eigs for sparse matrices.

8.4.5Exercises¶

✍ Let $\mathbf{A}=\displaystyle \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \end{bmatrix}.$
(a) Find the Krylov matrix $\mathbf{K}_3$ for the seed vector $\mathbf{u}=\mathbf{e}_1$ .
(b) Find $\mathbf{K}_3$ for the seed vector $\mathbf{u}=\begin{bmatrix}1; \: 1;\: 1; \: 1\end{bmatrix}.$
⌨ For each matrix, make a table of the 2-norm condition numbers $\kappa(\mathbf{K}_m)$ for $m=1,\ldots,10$ . Use a vector of all ones as the Krylov seed.
(a) Matrix from Demo 8.4.1
(b) $\begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ & & & 1 & -2 \end{bmatrix} \: (100\times 100)$
(c) $\begin{bmatrix} -2 & 1 & & & 1 \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ 1 & & & 1 & -2 \end{bmatrix} \: (200\times 200)$

✍ Show that if $\mathbf{x}\in\mathcal{K}_m$ , then $\mathbf{x}=p(\mathbf{A})\mathbf{u}$ for a polynomial $p$ of degree at most $m-1$ . (See (7.2.12) for applying a polynomial to a matrix.)
✍ Compute the asymptotic flop requirements for Function 8.4.2. Assume that due to sparsity, a matrix-vector multiplication $\mathbf{A}\mathbf{u}$ requires only $c n$ flops for a constant $c$ , rather than the usual $O(n^2)$ .
⌨ When Arnoldi iteration is performed on the Krylov subspace generated using the matrix $\mathbf{A}=\displaystyle \begin{bmatrix} 2& 1& 1& 0\\ 1 &3 &1& 0\\ 0& 1& 3& 1\\ 0& 1& 1& 2 \end{bmatrix}$ , the results can depend strongly on the initial vector $\mathbf{u}$ .
(a) Apply Function 8.4.2 for 3 iterations and output Q (which should be square) and H when using the following seed vectors.
(i) $\bigl[1,\,0,\,0,\,0\bigr]$ $\qquad$ (ii) $\bigl[1,\,1,\,1,\,1\bigr]$ $\qquad$ (iii) rand(4)
(b) Can you explain why case (ii) in part (a) cannot finish successfully? (Hint: What line(s) of the function can possibly return NaN when applied to finite values?)

✍ As mentioned in the text, Function 8.4.2 does not compute $H_{ij}$ as defined by (8.4.7), but rather
$S_{ij} = \mathbf{q}_i^* ( \mathbf{A}\mathbf{q}_j - S_{1j}\,\mathbf{q}_1 - \cdots - S_{i-1,j}\,\mathbf{q}_{i-1} )$
(8.4.11)
for $i=1,\ldots,j$ . Show that $S_{ij}=H_{ij}$ . (Hence the function is mathematically equivalent to our Arnoldi formulas.)

One way to approximate the eigenvalue problem $\mathbf{A}\mathbf{x}=\lambda\mathbf{x}$ over $\mathcal{K}_m$ is to restrict $\mathbf{x}$ to the low-dimensional spaces $\mathcal{K}_m$ .
(a) ✍ Show starting from (8.4.10) that
$\mathbf{Q}_m^* \mathbf{A} \mathbf{Q}_m = \tilde{\mathbf{H}}_m,$
(8.4.12)
where $\tilde{\mathbf{H}}_m$ is the upper Hessenberg matrix resulting from deleting the last row of $\mathbf{H}_m$ . What is the size of this matrix?
(b) ✍ Show the reasoning above leads to the approximate eigenvalue problem $\tilde{\mathbf{H}}_m\mathbf{z} \approx \lambda\mathbf{z}$ . (Hint: Start with $\mathbf{A}\mathbf{x} \approx \lambda\mathbf{x}$ , and let $\mathbf{x}=\mathbf{Q}_m\mathbf{z}$ before applying part (a).)
(c) ⌨ Apply Function 8.4.2 to the matrix of Demo 8.4.1 using a random seed vector. Compute eigenvalues of $\tilde{\mathbf{H}}_m$ for $m=1,\ldots,40$ , keeping track in each case of the error between the largest of those values (in magnitude) and the largest eigenvalue of $\mathbf{A}$ . Make a log-linear graph of the error as a function of $m$ .

Footnotes¶

The proper pronunciation of “Krylov” is something like “kree-luv,” but American English speakers often say “kreye-lahv.”
↩