Application to the eigenproblem

The Lanczos algorithm is most often brought up in the context of finding the eigenvalues and eigenvectors of a matrix, but whereas an ordinary diagonalization of a matrix would make eigenvectors and eigenvalues apparent from inspection, the same is not true for the tridiagonalization performed by the Lanczos algorithm; nontrivial additional steps are needed to compute even a single eigenvalue or eigenvector. Nonetheless, applying the Lanczos algorithm is often a significant step forward in computing the eigendecomposition. If ${\displaystyle \lambda }$ is an eigenvalue of ${\displaystyle A}$, and if ${\displaystyle Tx=\lambda x}$ (${\displaystyle x}$ is an eigenvector of ${\displaystyle T}$) then ${\displaystyle y=Vx}$ is the corresponding eigenvector of ${\displaystyle A}$ (since ${\displaystyle Ay=AVx=VTV^{*}Vx=VTIx=VTx=V(\lambda x)=\lambda Vx=\lambda y}$). Thus the Lanczos algorithm transforms the eigendecomposition problem for ${\displaystyle A}$ into the eigendecomposition problem for ${\displaystyle T}$.

1. For tridiagonal matrices, there exist a number of specialised algorithms, often with better computational complexity than general-purpose algorithms. For example, if ${\displaystyle T}$ is an ${\displaystyle m\times m}$ tridiagonal symmetric matrix then:
   • The continuant recursion allows computing the characteristic polynomial in ${\displaystyle O(m^{2})}$ operations, and evaluating it at a point in ${\displaystyle O(m)}$ operations.
   • The divide-and-conquer eigenvalue algorithm can be used to compute the entire eigendecomposition of ${\displaystyle T}$ in ${\displaystyle O(m^{2})}$ operations.
   • The Fast Multipole Method^[8] can compute all eigenvalues in just ${\displaystyle O(m\log m)}$ operations.
2. Some general eigendecomposition algorithms, notably the QR algorithm, are known to converge faster for tridiagonal matrices than for general matrices. Asymptotic complexity of tridiagonal QR is ${\displaystyle O(m^{2})}$ just as for the divide-and-conquer algorithm (though the constant factor may be different); since the eigenvectors together have ${\displaystyle m^{2}}$ elements, this is asymptotically optimal.
3. Even algorithms whose convergence rates are unaffected by unitary transformations, such as the power method and inverse iteration, may enjoy low-level performance benefits from being applied to the tridiagonal matrix ${\displaystyle T}$ rather than the original matrix ${\displaystyle A}$. Since ${\displaystyle T}$ is very sparse with all nonzero elements in highly predictable positions, it permits compact storage with excellent performance vis-à-vis caching. Likewise, ${\displaystyle T}$ is a real matrix with all eigenvectors and eigenvalues real, whereas ${\displaystyle A}$ in general may have complex elements and eigenvectors, so real arithmetic is sufficient for finding the eigenvectors and eigenvalues of ${\displaystyle T}$.
4. If ${\displaystyle n}$ is very large, then reducing ${\displaystyle m}$ so that ${\displaystyle T}$ is of a manageable size will still allow finding the more extreme eigenvalues and eigenvectors of ${\displaystyle A}$; in the ${\displaystyle m\ll n}$ region, the Lanczos algorithm can be viewed as a lossy compression scheme for Hermitian matrices, that emphasises preserving the extreme eigenvalues.

The combination of good performance for sparse matrices and the ability to compute several (without computing all) eigenvalues are the main reasons for choosing to use the Lanczos algorithm.

Application to tridiagonalization

Though the eigenproblem is often the motivation for applying the Lanczos algorithm, the operation the algorithm primarily performs is tridiagonalization of a matrix, for which numerically stable Householder transformations have been favoured since the 1950s. During the 1960s the Lanczos algorithm was disregarded. Interest in it was rejuvenated by the Kaniel–Paige convergence theory and the development of methods to prevent numerical instability, but the Lanczos algorithm remains the alternative algorithm that one tries only if Householder is not satisfactory.^[9]

