Lanczos algorithm explained

The Lanczos algorithm is an iterative method devised by Cornelius Lanczos that is an adaptation of power methods to find the

m

"most useful" (tending towards extreme highest/lowest) eigenvalues and eigenvectors of an

n x n

Hermitian matrix, where

m

is often but not necessarily much smaller than

n

.[1] 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.[2] 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

m

, 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.[3] [4] In 1988, Ojalvo produced a more detailed history of this algorithm and an efficient eigenvalue error test.[5]

The algorithm

A

of size

n x n

, and optionally a number of iterations

m

(as default, let

m=n

).

v\mapstoAv

that computes the product of the matrix by an arbitrary vector. This function is called at most

m

times.

Output an

n x m

matrix

V

with orthonormal columns and a tridiagonal real symmetric matrix

T=V*AV

of size

m x m

. If

m=n

, then

V

is unitary, and

A=VTV*

.

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.

  1. Let

v1\inCn

be an arbitrary vector with Euclidean norm

1

.
  1. Abbreviated initial iteration step:
    1. Let

w1'=Av1

.
    1. Let

\alpha1=

*
w
1'

v1

.
    1. Let

w1=w1'-\alpha1v1

.
  1. For

j=2,...,m

do:
    1. Let

\betaj=\|wj-1\|

(also Euclidean norm).
    1. If

\betaj0

, then let

vj=wj-1/\betaj

,

else pick as

vj

an arbitrary vector with Euclidean norm

1

that is orthogonal to all of

v1,...,vj-1

.
    1. Let

wj'=Avj

.
    1. Let

\alphaj=

*
w
j'

vj

.
    1. Let

wj=wj'-\alphajvj-\betajvj-1

.
  1. Let

V

be the matrix with columns

v1,...,vm

. Let

T=\begin{pmatrix} \alpha1&\beta2&&&&0\\ \beta2&\alpha2&\beta3&&&\\ &\beta3&\alpha3&\ddots&&\\ &&\ddots&\ddots&\betam-1&\\ &&&\betam-1&\alpham-1&\betam\\ 0&&&&\betam&\alpham\\ \end{pmatrix}

.

Note

Avj=wj'=\betaj+1vj+1+\alphajvj+\betajvj-1

for

2<j<m

.

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.[6] [7] In practice the initial vector

v1

may be taken as another argument of the procedure, with

\betaj=0

and indicators of numerical imprecision being included as additional loop termination conditions.

Not counting the matrix–vector multiplication, each iteration does

O(n)

arithmetical operations. The matrix–vector multiplication can be done in

O(dn)

arithmetical operations where

d

is the average number of nonzero elements in a row. The total complexity is thus

O(dmn)

, or

O(dn2)

if

m=n

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

The vectors

vj

are called Lanczos vectors. The vector

wj'

is not used after

wj

is computed, and the vector

wj

is not used after

vj+1

is computed. Hence one may use the same storage for all three. Likewise, if only the tridiagonal matrix

T

is sought, then the raw iteration does not need

vj-1

after having computed

wj

, although some schemes for improving the numerical stability would need it later on. Sometimes the subsequent Lanczos vectors are recomputed from

v1

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

T

, and

x

its eigenvector (

Txx

), then

y=Vx

is a corresponding eigenvector of

A

with the same eigenvalue:

\begin A y &= A V x \\ &= V T V^* V x \\ &= V T I x \\ &= V T x \\ &= V (\lambda x) \\ &= \lambda V x \\& = \lambda y.\end

Thus the Lanczos algorithm transforms the eigendecomposition problem for

A

into the eigendecomposition problem for

T

.
  1. For tridiagonal matrices, there exist a number of specialised algorithms, often with better computational complexity than general-purpose algorithms. For example, if

T

is an

m x m

tridiagonal symmetric matrix then:

O(m2)

operations, and evaluating it at a point in

O(m)

operations.

T

in

O(m2)

operations.
    • The Fast Multipole Method[8] can compute all eigenvalues in just

O(mlogm)

operations.
  1. 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

O(m2)

just as for the divide-and-conquer algorithm (though the constant factor may be different); since the eigenvectors together have

m2

elements, this is asymptotically optimal.
  1. 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

T

rather than the original matrix

