24.6  Krylov Subspaces

Sections 4 and 5 both consist of multiplying by A over and over: the power method forms b,Ab,A2b,… and keeps only the newest vector, and the QR algorithm, by Theorem 24.5.2, does the same thing to e1 inside its first column. In both cases everything computed before the last step is thrown away, and this section keeps it instead. The span of the vectors produced by k matrix-vector products is a subspace of dimension up to k, and the best approximations it contains are far better than its last vector alone. That subspace, the process that builds an orthonormal basis of it, and the eigenvalue estimates it yields are the subject here.

Throughout, F=R or F=C, and Fn carries the standard inner product, linear in the first slot, with ⟨x,y⟩=y∗x and ‖⋅‖=‖⋅‖2.

24.6.1. The space that k products span

A large matrix is often not available as a table of n2 numbers at all. What is available is a rule for computing Ax from x — a sparse pattern of neighbors on a grid, a graph, a differential operator discretized. For such a matrix, forming Ax is cheap and forming A−1, or even looking at an entry of A2, is not. So the honest question is: after k multiplications, what do we possess?

We possess b,Ab,…,Ak−1b, and hence every linear combination of them. The power method uses one of these vectors. That is wasteful, and the waste has a name only because the span deserves one.

The Krylov space is everything that k matrix-vector products can reach.

Definition 24.6.1 (Krylov Subspace)

Let A∈Mn(F), let b∈Fn and let k≥1 be an integer. The k-th Krylov subspace of the pair (A,b) is Kk(A,b):=span(b, Ab, A2b, …, Ak−1b)⊆Fn.

In words: we list the first k vectors of the sequence obtained by repeatedly applying A to b — starting at b itself, so that the last one is Ak−1b and not Akb — and take their span. It is a subspace for the usual reason, being a span. The index k counts vectors, not powers; the off-by-one is a convention, and it is the one that makes dim⁡Kk≤k with equality in the typical case.

Examples. Take F=R throughout.

  • The identity. For A=In and any b≠0, every Ajb equals b, so Kk(In,b)=span(b) for every k. One product buys nothing, and neither do a hundred.
  • A companion matrix. Let Aej=ej+1 for j≤n−1, as for the companion matrix of Definition 10.5.3, and let b=e1. Then Aje1=ej+1 for j≤n−1, so Kk(A,e1)=span(e1,…,ek) for k≤n, of dimension exactly k. Here every product buys a new dimension.
  • A diagonal matrix. Let A=diag(4,3,2,1) and b=(1,1,1,1). Then Ajb=(4j,3j,2j,1), and these four vectors for j=0,1,2,3 form a Vandermonde matrix with distinct nodes, hence are independent (Theorem 7.8.1). So dim⁡Kk=k for k≤4, and K4=R4.
  • Degenerate cases. K1(A,b)=span(b); and Kk(A,0)={0} for every k, which is why b≠0 is assumed from here on.

Non-example by minimal change. Drop the middle vector of the list. Let A∈M3(R) satisfy Ae1=e2, Ae2=e3, Ae3=0, and let b=e1. The set span(b,A2b)=span(e1,e3) is a perfectly good subspace, and it is not Kk(A,b) for any k, since those are span(e1), span(e1,e2) and R3. The clause that fails is that the powers be consecutive, and what it costs is everything: Aspan(e1,e3)=span(e2), which meets span(e1,e3) only at 0. The Krylov spaces, by contrast, are nested and almost invariant, which is the next proposition.

Why this definition. Chapter 10 already met the union of all of them. The cyclic subspace Z(b;A)=span(b,Ab,A2b,…) of Definition 10.5.1 is the smallest A-invariant subspace containing b (Proposition 10.5.1), and Theorem 10.5.2 says that the list stops growing exactly at the degree of the A-annihilator mA,b. The Krylov space is that construction truncated at k, and the truncation is the whole point: Z(b;A) is a theoretical object, while Kk(A,b) is what k units of work actually deliver.

Proposition 24.6.1 (Krylov Spaces: Nesting, Dimension, Invariance)

Let A∈Mn(F) and b∈Fn with b≠0, and let d=deg⁡mA,b≥1 be the degree of the A-annihilator of b (Definition 10.5.2). Then, for every k≥1:

  • (a) Kk(A,b)={p(A)b:p∈F[x], deg⁡p≤k−1}, and Kk⊆Kk+1 and AKk⊆Kk+1;
  • (b) dim⁡Kk(A,b)=min(k,d);
  • (c) Kk+1=Kk if and only if k≥d, and in that case Kk=Z(b;A) is A-invariant.

Idea. Everything is Theorem 10.5.2 read one index at a time. That theorem says (b,Ab,…,Ad−1b) is a basis of Z(b;A); so the first d vectors of our list are independent, and from Adb on they add nothing.

Proof.

  • (a) A linear combination ∑j<kajAjb is p(A)b for p=∑j<kajxj, and conversely, which is the first statement. The inclusion Kk⊆Kk+1 holds because the spanning list of Kk is part of that of Kk+1. And AAjb=Aj+1b∈Kk+1 for j≤k−1, so A maps a spanning list of Kk into Kk+1, hence maps Kk into Kk+1.

  • (b) By Theorem 10.5.2 (a), the list (b,Ab,…,Ad−1b) is a basis of Z(b;A). If k≤d, the first k entries of that list are a sublist of a basis, hence independent, so dim⁡Kk=k. If k>d, then Kk⊆Z(b;A), because every Ajb lies there (Proposition 10.5.1 (a)); and Kk⊇Kd=Z(b;A), since Kd already contains a basis of it. So Kk=Z(b;A) and dim⁡Kk=d. In both cases dim⁡Kk=min(k,d).

  • (c) (⇐) If k≥d, then (b) gives dim⁡Kk+1=d=dim⁡Kk, and with Kk⊆Kk+1 from (a) this forces equality; and the common space is Z(b;A), A-invariant by Proposition 10.5.1 (b). (⇒) If k<d, then (b) gives dim⁡Kk+1=k+1>k=dim⁡Kk, so the two differ. This proves the proposition.

