The Arnoldi Process#
The previous section introduced the projected eigenvalue problem on a Krylov subspace. Arnoldi constructs its orthonormal basis and the small Hessenberg matrix together. Each step uses one matrix-vector product and removes the components along the directions already found.
We first derive this construction in exact arithmetic for \(A\in\mathbb R^{n\times n}\). For complex matrices, use conjugate transposes and complex inner products. The implementation below supports both cases.
Building One New Direction#
Start with a unit vector \(q_1\). Suppose we have constructed orthonormal vectors \(q_1,\ldots,q_j\), and write \(Q_j=[q_1,\ldots,q_j]\). To extend the basis, multiply the newest vector by \(A\):
Part of \(v\) may already lie in the current subspace. Its coordinates along the existing basis vectors are
Remove these components:
The remainder is perpendicular to every existing basis vector. For \(\ell=1,\ldots,j\),
If \(w\ne0\), normalize it to obtain the next direction:
Rearranging the construction gives
If \(w=0\) instead, no new direction can be added and the process stops. This is called breakdown; its meaning for the eigenvalue problem is explained below.
This is Gram-Schmidt applied to the newly generated vector \(Aq_j\). Unlike orthogonal iteration, which updates a whole basis at each step, Arnoldi keeps the existing vectors and adds one new direction.
The Subspace Generated by Arnoldi#
The construction keeps the basis orthonormal. We also need to verify that it spans the Krylov subspace chosen on the previous page.
Theorem 30 (Arnoldi Produces a Krylov Basis)
If Arnoldi constructs \(q_1,\ldots,q_k\) without earlier breakdown, then these vectors form an orthonormal basis of \(\mathcal K_k(A,q_1)\).
Proof. The claim is immediate for \(k=1\). Suppose it holds for the first \(j\) vectors. The Arnoldi construction gives
The inclusion holds because multiplication by \(A\) raises the largest power by at most one. When \(h_{j+1,j}\ne0\), this is a new unit vector perpendicular to the first \(j\). Thus we have \(j+1\) orthonormal vectors in a space generated by \(j+1\) vectors. They must span that space, completing the induction.
There is also a QR interpretation. The Krylov matrix
has a factorization \(K_k=Q_kR_k\). The matrix \(R_k\) is upper triangular because column \(j\) of \(K_k\) lies in the span of the first \(j\) columns of \(Q_k\).
Arnoldi obtains this basis without forming \(K_k\) explicitly. Successive powers can become nearly parallel, making the raw Krylov matrix poorly conditioned. Generating each new direction from the existing orthonormal basis avoids explicitly storing and factoring these nearly parallel columns.
The Arnoldi Relation#
After \(k\) steps, collect the first \(k\) basis vectors in \(Q_k\) and the coefficients in a \(k\times k\) matrix \(H_k\). Column \(j\) has no entries below row \(j+1\), so \(H_k\) is upper Hessenberg.
For \(j<k\), the new direction \(q_{j+1}\) is already a column of \(Q_k\). Only the last product, \(Aq_k\), can have a component outside this basis. Let \(f_k\) be the remainder from this final step:
When \(f_k\ne0\), normalization gives \(f_k=h_{k+1,k}q_{k+1}\). If it is zero, there is no new direction to add.
Stacking the column equations gives the Arnoldi relation:
where \(e_k\) is the last coordinate vector in \(\mathbb R^k\). The matrix \(f_ke_k^T\) has zeros in its first \(k-1\) columns and \(f_k\) in its last column.
Since \(Q_k^TQ_k=I_k\) and \(Q_k^Tf_k=0\), multiplication on the left by \(Q_k^T\) gives
Thus the coefficients collected during orthogonalization are exactly the projected matrix from the previous section. We do not need to form \(AQ_k\) again to compute it.
When the Subspace Stops Growing#
At breakdown in step \(j\), the remainder is zero and \(h_{j+1,j}=0\). The Arnoldi relation becomes
The current space is invariant under \(A\). For any eigenpair \(H_jy=\theta y\),
Every Ritz pair from this space is therefore an exact eigenpair of \(A\). This is often called happy or lucky breakdown. It does not mean that we have found every eigenvalue of \(A\): the starting vector may generate only a proper invariant subspace.
In floating-point arithmetic, the remainder may be very small without being exactly zero. Dividing by its norm can magnify rounding errors, so implementations use a small threshold to decide when to stop expanding. Such a decision indicates approximate invariance; it does not make the computed Ritz pairs exact.
Implementation#
The implementation uses modified Gram-Schmidt: each component is subtracted immediately, and the next inner product uses the current remainder. In exact arithmetic, this gives the coefficients and remainder derived above.
It normalizes the starting vector and performs a second orthogonalization pass by default to help preserve orthogonality despite rounding. It returns the basis, the projected matrix, and the final remainder \(f_k\). The next section uses this remainder to check Ritz residuals.
Here \(k\) is the maximum number of Arnoldi steps. If the process stops earlier, the returned matrices have the smaller dimensions.
import numpy as np
def arnoldi_iteration(A, q1, k, reorthogonalize=True,
breakdown_tol=None):
"""Return Q, H, f with A @ Q approximately Q @ H + outer(f, e_last).
A is a square matrix or linear operator supporting A @ v.
Q has at most k orthonormal columns, and H is square Hessenberg.
f is a vector: the remainder from the last completed step.
breakdown_tol is relative to the norm of A @ q_j.
"""
q = np.asarray(q1)
if q.ndim != 1 or A.shape != (q.size, q.size):
raise ValueError("A must be square and q1 a matching vector.")
if not isinstance(k, (int, np.integer)) or not 1 <= k <= q.size:
raise ValueError("k must be an integer between 1 and n.")
dtype = np.result_type(A.dtype, q.dtype, np.float64)
q = q.astype(dtype, copy=True)
qnorm = np.linalg.norm(q)
if qnorm == 0:
raise ValueError("The starting vector must be nonzero.")
q /= qnorm
if breakdown_tol is None:
breakdown_tol = 10 * np.finfo(q.real.dtype).eps
if breakdown_tol < 0:
raise ValueError("breakdown_tol must be nonnegative.")
Q = [q]
H = np.zeros((k, k), dtype=dtype)
passes = 2 if reorthogonalize else 1
for j in range(k):
v = A @ Q[j]
scale = np.linalg.norm(v)
# Accumulate coefficients from all orthogonalization passes.
for _ in range(passes):
for i in range(j + 1):
h = np.vdot(Q[i], v)
H[i, j] += h
v = v - h * Q[i]
beta = np.linalg.norm(v)
if beta <= breakdown_tol * scale or j + 1 == k:
m = j + 1
return np.column_stack(Q), H[:m, :m].copy(), v
H[j + 1, j] = beta
Q.append(v / beta)
The last step still computes its remainder, even though no further basis vector is needed. If the requested number of steps is reached, \(f_k\) need not be small: stopping after \(k\) steps does not by itself imply an accurate eigenpair.
The threshold in this example is a simple relative test, not an eigenvalue accuracy tolerance. In finite precision, the basis is only approximately orthonormal and the Arnoldi relation holds up to rounding errors. Reorthogonalization helps preserve the orthogonality on which the projection relies.
Computational Cost and the Next Step#
For \(k\) steps without early termination, Arnoldi uses \(k\) matrix-vector products, \(O(nk^2)\) work for orthogonalization, and \(O(nk)\) storage for the basis. A second orthogonalization pass changes the constant, not these orders. Computing the eigenvalues of the small dense \(H_k\) adds work of order \(O(k^3)\).
Keeping \(k\) small makes the method useful for large matrices. The next section explains how to extract Ritz pairs, check their residuals, and understand their accuracy through polynomials of \(A\).