A

. Since

T

is very sparse with all nonzero elements in highly predictable positions, it permits compact storage with excellent performance vis-à-vis caching. Likewise,

T

is a real matrix with all eigenvectors and eigenvalues real, whereas

A

in general may have complex elements and eigenvectors, so real arithmetic is sufficient for finding the eigenvectors and eigenvalues of

T

.
  1. If

n

is very large, then reducing

m

so that

T

is of a manageable size will still allow finding the more extreme eigenvalues and eigenvectors of

A

; in the

m\lln

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:

A

being a sparse matrix, whereas Householder does not, and will generate fill-in.

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

V

, whereas an iteration of Householder produces another factor in a unitary factorisation

Q1Q2...Qn

of

V

. Each factor is however determined by a single vector, so the storage requirements are the same for both algorithms, and

V=Q1Q2...Qn

can be computed in

O(n3)

time.

O(n)

points of synchronisation (the computations of

\alphaj

and

\betaj

). Householder is less parallel, having a sequence of

O(n2)

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

See main article: Power iteration. The power method for finding the eigenvalue of largest magnitude and a corresponding eigenvector of a matrix

A

is roughly
  1. Pick a random vector

u10

.
  1. For

j\geqslant1

(until the direction of

uj

has converged) do:
    1. Let

uj+1'=Auj.

    1. Let

uj+1=uj+1'/\|uj+1'\|.

j

limit,

uj

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

A

, but pays attention only to the very last result; implementations typically use the same variable for all the vectors

uj

, 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

uj

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

\{vj\}

m
j=1
of a basis of

\Complexn

such that

Ax\in\operatorname{span}(v1,...c,vj+1)

for every

x\in\operatorname{span}(v1,...c,vj)

and all

1\leqslantj<m;

this is trivially satisfied by

vj=uj

as long as

uj

is linearly independent of

u1,...c,uj-1

(and in the case that there is such a dependence then one may continue the sequence by picking as

vj

an arbitrary vector linearly independent of

u1,...c,uj-1

). A basis containing the

uj

vectors is however likely to be numerically ill-conditioned, since this sequence of vectors is by design meant to converge to an eigenvector of

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

u1

of Euclidean norm

1

. Let

v1=u1

.
  1. For

j=1,...c,m-1

do:
    1. Let

uj+1'=Auj

.
    1. For all

k=1,...c,j

let

gk,j=

*
v
k

uj+1'

. (These are the coordinates of

Auj=uj+1'

with respect to the basis vectors

v1,...c,vj

.)
    1. Let

wj+1=uj+1'-

j
\sum
k=1

gk,jvk

. (Cancel the component of

uj+1'

that is in

\operatorname{span}(v1,...c,vj)

.)
    1. If

wj+10

then let

uj+1=uj+1'/\|uj+1'\|

and

vj+1=wj+1/\|wj+1\|

,

otherwise pick as

uj+1=vj+1

an arbitrary vector of Euclidean norm

1

that is orthogonal to all of

v1,...c,vj

.The relation between the power iteration vectors

uj

and the orthogonal vectors

vj

is that

Auj=\|uj+1'\|uj+1=uj+1'=wj+1+

j
\sum
k=1

gk,jvk=\|wj+1\|vj+1+

j
\sum
k=1

gk,jvk

.Here it may be observed that we do not actually need the

uj

vectors to compute these

vj

, because

uj-vj\in\operatorname{span}(v1,...c,vj-1)

and therefore the difference between

uj+1'=Auj

and

wj+1'=Avj

is in

\operatorname{span}(v1,...c,vj)

, 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

v1

of Euclidean norm

1

.
  1. For

j=1,...c,m-1

do:
    1. Let

wj+1'=Avj

.
    1. For all

k=1,...c,j

let

hk,j=

*
v
k

wj+1'

.
    1. Let

wj+1=wj+1'-

j
\sum
k=1

hk,jvk

.
    1. Let

hj+1,j=\|wj+1\|

.
    1. If

hj+1,j0

then let

vj+1=wj+1/hj+1,j

,

otherwise pick as

vj+1

an arbitrary vector of Euclidean norm

1

that is orthogonal to all of

v1,...c,vj

.A priori the coefficients

hk,j

satisfy