Aspects in which the two algorithms differ include:

• Lanczos takes advantage of ${\displaystyle A}$ being a sparse matrix, whereas Householder does not, and will generate fill-in.
• Lanczos works throughout with the original matrix ${\displaystyle A}$ (and has no problem with it being known only implicitly), whereas raw Householder wants to modify the matrix during the computation (although that can be avoided).
• Each iteration of the Lanczos algorithm produces another column of the final transformation matrix ${\displaystyle V}$, whereas an iteration of Householder produces another factor in a unitary factorisation ${\displaystyle Q_{1}Q_{2}\dots Q_{n}}$ of ${\displaystyle V}$. Each factor is however determined by a single vector, so the storage requirements are the same for both algorithms, and ${\displaystyle V=Q_{1}Q_{2}\dots Q_{n}}$ can be computed in ${\displaystyle O(n^{3})}$ time.
• Householder is numerically stable, whereas raw Lanczos is not.
• Lanczos is highly parallel, with only ${\displaystyle O(n)}$ points of synchronisation (the computations of ${\displaystyle \alpha _{j}}$ and ${\displaystyle \beta _{j}}$). Householder is less parallel, having a sequence of ${\displaystyle O(n^{2})}$ scalar quantities computed that each depend on the previous quantity in the sequence.

A more provident power method

The power method for finding the eigenvalue of largest magnitude and a corresponding eigenvector of a matrix ${\displaystyle A}$ is roughly

1. Pick a random vector ${\displaystyle u_{1}\neq 0}$.
2. For ${\displaystyle j\geqslant 1}$ (until the direction of ${\displaystyle u_{j}}$ has converged) do:
   1. Let ${\displaystyle u_{j+1}'=Au_{j}.}$
   2. Let ${\displaystyle u_{j+1}=u_{j+1}'/\|u_{j+1}'\|.}$

• In the large ${\displaystyle j}$ limit, ${\displaystyle u_{j}}$ approaches the normed eigenvector corresponding to the largest magnitude eigenvalue.

A critique that can be raised against this method is that it is wasteful: it spends a lot of work (the matrix–vector products in step 2.1) extracting information from the matrix ${\displaystyle A}$, but pays attention only to the very last result; implementations typically use the same variable for all the vectors ${\displaystyle u_{j}}$, having each new iteration overwrite the results from the previous one. It may be desirable to instead keep all the intermediate results and organise the data.

One piece of information that trivially is available from the vectors ${\displaystyle u_{j}}$ is a chain of Krylov subspaces. One way of stating that without introducing sets into the algorithm is to claim that it computes

a subset ${\displaystyle \{v_{j}\}_{j=1}^{m}}$ of a basis of ${\displaystyle \mathbb {C} ^{n}}$ such that ${\displaystyle Ax\in \operatorname {span} (v_{1},\dotsc ,v_{j+1})}$ for every ${\displaystyle x\in \operatorname {span} (v_{1},\dotsc ,v_{j})}$ and all ${\displaystyle 1\leqslant j<m;}$

this is trivially satisfied by ${\displaystyle v_{j}=u_{j}}$ as long as ${\displaystyle u_{j}}$ is linearly independent of ${\displaystyle u_{1},\dotsc ,u_{j-1}}$ (and in the case that there is such a dependence then one may continue the sequence by picking as ${\displaystyle v_{j}}$ an arbitrary vector linearly independent of ${\displaystyle u_{1},\dotsc ,u_{j-1}}$). A basis containing the ${\displaystyle u_{j}}$ vectors is however likely to be numerically ill-conditioned, since this sequence of vectors is by design meant to converge to an eigenvector of ${\displaystyle A}$. To avoid that, one can combine the power iteration with a Gram–Schmidt process, to instead produce an orthonormal basis of these Krylov subspaces.

1. Pick a random vector ${\displaystyle u_{1}}$ of Euclidean norm ${\displaystyle 1}$. Let ${\displaystyle v_{1}=u_{1}}$.
2. For ${\displaystyle j=1,\dotsc ,m-1}$ do:
   1. Let ${\displaystyle u_{j+1}'=Au_{j}}$.
   2. For all ${\displaystyle k=1,\dotsc ,j}$ let ${\displaystyle g_{k,j}=v_{k}^{*}u_{j+1}'}$. (These are the coordinates of ${\displaystyle Au_{j}=u_{j+1}'}$ with respect to the basis vectors ${\displaystyle v_{1},\dotsc ,v_{j}}$.)
   3. Let ${\displaystyle w_{j+1}=u_{j+1}'-\sum _{k=1}^{j}g_{k,j}v_{k}}$. (Cancel the component of ${\displaystyle u_{j+1}'}$ that is in ${\displaystyle \operatorname {span} (v_{1},\dotsc ,v_{j})}$.)
   4. If ${\displaystyle w_{j+1}\neq 0}$ then let ${\displaystyle u_{j+1}=u_{j+1}'/\|u_{j+1}'\|}$ and ${\displaystyle v_{j+1}=w_{j+1}/\|w_{j+1}\|}$,

      otherwise pick as ${\displaystyle u_{j+1}=v_{j+1}}$ an arbitrary vector of Euclidean norm ${\displaystyle 1}$ that is orthogonal to all of ${\displaystyle v_{1},\dotsc ,v_{j}}$.