Warning. The list b,Ab,…,Ak−1b is a basis of Kk and is useless as one. By Section 4’s reading of the power method (Theorem 24.4.1), Ajb turns towards the dominant eigenvector as j grows, so the later vectors of the list are nearly parallel to each other. Even for A=diag(4,3,2,1) and b=(1,1,1,1), with its modest ratio λ1/λ2=4/3, the angle between consecutive vectors of the list falls from 24 degrees to 8.9 degrees across those four vectors, and the matrix with those four columns has 2-norm condition number about 1200 — for a 4×4 problem with one-digit data. Independent they may be; in the arithmetic of Definition 24.1.1 the coordinates of a vector in that basis are computed with errors magnified by that condition number, which grows quickly with k. A Krylov method never uses this basis. It builds an orthonormal one instead, and that is the next subsection.

Quick check. Let A∈Mn(F) be invertible with A2=In, and let b≠0 be a vector that is not an eigenvector of A. What is dim⁡Kk(A,b) for each k≥1?

Solution. It is 1 for k=1 and 2 for every k≥2. The polynomial x2−1 kills b, so mA,b divides it and has degree at most 2 (Definition 10.5.2). Degree 1 would mean mA,b=x−λ, that is Ab=λb, making b an eigenvector, which is excluded; and deg⁡mA,b≥1 because b≠0. So d=2, and Proposition 24.6.1 (b) gives min(k,2).

24.6.2. Arnoldi: Gram–Schmidt with a purpose

We want an orthonormal basis q1,…,qk of Kk, with the nesting preserved: span(q1,…,qj)=Kj for every j≤k. Running Theorem 11.2.2 on the list b,Ab,…,Ak−1b would do it, and would inherit exactly the defect the warning describes, since it needs those vectors first.

The repair is one line long, and it is the idea of the whole subsection. Do not orthogonalize Ajb. Orthogonalize Aqj. The vector Aqj lies in Kj+1 by Proposition 24.6.1 (a), and subtracting its components along q1,…,qj leaves a vector in Kj+1 orthogonal to Kj — which is what the next basis vector must be. The powers of A are never formed.

Algorithm. The Arnoldi process. Given A∈Mn(F) and b≠0, set q1=b/‖b‖. For j=1,2,…: hij=⟨Aqj,qi⟩(1≤i≤j),wj=Aqj−∑i=1jhijqi,hj+1,j=‖wj‖,qj+1=wj/hj+1,jif hj+1,j≠0. If hj+1,j=0 the process stops; this is called a breakdown, and the theorem below says it is good news.

Write Qj=(q1∣⋯∣qj)∈Mn×j(F), let Hk∈Mk(F) have entries hij for i≤j+1 and zeros elsewhere, and let H~k∈M(k+1)×k(F) be Hk with the extra row hk+1,kek⊤ appended at the bottom.

Theorem 24.6.2 (The Arnoldi Process)

Let A∈Mn(F), let b∈Fn be non-zero, let d=deg⁡mA,b, and run the Arnoldi process. Then the Arnoldi relation

AQk=QkHk+hk+1,kqk+1ek⊤(1≤k≤d)
(24.6.1)

holds, the last term being 0 when k=d. Moreover:

  • (a) For every j≤d, the list (q1,…,qj) is an orthonormal basis of Kj(A,b); moreover hj+1,j>0 for j<d, and hd+1,d=0, so the process breaks down at step d and at no earlier step.
  • (b) For 1≤k≤d we have Qk∗Qk=Ik, and Hk=Qk∗AQk is upper Hessenberg. For k<d, Eq. 24.6.1 can also be written AQk=Qk+1H~k.
  • (c) A breakdown means an invariant subspace: if hk+1,k=0, then AKk⊆Kk, and every eigenvalue of Hk is an eigenvalue of A.

Idea. Only one thing has to be checked, and it is the nesting: if q1,…,qj is an orthonormal basis of Kj, then Aqj lies in Kj+1 but generally not in Kj, and what the step subtracts is exactly its Kj-component. The matrix identity is then just the definition of wj rearranged, column by column, and the Hessenberg shape is the statement that Aqj needs only q1,…,qj+1 — never qj+2 — because AKj⊆Kj+1.

Proof.

  • (a) Induction on j. For j=1, q1=b/‖b‖ is a unit vector spanning K1=span(b). Let 1≤j and suppose (q1,…,qj) is an orthonormal basis of Kj. The subtraction in the process makes wj orthogonal to each qi with i≤j: for such i, ⟨wj,qi⟩=⟨Aqj,qi⟩−∑l=1jhlj⟨ql,qi⟩=hij−hij=0, using orthonormality. Also Aqj∈AKj⊆Kj+1 by Proposition 24.6.1 (a), and each qi∈Kj⊆Kj+1, so wj∈Kj+1.

Claim. Kj+1=Kj+span(Aqj).

Proof. (⊇) Both summands lie in Kj+1, by Proposition 24.6.1 (a). (⊆) Since Kj+1=Kj+span(Ajb), it is enough to place Ajb on the right. Write Aj−1b=∑i≤jciqi, possible because Aj−1b∈Kj and the qi are a basis of it. Applying A, Ajb=∑i≤jciAqi, and for i≤j−1 we have Aqi∈AKi⊆Ki+1⊆Kj. So only the term cjAqj leaves Kj.

Now wj=0 if and only if Aqj∈Kj, and by the Claim that happens if and only if Kj+1=Kj. Suppose first that j<d. Then dim⁡Kj+1=j+1>j=dim⁡Kj by Proposition 24.6.1 (b), so Kj+1≠Kj and wj≠0. Hence hj+1,j>0, and qj+1=wj/hj+1,j is a unit vector in Kj+1 orthogonal to q1,…,qj. An orthogonal list of non-zero vectors is independent, so (q1,…,qj+1) is an orthonormal list of j+1 vectors in the (j+1)-dimensional space Kj+1, hence a basis of it.

If instead j=d, then Kd+1=Kd by Proposition 24.6.1 (c), so Aqd∈AKd⊆Kd+1=Kd and wd=0. This gives (a).