Avj=

j+1
\sum
k=1

hk,jvk

for all

j<m

;the definition

hj+1,j=\|wj+1\|

may seem a bit odd, but fits the general pattern

hk,j=

*
v
k

wj+1'

since
*
v
j+1

wj+1'=

*
v
j+1

wj+1=\|wj+1\|

*
v
j+1

vj+1=\|wj+1\|.

Because the power iteration vectors

uj

that were eliminated from this recursion satisfy

uj\in\operatorname{span}(v1,\ldots,vj),

the vectors

\{vj\}

m
j=1
and coefficients

hk,j

contain enough information from

A

that all of

u1,\ldots,um

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

A

is Hermitian—in particular most of the

hk,j

coefficients turn out to be zero.

Elementarily, if

A

is Hermitian then

hk,j=

*
v
k

wj+1'=

*
v
k

Avj=

*
v
k

A*vj=(A

*
v
k)

vj.

For

k<j-1

we know that

Avk\in\operatorname{span}(v1,\ldots,vj-1)

, and since

vj

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

k=j-1

one gets

hj-1,j=(Avj-1)*vj=

*
\overline{v
j

Avj-1}=\overline{hj,j-1}=hj,j-1

since the latter is real on account of being the norm of a vector. For

k=j

one gets

hj,j=(A

*
v
j)

vj=

*
\overline{v
j

Avj}=\overline{hj,j

},

meaning this is real too.

More abstractly, if

V

is the matrix with columns

v1,\ldots,vm

then the numbers

hk,j

can be identified as elements of the matrix

H=V*AV

, and

hk,j=0

for

k>j+1;

the matrix

H

is upper Hessenberg. Since

H*=\left(V*AV\right)*=V*A*V=V*AV=H

the matrix

H

is Hermitian. This implies that

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,

H

is a real, symmetric matrix—the matrix

T

of the Lanczos algorithm specification.

Simultaneous approximation of extreme eigenvalues

One way of characterising the eigenvectors of a Hermitian matrix

A

is as stationary points of the Rayleigh quotient

r(x)=

x*Ax
x*x

,    x\in\Complexn.

In particular, the largest eigenvalue

λmax

is the global maximum of

r

and the smallest eigenvalue

λmin

is the global minimum of

r

.

Within a low-dimensional subspace

l{L}

of

\Complexn

it can be feasible to locate the maximum

x

and minimum

y

of

r

. Repeating that for an increasing chain

l{L}1\subsetl{L}2\subset

produces two sequences of vectors:

x1,x2,\ldots

and

y1,y2,...c

such that

xj,yj\inl{L}j

and

\begin{align} r(x1)&\leqslantr(x2)\leqslant\leqslantλmax\\ r(y1)&\geqslantr(y2)\geqslant\geqslantλmin\end{align}

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

From

xj

, the optimal direction in which to seek larger values of

r

is that of the gradient

\nablar(xj)

, and likewise from

yj

the optimal direction in which to seek smaller values of

r

is that of the negative gradient

-\nablar(yj)

. In general

\nablar(x)=

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

xj

and

yj

then there are two new directions to take into account:

Axj

and

Ayj;

since

xj

and

yj

can be linearly independent vectors (indeed, are close to orthogonal), one cannot in general expect

Axj

and

Ayj

to be parallel. It is not necessary to increase the dimension of

l{L}j

by

2

on every step if

\{l{L}j\}

m
j=1
are taken to be Krylov subspaces, because then

Az\inl{L}j+1

for all

z\inl{L}j,

thus in particular for both

z=xj

and

z=yj

.

In other words, we can start with some arbitrary initial vector

x1=y1,

construct the vector spaces

l{L}j=\operatorname{span}(x1,Ax1,\ldots,Aj-1x1)

and then seek

xj,yj\inl{L}j

such that

r(xj)=maxzj}r(z)    and    r(yj)=minzj}r(z).

Since the

j

th power method iterate

uj

belongs to

l{L}j,

it follows that an iteration to produce the

xj

and

yj

cannot converge slower than that of the power method, and will achieve more by approximating both eigenvalue extremes. For the subproblem of optimising

r

on some

l{L}j

, it is convenient to have an orthonormal basis

\{v1,\ldots,vj\}

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