The relation between the power iteration vectors ${\displaystyle u_{j}}$ and the orthogonal vectors ${\displaystyle v_{j}}$ is that

${\displaystyle Au_{j}=\|u_{j+1}'\|u_{j+1}=u_{j+1}'=w_{j+1}+\sum _{k=1}^{j}g_{k,j}v_{k}=\|w_{j+1}\|v_{j+1}+\sum _{k=1}^{j}g_{k,j}v_{k}}$.

Here it may be observed that we do not actually need the ${\displaystyle u_{j}}$ vectors to compute these ${\displaystyle v_{j}}$, because ${\displaystyle u_{j}-v_{j}\in \operatorname {span} (v_{1},\dotsc ,v_{j-1})}$ and therefore the difference between ${\displaystyle u_{j+1}'=Au_{j}}$ and ${\displaystyle w_{j+1}'=Av_{j}}$ is in ${\displaystyle \operatorname {span} (v_{1},\dotsc ,v_{j})}$, which is cancelled out by the orthogonalisation process. Thus the same basis for the chain of Krylov subspaces is computed by

1. Pick a random vector ${\displaystyle v_{1}}$ of Euclidean norm ${\displaystyle 1}$.
2. For ${\displaystyle j=1,\dotsc ,m-1}$ do:
   1. Let ${\displaystyle w_{j+1}'=Av_{j}}$.
   2. For all ${\displaystyle k=1,\dotsc ,j}$ let ${\displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'}$.
   3. Let ${\displaystyle w_{j+1}=w_{j+1}'-\sum _{k=1}^{j}h_{k,j}v_{k}}$.
   4. Let ${\displaystyle h_{j+1,j}=\|w_{j+1}\|}$.
   5. If ${\displaystyle h_{j+1,j}\neq 0}$ then let ${\displaystyle v_{j+1}=w_{j+1}/h_{j+1,j}}$,

      otherwise pick as ${\displaystyle v_{j+1}}$ an arbitrary vector of Euclidean norm ${\displaystyle 1}$ that is orthogonal to all of ${\displaystyle v_{1},\dotsc ,v_{j}}$.

A priori the coefficients ${\displaystyle h_{k,j}}$ satisfy

${\displaystyle Av_{j}=\sum _{k=1}^{j+1}h_{k,j}v_{k}}$ for all ${\displaystyle j<m}$;

the definition ${\displaystyle h_{j+1,j}=\|w_{j+1}\|}$ may seem a bit odd, but fits the general pattern ${\displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'}$ since

${\displaystyle v_{j+1}^{*}w_{j+1}'=v_{j+1}^{*}w_{j+1}=\|w_{j+1}\|v_{j+1}^{*}v_{j+1}=\|w_{j+1}\|.}$

Because the power iteration vectors ${\displaystyle u_{j}}$ that were eliminated from this recursion satisfy ${\displaystyle u_{j}\in \operatorname {span} (v_{1},\ldots ,v_{j}),}$ the vectors ${\displaystyle \{v_{j}\}_{j=1}^{m}}$ and coefficients ${\displaystyle h_{k,j}}$ contain enough information from ${\displaystyle A}$ that all of ${\displaystyle u_{1},\ldots ,u_{m}}$ can be computed, so nothing was lost by switching vectors. (Indeed, it turns out that the data collected here give significantly better approximations of the largest eigenvalue than one gets from an equal number of iterations in the power method, although that is not necessarily obvious at this point.)

This last procedure is the Arnoldi iteration. The Lanczos algorithm then arises as the simplification one gets from eliminating calculation steps that turn out to be trivial when ${\displaystyle A}$ is Hermitian—in particular most of the ${\displaystyle h_{k,j}}$ coefficients turn out to be zero.

Elementarily, if ${\displaystyle A}$ is Hermitian then

${\displaystyle h_{k,j}=v_{k}^{*}w_{j+1}'=v_{k}^{*}Av_{j}=v_{k}^{*}A^{*}v_{j}=(Av_{k})^{*}v_{j}.}$

For ${\displaystyle k<j-1}$ we know that ${\displaystyle Av_{k}\in \operatorname {span} (v_{1},\ldots ,v_{j-1})}$, and since ${\displaystyle v_{j}}$ by construction is orthogonal to this subspace, this inner product must be zero. (This is essentially also the reason why sequences of orthogonal polynomials can always be given a three-term recurrence relation.) For ${\displaystyle k=j-1}$ one gets

${\displaystyle h_{j-1,j}=(Av_{j-1})^{*}v_{j}={\overline {v_{j}^{*}Av_{j-1}}}={\overline {h_{j,j-1}}}=h_{j,j-1}}$

since the latter is real on account of being the norm of a vector. For ${\displaystyle k=j}$ one gets

${\displaystyle h_{j,j}=(Av_{j})^{*}v_{j}={\overline {v_{j}^{*}Av_{j}}}={\overline {h_{j,j}}},}$

meaning this is real too.

More abstractly, if ${\displaystyle V}$ is the matrix with columns ${\displaystyle v_{1},\ldots ,v_{m}}$ then the numbers ${\displaystyle h_{k,j}}$ can be identified as elements of the matrix ${\displaystyle H=V^{*}AV}$, and ${\displaystyle h_{k,j}=0}$ for ${\displaystyle k>j+1;}$ the matrix ${\displaystyle H}$ is upper Hessenberg. Since

${\displaystyle H^{*}=\left(V^{*}AV\right)^{*}=V^{*}A^{*}V=V^{*}AV=H}$

the matrix ${\displaystyle H}$ is Hermitian. This implies that ${\displaystyle H}$ is also lower Hessenberg, so it must in fact be tridiagional. Being Hermitian, its main diagonal is real, and since its first subdiagonal is real by construction, the same is true for its first superdiagonal. Therefore, ${\displaystyle H}$ is a real, symmetric matrix—the matrix ${\displaystyle T}$ of the Lanczos algorithm specification.

Simultaneous approximation of extreme eigenvalues

One way of characterising the eigenvectors of a Hermitian matrix ${\displaystyle A}$ is as stationary points of the Rayleigh quotient

${\displaystyle r(x)={\frac {x^{*}Ax}{x^{*}x}},\qquad x\in \mathbb {C} ^{n}.}$

In particular, the largest eigenvalue ${\displaystyle \lambda _{\max }}$ is the global maximum of ${\displaystyle r}$ and the smallest eigenvalue ${\displaystyle \lambda _{\min }}$ is the global minimum of ${\displaystyle r}$.

Within a low-dimensional subspace ${\displaystyle {\mathcal {L}}}$ of ${\displaystyle \mathbb {C} ^{n}}$ it can be feasible to locate the maximum ${\displaystyle x}$ and minimum ${\displaystyle y}$ of ${\displaystyle r}$. Repeating that for an increasing chain ${\displaystyle {\mathcal {L}}_{1}\subset {\mathcal {L}}_{2}\subset \cdots }$ produces two sequences of vectors: ${\displaystyle x_{1},x_{2},\ldots }$ and ${\displaystyle y_{1},y_{2},\dotsc }$ such that ${\displaystyle x_{j},y_{j}\in {\mathcal {L}}_{j}}$ and

${\displaystyle {\begin{aligned}r(x_{1})&\leqslant r(x_{2})\leqslant \cdots \leqslant \lambda _{\max }\\r(y_{1})&\geqslant r(y_{2})\geqslant \cdots \geqslant \lambda _{\min }\end{aligned}}}$

The question then arises how to choose the subspaces so that these sequences converge at optimal rate.

From ${\displaystyle x_{j}}$, the optimal direction in which to seek larger values of ${\displaystyle r}$ is that of the gradient ${\displaystyle \nabla r(x_{j})}$, and likewise from ${\displaystyle y_{j}}$ the optimal direction in which to seek smaller values of ${\displaystyle r}$ is that of the negative gradient ${\displaystyle -\nabla r(y_{j})}$. In general

${\displaystyle \nabla r(x)={\frac {2}{x^{*}x}}(Ax-r(x)x),}$

so the directions of interest are easy enough to compute in matrix arithmetic, but if one wishes to improve on both ${\displaystyle x_{j}}$ and ${\displaystyle y_{j}}$ then there are two new directions to take into account: ${\displaystyle Ax_{j}}$ and ${\displaystyle Ay_{j};}$ since ${\displaystyle x_{j}}$ and ${\displaystyle y_{j}}$ can be linearly independent vectors (indeed, are close to orthogonal), one cannot in general expect ${\displaystyle Ax_{j}}$ and ${\displaystyle Ay_{j}}$ to be parallel. It is not necessary to increase the dimension of ${\displaystyle {\mathcal {L}}_{j}}$ by ${\displaystyle 2}$ on every step if ${\displaystyle \{{\mathcal {L}}_{j}\}_{j=1}^{m}}$ are taken to be Krylov subspaces, because then ${\displaystyle Az\in {\mathcal {L}}_{j+1}}$ for all ${\displaystyle z\in {\mathcal {L}}_{j},}$ thus in particular for both ${\displaystyle z=x_{j}}$ and ${\displaystyle z=y_{j}}$.

In other words, we can start with some arbitrary initial vector ${\displaystyle x_{1}=y_{1},}$ construct the vector spaces

${\displaystyle {\mathcal {L}}_{j}=\operatorname {span} (x_{1},Ax_{1},\ldots ,A^{j-1}x_{1})}$

and then seek ${\displaystyle x_{j},y_{j}\in {\mathcal {L}}_{j}}$ such that

${\displaystyle r(x_{j})=\max _{z\in {\mathcal {L}}_{j}}r(z)\qquad {\text{and}}\qquad r(y_{j})=\min _{z\in {\mathcal {L}}_{j}}r(z).}$

Since the ${\displaystyle j}$th power method iterate ${\displaystyle u_{j}}$ belongs to ${\displaystyle {\mathcal {L}}_{j},}$ it follows that an iteration to produce the ${\displaystyle x_{j}}$ and ${\displaystyle y_{j}}$ cannot converge slower than that of the power method, and will achieve more by approximating both eigenvalue extremes. For the subproblem of optimising ${\displaystyle r}$ on some ${\displaystyle {\mathcal {L}}_{j}}$, it is convenient to have an orthonormal basis ${\displaystyle \{v_{1},\ldots ,v_{j}\}}$ for this vector space. Thus we are again led to the problem of iteratively computing such a basis for the sequence of Krylov subspaces.

Kaniel–Paige convergence theory

After ${\displaystyle m}$ iteration steps of the Lanczos algorithm, ${\displaystyle T}$ is an ${\displaystyle m\times m}$ real symmetric matrix, that similarly to the above has ${\displaystyle m}$ eigenvalues ${\displaystyle \theta _{1}\geqslant \theta _{2}\geqslant \dots \geqslant \theta _{m}.}$ By convergence is primarily understood the convergence of ${\displaystyle \theta _{1}}$ to ${\displaystyle \lambda _{1}}$ (and the symmetrical convergence of ${\displaystyle \theta _{m}}$ to ${\displaystyle \lambda _{n}}$) as ${\displaystyle m}$ grows, and secondarily the convergence of some range ${\displaystyle \theta _{1},\ldots ,\theta _{k}}$ of eigenvalues of ${\displaystyle T}$ to their counterparts ${\displaystyle \lambda _{1},\ldots ,\lambda _{k}}$ of ${\displaystyle A}$. The convergence for the Lanczos algorithm is often orders of magnitude faster than that for the power iteration algorithm.^[9]: 477

The bounds for ${\displaystyle \theta _{1}}$ come from the above interpretation of eigenvalues as extreme values of the Rayleigh quotient ${\displaystyle r(x)}$. Since ${\displaystyle \lambda _{1}}$ is a priori the maximum of ${\displaystyle r}$ on the whole of ${\displaystyle \mathbb {C} ^{n},}$ whereas ${\displaystyle \theta _{1}}$ is merely the maximum on an ${\displaystyle m}$-dimensional Krylov subspace, we trivially get ${\displaystyle \lambda _{1}\geqslant \theta _{1}}$. Conversely, any point ${\displaystyle x}$ in that Krylov subspace provides a lower bound ${\displaystyle r(x)}$ for ${\displaystyle \theta _{1}}$, so if a point can be exhibited for which ${\displaystyle \lambda _{1}-r(x)}$ is small then this provides a tight bound on ${\displaystyle \theta _{1}}$.

The dimension ${\displaystyle m}$ Krylov subspace is

${\displaystyle \operatorname {span} \left\{v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}\right\},}$

so any element of it can be expressed as ${\displaystyle p(A)v_{1}}$ for some polynomial ${\displaystyle p}$ of degree at most ${\displaystyle m-1}$; the coefficients of that polynomial are simply the coefficients in the linear combination of the vectors ${\displaystyle v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}}$. The polynomial we want will turn out to have real coefficients, but for the moment we should allow also for complex coefficients, and we will write ${\displaystyle p^{*}}$ for the polynomial obtained by complex conjugating all coefficients of ${\displaystyle p}$. In this parametrisation of the Krylov subspace, we have

${\displaystyle r(p(A)v_{1})={\frac {(p(A)v_{1})^{*}Ap(A)v_{1}}{(p(A)v_{1})^{*}p(A)v_{1}}}={\frac {v_{1}^{*}p(A)^{*}Ap(A)v_{1}}{v_{1}^{*}p(A)^{*}p(A)v_{1}}}={\frac {v_{1}^{*}p^{*}(A^{*})Ap(A)v_{1}}{v_{1}^{*}p^{*}(A^{*})p(A)v_{1}}}={\frac {v_{1}^{*}p^{*}(A)Ap(A)v_{1}}{v_{1}^{*}p^{*}(A)p(A)v_{1}}}}$

Using now the expression for ${\displaystyle v_{1}}$ as a linear combination of eigenvectors, we get

${\displaystyle Av_{1}=A\sum _{k=1}^{n}d_{k}z_{k}=\sum _{k=1}^{n}d_{k}\lambda _{k}z_{k}}$

and more generally

${\displaystyle q(A)v_{1}=\sum _{k=1}^{n}d_{k}q(\lambda _{k})z_{k}}$

for any polynomial ${\displaystyle q}$.

Thus

${\displaystyle \lambda _{1}-r(p(A)v_{1})=\lambda _{1}-{\frac {v_{1}^{*}\sum _{k=1}^{n}d_{k}p^{*}(\lambda _{k})\lambda _{k}p(\lambda _{k})z_{k}}{v_{1}^{*}\sum _{k=1}^{n}d_{k}p^{*}(\lambda _{k})p(\lambda _{k})z_{k}}}=\lambda _{1}-{\frac {\sum _{k=1}^{n}|d_{k}|^{2}\lambda _{k}p(\lambda _{k})^{*}p(\lambda _{k})}{\sum _{k=1}^{n}|d_{k}|^{2}p(\lambda _{k})^{*}p(\lambda _{k})}}={\frac {\sum _{k=1}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})\left|p(\lambda _{k})\right|^{2}}{\sum _{k=1}^{n}|d_{k}|^{2}\left|p(\lambda _{k})\right|^{2}}}.}$

A key difference between numerator and denominator here is that the ${\displaystyle k=1}$ term vanishes in the numerator, but not in the denominator. Thus if one can pick ${\displaystyle p}$ to be large at ${\displaystyle \lambda _{1}}$ but small at all other eigenvalues, one will get a tight bound on the error ${\displaystyle \lambda _{1}-\theta _{1}}$.

Since ${\displaystyle A}$ has many more eigenvalues than ${\displaystyle p}$ has coefficients, this may seem a tall order, but one way to meet it is to use Chebyshev polynomials. Writing ${\displaystyle c_{k}}$ for the degree ${\displaystyle k}$ Chebyshev polynomial of the first kind (that which satisfies ${\displaystyle c_{k}(\cos x)=\cos(kx)}$ for all ${\displaystyle x}$), we have a polynomial which stays in the range ${\displaystyle [-1,1]}$ on the known interval ${\displaystyle [-1,1]}$ but grows rapidly outside it. With some scaling of the argument, we can have it map all eigenvalues except ${\displaystyle \lambda _{1}}$ into ${\displaystyle [-1,1]}$. Let

${\displaystyle p(x)=c_{m-1}\left({\frac {2x-\lambda _{2}-\lambda _{n}}{\lambda _{2}-\lambda _{n}}}\right)}$

(in case ${\displaystyle \lambda _{2}=\lambda _{1}}$, use instead the largest eigenvalue strictly less than ${\displaystyle \lambda _{1}}$), then the maximal value of ${\displaystyle |p(\lambda _{k})|^{2}}$ for ${\displaystyle k\geqslant 2}$ is ${\displaystyle 1}$ and the minimal value is ${\displaystyle 0}$, so

${\displaystyle \lambda _{1}-\theta _{1}\leqslant \lambda _{1}-r(p(A)v_{1})={\frac {\sum _{k=2}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})|p(\lambda _{k})|^{2}}{\sum _{k=1}^{n}|d_{k}|^{2}|p(\lambda _{k})|^{2}}}\leqslant {\frac {\sum _{k=2}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})}{|d_{1}|^{2}|p(\lambda _{1})|^{2}}}\leqslant {\frac {(\lambda _{1}-\lambda _{n})\sum _{k=2}^{n}|d_{k}|^{2}}{|p(\lambda _{1})|^{2}|d_{1}|^{2}}}.}$