(b), and Eq. 24.6.1. The columns of Qk are orthonormal by (a), which is Qk∗Qk=Ik. Rearranging the definition of wj, Aqj=∑i=1jhijqi+hj+1,jqj+1=∑i=1j+1hijqi, for each j≤k, where for j=d the last term is 0. For j<k the right-hand side is QkHkej, since column j of Hk has entries h1j,…,hj+1,j and zeros below and j+1≤k; for j=k it is QkHkek+hk+1,kqk+1. Collecting the columns gives Eq. 24.6.1, and the form AQk=Qk+1H~k is the same identity with the extra column of Qk+1 and the extra row of H~k absorbing the last term. The entries hij with i>j+1 are zero by definition, which is exactly upper Hessenberg form. Finally, multiplying Eq. 24.6.1 on the left by Qk∗ and using Qk∗Qk=Ik and Qk∗qk+1=0, Qk∗AQk=Hk+hk+1,kQk∗qk+1ek⊤=Hk.

  • (c) If hk+1,k=0, then Eq. 24.6.1 reads AQk=QkHk, so Aqj∈Kk for every j≤k and hence AKk⊆Kk. If Hky=θy with y≠0, put x=Qky, which is non-zero because the columns of Qk are independent; then Ax=AQky=QkHky=θQky=θx. This proves the theorem.

The cost of k steps is k matrix-vector products plus the orthogonalizations, which use about 4nj flops at step j — j inner products and j subtractions, about 2n flops each — so O(nk2) in total; and all k+1 vectors must be stored. Both grow with k, which is why Arnoldi is restarted in practice. For Hermitian A neither grows, and that is the next theorem.

24.6.3. Lanczos: the recurrence collapses to three terms

Theorem 24.6.3 (The Lanczos Recurrence)

Let A∈Mn(C) be Hermitian, let b≠0, and run the Arnoldi process. Then Hk is Hermitian for every k≤d, hence tridiagonal, with real diagonal entries; writing αj=hjj∈R and βj=hj+1,j, which is >0 for j<d and 0 for j=d, we have hj,j+1=βj and

Aqj=βj−1qj−1+αjqj+βjqj+1(1≤j≤d),
(24.6.2)

with β0q0 read as 0 and βdqd+1 as 0. Consequently Hk=Tk is the real symmetric tridiagonal matrix with diagonal (α1,…,αk) and off-diagonal (β1,…,βk−1), and each step of the process needs only qj−1 and qj.

Idea. Two facts collide. The matrix Hk=Qk∗AQk is Hermitian because A is, and it is upper Hessenberg by Theorem 24.6.2. A Hermitian matrix with zeros below the subdiagonal has zeros above the superdiagonal too, by conjugating. Everything else is bookkeeping. The reason is worth naming separately: hij=⟨Aqj,qi⟩ lets A move to the other slot for free, and Aqi lives one index further along than qi does.

Proof. By Theorem 24.6.2 (b), Hk=Qk∗AQk, so Hk∗=Qk∗A∗Qk=Qk∗AQk=Hk and Hk is Hermitian; in particular its diagonal entries are real.

Let i≤j−2. Since A∗=A is self-adjoint, hij=⟨Aqj,qi⟩=⟨qj,Aqi⟩. Now qi∈Ki, so Aqi∈Ki+1⊆Kj−1 by Proposition 24.6.1 (a) and i+1≤j−1, while qj is orthogonal to Kj−1 by Theorem 24.6.2 (a). Hence hij=0 whenever i≤j−2. Together with hij=0 for i>j+1, this leaves only i∈{j−1,j,j+1}, so Hk is tridiagonal. Being Hermitian, hj,j+1=hj+1,j―=βj, which is real and positive. The displayed recurrence Eq. 24.6.2 is then the expansion Aqj=∑ihijqi of Theorem 24.6.2 (b) with only three surviving terms; and Hk has real entries, is symmetric and is tridiagonal, as claimed. Since the step needs only hj−1,j=βj−1, hjj=αj and qj−1,qj, this proves the theorem.

This is Chapter 11’s three-term recurrence, and for the same reason. Chapter 11 §10 built the monic orthogonal polynomials for an inner product on R[x] and found (Theorem 11.10.3) that pk+1=(x−ak)pk−bkpk−1, with only three terms; the reason given there was Proposition 11.10.1, “multiplication by x moves from one slot to the other for free”. The two statements are the same statement. Given Hermitian A and b≠0, define on R[x] ⟨p,q⟩b:=⟨p(A)b,q(A)b⟩. This is symmetric and bilinear, and on R[x]≤d−1 it is positive definite: ⟨p,p⟩b=‖p(A)b‖2 vanishes only when p(A)b=0, which by Definition 10.5.2 means mA,b∣p, impossible for a non-zero p of degree below d. Multiplication by x is symmetric for it, since A is Hermitian: ⟨xp,q⟩b=⟨Ap(A)b,q(A)b⟩=⟨p(A)b,Aq(A)b⟩=⟨p,xq⟩b. Now Proposition 24.6.1 (a) says qj+1 is r(A)b for some r of degree at most j, and orthogonality to Kj says ⟨r,s⟩b=0 for every s of degree below j. So for j≤d−1 the polynomial r is a scalar multiple of the j-th monic orthogonal polynomial pj for ⟨⋅,⋅⟩b, and Eq. 24.6.2, at those indices, is the recurrence of Theorem 11.10.3 rewritten in terms of the vectors pj(A)b — read, as that theorem’s node inner product is read, only in the range of degrees where the form is positive definite. The last step j=d lies outside that range: there the monic polynomial of degree d orthogonal to everything of lower degree is mA,b itself, and ⟨mA,b,mA,b⟩b=0. A third appearance of the same three terms is Chapter 7’s: by Theorem 7.8.3, the characteristic polynomials Dk(x)=det(xIk−Tk) of the leading blocks satisfy Dk=(x−αk)Dk−1−βk−12Dk−2, which is Theorem 11.10.3 with ak−1=αk and bk−1=βk−12>0. The Lanczos coefficients, the orthogonal polynomials and the determinants of the tridiagonal blocks are one object seen three ways.

