Lanczos algorithm on arvau

Wikipedia / Lanczos algorithm

The Lanczos algorithm is an iterative method devised by Cornelius Lanczos that is an adaptation of power methods "most useful" (tending towards extreme highest/lowest) eigenvalues and eigenvectors Hermitian matrix. 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.

, 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).

times.
Output with orthonormal columns and a tridiagonal is unitary.

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.

be an arbitrary vector with Euclidean norm.

Abbreviated initial iteration step:

do:

(also Euclidean norm).
,
1. .
.
.
.
.
Note.

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. and indicators of numerical imprecision being included as additional loop termination conditions.

; the Lanczos algorithm can be very fast for sparse matrices. Schemes for improving numerical stability are typically judged against this high performance.

are called Lanczos vectors 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.

with the same eigenvalue:

tridiagonal symmetric matrix then:

The continuant recursion allows computing the characteristic polynomial operations.
The divide-and-conquer eigenvalue algorithm operations.
The Fast Multipole Method operations.
Some general eigendecomposition algorithms, notably the QR algorithm just as for the divide-and-conquer algorithm elements, this is asymptotically optimal.
Even algorithms whose convergence rates are unaffected by unitary transformations, such as the power method and inverse iteration is very sparse with all nonzero elements in highly predictable positions, it permits compact storage with excellent performance vis-à-vis caching is a real.
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:

being a sparse matrix, whereas Householder does not, and will generate fill-in.
(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 time.
Householder is numerically stable, whereas raw Lanczos is not.
points of synchronisation 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

is roughly

has converged) do:

approaches the normed eigenvector corresponding to the largest magnitude eigenvalue.

, 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.

is a chain of Krylov subspaces. One way of stating that without introducing sets into the algorithm is to claim that it computes

vectors is however likely to be numerically ill-conditioned. To avoid that, one can combine the power iteration with a Gram–Schmidt process, to instead produce an orthonormal basis of these Krylov subspaces.

do:

.
.)
.)
,
1. .

is that

, which is cancelled out by the orthogonalisation process. Thus the same basis for the chain of Krylov subspaces is computed by

do:

.
.
.
.
,
1. .

satisfy

;

since

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 coefficients turn out to be zero.

is Hermitian then

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 one gets

one gets

meaning this is real too.

is upper Hessenberg. Since

is also lower Hessenberg, so it must in fact be tridiagional. Being Hermitian, its main diagonal of the Lanczos algorithm specification.

Simultaneous approximation of extreme eigenvalues

is as stationary points of the Rayleigh quotient

and

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

is that of the gradient. In general

construct the vector spaces

such that

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

by first drawing the elements randomly according to the same normal distribution will have a uniform distribution on the unit sphere.

a diagonal matrix.

Kaniel–Paige convergence theory

. The convergence for the Lanczos algorithm is often orders of magnitude faster than that for the power iteration algorithm.

Krylov subspace is

; the coefficients of that polynomial are simply the coefficients in the linear combination. In this parametrisation of the Krylov subspace, we have

as a linear combination of eigenvectors, we get

and more generally

Thus

has coefficients, this may seem a tall order, but one way to meet it is to use Chebyshev polynomials. Let

, 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

for each extra iteration.

, so consider that. Late in the power method, the iteration vector:

The estimate of the largest eigenvalue is then

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

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 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:

Prevent the loss of orthogonality,
Recover the orthogonality after the basis is generated.
After the good and "spurious" eigenvalues are all identified, remove the spurious ones.

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

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.

ARPACK (FORTRAN 77, also available in MATLAB http://www.mathworks.com/help/techdoc/ref/eigs.html, GNU Octave (eigs), Julia https://github.com/JuliaLinearAlgebra/Arpack.jl, and Python via the SciPy https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html package) focuses on eigenvalue problems, and supports both stored and implicit matrices.

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.

Julia implementations of Lanczos and related Krylov methods can be found in Krylov.jl, KrylovKit.jl, IterativeSolvers.jl, and ArnoldiMethod.jl.

The PRIMME library also implements a Lanczos-like algorithm.

An open source Python package Leymosun implements Lanczos algorithm in the context of Krylov Complexity in pure Python: Its function called lanczos will generate Krylov bases and the coefficients.