Lanczos algorithm

The Lanczos algorithm is a direct algorithm devised by Cornelius Lanczos that is an adaptation of power methods to find the formula_1 "most useful" (tending towards extreme highest/lowest) eigenvalues and eigenvectors of an formula_2 Hermitian matrix, where formula_3 is often but not necessarily much smaller than formula_4. 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 formula_3, 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.
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 formula_39 may be taken as another argument of the procedure, with formula_40 and indicators of numerical imprecision being included as additional loop termination conditions.
Not counting the matrix–vector multiplication, each iteration does formula_41 arithmetical operations. The matrix–vector multiplication can be done in formula_42 arithmetical operations where formula_43 is the average number of nonzero elements in a row. The total complexity is thus formula_44, or formula_45 if formula_9; the Lanczos algorithm can be very fast for sparse matrices. Schemes for improving numerical stability are typically judged against this high performance.
The vectors formula_28 are called "Lanczos vectors". 
The vector formula_48 is not used after formula_49 is computed, and the vector formula_49 is not used after formula_51 is computed. Hence one may use the same storage for all three. Likewise, if only the tridiagonal matrix formula_52 is sought, then the raw iteration does not need formula_53 after having computed formula_49, although some schemes for improving the numerical stability would need it later on. Sometimes the subsequent Lanczos vectors are recomputed from formula_39 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 formula_56 is an eigenvalue of formula_6, and if formula_58 (formula_59 is an eigenvector of formula_52) then formula_61 is the corresponding eigenvector of formula_6 (since formula_63). Thus the Lanczos algorithm transforms the eigendecomposition problem for formula_6 into the eigendecomposition problem for formula_52.
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:
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 formula_6 is roughly
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 formula_6, but pays attention only to the very last result; implementations typically use the same variable for all the vectors formula_100, having each new iteration overwrite the results from the previous one. What if we instead kept all the intermediate results and organised their data?
One piece of information that trivially is available from the vectors formula_100 is a chain of Krylov subspaces. One way of stating that without introducing sets into the algorithm is to claim that it computes
this is trivially satisfied by formula_113 as long as formula_100 is linearly independent of formula_115 (and in the case that there is such a dependence then one may continue the sequence by picking as formula_28 an arbitrary vector linearly independent of formula_115). A basis containing the formula_100 vectors is however likely to be numerically ill-conditioned, since this sequence of vectors is by design meant to converge to an eigenvector of formula_6. To avoid that, one can combine the power iteration with a Gram–Schmidt process, to instead produce an orthonormal basis of these Krylov subspaces. 
The relation between the power iteration vectors formula_100 and the orthogonal vectors formula_28 is that
Here it may be observed that we do not actually need the formula_100 vectors to compute these formula_28, because formula_143 and therefore the difference between formula_124 and formula_145 is in formula_131, which is cancelled out by the orthogonalisation process. Thus the same basis for the chain of Krylov subspaces is computed by
A priori the coefficients formula_160 satisfy
the definition formula_154 may seem a bit odd, but fits the general pattern formula_164 since 
Because the power iteration vectors formula_100 that were eliminated from this recursion satisfy formula_167 the vectors formula_108 and coefficients formula_160 contain enough information from formula_6 that all of formula_171 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 formula_6 is Hermitian—in particular most of the formula_160 coefficients turn out to be zero.
Elementarily, if formula_6 is Hermitian then 
For formula_176 we know that formula_177, and since formula_178 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 formula_179 one gets 
since the latter is real on account of being the norm of a vector. For formula_181 one gets 
meaning this is real too.
More abstractly, if formula_13 is the matrix with columns formula_184 then the numbers formula_160 can be identified as elements of the matrix formula_186, and formula_187 for formula_188 the matrix formula_189 is upper Hessenberg. Since 
the matrix formula_189 is Hermitian. This implies that formula_189 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, formula_189 is a real, symmetric matrix—the matrix formula_52 of the Lanczos algorithm specification.
Simultaneous approximation of extreme eigenvalues.
One way of characterising the eigenvectors of a Hermitian matrix formula_6 is as stationary points of the Rayleigh quotient
In particular, the largest eigenvalue formula_197 is the global maximum of formula_198 and the smallest eigenvalue formula_199 is the global minimum of formula_198.
Within a low-dimensional subspace formula_201 of formula_109 it can be feasible to locate the maximum formula_59 and minimum formula_204 of formula_198. Repeating that for an increasing chain formula_206 produces two sequences of vectors: formula_207 and formula_208 such that formula_209 and 
The question then arises how to choose the subspaces so that these sequences converge at optimal rate.
From formula_211, the optimal direction in which to seek larger values of formula_198 is that of the gradient formula_213, and likewise from formula_214 the optimal direction in which to seek smaller values of formula_198 is that of the negative gradient formula_216. In general
so the directions of interest are easy enough to compute in matrix arithmetic, but if one wishes to improve on both formula_211 and formula_214 then there are two new directions to take into account: formula_220 and formula_221 since formula_211 and formula_214 can be linearly independent vectors (indeed, are close to orthogonal), one cannot in general expect formula_220 and formula_225 to be parallel. Is it therefore necessary to increase the dimension of formula_226 by formula_227 on every step? Not if formula_228 are taken to be Krylov subspaces, because then formula_229 for all formula_230 thus in particular for both formula_231 and formula_232.
In other words, we can start with some arbitrary initial vector formula_233 construct the vector spaces
and then seek formula_235 such that 
Since the formula_103th power method iterate formula_100 belongs to formula_239 it follows that an iteration to produce the formula_211 and formula_214 cannot converge slower than that of the power method, and will achieve more by approximating both eigenvalue extremes. For the subproblem of optimising formula_198 on some formula_243, it is convenient to have an orthonormal basis formula_244 for this vector space. Thus we are again led to the problem of iteratively computing such a basis for the sequence of Krylov subspaces.
Convergence and other dynamics.
When analysing the dynamics of the algorithm, it is convenient to take the eigenvalues and eigenvectors of formula_6 as given, even though they are not explicitly known to the user. To fix notation, let formula_246 be the eigenvalues (these are known to all be real, and thus possible to order) and let formula_247 be an orthonormal set of eigenvectors such that formula_248 for all formula_249.
It is also convenient to fix a notation for the coefficients of the initial Lanczos vector formula_39 with respect to this eigenbasis; let formula_251 for all formula_249, so that formula_253. A starting vector formula_39 depleted of some eigenvalue will delay convergence to the corresponding eigenvalue, and even though this just comes out as a constant factor in the error bounds, depletion remains undesirable. One common technique for avoiding being consistently hit by it is to pick formula_39 by first drawing the elements randomly according to the same normal distribution with mean formula_256 and then rescale the vector to norm formula_20. Prior to the rescaling, this causes the coefficients formula_258 to also be independent normally distributed stochastic variables from the same normal distribution (since the change of coordinates is unitary), and after rescaling the vector formula_259 will have a uniform distribution on the unit sphere in formula_260. This makes it possible to bound the probability that for example formula_261.
The fact that the Lanczos algorithm is coordinate-agnostic – operations only look at inner products of vectors, never at individual elements of vectors – makes it easy to construct examples with known eigenstructure to run the algorithm on: make formula_6 a diagonal matrix with the desired eigenvalues on the diagonal; as long as the starting vector formula_39 has enough nonzero elements, the algorithm will output a general tridiagonal symmetric matrix as formula_52.
Kaniel–Paige convergence theory.
After formula_1 iteration steps of the Lanczos algorithm, formula_52 is an formula_15 real symmetric matrix, that similarly to the above has formula_1 eigenvalues formula_269 By convergence is primarily understood the convergence of formula_270 to formula_271 (and the symmetrical convergence of formula_272 to formula_273) as formula_1 grows, and secondarily the convergence of some range formula_275 of eigenvalues of formula_52 to their counterparts formula_277 of formula_6. The convergence for the Lanczos algorithm is often orders of magnitude faster than that for the power iteration algorithm.
The bounds for formula_270 come from the above interpretation of eigenvalues as extreme values of the Rayleigh quotient formula_280. Since formula_271 is a priori the maximum of formula_198 on the whole of formula_283 whereas formula_270 is merely the maximum on an formula_1-dimensional Krylov subspace, we trivially get formula_286. Conversely, any point formula_59 in that Krylov subspace provides a lower bound formula_280 for formula_270, so if a point can be exhibited for which formula_290 is small then this provides a tight bound on formula_270.
The dimension formula_1 Krylov subspace is 
so any element of it can be expressed as formula_294 for some polynomial formula_295 of degree at most formula_296; the coefficients of that polynomial are simply the coefficients in the linear combination of the vectors formula_297. 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 formula_298 for the polynomial obtained by complex conjugating all coefficients of formula_295. In this parametrisation of the Krylov subspace, we have
Using now the expression for formula_39 as a linear combination of eigenvectors, we get
and more generally 
for any polynomial formula_304.
Thus
A key difference between numerator and denominator here is that the formula_306 term vanishes in the numerator, but not in the denominator. Thus if one can pick formula_295 to be large at formula_271 but small at all other eigenvalues, one will get a tight bound on the error formula_309.
Since formula_6 has many more eigenvalues than formula_295 has coefficients, this may seem a tall order, but one way to meet it is to use Chebyshev polynomials. Writing formula_312 for the degree formula_313 Chebyshev polynomial of the first kind (that which satisfies formula_314 for all formula_59), we have a polynomial which stays in the range formula_316 on the known interval formula_316 but grows rapidly outside it. With some scaling of the argument, we can have it map all eigenvalues except formula_271 into formula_316. Let
(in case formula_321, use instead the largest eigenvalue strictly less than formula_271), then the maximal value of formula_323 for formula_324 is formula_20 and the minimal value is formula_256, so
Furthermore
the quantity 
(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 
we may conclude that
The convergence rate is thus controlled chiefly by formula_332, since this bound shrinks by a factor formula_333 for each extra iteration.
For comparison, one may consider how the convergence rate of the power method depends on formula_334, but since the power method primarily is sensitive to the quotient between absolute values of the eigenvalues, we need formula_335 for the eigengap between formula_271 and formula_337 to be the dominant one. Under that constraint, the case that most favours the power method is that formula_338, so consider that. Late in the power method, the iteration vector:
where each new iteration effectively multiplies the formula_340-amplitude formula_341 by
The estimate of the largest eigenvalue is then 
so the above bound for the Lanczos algorithm convergence rate should be compared to 
which shrinks by a factor of formula_345 for each iteration. The difference thus boils down to that between formula_346 and formula_347. In the formula_348 region, the latter is more like formula_349, and performs like the power method would with an eigengap twice as large; a notable improvement. The more challenging case is however that of formula_350 in which formula_351 is an even larger improvement on the eigengap; the formula_352 region is where the Lanczos algorithm convergence-wise makes the "smallest" improvement on the power method.
Numerical stability.
Stability means how much the algorithm will be affected (i.e. will it produce the approximate result close to the original one) if there are small numerical errors introduced and accumulated. Numerical stability is the central criterion for judging the usefulness of implementing an algorithm on a computer with roundoff.
For the Lanczos algorithm, it can be proved that with "exact arithmetic", the set of vectors formula_353 constructs an "orthonormal" basis, and the eigenvalues/vectors solved are good approximations to those of the original matrix. However, in practice (as the calculations are performed in floating point arithmetic where inaccuracy is inevitable), the orthogonality is quickly lost and in some cases the new vector could even be linearly dependent on the set that is already constructed. As a result, some of the eigenvalues of the resultant tridiagonal matrix may not be approximations to the original matrix. Therefore, the Lanczos algorithm is not very stable.
Users of this algorithm must be able to find and remove those "spurious" eigenvalues. Practical implementations of the Lanczos algorithm go in three directions to fight this stability issue:
Variations.
Variations on the Lanczos algorithm exist where the vectors involved are tall, narrow matrices instead of vectors and the normalizing constants are small square matrices. These are called "block" Lanczos algorithms and can be much faster on computers with large numbers of registers and long memory-fetch times.
Many implementations of the Lanczos algorithm restart after a certain number of iterations. One of the most influential restarted variations is the implicitly restarted Lanczos method, which is implemented in ARPACK. This has led into a number of other restarted variations such as restarted Lanczos bidiagonalization. Another successful restarted variation is the Thick-Restart Lanczos method, which has been implemented in a software package called TRLan.
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.
Applications.
Lanczos algorithms are very attractive because the multiplication by formula_354 is the only large-scale linear operation. Since weighted-term text retrieval engines implement just this operation, the Lanczos algorithm can be applied efficiently to text documents (see Latent Semantic Indexing). Eigenvectors are also important for large-scale ranking methods such as the HITS algorithm developed by Jon Kleinberg, or the PageRank algorithm used by Google.
Lanczos algorithms are also used in Condensed Matter Physics as a method for solving Hamiltonians of strongly correlated electron systems, as well as in shell model codes in nuclear physics.
Implementations.
The NAG Library contains several routines for the solution of large scale linear systems and eigenproblems which use the Lanczos algorithm.
MATLAB and GNU Octave come with ARPACK built-in. Both stored and implicit matrices can be analyzed through the "eigs()" function (Matlab/Octave).
A Matlab implementation of the Lanczos algorithm (note precision issues) is available as a part of the Gaussian Belief Propagation Matlab Package. The GraphLab collaborative filtering library incorporates a large scale parallel implementation of the Lanczos algorithm (in C++) for multicore.
The PRIMME library also implements a Lanczos like algorithm.