Example (A breakdown at step two)

Let A=(311131113) and b=e1. Run the Lanczos recurrence of Theorem 24.6.3, and identify what it finds.

Solution. q1=e1, and Aq1=(3,1,1). So α1=⟨Aq1,q1⟩=3 and w1=(3,1,1)−3(1,0,0)=(0,1,1),β1=2,q2=12(0,1,1). Next, Aq2=12(2,4,4), so α2=⟨Aq2,q2⟩=12(0+4+4)=4 and w2=12(2,4,4)−4⋅12(0,1,1)−2(1,0,0)=12(2,0,0)−2(1,0,0)=0. The process breaks down at step 2, so d=2 and K2 is A-invariant, by Theorem 24.6.2 (a), (c). Then T2=(3224),pT2(x)=(x−3)(x−4)−2=x2−7x+10, with roots 5 and 2. And indeed spec(A)={5,2}, with 5 simple and 2 of multiplicity 2: the two steps have produced both eigenvalues exactly. This is Theorem 24.6.2 (c) in action, and it is also Proposition 24.6.1 (b) with d=deg⁡mA,e1=2, since mA=(x−5)(x−2) and neither factor alone kills e1.

24.6.4. Ritz values

An invariant subspace would give eigenvalues exactly. A Krylov space is an almost invariant subspace — by Eq. 24.6.1 it fails to be invariant only through the single vector hk+1,kqk+1 — so it should give almost-eigenvalues, and it should be possible to say how almost.

Definition 24.6.2 (Ritz Values and Ritz Vectors)

Let A∈Mn(F), let Qk∈Mn×k(F) have orthonormal columns, and put Hk=Qk∗AQk. The eigenvalues θ1,…,θk of Hk are the Ritz values of A from the subspace span(q1,…,qk), and for an eigenvector y of Hk with Hky=θy and ‖y‖=1, the vector x=Qky is the corresponding Ritz vector; the pair (θ,x) is a Ritz pair.

In words: project the problem onto the subspace, solve the small problem exactly, and lift the answer back. The Ritz values depend on the subspace and not on which orthonormal basis of it we chose: replacing Qk by QkV with V unitary replaces Hk by V∗HkV, which has the same eigenvalues. Note that ‖x‖=‖Qky‖=‖y‖=1, since Qk has orthonormal columns.

Theorem 24.6.4 (Ritz Values Interlace, and Their Residual Is One Number)

Let A∈Mn(C) be Hermitian with eigenvalues λ1(A)≥⋯≥λn(A), let 1≤k≤d, and let θ1≥⋯≥θk be the Ritz values from Kk(A,b), that is, the eigenvalues of Tk.

  • (a) λi(A) ≥ θi ≥ λi+n−k(A) for 1≤i≤k.
  • (b) If Tky=θy with ‖y‖=1 and x=Qky, then Ax−θx=βkykqk+1,so‖Ax−θx‖=βk|yk|, where yk is the last coordinate of y and βk=hk+1,k, with βdqd+1 read as 0.

Idea. For (a), the obstacle is that Tk=Qk∗AQk is a compression, not a principal submatrix, and Chapter 17’s interlacing is about principal submatrices. The repair is to complete Qk to a unitary matrix: then Tk is the leading principal submatrix of a matrix unitarily similar to A, and Corollary 17.4.3 applies unchanged. For (b), just apply Eq. 24.6.1 to y and watch the two QkTky terms cancel.

Proof.

  • (a) Extend q1,…,qk to an orthonormal basis q1,…,qn of Cn, which is possible by Theorem 11.2.2 applied to any extension of the list to a basis (Theorem 1.6.5), and let U=(q1∣⋯∣qn), which is unitary. Then B=U∗AU is Hermitian and unitarily similar to A, so λi(B)=λi(A) for every i. Its entries are bij=⟨Aqj,qi⟩, and for i,j≤k these are exactly the entries of Qk∗AQk=Tk. So Tk=BI,I with I={1,…,k}, a principal submatrix of B obtained by deleting m=n−k rows and the same columns. Now Corollary 17.4.3 gives λi(B) ≥ λi(Tk) ≥ λi+m(B)(1≤i≤k), which is the assertion, since λi(Tk)=θi and λi(B)=λi(A).

  • (b) By Eq. 24.6.1, AQk=QkTk+βkqk+1ek⊤. Applying both sides to y and using Tky=θy and ek⊤y=yk, Ax=AQky=QkTky+βkykqk+1=θQky+βkykqk+1=θx+βkykqk+1. Subtracting θx gives the displayed identity, and taking norms gives the second form, since ‖qk+1‖=1 and βk≥0. (For k=d the right-hand side is 0, and the Ritz pair is an exact eigenpair, as Theorem 24.6.2 (c) already said.) This proves the theorem.

Part (b) is where the chapter’s perturbation theory pays off, and it pays off for free. The residual ‖Ax−θx‖ is the one quantity Chapter 20 §06 asks for, and here it is available without forming Ax: it is the product of the number βk the process has just computed and one coordinate of the small eigenvector y. So Theorem 20.6.1 applies verbatim, and says that A has an eigenvalue within βk|yk| of θ. Better still, θ=⟨Ax,x⟩ is the Rayleigh quotient of x — because ⟨Ax,x⟩=y∗Qk∗AQky=y∗Tky=θ — so Theorem 20.6.2 applies too, and turns the same number into a bound of size (βk|yk|)2 divided by the distance to the rest of the spectrum. Nothing in this section had to be proved twice.

Warning. A Ritz value that happens to equal an eigenvalue is not a converged one, and for a non-Hermitian A there is no interlacing to fall back on. Take A=(211021001),b=e3. Two Arnoldi steps give q1=e3, q2=12(1,1,0) and H2=(10252), hence the Ritz values 52 and 1; Exercise (B2: Two Arnoldi steps) works the details. The second is an exact eigenvalue of A — and its Ritz vector 117(−2,−2,3) is not an eigenvector of A, the residual being 34/17≈0.343, nowhere near zero. What certifies convergence is the residual, never the agreement of a number with a number. Part (a) of Theorem 24.6.4 needs A Hermitian; without that hypothesis there is no interlacing and, by Theorem 20.4.4, a residual bound carries the factor κ2(X).