Furthermore

${\displaystyle p(\lambda _{1})=c_{m-1}\left({\frac {2\lambda _{1}-\lambda _{2}-\lambda _{n}}{\lambda _{2}-\lambda _{n}}}\right)=c_{m-1}\left(2{\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}-\lambda _{n}}}+1\right);}$

the quantity

${\displaystyle \rho ={\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}-\lambda _{n}}}}$

(i.e., the ratio of the first eigengap to the diameter of the rest of the spectrum) is thus of key importance for the convergence rate here. Also writing

${\displaystyle R=e^{\operatorname {arcosh} (1+2\rho )}=1+2\rho +2{\sqrt {\rho ^{2}+\rho }},}$

we may conclude that

${\displaystyle {\begin{aligned}\lambda _{1}-\theta _{1}&\leqslant {\frac {(\lambda _{1}-\lambda _{n})\left(1-|d_{1}|^{2}\right)}{c_{m-1}(2\rho +1)^{2}|d_{1}|^{2}}}\\[6pt]&={\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {1}{\cosh ^{2}((m-1)\operatorname {arcosh} (1+2\rho ))}}\\[6pt]&={\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {4}{\left(R^{m-1}+R^{-(m-1)}\right)^{2}}}\\[6pt]&\leqslant 4{\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n})R^{-2(m-1)}\end{aligned}}}$

