The Lanczos algorithm is an iterative method devised by Cornelius Lanczos that is an adaptation of power methods to find the "most useful" (tending towards extreme highest/lowest) eigenvalues and eigenvectors of an Hermitian matrix, where is often but not necessarily much smaller than . Although computationally efficient in principle, the method as initially formulated was not useful, due to its numerical instability.
In 1970, Ojalvo and Newman showed how to make the method numerically stable and applied it to the solution of very large engineering structures subjected to dynamic loading. This was achieved using a method for purifying the Lanczos vectors (i.e. by repeatedly reorthogonalizing each newly generated vector with all previously generated ones) to any degree of accuracy, which when not performed, produced a series of vectors that were highly contaminated by those associated with the lowest natural frequencies.
In their original work, these authors also suggested how to select a starting vector (i.e. use a random-number generator to select each element of the starting vector) and suggested an empirically determined method for determining , the reduced number of vectors (i.e. it should be selected to be approximately 1.5 times the number of accurate eigenvalues desired). Soon thereafter their work was followed by Paige, who also provided an error analysis. In 1988, Ojalvo produced a more detailed history of this algorithm and an efficient eigenvalue error test.
The algorithm Input a Hermitian matrix of size , and optionally a number of iterations (as default, let ). Strictly speaking, the algorithm does not need access to the explicit matrix, but only a function that computes the product of the matrix by an arbitrary vector. This function is called at most times. Output an matrix with orthonormal columns and a tridiagonal real symmetric matrix of size . If , then is unitary, and . Warning The Lanczos iteration is prone to numerical instability. When executed in non-exact arithmetic, additional measures (as outlined in later sections) should be taken to ensure validity of the results. Let be an arbitrary vector with Euclidean norm . Abbreviated initial iteration step: Let . Let . Let . For do: Let (also Euclidean norm). If , then let , else pick as an arbitrary vector with Euclidean norm that is orthogonal to all of . Let . Let . Let . Let be the matrix with columns . Let . Note for .
There are in principle four ways to write the iteration procedure. Paige and other works show that the above order of operations is the most numerically stable. In practice the initial vector may be taken as another argument of the procedure, with and indicators of numerical imprecision being included as additional loop termination conditions.
Not counting the matrix–vector multiplication, each iteration does arithmetical operations. The matrix–vector multiplication can be done in arithmetical operations where is the average number of nonzero elements in a row. The total complexity is thus , or if ; the Lanczos algorithm can be very fast for sparse matrices. Schemes for improving numerical stability are typically judged against this high performance.
The vectors are called Lanczos vectors. The vector is not used after is computed, and the vector is not used after is computed. Hence one may use the same storage for all three. Likewise, if only the tridiagonal matrix is sought, then the raw iteration does not need after having computed , although some schemes for improving the numerical stability would need it later on. Sometimes the subsequent Lanczos vectors are recomputed from when needed.
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 is an eigenvalue of , and its eigenvector (), then is a corresponding eigenvector of with the same eigenvalue:
Thus the Lanczos algorithm transforms the eigendecomposition problem for into the eigendecomposition problem for . For tridiagonal matrices, there exist a number of specialised algorithms, often with better computational complexity than general-purpose algorithms. For example, if is an tridiagonal symmetric matrix then: The continuant recursion allows computing the characteristic polynomial in operations, and evaluating it at a point in operations. The divide-and-conquer eigenvalue algorithm can be used to compute the entire eigendecomposition of in operations. The Fast Multipole Method can compute all eigenvalues in just operations. 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 just as for the divide-and-conquer algorithm (though the constant factor may be different); since the eigenvectors together have elements, this is asymptotically optimal. 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 rather than the original matrix . Since is very sparse with all nonzero elements in highly predictable positions, it permits compact storage with excellent performance vis-à-vis caching. Likewise, is a real matrix with all eigenvectors and eigenvalues real, whereas in general may have complex elements and eigenvectors, so real arithmetic is sufficient for finding the eigenvectors and eigenvalues of . If is very large, then reducing so that is of a manageable size will still allow finding the more extreme eigenvalues and eigenvectors of ; in the 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.
Aspects in which the two algorithms differ include: Lanczos takes advantage of being a sparse matrix, whereas Householder does not, and will generate fill-in. Lanczos works throughout with the original matrix (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 , whereas an iteration of Householder produces another factor in a unitary factorisation of . Each factor is however determined by a single vector, so the storage requirements are the same for both algorithms, and can be computed in time. Householder is numerically stable, whereas raw Lanczos is not. Lanczos is highly parallel, with only points of synchronisation (the computations of and ). Householder is less parallel, having a sequence of scalar quantities computed that each depend on the previous quantity in the sequence.
Derivation of the algorithm There are several lines of reasoning which lead to the Lanczos algorithm.
A more provident power method
The power method for finding the eigenvalue of largest magnitude and a corresponding eigenvector of a matrix is roughly Pick a random vector . For (until the direction of has converged) do: Let Let In the large limit, 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 , but pays attention only to the very last result; implementations typically use the same variable for all the vectors , 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 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 of a basis of such that for every and all
this is trivially satisfied by as long as is linearly independent of (and in the case that there is such a dependence then one may continue the sequence by picking as an arbitrary vector linearly independent of ). A basis containing the vectors is however likely to be numerically ill-conditioned, since this sequence of vectors is by design meant to converge to an eigenvector of . To avoid that, one can combine the power iteration with a Gram–Schmidt process, to instead produce an orthonormal basis of these Krylov subspaces. Pick a random vector of Euclidean norm . Let . For do: Let . For all let . (These are the coordinates of with respect to the basis vectors .) Let . (Cancel the component of that is in .) If then let and , otherwise pick as an arbitrary vector of Euclidean norm that is orthogonal to all of . The relation between the power iteration vectors and the orthogonal vectors is that . Here it may be observed that we do not actually need the vectors to compute these , because and therefore the difference between and is in , which is cancelled out by the orthogonalisation process. Thus the same basis for the chain of Krylov subspaces is computed by Pick a random vector of Euclidean norm . For do: Let . For all let . Let . Let . If then let , otherwise pick as an arbitrary vector of Euclidean norm that is orthogonal to all of . A priori the coefficients satisfy for all ; the definition may seem a bit odd, but fits the general pattern since
Because the power iteration vectors that were eliminated from this recursion satisfy the vectors and coefficients contain enough information from that all of 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 is Hermitian—in particular most of the coefficients turn out to be zero.
Elementarily, if is Hermitian then
For we know that , and since 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 one gets
since the latter is real on account of being the norm of a vector. For one gets
meaning this is real too.