Example (Ritz values closing in on the extremes)

Let A=diag(4,3,2,1) and b=(1,1,1,1). Compute the Lanczos coefficients, the Ritz values for k=1,2,3, and the residual of the largest Ritz pair at k=3.

Solution. Here q1=12(1,1,1,1) and Aq1=12(4,3,2,1), so α1=14(4+3+2+1)=52. Then w1=12(4,3,2,1)−52⋅12(1,1,1,1)=14(3,1,−1,−3),β1=204=52, and q2=120(3,1,−1,−3). Continuing in the same way gives α2=α3=α4=52 as well, together with β2=25 and β3=325.

The repetition of 52 is a symmetry, not an accident. Let S be the reversal permutation matrix, Sej=e5−j, and put B=A−52I=diag(32,12,−12,−32). Then S is orthogonal with S2=I, and SBS=−B, Sb=b. We show by induction, simultaneously, that Sqj=εjqj with εj=(−1)j−1 and that αj=52. The first holds at j=1. Given it at j, αj−52=⟨Bqj,qj⟩=⟨SBqj,Sqj⟩=⟨SBS(Sqj),Sqj⟩=−εj2⟨Bqj,qj⟩, and εj2=1, so this number equals its own negative and αj=52. The Lanczos step then reads wj=Bqj−βj−1qj−1, and applying S with εj−1=−εj gives Swj=−εjwj, so Sqj+1=εj+1qj+1. The Ritz values are the eigenvalues of the leading blocks of T4:

k Ritz values check against 4,3,2,1
1 52 4≥2.5≥1
2 5±52≈3.6180, 1.3820 4≥3.618≥2; 3≥1.382≥1
3 52±20510≈3.9318, 1.0682, and 52 4≥3.932≥3; 3≥2.5≥2; 2≥1.068≥1
4 4,3,2,1 exact

Every row obeys Theorem 24.6.4 (a), with m=n−k shifting the index on the right. The extreme Ritz values close in on 4 and 1 fast, and the interior ones lag: at k=3 the outer two are within 0.07 of 4 and 1, while the single interior value 52 is half a unit from each of 3 and 2.

The residual at k=3. The eigenvector of T3 for θ1=52+20510 is, normalized, y=(582, 12, 482), so Theorem 24.6.4 (b) gives ‖Ax−θ1x‖=β3|y3|=325⋅482=3410205≈0.2963. By Theorem 20.6.1, A has an eigenvalue within 0.2963 of 3.9318; the truth is 4, at distance 0.0682, so the certificate costs a factor of about 4. Suppose we also know that A has no eigenvalue in the open interval (3,4). Then Theorem 20.6.2 (b) with α=3 improves the upper end to θ1+(0.2963)23.9318−3=3.9318+0.0942=4.0260, and the lower end is θ1 itself by interlacing. The bracket [3.9318, 4.0260] has width 0.094 where the residual bound gave 0.593. That is the quadratic gain of Chapter 20 §06, obtained here from one number the process had already produced.

24.6.5. How fast the extreme Ritz values converge

The example suggests the extremes converge quickly. Here is why, and how much of it we can prove.

Theorem 24.6.5 (A Rate for the Largest Ritz Value)

Let A∈Mn(C) be Hermitian with eigenvalues λ1>λ2≥⋯≥λn and an orthonormal eigenbasis v1,…,vn. Let b be a unit vector with b=∑icivi and c1≠0, and put τ2=∑i≥2|ci|2|c1|2,ρ=λ2−λn2λ1−λ2−λn∈[0,1). Then the largest Ritz value θ1(k) from Kk(A,b) satisfies, for 1≤k≤d, 0 ≤ λ1−θ1(k) ≤ (λ1−λn)τ2ρ2(k−1), and θ1(k) is non-decreasing in k.

Idea. Three moves. ① The largest Ritz value is the largest Rayleigh quotient over the Krylov space, by the Rayleigh principle transported through Qk. ② Every vector of the Krylov space is p(A)b for a polynomial of degree below k, so we may choose p and get a bound for free: whichever p makes p large at λ1 and small at the other eigenvalues gives a good one. ③ Take the simplest such p, the (k−1)-st power of x minus the midpoint of [λn,λ2]. That is not the best choice, and the best one is the subject of the next section.

Proof. Step 1. θ1(k)=max{RA(x):x∈Kk(A,b), x≠0}. Indeed θ1(k)=λ1(Tk)=maxy≠0RTk(y) by Proposition 17.1.1, and the map y↦Qky is a bijection from Ck∖{0} onto Kk∖{0}, because the columns of Qk are a basis of Kk (Theorem 24.6.2 (a)). Under it, with x=Qky, ⟨Ax,x⟩=y∗Qk∗AQky=y∗Tky,⟨x,x⟩=y∗Qk∗Qky=y∗y, so RA(x)=RTk(y) and the two maxima agree. Monotonicity in k follows at once from Kk⊆Kk+1 (Proposition 24.6.1 (a)), and θ1(k)≤λ1 from Proposition 17.1.1 applied on all of Cn — which is also Theorem 24.6.4 (a) at i=1.

Step 2. Let p∈R[x] have deg⁡p≤k−1 and p(λ1)≠0, and put x=p(A)b. Then x∈Kk by Proposition 24.6.1 (a), and x=∑icip(λi)vi since Avi=λivi. The list (vi) being orthonormal, Theorem 11.2.1 gives RA(x)=∑i|ci|2p(λi)2λi∑i|ci|2p(λi)2, and x≠0 because its first coefficient c1p(λ1) is non-zero. Subtracting from λ1, λ1−RA(x)=∑i≥2|ci|2p(λi)2(λ1−λi)∑i|ci|2p(λi)2, the term i=1 having vanished from the numerator. Every factor λ1−λi is at most λ1−λn; and the denominator is at least its single term |c1|2p(λ1)2, all terms being non-negative. Hence (∗)λ1−RA(x) ≤ (λ1−λn)∑i≥2|ci|2p(λi)2|c1|2p(λ1)2 ≤ (λ1−λn)τ2maxi≥2(p(λi)p(λ1))2. By Step 1, λ1−θ1(k)≤λ1−RA(x).