The convergence rate is thus controlled chiefly by ${\displaystyle R}$, since this bound shrinks by a factor ${\displaystyle R^{-2}}$ for each extra iteration.

For comparison, one may consider how the convergence rate of the power method depends on ${\displaystyle \rho }$, but since the power method primarily is sensitive to the quotient between absolute values of the eigenvalues, we need ${\displaystyle |\lambda _{n}|\leqslant |\lambda _{2}|}$ for the eigengap between ${\displaystyle \lambda _{1}}$ and ${\displaystyle \lambda _{2}}$ to be the dominant one. Under that constraint, the case that most favours the power method is that ${\displaystyle \lambda _{n}=-\lambda _{2}}$, so consider that. Late in the power method, the iteration vector:

${\displaystyle u=(1-t^{2})^{1/2}z_{1}+tz_{2}\approx z_{1}+tz_{2},}$^{[note 1]}

where each new iteration effectively multiplies the ${\displaystyle z_{2}}$-amplitude ${\displaystyle t}$ by

${\displaystyle {\frac {\lambda _{2}}{\lambda _{1}}}={\frac {\lambda _{2}}{\lambda _{2}+(\lambda _{1}-\lambda _{2})}}={\frac {1}{1+{\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}}}}}={\frac {1}{1+2\rho }}.}$