A

as given, even though they are not explicitly known to the user. To fix notation, let

λ1\geqslantλ2\geqslant...b\geqslantλn

be the eigenvalues (these are known to all be real, and thus possible to order) and let

z1,...c,zn

be an orthonormal set of eigenvectors such that

Azk=λkzk

for all

k=1,...c,n

.

It is also convenient to fix a notation for the coefficients of the initial Lanczos vector

v1

with respect to this eigenbasis; let

dk=

*
z
k

v1

for all

k=1,...c,n

, so that

stylev1=

n
\sum
k=1

dkzk

. A starting vector

v1

depleted of some eigencomponent 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

v1

by first drawing the elements randomly according to the same normal distribution with mean

0

and then rescale the vector to norm

1

. Prior to the rescaling, this causes the coefficients

dk

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

(d1,...c,dn)

will have a uniform distribution on the unit sphere in

Cn

. This makes it possible to bound the probability that for example

|d1|<\varepsilon

.

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

A

a diagonal matrix with the desired eigenvalues on the diagonal; as long as the starting vector

v1

has enough nonzero elements, the algorithm will output a general tridiagonal symmetric matrix as

T

.

Kaniel–Paige convergence theory

After

m

iteration steps of the Lanczos algorithm,

T

is an

m x m

real symmetric matrix, that similarly to the above has

m

eigenvalues

\theta1\geqslant\theta2\geqslant...\geqslant\thetam.

By convergence is primarily understood the convergence of

\theta1

to

λ1

(and the symmetrical convergence of

\thetam

to

λn

) as

m

grows, and secondarily the convergence of some range

\theta1,\ldots,\thetak

of eigenvalues of

T

to their counterparts

λ1,\ldots,λk

of

A

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

The bounds for

\theta1

come from the above interpretation of eigenvalues as extreme values of the Rayleigh quotient

r(x)

. Since

λ1

is a priori the maximum of

r

on the whole of

\Complexn,

whereas

\theta1

is merely the maximum on an

m

-dimensional Krylov subspace, we trivially get

λ1\geqslant\theta1

. Conversely, any point

x

in that Krylov subspace provides a lower bound

r(x)

for

\theta1

, so if a point can be exhibited for which

λ1-r(x)

is small then this provides a tight bound on

\theta1

.

The dimension

m

Krylov subspace is

\operatorname{span}\left\{v1,Av1,A2v1,\ldots,Am-1v1\right\},

so any element of it can be expressed as

p(A)v1

for some polynomial

p

of degree at most

m-1

; the coefficients of that polynomial are simply the coefficients in the linear combination of the vectors

v1,Av1,A2v1,\ldots,Am-1v1

. 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

p*

for the polynomial obtained by complex conjugating all coefficients of

p

. In this parametrisation of the Krylov subspace, we have

r(p(A)v1)=

(p(A)
*
v
1)
Ap(A)v1
(p(A)
*
v
1)
p(A)v1

=

*
vp(A)*Ap(A)v1
1
*
vp(A)*p(A)v1
1

=

*
vp*(A*)Ap(A)v1
1
*
vp*(A*)p(A)v1
1

=

*
v
1
p*(A)Ap(A)v1
*
vp*(A)p(A)v1
1

Using now the expression for

v1

as a linear combination of eigenvectors, we get

Av1=A

n
\sum
k=1

dkzk=

n
\sum
k=1

dkλkzk

and more generally

q(A)v1=

n
\sum
k=1

dkq(λk)zk

for any polynomial

q

.

Thus

λ1-r(p(A)v1)=λ1-

*
v
n
\sum
k=1
dk
*(λ
p
k)
λkp(λk)zk
1
*
v
n
\sum
k=1
dk
*(λ
p
k)
p(λk)zk
1

=λ1-

n
\sum
2
|d
k|
λk
*
p(λ
k)
p(λk)
k=1
n
\sum
2
|d
k|
*
p(λ
k)
p(λk)
k=1

=

n
\sum
2
|d
k|
(λ1-λk)\left|p(λk)\right|2
k=1
n
\sum
2
|d
k|
\left|p(λk)\right|2
k=1

.

A key difference between numerator and denominator here is that the

k=1

term vanishes in the numerator, but not in the denominator. Thus if one can pick