Step 3. Take μ=12(λ2+λn), the midpoint of the interval holding the unwanted eigenvalues, and p(x)=(x−μ)k−1, of degree k−1. For i≥2 we have λn≤λi≤λ2, so |λi−μ|≤12(λ2−λn). And λ1−μ=12(λ2−λn)+(λ1−λ2)>12(λ2−λn)≥0, using λ1>λ2; in particular p(λ1)=(λ1−μ)k−1≠0. Therefore maxi≥2|p(λi)p(λ1)|≤((λ2−λn)/2λ1−μ)k−1=ρk−1, since 2(λ1−μ)=2λ1−λ2−λn. Substituting into (∗) gives the bound, and ρ<1 because λ1>λ2. This proves the theorem.

For A=diag(4,3,2,1) and b=(1,1,1,1)/2 the constants are τ2=3, λ1−λn=3 and ρ=12, so the bound reads 9⋅4−(k−1): the values 9,2.25,0.5625 for k=1,2,3, against the true gaps 1.5, 0.382, 0.068 of Example (Ritz values closing in on the extremes). The bound is correct and loose, and the looseness is in Step 3.

It is worth naming what ρ is. Assume for the rest of this subsection that λ2>λn, which the theorem does not require; if instead λ2=λn, then ρ=0 and the bound of Theorem 24.6.5 is already as strong as this route can make it. Put γ=λ1−λ2λ2−λn, the size of the gap above λ2 measured against the spread of the rest. Then 2λ1−λ2−λn=2(λ1−λ2)+(λ2−λn)=(λ2−λn)(1+2γ), so ρ=11+2γ. (For A=diag(4,3,2,1), γ=12 and ρ=12.) So Step 3 buys a factor (1+2γ)−(k−1) per unit of k, and a small gap makes it nearly worthless.

What is not proved here. The polynomial of Step 3 is not the best one. Minimizing maxi≥2|p(λi)/p(λ1)| over the polynomials of degree at most k−1 with p(λ1)=1 is a question about how small a polynomial can be kept on an interval while staying large at a point outside it, and its answer is the Chebyshev polynomial Tk−1, shifted and scaled to that interval. The resulting estimate, due to Kaniel, Paige and Saad, is 0≤λ1−θ1(k)≤(λ1−λn)τ2(Tk−1(1+2γ))2. The next section defines the Chebyshev polynomials and proves the two facts this bound rests on: that Tk stays within 1 on [−1,1], and that it grows geometrically just outside it; the substitution into (∗) is then immediate, but we do not carry it out here, and nothing in this section depends on it. The reason it is worth stating is the shape of the improvement. For small γ, Tk−1(1+2γ) grows roughly like (1+2γ)k−1, where Step 3 managed only (1+2γ)k−1. A small gap costs a square root, not a full factor. That is the same square root that the next section extracts for conjugate gradients, and it is the reason Krylov methods are used at all.

24.6.6. Losing orthogonality

Everything above is exact arithmetic. In the floating-point model of Definition 24.1.1, one feature of the Lanczos recurrence fails, and it fails in a way that is worth describing honestly rather than proving badly.

The three-term recurrence never re-examines q1,…,qj−2. In exact arithmetic it does not need to: Theorem 24.6.3 proves that wj is automatically orthogonal to them. In floating-point arithmetic the computed vectors satisfy the recurrence to within a small multiple of u‖A‖, and nothing removes the small components along the earlier vectors that rounding introduces. Those components do not stay small.

The mechanism is specific, and it is not random drift. The component of the computed qk+1 along a Ritz vector x grows in inverse proportion to that Ritz pair’s residual βk|yk|. So orthogonality is lost first, and most, in exactly the direction of a Ritz pair that has converged — the better the approximation, the faster the basis is corrupted in its direction. The visible symptom is that the computed Tk acquires a second copy of a Ritz value that converged several steps earlier, and later a third: the spurious duplicates, which correspond to no extra multiplicity in A. An algorithm that counts eigenvalues by counting Ritz values will count wrong.

This paragraph is a description of observed behavior, not a theorem, and this book does not prove it. The analysis that makes it precise — that the loss of orthogonality is confined to the converged directions, and that the computed Tk is the exact tridiagonal matrix of a larger problem — is due to Paige, and it is outside the scope of this chapter. Nothing above depends on it. The standard remedies are to reorthogonalize each new vector against all the earlier ones, which restores the guarantee and abandons the three-term economy, or to reorthogonalize only against the converged Ritz vectors, which is cheaper and rests on precisely the analysis we are not giving. Both are engineering, not theorems, and we label them as such.

24.6.7. Exercises

A. Check your understanding

Exercise (A1)
  • (a) Define Kk(A,b) and state its dimension in terms of mA,b.
  • (b) State the Arnoldi identity of Theorem 24.6.2 (b), naming the shape of every matrix in it.
  • (c) True or false: for Hermitian A, the Arnoldi process needs all of q1,…,qj to produce qj+1. Justify your answer.
  • (d) What does a breakdown hk+1,k=0 mean, and why is it good news?
  • (e) Give the formula for the residual of a Ritz pair, and say why it costs nothing to evaluate.

Solution.

  • (a) Kk(A,b)=span(b,Ab,…,Ak−1b), and dim⁡Kk=min(k,d) with d=deg⁡mA,b (Proposition 24.6.1).

  • (b) AQk=Qk+1H~k, where Qk∈Mn×k and Qk+1∈Mn×(k+1) have orthonormal columns and H~k∈M(k+1)×k is upper Hessenberg; equivalently AQk=QkHk+hk+1,kqk+1ek⊤ with Hk=Qk∗AQk∈Mk upper Hessenberg.

  • (c) False. By Theorem 24.6.3 the recurrence has three terms, so only qj−1 and qj are needed. The reason is that Hk is both Hermitian and upper Hessenberg, hence tridiagonal.

  • (d) It means Kk is A-invariant (Theorem 24.6.2 (c)), equivalently k=d. It is good news because every eigenvalue of Hk is then an exact eigenvalue of A: the method has finished, not failed.

  • (e) ‖Ax−θx‖=hk+1,k|yk| (Theorem 24.6.4 (b)). Both factors are already known — hk+1,k from the last step of the process, yk from the small eigenvector — so no multiplication by A is needed.