The estimate of the largest eigenvalue is then

${\displaystyle u^{*}Au=(1-t^{2})\lambda _{1}+t^{2}\lambda _{2},}$

so the above bound for the Lanczos algorithm convergence rate should be compared to

${\displaystyle \lambda _{1}-u^{*}Au=(\lambda _{1}-\lambda _{2})t^{2},}$

which shrinks by a factor of ${\displaystyle (1+2\rho )^{-2}}$ for each iteration. The difference thus boils down to that between ${\displaystyle 1+2\rho }$ and ${\displaystyle R=1+2\rho +2{\sqrt {\rho ^{2}+\rho }}}$. In the ${\displaystyle \rho \gg 1}$ region, the latter is more like ${\displaystyle 1+4\rho }$, and performs like the power method would with an eigengap twice as large; a notable improvement. The more challenging case is however that of ${\displaystyle \rho \ll 1,}$ in which ${\displaystyle R\approx 1+2{\sqrt {\rho }}}$ is an even larger improvement on the eigengap; the ${\displaystyle \rho \gg 1}$ region is where the Lanczos algorithm convergence-wise makes the smallest improvement on the power method.

Nullspace over a finite field

In 1995, Peter Montgomery published an algorithm, based on the Lanczos algorithm, for finding elements of the nullspace of a large sparse matrix over GF(2); since the set of people interested in large sparse matrices over finite fields and the set of people interested in large eigenvalue problems scarcely overlap, this is often also called the block Lanczos algorithm without causing unreasonable confusion.