Eigenpairs from Matrix–Vector Products: Power Iteration, Deflation, Lanczos, and Arnoldi
Linear AlgebraEigenvaluesPower IterationLanczos MethodArnoldi MethodMatrix-Free Methods
This document rigorously studies what can be learned about the A nonzero vector u is an eigenvector of A, with eigenvalueλ, when Au=λu. Multiplying u by A only rescales it, by the factor λ, and the pair (λ,u) is called an eigenpair. For example, A=[2112] maps [1,1]⊤ to 3[1,1]⊤ and [1,−1]⊤ to itself, so these two vectors are eigenvectors, with the eigenvalues 3 and 1. Rewriting the definition as (A−λId)u=0 shows that λ is an eigenvalue exactly when the matrix A−λId is singular, that is, when λ is a root of the polynomial λ↦det(λId−A), of degree d, called the characteristic polynomial of A. The eigenvectors of λ are then the nonzero vectors of the null space ker(A−λId), called the eigenspace of λ. A real polynomial can have complex roots, so a real matrix can have complex eigenvalues. The rotation [01−10] turns every real vector by a quarter turn, so it has no real eigenvector, and its characteristic polynomial λ2+1 has the roots ±i. A real symmetric matrix, in contrast, has only real eigenvalues and an orthonormal basis of eigenvectors, by the spectral theorem used in Sections 1 to 4. Computing the characteristic polynomial requires every entry of A, which the setting of this post does not provide. The definition itself, however, involves A only through the product Au. One product is enough to test a candidate pair (λ,u), because the residual Au−λu=(A−λId)u vanishes exactly when the pair is an eigenpair. The methods of this post use A in no other way. For a symmetric matrix, Lemma 2 and Remark 11 also turn a small residual into a guarantee: if ∥u∥=1, then λ is within ∥Au−λu∥ of an eigenvalue. of a real matrix A∈Rd×d that is known only through the map x↦Ax. Storing A takes d2 numbers, and forming it can cost much more. In many applications, however, the product Ax is cheap, because it follows directly from how A is defined, whether A is built from a data matrix, a sparse graph or a discretized differential operator, or is the Hessian of a function whose gradient can be computed (Section 2 gives more examples). Two families of methods need nothing more than such products: power iteration (Section 3), which extracts one eigenpair at a time and reaches further eigenpairs through shifts and deflation, and the Lanczos method (Section 4), which builds an orthonormal basis of a Krylov subspace and solves a small tridiagonal eigenvalue problem inside it. Sections 1 to 4 assume that A is symmetric, as Hessians, covariance matrices and the adjacency matrices of undirected graphs are. The Arnoldi method (Section 5) extends the Lanczos method to nonsymmetric matrices, such as the transition matrix of a random walk on a directed graph.
The results assume neither positive definiteness nor distinct eigenvalues. Each statement says precisely which eigenpairs are recovered, and the remarks say which are not. All the delicate cases are treated explicitly: a dominant eigenvalue of either sign, opposite eigenvalues of equal magnitude, repeated eigenvalues, and zero eigenvalues. Sections 3, 4 and 5 each end with an implementable procedure: Algorithm 1 for power iteration, Algorithm 2 for the Lanczos method and Algorithm 3 for the Arnoldi method. All three stop when a residual computed with the original matrix falls below a tolerance. The figures run these procedures on matrices with prescribed spectra and on five matrices from applications: the adjacency matrix of a graph, a covariance matrix and the Hessian of a neural network (Section 2), and the transition matrix of a random walk on a directed graph and a discretized convection–diffusion operator (Section 5.1). The code that produces them, with implementations of the three algorithms, is on GitHub.
For a symmetric A, the central object is the decomposition of a vector along the eigenspaces of A. Let Λ(A) be the set of distinct eigenvalues of A, and let Πλ be the orthogonal projection onto the eigenspace of λ∈Λ(A). Then, for every vector v∈Rd and every integer k≥0,
v=λ∈Λ(A)∑Πλv,Akv=λ∈Λ(A)∑λkΠλv.
Started from v, power iteration keeps the last vector of the sequence v,Av,A2v,… and lets the factors λk select the eigenvalue of largest magnitude. The Lanczos method keeps the whole sequence and optimizes over its span. Every proof of Sections 3 and 4 follows the components Πλv through these operations, and Section 5 shows what changes when the eigenspaces are no longer orthogonal.
1. Notation and preliminaries
Notation
∥⋅∥ is the Euclidean norm, Id the d×d identity matrix, δa,b the Kronecker delta, equal to 1 if a=b and to 0 otherwise, for a and b in any set (we write δij, without the comma, when i and j are integers), and em=[0,…,0,1]⊤∈Rm the last standard basis vector of Rm.
For a subspace S⊆Rd, S⊥={x∈Rd:x⊤s=0∀s∈S} is its orthogonal complement.
For a matrix M with n columns, ker(M)={x∈Rn:Mx=0} is its null space and range(M)={Mx:x∈Rn} is its column space.
In Sections 1 to 4, A∈Rd×d is symmetric, and Section 5 drops this assumption. We let Λ(A) be the set of its distinct eigenvalues, so that a repeated eigenvalue appears once in Λ(A), and we write λmin(A)=minΛ(A) and λmax(A)=maxΛ(A).
For each λ∈Λ(A), Eλ=ker(A−λId) is the eigenspace of λ, and Πλ is the orthogonal projection onto Eλ. If the columns u1,…,umλ of Uλ∈Rd×mλ form an orthonormal basis of Eλ, where mλ is the multiplicity of λ, then Πλ=UλUλ⊤=∑i=1mλuiui⊤, that is, Πλx=∑i=1mλ(ui⊤x)ui keeps the components of x along these basis vectors. An expression in terms of A alone is given below.
For vectors u1,…,um∈Rd, [u1,…,um]∈Rd×m is the matrix with these columns, and span{u1,…,um} is the set of their linear combinations
diag(c1,…,cm) is the diagonal matrix with diagonal entries c1,…,cm∈R.
N(0,Id) is the standard Gaussian distribution on Rd, and sgn(t)∈{−1,1} is the sign of a real number t=0.
Spectral decomposition
When A is symmetric, the Every real symmetric matrix has an orthonormal basis of eigenvectors. Collecting the basis vectors that share the eigenvalue λ as the columns of a matrix Uλ gives Πλ=UλUλ⊤, and the three identities follow from the orthonormality of the basis. gives
A=λ∈Λ(A)∑λΠλ,Id=λ∈Λ(A)∑Πλ,ΠλΠμ=δλ,μΠλ.
The powers of A act on each eigenspace separately, and so does every polynomial:
for every integer k≥0 and every polynomial p,
Ak=λ∈Λ(A)∑λkΠλ,p(A)=λ∈Λ(A)∑p(λ)Πλ.
In fact, each column ui of Uλ is an eigenvector, Aui=λui, so multiplying Πλ=∑i=1mλuiui⊤ by A gives AΠλ=∑i=1mλ(Aui)ui⊤=λΠλ.
The spectral theorem gives the same identity, since AΠλ=∑μμΠμΠλ=λΠλ2 and Πλ2=Πλ. Applying it k times gives AkΠλ=λkΠλ, and multiplying Id=∑λΠλ by Ak then gives the first identity above.
For a polynomial p(t)=∑k=0ncktk, we have p(A)=∑k=0nckAk=∑λ(∑k=0nckλk)Πλ=∑λp(λ)Πλ, which is the second identity.
Taking for p the Given finitely many distinct numbers, the Lagrange polynomial attached to one of them equals 1 at that number and 0 at all the others. For λ∈Λ(A) it is ℓλ(t)=μ∈Λ(A),μ=λ∏λ−μt−μ. Each factor vanishes at its own μ and equals 1 at t=λ, so ℓλ(λ)=1 and ℓλ(μ)=0 for every other μ∈Λ(A). Its degree is ∣Λ(A)∣−1. Substituting it for p in p(A)=∑νp(ν)Πν keeps only the term ν=λ, so ℓλ(A)=Πλ. For example, if Λ(A)={1,3}, then ℓ1(t)=(3−t)/2 and Π1=(3Id−A)/2. that equals 1 at λ and 0 at the other distinct eigenvalues expresses each projection through A alone:
Πλ=μ∈Λ(A)μ=λ∏λ−μA−μId,
where the empty product, when A has a single distinct eigenvalue, equals Id.
These projections are used in the analysis. The algorithms below do not require computing them.
2. Matrices known only through their products
The results below use A only through the products x↦Ax, so they hold whatever routine computes those products. The following cases are typical.
2.1 Hessians
In many problems, A=∇θ2f(θ) is the Hessian of a function f:Rd→R that is twice continuously differentiable in a neighborhood of a point θ. Examples include the following.
Training losses. The training loss of a (deep) neural network with parameters θ is f(θ)=n1∑i=1nℓ(zi(θ),yi), where zi(θ) is the output of the network on the i-th of n training inputs, yi is the target for that input, and ℓ is a cost function such as ℓ(z,y)=21∥z−y∥2. Its curvature is studied to understand optimization (how training progresses) and generalization (how the network performs on new inputs).
Likelihoods. The negative log-likelihood of observations y1,…,yn under a statistical model with density p(⋅;θ) is f(θ)=−∑i=1nlogp(yi;θ). Its minimizer θ^ is the maximum-likelihood estimate, and the Hessian ∇2f(θ^) is the observed Fisher information. For large samples, its inverse approximates the covariance of θ^, that is, how precisely the data determine the parameters.
Newton-type optimization. To minimize f(θ), Near θ, Newton’s method replaces f by its second-order Taylor expansion, a quadratic function of the step: f(θ+Δθ)≈f(θ)+∇f(θ)⊤Δθ+21Δθ⊤∇2f(θ)Δθ. The gradient of this model with respect to Δθ is ∇f(θ)+∇2f(θ)Δθ, and setting it to zero gives the linear system ∇2f(θ)Δθ=−∇f(θ). When ∇2f(θ) is positive definite, the model is a bowl and Δθ goes to its bottom. Otherwise, the model has no minimum, and practical methods modify the step. In one dimension, the step Δθ=−f′(θ)/f′′(θ) moves to the point where the tangent line of f′ at θ crosses zero, as Newton’s method does to solve f′=0. Write the Hessian in an orthonormal eigenbasis, ∇2f(θ)=∑j=1dλjujuj⊤, with every λj>0. Since ui⊤uj=δij, multiplying the sum on the left by ui⊤ gives ui⊤∇2f(θ)=λiui⊤. Multiplying the linear system on the left by ui⊤ therefore gives λiui⊤Δθ=−ui⊤∇f(θ), that is, ui⊤Δθ=−ui⊤∇f(θ)/λi. Expanding the step in the same basis, Δθ=∑i=1d(ui⊤Δθ)ui, and substituting each coordinate gives Δθ=−i=1∑dλiui⊤∇f(θ)ui. Each component of the gradient is divided by the curvature λi along ui, so the step is short where f is steep and long where f is flat. Gradient descent instead takes the step Δθ=−α∇f(θ)=−i=1∑dα(ui⊤∇f(θ))ui, which multiplies every coordinate by the same factor α. takes the step Δθ that solves the linear system ∇2f(θ)Δθ=−∇f(θ). Solving it by forming and factorizing the d×d Hessian would take d2 numbers of storage and about d3 operations. Truncated Newton methods (Dembo and Steihaug, 1983) instead run a few iterations of the conjugate gradient method (Hestenes and Stiefel, 1952), which uses the Hessian only through products ∇2f(θ)x, computed as explained below, and stop once the step is accurate enough. The conjugate gradient method builds the same Krylov subspaces as the Lanczos method.
Molecular vibrations. The potential energy f(θ) of a molecule depends on the positions θ∈R3N of its N atoms. Near an equilibrium θ⋆, where ∇f(θ⋆)=0, small displacements obey Newton’s second law says that mass times acceleration equals force, and the force on the atoms is minus the gradient of the potential energy, −∇f(θ): the atoms are pushed downhill. Near θ⋆, the first-order Taylor expansion of the gradient is ∇f(θ)≈∇f(θ⋆)+∇2f(θ⋆)(θ−θ⋆)=∇2f(θ⋆)(θ−θ⋆), because ∇f(θ⋆)=0 at an equilibrium. The force is therefore a linear restoring force, as for a system of springs, and Mθ¨=−∇f(θ) becomes the stated equation, where θ¨ is the acceleration, the second time derivative of the positions.Mθ¨=−∇2f(θ⋆)(θ−θ⋆), where M is the diagonal matrix of atomic masses, each repeated for the three coordinates of its atom. The eigenvectors of the In the mass-weighted coordinates x=M1/2(θ−θ⋆), substituting θ−θ⋆=M−1/2x into the equation of motion and multiplying it on the left by M−1/2 gives x¨=−Kx,K=M−1/2∇2f(θ⋆)M−1/2. A motion in which every coordinate oscillates at the same frequency ω in a fixed pattern v, x(t)=cos(ωt)v, has x¨=−ω2x, so it solves the equation exactly when Kv=ω2v. Each eigenvector of K is thus a vibration pattern, a normal mode, and its eigenvalue is the squared frequency of that mode. Because K is symmetric, its orthonormal eigenvectors split every small motion into independent normal modes, each vibrating at its own frequency, and the weighting by M−1/2 on both sides is what keeps the matrix symmetric. At a minimum of f, the Hessian is positive semidefinite, and so is K, since x⊤Kx=(M−1/2x)⊤∇2f(θ⋆)(M−1/2x)≥0. Its eigenvalues are therefore nonnegative, and each one is the square of a real frequency ω. At a saddle point of f, an eigenvalue −μ<0 would have no real square root, and the displacement along its mode would grow like eμt instead of oscillating. For example, take two atoms of masses m1 and m2 at positions θ1 and θ2 on a line, joined by a spring of stiffness k and rest length a. The energy of the spring is f(θ)=2k(θ2−θ1−a)2, and ∇2f(θ⋆)=k[1−1−11],K=k[1/m1−1/m1m2−1/m1m21/m2]. The eigenvalues of K are 0, whose mode moves both atoms together and stretches nothing, and k(1/m1+1/m2), the trace of K. The second mode moves the atoms in opposite directions, the heavier one less, and it is the stretching vibration, with frequency ω=k(1/m1+1/m2).M−1/2∇2f(θ⋆)M−1/2 are the vibration patterns of the molecule (its normal modes), and its eigenvalues are the squared vibration frequencies ω2. A product with this matrix is a Hessian–vector product between two diagonal scalings, M−1/2(∇2f(θ⋆)(M−1/2x)), so the lowest-frequency modes of a large molecule, its slowest collective motions, can be computed without forming the 3N×3N Hessian.
Such a Hessian is symmetric by When the second partial derivatives of f are continuous near θ, they do not depend on the order of differentiation: ∂2f/∂θi∂θj=∂2f/∂θj∂θi. The Hessian is therefore a real symmetric matrix.. In that case, each product can be evaluated as
Ax=∇θ(∇θf(θ)⊤x),
with x held constant during differentiation. The scalar ∇θf(θ)⊤x is the directional derivative of f along x. Its gradient with respect to θ is ∇θ2f(θ)⊤x, which equals Ax because the Hessian is symmetric. Automatic differentiation therefore returns the product at a The product is obtained by differentiating the gradient computation once more, in forward or reverse mode. Each pass costs a small constant multiple of the original computation, so the product costs about as much as a few gradient evaluations. of the cost of a gradient (Pearlmutter, 1994), and the d×d matrix is never formed.
The point θ and everything that defines f must remain fixed during an eigensolver run. For a training loss, evaluating it on a different random subset of the data (a minibatch) between two products changes the matrix being analyzed, and the statements below then no longer describe the computation.
Every symmetric matrix A is in fact a Hessian: the quadratic x↦21x⊤Ax has the Hessian A at every point. This does not help to compute products, since the gradient of that quadratic is Ax itself. The Hessian route pays off when f is given by a program whose gradient is cheap, so that curvature is obtained from first-order code. For the matrices of the following section, the products come instead from the structure that defines A.
2.2 Other matrices known through their products
Products are just as cheap in many cases where A is not presented as a Hessian.
Gram and covariance matrices. For a data matrix X∈Rn×d, whose n rows are observations of d variables, the product X⊤Xx=X⊤(Xx)The parentheses fix the order: two matrix–vector products, and X⊤X is never formed. First, each of the n entries of y=Xx is the inner product of a row of X with x, about d multiplications and additions. Then each of the d entries of X⊤y is the inner product of a column of X with y, about n of them. Each pass costs about nd operations, and only y is stored in between. Forming X⊤X instead takes d2 inner products of columns of length n, about nd2 operations, and d2 numbers of storage. For n=106 observations of d=104 variables, one product costs about 2⋅1010 operations, whereas forming X⊤X costs about 1014 operations and 108 numbers of storage. If X is sparse, a product costs on the order of its number of nonzero entries., whereas forming X⊤X costs O(nd2) and storing it costs O(d2). When each column of X has zero mean, the sample covariance matrix of the d variables is the average of the outer products of the observations, X⊤X/n=(1/n)∑i=1nXiXi⊤, where Xi∈Rd is the i-th row of X. Its leading eigenvector, the first principal component, is the direction of largest variance in the data: it maximizes the variance n1∥Xu∥2 of the data projected on u over unit vectors u.
Gauss–Newton matrices. Many objectives are compositions f(θ)=ψ(F(θ)) of a map F:Rd→Rp and a function ψ:Rp→R, both twice differentiable. For the training loss of Section 2.1, F(θ) stacks the outputs z1(θ),…,zn(θ) into one vector of Rp, and ψ(z)=n1∑i=1nℓ(zi,yi) for z=[z1;…;zn]. Let J∈Rp×d be the Jacobian of F at θ, with entries Jkj=∂Fk(θ)/∂θj, and let H=∇2ψ(F(θ))∈Rp×p. The chain rule gives ∇2f(θ)=J⊤HJ+∑k=1p∂kψ(F(θ))∇2Fk(θ), where ∂kψ is the partial derivative of ψ with respect to its k-th argument. The generalized Gauss–Newton matrix G=J⊤HJ keeps the first term. It is symmetric, and positive semidefinite when ψ is convex. It is applied as Gx=J⊤(H(Jx)), where automatic differentiation computes Jx in forward mode and J⊤u in reverse mode without forming J. For least squares, ψ(z)=21∥z−y∥2 for a target y∈Rp, so H=Ip and G=J⊤J. When F is linear, F(θ)=Jθ, every ∇2Fk vanishes and G is the exact Hessian. When F is nonlinear, G=J⊤J is the classical Gauss–Newton approximation of the Hessian. For the squared loss and the softmax cross-entropy of a neural network, G is also the Suppose the model gives each target a probability density p(y∣z) that depends on the output z=F(θ), and that ℓ is the corresponding negative log-likelihood. The Fisher information matrix is Ey∼p(⋅∣F(θ))[∇θlogp(y∣F(θ))∇θlogp(y∣F(θ))⊤]. By the chain rule it equals J⊤IzJ at z=F(θ), where Iz=Ey∼p(⋅∣z)[∇zlogp(y∣z)∇zlogp(y∣z)⊤]=Ey∼p(⋅∣z)[−∇z2logp(y∣z)]. When −∇z2logp(y∣z) does not depend on y, the last expectation is H, so the Fisher matrix is G. This holds for a Gaussian with mean z (the squared loss, H=I) and for a categorical distribution with probabilities sk=ezk/∑jezj (the softmax cross-entropy, H=diag(s)−ss⊤). of the model.
Graphs. For an undirected graph on d nodes, the adjacency matrix W has the entries Wij=Wji≥0, the weight of the edge between nodes i and j, with Wij=0 when there is no edge. It is symmetric, and so is the Laplacian L=diag(W1)−W, where 1 is the all-ones vector. A product with either matrix costs one pass over the edges. Since x⊤Lx=∑i<jWij(xi−xj)2, the eigenvectors of L with the smallest eigenvalues vary little across heavily weighted edges. They are used to split the graph into well-connected groups, in spectral clustering and graph partitioning (von Luxburg, 2007).
Discretized differential operators. Finite-difference and finite-element discretizations of self-adjoint differential operators, such as those of heat conduction, elasticity or quantum mechanics, are applied point by point or element by element, in time proportional to the number of unknowns and without assembling the matrix. For instance, the operator u↦−u′′ maps a twice-differentiable function u on [0,1] with u(0)=u(1)=0 to the function −u′′. To discretize it, we use the grid ti=ih for i∈{0,…,d+1}, with step h=1/(d+1). A function u is represented by the vector x∈Rd of its values xi=u(ti) at the d interior points ti for i∈{1,…,d}, and the boundary conditions fix x0=xd+1=0, which are not entries of x. At each interior point, the second derivative is replaced by the If u has four continuous derivatives on [0,1], Taylor’s theorem at ti with the increment ±h and the remainder in Lagrange form gives for some ξ+∈(ti,ti+1) and some ξ−∈(ti−1,ti), u(ti±1)=u(ti±h)=u(ti)±hu′(ti)+2h2u′′(ti)±6h3u′′′(ti)+24h4u′′′′(ξ±).Adding the two expansions, rearranging the terms, dividing by h2, and writing xj=u(tj), with x0=xd+1=0 because u(0)=u(1)=0, gives −u′′(ti)=h22xi−xi−1−xi+1+24h2(u′′′′(ξ+)+u′′′′(ξ−)), so the error of the centered difference is at most 12h2max[0,1]∣u′′′′∣.
−u′′(ti)≈h22xi−xi−1−xi+1,i∈{1,…,d},
with an error O(h2). The right-hand side is a linear function of x, so it defines a matrix A∈Rd×d by (Ax)i=(2xi−xi−1−xi+1)/h2 for every x∈Rd and i∈{1,…,d}. When x holds the values of u, the vector Ax therefore holds approximations of the values of −u′′ at the interior points (Figure 1). Row i of A holds the coefficients of xi−1, xi and xi+1 in this formula, except that rows 1 and d have no entry for x0 and xd+1, which are zero:
A=h212−1−12⋱−1⋱−1⋱2−1−12.Figure 1. The discretization of u↦−u′′ of Section 2.2, with d=7 interior grid points. Left: a function u with u(0)=u(1)=0 is represented by its values xi=u(ti) at the interior points (blue). The boundary values x0=x8=0 (orange) are fixed and are not entries of x. The centered difference at t3 combines x2, x3 and x4 with the weights −1, 2 and −1, divided by h2. Right: the same weights form row 3 of the tridiagonal matrix A, so that the product Ax gives the approximations of −u′′ at all the interior points at once. In rows 1 and d, the terms in x0 and xd+1 vanish and have no column.
It is tridiagonal and symmetric, and the formula computes Ax in O(d) operations without storing A. Its eigenpairs approximate those of u↦−u′′. For every integer k≥1, the function u(t)=sin(kπt) satisfies the boundary conditions and −u′′=k2π2u. For k∈{1,…,d}, the vector with entries sin(kπti) is Fix k∈{1,…,d} and set vi=sin(kπti) for i∈{0,…,d+1}. The end values v0=sin(0) and vd+1=sin(kπ) are zero, so the formula for (Av)i holds in every row for v=(v1,…,vd). Since ti±1=ti±h, the identity sin(a−b)+sin(a+b)=2sin(a)cos(b) with a=kπti and b=kπh gives vi−1+vi+1=2cos(kπh)vi, hence (Av)i=h22vi−vi−1−vi+1=h22−2cos(kπh)vi,i∈{1,…,d}. The vector v is nonzero, since 0<kπh<π gives v1=sin(kπh)>0. As cos is decreasing on [0,π], the d eigenvalues (2−2cos(kπh))/h2 obtained for k=1,…,d are distinct, so these are all the eigenpairs of A. Finally, 2−2cos(z)=z2+O(z4), so the eigenvalue equals k2π2+O(k4h2). with eigenvalue (2−2cos(kπh))/h2, which tends to k2π2 as h→0 for fixed k. Lanczos (1950) designed his method for linear differential and integral operators.
Structured matrices. A symmetric Toeplitz matrix is constant along each diagonal, Aij=a∣i−j∣ for numbers a0,…,ad−1. The covariance matrix of d consecutive values X1,…,Xd of a stationary time series has this form, with ak=Cov(Xi,Xi+k), because stationarity means that Cov(Xi,Xj) depends only on ∣i−j∣. A product with A costs O(dlog(d)) operations instead of O(d2): A is the upper-left block of a A circulant matrix C∈Rn×n is determined by a vector c=(c0,…,cn−1) through Cij=c(i−j)modn, so each row is the previous one shifted cyclically by one position. Take n=2d and c=(a0,a1,…,ad−1,0,ad−1,…,a1). For i,j∈{1,…,d}, the index (i−j)mod2d equals i−j when i≥j and 2d−(j−i) when i<j, and in both cases the entry of c at that index is a∣i−j∣. Hence A is the upper-left d×d block of C, and Ax consists of the first d entries of C[x;0], with 0∈Rd. The product (Cy)i=∑jc(i−j)modnyj is the cyclic convolution of c and y, and the discrete Fourier transform F of length n turns a cyclic convolution into an entrywise product, F(Cy)=(Fc)⊙(Fy). Therefore Ax=first d entries of F−1((Fc)⊙F[x;0]). The vector Fc is computed once. Each product then takes one transform and one inverse transform of length 2d, each O(dlog(d)) operations with the fast Fourier transform, and 2d multiplications., whose rows are cyclic shifts of one another, and the fast Fourier transform (Cooley and Tukey, 1965) multiplies such a matrix by a vector in O(dlog(d)) operations. Kronecker products arise on grids. For M1∈Rd1×d1 and M2∈Rd2×d2, the Kronecker product M1⊗M2 is the d1d2×d1d2 matrix made of d1×d1 blocks of size d2×d2, with block (k,l) equal to (M1)klM2. It is symmetric when M1 and M2 are, since (M1⊗M2)⊤=M1⊤⊗M2⊤. It is applied without being formed, through the Write yl∈Rd2 for column l of Y, so that vec(Y)=[y1;…;yd1]. Block row k of M1⊗M2 is [(M1)k1M2,…,(M1)kd1M2], so block k of the product is ∑l=1d1(M1)klM2yl=M2∑l=1d1yl(M1⊤)lk. The last sum is column k of YM1⊤, so this block is column k of M2YM1⊤. Computing M2Y costs about d22d1 multiplications, and multiplying the result by M1⊤ about d2d12, whereas the formed matrix has d12d22 entries.
(M1⊗M2)vec(Y)=vec(M2YM1⊤),Y∈Rd2×d1,
where vec(⋅) stacks the columns of a matrix into one vector. The right-hand side costs O(d1d2(d1+d2)) operations, whereas a product with the formed matrix costs O(d12d22). When Y holds the values of a function on a grid with d2 points along the first axis and d1 along the second, M2Y applies M2 to each column, along the first axis, and YM1⊤ applies M1 to each row, along the second axis. For instance, take the operator u↦−∂2u/∂s2−∂2u/∂t2 on the unit square, with u=0 on its boundary, and the grid points (si,tj)=(ih,jh) for i,j∈{1,…,d}, with h=1/(d+1). With Yij=u(si,tj), the centered difference of the previous item along each axis gives −∂2u/∂s2(si,tj)≈(AY)ij and −∂2u/∂t2(si,tj)≈(YA)ij, where A is the matrix of the previous item. The discretized operator is therefore vec(Y)↦vec(AY+YA), that is, the matrix Id⊗A+A⊗Id of size d2, whose product costs O(d2) operations because A is tridiagonal.
Shifted inverses. Let σ be a number that is not an eigenvalue of A, so that A−σId is invertible. Suppose that linear systems (A−σId)y=x can be solved, for instance with a sparse factorization of A−σId into triangular factors, computed once from the entries of A, after which each solve consists of two triangular solves. The solution is y=Bx with B=(A−σId)−1, so each solve is a product with B, a symmetric matrix as the inverse of a symmetric one. It has the eigenvectors of A. The eigenvalue μ=1/(λ−σ) of B is largest in absolute value for the eigenvalue λ of A closest to σ. Power iteration on B therefore returns this μ (Remark 4), unless another eigenvalue of A lies at the same distance on the other side of σ (Remark 5), and λ=σ+1/μ. When λmin<σ<λmax, the extreme eigenvalues of B, which the extreme Ritz values of the Lanczos method approach from inside (Remark 14), are 1/(λ+−σ)>0 and 1/(λ−−σ)<0, where λ+ and λ− are the eigenvalues of A closest to σ above and below it. This is the shift-and-invert strategy. It reaches eigenvalues inside the spectrum, which power iteration on a shifted matrix A−τId never returns, since it returns λmax or λmin, whichever is farther from τ (Corollary 1).
All these matrices are symmetric, and the results of Sections 3 and 4 rely on symmetry through the spectral theorem. For a nonsymmetric matrix accessed through products, such as the transition matrix of a random walk on a directed graph, whose entry (i,j) is the probability of stepping from node i to node j, the Arnoldi method (Arnoldi, 1951) replaces the Lanczos method. Section 5 presents it, together with nonsymmetric examples and what power iteration retains without symmetry.
3. Power iteration
In this part, we first establish the convergence conditions for ordinary power iteration, then describe how shifting and deflation recover additional eigenpairs. Starting from a unit vector v0, power iteration alternates a product with A and a normalization:
vk+1=∥Avk∥Avk,k≥0.
If Av0=0, then v0 is already an eigenvector with eigenvalue zero. In this case, the algorithm should return (0,v0) without attempting to normalize Av0.
The numerical analyst Alston S. Householder attributed the first treatment of the method to Müntz (1913) in his book The Theory of Matrices in Numerical Analysis (1964), as recalled in the historical survey of Golub and van der Vorst (2000). Deflation schemes built on it were analyzed early on, for example by Osborne (1958). The eigenvalue estimate of power iteration is the For a unit vector v, the Rayleigh quotient v⊤Av is the curvature of x↦21x⊤Ax along v. It lies between λmin(A) and λmax(A). It equals λ for every unit vector v of the eigenspace Eλ, since Av=λv gives v⊤Av=λ∥v∥2=λ. The converse fails: v⊤Av=λ does not make v an eigenvector, except when λ is λmin(A) or λmax(A).vk⊤Avk.
3.1 Proposition 1 [Convergence of power iteration]
Let v0∈Rd be a unit vector, and define
Λ0={λ∈Λ(A):Πλv0=0}.
Thus, Λ0 contains precisely the eigenvalues whose eigenspaces have a nonzero component in the initialization. Suppose that Λ0 contains a nonzero eigenvalue λ⋆ satisfying ∣λ⋆∣>∣λ∣ for every λ∈Λ0∖{λ⋆}. The eigenvalue λ⋆ may have arbitrary multiplicity. Then the iteration
vk+1:=∥Avk∥Avk=∥Ak+1v0∥Ak+1v0,k≥0,
is well defined, and
(sgn(λ⋆))kvk⟶q⋆:=∥Πλ⋆v0∥Πλ⋆v0.
In particular,
Aq⋆=λ⋆q⋆,vk⊤Avk⟶λ⋆.
Remark 1 [Why the dominant component wins]
Recall that for every vector v∈Rd and every integer k≥0,
v=∑λ∈Λ(A)Πλv and
Akv=∑λ∈Λ(A)λkΠλv.
So, for the iteration above, each application of A multiplies the component of the current vector in Eλ by λ, so after k steps the component in Eλ carries the factor λk:
Among the components present in v0, the blue one, with the largest ∣λ∣ (we have ∣λ⋆∣>∣λ∣ for every λ∈Λ0∖{λ⋆} by assumption), grows fastest, and the normalization removes the common factor ∣λ⋆∣k, so the orange components fade relative to it. Only the sign of λ⋆k survives, which is why the statement multiplies vk by (sgn(λ⋆))k. The limit is the normalized projection of v0 onto the whole eigenspace Eλ⋆: inside a repeated eigenspace, the direction is chosen by the initialization, not by the method.
Proof
Set pλ=Πλv0. The spectral decomposition gives
v0=λ∈Λ0∑pλ,Akv0=λ∈Λ0∑λkpλ(k≥1).
Well-definedness. Because pλ⋆=0 and λ⋆=0, the component
Πλ⋆Akv0=λ⋆kpλ⋆
is nonzero for every k≥0. Consequently, Akv0=0 for every k≥0. Since ∥v0∥=1, repeated normalization satisfies
vk=∥Akv0∥Akv0.
Indeed, this identity holds at k=0. If it holds at iteration k, then
The blue term is the component along the dominant eigenspace, and it keeps its size after rescaling. The orange term collects every other component, and it shrinks geometrically. Since pλ⋆∈Eλ⋆, the limit q⋆ satisfies Aq⋆=λ⋆q⋆.
Rayleigh quotient. For k≥1, the symmetry of A and its spectral decomposition A=∑λ∈Λ(A)λΠλ give
The second line uses (Akv0)⊤=v0⊤Ak, since A is symmetric. The third expands A2k+1 and A2k with Aj=∑λλjΠλ. The fourth uses v0⊤Πλv0=v0⊤Πλ⊤Πλv0=∥Πλv0∥2=∥pλ∥2, which holds because Πλ is symmetric and Πλ2=Πλ, and only the eigenvalues of Λ0 remain, since pλ=0 for the others. The fifth divides the numerator and the denominator by λ⋆2k>0.
■□
Remark 2 [A negative dominant eigenvalue]
If λ⋆>0, Proposition 1 gives vk→q⋆. If λ⋆<0, it gives (−1)kvk→q⋆, or, equivalently,
v2k⟶q⋆ and v2k+1⟶−q⋆.
Thus, in the latter case, the even and odd subsequences converge to opposite unit vectors. The raw sequence {vk}k≥0 does not converge, but min{∥vk−q⋆∥,∥vk+q⋆∥}⟶0.
A stopping rule that compares consecutive iterates, such as ∥vk+1−vk∥≤ε with ε<2, is therefore eventually never satisfied when λ⋆<0, because ∥vk+1−vk∥→2∥q⋆∥=2, even though the eigendirection has converged. The Rayleigh quotient converges in both cases. In panel (a) of Figure 2, where λ⋆=−1, the raw distance ∥vk−q⋆∥ (grey) alternates between values close to 0 and values close to 2, while ∥(−1)kvk−q⋆∥ (blue) tends to zero.
3.2 Theorem 1 [Rate of convergence]
Assume the hypotheses of Proposition 1, and write pλ⋆=Πλ⋆v0. Let ρ2=max{∣λ∣:λ∈Λ0∖{λ⋆}},
with ρ2=0 when Λ0={λ⋆}, and let ϕ0∈[0,π/2) be Write v0=pλ⋆+(v0−pλ⋆), where pλ⋆=Πλ⋆v0 lies in Eλ⋆ and the remainder is orthogonal to it. The three vectors form a right triangle: the hypotenuse v0 (black), a side pλ⋆ along Eλ⋆ (blue), and a side v0−pλ⋆ perpendicular to it (orange). The angle between a vector and a subspace is the angle between the vector and its orthogonal projection, the smallest angle it makes with a vector of the subspace, so ϕ0 is the angle between v0 and pλ⋆. The tangent is the opposite side over the adjacent side, and since ∥v0∥=1, also cos(ϕ0)=∥pλ⋆∥ and sin(ϕ0)=∥v0−pλ⋆∥. The hypothesis pλ⋆=0 keeps ϕ0 below π/2., so that tan(ϕ0)=∥v0−pλ⋆∥/∥pλ⋆∥. Then, for every k≥1,
Relative to the dominant component, each product multiplies the component in Eλ by λ/λ⋆ (Remark 1), so the slowest of them, with ∣λ∣=ρ2, sets the pace. The eigenvector error therefore shrinks by the factorρ2/∣λ⋆∣ per product, and the eigenvalue error by its square, because the Rayleigh quotient is stationary at an eigenvector. In panel (a) of Figure 2, the factor is 0.8, and the two errors decay like 0.8k and 0.82k. The factor can be arbitrarily close to one. Remark 5 below is the limiting case ρ2=∣λ⋆∣, in which convergence can fail altogether.
This factor is related to the condition number κ(A) of A in the following way.
For an invertible A, the For a matrix M, ∥M∥=max∥z∥=1∥Mz∥ is the largest factor by which M stretches a vector, and the condition number of an invertible A is κ(A)=∥A∥∥A−1∥. For a symmetric A, write a unit vector along the eigenspaces, z=∑λΠλz, so that Az=∑λλΠλz. The eigenspaces are orthogonal, so Pythagoras gives ∥Az∥2=∑λλ2∥Πλz∥2≤(maxλλ2)∑λ∥Πλz∥2=maxλλ2, with equality when z is a unit eigenvector of an eigenvalue of largest magnitude. Hence ∥A∥=maxλ∈Λ(A)∣λ∣. If A is invertible, 0 is not an eigenvalue, since Au=0 forces u=0, and A−1=∑λλ−1Πλ, because multiplying it by A=∑μμΠμ gives ∑λΠλ=Id. So A−1 is symmetric, with the same eigenspaces and the eigenvalues 1/λ, and the same argument gives ∥A−1∥=maxλ1/∣λ∣=1/minλ∈Λ(A)∣λ∣. Multiplying the two norms, κ(A)=minλ∈Λ(A)∣λ∣maxλ∈Λ(A)∣λ∣. Without symmetry, ∥A∥ and 1/∥A−1∥ are the largest and smallest singular values of A, which can differ greatly from the magnitudes of its eigenvalues.κ(A)=maxλ∈Λ(A)∣λ∣/minλ∈Λ(A)∣λ∣ compares the largest magnitude with the smallest, while the factor compares the largest with the next one. The two are linked only by a lower bound. When Λ0 contains an eigenvalue other than λ⋆, ∣λ⋆∣≤∥A∥:=max∥z∥=1∥Az∥ and ρ2≥minλ∈Λ(A)∣λ∣=1/∥A−1∥, so
ρ2/∣λ⋆∣≥1/κ(A),
with equality when Λ0=Λ(A) has exactly two elements. A well-conditioned matrix is therefore a slow one for power iteration. When κ(A) is close to one, all the eigenvalues have nearly the same magnitude, and the factor is close to one as well. For A=diag(1.01,1), κ(A)=1.01, and each decimal digit of the eigenvector costs about One decimal digit of accuracy means dividing the error by 10. By Theorem 1, after k products the eigenvector error is at most (ρ2/∣λ⋆∣)ktan(ϕ0), here (1/1.01)ktan(ϕ0). Each product multiplies this bound by 1/1.01, so k products multiply it by (1/1.01)k. The bound has thus gained one digit when (1.011)k=101⟺k=ln1.01ln10≈0.009952.3026≈231. The logarithm solves for the unknown exponent, and the ratio does not depend on its base, since ln10/ln1.01=1/log101.01. For a factor 1−δ, where δ=1−ρ2/∣λ⋆∣ is the relative gap, the same computation gives k=−ln(1−δ)ln10≈δ2.3 products per digit when δ is small, because −ln(1−δ)≈δ. The cost grows like 1/δ. Here δ≈0.0099, and the approximation gives about 232.. This is the reverse of the conjugate gradient method for Ax=b (Hestenes and Stiefel, 1952), which converges fast when κ(A) is small. A large condition number, on the other hand, implies nothing. diag(1,10−6) has κ=106 and the factor 10−6, so each product gains six digits, whereas diag(1.01,1,10−6) also has κ≈106 but is as slow as diag(1.01,1).
The pace is set by the relative gap 1−ρ2/∣λ⋆∣ at the top of the spectrum, which depends only on the two largest magnitudes present.
For a singular A, κ(A)=∞ by convention, and the bound above says nothing. It still holds with the zero eigenvalue left out. The pseudo-inverse A+=∑λ∈Λ(A)∖{0}λ−1Πλ inverts A on range(A) and vanishes on E0=ker(A), and
κ+(A)=∥A∥∥A+∥=maxλ∈Λ(A)∣λ∣/minλ∈Λ(A)∖{0}∣λ∣
compares the largest magnitude with the smallest nonzero one. When Λ0 contains a nonzero eigenvalue other than λ⋆, ρ2≥minλ∈Λ(A)∖{0}∣λ∣=1/∥A+∥, so ρ2/∣λ⋆∣≥1/κ+(A), and everything said above about κ(A) applies to κ+(A). For instance, diag(1.01,1,0) has κ+=1.01 and is as slow as diag(1.01,1). Otherwise, the only other eigenvalue that can be present is 0, so ρ2=0, and a single product gives the eigenvector exactly. In both cases, a zero eigenvalue never slows power iteration. Its magnitude is below every other, and the first product removes the component along E0.
Proof
Recall from the proof of Proposition 1 the scaled remainder and the normalized iterate:
The vectors pλ lie in mutually orthogonal eigenspaces, so the squared norm of a combination of them is the sum of the squared norms of its terms (Pythagoras). Every ∣λ∣ in the sum defining rk is at most ρ2, hence
where the last equality is Pythagoras again, applied to v0−pλ⋆=∑λ=λ⋆pλ.
Write tk=∥rk∥/∥pλ⋆∥, so that tk≤(ρ2/∣λ⋆∣)ktan(ϕ0).
Since rk is orthogonal to pλ⋆ and q⋆=pλ⋆/∥pλ⋆∥, the inner product of the unit vector (sgn(λ⋆))kvk with q⋆ is
For unit vectors w and q, ∥w−q∥2=∥w∥2−2w⊤q+∥q∥2=2−2w⊤q, so
(sgn(λ⋆))kvk−q⋆2=2−1+tk22≤tk2,
because 1/1+t≥1−t/2 for all t≥0 (the convex function t↦1/1+t lies above its tangent line at 0). Taking square roots and using the bound on tk gives the first inequality. For the Rayleigh quotient, recall equation (1) from the proof of Proposition 1, whose denominator equals ∥pλ⋆∥2+∥rk∥2 by the expression of ∥rk∥2 above:
The first inequality moves the absolute value inside the sum (triangle inequality). The second uses ∣λ−λ⋆∣≤∣λ∣+∣λ⋆∣<2∣λ⋆∣, since every other eigenvalue present satisfies ∣λ∣<∣λ⋆∣, and makes the denominator smaller by dropping ∥rk∥2. The two equalities are the expression of ∥rk∥2 above and the definition of tk, and the last inequality squares the bound on tk. This is the second inequality of the theorem.
■□
Remark 4 [Which eigenvalue is found]
The hypothesis of Proposition 1 concerns distinct eigenvalues present in the initialization. It does not require the dominant eigenvalue to be A simple eigenvalue is an eigenvalue that appears only once as a root of a matrix’s characteristic polynomial λ↦det(λId−A), that is, with algebraic multiplicity one. For a symmetric matrix, the number of times λ appears as a root equals the dimension mλ of its eigenspace (geometric multiplicity), so λ is simple exactly when mλ=1. Its eigenspace Eλ is then a line, and its unit eigenvector is unique up to sign. For example, in diag(3,3,1), the eigenvalue 1 is simple, with E1=span{[0,0,1]⊤}, whereas 3 is a double root, with E3=span{[1,0,0]⊤,[0,1,0]⊤}. Power iteration started from v0=[a,b,c]⊤, with (a,b)=(0,0), converges to a2+b2[a,b,0]⊤, an eigenvector of 3 that depends on v0, so different initializations give different eigenvectors. When the dominant eigenvalue is simple, every admissible initialization gives the same eigenvector, up to sign., since the limit is the normalized projection of v0 onto its entire eigenspace. To obtain an eigenvalue of globally largest absolute value, the initialization must have a nonzero projection onto a corresponding eigenspace. For example, if
v0=g/∥g∥ with
g∼N(0,Id),
then Πλv0=0 for every λ∈Λ(A) with probability one. Indeed, for each λ, the event Πλg=0 means that g belongs to the A subspace of Rd is proper when it is not all of Rd, that is, when its dimension is less than d. Here Eλ contains an eigenvector u=0 of λ, and u⊤u=∥u∥2>0, so u is not orthogonal to itself and u∈/Eλ⊥. Hence Eλ⊥=Rd, and indeed dimEλ⊥=d−mλ≤d−1<d.Eλ⊥, which has A proper subspace S lies in a hyperplane {x:c⊤x=0} with c=0. For g∼N(0,Id), the number c⊤g is Gaussian with variance ∥c∥2>0, so it equals zero with probability zero.. There are only finitely many eigenspaces. Consequently, such an initialization satisfies the hypotheses of Proposition 1 whenever A=0 and exactly one distinct eigenvalue of A has maximal absolute value.
Remark 5 [Failure for opposite eigenvalues of equal dominant magnitude]
The condition examined here is the strict dominance assumption of Proposition 1. In words, some nonzero eigenvalue present in the initialization is strictly larger in absolute value than every other eigenvalue present. With Λ0={λ∈Λ(A):Πλv0=0}, the set of eigenvalues present, it reads
∃λ⋆∈Λ0∖{0}∀λ∈Λ0∖{λ⋆}:∣λ⋆∣>∣λ∣.
The condition The matrix A=[0110] swaps the two coordinates of a vector. Its eigenvalues are 1, with the eigenvector [1,1]⊤, and −1, with the eigenvector [1,−1]⊤. From v0=[1,0]⊤=21[1,1]⊤+21[1,−1]⊤, both are present, and each is matched in absolute value by the other. Power iteration gives v1=[0,1]⊤,v2=[1,0]⊤,v3=[0,1]⊤,… The iterates alternate between two vectors, neither of which is an eigenvector, and the Rayleigh quotient vk⊤Avk=0 is constant but is not an eigenvalue. The condition also fails for diag(2,−2,1) when the first two coordinates of v0 are nonzero: 1 is exceeded by 2, and 2 and −2 match each other. It holds for the same matrix when v0=[1,0,1]⊤/2, because −2 is then absent and 2 dominates 1. It also holds for diag(2,2,1), because 2 counts as one distinct eigenvalue, however many times it is repeated. when no eigenvalue present stands out in this way: every nonzero eigenvalue present is matched or exceeded in absolute value by another eigenvalue present, i.e.
∀μ∈Λ0∖{0}∃λ∈Λ0∖{μ}:∣λ∣≥∣μ∣.
Suppose that the condition fails, while some nonzero eigenvalue is present in the initialization. Let ρ=max{∣λ∣:λ∈Λ0} be the largest absolute value among the eigenvalues present, so that ρ>0. If only one element of Λ0 had absolute value ρ, it would be strictly larger in absolute value than every other eigenvalue present, and the condition would hold with λ⋆ equal to that element. Since the condition fails, at least two distinct elements of Λ0 have absolute value ρ, +ρ and −ρ, and both are present in the initialization. We show that ordinary power iteration can then fail to converge to any eigendirection, even while its Rayleigh quotient converges.
Suppose that
p+=Πρv0=0,p−=Π−ρv0=0;
and that every other eigenvalue present in v0 has absolute value smaller than ρ. For k≥1, define the scaled remainder
rk=λ∈Λ0∣λ∣<ρ∑(ρλ)kpλ.
Then
∥rk∥2=λ∈Λ0∣λ∣<ρ∑ρλ2k∥pλ∥2⟶0,
since every ratio ∣λ/ρ∣ in this sum is smaller than one, and
The second line uses the spectral decomposition Akv0=ρkp++(−ρ)kp−+∑∣λ∣<ρλkpλ, and its denominator follows from Pythagoras, since the three terms in the parentheses lie in Eρ, in E−ρ and in the sum of the other eigenspaces, which are mutually orthogonal. It is positive because p+=0.
The blue component keeps its sign, the orange component changes sign at every step, and the remainder vanishes. Hence
v2k⟶∥p+∥2+∥p−∥2p++p−,
whereas
v2k+1⟶∥p+∥2+∥p−∥2p+−p−.
These two limits are not collinear. Indeed, if p++p−=c(p+−p−) for some scalar c, projecting onto Eρ gives c=1, whereas projecting onto E−ρ gives c=−1. Thus no choice of an overall factor in {−1,1} can make the iterates converge to a single unit eigenvector.
Nevertheless, the Rayleigh quotient satisfies
This limit need not be an eigenvalue. For example, let
A=[100−1],v0=21[11].
Then, at every iteration,
vk=21[1(−1)k],vk⊤Avk=0.
However,
0∈/Λ(A),Avk−(vk⊤Avk)vk=∥Avk∥=1.
Thus a constant or convergent Rayleigh quotient does not by itself establish convergence to an eigenpair. This obstruction motivates the squared-power construction of Proposition 2 below and the use of an eigenpair residual as a stopping criterion. Panel (b) of Figure 2, below, shows this obstruction on a matrix whose eigenvalues +1 and −1 tie for the largest magnitude: the residual ∥Avk−(vk⊤Avk)vk∥ of plain power iteration (orange) stays of order one.
Figure 2. Power iteration on two symmetric matrices of size d=60 with random eigenvectors. Each curve is the median over 200 Gaussian initializations, and its band contains the middle half of the runs (the interquartile range). (a) The dominant eigenvalue λ⋆=−1 has multiplicity 3 and the next largest magnitude is ρ2=0.8. The sign-corrected error ∥(−1)kvk−q⋆∥ (blue) decays like 0.8k and the Rayleigh-quotient error (orange) like 0.82k, below the bounds of Theorem 1 (dashed), while the raw distance ∥vk−q⋆∥ (grey) alternates between those small values and nearly 2 (Remark 2). (b) The eigenvalues +1 and −1 tie for the largest magnitude. The residual ∥Av−(v⊤Av)v∥ of plain power iteration (orange) stalls at order one (Remark 5), whereas squared power iteration with sign separation (blue; Proposition 2, run as in Algorithm 1) drives the residual of its candidate to zero, by a factor of about 0.8 per product. Both horizontal axes count matrix–vector products: vk costs k products and its Rayleigh quotient one more, and each squared step costs two.
3.3 Corollary 1 [Recovering algebraic extrema by shifting]
Let λmin=minΛ(A) and λmax=maxΛ(A).
If α≥λmax and α>λmin, then power iteration applied to B=αId−A converges to a unit vector in Eλmin, provided its initialization has a nonzero projection onto that eigenspace. Its Rayleigh quotient converges to μmax=α−λmin, so λmin=α−μmax.
Similarly, if β≤λmin and β<λmax, then power iteration applied to C=A−βId converges to a unit vector in Eλmax, provided its initialization has a nonzero projection onto that eigenspace. Its Rayleigh quotient converges to νmax=λmax−β, so λmax=β+νmax.
Remark 6 [What a shift does to the spectrum]
Power iteration returns the eigenvalue of largest magnitude, which is λmax or λmin depending on the signs. A scalar shift moves the whole spectrum without changing the eigenvectors, and the shift αId−A also reverses its order. With α≥λmax, the spectrum of B lies in [0,α−λmin], so the eigenvalue of A farthest below α becomes the dominant eigenvalue of B.
Proof
For any eigenpair (λ,u) of A,
Bu=(αId−A)u=(α−λ)u.
Since an orthonormal eigenbasis of A is also an eigenbasis of B,
Λ(B)={α−λ:λ∈Λ(A)}.
The assumptions on α give
0≤α−λ≤α−λmin and
α−λmin>0.
Whenever λ>λmin, α−λ<α−λmin.
Therefore, B has a positive dominant eigenvalue α−λmin, whose eigenspace is exactly Eλmin. Proposition 1 applies, including when this eigenvalue is repeated.
For C,
Cu=(A−βId)u=(λ−β)u
and Λ(C)={λ−β:λ∈Λ(A)}.
The assumptions on β give
0≤λ−β≤λmax−β
and λmax−β>0.
Whenever λ<λmax, λ−β<λmax−β.
Thus C has a positive dominant eigenvalue λmax−β, with eigenspace Eλmax. Proposition 1 again applies.
■□
In particular, choosing α=λmax recovers the minimum eigenvalue from the maximum, unless A=λmaxId. In that exceptional case, every eigenvalue already equals λmax. The analogous statement holds for β=λmin.
Remark 7 [The shift only needs to pass the midpoint]
The condition α≥λmax is convenient but stronger than needed. Proposition 1 only requires that α−λmin strictly dominate ∣α−λ∣ for every other eigenvalue λ present in the initialization. For λ≤α, this holds because λ>λmin. For λ>α, it reads λ−α<α−λmin, that is, λ+λmin<2α. Hence every shift α>(λmin+λmax)/2 works. In particular, an underestimate of λmax, such as a Rayleigh quotient
is a valid shift as long as it exceeds the midpoint of the spectrum. Symmetrically, every β<(λmin+λmax)/2 works for C.
The midpoint is unknown when λmin and λmax are, so x must be chosen without it. A random vector is not enough. For g∼N(0,Id), the unit vector u=g/∥g∥ is uniformly distributed on the sphere, so E[uu⊤]=Id/d and
E[g⊤Ag/∥g∥2]=tr(AE[uu⊤])=tr(A)/d,
the mean of the eigenvalues counted with multiplicity, which can lie far below the midpoint. For the spectrum {10,0,…,0}, the midpoint is 5, while the mean is 10/d.
Power iteration itself supplies a suitable x. Assume that one of the two extremes is strictly larger in absolute value than the other, so that, with a random initialization (Remark 4), power iteration on A converges to it. Take for x the iterate vk, and for the shift its Rayleigh quotient θk=vk⊤Avk. By Theorem 1, θk converges to the dominant eigenvalue λ⋆∈{λmin,λmax}, with
∣θk−λ⋆∣≤2∣λ⋆∣(ρ2/∣λ⋆∣)2ktan2(ϕ0).
Both extremes lie at the distance (λmax−λmin)/2 from the midpoint. As soon as ∣θk−λ⋆∣ is smaller than this distance, θk lies on the same side of the midpoint as λ⋆, and θk is a valid shift. The sign of λ⋆ tells which extreme was found. If λ⋆>0, every other eigenvalue satisfies λ≤∣λ∣<λ⋆, so λ⋆=λmax. If λ⋆<0, every other eigenvalue satisfies λ≥−∣λ∣>λ⋆, so λ⋆=λmin. This sign can be read from θk. In both cases, the midpoint lies strictly between 0 and λ⋆. For instance, if λ⋆=λmax>0, then ∣λmin∣<λmax gives λmin+λmax>0. Hence a θk on the same side of the midpoint as λ⋆ has the sign of λ⋆. The sign also shows in the iterates: vk⊤vk+1=θk/∥Avk∥, so consecutive iterates point to opposite sides exactly when θk<0, which is the alternation of Remark 2. If θk>0, the matrix to use is B=θkId−A of Corollary 1, with α=θk. If θk<0, it is C=A−θkId, with β=θk. Even this distinction is a convenience. The two matrices differ only by their sign, so they have the same eigenvectors, and power iteration on either, from a random initialization, converges, up to sign (Remark 2), to an eigenvector of the eigenvalue of A farthest from θk, which is the other extreme.
Products alone cannot certify that θk has passed the midpoint, because the spread λmax−λmin is unknown. A failure is nevertheless visible. If θk is still on the wrong side of the midpoint, the eigenvalue of A farthest from θk is λ⋆ itself, so the second run converges to an eigenvector of λ⋆, the eigenvalue already found by the first run, instead of the other extreme, and comparing the eigenvalues returned by the two runs reveals it. The remedy is to continue the first run, which brings θk closer to λ⋆, and to repeat the second. Once ∣θk−λ⋆∣<(λmax−λmin)/2, the second run reaches the other extreme. This can only fail if the half-spread is below the accuracy reached, that is, if every eigenvalue of A is within twice that accuracy of λ⋆, and then the other extreme coincides with λ⋆ to that accuracy anyway. Algorithm 1 avoids the question altogether: it uses no shift, reaches the remaining eigenpairs by deflation, and separates eigenvalues of opposite signs by sign separation.
Beyond the two extremes, however, a scalar shift does not help. It changes every eigenvalue, and the dominant eigenvalue of A−σId comes from the eigenvalue of A farthest from σ, which is λmin or λmax whatever σ is. No shift therefore makes an interior eigenvalue dominant. Deflation does. Power iteration finds one eigenpair, the dominant one, and running it again on A would only find the same pair. To reach the next pair, deflation removes from A the eigendirections already found and leaves every other eigenpair untouched. In an orthonormal eigenbasis q1,…,qd of A, with Aqi=λiqi, the spectral decomposition writes A as a sum of rank-one layers, one per eigenvector, and deflating the first r eigenvectors means deleting their layers, which gives a new matrix, the deflated matrix B:
This matrix has the same eigenvectors as A. The ones already found (orange) now have the eigenvalue 0, and all the others (blue) keep their eigenvalues. The largest remaining eigenvalue of A in absolute value is therefore the dominant eigenvalue of B, and power iteration on B finds it. The layer picture uses every eigenvector, but computing B does not: the lemma below shows that it only needs the eigenvectors already found, through the projection P=Id−QQ⊤ that removes the components along them.
3.4 Lemma 1 [Orthogonal deflation]
Let q1,…,qr be exact orthonormal eigenvectors of A, with associated eigenvalues λ1,…,λr:
Aqi=λiqi,qi⊤qℓ=δiℓ,i,ℓ∈{1,…,r}.
Define
Q=[q1,…,qr],D=diag(λ1,…,λr),P=Id−QQ⊤.
Then
B=PAP=A−QDQ⊤=A−i=1∑rλiqiqi⊤.
Moreover,
The eigenvectors already found have the eigenvalue 0 for B. For each column qi of Q, Bqi=0, for i∈{1,…,r}.
On the vectors orthogonal to the eigenvectors already found, B acts exactly as A. That is, Bx=Ax for every x with Q⊤x=0, which is to say for every x∈range(P)=span{q1,…,qr}⊥.
The other eigenpairs of A are eigenpairs of B. Every eigenvector of A orthogonal to q1,…,qr is an eigenvector of B with the same eigenvalue. Hence the restriction of B to range(P)=span{q1,…,qr}⊥ has the remaining eigenvalues of A, counted with multiplicity, and B has the eigenvalue 0 on range(P)⊥=span{q1,…,qr}.
Remark 8 [Two views of the deflated matrix]
The two formulas for B are two ways of computing a product Bx with one product by A.
Projection,Bx=P(A(Px)): remove from x its components along the stored eigenvectors, apply A, and remove those components again. A vector orthogonal to every stored eigenvector passes through unchanged (item 2), and a stored eigenvector is sent to 0 (item 1).
Subtraction,Bx=Ax−∑i=1rλi(qi⊤x)qi: apply A, then subtract the layers of the stored eigenpairs, exactly as in the picture above.
The two agree because the qi are exact eigenvectors, and the proof below uses AQ=QD. For computed, approximate eigenvectors they differ. The projection still returns vectors exactly orthogonal to the stored ones, so the eigenvectors found later stay orthogonal to them, which is why the implementation below applies P explicitly.
Extend q1,…,qr to an orthonormal eigenbasis q1,…,qd of A. For every i>r,
Q⊤qi=0,Bqi=Aqi=λiqi.
Thus the eigenpairs associated with qr+1,…,qd are preserved on range(P).
■□
When an eigenvalue is repeated, deflating one unit eigenvector removes only that direction. Other directions in the same eigenspace that are orthogonal to it retain their eigenvalue.
Ordinary power iteration can be applied repeatedly to the deflated operator whenever the hypotheses of Proposition 1 hold at each stage. The following modification also handles opposite eigenvalues of equal magnitude, without requiring a spectral bound, such as a lower bound β≤λmin that would make A−βId positive semidefinite and thus free of opposite eigenvalues. It squares the operator instead, which makes it positive semidefinite without any bound, and then separates the signs.
3.5 Proposition 2 [Power iteration with sign separation]
Let B∈Rd×d be symmetric, and let w0 be any unit vector. If Bw0=0, then (0,w0) is an eigenpair. Otherwise, the iteration
wk+1=∥B2wk∥B2wk
is well defined and converges to a unit vector w satisfying
B2w=ρ2w
for some ρ>0. This value can be recovered as
ρ=∥Bw∥.
Define
u+=21(w+ρBw),u−=21(w−ρBw).
Then
Bu+=ρu+,Bu−=−ρu−.
At least one of u+ and u− is nonzero. Normalizing either nonzero vector gives an eigenvector of B.
Remark 9 [Sign separation is a reflection]
Squaring maps +ρ and −ρ to the same eigenvalue ρ2 of B2, so the obstruction of Remark 5 disappears: B2 is positive semidefinite, and its largest eigenvalue present in w0 is strictly dominant. The price is that the limit w may mix the two eigenspaces of B for +ρ and −ρ. One more product separates them. The limit lies in the eigenspace of B2 for ρ2, which is the sum of these two eigenspaces, so we can write w=w++w− with Bw±=±ρw±. Then Bw/ρ=w+−w− is the On the sum of the two eigenspaces, B/ρ fixes every vector of the +ρ eigenspace and reverses every vector of the −ρ eigenspace: it is the mirror reflection across the +ρ eigenspace. Averaging a vector with its mirror image keeps its part on the mirror; averaging it with minus its mirror image keeps the other part. of w that fixes the blue part and reverses the orange part, so u+=w+ and u−=w− are exactly the orthogonal projections of w onto the two eigenspaces. The orthogonality and norm identities in the proof below are this observation written out. Panel (b) of Figure 2 runs this construction on a matrix whose eigenvalues +1 and −1 tie for the largest magnitude: the residual of its candidate tends to zero, by a factor of about 0.8 per product, where plain power iteration stalls.
Proof
Let ΠλB denote the orthogonal projection onto the eigenspace of B associated with λ. Since Bw0=0, the initialization has a nonzero component in at least one nonzero eigenspace of B. Define
ρ=max{∣λ∣:λ∈Λ(B),ΠλBw0=0}>0.
The eigenvalues of B2 are the numbers λ2. Its largest eigenvalue present in w0 is therefore ρ2>0, and every other distinct eigenvalue present is strictly smaller than ρ2.
Proposition 1, applied to B2, shows that the iterates are well defined and converge to a unit vector w in the eigenspace associated with ρ2:
B2w=ρ2w.
Also,
u++u−=w=0,
so the two vectors cannot both be zero. More precisely,
u+⊤u−=41(∥w∥2−ρ2∥Bw∥2)=0,
and
∥u+∥2+∥u−∥2=1.
Therefore,
max{∥u+∥,∥u−∥}≥1/2.
Choosing the vector with larger norm avoids normalizing a potentially very small component.
■□
This construction does not require forming B2. Its action is obtained through B2x=B(Bx), using two applications of B. Squaring removes the distinction between the eigenvalues +ρ and −ρ during power iteration. The subsequent construction of u+ and u− separates their eigenspaces.
3.6 Corollary 2 [Recovery of a complete eigenbasis]
For every real symmetric matrix A, orthogonal deflation combined with Proposition 2 recovers a complete orthonormal eigenbasis, including repeated and zero eigenvalues, when each inner power iteration is taken to its limit. The procedure is as follows.
Initialize Q0∈Rd×0 as an empty matrix, that is, a matrix with no columns, so that Q0Q0⊤=0d×d (each entry is an empty sum).
At stage j∈{1,…,d}, define Pj−1=Id−Qj−1Qj−1⊤ and Bj=Pj−1APj−1. Choose z such that Pj−1z=0, and set w0=Pj−1z/∥Pj−1z∥. For example, a fresh draw z∼N(0,Id) gives a nonzero denominator with probability one.
If Bjw0=0, set qj=w0 and λj=0. Otherwise, apply Proposition 2 to Bj, choose the larger of the two sign-separated components in norm, and normalize it to obtain qj. Recover the eigenvalue as λj=qj⊤Aqj.
Append the new eigenvector, Qj=[Qj−1,qj], and continue to the next stage.
Proof
The proof is by induction on the stage j. Assume that the columns of Qj−1 are orthonormal eigenvectors of A, which holds trivially for j=1, since Q0 has no columns. Let
Sj={x∈Rd:Qj−1⊤x=0}
be the orthogonal complement of the eigenvectors already found, with S1=Rd. We show that qj is a unit eigenvector of A that lies in Sj, so that the columns of Qj=[Qj−1,qj] are again orthonormal eigenvectors of A.
On Sj, Bj acts as A. By Lemma 1, Bj=A−Qj−1Dj−1Qj−1⊤, where Dj−1 is the diagonal matrix of the eigenvalues already found. Hence Bjx=Ax for every x∈Sj.
Bj maps every vector into Sj. Since Qj−1⊤Qj−1=Ij−1,
Qj−1⊤Pj−1=Qj−1⊤−Qj−1⊤Qj−1Qj−1⊤=0,
so Pj−1y∈Sj for every y∈Rd. In particular, Bjx=Pj−1(APj−1x)∈Sj for every x.
Every vector computed at stage j lies in Sj. The initialization w0=Pj−1z/∥Pj−1z∥ does, by the previous step. Each iterate wk+1 is a multiple of Bj(Bjwk), so it does too. The limit w does, because a subspace is closed. Finally, u±=21(w±Bjw/ρ) are linear combinations of the two vectors w and Bjw of Sj.
Conclusion. The matrix Bj is symmetric, so Proposition 2 applies to it. If Bjw0=0, the vector qj is a normalized u+ or u−. It lies in Sj, that is, Qj−1⊤qj=0, and it satisfies Bjqj=μqj with μ=ρ or μ=−ρ. The first step turns this into
Aqj=Bjqj=μqj,
so qj is an eigenvector of A, and λj=qj⊤Aqj=μ. If Bjw0=0, then qj=w0 lies in Sj, and the first step gives
Aw0=Bjw0=0,
so w0 is an eigenvector of A with the eigenvalue λj=0.
Thus each stage adds one unit eigenvector orthogonal to all previously stored eigenvectors.
After d stages,
Qd⊤Qd=Id and
AQd=Qddiag(λ1,…,λd).
Since the square matrix Qd is orthogonal, QdQd⊤=Id, and hence
A=Qddiag(λ1,…,λd)Qd⊤.
■□
Any admissible nonzero initialization suffices to recover some eigenpair at each stage. With fresh Gaussian initializations, the procedure selects an eigenvalue of largest remaining absolute value with probability one. Consequently, the eigenvalues are obtained in nonincreasing absolute value, with arbitrary ordering within ties. Algebraic ordering can be obtained by sorting afterward.
The exact test Bjw0=0 identifies a zero eigenvector. It does not require concluding that the entire remaining operator is zero.
Figure 3 below shows the adjacency matrix of a bipartite graph, whose nonzero eigenvalues come in opposite pairs and whose eigenvalue 0 is repeated. Figure 4 runs this procedure on it and compares it with deflation combined with ordinary power iteration, with and without a shift.
Figure 3. The bipartite graph used in Figures 4 and 11. (a) The graph has 30 nodes on one side and 20 on the other, each side split into two planted communities (teal and gold). Each of the 30×20 possible edges between the two sides is present with probability 0.45 when the communities match and 0.08 otherwise, and no edge joins two nodes of the same side. The draw has 148 edges, 128 within communities (colored) and 20 across (grey, dashed). (b) Its adjacency matrix W, with the nodes ordered side by side and community by community. The diagonal blocks are zero, and the off-diagonal blocks are the 30×20 biadjacency matrix B and its transpose, so a product Wx is one pass over the 148 edges. (c) The 50 eigenvalues of W in decreasing order. They are symmetric about 0. Writing a vector of R50 as (x,y) with x∈R30 and y∈R20, W(x,y)=(By,B⊤x), so a pair of singular vectors with Bv=σu and B⊤u=σv gives the eigenvector (u,v) for σ and (u,−v) for −σ. The 20 singular values of B give the 40 nonzero eigenvalues, and 0 is repeated 30−20=10 times.Figure 4. A complete eigenbasis of the adjacency matrix of the bipartite graph of Figure 3, applied through its edge list (Section 2.2). Because the graph is bipartite, every nonzero eigenvalue σ comes with −σ, so plain power iteration meets the tie of Remark 5 at once, and 0 is an eigenvalue of multiplicity 10. Every stage starts from a fresh Gaussian vector and, as in Corollary 2, is taken to its limit: it runs until the residual of its candidate for the deflated matrix Bj is below 10−13. (a) Eigenvalue found at each stage. Algorithm 1 (blue) finds the eigenvalues in nonincreasing absolute value, the two members of each pair at consecutive stages, and the ten zero eigenvalues last. Power iteration on A−λmin(A)Id with deflation (purple; Corollary 1 and Lemma 1) finds them in decreasing order. Grey lines mark the 41 distinct eigenvalues. (b) Number of eigenpairs found against the number of products. Algorithm 1 uses 15,296 products in all, two for each zero eigenvalue. The shifted iteration uses 61,453, even with the most favorable shift λmin(A)=−6.48: adding 6.48 to every eigenvalue brings the ratios that set the speed of convergence (Theorem 1) closer to one, and each zero eigenvalue takes a full stage. Plain deflated power iteration (orange) never completes its first stage, and after 20,000 products its residual is still above 1.9. (c) Accuracy of the eigenbasis of Algorithm 1, far below the tolerance 10−10∥A∥ (grey, dashed): residuals (blue) below 3.3⋅10−13, eigenvalues within 6.3⋅10−15 of the spectrum (purple; exact zeros drawn at 10−17), inner products between distinct eigenvectors (grey) below 7.3⋅10−14, and ∥A−Q50diag(λ1,…,λ50)Q50⊤∥=2.2⋅10−13.
3.7 Lemma 2 [Residual certificate for an approximate eigenpair]
For any unit vector q, define
λ=q⊤Aq
and
r=Aq−λq.
Then minλ∈Λ(A)∣λ−λ∣≤∥r∥.
Proof
Expand q in an orthonormal eigenbasis u1,…,ud of A:
The squared residual is thus an average of the squared distances from λ to the eigenvalues, with weightsci2 that sum to one, so it is at least the smallest of them. Taking square roots gives the result.
■□
Remark 10 [A computable bound on an unknown distance]
The residual r=Aq−λq is computable with one product, and it bounds an unknown quantity: the distance from the estimate λ=q⊤Aq to the nearest true eigenvalue. If q is exactly an eigenvector, the residual is zero, and the bound says that λ is exactly an eigenvalue. The lemma has three consequences.
A stopping rule with a guarantee. The eigenvalues of A are unknown, so the error of λ cannot be measured directly. The lemma converts the computable residual into a guarantee: the interval [λ−∥r∥,λ+∥r∥] contains an eigenvalue of A. Stopping when ∥r∥≤ε therefore ensures that λ is within ε of an eigenvalue, and this is how Algorithm 1 and Algorithm 2 stop. Lemma 7 is the counterpart for nonsymmetric matrices.
Independence from the method. The lemma involves only A, q and one product, not the way q was obtained. It applies equally to an unconverged power iterate, an approximate eigenvector produced by the Lanczos method of Section 4, or a vector from a late deflation stage whose earlier eigenvectors carry rounding errors. This is why the algorithms compute the residual with the original matrix A rather than with the matrix iterated at stage j of the deflation in Corollary 2, Bj=Pj−1APj−1, where Pj−1=Id−Qj−1Qj−1⊤ removes the directions of the eigenvectors q1,…,qj−1 already found. Errors in these eigenvectors change Bj, but not the certificate.
What it does not say. It certifies some eigenvalue, not a particular one, such as the largest remaining one. Nor does a small residual alone make q close to an eigenvector. The expansion in the proof above also bounds the angle between q and an eigenspace. Fix an eigenvalue λ, and let δ=min{∣μ−λ∣:μ∈Λ(A)∖{λ}} be the distance from λ to the other eigenvalues, assumed positive. Then
The first line is the expansion of the proof. The second drops the terms with λi=λ, which are nonnegative. The third uses ∣λi−λ∣≥δ for every λi=λ, which is the definition of δ. For the fourth, the ui with λi=λ form an orthonormal basis of Eλ, so Πλq=∑i:λi=λciui and q−Πλq=∑i:λi=λciui, whose squared norm is ∑i:λi=λci2 by Pythagoras. The last is the sine of the angle between the unit vector q and Eλ, as in Theorem 1: the side of the right triangle opposite the angle, over the hypotenuse ∥q∥=1. Dividing by δ2 and taking square roots gives
sin∠(q,Eλ)≤δ∥r∥=min{∣μ−λ∣:μ∈Λ(A)∖{λ}}∥Aq−λq∥.
The eigenvector is accurate only when λ=q⊤Aq is well separated from the rest of the spectrum. When another eigenvalue is close, a small residual still allows q to mix the two eigenvectors.
For example, take A=diag(1,3) and the unit vector q=[cos(t),sin(t)]⊤, at the angle t from the eigenvector [1,0]⊤ of 1. Then
For sin(t)=0.1, λ=1.02 and ∥r∥≈0.199, so the lemma guarantees an eigenvalue in [0.821,1.219], and 1 is there. The eigenvector bound gives sin(t)≤0.199/1.98≈0.1005, close to the true 0.1. The actual eigenvalue error, 0.02, is ten times smaller than the residual: it is quadratic in the angle t, while the residual is linear in it.
This holds in general. Since the ci2 sum to one, λ=∑ici2λi, so λ−λ=∑ici2(λi−λ), where the terms with λi=λ vanish. Hence
where the equality combines the last two lines of the chain above and the last step is the eigenvector bound. In the example, the right-hand side is 2⋅0.1992/1.982≈0.0202, against the true error 0.02. Once the residual is small compared with the gap δ, the eigenvalue error is of the order of ∥r∥2, far below the certificate ∥r∥. The certificate is therefore valid but pessimistic: it guarantees roughly the square root of the accuracy actually reached, as Figure 5 shows.
Remark 11 [Any scalar, and the best one]
The proof never uses the definition of λ. Replacing λ by any scalar μ, the same expansion gives ∥Aq−μq∥2=∑i=1dci2(λi−μ)2≥mini(λi−μ)2, that is,
λ∈Λ(A)min∣μ−λ∣≤∥Aq−μq∥for every unit vector q and every scalar μ.
No relation between q and μ is required. They can be chosen independently, and a poor match only makes the bound loose, never wrong. For A=diag(1,3), q=[1,0]⊤ and μ=10, it reads 7≤9. Hence the certificate applies to an eigenvalue estimate computed by any means, as long as it is paired with a unit vector. The Rayleigh quotient is the best choice of μ. Indeed, expanding the square with ∥q∥=1 and λ=q⊤Aq gives
where the last line is the second one at μ=λ. The blue term does not depend on μ, so the certificate is smallest at μ=λ, where it equals the residual of Lemma 2, and the orange term is the price of any other estimate. As a function of μ, the certificate (μ−λ)2+∥r∥2 is a hyperbola centered at the Rayleigh quotient λ, and it lies above the distance minλ∈Λ(A)∣μ−λ∣ from μ to the spectrum, as panel (a) of Figure 5 below shows.
This distance is the pointwise minimum of the functions μ↦∣μ−λ∣, one for each eigenvalue λ, whose graphs are V shapes with their vertex on the axis at λ. Its graph is therefore the lower envelope of these V shapes: a zigzag, piecewise linear with slopes ±1, that vanishes at each eigenvalue. Between two consecutive eigenvalues λ<λ′, the nearest eigenvalue is λ up to the midpoint (λ+λ′)/2 and λ′ after it, so the distance equals μ−λ, which rises with slope 1, and then λ′−μ, which falls with slope −1. Each gap of the spectrum thus carries a triangle whose height (λ′−λ)/2 is half the gap, and beyond the extreme eigenvalues the distance grows along half-lines of slope ±1. The hyperbola has the same slopes far from its center, since it approaches ∣μ−λ∣ as ∣μ−λ∣ grows, and it is lowest at μ=λ, where its height is ∥r∥. The bound says that the hyperbola never goes below this envelope. At its lowest point, this is Lemma 2: the height ∥r∥ is at least the distance from λ to the spectrum, so a small residual places λ near a zero of the envelope, that is, near an eigenvalue.
The remark has three consequences.
Certifying an estimate obtained elsewhere. An eigenvalue estimate μ that does not come from a Rayleigh quotient, such as a value known from theory, the value α−μmax recovered from a shifted run in Corollary 1, or the value ±ρ of Proposition 2, is certified by pairing it with a unit vector q and computing one product: some eigenvalue of A lies within ∥Aq−μq∥=(μ−λ)2+∥r∥2 of μ, r=Aq−λq.
Reporting the Rayleigh quotient. For a given q, replacing μ by λ never weakens the certificate, and strictly improves it when μ=λ. It costs no extra product, since λ=q⊤(Aq) reuses the product already needed for the residual. This is why Corollary 2 recovers each eigenvalue as λj=qj⊤Aqj rather than reporting ±ρ: the two coincide in exact arithmetic, but in practice qj is only approximately an eigenvector. Geometrically, q⊤r=λ−λ=0, so λq is the orthogonal projection of Aq onto the line spanned by q, the point of that line closest to Aq. The certificate measures how well the pair (μ,q) fits together, not the error of μ alone. In the example of Remark 10, the exact eigenvalue μ=1 receives the certificate ∥Aq−q∥=∥[0,0.2]⊤∥=0.2, slightly worse than the 0.199 of λ=1.02, because q is not its eigenvector.
Testing a given value. For a fixed μ, every unit vector q gives an upper bound on the distance from μ to the spectrum, so a single well-chosen q suffices to show that μ is nearly an eigenvalue. With μ=0, the bound reads minλ∈Λ(A)∣λ∣≤∥Aq∥.
Any unit vector q that A nearly annihilates proves that A has an eigenvalue near zero, that is, that A is nearly singular. Combined with ∥A∥=max∥z∥=1∥Az∥≥∥Aq′∥ for any other unit vector q′, it gives a lower bound on the condition number computed from two products,
κ(A)=maxλ∈Λ(A)∣λ∣/minλ∈Λ(A)∣λ∣=∥A∥/minλ∈Λ(A)∣λ∣≥∥Aq′∥/∥Aq∥. The inequality bounds the two parts of the fraction separately. The numerator can only decrease, from ∥A∥ to ∥Aq′∥. The denominator can only increase, from minλ∣λ∣ to ∥Aq∥, by the bound with μ=0. Both changes make the fraction smaller. A vector q′ that A stretches and a vector q that it nearly annihilates therefore give a large ratio, which proves that A is ill conditioned.
Figure 5. The residual certificate on the sample covariance matrix A=n1X⊤X of n=400 centered observations of d=60 variables, applied as X⊤(Xv)/n (Section 2.2). Its three largest eigenvalues, the variances along the first three principal components, are 9.72, 4.92 and 3.20, and the other 57 lie in [0.40,1.80]. (a) For the iterates vk of power iteration from one Gaussian start, with Rayleigh quotients λk=vk⊤Avk and residuals rk=Avk−λkvk, the certificate μ↦∥Avk−μvk∥=(μ−λk)2+∥rk∥2 of Remark 11 (light to dark as k grows) lies above the distance minλ∈Λ(A)∣μ−λ∣ from μ to the spectrum (black, with the eigenvalues marked in orange on the axis) for every scalar μ. This distance is the lower envelope of the V shapes ∣μ−λ∣ (Remark 11), a zigzag with one triangle over each gap of the spectrum, of height half the gap. The tallest, of height 2.40, spans the gap from 4.92 to 9.72, and those over the 57 closely spaced eigenvalues in [0.40,1.80] are too small to see. The minimum of each hyperbola, marked by a dot, is attained at μ=λk, whose value is written under the axis at the foot of the dashed line, and equals the residual norm ∥rk∥ of Lemma 2. (b) The same quantities along the whole run rather than for a few iterates. Each curve is the median over 200 Gaussian starts of a quantity computed from vk, and its band contains the middle half of the runs. The residual ∥rk∥ (blue) shrinks by the factor λ2/λ1≈0.51 per product, and the distance from λk to the nearest eigenvalue (orange) by its square (Theorem 1), until rounding errors of order 10−15 stop it. The certificate is valid but pessimistic: it bounds the error by the residual, whereas the error of a Rayleigh quotient is of the order of the squared residual (Remark 10). After 20 products, the median residual is 1.2⋅10−5 and the median distance 3.2⋅10−11, about its square divided by the gap 4.81 between the two largest eigenvalues. In the first steps, λk lies among the small eigenvalues, and its distance to the nearest of them (orange) is much smaller than its distance to λmax(A) (grey, dotted): the certificate concerns some eigenvalue, not necessarily the largest. (c) The pairs formed by the residual and the distance, one point for each iterate vk, 0≤k≤39, of 50 runs, colored by k from light to dark. None lies above the diagonal, as Lemma 2 requires, and over all 200 runs and 40 steps the largest ratio of the distance to the residual is 0.995. Distances below 10−15, at rounding level, are drawn at 10−15.
3.8 Algorithm 1 [Deflated squared power iteration with sign separation]
Corollary 2 recovers a complete eigenbasis, but each of its stages runs the iteration of Proposition 2 to its limit, which takes infinitely many products. The algorithm below is its finite version. It keeps the two loops of Corollary 2, an outer loop over deflation stages (Lemma 1) and an inner loop of squared power iteration on the deflated matrix B=PAP of the current stage (Proposition 2), with three changes.
No matrix is formed. Neither A, P nor B is built. A product with B is a product with A followed by the projection Px=x−Q(Q⊤x), as in the projection view of Remark 8.
The signs are separated at every step. Each inner step turns the current iterate into a candidate eigenvector by the construction of u± in Proposition 2, without waiting for the limit.
The inner loop stops on a certificate. The candidate is accepted as soon as its residual with A is below a tolerance, which by Lemma 2 places its Rayleigh quotient near an eigenvalue of A.
At a given deflation stage, let Q contain the previously computed orthonormal vectors. At the first stage, no vector has been computed yet: Q is the empty matrix with d rows and no columns, as Q0 in Corollary 2, so that QQ⊤=0d×d and P=Id, and the steps below are plain squared power iteration with sign separation. Each stage ends by appending the vector it accepts to Q (step 4), and the next stage starts with the enlarged Q, until the desired number of eigenpairs has been found. Figure 6 draws these two loops.
Choose tolerances εabs>0 and εrel≥0. Step 4 accepts a candidate when ∥r∥≤εabs+εrel∥Aq∥, and by Lemma 2 the accepted λ then lies within this threshold of an eigenvalue of A.
In the simplest case εrel=0, every eigenvalue is found to within the same distance εabs. This is the recommended choice, used in Figures 4 and 11 and in the constant-threshold runs of Figure 7 with εabs=10−10∥A∥. This tolerance must be positive, since rounding errors keep the computed residual of the order of the unit roundoff times ∥A∥, never exactly zero.
A positive εrel asks instead for a number of correct digits. The single tolerance εabs of the first case treats large and small eigenvalues unequally, since it accepts the same error on all of them. With εabs=10−10, an error of 10−10 on the eigenvalue λ=10 gives the estimate 10.0000000001, which agrees with λ in its first eleven significant digits, but the same error on λ=0.001 gives 0.0010000001, which agrees in only seven. The relative part removes this imbalance by accepting an error proportional to the eigenvalue. With a negligible εabs, the test accepts a residual up to εrel∥Aq∥, and ∥Aq∥ is nearly ∣λ∣: since Aq=λq+r with q⊤r=0 and ∥q∥=1, we have ∥Aq∥2=λ2+∥r∥2, and so ∣λ∣≤∥Aq∥≤∣λ∣+∥r∥. With εrel=10−10, the accepted residual is thus about 10−10∣λ∣, that is, about 10−9 near λ=10 and about 10−13 near λ=0.001. By Lemma 2, the estimates are then within these distances of λ, for instance 10.000000001 and 0.0010000000001, and both agree with λ in about ten significant digits. In general, the upper bound ∣λ∣+∥r∥ on ∥Aq∥ turns the test ∥r∥≤εrel∥Aq∥ into ∥r∥≤εrel∣λ∣/(1−εrel), and Lemma 2 bounds the relative error by ∣λ−λ∣/∣λ∣≤∥r∥/∣λ∣≤εrel/(1−εrel)≈εrel, whatever the eigenvalue. Using ∥Aq∥ rather than ∣λ∣ in the test thus changes little, and it gives the test a geometric meaning: r is the part of Aq orthogonal to q, so ∥r∥/∥Aq∥ is the sine of the angle between q and Aq, and the test asks that A map q nearly onto its own line. The absolute part remains as a floor. For the eigenvalue λ=0, the relative part accepts no residual at all, rounding errors never bring the residual to zero, and without εabs the stage would never stop. In deflation, this balance has a price. The eigenvalues are found in decreasing absolute value, so the accepted residual shrinks from stage to stage, from about 10−9 to about 10−13 in the example, and later stages can stall, for a reason explained below.
At each stage, start from a unit vector w satisfying Q⊤w=0, for instance w=Pz/∥Pz∥ for a fresh draw z∼N(0,Id) (at the first stage, any unit vector qualifies), and repeat the following steps.
First product. Compute h1=Aw, a=Ph1 and ρ=∥a∥. If ρ=0, go to step 5.
Second product. Set y=a/ρ, then compute h2=Ay and b=Ph2.
Candidate. Choose the scalar sign s=+1 if w⊤y≥0 and s=−1 if w⊤y<0. This selects whichever of w+y and w−y has larger norm. Form n=∥w+sy∥ and q=(w+sy)/n. The product of the candidate with the original matrix is available by linearity: Aq=(h1+sh2)/n.
Check or update. Compute λ=q⊤Aq and r=Aq−λq. If ∥r∥≤εabs+εrel∥Aq∥, accept (λ,q), append q to Q, and end the stage. Otherwise, update w←b/∥b∥ and return to step 1.
Zero projected product. If ρ=0, use q=w as the candidate, compute λ=q⊤Aq and r=Aq−λq, where Aq=h1, and apply the test of step 4. Either way, the stage ends here and does not return to step 1. The vector w satisfies Q⊤w=0: the starting vector is chosen so, and each update b/∥b∥ is a projected vector, with Q⊤P=Q⊤−(Q⊤Q)Q⊤=0 since Q⊤Q=I. Hence Pw=w−Q(Q⊤w)=w, and the condition ρ=0 means Bw=PAPw=0, so there is no direction y to iterate on. If the test passes, accept (λ,q) and append q to Q, as in step 4. In exact arithmetic it always passes, with λ=0, by the zero case of Corollary 2. If it fails, which can only come from errors in the stored vectors or from rounding, report nonconvergence or restart the stage from a new w.
■□Figure 6. The execution of Algorithm 1. Each stage (outer frame) computes one eigenpair and appends its vector to Q, and the green arrow starts the next stage, on the matrix deflated further. Within a stage, the inner loop (inner frame) is squared power iteration on B=PAP. Each inner step costs the two products of steps 1 and 2, tinted with the colours of h1 and h2, and the blue arrow returns to step 1 when the test fails. Step 5 handles the case ρ=0, where Bw=0. If its test fails, which can only come from rounding, the stage is restarted from a new w or nonconvergence is reported.
How the steps implement the theory
Let B=PAP be the deflated matrix of the current stage. The vectors w and y=Ph1/∥Ph1∥ are orthogonal to the stored vectors (the columns of Q): w by the choice of the start and because each update is a projected vector, and y because Q⊤P=0. Hence Pw=w−QQ⊤w=w and Py=y−QQ⊤y=y, and the quantities computed in steps 1 and 2 are products with B:
The update is the iteration of Proposition 2. Normalizing b gives w←B2w/∥B2w∥, the iteration of Proposition 2 on B, at the cost of two products with A per step. It is well defined: ρ>0 implies b=0, because
w⊤B2w=∥Bw∥2=ρ2>0.
By Proposition 2, the iterates converge to a unit vector w with B2w=ρ2w, where ρ is the largest absolute value among the eigenvalues of B present in the starting vector. Squaring turns the eigenvalues +ρ and −ρ of B into the same eigenvalue ρ2 of B2, so the tie of Remark 5 cannot stall the iteration.
The candidate is the sign separation of Proposition 2, before the limit. Since w+sy=w+sBw/ρ=2us, the candidate q=us/∥us∥ is the normalized u+ when s=+1 and the normalized u− when s=−1, computed from the current iterate instead of the limit. The sign picks the longer of the two, since ∥u+∥2−∥u−∥2=w⊤Bw/ρ=w⊤y, and the normalization is never small:
n2=∥w+sy∥2=∥w∥2+∥y∥2+2sw⊤y=2+2∣w⊤y∣≥2,
which is the bound max{∥u+∥,∥u−∥}≥1/2 of the proof of Proposition 2. Forming the candidate before the limit is justified by a direct calculation:
(B−sρId)q=nBw+sB2w/ρ−sρw−Bw=ρns(B2w−ρ2w).
The blue terms cancel because s2=1. Since n≥2, the residual of (sρ,q) as an eigenpair of B is at most ∥B2w−ρ2w∥/(2ρ), which measures how far the current iterate is from an eigenvector of B2. As squared power iteration converges, the candidate therefore becomes an eigenvector of B, for the eigenvalue ρ or −ρ, and nothing is lost by separating the signs at every step.
The check uses A, not B. With exact stored eigenvectors, an eigenvector of B orthogonal to them is an eigenvector of A with the same eigenvalue, by the first step of the proof of Corollary 2, and its Rayleigh quotient λ=q⊤Aq is then sρ. Computed eigenvectors are only approximate, so B is not exactly the deflated matrix of Lemma 1, and an eigenvector of B need not be one of A. The test therefore computes the residual with A itself, which Lemma 2 turns into a guarantee however q was obtained (Remark 10), and it reports the Rayleigh quotient, the scalar with the smallest certificate (Remark 11). For the same reason, the exact recovery of Corollary 2 does not carry over to finite computations: convergence for B need not drive the residual for A below an arbitrary tolerance, as the end of this section shows.
The zero case. If ρ=0, then Bw=0, and w is an eigenvector of B for the eigenvalue 0, the first case of Proposition 2. In exact arithmetic, this can only happen at the first step of a stage: if Bw=0, then w has a component in an eigenspace of B with a nonzero eigenvalue, and so does every later iterate, since it is a multiple of B2kw. A small computed ρ does not establish that w is an eigenvector for 0, still less that every remaining eigenvalue is zero, which is why the zero case is also accepted only through the residual test of step 4.
Floating-point safeguards. The implementation should reorthogonalize projected vectors, impose a maximum number of inner iterations, and report nonconvergence if the original residual does not reach the requested tolerance. No vector with computed norm zero should be normalized. A numerical breakdown should instead trigger a residual check, a restart, or a report of nonconvergence.
Why a stage can stall. The residual for A splits exactly into a part that the inner iteration reduces and a part that it inherits. Let λ1,…,λr be the eigenvalue estimates of the stored pairs, and let the columns of R=AQ−Qdiag(λ1,…,λr) be their residuals. For a unit candidate q with Q⊤q=0 and λ=q⊤Aq, we have Pq=q, hence Aq=Bq+QQ⊤Aq, where Bq−λq lies in range(P) and QQ⊤Aq in range(Q). Since Q⊤q=0, we also have Q⊤Aq=(AQ)⊤q=R⊤q. The two parts are orthogonal, so
∥Aq−λq∥2=∥Bq−λq∥2+∥R⊤q∥2.
The inner iteration drives the orange term to zero. The purple term is inherited from the stored pairs: each entry of R⊤q is the inner product of q with a stored residual, so the term is at most the sum of their squared norms, but nothing in the inner iteration makes it small. A stage therefore stalls when its threshold εabs+εrel∥Aq∥ is below the inherited part. This happens readily with εrel>0, because the threshold shrinks with ∣λ∣, whereas the inherited part comes from earlier stages, whose eigenvalues are larger in absolute value and whose accepted residuals can be nearly as large as their own thresholds. A threshold that is the same at every stage (εrel=0) removes this systematic mismatch, but over many stages the inherited parts can still exceed it. Continuing each stage until the residual for B reaches rounding level, and only then checking the residual for A, keeps every inherited part far below the threshold, at the price of more products. Figure 7 below shows the stall on the Hessian of a neural network. The complete eigenbasis in Figure 4 was computed in the last way. With the same constant threshold but without the continuation, that run stalls at stage 20 of 50, where the inherited part is 6.89⋅10−10 and the threshold 6.48⋅10−10.
Figure 7.Algorithm 1 on the Hessian A=∇2f(θ) of the training loss f of a small neural network, applied by automatic differentiation (Section 2.1). The network is a fully connected tanh network with layer widths 4, 16, 16 and 1, so that θ∈R369, trained by 8,000 steps of gradient descent on n=256 Gaussian inputs labelled by a random teacher network. The top row describes the matrix, and the bottom row runs the algorithm on it. The 8 eigenvalues of largest absolute value at the last step, which the algorithm computes, are colored from light to dark in panels (b), (c) and (d). (a) The squared losses along gradient descent, on the training inputs (blue) and on 1,024 fresh inputs labelled by the same teacher (orange, dashed). The Hessian A of panels (c) to (f) is taken at the last step (black dot), where the training loss is 4.8⋅10−4 and the test loss reaches its lowest value, 1.2⋅10−3. (b) The eigenvalues of the Hessian along the same run, computed from the dense 369×369 matrix every 80 steps, on a vertical scale that is logarithmic beyond ±0.01 and linear in between. The 8 largest are drawn as lines, and the other 361 lie in the grey band, whose lower edge is the smallest eigenvalue. Training sharpens the top of the spectrum: the five largest eigenvalues rise from between 1.07 and 2.32 at initialization to between 3.04 and 6.50, and they separate from the next three, which settle between 0.25 and 0.34. The negative curvature nearly disappears: the smallest eigenvalue rises from −0.83 to −0.0014, and the number of negative eigenvalues falls from 191 to 139. (c) The 369 eigenvalues at the last step, counted by absolute value in bins of half a decade, the positive ones above the axis (blue) and the negative ones below it (orange). Their absolute values span more than eight orders of magnitude, from 1.1⋅10−8 to 6.50, and 269 of them lie between 10−5 and 10−2. The five largest stand far above the rest, and every negative eigenvalue lies above −0.0015. The ticks at the top mark the 8 largest, which the bottom row computes. The bottom row uses two threshold rules (Section 3.8): a constant one, εabs=10−10∥A∥ and εrel=0, where ∥A∥=6.50 is the absolute value of the first eigenvalue found, and a relative one, εabs=εrel=10−10. (d) One run with the constant threshold: the residual of the candidate against the number of Hessian–vector products, counted over all stages, with each stage labelled by the eigenvalue it finds. Each stage decreases geometrically, by a factor per product equal to the ratio of the next largest remaining absolute value to its own, until it reaches the threshold (grey, dashed), and the eight stages take 1,072 products. (e) A run with the relative threshold, at its second stage. By the splitting of Section 3.8, the squared residual for A (blue), which the algorithm checks, is the sum of the squares of the residual for the deflated matrix B (orange, dashed) and of the inherited part ∥Q⊤Aq∥ (purple, dotted). The inner iteration drives the first down to 5⋅10−16, but the second stays at 7.31⋅10−10, the residual accepted at the first stage. The threshold of the second stage, 10−10(1+5.90)=6.90⋅10−10, is just below it, so the stage never stops. Stage 2 is where stalls are most common, and the run shown is the first of the 20 to stall there. (f) Fraction of 20 runs, each with its own Gaussian starting vectors, that certify the j largest pairs. With the relative threshold (orange), 14 runs stall: 8 at stage 2, 3 at stage 4, and one each at stages 5, 6 and 7. With the constant threshold (blue), all 20 runs certify the 8 pairs, and in the run of panel (d) every eigenvalue is within 5⋅10−15 of the spectrum.
The operator must stay fixed for the whole run. For a Hessian A=∇θ2f(θ), each product is computed as in Section 2.1, with θ and everything that defines f held fixed.
Each nondegenerate inner iteration above uses two matrix–vector products, h1 and h2, since Aq comes for free by linearity. After j−1 eigenvectors have been stored, projecting a vector costs O(d(j−1)) arithmetic operations. Retaining m eigenvectors requires O(dm) storage. Thus computing and storing a complete eigenbasis still requires O(d2) storage, even though A itself is never formed.
4. The Lanczos method
The Lanczos method constructs an orthonormal basis of a Krylov subspace and projects the original eigenvalue problem onto that subspace. This produces a smaller symmetric tridiagonal matrix whose eigenpairs provide approximations to those of the original matrix. The construction requires only matrix–vector products (Lanczos, 1950). Unless stated otherwise, all identities and termination results below concern exact arithmetic and a fixed matrix A.
4.1 Definition 1 [Krylov subspace]
Given a unit vector q1, the Krylov subspace of order m is
Km(A,q1)=span{q1,Aq1,…,Am−1q1}.
Equivalently,
Km(A,q1)={p(A)q1:p is a polynomial of degree at most m−1}.
These subspaces are nested:
Km(A,q1)⊆Km+1(A,q1).■□
Power iteration started at q1 keeps only the last vector Am−1q1 of this list. The Krylov subspace keeps all of them, so it contains every polynomial filter of degree at most m−1 applied to q1. By the spectral decomposition, p(A)q1=∑λp(λ)Πλq1, so a polynomial that is large at one eigenvalue and small at the others can isolate that eigencomponent much better than the monomial tm−1 used by power iteration.
4.2 The Gram–Schmidt process
The Gram–Schmidt process turns linearly independent vectors v1,…,vm of Rd into orthonormal vectors q1,…,qm, that is, vectors with
qi⊤qj=δij
that span the same subspaces: span{q1,…,qj}=span{v1,…,vj} for every j. It builds them one at a time. The first is q1=v1/∥v1∥. At step j≥2, it subtracts from vj its components along the vectors already built, and normalizes the remainder:
rj=vj−i=1∑j−1(qi⊤vj)qi,qj=∥rj∥rj.
Both properties are proved together, by induction on j. Suppose that q1,…,qj−1 are orthonormal and span the same subspace as v1,…,vj−1. This holds for j=2, since q1 is a unit multiple of v1. Three facts carry it over to j.
First, the remainder rj is not zero, so the division by ∥rj∥ is allowed. If rj were zero, its definition would give vj=∑i<j(qi⊤vj)qi, a vector of span{q1,…,qj−1}. By the hypothesis, this subspace is span{v1,…,vj−1}, so vj would be a combination of v1,…,vj−1, which linear independence forbids. Without independence, rj=0 is exactly the signal that vj adds no new direction.
Second, qj is orthogonal to q1,…,qj−1. For k<j, every term of the sum but i=k vanishes in the inner product with qk since qk⊤qi=δki, that is
qk⊤rj=qk⊤vj−∑i=1j−1(qi⊤vj)qk⊤qi=qk⊤vj−qk⊤vj=0.
Dividing by ∥rj∥ keeps this orthogonality and gives qj=rj/∥rj∥ norm one, so q1,…,qj are orthonormal.
Third, the spans agree. The vector qj is a combination of vj and q1,…,qj−1, hence of v1,…,vj by the hypothesis. Conversely, vj=∥rj∥qj+∑i<j(qi⊤vj)qi is a combination of q1,…,qj.
For example, v1=[1,1,0]⊤ and v2=[1,0,1]⊤ give q1=[1,1,0]⊤/2, then q1⊤v2=1/2, r2=v2−q1/2=[21,−21,1]⊤ and q2=[1,−1,2]⊤/6, and indeed q1⊤q2=(1−1+0)/12=0. Step j costs j−1 inner products and subtractions, so the cost per step grows with j. The following section applies the process to the spanning vectors q1,Aq1,A2q1,… of the Krylov subspace.
4.3 How the Lanczos method works
The idea The Lanczos method runs the Gram–Schmidt process on the Krylov subspace Km(A,q1)=span{q1,Aq1,…,Am−1q1}, one matrix–vector product per step.
The numbers that Gram–Schmidt computes along the way form a small symmetric matrix, and the eigenvalues of that matrix approximate eigenvalues of A.
This section unpacks these sentences. It first explains why the search for eigenvalues needs an orthonormal basis of Km(A,q1), then runs Gram–Schmidt on the Krylov subspace by hand for four steps, states the general step, and finally shows how eigenvalues and eigenvectors come out of the coefficients. The proofs are in Lemmas 3 to 5 and Proposition 3.
After m−1 products, every vector of Km(A,q1) is within reach. It is a combination c0q1+c1Aq1+⋯+cm−1Am−1q1=p(A)q1, and by the spectral decomposition, p(A)q1=∑λp(λ)Πλq1 reweights each eigencomponent of q1 by the value p(λ). Power iteration only ever uses the polynomial t↦tm−1. The Lanczos method considers all of them at once.
To turn these vectors into eigenvalue estimates, recall that the Rayleigh quotient of a unit vector never exceeds λmax(A)=max∥x∥=1x⊤Ax and never falls below λmin(A)=min∥x∥=1x⊤Ax (Remark 7), with equality at the corresponding eigenvectors.
These optimizations run over all unit vectors of Rd, which are out of reach. Running them over the unit vectors of Km(A,q1) only gives two numbers ϑmax(m)≤λmax(A) and ϑmin(m)≥λmin(A). They are the best estimates of the extreme eigenvalues that Rayleigh quotients of vectors of Km(A,q1) can give, and they improve as m grows, because the subspaces are nested. The power iterate v=Am−1q1/∥Am−1q1∥ is a unit vector of Km(A,q1), so the estimate of power iteration, v⊤Av, lies between ϑmin(m) and ϑmax(m): whichever extreme eigenvalue power iteration approaches, the corresponding estimate of the Lanczos method is at least as close. Lemma 5 states these facts precisely.
Searching a subspace needs coordinates, and the convenient ones come from an orthonormal basis q1,…,qm of Km(A,q1), the columns of Qm=[q1,…,qm]. Every vector of the subspace is then x=Qmy=y1q1+⋯+ymqm for a unique y∈Rm, lengths are preserved, ∥x∥=∥y∥, and
x⊤Ax=y⊤Tmy,Tm=Qm⊤AQm∈Rm×m.
The search over Km(A,q1) thus becomes the eigenvalue problem of the m×m matrix Tm, whose entry (Tm)ij=qi⊤Aqj is the coordinate of Aqj along qi. Gram–Schmidt applied to the spanning vectors q1,Aq1,A2q1,… of the Krylov subspace produces the orthonormal basis, and, as the following steps show, it computes these coordinates on the way.
Step 0: Gram–Schmidt on K1(A,q1)=span{q1}={c0q1:c0∈R}. There is nothing to do. The starting vector q1 is a unit vector, so on its own it is already an orthonormal basis of K1(A,q1): it is the first basis vector, obtained without any product.
Step 1: Gram–Schmidt on K2(A,q1)=span{q1,Aq1}={(c0Id+c1A)q1:c0,c1∈R}. The first product gives Aq1, the vector that K2(A,q1) adds. Gram–Schmidt removes its component along q1, whose coefficient α1=q1⊤Aq1 is the Rayleigh quotient of q1. It is the first eigenvalue estimate, the same as that of power iteration. The remainder r1=Aq1−α1q1 is orthogonal to q1, and dividing it by its length β1=∥r1∥ gives the second basis vector q2=r1/β1. Rearranged, the step reads
So q1,q2 is an orthonormal basis of K2(A,q1), and q2=p1(A)q1 for the polynomial p1(t)=(t−α1)/β1 of degree one. In the colors used throughout this section, a purple coefficient is computed at the current step, and a green one is the length of the new direction.
Step 2: Gram–Schmidt on K3(A,q1)=span{q1,Aq1,A2q1}. The vector that K3(A,q1) adds is A2q1, and textbook Gram–Schmidt would orthogonalize it next. The Lanczos method applies A to the newest basis vector q2 instead, which also costs one product. Indeed, by Step 1,
Aq2=β11A2q1−β1α1Aq1.
The orange vector already lies in K2(A,q1)=span{q1,q2}, so removing the components along q1 and q2 deletes it, and only the blue power contributes. Hence Aq2 brings the same new direction as A2q1. Gram–Schmidt needs two coefficients. The one on q1 comes for free from symmetry and Step 1:
q1⊤Aq2=(Aq1)⊤q2=(α1q1+β1q2)⊤q2=β1.
It is the length found in Step 1, now reused, and orange marks such known coefficients. The one on q2 is the only new number, α2=q2⊤Aq2. Subtracting both components and normalizing gives
that is, Aq2=β1q1+α2q2+β2q3. Now q1,q2,q3 is an orthonormal basis of K3(A,q1), and q3=p2(A)q1 for the polynomial p2(t)=((t−α2)p1(t)−β1)/β2 of degree two. The coefficients found so far already give the matrix of the Rayleigh quotient on K2(A,q1),
T2=Q2⊤AQ2=[α1β1β1α2],
whose two eigenvalues are the smallest and largest values of x⊤Ax over the unit vectors of K2(A,q1).
Step 3: Gram–Schmidt on K4(A,q1)=span{q1,Aq1,A2q1,A3q1}, where a coefficient vanishes. The vector that K4(A,q1) adds is A3q1, and, as in Step 2, the third product is applied to the newest basis vector instead, giving Aq3. The argument of Step 2 carries over. Since q3=p2(A)q1, where p2 has degree two and leading coefficient 1/(β1β2),
where the orange vector collects the terms of lower degree. Removing the components along q1,q2,q3 deletes it, so Aq3 brings the same new direction as A3q1. Its coefficient on q1 is
q1⊤Aq3=(Aq1)⊤q3=(α1q1+β1q2)⊤q3=0,
because Aq1 lies in K2(A,q1), to which q3 is orthogonal. Gram–Schmidt has nothing to subtract along q1, and grey marks such vanishing coefficients. The coefficient on q2 is (Aq2)⊤q3=β2, known from Step 2, and the one on q3 is the new number α3=q3⊤Aq3. Two subtractions and a normalization give q4:
Aq3=β2q2+α3q3+β3q4,
and q1,…,q4 is an orthonormal basis of K4(A,q1). Solving this relation for q4, with q2=p1(A)q1 and q3=p2(A)q1, gives q4=p3(A)q1 for the polynomial p3(t)=((t−α3)p2(t)−β2p1(t))/β3 of degree three. The coefficients found so far give the matrix of the Rayleigh quotient on K3(A,q1),
T3=Q3⊤AQ3=α1β10β1α2β20β2α3,
where the grey zeros are the vanishing coefficient q1⊤Aq3 and its symmetric copy. Its smallest and largest eigenvalues are the smallest and largest values of x⊤Ax over the unit vectors of K3(A,q1).
Step 4: Gram–Schmidt on K5(A,q1)=span{q1,Aq1,A2q1,A3q1,A4q1}, and the pattern. The fourth product gives Aq4, which, by the argument of Steps 2 and 3, brings the same new direction as A4q1. Its coefficients on q1 and q2 vanish, since Aq1 and Aq2 lie in K3(A,q1), to which q4 is orthogonal. Its coefficient on q3 is β3, known from Step 3, and the new number is α4=q4⊤Aq4, so that Aq4=β3q3+α4q4+β4q5, and q5=p4(A)q1 for the polynomial p4(t)=((t−α4)p3(t)−β3p2(t))/β4 of degree four. The four steps have computed the coordinates of Aq1,…,Aq4 in the basis:
on q1
on q2
on q3
on q4
on q5
Aq1 (step 1)
α1
β1
Aq2 (step 2)
β1
α2
β2
Aq3 (step 3)
0
β2
α3
β3
Aq4 (step 4)
0
0
β3
α4
β4
Each green length reappears in orange in the next row, by symmetry. The first four columns are the matrix of the Rayleigh quotient on K4(A,q1), which Gram–Schmidt has therefore filled for free:
It is tridiagonal: only the diagonal and its two neighbors are nonzero. The last column of the table, β4 on q5, lies outside K4(A,q1), and it plays a central role below.
The general step. The same reasoning holds at every step j∈{1,2,3,…}, as long as the previous steps have produced qj, which for j≥2 means that the length βj−1=∥rj−1∥ found at step j−1 is positive. For j=1, there is no earlier basis vector, and the conventions q0=0 and β0=0 of Definition 2 let the formulas below cover this case too. By symmetry, the coefficient of Aqj on an earlier basis vector qi is (Aqi)⊤qj, and Aqi lies in Ki+1(A,q1)=span{q1,…,qi+1}, to which qj is orthogonal as soon as i+2≤j. Only two coefficients then survive, i=j−1 and i=j:
The one on qj−1 is the length βj−1=∥rj−1∥ already computed at the previous step, as for β1 in Step 2:
Subtracting both components from Aqj leaves the remainder rj, orthogonal to q1,…,qj, and when its length βj=∥rj∥ is positive, dividing by it gives the next basis vector qj+1:
rj=Aqj−βj−1qj−1−αjqj=βjqj+1.
Step j of Definition 2 below is therefore the following Gram–Schmidt step, which makes the two subtractions one after the other.
qjAAqj−βj−1qj−1zj−αjqjrj÷βjqj+1.
Quantity
What it is
Why it is needed
Aqj
the only product of step j
temporary candidate that carries the new Krylov direction
zj=Aqj−βj−1qj−1
the candidate without its component on qj−1
that coefficient is βj−1, already known from step j−1 by symmetry
αj=qj⊤zj
the coefficient on qj, equal to qj⊤Aqj
removes the last old component and fills the diagonal of Tm
rj=zj−αjqj
the part of Aqj orthogonal to q1,…,qj
the new direction, not yet normalized
βj=∥rj∥
the length of the new direction
normalizes it and fills the off-diagonal of Tm. The value βj=0 means that no new direction exists (exact breakdown).
qj+1=rj/βj
the new basis vector
extends the orthonormal basis by one vector
Rearranging the definitions of zj and rj gives the three-term recurrence
Aqj=βj−1qj−1+αjqj+βjqj+1.
The orange term was found at the previous step, the purple term is the only coefficient computed at step j, and the green term is the new direction. In terms of polynomials, qj+1=pj(A)q1 with p0=1 and
βjpj(t)=(t−αj)pj−1(t)−βj−1pj−2(t) for all j≥1, where p−1=0.
Gram–Schmidt on Km(A,q1) is thus Gram–Schmidt on the polynomials 1,t,t2,… for the inner product ⟨p,r⟩=(p(A)q1)⊤r(A)q1=∑λp(λ)r(λ)∥Πλq1∥2, and the polynomials pj are Everything follows from qj+1=pj(A)q1. First, the polynomials are orthonormal because the basis vectors are: for all i,j≥0, ⟨pi,pj⟩=(pi(A)q1)⊤pj(A)q1=qi+1⊤qj+1=δij. Second, the sum over eigenvalues comes from the expansion Akq1=∑λλkΠλq1 of the introduction. By linearity, p(A)q1=∑λp(λ)Πλq1 for every polynomial p, and since the eigenspaces of distinct eigenvalues are orthogonal, only the terms with μ=λ survive in the product: (p(A)q1)⊤r(A)q1=∑λ,μp(λ)r(μ)(Πλq1)⊤Πμq1=∑λp(λ)r(λ)∥Πλq1∥2. Third, this is an inner product on the polynomials that matter. It is symmetric and bilinear, and ⟨p,p⟩=∑λp(λ)2∥Πλq1∥2 vanishes only if p is zero at every eigenvalue λ with Πλq1=0. If there are s such eigenvalues, a nonzero polynomial of degree less than s cannot vanish at all of them, and these degrees cover p0,…,ps−1, all the polynomials the method produces, since it stops at step s by Proposition 3. Finally, pj has degree j and the positive leading coefficient 1/(β1⋯βj), so p0,…,pj are exactly what Gram–Schmidt produces from 1,t,…,tj for this inner product..
Compared with textbook Gram–Schmidt on q1,Aq1,A2q1,…, the method makes two changes, and neither changes the subspaces.
First, at step j it multiplies the newest basis vector qj by A instead of forming the next power Ajq1. In general, qj has a nonzero coefficient on Aj−1q1 (recall that qj exists only when βj−1=∥rj−1∥>0), since otherwise it would lie in Kj−1(A,q1)=span{q1,Aq1,…,Aj−2q1}=span{q1,…,qj−1}, to which it is orthogonal. By the recurrence for the polynomials pj above, this coefficient is 1/(β1β2⋯βj−1), the leading coefficient of pj−1. Hence Aqj is a nonzero multiple of Ajq1 plus a vector of Kj(A,q1), and the two bring the same new direction, the direction qj+1 orthogonal to Kj(A,q1), as in Step 2, where Aq2 and A2q1 both bring q3. In exact arithmetic, both choices therefore give the same basis. In floating-point arithmetic they do not: the basis computed from the powers Ajq1 loses its orthogonality within a few steps, as the example below shows, and the reason is the following. Since pj has leading coefficient 1/(β1⋯βj), and by the three-term recurrence,
for some vector kj. In each of the two vectors, the green term is orthogonal to Kj(A,q1)=span{q1,…,qj} and the other terms lie in it, so Gram–Schmidt removes the other terms and keeps only the green term. Let x be either vector, Ajq1 or Aqj, and g its green term. What decides the accuracy is the ratio f=∥g∥/∥x∥.
Since g is orthogonal to Kj(A,q1) and x−g lies in it, f is the sine of the angle between x and Kj(A,q1). It lies in [0,1]. It equals 1 when x is orthogonal to Kj(A,q1) and there is nothing to remove, the best case. It is close to 0 when x is close to Kj(A,q1), the bad case, because the computed qj+1 then loses orthogonality to the earlier basis vectors, by an amount that grows as f shrinks, as shown next. Gram–Schmidt removes the terms x−g from x. In double precision, whose unit roundoff is u≈1.1⋅10−16, it does so with rounding errors Gram–Schmidt computes the coefficients ci=qi⊤x of x on the basis vectors of Kj(A,q1) and subtracts the nonzero terms ciqi from x one at a time. There are n=2 of them for Aqj and n=j for Ajq1. Write yk for the intermediate vector left after the first k subtractions, so that y0=x and yn=g. In floating-point arithmetic, every addition, subtraction and multiplication is rounded: its computed result is the exact one times 1+δ, with ∣δ∣≤u. Each operation therefore errs by at most u times the size of its result, and the results here have the size of x, not of g. Indeed, ∣ci∣≤∥qi∥∥x∥=∥x∥ by the Cauchy–Schwarz inequality, and the vectors yk are no longer than x. The difference x−yk, the sum of the first k terms, is the orthogonal projection of x onto the span of the first k basis vectors used, and yk is orthogonal to this span. Pythagoras’ theorem therefore gives ∥x∥2=∥x−yk∥2+∥yk∥2, so ∥yk∥≤∥x∥ (the triangle inequality would only give ∥yk∥≤(k+1)∥x∥ because it ignores that the terms ciqi are orthogonal). An inner product of two vectors of Rd sums d rounded products, so each ci is computed with an error of at most about du∥x∥, and each subtraction of a term adds at most about 2u∥x∥, one u∥x∥ for forming the term and one for subtracting it from yk−1. Up to terms of order u2, the computed green term g therefore satisfies ∥g−g∥≤n(d+2)u∥x∥. The factor n(d+2) is a worst case, and actual errors are usually much smaller, because rounding errors of different signs partly cancel. The number of terms makes this bound j/2 times larger for Ajq1 than for Aqj, which costs only log10(j/2) digits, one digit at j=20. The large difference between the two vectors comes from the ratio f=∥g∥/∥x∥, not from n: at j=20, in the example below, f is 8⋅10−12 for A20q1 and 0.4 for Aq20. What matters here is that the bound is proportional to ∥x∥, not to ∥g∥: a small g does not come with a small error., and these errors stay in the computed g. Relative to ∥g∥=f∥x∥, they are therefore of the order of u/f, and so are the components that the computed qj+1 keeps on the earlier basis vectors. This loss of accuracy is called cancellation. The same thing happens with a number in place of the vector x. Subtracting 1 from x=1.000001 leaves g=10−6, so f=g/x≈10−6, and the error of about 10−16 with which x is stored stays in g, where it is a relative error of about u/f=10−10.
The two ratios behave very differently. For x=Aqj, the three terms are orthogonal, so ∥x∥2=βj−12+αj2+βj2 and
f=βj−12+αj2+βj2βj,
which depends only on the coefficients of step j, and is close to 0 only if βj is small compared with βj−12+αj2. For x=Ajq1, the green term has length β1β2⋯βj, so
f=∥Ajq1∥β1β2⋯βj.
Up to scaling, the powers q1,Aq1,A2q1,… are the iterates of power iteration started from q1 (Section 3), which, under the hypotheses of Theorem 1, converge up to sign to the dominant eigenvector, so the angle between Ajq1 and the previous power Aj−1q1, which lies in Kj(A,q1), tends to 0. As a result, this ratio shrinks Write vk=Akq1/∥Akq1∥ for the iterates of power iteration started from q1, and s=sgn(λ⋆). The green term of vj has length f and is orthogonal to Kj(A,q1), so f is the distance from vj to this subspace, which contains svj−1. Under the hypotheses of Theorem 1, skvk lies within (ρ2/∣λ⋆∣)ktan(ϕ0) of q⋆ for every k, and since s2=1, f≤∥vj−svj−1∥=∥sjvj−sj−1vj−1∥≤∥sjvj−q⋆∥+∥q⋆−sj−1vj−1∥≤2(∣λ⋆∣ρ2)j−1tan(ϕ0). This bound uses only the last power in Kj(A,q1). Using all of them gives the much faster decrease seen in practice. with j, and in practice much faster, since Kj(A,q1) contains all the earlier powers, not only the last one. For example, we ran both methods in double precision on the diagonal matrix of size 1000 with eigenvalues drawn uniformly from [0,1], from a random unit vector q1, and measured the ratios:
step j
1
5
10
20
f for Ajq1
0.5
4⋅10−3
6⋅10−6
8⋅10−12
f for Aqj
0.5
0.4
0.4
0.4
The same runs also measured how far each computed basis is from orthonormal. Let Qm be the matrix whose columns are the first m computed basis vectors. The entries of Qm⊤Qm−Im are the numbers qi⊤qk−δik, so its norm is 0 for an exactly orthonormal basis, and close to 1 when two of the basis vectors are nearly parallel. The measured norms are:
number m of basis vectors
10
20
25
30
basis from the powers Ajq1
9⋅10−11
4⋅10−3
0.9
0.9
basis from the Lanczos recurrence, without reorthogonalization
3⋅10−15
5⋅10−15
5⋅10−15
6⋅10−15
The basis built from the powers thus loses its orthogonality step by step and is far from orthonormal after 25 vectors, while the Lanczos basis stays orthonormal to within 10−14.
Second, it skips the coefficients that vanish. Full Gram–Schmidt on Aqj computes the j coefficients qi⊤Aqj, i=1,…,j, and subtracts the j terms, which costs O(jd) operations at step j, a cost that grows with j. By symmetry, the coefficients with i≤j−2 are zero, as shown in the general step above. The Lanczos method therefore computes only αj, reuses βj−1, and subtracts two terms, so that every step costs one product and O(d) further operations, whatever j.
This saving has a price in floating-point arithmetic. The proof that qi⊤Aqj=0 for i≤j−2 uses that qj is exactly orthogonal to q1,…,qj−1, and the computed basis vectors are orthogonal only up to rounding errors. The skipped coefficients are therefore small but not zero, and since the method does not remove them, they stay in the computed rj. For i≤j−2, the new basis vector qj+1=rj/βj then has, up to the rounding errors of step j itself, the component
qi⊤qj+1=βjqi⊤Aqj−βj−1qi⊤qj−1−αjqi⊤qj
on the earlier vector qi, which is zero in exact arithmetic. Its numerator is made of inner products between earlier basis vectors: qi⊤qj−1 and qi⊤qj directly, and qi⊤Aqj=(Aqi)⊤qj through the recurrence Aqi=βi−1qi−1+αiqi+βiqi+1, up to rounding. The loss of orthogonality of the previous steps thus feeds into every new basis vector, and it accumulates. The same runs as in the first item measured ∥Qm⊤Qm−Im∥ for the Lanczos basis, without reorthogonalization and with the reorthogonalization of Algorithm 2 described below:
number m of basis vectors
30
100
150
200
without reorthogonalization
6⋅10−15
3⋅10−12
6⋅10−10
1⋅10−6
with reorthogonalization
9⋅10−16
1⋅10−15
1⋅10−15
1⋅10−15
Without reorthogonalization, the loss is much slower than for the powers, where it reached 0.9 after 25 vectors, but it does not stop. This is why Algorithm 2 below reorthogonalizes: after the three-term step, it removes from rj its components on all the stored basis vectors q1,…,qj once more, here in two passes, and the basis then stays orthonormal to about 10−15 at every m of the table. This gives back the O(jd) cost per step that symmetry had saved. Unlike Gram–Schmidt on the powers, however, this extra pass is accurate. The components it removes have the size of the loss of orthogonality, tiny compared with ∥rj∥, so the ratio f of the first item is close to 1 for it.
Eigenvalues and eigenvectors from the coefficients. After m steps, the coefficients fill the m×m tridiagonal matrix
and the m recurrences Aqj=βj−1qj−1+αjqj+βjqj+1, for j=1,…,m, combine into Compare the columns. Column j of AQm is Aqj. Column j of QmTm is Qm times column j of Tm, whose only nonzero entries are βj−1, αj and βj, in rows j−1, j and j+1. It is therefore βj−1qj−1+αjqj+βjqj+1, without the first term when j=1 since β0=0, which is Aqj by the recurrence. The last column is the exception: Tm has no row m+1, so column m of QmTm stops at βm−1qm−1+αmqm and misses the term βmqm+1 of the recurrence. The green matrixβmqm+1em⊤ adds it back: its columns are zero except the last, which is βmqm+1.:
AQm=QmTm+βmqm+1em⊤,
where em is the last standard basis vector of Rm. Inside Km(A,q1)=span{q1,…,qm}, the matrix A acts on coordinates as the small matrix Tm, except for the green term, through which the last basis vector leaks out of the subspace. Every eigenpair (ϑ,y) of Tm, with ∥y∥=1, therefore gives a candidate eigenpair (ϑ,Qmy) of A. Such a candidate is called a Ritz pair: the Ritz valueϑ estimates an eigenvalue of A, and the Ritz vectorQmy, a unit vector of Km(A,q1), estimates a corresponding eigenvector. Lemma 4 states their properties. In short,
d×dA⟶m×mTm=Qm⊤AQm⟶(ϑ,y)⟶(ϑ,Qmy).
Multiplying the matrix identity AQm=QmTm+βmqm+1em⊤ above on the right by y gives
AQmy=QmTmy+βmqm+1(em⊤y)=ϑQmy+βmymqm+1.
The residual of the Ritz pair is therefore exactly the leak, and since qm+1 is a unit vector,
AQmy−ϑQmy=βm∣ym∣.
The residual is available without any further product, and Lemma 2 turns it into a certificate. The basis vectors are building blocks, not approximate eigenvectors. An approximate eigenvector Qmy combines all of them, so the basis must be stored to form it, whereas the Ritz values need only the numbers αj and βj.
Gram–Schmidt also explains why the Ritz values approach eigenvalues. The argument has three steps.
First, the green term of Amq1 has the smallest Euclidean norm among the vectors of a certain form. As shown above with j=m,
Amq1=β1β2⋯βmqm+1+km with km∈Km(A,q1).
Every k∈Km(A,q1) is a combination of q1,Aq1,…,Am−1q1, say k=∑i=0m−1aiAiq1, so
Amq1−k=π(A)q1
with π(t)=tm−∑i=0m−1aiti a A polynomial of degree m is called monic if its coefficient of tm is 1, i.e. if it has the form t↦tm+∑i=m−10citi. of degree m. Conversely, any monic polynomial tm+∑i=0m−1citi of degree m arises in this way from ai=−ci. The vectors Amq1−k with k∈Km(A,q1) are therefore exactly the vectors π(A)q1 with π monic of degree m. Each of them is the green term β1β2⋯βmqm+1 plus km−k∈Km(A,q1), two orthogonal vectors, so Pythagoras’ theorem gives ∥Amq1−k∥2=(β1⋯βm)2+∥km−k∥2, which is smallest for k=km, where the vector is the green term itself. Among the vectors π(A)q1 with π monic of degree m, the green term thus has the smallest Euclidean norm ∥π(A)q1∥. Its own polynomial follows from qm+1=pm(A)q1: the green term is πm(A)q1 for πm(t)=β1⋯βmpm(t), which is monic because pm has leading coefficient 1/(β1⋯βm). Hence πm minimizes ∥π(A)q1∥ over the monic polynomials of degree m. Recall the inner product on polynomials defined with the three-term recurrence, ⟨p,r⟩=(p(A)q1)⊤r(A)q1=∑λp(λ)r(λ)∥Πλq1∥2. It gives ∥π(A)q1∥2=⟨π,π⟩, so this Euclidean norm of the vector π(A)q1 is also the norm of the polynomial π for that inner product, and with its expression as a sum over the eigenvalues, the minimization reads
In words, πm is the monic polynomial of degree m that is smallest on the eigenvalues of A, each eigenvalue weighted by ∥Πλq1∥2, the squared length of the component of q1 in its eigenspace.
Second, the Ritz values at step m are exactly the roots of πm, because Let Dk(t)=det(tIk−Tk), with D0=1 and D−1=0. The matrix tIk−Tk is tridiagonal, and its last row has only two nonzero entries, −βk−1 and t−αk. Expanding the determinant along this row, and the remaining minor along its last column, gives Dk(t)=(t−αk)Dk−1(t)−βk−12Dk−2(t). Multiplying the recurrence βkpk=(t−αk)pk−1−βk−1pk−2 by β1⋯βk−1 shows that Pk=β1⋯βkpk satisfies the same recurrence, Pk(t)=(t−αk)Pk−1(t)−βk−12Pk−2(t), with the same starting values P0=1 and P−1=0. Hence Dk=Pk for every k, and in particular det(tIm−Tm)=Pm=πm.. If ϑ1,…,ϑm are the Ritz values, then πm(t)=∏i=1m(t−ϑi).
Third, the two steps combine into one statement. Any m real numbers z1,…,zm define the monic polynomial ∏i=1m(t−zi) of degree m, and πm is of this form with zi=ϑi, since Tm is symmetric and its eigenvalues are real. Since πm minimizes over all monic polynomials of degree m, the Ritz values minimize in particular, over all (z1,…,zm)∈Rm, the sum
λ∈Λ(A)∑∥Πλq1∥2i=1∏m(λ−zi)2.
Each term is the weight of an eigenvalue λ times the product of its squared distances to z1,…,zm. A term with a large weight is reduced by placing one of the zi close to its eigenvalue, and this is why the Ritz values approach the eigenvalues that q1 contains. Two statements make this exact. When m equals the number s of eigenvalues with Πλq1=0, choosing the zi equal to these eigenvalues makes the sum zero, so the minimum is zero, β1⋯βs=0, hence βs=0, and the Ritz values are exactly these eigenvalues (Proposition 3). Before that, Lemma 5 shows that the extreme Ritz values move monotonically toward λmin(A) and λmax(A).
4.4 A three-dimensional example
The steps of Section 4.3 can be followed with numbers on a small matrix. Take
A=100020003,q1=31111.
Step 1. The candidate Aq1=31[1,2,3]⊤ has the coefficient α1=q1⊤Aq1=2 on q1, and there is no q0 to remove. The rest, r1=31[−1,0,1]⊤, has length β1=2/3 and gives q2=21[−1,0,1]⊤.
Step 2. The candidate Aq2=21[−1,0,3]⊤ has the coefficient q1⊤Aq2=2/3=β1 on q1, as symmetry predicts, and the coefficient α2=2 on q2. Removing both leaves r2=321[1,−2,1]⊤, of length β2=1/3, which gives q3=61[1,−2,1]⊤.
Step 3. The candidate Aq3=61[1,−4,3]⊤ minus β2q2 equals 61[2,−4,2]⊤=2q3. Hence α3=2 and r3=0, an exact breakdown: K3(A,q1) is already all of R3, so no new direction exists.
and their eigenvalues, the Ritz values, are {2}, then {2−2/3,2+2/3}≈{1.18,2.82}, then {1,2,3}. The extreme Ritz values move outward toward the extreme eigenvalues 1 and 3, as Lemma 5 guarantees. At m=3, they are exactly the eigenvalues of A, as Proposition 3 guarantees for a starting vector with a nonzero component on every eigenvector. At m=2, the eigenvectors of T2 are 21[1,±1]⊤, so both Ritz pairs have residual β2∣y2∣=31⋅21≈0.41 (Lemma 4). This bounds their actual distances to the spectrum, 0.18 in both cases (Lemma 2). For comparison, power iteration from the same q1 estimates λmax(A)=3 after the same two products by v1⊤Av1=36/14≈2.57, with v1=Aq1/∥Aq1∥, which is farther from 3 than the Ritz value 2.82. In the polynomial reading of Section 4.3, the Ritz values 2±2/3 are the roots of π2(t)=(t−2)2−2/3, the monic quadratic that makes ∥π(A)q1∥ smallest, and at m=3 the polynomial (t−1)(t−2)(t−3) makes it zero, which is the breakdown.
4.5 Definition 2 [Lanczos recurrence and exact breakdown]
If βj>0, define qj+1=rj/βj and continue. If βj=0, stop the current run. We refer to this event as an exact breakdown. As shown below, it means that the generated subspace is invariant under A.
■□
The vectors qj are basis vectors, not generally eigenvectors of A. Each step orthogonalizes Aqj against the two most recent basis vectors only, for the reasons given in Section 4.3. This Orthogonalizing Aqj against every previous basis vector, as in Gram–Schmidt, would cost more at every step. For symmetric A, the coefficients against q1,…,qj−2 vanish (Lemma 3), so the two coefficients αj and βj−1 suffice. is enough because A is symmetric. The following lemma shows that, in exact arithmetic, Aqj is automatically orthogonal to q1,…,qj−2, so the projected matrix is tridiagonal and each step costs one product plus O(d) further operations.
4.6 Lemma 3 [Orthogonality and tridiagonal projection]
Before and including the first exact breakdown, the Lanczos vectors satisfy
qi⊤qj=δij,span{q1,…,qm}=Km(A,q1).
Let Qm=[q1,…,qm], and let Tm∈Rm×m be the symmetric tridiagonal matrix
that is, (Tm)ii=αi for 1≤i≤m, (Tm)i,i+1=(Tm)i+1,i=βi for 1≤i<m, and all other entries equal to zero. Then
AQm=QmTm+rmem⊤.
Moreover,
Qm⊤rm=0,Tm=Qm⊤AQm.
Remark 12 [What the tridiagonal matrix represents]
The identity Tm=Qm⊤AQm says that Tm represents A compressed to the Krylov subspace, in the orthonormal basis q1,…,qm. The rank-one term rmem⊤ measures how far that subspace is from being invariant under A: only the image of the last basis vector leaves it.
Proof
Orthogonality. We argue by induction. The initial vector is normalized by assumption. Suppose that q1,…,qj are orthonormal and that the previous iterations have been completed. We show that rj is orthogonal to each qi with i≤j.
The case i=j. Since qj⊤qj−1=0, with q0=0 when j=1, we have
αj=qj⊤zj=qj⊤(Aqj−βj−1qj−1)=qj⊤Aqj.
Therefore,
qj⊤rj=qj⊤zj−αjqj⊤qj=αj−αj=0.
The case i=j−1, for j≥2. The previous recurrence gives
Aqj−1=βj−2qj−2+αj−1qj−1+βj−1qj.
Using symmetry,
The first line expands zj=Aqj−βj−1qj−1, and the two inner products vanish by the induction hypothesis, since i differs from j−1 and j. The second line moves A onto qi by symmetry and expands Aqi=βi−1qi−1+αiqi+βiqi+1 by the recurrence at step i, with q0=0 when i=1, and its three terms are orthogonal to qj, as noted above.
The blue terms are where symmetry enters. Moving A from qj onto an older basis vector qi turns the inner product into one with Aqi, and an earlier step of the recurrence has already expanded Aqi in the basis. This is why orthogonalizing against the two most recent vectors suffices.
Thus rj is orthogonal to every column of Qj. If βj>0, normalizing rj extends the orthonormal family.
Krylov subspaces. We prove span{q1,…,qm}=Km(A,q1) by induction on m. It holds at m=1, since both sides equal span{q1}. Suppose it holds at m=j and βj>0. Then
βjqj+1=Aqj−βj−1qj−1−αjqj∈Kj+1(A,q1),
because qj−1 and qj belong to Kj(A,q1)⊆Kj+1(A,q1), and AqjSince qj∈Kj(A,q1), there are numbers c0,…,cj−1 with qj=∑k=0j−1ckAkq1, so Aqj=∑k=0j−1ckAk+1q1 is a combination of Aq1,…,Ajq1, which belong to Kj+1(A,q1) by Definition 1.. Consequently,
span{q1,…,qj+1}⊆Kj+1(A,q1).
The left-hand side has dimension j+1, because its vectors are orthonormal, so j+1≤dimKj+1(A,q1). We also have dimKj+1(A,q1)≤j+1, since Kj+1(A,q1) is generated by j+1 vectors. We thus have dimKj+1(A,q1)=dimspan{q1,…,qj+1}=j+1. The two spaces are therefore equal, since span{q1,…,qj+1} is contained in Kj+1(A,q1) and both have the same finite dimension.
Matrix identity. For each completed iteration j<m,
Aqj=βj−1qj−1+αjqj+βjqj+1.
For the final column,
Aqm=βm−1qm−1+αmqm+rm.
Collecting these identities column by column gives
AQm=QmTm+rmem⊤.
Multiplying on the left by Qm⊤, and using
Qm⊤Qm=Im and Qm⊤rm=0 gives Qm⊤AQm=Tm.
■□
The residual form of this identity remains valid when βm=0. When βm>0, it can also be written as
AQm=QmTm+βmqm+1em⊤.
No vector qm+1 needs to be defined at an exact breakdown.
4.7 Lemma 4 [Ritz pairs and their residuals]
Let
Tmy=ϑy,∥y∥=1,
and define
x=Qmy.
Then
∥x∥=1,x⊤Ax=ϑ,
and
Ax−ϑx=rmym,∥Ax−ϑx∥=βm∣ym∣.
Here ym=em⊤y is the last coordinate of y. In particular,
λ∈Λ(A)min∣ϑ−λ∣≤βm∣ym∣.
The scalar ϑ is called a Ritz value, and x a Ritz vector.
Remark 13 [Residuals at no extra cost]
Ritz pairs are the eigenpairs of the compression of A to the Krylov subspace, lifted back to Rd. Computing them is the Approximating the eigenpairs of A by those of its compression Q⊤AQ to a subspace with orthonormal basis Q is the Rayleigh–Ritz procedure. Lanczos applies it to Krylov subspaces, where the compression is tridiagonal. on a Krylov subspace. Their residuals cost nothing extra: βm and the last coordinate of y are available from the small problem. A Ritz pair can therefore be certified as accurate even when βm is not small, as long as ∣ym∣ is.
Proof
Orthonormality of the columns of Qm gives
∥x∥2=y⊤Qm⊤Qmy=y⊤y=1.
Also,
x⊤Ax=y⊤Qm⊤AQmy=y⊤Tmy=ϑy⊤y=ϑ.
Using Lemma 3,
Ax−ϑx=AQmy−ϑQmy=QmTmy+rmem⊤y−ϑQmy=Qm(Tmy−ϑy)+rmym=rmym.
The blue factor vanishes because y is an eigenvector of Tm. Taking norms gives
∥Ax−ϑx∥=∥rm∥∣ym∣=βm∣ym∣.
The eigenvalue bound minλ∈Λ(A)∣ϑ−λ∣≤βm∣ym∣ is Lemma 2 applied to the unit vector x, whose Rayleigh quotient is ϑ.
The eigenbasis expansion, written out for the Ritz vector
Lemma 2 states that every unit vector q, with Rayleigh quotient λ=q⊤Aq and residual r=Aq−λq, satisfies minλ∈Λ(A)∣λ−λ∣≤∥r∥. Take q=x. It is a unit vector, its Rayleigh quotient is λ=x⊤Ax=ϑ, and its residual is r=Ax−ϑx, all by the computations above. Lemma 2 therefore gives
minλ∈Λ(A)∣ϑ−λ∣≤∥Ax−ϑx∥=βm∣ym∣.
More specifically, expand x in an orthonormal eigenbasis u1,…,ud of A: x=∑i=1dciui, Aui=λiui, ∑i=1dci2=1.
Then Ax−ϑx=∑i=1dci(λi−ϑ)ui,
and therefore
This bound certifies proximity to some eigenvalue of A. By itself, a small residual does not certify that the eigenvalue is a global extremum or that all eigenvalues have been found.
4.8 Lemma 5 [Variational characterization and monotonicity]
The extreme Ritz values are the extreme Rayleigh quotients over a growing family of subspaces. They can therefore only move outward, toward the extreme eigenvalues of A, and never overshoot them, as the green curves of panel (a) in Figure 8 show.
Proof
Every vector in Km(A,q1) has a unique representation
x=Qmy.
Moreover,
∥x∥2=y⊤Qm⊤Qmy=∥y∥2,
and
x⊤Ax=y⊤Tmy.
Let s1,…,sm be an orthonormal eigenbasis of Tm, with eigenvalues ordered as ϑ1≤ϑ2≤⋯≤ϑm, so that ϑmin(m)=ϑ1 and ϑmax(m)=ϑm. For any unit vector
y=∑i=1mcisi, ∑i=1mci2=1,
we have y⊤Tmy=∑i=1mϑici2.
Thus ϑ1≤y⊤Tmy≤ϑm, with equality for y=s1 and y=sm respectively. This proves the two variational identities.
Since Km(A,q1)⊆Km+1(A,q1),
minimizing over the larger set can only decrease the minimum, while maximizing over the larger set can only increase the maximum. Therefore,
ϑmin(m+1)≤ϑmin(m) and ϑmax(m)≤ϑmax(m+1).
Finally, every unit vector x∈Rd satisfies
λmin(A)≤x⊤Ax≤λmax(A),
as follows from the same eigenbasis expansion for A. Apply this to the Ritz vectors x=Qms1 and x=Qmsm. They are unit vectors, since ∥Qmsi∥=∥si∥=1, and their Rayleigh quotients are s1⊤Tms1=ϑ1=ϑmin(m) and sm⊤Tmsm=ϑm=ϑmax(m). Hence λmin(A)≤ϑmin(m) and ϑmax(m)≤λmax(A) for every m, and with m+1 in place of m, these are the remaining bounds.
■□
This characterization provides a precise connection with power iteration. Started from v0=q1, and whenever Am−1q1=0, the normalized power iterate
vm−1=Am−1q1/∥Am−1q1∥
belongs to Km(A,q1). Consequently, ϑmin(m)≤vm−1⊤Avm−1≤ϑmax(m).
Lanczos optimizes the Rayleigh quotient over all unit vectors in the available Krylov subspace. The normalized power vector is one particular member of that subspace.
Both estimates cost m matrix–vector products: Tm needs Aq1,…,Aqm, and the Rayleigh quotient of vm−1 needs Av0,…,Avm−1. The eigenvalue of largest magnitude is λmax(A) or λmin(A). In the first case, vm−1⊤Avm−1≤ϑmax(m)≤λmax(A), and in the second, λmin(A)≤ϑmin(m)≤vm−1⊤Avm−1. Either way, at an equal number of products, the extreme Ritz value on the side of the dominant eigenvalue is at least as close to it as the power-iteration estimate (middle panel of the figure in Section 4.10).
This variational statement does not imply that Lanczos has a lower total runtime for every matrix. Constructing and retaining the basis introduces additional work and storage.
4.9 Proposition 3 [Exact termination and the eigenvalues recovered by one run]
For each λ∈Λ(A), let Πλ be the orthogonal projection onto the eigenspace Eλ=ker(A−λId). Define
Λ1={λ∈Λ(A):Πλq1=0} and let s=∣Λ1∣.
Here Λ1 contains distinct eigenvalues, irrespective of their multiplicities.
Then, a single run of the Lanczos recurrence of Definition 2 terminates at exactly iteration s:
βj>0(1≤j<s),βs=0.
Moreover, at termination, AQs=QsTs and Λ(Ts)=Λ1. Since Ts has size s=∣Λ1∣, each of these eigenvalues is simple for Ts.
Remark 15 [One direction per distinct eigenvalue]
A run can only discover the eigenvalues present in q1, and it discovers each of them once, however large its multiplicity. Inside the Krylov subspace, the whole eigenspace Eλ is represented by the single vector Πλq1. The number of steps before an exact breakdown is therefore the number of distinct eigenvalues present in q1, not d. In Figure 8, below, the matrix of size 60 has 12 distinct eigenvalues: βm falls from about 0.3 at m=11 to about 10−12 at m=12 (panel (c)), and the 12 Ritz values are then the 12 distinct eigenvalues, each appearing once (panel (a)).
Figure 8. The Lanczos method on a symmetric matrix A=Udiag(λ1,…,λ60)U⊤ of size d=60 with only s=12 distinct eigenvalues, of multiplicities 1 to 9 (labels on the right of panel a), including the eigenvalue 0 with multiplicity 8. The numbers λ1,…,λ60 are these eigenvalues, each repeated according to its multiplicity, and the eigenvectors, the columns of U, are random: U is uniformly distributed over the orthogonal matrices, and is sampled as the orthogonal factor of the QR decomposition of a 60×60 matrix with independent standard Gaussian entries (with the signs chosen so that the triangular factor has a positive diagonal). (a) One run: Ritz values at each step m (black dots) against the distinct eigenvalues (grey lines). The extreme Ritz values (green) move outward monotonically (Lemma 5), the power-iteration estimate after the same number of products (orange) stays between them, and at m=s=12 the Ritz values are exactly the 12 distinct eigenvalues, each appearing once (Proposition 3). (b) Distance from λmax(A)=2 of the largest Ritz value (green) and of the power-iteration estimate (orange), at equal numbers of matrix–vector products. (c) The coefficient βm falls from about 0.3 at m=11 to about 10−12 at m=12: the exact breakdown of Proposition 3, up to floating-point error. For every starting vector, the 12 Ritz values at m=12 match the 12 distinct eigenvalues to within 3⋅10−15. Panels (b) and (c) show medians over 200 Gaussian starting vectors, with bands containing the middle half of the runs (the interquartile range), and all runs use full reorthogonalization.
Proof
Write pλ=Πλq1 for λ∈Λ1,
and define W=span{pλ:λ∈Λ1}.
The vectors pλ are nonzero and mutually orthogonal, so dimW=s.
The Krylov subspaces lie in W. Every power of A applied to q1 belongs to W, because q1=∑λ∈Λ1pλ and, for k≥1,
Akq1=∑λ∈Λ1λkpλ.
Therefore,
Km(A,q1)⊆W for every m.
The Krylov subspace of order s equals W. For each λ∈Λ1, define the ℓλ is the unique polynomial of degree at most s−1 that equals 1 at λ and 0 at the other s−1 points of Λ1. Applied to A, it keeps the component of q1 along Eλ and removes every other component present in q1.
ℓλ(t)=μ∈Λ1μ=λ∏λ−μt−μ.
When s=1, the empty product is defined as 1. This polynomial has degree at most s−1 and satisfies
ℓλ(μ)=δλ,μ for all μ∈Λ1, with the Kronecker delta of the notation.
Consequently,
ℓλ(A)q1=∑μ∈Λ1ℓλ(μ)pμ=pλ.
Since ℓλ has degree at most s−1, pλ∈Ks(A,q1).
Thus W⊆Ks(A,q1), and hence Ks(A,q1)=W.
No breakdown before iteration s. Suppose a breakdown occurred at some m<s. Lemma 3 would give
AQm=QmTm,
so
Km(A,q1)=range(Qm)
would be invariant under A. Since it contains q1, it would then contain every power Akq1 with k∈{0,1,2,…}, not only the powers with k≤m−1 that span it by definition. In particular, it would contain q1,Aq1,…,As−1q1, and therefore
Ks(A,q1)⊆Km(A,q1).
This is impossible, because
dimKs(A,q1)=s>m=dimKm(A,q1).
Thus no breakdown occurs before iteration s.
Breakdown at iteration s. At iteration s, W=Ks(A,q1)=range(Qs). The first equality was shown above, and the second is the identity span{q1,…,qs}=Ks(A,q1) of Lemma 3, which applies because no breakdown occurs before iteration s. Since W is spanned by eigenvectors, it is invariant under A. Therefore,
rs=Aqs−βs−1qs−1−αsqs∈W.
But Lemma 3 also gives
rs⊥range(Qs)=W.
Hence rs=0, and βs=0. The identity AQs=QsTs+rses⊤ of Lemma 3 then becomes AQs=QsTs.
The eigenvalues of Ts. The normalized vectors
uλ=pλ/∥pλ∥ for
λ∈Λ1 form an orthonormal eigenbasis of W.
Set yλ=Qs⊤uλ.
Transposing AQs=QsTs, with A and Ts symmetric, gives Qs⊤A=TsQs⊤, so
Tsyλ=TsQs⊤uλ=Qs⊤Auλ=λQs⊤uλ=λyλ.
Moreover, uμ∈W=range(Qs), so QsQs⊤uμ=uμ, and for λ,μ∈Λ1,
yλ⊤yμ=uλ⊤QsQs⊤uμ=uλ⊤uμ=δλ,μ.
Thus the s vectors yλ, λ∈Λ1, are orthonormal, hence an orthonormal basis of Rs, and each of them is an eigenvector of Ts. Every λ∈Λ1 is therefore an eigenvalue of Ts, with the nonzero eigenvector yλ. Conversely, let Tsz=ξz with z=0, and expand z=∑λ∈Λ1cλyλ in this basis. Then
0=Tsz−ξz=∑λ∈Λ1cλ(λ−ξ)yλ,
so cλ(λ−ξ)=0∀λ∈Λ1, by linear independence. Some cλ is nonzero, since z=0, so ξ equals that λ. The elements of Λ1 being distinct, cμ=0∀μ=ξ, and z is a multiple of yξ. Hence Λ(Ts)=Λ1, and each eigenvalue of Ts is simple, with eigenvectors the multiples of a single yλ.
■□
4.10 Corollary 3 [Random initialization and repeated eigenvalues]
If g∼N(0,Id) and q1=g/∥g∥, then, with probability one, a terminating Lanczos run recovers every distinct eigenvalue of A. If A has d distinct eigenvalues, the run reaches s=d and yields a complete eigenbasis. If eigenvalues are repeated, a single run yields only one eigenvector from each corresponding eigenspace. With probability one, the eigenvector found in Eλ at termination is, up to sign, the normalized projection of q1 onto Eλ.
Proof
For each λ∈Λ(A), the event Πλg=0 means that g lies in the proper linear subspace Eλ⊥, which has Gaussian probability zero. Since there are finitely many distinct eigenvalues, the For finitely many events, the probability that at least one of them occurs is at most the sum of their probabilities. A finite union of probability-zero events therefore has probability zero. gives
P(∃λ∈Λ(A):Πλg=0)≤∑λ∈Λ(A)P(Πλg=0)=0.
Thus Λ1=Λ(A) with probability one, and Proposition 3 applies.
For any polynomial p, Πλp(A)q1=p(λ)Πλq1. Therefore, the projection of the entire Krylov subspace onto Eλ is contained in span{Πλq1}, which is one-dimensional when Πλq1=0. A single run cannot recover additional independent directions in that eigenspace. At termination, the proof of Proposition 3 gives, for each λ∈Λ1, the unit eigenvector yλ=Qs⊤uλ of Ts, where uλ=Πλq1/∥Πλq1∥ is the normalized projection of q1 onto Eλ. Its Ritz vector is
Qsyλ=QsQs⊤uλ=uλ,
where the last equality holds because QsQs⊤ is the orthogonal projection onto range(Qs), which contains uλ. Since λ is a simple eigenvalue of Ts, every unit eigenvector of Ts for λ is ±yλ, so the eigenvector found in Eλ is ±uλ.
When the spectrum is simple and every eigenvalue is present in the initialization, s=d. In this case, Qd is an orthogonal matrix, and
A=QdTdQd⊤.
An orthonormal eigenbasis of Td therefore yields a complete orthonormal eigenbasis of A after multiplication by Qd.
■□
An early exact breakdown does not necessarily mean that a complete eigenbasis of A has been found. It means that the current Krylov subspace has become invariant under A. For example, when A=cId, every initialization produces a breakdown after one iteration, even though the eigenspace has dimension d.
To recover all eigendirections, the construction can continue with new starting vectors orthogonal to the previously generated basis.
4.11 Proposition 4 [Complete decomposition by continuation with new starting vectors]
Let A∈Rd×d be symmetric. Initialize q0=0, β0=0, and any unit vector q1. Run the Lanczos recurrence of Definition 2 for j=1,…,d, with the following rule whenever j<d:
qj+1={rj/βj,any unit vector orthogonal to q1,…,qj,βj>0,βj=0.
At a breakdown, retain βj=0 and all previously generated basis vectors.
Then the following hold.
Orthogonal basis. The matrix Qd=[q1,…,qd] is orthogonal, Qd⊤Qd=Id.
Tridiagonal form.βd=0 and AQd=QdTd, hence A=QdTdQd⊤. The matrix Td is symmetric and tridiagonal, with a zero off-diagonal entry at every boundary between runs.
Complete eigenbasis. If Td=ZDZ⊤ with Z⊤Z=Id and D∈Rd×d the diagonal matrix of the eigenvalues of Td, then V=QdZ is a complete orthogonal eigenvector matrix of A:
V⊤V=Id and AV=VD.
Remark 16 [Why continuation completes the basis]
Each run discovers one direction in every eigenspace that its starting vector reaches, then stops at a subspace invariant under A (Proposition 3). Continuing with a new starting vector in the orthogonal complement, which is invariant under A as well, lets the next run find the directions that the previous runs could not see. Figure 9 shows the resulting basis and tridiagonal matrix on an example.
Figure 9. The structure built by Proposition 4, on an example with d=9 and runs of lengths 4, 3 and 2. Left: the columns of Q9, coloured by run. Each run ends with an exact breakdown, β4=0, β7=0 and β9=0, and the next run starts from Pjg/∥Pjg∥, a unit vector orthogonal to q1,…,qj, chosen as described after the proof. Right: the matrix T9 of the identity AQ9=Q9T9. Its diagonal blocks, framed in the colour of their run, are the tridiagonal matrices of the three runs, and the entries that would couple two runs are β4=β7=0 (grey).
Proof
We follow the runs one after the other and prove, by induction over the runs, the following claim: if a run ends with an exact breakdown after iteration j, then q1,…,qj are orthonormal and AQj=QjTj. The three items then follow from it.
The first run. At its exact breakdown, after iteration j, Lemma 3 gives orthonormal vectors q1,…,qj and AQj=QjTj+rjej⊤ with rj=0, since βj=∥rj∥=0, that is, AQj=QjTj.
From one run to the next. Suppose that a run ends with an exact breakdown after iteration j<d, with q1,…,qj orthonormal and AQj=QjTj. Then range(Qj) and its orthogonal complement are both A subspace S is invariant under A if Ay∈S for every y∈S. The range. If y∈range(Qj), there is x∈Rj with y=Qjx, so Ay=AQjx=QjTjx∈range(Qj).Its orthogonal complement. The orthogonal complement of range(Qj) is ker(Qj⊤). In fact, y∈range(Qj)⊥⟺y⊤Qjx=(Qj⊤y)⊤x=0∀x∈Rj⟺Qj⊤y=0. Now let y∈ker(Qj⊤). The symmetry of A and of Tj gives Qj⊤Ay=(AQj)⊤y=(QjTj)⊤y=TjQj⊤y=0, so Ay∈ker(Qj⊤)..
The next run starts from a unit vector qj+1 of this complement, which exists because the complement has dimension d−j>0. Since βj=0, its first step computes zj+1=Aqj+1, with no term in qj, so the run follows Definition 2 from qj+1, and Lemma 3 applies to it. Its vectors are therefore orthonormal and lie in the Krylov subspaces Km(A,qj+1). These subspaces lie in the complement, because it contains qj+1 and is invariant under A, so the vectors of the new run are orthogonal to q1,…,qj. If the new run ends with an exact breakdown after iteration j′, Lemma 3 also gives A[qj+1,…,qj′]=[qj+1,…,qj′]T′, where T′ is the tridiagonal matrix of this run, with no residual term since βj′=0. Together with AQj=QjTj, this gives AQj′=Qj′Tj′, because the entries of Tj′ that couple the two runs are βj=0. This proves the claim for the new run.
1. Orthogonal basis. At every iteration j<d, the rule produces a unit vector qj+1, either rj/βj or a new starting vector, so the construction reaches qd. By the previous paragraph, the vectors of each run are orthonormal and orthogonal to those of the earlier runs. Hence q1,…,qd are orthonormal, that is, Qd⊤Qd=Id.
2. Tridiagonal form. At iteration d, the residual rd=Aqd−βd−1qd−1−αdqd of the last run is orthogonal to the vectors of that run, by Lemma 3. It also lies in the orthogonal complement of the earlier runs, which is invariant under A and contains the vectors of the last run, since rd is a combination of Aqd and of vectors of the last run (if qd−1 belongs to an earlier run, then βd−1=0). Thus rd is orthogonal to every column of Qd. These columns form a basis of Rd, so rd=0 and βd=0. The last run therefore ends with an exact breakdown after iteration d, and the claim gives AQd=QdTd. Since Qd is square with Qd⊤Qd=Id, also QdQd⊤=Id, and multiplying on the right by Qd⊤ gives A=QdTdQd⊤. By construction, Td is symmetric and tridiagonal, and its off-diagonal entry βj is zero at every boundary between runs.
3. Complete eigenbasis. Since TdZ=ZDZ⊤Z=ZD,
A(QdZ)=QdTdZ=QdZD
and
(QdZ)⊤(QdZ)=Z⊤Qd⊤QdZ=Z⊤Z=Id.
Therefore, V=QdZ is a complete orthonormal eigenbasis.
■□
Proposition 4 accepts any unit vector orthogonal to q1,…,qj as the new starting vector. One way to choose it is to draw a Gaussian vector g, independently of q1,…,qj, and to remove its components along q1,…,qj by applying Pj=Id−QjQj⊤. Since Pj is symmetric and Pj2=Id−2QjQj⊤+Qj(Qj⊤Qj)Qj⊤=Pj, the vector u=Pjg is Gaussian with covariance PjPj⊤=Pj:
The matrix Pj is the orthogonal projection onto range(Qj)⊥=ker(Qj⊤) since QjQj⊤ is the one onto range(Qj), so Qj⊤qj+1=0, that is, qj+1 is orthogonal to q1,…,qj. The division by ∥u∥ requires u=0, which holds with probability one when j<d. Indeed, u=0 means that g∈range(Qj), a subspace of dimension j<d, which has Gaussian probability zero, as in the proof of Corollary 3.
This continuation keeps all the vectors q1,…,qj already computed, and the new run appends its vectors to them, so that Qd contains the vectors of every run. A restart, in contrast, limits memory by discarding the basis and running the method again from a vector that combines a few Ritz vectors, as described after Algorithm 3.
4.12 Algorithm 2 [Lanczos with reorthogonalization and residual checks]
Proposition 3 and Proposition 4 recover eigenpairs exactly, but only at exact breakdowns and in exact arithmetic. The algorithm below is their finite version. It keeps the recurrence of Definition 2 and the Ritz pairs of Lemma 4, with three changes.
The basis is reorthogonalized. In floating-point arithmetic, the short recurrence gradually loses orthogonality, as measured in Section 4.3. Each new remainder is therefore projected once more against all the stored basis vectors, rj←rj−Qj(Qj⊤rj),
before βj is computed. Such See the discussion of full reorthogonalization in Bai et al., Templates for the Solution of Algebraic Eigenvalue Problems. The projection is carried out with classical or modified Gram–Schmidt, and it is repeated when the norm of the new vector drops below about 1/2 of its previous value, which rarely happens. is a standard part of practical Lanczos implementations (Bai et al., 2000).
Ritz pairs are accepted on a certificate. Instead of waiting for an exact breakdown, the algorithm tests the Ritz pairs as the basis grows, and accepts a pair as soon as its residual with A is below a tolerance, which by Lemma 2 places its Rayleigh quotient near an eigenvalue of A.
Breakdowns are detected with a tolerance. A computed βj is almost never exactly zero. A step is treated as a breakdown when βj falls below a small threshold, and for a complete decomposition the basis is then continued from a new starting vector, as in Proposition 4.
The algorithm uses three tolerances.
The tolerances εabs>0 and εrel≥0 play the same roles as in Algorithm 1. Step 5 accepts a candidate (λ,x) when ∥Ax−λx∥≤εabs+εrel∥Ax∥, and by Lemma 2 the accepted λ then lies within this threshold of an eigenvalue of A. With εrel=0, every eigenvalue is found to within the same distance εabs, and a positive εrel asks instead for about the same number of correct digits on every eigenvalue, as explained for Algorithm 1. One difference makes the relative rule safer here. Algorithm 2 does not deflate, so the residual of an accepted pair does not enter the test of any other pair, and the stall of Algorithm 1, where residuals inherited from earlier stages exceed the threshold, cannot occur. The floor remains: rounding errors keep computed residuals of the order of the unit roundoff times ∥A∥, so εabs must exceed that level. Figure 11 uses the constant threshold εabs=10−10∥A∥ and εrel=0, the same as for Algorithm 1, so that the two algorithms are compared on equal terms.
The breakdown tolerance τ>0 decides when βj counts as zero in step 3. Even when the Krylov subspace is invariant under A, rounding errors leave a computed βj of the order of the unit roundoff times ∥A∥, so τ must exceed that level. The runs of Figure 11 use τ=10−10 times the largest ∣αi∣ or βi computed so far, an estimate of ∥A∥ that costs nothing.
Choose the number k of requested eigenpairs (with k=d for a complete decomposition), a period p≥1 for the checks of steps 4 and 5, and a maximum number of iterations, at most d. Steps 1 and 2 are done once. Steps 3 to 6 form one iteration, which adds one vector to the basis. The iterations are numbered j=1,2,…, and step 6 either stops the algorithm or increases j by one and starts the next iteration at step 3. Figure 10 draws this execution.
Setup.
Fix the operator. Provide a routine computing x↦Ax, and keep everything that defines A fixed throughout the run, so that every product applies the same matrix. For a Hessian ∇θ2f(θ), this means the point θ and the function f, as noted in Section 2.1.
Initialize the basis. Choose a unit vector q1, for instance q1=g/∥g∥ with g∼N(0,Id), which by Corollary 3 has a component in every eigenspace of A with probability one. Set q0=0, β0=0, Q1=[q1] and j=1.
Iteration j≥1.
Extend the basis (one product). Compute zj=Aqj−βj−1qj−1, αj=qj⊤zj and rj=zj−αjqj, as in Definition 2. Reorthogonalize, rj←rj−Qj(Qj⊤rj), repeating the projection once when needed, and set βj=∥rj∥. If βj>τ, set qj+1=rj/βj. Otherwise, the run has broken down: set βj=0, and step 6 draws qj+1.
Solve the small eigenvalue problem. Do this step and step 5 only when j is a multiple of p or the run has just broken down, and otherwise go to step 6. Assemble Tm with m=j, and compute the eigenpairs Tmyi=ϑiyi with ∥yi∥=1 that correspond to the requested ones, for instance the min(j,k) whose eigenvalues ϑi have the largest absolute values, since Tm has only m=j eigenpairs. The candidates for step 5 are the pairs with a small residual estimate, βm∣(yi)m∣≤εabs+εrel∣ϑi∣, so there are at most min(j,k) of them. A pair with a larger estimate is not accurate yet: it is skipped at this iteration and examined again at the next check. By Lemma 4, this estimate equals the residual norm ∥AQmyi−ϑiQmyi∥ in exact arithmetic, and it costs no product, whereas step 5 costs one product per pair. After a breakdown, βm=0, so every pair becomes a candidate. If there are fewer than k candidates, which is always the case when j<k, go to step 6: the algorithm stops only when all k pairs pass step 5 at the same iteration, so step 5 would spend products in vain.
Check the k candidates with the original operator (one product each). For each candidate of step 4, form xi=Qmyi/∥Qmyi∥, compute λi=xi⊤Axi and εi=∥Axi−λixi∥, and accept (λi,xi) when εi≤εabs+εrel∥Axi∥.
Stop or continue. Apply the first of the following three cases that holds.
Success. If step 5 has run at this iteration and accepted all k candidates, stop and return them. Pairs accepted at an earlier iteration do not count, since step 4 recomputes all the Ritz pairs at each check, and the same eigenpair would be counted twice.
Failure. If j has reached the maximum number of iterations, stop and report nonconvergence.
Next iteration. Otherwise, obtain qj+1. If the run has not broken down in step 3, qj+1=rj/βj was computed there. If it has broken down, start a new run as in Proposition 4, from the new starting vector described after it: draw g∼N(0,Id) and set qj+1=Pjg/∥Pjg∥ with Pjg=g−Qj(Qj⊤g), keeping βj=0. In both cases, set Qj+1=[Qj,qj+1], increase j by one, and go to step 3.
■□Figure 10. The execution of Algorithm 2. Steps 1 and 2 run once. Steps 3 to 6 form iteration j, which adds one vector to the basis, and the blue arrow is the passage to iteration j+1. The blue-tinted boxes cost products with A: one in step 3 and k in step 5. Steps 4 and 5 are skipped when j is not a multiple of p and the run has not broken down, and step 5 is skipped when fewer than k Ritz pairs are candidates, since the algorithm stops only when all k pass step 5 at the same iteration. Step 6 then stops, or obtains qj+1: from step 3 while the run continues, or as a new starting vector after a breakdown, as in Proposition 4.
How the steps implement the theory
The table lists what each step computes and the result it implements. The paragraphs after it justify the three points where the finite computation departs from the exact theory.
Reorthogonalization restores what the theory assumes.Lemma 3, and through it Lemma 4 and Propositions 3 and 4, assume that the basis vectors are exactly orthonormal. The short recurrence keeps them orthogonal only up to errors that accumulate from step to step, and the projection of step 3 removes these errors at every step. In the example of Section 4.3, it keeps ∥Qm⊤Qm−Im∥ at about 10−15 up to m=200, against 10−6 without it.
The inexpensive estimate selects, the explicit residual certifies. In exact arithmetic, βm∣(yi)m∣ equals the residual norm ∥Axi−ϑixi∥ (Lemma 4). In floating-point arithmetic, the remaining loss of orthogonality and the corrections of the reorthogonalization perturb the identity AQm=QmTm+rmem⊤ on which this equality rests. Step 4 therefore uses the estimate only to decide which pairs are worth a product, and step 5 computes the residual with A itself, which Lemma 2 turns into a guarantee however xi was obtained. The normalization by ∥Qmyi∥ in step 5 matters for the same reason: when the computed basis is not exactly orthonormal, Qmyi need not have unit norm, and the inexpensive estimate can then be too optimistic (see the implementation discussion of the Lanczos algorithm in Bai et al., 2000).
A small βj is an approximate breakdown. By Lemma 3, ∥AQj−QjTj∥=∥rjej⊤∥=βj, so βj≤τ means that A maps range(Qj) into itself up to an error τ. It is not an exact proof of invariance or spectral completeness. Likewise, satisfying a residual tolerance certifies an approximate eigenpair, not that no more extreme eigenvalue has been missed.
Each iteration requires one matrix–vector product, in step 3. If CA denotes its cost, m iterations require m such products, with total cost mCA. Step 5 adds k products each time it runs, which happens only when all k pairs are candidates. Retaining Qm requires O(dm) storage. Full reorthogonalization costs O(dj) operations at iteration j, and therefore
∑j=1mO(dj)=O(dm2)
operations over the run, in addition to the eigenvalue problems of Tj that step 4 solves every p iterations and after each breakdown.
When A is available only through products, as in the examples of Section 2, the run never forms A. Nevertheless, a complete decomposition with m=d still requires O(d2) storage for the basis, and full reorthogonalization then costs O(d3) operations. The practical benefit is greatest when the required eigenpairs reach the desired accuracy with m≪d.
Figure 11 below compares the two algorithms on the network Hessian of Figure 7 and on the bipartite graph of Figure 3.
Figure 11.Algorithm 2 (green) against Algorithm 1 (blue), with the same constant threshold 10−10∥A∥ on the residuals for A. (a) The network Hessian of Figure 7: Hessian–vector products used until the j eigenpairs of largest absolute value are certified, medians over 20 starting vectors with bands containing the middle half of the runs. For Algorithm 1, the count runs over stages 1 to j. For Algorithm 2, with checks at every iteration (p=1), it is the number of Lanczos steps until the j Ritz pairs of largest absolute value pass the inexpensive test βm∣(yi)m∣≤10−10∥A∥ of Lemma 4, plus one product per pair for the explicit check of step 5, whose residuals were all below 5.6⋅10−10. The eight pairs cost 35 products with Algorithm 2 and about 1,100 with Algorithm 1. (b) The complete eigenbasis of the bipartite graph of Figure 3. The first Lanczos run breaks down at step 41, the number of distinct eigenvalues (Proposition 3), where its 41 Ritz pairs all pass the test. Each of the nine continuation runs of Proposition 4 starts in the eigenspace of 0, orthogonally to the basis built so far, and breaks down after one step. The 50 products, plus 50 explicit checks (green square), give all 50 eigenpairs with residuals below 1.2⋅10−11, against 15,296 products for Algorithm 1 with every stage taken to its limit, as in Figure 4.
5. The Arnoldi method
All the matrices of Section 2 are symmetric, and every proof so far uses the spectral theorem. Many matrices known through their products are not symmetric, and for them three facts of Section 1 fail. First, the eigenvalues of a real matrix can be complex, and then they come in conjugate pairs λ,λˉ. Second, eigenvectors of distinct eigenvalues need not be orthogonal. Third, the eigenvectors need not span the whole space: the matrix [0010] has the double eigenvalue 0, but its only eigenvectors are the multiples of [1,0]⊤.
The Krylov subspace of Definition 1 and the Gram–Schmidt process need no symmetry. The Arnoldi method (Arnoldi, 1951) runs Gram–Schmidt on Km(A,q1) for any square matrix, as the Lanczos method does for a symmetric one, and Saad (1980) developed it into a method for approximating a few eigenvalues of a large nonsymmetric matrix. This section follows the plan of Section 4: examples first, then the method worked by hand, its exact properties, the certificate that replaces Lemma 2, and an implementable procedure.
From here on, A∈Rd×d is any real square matrix, and Λ(A)⊂C is the set of its distinct eigenvalues. For a complex vector y, y∗ is its conjugate transpose and ∥y∥2=y∗y.
5.1 Nonsymmetric matrices known through their products
Convection–diffusion. The operator u↦−u′′+cu′ on [0,1], with u(0)=u(1)=0, describes a quantity that diffuses while a flow of speed c carries it, such as heat or a pollutant in a moving fluid. When c=0, it is the operator u↦−u′′ discretized in Section 2.2, and we discretize it in the same way. The grid points are ti=ih for i∈{0,…,d+1}, with step h=1/(d+1). A function u is represented by the vector x∈Rd of its values xi=u(ti) at the d interior points ti for i∈{1,…,d}, and the boundary conditions fix x0=xd+1=0. Replacing both derivatives at ti by If u has four continuous derivatives on [0,1], Taylor’s theorem gives u(ti±1)=u(ti±h)=u(ti)±hu′(ti)+2h2u′′(ti)±6h3u′′′(ti)+O(h4). Adding the two expansions cancels the odd powers of h and gives the difference for −u′′(ti) of Section 2.2, with error O(h2). Subtracting them cancels the even powers: u(ti+1)−u(ti−1)=2hu′(ti)+3h3u′′′(ti)+O(h4), so (xi+1−xi−1)/(2h)=u′(ti)+O(h2) when xj=u(tj). gives
−u′′(ti)+cu′(ti)≈h22xi−xi−1−xi+1+c2hxi+1−xi−1=h22xi−(1+γ)xi−1−(1−γ)xi+1,γ=2ch,
where the equality writes c/(2h) as γ/h2 and collects the terms in xi−1 and in xi+1. As in Section 2.2, the last expression is a linear function of x, so it defines a matrix A∈Rd×d: for every x∈Rd and i∈{1,…,d}, (Ax)i is this expression. When x holds the values of u, the vector Ax therefore holds approximations of the values of −u′′+cu′ at the interior points. The diagonal entries of A are 2/h2, its entries below the diagonal are Ai,i−1=−(1+γ)/h2, and its entries above the diagonal are Ai−1,i=−(1−γ)/h2. When c=0, γ=0 and A is the symmetric matrix of Section 2.2. When c=0, Ai,i−1=Ai−1,i, so A is not symmetric. For 0<γ<1, it is similar to a symmetric matrix. With s=(1+γ)/(1−γ) and S=diag(1,s,s2,…,sd−1), the matrix A=S−1AS is Since S is diagonal with Sii=si−1, the entries of A=S−1AS are Aik=Sii−1AikSkk=sk−iAik. We check that Aik=Aki in three cases. If ∣i−k∣≥2, then Aik=Aki=0, because A is tridiagonal, so Aik=Aki=0. If k=i, there is nothing to check, and Aii=Aii=2/h2. If ∣i−k∣=1, the two entries are, for i∈{2,…,d}, Ai,i−1Ai−1,i=sAi,i−1=−h21+γ1+γ1−γ=−h21−γ2,=sAi−1,i=−h21−γ1−γ1+γ=−h21−γ2. They are equal, so A is symmetric, and it is tridiagonal by the first case. The rescaling makes the weight 1+γ of the upstream neighbor and the weight 1−γ of the downstream one meet at their geometric mean., with 2/h2 on its diagonal and −1−γ2/h2 on both neighboring diagonals. The eigenvalues of A are therefore real. They are (2−21−γ2cos(kπh))/h2 for k=1,…,d, because the computation of Section 2.2, with −1−γ2/h2 in place of −1/h2 on the neighboring diagonals, shows that the vector with entries sin(kπti) is an eigenvector of A with this eigenvalue. The eigenvectors of A, however, are not orthogonal. Since A=SAS−1, an eigenvector v of A gives the eigenvector Sv of A for the same eigenvalue, so the eigenvector of A for the k-th eigenvalue has the entries si−1sin(kπti). For c>0, the fluid that carries the heat or pollutant moves at speed c in the direction of increasing t, toward t=1, and γ=ch/2>0 gives s>1. The factor si−1 then grows geometrically along the grid, from 1 at t1 to sd−1 at td, so every eigenvector has most of its norm near t=1. For d=40 and γ=0.2, the matrix of Figure 13, sd−1≈2.7⋅103, every unit eigenvector has at least 85% of its squared norm at the grid points with ti≥3/4, and two of them have an inner product of 0.98 in absolute value. The matrix V=SV, where the columns of V are the orthonormal eigenvectors of A, diagonalizes A with the condition number κ(V)=κ(S)=sd−1, the factor of Lemma 7. A product costs O(d) operations, as for −u′′.
Jacobians of dynamical systems. A system that evolves by x˙=F(x), for a map F:Rd→Rd, is at rest at a point x⋆ where F(x⋆)=0. A small perturbation y=x−x⋆ then evolves approximately by y˙=Jy, where J is the Jacobian of F at x⋆, with entries Jij=∂Fi(x⋆)/∂xj, because F(x⋆+y)≈Jy to first order. If J has a basis of eigenvectors v1,…,vd, with Jvj=λjvj, and the perturbation starts at y(0)=∑j=1dcjvj, then Each term solves the equation, dtd(cjeλjtvj)=λjcjeλjtvj=J(cjeλjtvj), so their sum does too, because the equation is linear, and at t=0 the sum equals y(0). A linear system with a constant matrix has exactly one solution for each starting point, so this sum is the solution. of y˙=Jy at time t is y(t)=∑j=1dcjeλjtvj. Since ∣eλjt∣=eRe(λj)t, the term of λj decays when the real part of λj is negative and grows when it is positive. The largest real part, α=maxjRe(λj), therefore decides whether small perturbations die out, since ∥y(t)∥≤(∑j∣cj∣∥vj∥)eαt. When α<0, every perturbation decays at least like eαt, and the equilibrium is stable. When α>0, a perturbation with cj=0 for some λj of real part α grows like eαt, and the equilibrium is unstable. The matrix J is not symmetric in general. A product with it is a directional derivative, Jx=dεdF(x⋆+εx)ε=0, which forward-mode automatic differentiation computes at the cost of a few evaluations of F. In machine learning, simultaneous gradient steps on the two players of a game, such as the generator and the discriminator of a generative adversarial network (GAN), follow such a vector field, and eigenvalues of its Jacobian with large imaginary parts make the iterates circle around an equilibrium rather than approach it (Mescheder et al., 2017). The simplest case is minθmaxϕθϕ, where θ˙=−ϕ and ϕ˙=θ, so that J=[01−10], whose eigenvalues ±i describe a pure rotation.
Random walks on directed graphs. On a directed graph with d nodes, let Wij≥0 be the weight of the edge from node i to node j, with Wij=0 when there is no such edge. An edge from i to j need not come with an edge from j to i, so W is not symmetric. A random walk steps from node i to node j with probability Pij=Wij/∑kWik, assuming that every node has an outgoing edge. The rows of the transition matrix P=diag(W1)−1W sum to one, P1=1. If (pk)i is the probability that the walker is at node i after k steps, then (pk+1)j=∑i(pk)iPij, that is, pk+1=P⊤pk. A stationary distribution, which one step of the walk leaves unchanged, is therefore an eigenvector of P⊤ with eigenvalue 1, P⊤π=π. It gives the long-run fraction of time that the walker spends at each node, and on the graph of hyperlinks between web pages, with random jumps added, it is the PageRank of Brin and Page (1998). The other eigenvalues of P⊤ set how Assume that P⊤ is diagonalizable and that 1 is a simple eigenvalue, and scale π so that its entries sum to one. Let λ2,…,λd be the other eigenvalues of P⊤, repeated according to multiplicity, with eigenvectors w2,…,wd such that π,w2,…,wd is a basis of Cd, which exists because P⊤ is diagonalizable. Write p0=c1π+∑j=2dcjwj. Since P1=1, we have 1⊤P⊤=1⊤, so 1⊤wj=1⊤P⊤wj=λj1⊤wj,j∈{2,…,d}, and λj=1, because 1 is simple, gives 1⊤wj=0. Summing the entries of p0, a distribution, therefore gives c1=1⊤p0=1. Since pk=(P⊤)kp0, and (P⊤)k multiplies each eigenvector by the k-th power of its eigenvalue, the term of π stays and, as c1=1, pk−π=(P⊤)kp0−π=c1(P⊤)kπ−π+j=2∑dcj(P⊤)kwj=c1π−π+j=2∑dcjλjkwj=j=2∑dcjλjkwj. The term of λj shrinks by the factor ∣λj∣ at each step, so the slowest terms, those of largest modulus, set the pace. Precisely, if every λj has modulus at most r<1, then ∥pk−π∥≤Crk,C=∑j=2d∣cj∣∥wj∥. In Theorem 1, the error shrinks by the factor ρ2/∣λ⋆∣ per step, where λ⋆ is the eigenvalue of largest modulus present in the starting vector and ρ2 is the largest modulus among the other eigenvalues present. Here that eigenvalue is 1, and r plays the role of ρ2. That theorem does not apply here, because it assumes a symmetric matrix and P⊤ is not one: the wj need not be orthogonal, and C can be much larger than ∥p0−π∥. Power iteration without symmetry, below, states what carries over., and complex ones make it Assume that P⊤ is diagonalizable, and expand p0 on a basis of its eigenvectors. Since P⊤ is real, conjugating P⊤w=λw gives P⊤wˉ=λˉwˉ, so the basis can be chosen so that the eigenvectors for λˉ are the conjugates of those for λ. Since p0 is real, conjugating its expansion gives an expansion of p0 in the same basis, and by uniqueness the coefficient of wˉ is the conjugate cˉ of the coefficient c of w. Write λ=ρeiω with ρ>0 and ω∈(0,π), and cw=a+ib with a,b∈Rd. After k steps, the two terms of the pair add up to the real vector cλkw+cˉλˉkwˉ=2Re(cλkw)=2ρk(cos(kω)a−sin(kω)b). Writing ai+ibi=mieiϕi with mi≥0, its entry i is 2ρkmicos(kω+ϕi), which changes sign every π/ω steps on average while its size shrinks like ρk. When this pair has the largest modulus among the eigenvalues other than 1 and c=0, it dominates pk−π for large k, so at every node with mi>0 the probability alternately exceeds and falls short of its limit. In Figure 12, ω≈119∘, so the pattern repeats about every three steps. A negative real eigenvalue is the case ω=π: its term changes sign at every step. A positive one keeps its sign and decays without oscillating. on the way, as Figure 12 below shows. A product with P⊤ costs one pass over the edges, (P⊤x)j=∑iPijxi over the edges i→j.
Figure 12. The random walk on a directed graph used in this section (Section 5.1). (a) The graph: 60 nodes in three groups of 20, teal, gold and pink, arranged in the cycle 1 → 2 → 3 → 1 that the large arrows show. Every node has five outgoing edges: four to distinct nodes of the next group, and one to a node anywhere else, other than itself and these four, all drawn uniformly at random. The 259 edges that end in the next group are drawn in the color of the group they leave, and the other 41 in black (dashed). (b) The transition matrix P, with rows and columns ordered by group and Pij=1/5 for every edge i→j: the edges into the next group fill three colored blocks, and the other edges are the scattered black entries. (c) The probability of each group under the distribution pk=(P⊤)kp0 of the walker, started uniformly on group 1. It circulates through the groups with period three and settles at the stationary values 0.335, 0.331 and 0.334, close to 1/3 (grey, dashed). Between steps 10 and 30, the largest gap to these values shrinks by a factor of 0.80 per step on average, measured as the 20th root of the ratio of the gaps at steps 30 and 10. This is close to the modulus 0.81 of the pair λ2,λˉ2 of panel (d), whose terms dominate the gap and shrink by that factor at each step. The two numbers differ slightly because these terms oscillate. The slider above the figure chooses a step k from 0 to 30: panel (c) marks it (black, dashed), and panel (a) draws at each node a disc, in the color of its group, whose area is proportional to the probability (pk)i that the walker is at that node after k steps. (d) The eigenvalues of P in the complex plane, all inside the unit circle (grey, dashed). The eigenvalue 1 has the eigenvector 1 for P, and the stationary distribution for P⊤. The pair λ2,λˉ2=0.81e±i119∘ produces the circulation of panel (c), a third of a turn per step damped by the factor 0.81. The other 57 eigenvalues have modulus below 0.45.
The figures of this section use a random walk on a graph with three groups of 20 nodes arranged in a cycle, shown in Figure 12 above. Every node has five outgoing edges, four to distinct nodes of the next group and one to another node anywhere, all drawn uniformly at random, and the walk follows one of them, chosen uniformly. The walk therefore tends to move from group 1 to group 2, then to group 3, and back to group 1. Started in group 1, its distribution circulates through the groups while it settles (panel c), and the eigenvalues of P⊤ account for both behaviors. The 60 eigenvalues of P⊤ are distinct, so it has a basis of eigenvectors. Let w be the eigenvector for λ2, the eigenvalue of largest modulus after 1 with a positive imaginary part, and let c be the coefficient of w when p0 is expanded in this basis. As shown in the item on random walks above, the coefficient of π is 1, the coefficient of wˉ, the eigenvector for λˉ2, is cˉ, and each step multiplies each term by its eigenvalue, so
pk=π+cλ2kw+cˉλˉ2kwˉ+rk=π+2Re(cλ2kw)+rk,
where rk is the sum of the terms of the other 57 eigenvalues. These have modulus below 0.45, so rk shrinks by at least that factor per step and soon becomes negligible against the middle term. Write λ2=ρeiω, with ρ≈0.81 and ω≈119∘ (panel d). The middle term equals 2ρkRe(eikωcw): it shrinks by the factor ρ at each step, which is the settling, and its phase advances by ω, about a third of a turn, which is the circulation. In other words, the pair λ2,λˉ2 is the cycle 1→2→3→1 of the graph, seen in the spectrum. The items below show this by following the probability of each group g∈{1,2,3}, the sum of the entries of pk over the nodes of g, where the group after group 3 is group 1.
The deviation of a group. Let sg=c∑i∈gwi, and write sg=∣sg∣eiϕg. After k steps, the probability of group gRecall the expansion pk=π+2Re(cλ2kw)+rk. Sum it over the nodes i of g. A sum of real parts is the real part of the sum, so ∑i∈g(pk)i−∑i∈gπi=2Re(λ2kc∑i∈gwi)+∑i∈g(rk)i=2Re(λ2ksg)+∑i∈g(rk)i. Since λ2ksg=ρkeikω∣sg∣eiϕg=ρk∣sg∣ei(kω+ϕg), its real part is ρk∣sg∣cos(kω+ϕg), which gives δg(k). by δg(k)=2ρk∣sg∣cos(kω+ϕg),
up to the sum of the entries of rk over g.
The phases. For this walk, ∣s1∣≈0.32, ∣s2∣≈∣s3∣≈0.33, and ϕ1,ϕ2,ϕ3≈−1∘,−120∘,118∘. Since 118∘=−242∘+360∘, each phase is about 120∘ below the previous one: ϕg+1≈ϕg−120∘, also from group 3 to group 1.
The deviation moves to the next group. For group g+1 at step k+1, the angle in the cosine is (k+1)ω+ϕg+1≈kω+ϕg+(119∘−120∘)≈kω+ϕg, the angle for group g at step k, and the factor in front is 2ρk+1∣sg+1∣≈ρ⋅2ρk∣sg∣. Therefore δg+1(k+1)≈0.81δg(k). For instance, the deviations −0.26, 0.54 and −0.28 of groups 1, 2 and 3 after one step become −0.20, 0.43 and −0.23 in groups 2, 3 and 1 after two steps.
Two numbers of one eigenvalue therefore describe panel (c) without running the walk. The angle ω≈119∘, about a third of a turn, is the part of the cycle that the deviation travels in one step, so it returns to the same group every 360∘/ω≈3 steps. The modulus ρ≈0.81 is the part of the deviation that survives each step. Power iteration cannot compute λ2, as the next paragraph shows, because λ2 and λˉ2 have the same modulus. The Arnoldi method can.
Power iteration without symmetry. The iteration vk+1=Avk/∥Avk∥ of Section 3 needs no symmetry. For A=P⊤ started from a distribution, it is the walk itself, rescaled: vk=pk/∥pk∥. Suppose that A is diagonalizable, that a single real eigenvalue λ1 has the largest modulus, as 1 does here, and that v0 has a nonzero component on its eigenvector. Let ∣λ2∣ be the largest modulus among the other eigenvalues. Then the iterates approach the eigenvector of λ1, with an error that shrinks like Expand v0=∑jcjwj in a basis of eigenvectors, with Awj=λjwj and c1=0. Since A and λ1 are real, w1 can be taken real, and then c1 is real. Then Akv0ek∥ek∥=λ1k(c1w1+ek),=j≥2∑cj(λ1λj)kwj,≤(j≥2∑∣cj∣∥wj∥)(∣λ1∣∣λ2∣)k. Write xk=c1w1+ek, a real vector since it equals Akv0/λ1k. Then sgn(λ1)kvk=sgn(λ1)k∥Akv0∥Akv0=sgn(λ1)k∣λ1∣k∥xk∥λ1kxk=∥xk∥xk, because λ1k/∣λ1∣k=sgn(λ1)k and sgn(λ1)2k=1. When λ1<0, the iterates therefore alternate in sign, as in Remark 2. With u=c1w1/∥c1w1∥ and the inequality x/∥x∥−a/∥a∥≤2∥x−a∥/∥a∥ for nonzero x and a, applied to x=xk and a=c1w1, this gives sgn(λ1)kvk−u=∥xk∥xk−∥c1w1∥c1w1≤∥c1w1∥2∥ek∥≤∥c1w1∥2∑j≥2∣cj∣∥wj∥(∣λ1∣∣λ2∣)k. This is the rate of Theorem 1, with ∣λ2∣ in the role of ρ2, but that theorem does not prove it here: its proof uses orthogonal eigenvectors. Without orthogonality, the constant 2∑j≥2∣cj∣∥wj∥/∥c1w1∥ can be much larger than tan(ϕ0).. Two features of Section 3 are lost. The Rayleigh quotient vk⊤Avk converges at that same rate, not at its square as in Theorem 1, because the eigenvectors of the other eigenvalues, which make up the error of the iterate, are Recall Awj=λjwj, v0=∑jcjwj with c1=0 and c1w1 real, α=∥c1w1∥ and u=c1w1/α, so that Au=λ1u and ∥u∥=1. Then sgn(λ1)kvkek=∥αu+ek∥αu+ek=u+dk,=j≥2∑cj(λ1λj)kwj, which defines dk. Here ∥ek∥ shrinks like (∣λ2∣/∣λ1∣)k, and ∥dk∥≤2∥ek∥/α by x/∥x∥−a/∥a∥≤2∥x−a∥/∥a∥. Since ∥u+dk∥=∥u∥=1 gives dk⊤u=−∥dk∥2/2, vk⊤Avk−λ1=(u+dk)⊤A(u+dk)−λ1=u⊤Adk+λ1dk⊤u+dk⊤Adk=u⊤Adk−2λ1∥dk∥2+dk⊤Adk. Hence, using ∣dk⊤Adk∣≤∥A∥∥dk∥2 and ∥dk∥≤2∥ek∥/α, vk⊤Avk−λ1−u⊤Adk=−2λ1∥dk∥2+dk⊤Adk≤(2∣λ1∣+∥A∥)∥dk∥2≤2(∣λ1∣+2∥A∥)α2∥ek∥2. So the error equals u⊤Adk up to a term of the order of ∥ek∥2. Let N=∥αu+ek∥. Then ∣N−α∣≤∥ek∥, because the triangle inequality ∥a+b∥≤∥a∥+∥b∥ gives both N−α≤∥αu∥+∥ek∥−α=∥ek∥ and α−N=∥(αu+ek)−ek∥−N≤∥ek∥. So N≥α/2 once ∥ek∥≤α/2, which holds for large enough k. More precisely, with C=∑j≥2∣cj∣∥wj∥, the triangle inequality gives ∥ek∥≤C(∣λ2∣/∣λ1∣)k, so ∥ek∥≤α/2, and hence N≥α/2, for every k≥ln(∣λ1∣/∣λ2∣)ln(2C/α). With u⊤Au=λ1, exactly, dku⊤AdkN−α=Nαu+ek−u=N(α−N)u+ek,=Nu⊤Aek−λ1(N−α),=N+αN2−α2=N+α2αu⊤ek+∥ek∥2. In the formula for u⊤Adk, replace N by α in the denominator and N−α by u⊤ek, which is its first-order part. This gives Sk=αu⊤Aek−λ1u⊤ek=αu⊤(A−λ1Id)ek. Then u⊤Adk−Sk=u⊤Aek(N1−α1)−λ1(NN−α−αu⊤ek). In the first term, N1−α1=αNα−N. In the second, substitute the formula for N−α, put both fractions over αN(N+α), and factor 2α2−N(N+α)=(α−N)(2α+N): NN−α−αu⊤ek=N(N+α)2αu⊤ek+∥ek∥2−αu⊤ek=αN(N+α)2α2u⊤ek+α∥ek∥2−N(N+α)u⊤ek=αN(N+α)(2α2−N(N+α))u⊤ek+N(N+α)∥ek∥2=αN(N+α)(α−N)(2α+N)u⊤ek+N(N+α)∥ek∥2. Both terms of u⊤Adk−Sk now carry a factor α−N or ∥ek∥2. Bound each factor with ∣u⊤Aek∣≤∥A∥∥ek∥, ∣α−N∣≤∥ek∥, ∣u⊤ek∣≤∥ek∥, (2α+N)/(N+α)≤2 and 1/(N+α)≤1/α, then use N≥α/2: u⊤Adk−Sk≤∣u⊤Aek∣αN∣α−N∣≤+∣λ1∣αN(N+α)∣α−N∣(2α+N)∣u⊤ek∣≤+∣λ1∣N(N+α)∥ek∥2≤∥A∥αN∥ek∥2+2∣λ1∣αN∥ek∥2+∣λ1∣αN∥ek∥2=(∥A∥+3∣λ1∣)αN∥ek∥2≤2(∥A∥+3∣λ1∣)α2∥ek∥2. Adding this bound to the bound on ∣vk⊤Avk−λ1−u⊤Adk∣ above, vk⊤Avk−λ1−Sk≤vk⊤Avk−λ1−u⊤Adk≤+u⊤Adk−Sk≤2(∣λ1∣+2∥A∥)α2∥ek∥2≤+2(∥A∥+3∣λ1∣)α2∥ek∥2=2(4∣λ1∣+3∥A∥)α2∥ek∥2, where Sk=α1u⊤(A−λ1Id)ek=α1u⊤(A−λ1Id)j≥2∑cj(λ1λj)kwj=α1j≥2∑cj(λ1λj)ku⊤(A−λ1Id)wj=α1j≥2∑cj(λ1λj)k(λj−λ1)u⊤wj. If A is symmetric, u⊤wj=(c1/α)w1⊤wj=0 for j≥2, so Sk=0 and the error is of the order of ∥ek∥2, the square of the rate, as in Theorem 1. Otherwise u⊤wj=0 in general, and Sk, hence the error, is of the order of (∣λ2∣/∣λ1∣)k, the rate itself. to the eigenvector of λ1. And the next eigenvalues can form a conjugate pair of equal modulus, which defeats power iteration in the way the opposite eigenvalues of Remark 5 do. Once the eigenvector of 1 is deflated, the pair dominates, and the iterates keep turning in the real plane that carries λ2 and λˉ2, without settling on a direction. Panel (b) of Figure 13 below shows both effects, and the Arnoldi method removes them.
5.2 How the Arnoldi method works
The idea The Arnoldi method runs the Gram–Schmidt process on the Krylov subspace
Km(A,q1)=span{q1,Aq1,…,Am−1q1}={p(A)q1:p is a polynomial of degree at most m−1}, one matrix–vector product per step, exactly as the Lanczos method does, but for a matrix that need not be symmetric.
Without symmetry, no coefficient vanishes for free: each step removes the components of the new vector along all the previous basis vectors. The coefficients fill a small upper Hessenberg matrix, and its eigenvalues, real or complex, approximate eigenvalues of A.
After m−1 products, the Krylov subspace contains every vector p(A)q1 with p of degree at most m−1, whatever A is. To search it, we again need an orthonormal basis q1,…,qm, the columns of Qm, and the coordinates of Aqj in that basis, which form the matrix Hm=Qm⊤AQm with entries (Hm)ij=qi⊤Aqj. The approximations are defined as in Lemma 4: the eigenpairs (ϑ,y) of Hm, lifted to x=Qmy, now possibly complex. What changes is why they are good approximations. For a symmetric matrix, Lemma 5 locates the extreme Ritz values through the extreme values of the Rayleigh quotient. Without symmetry, a Ritz value is still the Rayleigh quotient x∗Ax of its unit Ritz vector (Proposition 5), but these quotients are complex numbers, which have no order, and they fill a region that can be much larger than the spectrum (Remark 18). There are no extreme values to locate the eigenvalues. What remains is a condition: we look for a vector x=Qmy of the Krylov subspace Km(A,q1)=range(Qm) and a number ϑ whose residual Ax−ϑx is orthogonal to Km(A,q1). This choice has four consequences.
It is the eigenvalue problem of Hm. The columns of Qm are an orthonormal basis of Km(A,q1) (Lemma 6), so a vector is orthogonal to Km(A,q1) exactly when Qm⊤ maps it to 0, and Qm⊤Qm=Im. Hence
Qm⊤(Ax−ϑx)=Qm⊤AQmy−ϑQm⊤Qmy=Hmy−ϑy.
The condition Qm⊤(Ax−ϑx)=0 is therefore Hmy=ϑy, which has m solutions ϑ, counted with multiplicity, found by a small dense eigenvalue solver.
It is the most that can be asked of a vector of the subspace. The residual Ax−ϑx has d entries, while x has only the m free coordinates of y, so the residual cannot be made zero in general. The condition asks that its m coordinates in the basis q1,…,qm vanish, which gives m equations for the m entries of y.
No eigenpair of the subspace is lost. If an eigenvector of A lies in the subspace, its residual is zero, so it satisfies the condition, and the method returns it exactly.
What is left of the residual is visible. The residual Ax−ϑx has no component in the subspace, and Proposition 5 shows that it is a multiple of the next basis direction, with norm hm+1,m∣ym∣. As in the symmetric case, its size costs no extra product.
Gram–Schmidt on q1,Aq1,A2q1,… computes the basis and the matrix Hm together, and its first steps, compared with those of Section 4.3, show what symmetry was providing. The colors are those of Section 4.3: a purple coefficient is computed by an inner product at the current step, a green one is the length of the new direction, and a grey one vanishes.
Step 0: Gram–Schmidt on K1(A,q1)=span{q1}. There is nothing to do, as in Step 0 of Section 4.3.
Step 1: Gram–Schmidt on K2(A,q1). The first product gives Aq1, and Gram–Schmidt removes its component h11=q1⊤Aq1 along q1 and divides the remainder r1=Aq1−h11q1 by its length h21=∥r1∥. This is Step 1 of Section 4.3, where these numbers were called α1 and β1:
As in Section 4.3, write each basis vector as a polynomial in A applied to q1, qj+1=pj(A)q1, with p0=1. Step 1 gives q2=p1(A)q1 for the polynomial of degree one
Each later step solves its relation for the new basis vector and gives the next polynomial in the same way.
Step 2: Gram–Schmidt on K3(A,q1). The second product gives Aq2, which brings the same new direction as A2q1, by the argument of Step 2 of Section 4.3, which used no symmetry. Its component along q1 is where the two methods part. For a symmetric matrix, this component was (Aq1)⊤q2=β1, the length found in Step 1. In general, q1⊤Aq2=(A⊤q1)⊤q2, and nothing relates A⊤q1 to the Krylov subspace, so this component is a new number h12, computed by an inner product. Removing it and the component h22=q2⊤Aq2 along q2, then normalizing, gives
Aq2=h12q1+h22q2+h32q3.
Solving for q3, with q1=p0(A)q1 and q2=p1(A)q1, gives q3=p2(A)q1 for the polynomial of degree two defined by
h32p2(t)=tp1(t)−h12p0(t)−h22p1(t).
The coefficients found so far give H2=Q2⊤AQ2=[h11h21h12h22], which has the shape of T2 in Section 4.3, except that its entry h12 above the diagonal is a new number, not the length h21 found in Step 1.
Step 3: Gram–Schmidt on K4(A,q1). The third product gives Aq3. In Section 4.3, its component along q1 vanished because q1⊤Aq3=(Aq1)⊤q3=0. Here this component is (A⊤q1)⊤q3, nonzero in general, and all three components have to be computed and removed:
Aq3=h13q1+h23q2+h33q3+h43q4.
Solving for q4 in the same way gives q4=p3(A)q1 for the polynomial of degree three defined by
Some coefficients still vanish. The vector Aq1 lies in K2(A,q1)=span{q1,q2}, so it has no component along q3, and Aq2, which lies in K3(A,q1), has none along q4. These zeros use only the nesting of the Krylov subspaces, not symmetry. Writing the coefficients of Aqj in column j, as the identity (H3)ij=qi⊤Aqj requires, and doing the same for the Lanczos coefficients, gives
In T3, symmetry made each entry above the diagonal either a length already found (orange) or zero (grey). In H3, these entries are new numbers (purple). The zeros below the first subdiagonal remain in both matrices, because they come from the nesting of Krylov subspaces. A matrix whose entries below the first subdiagonal are zero is called upper Hessenberg.
Step 4: Gram–Schmidt on K5(A,q1). The fourth product gives Aq4. By the argument of Step 3, its components along q1, q2 and q3 are nonzero in general, and the one along q4 is a new number too. Since Aq4 lies in K5(A,q1)=span{q1,…,q5}, it has no component beyond q5:
Aq4=h14q1+h24q2+h34q3+h44q4+h54q5.
The number of inner products grows by one at each step: step j computes j of them, where the Lanczos method computes one. Solving the relation above for q5, as in Steps 1 to 3, gives q5=p4(A)q1 for the polynomial of degree four defined by
Each new polynomial involves all the previous ones, and not only the last two as in the three-term recurrence of Section 4.3. The coefficients found so far fill the upper Hessenberg matrix
whose grey zeros are the three entries below the first subdiagonal. The last coefficient, h54 on q5, lies outside K4(A,q1), as β4 did in Section 4.3.
The general step. The same reasoning holds at every step j∈{1,2,3,…}, as long as the previous steps have produced qj, which for j≥2 means that hj,j−1>0. The product Aqj lies in Kj+1(A,q1)=span{q1,…,qj+1}, because qj∈Kj(A,q1), so it has no component beyond qj+1. Its components hij=qi⊤Aqj along all of q1,…,qj, on the other hand, have to be computed and removed. What remains is orthogonal to q1,…,qj, and when its length hj+1,j is positive, dividing by it gives the next basis vector:
Step j of Definition 3 below is therefore the following Gram–Schmidt step.
qjAAqj−h1jq1−⋯−hjjqjrj÷hj+1,jqj+1.
Step j costs one product, j inner products and j vector updates, that is, O(dj) operations besides the product, and every basis vector has to be kept. After m steps, the method has used m products, O(dm2) further operations and O(dm) storage. This is the price of nonsymmetry. It equals the cost of the Lanczos method with full reorthogonalization (Algorithm 2), which subtracts all the components anyway, to protect orthogonality in floating-point arithmetic.
Eigenvalues and eigenvectors from the coefficients. Writing the relations Aqj=∑i=1j+1hijqi for j=1,…,m as the columns of Compare the columns. Column j of AQm is Aqj. Column j of QmHm is Qm times column j of Hm, whose nonzero entries are h1j,…,hj+1,j, in rows 1 to j+1. For j<m, it is therefore ∑i=1j+1hijqi=Aqj. The last column is the exception: Hm has no row m+1, so column m of QmHm stops at ∑i=1mhimqi and misses the term hm+1,mqm+1. The green matrixhm+1,mqm+1em⊤ adds it back: its columns are zero except the last, which is hm+1,mqm+1. gives
AQm=QmHm+hm+1,mqm+1em⊤,
the analogue of the Lanczos identity, with the same leak through the last basis vector. Every eigenpair (ϑ,y) of Hm, with y∈Cm and ∥y∥=1, gives a Ritz pair (ϑ,Qmy), as in Section 4.3. Multiplying the identity by y gives
so the residual of the Ritz pair is again the leak, of norm hm+1,m∣ym∣. The matrix Hm is real but, in general, not symmetric, so its eigenvalues can be complex, in conjugate pairs, which is how the method finds the complex eigenvalues of A.
The polynomial reading of Section 4.3 carries over, with one change. Its first two steps used no symmetry: the green term of Amq1 is now h21h32⋯hm+1,mqm+1, it has the smallest norm among the vectors π(A)q1 with π monic of degree m, and the Ritz values at step m are the roots of the minimizing polynomial πm (Proposition 5). What changes is how ∥π(A)q1∥ depends on the eigenvalues. Suppose that A=VDV−1 as in Lemma 7 below, with the eigenvalues λ1,…,λd on the diagonal of D, and write q1=Vc. Then π(D)c=V−1Vπ(D)V−1q1=V−1π(A)q1, and so Vπ(D)c=π(A)q1. Since ∥z∥=∥V−1Vz∥≤∥V−1∥∥Vz∥ for every z,
For a symmetric matrix, V can be chosen orthogonal, both factors equal 1, and ∥π(A)q1∥ is exactly the weighted sum over the eigenvalues of Section 4.3. Otherwise, the two bounds differ by the factor κ(V)=∥V∥∥V−1∥, which is large when the eigenvectors are nearly parallel. The terms π(λi)civi of π(A)q1=∑i=1dπ(λi)civi, where vi is column i of V, can then cancel, and a polynomial can make ∥π(A)q1∥ small without being small at the eigenvalues that q1 contains. The Ritz values can then stay far from every eigenvalue until m is large, and in panel (c) of Figure 13 below they leave the real axis although every eigenvalue is real.
5.3 Definition 3 [Arnoldi recurrence and exact breakdown]
Start with a unit vector q1. At iteration j≥1, compute
If hj+1,j>0, define qj+1=rj/hj+1,j and continue. If hj+1,j=0, stop the run, an exact breakdown as in Definition 2. For m≥1, Hm∈Rm×m is the matrix with entries (Hm)ij=hij for i≤j+1 and (Hm)ij=0 for i>j+1:
The coefficients can equally be computed one at a time: starting from w=Aqj, compute hij=qi⊤w and update w←w−hijqi for i=1,…,j. This is the modified Gram–Schmidt process. Each update leaves the inner products of w with the other basis vectors unchanged, since the basis vectors are orthogonal, so in exact arithmetic both orders give the same numbers. In floating-point arithmetic, the one-at-a-time order loses less orthogonality, and implementations repeat the subtraction when needed, as in Algorithm 2.
5.4 Lemma 6 [Orthogonality and Hessenberg projection]
Before and including the first exact breakdown, the Arnoldi vectors satisfy
qi⊤qj=δij,span{q1,…,qm}=Km(A,q1).
With Qm=[q1,…,qm] and the upper Hessenberg matrix Hm of Definition 3,
AQm=QmHm+rmem⊤,Qm⊤rm=0,Hm=Qm⊤AQm.
If A is symmetric, then Hm is the symmetric tridiagonal matrix Tm of Lemma 3, and the Arnoldi vectors are the Lanczos vectors.
Remark 17 [Where the zeros of the projected matrix come from]
The zeros of Tm have two sources. Those below the first subdiagonal come from the nesting of Krylov subspaces: Aqj lies in Kj+1(A,q1), so it has no component along qj+2,qj+3,…. Those above the first superdiagonal come from symmetry, which turns the component qi⊤Aqj into (Aqi)⊤qj. The Arnoldi method keeps the first source and loses the second. This is why Hm is upper Hessenberg, and why each step removes every previous component.
Proof
Orthogonality. We argue by induction, as for Lemma 3. The initial vector is normalized by assumption. Suppose that q1,…,qj are orthonormal. For each i≤j, by the definition of hij and since qi⊤ql=δil,
qi⊤rj=qi⊤Aqj−∑l=1jhljqi⊤ql=hij−hij=0.
Thus rj is orthogonal to q1,…,qj, and if hj+1,j>0, normalizing it extends the orthonormal family. Unlike the proof of Lemma 3, this one needs no symmetry, because every component is removed explicitly.
Krylov subspaces. We prove span{q1,…,qm}=Km(A,q1) by induction on m. It holds at m=1, since both sides equal span{q1}. Suppose that it holds at m=j and that hj+1,j>0. The vector qj+1 is a combination of Aqj and q1,…,qj, all of which lie in Kj+1(A,q1). Indeed, q1,…,qj belong to Kj(A,q1)⊆Kj+1(A,q1), and AqjSince qj∈Kj(A,q1), there are numbers c0,…,cj−1 with qj=∑k=0j−1ckAkq1, so Aqj=∑k=0j−1ckAk+1q1 is a combination of Aq1,…,Ajq1, which belong to Kj+1(A,q1) by Definition 1.. Hence span{q1,…,qj+1}⊆Kj+1(A,q1). The left-hand side has dimension j+1, because its vectors are orthonormal, so j+1≤dimKj+1(A,q1). We also have dimKj+1(A,q1)≤j+1, since Kj+1(A,q1) is generated by the j+1 vectors q1,Aq1,…,Ajq1. We thus have dimKj+1(A,q1)=dimspan{q1,…,qj+1}=j+1. The two spaces are therefore equal, since span{q1,…,qj+1} is contained in Kj+1(A,q1) and both have the same finite dimension.
Matrix identity. The definition of rj gives Aqj=∑i=1jhijqi+rj, where rj=hj+1,jqj+1 for j<m. These identities, for j=1,…,m, are the columns of AQm=QmHm+rmem⊤. Multiplying on the left by Qm⊤, and using Qm⊤Qm=Im and Qm⊤rm=0 from the orthogonality part, gives Hm=Qm⊤AQm.
The symmetric case. If A is symmetric, then so is Hm=Qm⊤AQm. Its entries below the first subdiagonal vanish, so by symmetry its entries above the first superdiagonal vanish too, and hj−1,j=hj,j−1. Step j therefore subtracts only hj,j−1qj−1 and hjjqj from Aqj, which is the Lanczos step of Definition 2 with αj=hjj and βj−1=hj,j−1. By induction on j, both methods produce the same vectors and the same coefficients, and Hm=Tm.
■□
5.5 Proposition 5 [Ritz pairs, the minimizing polynomial, and exact termination]
Suppose that the run has reached iteration m, that is, hj+1,j>0 for 1≤j<m.
Ritz pairs. Let Hmy=ϑy with y∈Cm and ∥y∥=1, and let x=Qmy. Then
∥x∥=1,x∗Ax=ϑ,Ax−ϑx=rmym,∥Ax−ϑx∥=hm+1,m∣ym∣.
The minimizing polynomial. The characteristic polynomial πm(t)=det(tIm−Hm) satisfies
πm(A)q1=h21h32⋯hm,m−1rm. Moreover, πm is the unique monic polynomial of degree m that minimizes ∥π(A)q1∥, and the Ritz values at step m are its roots.
Exact termination. Let μ be the monic polynomial of smallest degree s such that μ(A)q1=0, called the minimal polynomial of q1. A run reaches iteration s and breaks down there: hj+1,j>0 for j<s, and hs+1,s=0. Moreover, at termination,
det(tIs−Hs)=μ(t),AQs=QsHs,
so the Ritz values are exactly the roots of μ, and every one of them is an eigenvalue of A. If A has a basis of eigenvectors and q1=∑λ∈Λ(A)qλ with Aqλ=λqλ, then, as in Proposition 3,
μ(t)=λ∈Λ1∏(t−λ),Λ1={λ∈Λ(A):qλ=0},s=∣Λ1∣.
Remark 18 [Ritz values lie in the field of values]
Item 1 says that every Ritz value is the value of x∗Ax at a unit vector. It therefore lies in the field of values
W(A)={x∗Ax:x∈Cd,∥x∥=1},
which also contains every eigenvalue, since x∗Ax=λ for a unit eigenvector x. For a symmetric matrix, W(A) is the interval [λmin(A),λmax(A)], and the Ritz values stay between the extreme eigenvalues (Lemma 5). For a nonsymmetric matrix, W(A) is a convex region of the complex plane that can be much larger than the smallest convex set containing the eigenvalues, and the Ritz values can lie anywhere in it. In panel (c) of Figure 13 below, every eigenvalue is real, but the field of values is an oval around them, and the Ritz values spread over it.
Proof
Ritz pairs. As in the proof of Lemma 4, the orthonormality of the columns of Qm gives ∥x∥2=y∗Qm⊤Qmy=∥y∥2=1, and Lemma 6 gives x∗Ax=y∗Qm⊤AQmy=y∗Hmy=ϑ and
Because Hm is upper Hessenberg, hil=0 for i>l+1, so (Hmv)i=∑l=1mhilvl=∑l=i−1mhilvl for every vector v∈Cm. If the last nonzero entry of v is in position p, that is, vp=0 and vl=0 for all l>p, this sum is zero for i>p+1, since it then runs over l≥i−1≥p+1 only, and for i=p+1 it reduces to (Hmv)p+1=hp+1,pvp, which is nonzero because hp+1,p>0 for p<m when the run has reached iteration m. The entries in positions 1 to p can be any numbers and play no role. So, for p<m, multiplying by Hm moves the last nonzero entry down by one position and multiplies it by a subdiagonal entry. Applying this twice, Hm2v=Hm(Hmv) has its last nonzero entry in position p+2, equal to hp+2,p+1hp+1,pvp, if p+2≤m. By induction on k, for p+k≤m, the vector Hmkv has nonzero entries only in its first p+k positions, and its last nonzero entry is
(Hmkv)p+k=hp+k,p+k−1⋯hp+2,p+1hp+1,pvp.
Applying this to v=e1, for which p=1 and vp=1, shows that for k≤m−1 the vector Hmke1 has nonzero entries only in its first k+1 positions, and its entry in position k+1 is the product h21h32⋯hk+1,k of the subdiagonal entries met on the way.
We now prove equations (2) and (3) by induction on k. For k=0, equation (2) holds, since A0q1=q1=Qme1=QmHm0e1. Suppose, as induction hypothesis, that Akq1=QmHmke1 for some k≤m−1. Multiplying by A and using the identity AQm=QmHm+rmem⊤ of Lemma 6 gives
The last factor (Hmke1)m is given by the previous paragraph. For k≤m−2, i.e. m>k+1, we have (Hmke1)m=0 since position m lies beyond position k+1, and equation (4) is equation (2) for k+1. This completes the induction, so equation (2) holds for 0≤k≤m−1. For k=m−1, i.e. m=k+1, (Hmke1)m=h21h32⋯hm,m−1, and equation (4) is equation (3).
Now let π(t)=tm+∑k=0m−1cktk be a monic polynomial of degree m. Applying equation (3) to Amq1 and equation (2) to the lower powers gives
For π=πm=det(tIm−Hm), the Cayley–Hamilton theorem gives πm(Hm)=0, hence the identity of item 2:
πm(A)q1=h21h32⋯hm,m−1rm.
By Lemma 6, the vectors q1,…,qm span Km(A,q1) and Qm⊤rm=0, so this vector πm(A)q1 is orthogonal to Km(A,q1). Now let π be any monic polynomial of degree m. The leading terms tm cancel in π−πm, so π−πm=∑k=0m−1bktk for some numbers bk, and (π−πm)(A)q1=∑k=0m−1bkAkq1 lies in Km(A,q1). Hence π(A)q1=πm(A)q1+(π−πm)(A)q1 is a sum of two orthogonal vectors, and Pythagoras’ theorem gives
So πm minimizes ∥π(A)q1∥ over all monic polynomials π of degree m. It is the only minimizer. Indeed, if ∥π(A)q1∥=∥πm(A)q1∥, equation (5) gives (π−πm)(A)q1=0, that is, ∑k=0m−1bkAkq1=0. The m vectors q1,Aq1,…,Am−1q1 span Km(A,q1), which has dimension m, since by Lemma 6 it is spanned by the m orthonormal vectors q1,…,qm. They are therefore linearly independent, so b0=⋯=bm−1=0, that is, π=πm. Finally, the Ritz values are the eigenvalues of Hm, which are the roots of its characteristic polynomial πm.
Exact termination. For j<s, no monic polynomial of degree j annihilates q1, by the minimality of s, so πj(A)q1=0. If the run has reached iteration j<s, item 2 at step j gives πj(A)q1=h21h32⋯hj,j−1rj, where the factors h21,…,hj,j−1 are positive and ∥rj∥=hj+1,j by Definition 3. Taking norms,
0<∥πj(A)q1∥=h21h32⋯hj,j−1∥rj∥=h21h32⋯hj,j−1hj+1,j,
hence hj+1,j>0, and the run reaches iteration j+1. Starting from j=1, this shows that the run reaches iteration s. There, item 2 at step s says that πs is the unique monic polynomial of degree s that minimizes ∥π(A)q1∥. The polynomial μ is monic of degree s and gives ∥μ(A)q1∥=0, the smallest value a norm can take, so μ is a minimizer, and by uniqueness πs=μ. The same computation at step s gives h21⋯hs,s−1hs+1,s=∥πs(A)q1∥=∥μ(A)q1∥=0, which forces 0=hs+1,s=∥rs∥, since the other factors are positive. So rs=0, and the identity of Lemma 6 becomes AQs=QsHs. If Hsy=ϑy with y=0, then A(Qsy)=ϑQsy with Qsy=0, so every Ritz value at termination is an eigenvalue of A.
The case of a basis of eigenvectors. If q1=∑λqλ with Aqλ=λqλ, then p(A)q1=∑λ∈Λ1p(λ)qλ for every polynomial p. Nonzero eigenvectors of distinct eigenvalues are linearly independent, so p(A)q1=0 exactly when p(λ)=0 for every λ∈Λ1. The monic polynomial of smallest degree with this property is ∏λ∈Λ1(t−λ), of degree ∣Λ1∣.
■□
Item 2 also proves, for the Lanczos method, the statement of Section 4.3 that the Ritz values are the roots of the minimizing polynomial, since Hm=Tm for a symmetric matrix (Lemma 6). Without a basis of eigenvectors, μ can have repeated roots, and Hs then has them as repeated eigenvalues.
5.6 Lemma 7 [Residual certificate without symmetry]
Suppose that A is diagonalizable, that is, Cd has a basis of eigenvectors of A, and write A=VDV−1, where the columns of the invertible matrix V∈Cd×d are these eigenvectors and the diagonal matrix D holds the eigenvalues. Let κ(V)=∥V∥∥V−1∥, where ∥M∥=max∥z∥=1∥Mz∥.
Then, for every unit vector x∈Cd and every z∈C,
λ∈Λ(A)min∣z−λ∣≤κ(V)∥Ax−zx∥.
Proof
If z is an eigenvalue of A, the left-hand side is zero and there is nothing to prove. Otherwise, let δ=minλ∈Λ(A)∣λ−z∣>0. The diagonal entries of D−zId are the numbers λ−z, all nonzero, so this matrix is invertible, and so is A−zId=V(D−zId)V−1, with inverse V(D−zId)−1V−1. We have ∥(D−zId)−1∥=maxλ∈Λ(A)1/∣λ−z∣=1/δ, so
The theorem of Bauer and Fike (1960), stated for the norm ∥⋅∥ of Lemma 7, is about perturbations of the matrix A=VDV−1 of the lemma: for every matrix E∈Cd×d, every eigenvalue μ of the perturbed matrix A+E satisfies minλ∈Λ(A)∣μ−λ∣≤κ(V)∥E∥. This theorem and Lemma 7 imply each other.
From the theorem to Lemma 7. Let x be a unit vector and z a number, with residual r=Ax−zx∈Cd. The matrix E=−rx∗∈Cd×d makes z an eigenvalue of A+E, since
(A+E)x=Ax−r(x∗x)=Ax−r=zx,
and it has the norm of the residual, since
∥E∥=max∥w∥=1∥Ew∥=max∥w∥=1∥r(x∗w)∥=∥r∥max∥w∥=1∣x∗w∣=∥r∥,
where we use the fact that For the unit vector x and every unit vector w, the Cauchy–Schwarz inequality gives ∣x∗w∣≤∥x∥∥w∥=1, and w=x reaches this bound, since ∣x∗x∣=∥x∥2=1. Hence max∥w∥=1∣x∗w∣=1.. The theorem, applied to the eigenvalue μ=z of A+E, therefore gives minλ∈Λ(A)∣z−λ∣≤κ(V)∥E∥=κ(V)∥r∥, the bound of Lemma 7.
From Lemma 7 to the theorem. Given a matrix E and an eigenvalue μ of A+E with a unit eigenvector x, the residual of (μ,x) for A is Ax−μx=−Ex, of norm at most ∥E∥. Lemma 7 with z=μ then gives the bound of the theorem.
Remark 19 [A small residual far from the spectrum]
When A is symmetric, V can be chosen orthogonal, so κ(V)=1 and the lemma is Remark 11. Otherwise κ(V) can be large, and it cannot be dropped. Take
A=[0δ10],x=[10],
with a small δ>0. The eigenvalues are ±δ, the Rayleigh quotient is z=x⊤Ax=0, and the norm of the residual r=Ax−zx=[0,δ]⊤ is ∥r∥=δ, whereas the distance from z to the spectrum is δ. For δ=10−6, the norm of the residual is 10−6 and the distance 10−3, a thousand times larger. The eigenvectors [1,±δ]⊤ are nearly parallel, κ(V)=1/δ=103, and the lemma holds with equality. What a small residual r=Ax−zx always certifies is that (z,x) is an exact eigenpair of a nearby matrix: (A−rx∗)x=Ax−r=zx, and the perturbation has norm ∥rx∗∥=∥r∥. In the example, A−rx∗=[0010], whose eigenvalue 0 is double. A perturbation of size δ has moved the eigenvalues by δ, and the eigenvalues of a nonsymmetric matrix can in general move much farther than the size of a perturbation. The set of numbers z with ∥Ax−zx∥<ε for some unit vector x is called the ε-pseudospectrum of A, and for a nonsymmetric matrix it can extend far beyond the points within ε of an eigenvalue (Trefethen and Embree, 2005).
For the convection–diffusion matrix of Section 5.1, A=SAS−1 and A=UDU⊤ with U orthogonal, so V=SU diagonalizes A, with κ(V)=∥S∥∥S−1∥=sd−1. Panel (d) of Figure 13 tests Lemma 2 and Lemma 7 on its Ritz pairs.
Figure 13. The Arnoldi method on the two nonsymmetric matrices of Section 5.1. (a) The random walk of Figure 12, through products with P⊤: its eigenvalues (grey dots), its field of values (grey line, Remark 18), the unit circle (grey, dashed), and the Ritz values of one run, from a Gaussian starting vector, at the step m chosen with the slider above the figure (dark purple circles) and at step m−1 (light green circles), for m=1,…,20. Without JavaScript, the panel shows the steps m=5, 10 and 20 together, colored from light to dark. At m=5, the Ritz values 0.75 and −0.23±0.51i already point to the eigenvalue 1 and to the pair λ2,λˉ2. They are within 6⋅10−4 and 3⋅10−4 of them at m=10, and within 1.1⋅10−8 and 5.7⋅10−7 at m=20, while the other Ritz values lie among the remaining eigenvalues. (b) Distance to an eigenvalue against the number of products, as medians over 200 Gaussian starting vectors, with bands containing the middle half of the runs. The Rayleigh quotient of power iteration on P⊤ (orange) approaches 1 at the rate ∣λ2∣=0.81 per product (grey, dotted), with a median error of 1.2⋅10−3 after 20 products and 1.2⋅10−5 after 40, and power iteration never approaches λ2. After 20 products, the Ritz value nearest to 1 (green) is within 4.4⋅10−9 of it, and the one nearest to λ2 (green, dashed) within 4.4⋅10−7. Both reach rounding level, about 10−15, by 40 products. (c) The convection–diffusion matrix with d=40 and γ=0.2, multiplied by h2. Its eigenvalues are real, between 0.046 and 3.954 (black ticks), and its field of values (grey) is an oval that reaches ±0.40i. The Ritz values of one run, from a Gaussian starting vector, at m=10, 20 and 30 (light to dark) are all real at m=10, but 14 of them at m=20 and 24 at m=30 are not, with imaginary parts up to 0.32. (d) Distance from each Ritz value to the spectrum against the norm of its residual, ∥Ax−ϑx∥, for all the Ritz pairs of steps m=1,…,30 from 5 Gaussian starting vectors, 2,325 pairs per matrix. For the symmetric matrix of −u′′ (γ=0, blue), every point lies below the diagonal, as Lemma 2 requires, and the largest ratio of the distance to the norm of the residual is 0.49. For convection–diffusion (orange), 311 points lie above the diagonal, with distances up to 5.2 times the norm of the residual. Lemma 7 holds, with κ(V)=sd−1=2.7⋅103, but it overestimates every distance by a factor of at least 520. The inexpensive norm hm+1,m∣ym∣ of the residual, from Proposition 5, agrees with the norm computed explicitly to a relative 2⋅10−14.
Remark 20 [Matrices without a basis of eigenvectors]
The hypothesis of Lemma 7 is that A has a basis of eigenvectors, not that A is invertible: D may have zeros on its diagonal. Lemma 6 and Proposition 5 hold for every real square matrix, apart from the formula for μ at the end of item 3 of Proposition 5, but Lemma 7 fails without its hypothesis. Take the matrix J=[0010] of Remark 19, whose eigenvalue 0 is double but whose eigenvectors are only the multiples of [1,0]⊤. For z=ε>0 and the unit vector x=[1,ε]⊤/1+ε2,
Jx−εx=1+ε21([ε0]−[εε2])=1+ε21[0−ε2],
so the norm of the residual is less than ε2, while the distance from z to the spectrum {0} is ε, more than 1/ε times this norm. This ratio is unbounded as ε→0, so no constant can play the role of κ(V).
What remains is a bound with a root. Every square matrix has a Jordan form A=VJV−1, where V∈Cd×d is invertible and J is block diagonal, each block having one eigenvalue on its diagonal, ones just above the diagonal, and zeros elsewhere. Let p be the size of the largest block, and η=pκ(V)∥Ax−zx∥. Then, for every unit vector x∈Cd and every z∈C, the distance from z to the spectrum is The setting: A=VJV−1 with V invertible and J in Jordan form, p is the size of the largest Jordan block, ∥x∥=1, r=Ax−zx, η=pκ(V)∥r∥ and δ=minλ∈Λ(A)∣z−λ∣. If δ=0, there is nothing to prove. Otherwise J−zId is invertible and x=V(J−zId)−1V−1r, so 1=∥x∥≤κ(V)∥(J−zId)−1∥∥r∥. A Jordan block of size k≤p with eigenvalue λ becomes the block (λ−z)Ik+N of J−zId, where N has ones just above its diagonal, so that ∥Ni∥≤1 and Nk=0. Because Nk=0, multiplying by (λ−z)Ik+N shows that its inverse is the finite geometric series below, and since ∣λ−z∣≥δ and k≤p, ((λ−z)Ik+N)−1=i=0∑k−1(λ−z)i+1(−1)iNi≤i=0∑k−1∣λ−z∣i+11≤i=1∑pδi1. The norm of a block diagonal matrix is the largest norm of its blocks, so 1≤κ(V)∥r∥i=1∑pδi1. If δ≤1, each of the p terms is at most δ−p, so 1≤pκ(V)∥r∥δ−p=ηδ−p,that is,δ≤η1/p. If δ>1, each term is at most δ−1, so 1≤pκ(V)∥r∥δ−1=ηδ−1,that is,δ≤η. In both cases, δ≤max{η,η1/p}. of η when η<1:
λ∈Λ(A)min∣z−λ∣≤max{η,η1/p}.
When A has a basis of eigenvectors, all the blocks have size p=1, J=D, and this is Lemma 7. For the matrix J above, p=2 and V=I2, so a residual of norm less than ε2 gives η<2ε2 and the bound 2ε, against the actual distance ε: the square root cannot be avoided. It is the square root of Remark 19, where a perturbation of size δ moved the double eigenvalue 0 of J to ±δ.
5.7 Algorithm 3 [Arnoldi with reorthogonalization and residual checks]
Proposition 5 recovers eigenvalues exactly, but only at an exact breakdown and in exact arithmetic. The algorithm below is its finite version. It follows Algorithm 2, with three differences.
Every component is removed. Step 3 subtracts the components along all the stored basis vectors, as Definition 3 requires, and repeating this subtraction once is the reorthogonalization. A step therefore costs the same as a step of Algorithm 2.
The small eigenvalue problem is not symmetric. The matrix Hm needs a dense eigensolver for nonsymmetric matrices, and its eigenpairs can be complex.
The residual certifies less. By Remark 19, an accepted pair is an exact eigenpair of a matrix close to A, but by Lemma 7, its distance to the spectrum of A can exceed the residual by the factor κ(V), which the run does not compute.
The tolerances εabs>0, εrel≥0 and τ>0 play the same roles as in Algorithm 2. Step 5 accepts a candidate (λ,x) when its residual r=Ax−λx satisfies ∥r∥≤εabs+εrel∥Ax∥, and step 3 treats hj+1,j≤τ as a breakdown. Only the meaning of an accepted pair changes: (λ,x) is then an exact eigenpair of A−rx∗, a matrix within ∥r∥ of A, and when A=VDV−1, its λ lies within κ(V)∥r∥ of an eigenvalue of A.
Choose the number k of wanted eigenpairs, a period p≥1 for the checks of steps 4 and 5, and a maximum number of iterations, at most d. As in Algorithm 2, steps 1 and 2 are done once, and steps 3 to 6 form iteration j=1,2,…, which adds one vector to the basis. Figure 14 draws this execution.
Setup.
Fix the operator. Provide a routine computing x↦Ax, and keep everything that defines A fixed throughout the run.
Initialize the basis. Choose a unit vector q1, for instance q1=g/∥g∥ with g∼N(0,Id). Set Q1=[q1] and j=1.
Iteration j≥1.
Extend the basis (one product). Compute w=Aqj, the coefficients h=Qj⊤w, whose entries are h1j,…,hjj, and w←w−Qjh. Repeat the projection once, with c=Qj⊤w, w←w−Qjc and h←h+c, and store h in column j of H. Set hj+1,j=∥w∥. If hj+1,j>τ, set qj+1=w/hj+1,j. Otherwise, the run has broken down: set hj+1,j=0, and step 6 draws qj+1.
Solve the small eigenvalue problem. Do this step and step 5 only when j is a multiple of p or the run has just broken down, and otherwise go to step 6. With m=j, compute the eigenpairs Hmyi=ϑiyi, with yi∈Cm and ∥yi∥=1, by a dense eigensolver for nonsymmetric matrices, and select the min(j,k) wanted ones, for instance those of largest modulus or of largest real part. Nonreal Ritz values come in conjugate pairs, and one member of each pair suffices: since Hm is real, Hmyˉi=ϑˉiyˉi, and the conjugate pair has the same residual norm. The candidates for step 5 are the wanted pairs with a small residual estimate, hm+1,m∣(yi)m∣≤εabs+εrel∣ϑi∣. A pair with a larger estimate is not accurate yet: it is skipped at this iteration and examined again at the next check. By Proposition 5, this estimate equals the residual norm ∥AQmyi−ϑiQmyi∥ in exact arithmetic, and it costs no product, whereas step 5 costs one or two products per pair. After a breakdown, hm+1,m=0, so every wanted pair becomes a candidate. If there are fewer than k candidates, go to step 6, as in Algorithm 2.
Check the k candidates with the original operator. For each candidate of step 4, form xi=Qmyi/∥Qmyi∥, compute λi=xi∗Axi and εi=∥Axi−λixi∥, and accept (λi,xi) when εi≤εabs+εrel∥Axi∥. This costs one product for a real xi, and two for a complex one, since Axi=ARe(xi)+iAIm(xi).
Stop, restart or continue. Apply the first of the following four cases that holds.
Success. If step 5 has run at this iteration and accepted all k candidates, stop and return them.
Failure. If j has reached the maximum number of iterations, stop and report nonconvergence.
Restart. If storing Qj or orthogonalizing against it has become too expensive, restart as described below.
Next iteration. Otherwise, obtain qj+1. If the run has not broken down in step 3, qj+1=w/hj+1,j was computed there. If it has broken down, range(Qj) is invariant under A up to τ, and the other eigenvalues require a new starting vector: draw qj+1 orthogonal to q1,…,qj as after Proposition 4, keeping hj+1,j=0. Since step 3 orthogonalizes against all the stored vectors, the identity AQm=QmHm+hm+1,mqm+1em⊤ of Lemma 6 continues to hold, with a zero subdiagonal entry in Hm at each breakdown. In both cases, set Qj+1=[Qj,qj+1], increase j by one, and go to step 3.
■□Figure 14. The execution of Algorithm 3, drawn as Figure 10 is for Algorithm 2. Step 3 projects Aqj against all the stored vectors, and step 5 costs one product for a real Ritz vector and two for a complex one. Step 6 has a fourth case: when storing Qj or orthogonalizing against it has become too expensive, the method restarts (blue arrow on the right) from a new starting vector built from the wanted Ritz vectors, as described after the algorithm.
a new starting vector built from the wanted Ritz vectors
a polynomial filter, in the sense of Section 4.1, applied to the starting vector
The repeated projection. In exact arithmetic, the second projection of step 3 finds c=0, since w is already orthogonal to Qj, and h is unchanged. In floating-point arithmetic, the first projection leaves in w components along the stored vectors of the order of the unit roundoff u times ∥Aqj∥. Relative to ∥w∥=hj+1,j, this is of the order of u/f with f=hj+1,j/∥Aqj∥, by the cancellation explained in Section 4.3, and the second projection removes these components. Adding c to h keeps the identity Aqj=Qjh+hj+1,jqj+1 exact up to rounding, so that Hm is still the matrix of Lemma 6.
The residual certifies a nearby matrix. An accepted pair (λi,xi) is an exact eigenpair of a matrix within εi of A (Remark 19). Whether λi is within εi of an eigenvalue of A depends on the sensitivity of the eigenvalues, such as the factor κ(V) of Lemma 7, which the run does not compute. When the eigenvalue itself must be accurate, for instance when the sign of a real part decides the stability of an equilibrium, the residual should be read together with an estimate of that sensitivity, as discussed by Saad (2011) and by Trefethen and Embree (2005).
Restarts bound the cost. Since step j costs O(dj) operations and every basis vector is stored, m cannot grow indefinitely. A restart runs the method again from a new starting vector that combines the wanted Ritz vectors, so that the next Krylov subspace concentrates on them. The implicitly restarted Arnoldi method of Sorensen (1992) restarts without any new product. It applies to Hm a few shifted QR steps, with the unwanted Ritz values as shifts, which amounts to multiplying the starting vector by the polynomial ∏(t−ϑ) over those values, a polynomial filter in the sense of Section 4.1 that damps the unwanted eigencomponents. It then keeps a basis of k<m vectors and extends it again. This is the method of the ARPACK software (Lehoucq, Sorensen and Yang, 1998), on which the sparse eigenvalue routines of several scientific computing environments are built. For a symmetric matrix, the procedure computes, in exact arithmetic, the same vectors and coefficients as the Lanczos method (Lemma 6), at the cost of Algorithm 2.
Lars Mescheder, Sebastian Nowozin, and Andreas Geiger, “The Numerics of GANs”, Advances in Neural Information Processing Systems 30, 2017.
Ch. Müntz, “Solution directe de l’équation séculaire et de quelques problèmes analogues transcendants”, Comptes Rendus de l’Académie des Sciences, Paris, 1913.