B. Practice

Exercise (B1: Dimensions)

For each pair (A,b), determine dim⁡Kk(A,b) for every k≥1. Justify your answer.

  • (a) A=diag(2,2,5), b=(1,1,1).
  • (b) A=diag(2,2,5), b=(1,−1,0).
  • (c) A=(0100), b=e2.
  • (d) A the 4×4 matrix with Aej=ej+1 for j≤3 and Ae4=0, b=e1.

Solution. By Proposition 24.6.1 (b) the answer is always min(k,d) with d=deg⁡mA,b.

  • (a) b has non-zero components in both eigenspaces, so no single factor of (x−2)(x−5) kills it, while the product does: d=2, and the dimensions are 1,2,2,2,….

  • (b) Ab=(2,−2,0)=2b, so b is an eigenvector and mA,b=x−2, d=1. The dimensions are 1,1,1,…: every further product is wasted.

  • (c) Ae2=e1≠0 and A2e2=0, so mA,e2=x2 and d=2: dimensions 1,2,2,….

  • (d) Aje1=ej+1 for j≤3 and A4e1=0, so mA,e1=x4 and d=4: dimensions 1,2,3,4,4,…. Here the Krylov spaces grow as long as they possibly can.

Exercise (B2: Two Arnoldi steps)

Let A=(211021001),b=e3. Carry out two Arnoldi steps. Write down Q2, H~2 and H2, verify Eq. 24.6.1, and compute the two Ritz values with their residuals. Hence say which, if either, is a converged eigenvalue estimate.

Solution. q1=e3 and Aq1=(1,1,1), so h11=⟨Aq1,q1⟩=1 and w1=(1,1,1)−(0,0,1)=(1,1,0),h21=2,q2=12(1,1,0). Next Aq2=12(3,2,0), so h12=⟨Aq2,q1⟩=0 and h22=12(3+2)=52, and w2=12(3,2,0)−52⋅12(1,1,0)=122(1,−1,0),h32=12, with q3=12(1,−1,0). Hence Q2=(01201210),H~2=(10252012),H2=(10252). Check Eq. 24.6.1 column by column: Aq1=(1,1,1)=1⋅q1+2q2, and Aq2=12(3,2,0)=0⋅q1+52q2+12q3, both correct.

H2 is lower triangular as well as upper Hessenberg, so its eigenvalues are θ=1 and θ=52. For θ=52: solving H2y=52y gives y1=0, so y=(0,1), y2=1, and the residual is h32|y2|=12. For θ=1: 2y1+32y2=0, so y=117(3,−22) and the residual is 12⋅2217=3417≈0.343.

Neither is converged. The eigenvalues of A are 2,2,1, so the Ritz value 1 is exactly an eigenvalue — yet its residual is 0.343, not 0, and its Ritz vector Q2y=117(−2,−2,3) is not an eigenvector: the eigenvectors for 1 are the multiples of (0,1,−1), which is not in K2=span(e3,(1,1,0)). The agreement is a coincidence of two numbers, and the residual is what says so.

Exercise (B3: Lanczos on a tridiagonal matrix)

Let A∈Mn(R) be symmetric tridiagonal with every subdiagonal entry non-zero, and let b=e1. Prove that the Lanczos process produces qj=±ej for every j, and hence that Tn=A up to signs. What does this say about the cost of Lanczos on a matrix that is already tridiagonal?

Solution. We show by induction on j that Aje1=∑i≤j+1ci(j)eiwithcj+1(j)=a21a32⋯aj+1,j≠0, for 0≤j≤n−1. At j=0 the sum is e1 and the coefficient is the empty product 1. Suppose it holds at j−1. Since A is tridiagonal, Aei∈span(ei−1,ei,ei+1), so applying A to Aj−1e1=∑i≤jci(j−1)ei lands in span(e1,…,ej+1); and the only term contributing to ej+1 is cj(j−1)Aej, whose ej+1-coefficient is cj(j−1)aj+1,j≠0, both factors being non-zero. This proves the induction.

Consequently the change-of-basis matrix from (e1,…,ej) to (e1,Ae1,…,Aj−1e1) is triangular with non-zero diagonal, so Kj(A,e1)=span(e1,…,ej) for 1≤j≤n.

By Theorem 24.6.2 (a), qj is a unit vector in Kj orthogonal to Kj−1=span(e1,…,ej−1), hence a unit multiple of ej; the sign is fixed by βj−1>0. Therefore (Tn)ij=⟨Aqj,qi⟩=±aij, and since qj=εjej with εj=±1, Tn=DAD with D=diag(ε1,…,εn) — a similarity, so the two have the same eigenvalues.

The cost lesson is that Lanczos does no work at all here beyond rediscovering the matrix: on a tridiagonal matrix with b=e1 it reproduces A and finds nothing early. The method earns its keep when n is large and only the first few extreme eigenvalues are wanted, so that k≪n.

C. Going deeper

Exercise (C1: Ritz values from the other end)

Let A∈Mn(C) be Hermitian and let θ1(k)≥⋯≥θk(k) be the Ritz values from Kk(A,b), 1≤k≤d.

  • (a) Prove that θk(k) is non-increasing in k and that θk(k)≥λn(A).
  • (b) Deduce a bound for θk(k)−λn(A) from Theorem 24.6.5, with no new argument.
  • (c) Show that at k=d the Ritz values are d of the eigenvalues of A, and say which.