p

to be large at

λ1

but small at all other eigenvalues, one will get a tight bound on the error

λ1-\theta1

.

Since

A

has many more eigenvalues than

p

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

ck

for the degree

k

Chebyshev polynomial of the first kind (that which satisfies

ck(\cosx)=\cos(kx)

for all

x

), we have a polynomial which stays in the range

[-1,1]

on the known interval

[-1,1]

but grows rapidly outside it. With some scaling of the argument, we can have it map all eigenvalues except

λ1

into

[-1,1]

. Let

p(x)=cm-1\left(

2x2n
λ2n

\right)

(in case

λ2=λ1

, use instead the largest eigenvalue strictly less than

λ1

), then the maximal value of

|

2
p(λ
k)|
for

k\geqslant2

is

1

and the minimal value is

0

, so

λ1-\theta1\leqslantλ1-r(p(A)v1)=

n
\sum
2
|d
k|
(λ1-λk)
2
|p(λ
k)|
k=2
n
\sum
k=1
2
|d
k|
2
|p(λ
k)|

\leqslant

n
\sum
2
|d
k|
(λ1-λk)
k=2
2
|d
2
|p(λ
1)|
1|

\leqslant

(λn)
n
\sum
k=2
2
|d
k|
1-λ
2
|p(λ
2
|d
1|
1)|

.

Furthermore

p(λ1)=cm-1\left(

12n
λ2n

\right)=cm-1\left(2

λ12
λ2n

+1\right);

the quantity

\rho=

λ12
λ2n

(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

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

we may conclude that

\begin{align} λ1-\theta1&\leqslant

(λn)\left(1-
2
|d
1|
\right)
1-λ
c(2\rho+1)2
2
|d
1|
m-1

\\[6pt] &=

1-
2
|d
1|
2
|d
1|

(λ1-λn)

1
\cosh2((m-1)\operatorname{arcosh

(1+2\rho))}\\[6pt] &=

1-
2
|d
1|
2
|d
1|

(λ1-λn)

4
\left(Rm-1+R-(m-1)\right)2

\\[6pt] &\leqslant4

1-
2
|d
1|
2
|d
1|

(λ1-λn)R-2(m-1)\end{align}

The convergence rate is thus controlled chiefly by

R

, since this bound shrinks by a factor

R-2

for each extra iteration.

For comparison, one may consider how the convergence rate of the power method depends on

\rho

, but since the power method primarily is sensitive to the quotient between absolute values of the eigenvalues, we need

|λn|\leqslant|λ2|

for the eigengap between

λ1

and

λ2

to be the dominant one. Under that constraint, the case that most favours the power method is that

λn=2

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

u=(1-t2)1/2z1+tz2z1+tz2,

where each new iteration effectively multiplies the

z2

-amplitude

t

by
λ2
λ1

=

λ2
λ2+(λ1-λ2)

=

1
1+
λ1-λ2
λ2

=

1
1+2\rho

.

The estimate of the largest eigenvalue is then

u*Au=

2)λ
(1-t
1

+

t
2,

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

λ1-u*Au=(λ1-λ2)t2,

which shrinks by a factor of

(1+2\rho)-2

for each iteration. The difference thus boils down to that between

1+2\rho

and

R=1+2\rho+2\sqrt{\rho2+\rho}

. In the

\rho\gg1

region, the latter is more like

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

\rho\ll1,

in which

R1+2\sqrt{\rho}

is an even larger improvement on the eigengap; the

\rho\gg1

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

v1,v2,,vm+1

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:[6] [7]

  1. Prevent the loss of orthogonality,
  2. Recover the orthogonality after the basis is generated.
  3. 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,[10] which is implemented in ARPACK.[11] This has led into a number of other restarted variations such as restarted Lanczos bidiagonalization.[12] Another successful restarted variation is the Thick-Restart Lanczos method,[13] which has been implemented in a software package called TRLan.[14]

Nullspace over a finite field

See main article: Block Lanczos algorithm.

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

A

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,[15] as well as in shell model codes in nuclear physics.[16]

Implementations

The NAG Library contains several routines[17] 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).