Solution.

  • (a) The argument of Step 1 in the proof of Theorem 24.6.5 applies verbatim with minima in place of maxima, using the second half of Proposition 17.1.1: θk(k)=min{RA(x):x∈Kk, x≠0}. A minimum over a larger set is no larger, and Kk⊆Kk+1, so the sequence is non-increasing; and the minimum over all of Cn is λn(A), so it is bounded below by that. (Alternatively, both statements are Theorem 24.6.4 (a) at i=k.)

  • (b) Apply Theorem 24.6.5 to −A. Its eigenvalues are −λn≥−λn−1≥⋯≥−λ1, with the same orthonormal eigenbasis in reversed order. Its Krylov spaces from b are the same subspaces, because (−A)jb=(−1)jAjb and scaling the vectors of a list does not change its span. And R−A=−RA, so by Step 1 of the proof the largest Ritz value of −A from Kk is −θk(k). The theorem gives θk(k)−λn≤(λ1−λn)τ~2ρ~2(k−1),ρ~=λ1−λn−1λ1+λn−1−2λn, with τ~2=∑i≤n−1|ci|2/|cn|2, provided λn−1>λn and cn≠0.

  • (c) At k=d the process breaks down, so Kd is A-invariant and every Ritz value is an eigenvalue of A by Theorem 24.6.2 (c). Which ones: writing b=∑icivi in an orthonormal eigenbasis and grouping the eigenvectors by eigenvalue, p(A)b=∑icip(λi)vi, so Kd is contained in the span of those vi with ci≠0. The Ritz values are therefore the distinct eigenvalues λ of A for which b has a non-zero component in Eλ(A), each occurring once; and d is their number, as Example (A breakdown at step two) illustrates with d=2 for a matrix with eigenvalues 5,2,2.

Exercise (C2: Shift invariance, and what it costs)

Let A∈Mn(F), b≠0, and let σ,c∈F with c≠0.

  • (a) Prove that Kk(cA+σI,b)=Kk(A,b) for every k.
  • (b) Deduce that the Arnoldi vectors qj are the same for A and for cA+σI, up to a scalar of modulus one, and that the Ritz values transform by θ↦cθ+σ.
  • (c) Contrast this with Section 5: explain in two sentences why shifting accelerates the QR algorithm but does not accelerate a Krylov method, and what one must do instead to get the effect.

Solution.

  • (a) Put B=cA+σI. By Proposition 24.6.1 (a), Kk(B,b)={p(B)b:deg⁡p≤k−1} and p(B)=p(cA+σI)=p~(A) where p~(x)=p(cx+σ) has the same degree, c≠0. The map p↦p~ is a bijection of F[x]≤k−1 onto itself, with inverse q↦q((x−σ)/c). So the two sets of vectors coincide.

  • (b) The Arnoldi vectors are determined by the flag K1⊆K2⊆… alone, up to a scalar of modulus one: qj is a unit vector in Kj orthogonal to Kj−1 (Theorem 24.6.2 (a)), and the unit vectors of Kj∩Kj−1⊥ are exactly the modulus-one multiples of one of them. By (a) the flag is the same for B as for A, so the run on B produces QkB=QkD with D=diag(d1,…,dk) and |dj|=1. The scalars need not all be 1: each run fixes its own by hj,j−1>0, and for c=−1, σ=0, A=diag(4,3,2,1) and b=(1,1,1,1) the two runs give qjB=(−1)j−1qj. Consequently (QkB)∗BQkB=D∗(cQk∗AQk+σQk∗Qk)D=D∗(cHk+σIk)D, whose eigenvalues, D being unitary, are those of cHk+σIk, namely cθi+σ.

  • (c) The QR algorithm converges at a rate governed by ratios |λi+1−μ|/|λi−μ|, and those ratios genuinely change with μ; a Krylov space does not change at all under a shift, by (a), so nothing it produces can change either. To get the analogous effect one must change the space, not the matrix: apply the method to (A−μI)−1, whose Krylov spaces are different, which is the Krylov version of the inverse iteration of Theorem 24.4.3 and costs a factorization of A−μI rather than a matrix-vector product.

Exercise (C3: When the Krylov space is everything)

Let A∈Mn(F).

  • (a) Prove that there exists b with Kn(A,b)=Fn if and only if mA=pA.
  • (b) Prove that in that case the Arnoldi process run to k=n produces a unitary Qn with Qn∗AQn upper Hessenberg, all of whose subdiagonal entries are non-zero.
  • (c) Hence give a second proof of Theorem 11.8.5 over C for such A, and explain why it does not replace Chapter 11’s.

Hint for (a): Chapter 10 §05 calls such a b a cyclic vector.

Solution.

  • (a) Kn(A,b)=Fn means dim⁡Kn=n, which by Proposition 24.6.1 (b) means deg⁡mA,b=n, that is, Z(b;A)=Fn and b is a cyclic vector for A (Definition 10.5.4). By Theorem 10.5.6, such a vector exists exactly when mA=pA.

  • (b) With d=n, Theorem 24.6.2 (a) says that hj+1,j≠0 for every j<n and that (q1,…,qn) is an orthonormal basis of Fn; so Qn is unitary, and Qn∗AQn=Hn is upper Hessenberg by Theorem 24.6.2 (b). Its subdiagonal entries are the hj+1,j, all non-zero. (A Hessenberg matrix with this property is called unreduced, and Proposition 24.5.5 is about exactly these.)

  • (c) Part (b) exhibits a unitary Qn with Qn∗AQn upper Hessenberg, which is the conclusion of Theorem 11.8.5. It does not replace that theorem for two reasons. First, it needs mA=pA by (a), while Theorem 11.8.5 holds for every square matrix; when no cyclic vector exists, the process breaks down at some d<n and delivers an invariant subspace instead of a full basis. Second, Theorem 11.8.5 is proved by reflections, each applied to the whole matrix, and its O(n3) cost is what Section 5 assumed; the Arnoldi route needs n matrix-vector products and O(n3) orthogonalization work as well, and in floating-point arithmetic the orthogonality it produces is exactly what the last subsection says can be lost. The two constructions agree in exact arithmetic when both apply, and the reflection route is the one that is run.

Found a typo or a mistake? Report it, or read this page’s source

The text of this book is licensed under CC BY-NC-SA 4.0.