Similarly, in Python, the SciPy package has scipy.sparse.linalg.eigsh which is also a wrapper for the SSEUPD and DSEUPD functions functions from ARPACK which use the Implicitly Restarted Lanczos Method.

A Matlab implementation of the Lanczos algorithm (note precision issues) is available as a part of the Gaussian Belief Propagation Matlab Package. The GraphLab[18] 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.

Further reading

Notes and References

  1. Lanczos . C. . An iteration method for the solution of the eigenvalue problem of linear differential and integral operators . Journal of Research of the National Bureau of Standards. 45 . 4. 255–282 . 1950 . 10.6028/jres.045.026 . free .
  2. Ojalvo . I. U. . Newman . M. . Vibration modes of large structures by an automatic matrix-reduction method . . 8 . 7 . 1234–1239 . 1970 . 10.2514/3.5878 . 1970AIAAJ...8.1234N .
  3. Paige . C. C. . The computation of eigenvalues and eigenvectors of very large sparse matrices . U. of London . Ph.D. thesis . 1971 . 654214109 .
  4. Paige . C. C. . Computational Variants of the Lanczos Method for the Eigenproblem . J. Inst. Maths Applics . 10 . 3 . 373–381 . 1972 . 10.1093/imamat/10.3.373 .
  5. Book: Ojalvo, I. U. . Origins and advantages of Lanczos vectors for large dynamic systems . Proc. 6th Modal Analysis Conference (IMAC), Kissimmee, FL . 489–494 . 1988 .
  6. Book: Cullum . Willoughby. Lanczos Algorithms for Large Symmetric Eigenvalue Computations. 1985. 1. 0-8176-3058-9.
  7. Book: Yousef Saad. Numerical Methods for Large Eigenvalue Problems. 0-470-21820-7. Yousef Saad. 1992-06-22.
  8. Coakley. Ed S.. Rokhlin. Vladimir. A fast divide-and-conquer algorithm for computing the spectra of real symmetric tridiagonal matrices. Applied and Computational Harmonic Analysis. 2013. 34. 3. 379–414. 10.1016/j.acha.2012.06.003.
  9. Book: Golub. Gene H.. Van Loan. Charles F.. Matrix computations. 1996. Johns Hopkins Univ. Press. Baltimore. 0-8018-5413-X. 3..
  10. D. Calvetti . Daniela Calvetti . L. Reichel . D.C. Sorensen . 1994 . An Implicitly Restarted Lanczos Method for Large Symmetric Eigenvalue Problems. Electronic Transactions on Numerical Analysis. 2. 1–21.
  11. Book: R. B. Lehoucq . D. C. Sorensen . C. Yang . 1998 . ARPACK Users Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods . SIAM . 10.1137/1.9780898719628 . 978-0-89871-407-4 .
  12. E. Kokiopoulou . C. Bekas . E. Gallopoulos . 2004 . Computing smallest singular triplets with implicitly restarted Lanczos bidiagonalization . Appl. Numer. Math. . 10.1016/j.apnum.2003.11.011 . 49 . 39–61.
  13. Kesheng Wu . Horst Simon . 2000 . Thick-Restart Lanczos Method for Large Symmetric Eigenvalue Problems . SIAM . 10.1137/S0895479898334605 . SIAM Journal on Matrix Analysis and Applications . 22 . 2 . 602–616 .
  14. Web site: Kesheng Wu . Horst Simon . 2001 . TRLan software package . 2007-06-30 . https://web.archive.org/web/20070701163755/http://crd.lbl.gov/~kewu/trlan.html . 2007-07-01 . dead .
  15. Chen. HY. Atkinson, W.A. . Wortis, R. . Disorder-induced zero-bias anomaly in the Anderson-Hubbard model: Numerical and analytical calculations. Physical Review B. July 2011. 84. 4. 045113. 10.1103/PhysRevB.84.045113. 1012.1031. 2011PhRvB..84d5113C. 118722138.
  16. Shimizu. Noritaka. Nuclear shell-model code for massive parallel computation, "KSHELL". 1310.5431. 21 October 2013. nucl-th.
  17. Web site: The Numerical Algorithms Group . Keyword Index: Lanczos . NAG Library Manual, Mark 23 . 2012-02-09 .
  18. http://www.graphlab.ml.cmu.edu/pmf.html GraphLab