# The Lanczos algorithm for self-adjoint operators

## The Lanczos algorithm for self-adjoint operators

Volume 81A, number 8 PHYSICS LETTERS 16 February 1981 THE LANCZOS ALGORITHM FOR SELF-ADJOINT OPERATORS K.G. KREUZER and H.G. MILLER Institut fur Th...

Volume 81A, number 8

PHYSICS LETTERS

16 February 1981

THE LANCZOS ALGORITHM FOR SELF-ADJOINT OPERATORS K.G. KREUZER and H.G. MILLER Institut fur Theoretische Physik der Universitöt Frankfurt, 6000 Frankfurt am Main, Fed. Rep. Germany and

W.A. BERGER Metallgesellschaft Frankfurt, 6000 Frankfurt am Main, Fed. Rep. Germany Received 30 September 1 980 Revised manuscript received 1 December 1980

The method of moments for self-adjoint operators provides an algorithm for determining its low lying eigenvalues and corresponding eigenvectors. The formal structure of this algorithm is completely identical to the Lanczos algorithm when it is formally extended to self-adjoint operators. The algorithm can be used to obtain the exact low lying energy levels of a quantum-mechanical system.

1. Introduction. As is well known the description of the bound states of a system of interacting particles consists of determining the eigenvalues and eigenvectors of a (e.g. unbounded) seif-adjoint operator, i.e. the hamiltonian of the system. A traditional way of making this problem tractable is to project the operator onto a large but finite-dimensional subspace and then attempt to diagonalize the resulting matrix. On the other hand, however, it is possible to expand a seif-adjoint operator, bounded as well as unbounded, in a convergent manner via the method of moments . In the present note we want to demonstrate that in the particular case where the self-adjoint operator is explicitly assumed to possess eigenvalues the method of moments provides a means of constructing an algorithm which yields sequences of numbers and vectors which converge to the eigenvalues and eigenvectors of the operator. The formal structure and the convergence behaviour of this algorithm are completely identical to the ordinary Lanczos algorithm which is well known from numerical linear algebra [2J. The essential difference however is that one iterates with a self-adjoint (not necessarily bounded) operator ~ The advantage of this extended Lanczos algorithm is -

.

.

.

.

that by iterating with the full (unprojected) hamiltonian the exact eigenstates of the latter are obtainable in a systematic manner. In the subsequent sections we formulate the algorithm, prove its convergence properties *2 and illustrate in a numerical example how the algorithm may be applied to quantum-mechanical eigenproblems. 2. The Lanczos algorithm for self-adjoint operators. Let ~?(be a separable Hilbert space and H a self-adjoint operator (not necessarily bounded), which is supposed to possess a number (not necessarily Imite) of eigenvalues En (not necessarily bounded from above) ascending from a minimal one E 1 Furthermore let .

IT) be a vector normalized to unity and having the properties that Hk T) exists for all non-negative integers k and that T) can be expanded in terms of eigenvectors only: IT) =

~ J~,IT)

~‘

*2

tn IE~).

=

n

(1)

n the

The case where operator is supposed to be bounded has already been well discussed by Vorobyev [31. . For the sake of brevity those details of the proof which are entirely analogous to the ordinary Lanczos algorithm will not be repeated.

0 031—9163/81/0000—0000/s 02.50 © North-Holland Publishing Company

429

Volume 81A, number 8

PHYSICS LETTERS

Here P

16 February 1981

0 denotes the projector onto the eigenspace~,~ corresponding to the eigenvalue F0 and iE~>is a corresponding normalized eigenvector. Now let!: ill <.... Eq. (1) then reduces to:

Eq. (4) indicates that the vectors (l~>are obtained by orthonormalizing the set {I T>. HIT), 122 IT>. ~via the Schmidt procedure. The self-adjointness of H guar. antees that all matrix elements (F~IIIVl~> are real atid vanish for IX X’I > 1. Therefore the vectors VP~>ful-

IT> = ~ tnyIEn >.

kb2> = i.JI ~HI~1>

(2)

The trial vector singles out a specific part of the spec tral decomposition of H, i.e. the operator:

...

--

fill the following recursion relations:

i

HI~1>),

~l >(k1

(7a)

(HI~>

‘~‘~x+1> = (7h)

Hr

:

~ EnyIEn~XEnyI.

(3)

and specifies a subspace ‘?~‘Twhich is defined to be the closed linear manifold spanned by the eigenvectors IE~~> and is invariant with respect to H and HT, respectively *3 The algorithm can now be stated as follows. Calculate successively the vectors:

IT>,

I~l>:=

(~T) ~

I~x’>(~x,IHxIT>),

X = 2, 3,4

and fII~)> Icl1>K~1jVb1>+(b7>K(1~2H4)1>, ~I~> +

=

(8a)

‘~xl>~~xl lPi~~> +

I+1>(~~+1IH~~>, X = 2: 3,4

(Sb)

and

(4a)

:=

1HI~,

~

=

p~~)IT> X

=

1, 2

(~)

I lere the polynomials p~(X)are defined in the following way:

X1

(4b) Xl,2,3

p1(x):=l,

1[(x

and diagonalize successively the operators: :=

~

~

p2(x):=wj~(x~v1).

1=1,2.3,...

(5)

--

v~)p~(x)

(lOa)

(lob)

Px+i(x) := w~ X 2,3,4 where

The resulting sequences of eigenvalues (e 11), (e12), (~I~lHkF~>,w~ ~+1iHi~~>. (11) The orthogonal and normalized vectors Vk~>are corn-

and the corresponding sequences of eigenvectors (Ic11 >), (Ic12>), ... of H1 then possess the following convergence properties:

Vy

e~y~En~:

for I~>~ ~r with I~>I~it follows that I~>= 0>. Here P denotes the closed linear manifold spanned by the vectors VI~>.Obviously .1? is contained in ~T• Let I~>E ~T with I~>I P. Then in particular K~IHkIT> = 0 for all non-negative integers k and therefore

~_=~

X

1,2,3

IE~), X =

1,

2, 3

(6a)

(6b)

In the case where Ii possesses only a finite number of eigenvalues or if IT) contains only a finite number of eigenvectors then ~Tis of finite dimension and HT can be considered as an ordinary hermitian operator defined in this space and its eigenvalues and eigenvectors can be obtained via the ordinary Lanczos algorithm. Our purpose of course is to investigate the case where the trial vector contains an infinite number of eigenvectors and this is assumed in the following, 430

plete in ~

To show this it is sufficient to show that

(~Iq(~)IT> = 0 for every operator-valued polynomial q, since ~klT) and q(H)IT> are both elements ofP. Evaluating this explicitly yields ~l q(E~~)t0~c~ = 0 where c~:= (i/iIE~~. Since this equation must hold for every polynomial q and because all t~ ~ 0, it can he inferred that all c~= 0 and thus I/i> = ID>. Expanding the eigenvectors vectors kI~>one obtains:

IE0~>in

terms

of the

Volume 81A, number 8

PHYSICS LETTERS

and therefore IIIRx~(Hi)III= which contradicts the vanishing of the norm of the operator Rx~(H,)

(~

2 IpX~(EnX)I2)I

lEn~>=

p~(E,1~)I~>~ (12)

where the normalization constant is just equal to the scalar product between En) and IT>. For the characteristic polynomial which corresponds to the eigenproblem of H1 one obtains det(H1

x”l) = (—l)’w1

...

w,p,+1(x).

(13)

For the zeros of these polynomials it is guaranteed that they separate one another, i.e. for the eigenvalues clx one obtains:

el+lx
(14)

.

For the decomposition of the eigenvectors elx) in terms of the vectors ‘~x>one obtains:

Clx>

=

(

Ipx(elx)12) 1/2

16 February 1981

px’(eix)I~x’>, (15)

where again the normalization constant is equal to the scalar product between clx) and IT). Relation (14) guarantees that the sequences (elx) decrease monotonically and therefore must be convergent to a limit x~since they all are bounded from below. The latter follows from the fact that e,x is the expectation value (e,~IHTIeIx) which is always greater than E01 , the smallest eigenvalue of HT. From eqs. (12) and (15) it follows Ielx> —~ofIEnx> as i provided that Xx = Enx. that The proof the latter

Rxx(~T)in the limit 1 ~. Hence the assumption must be wrong, i.e. x~must be an eigenvalueofH~. The convergence property of the sequence (H,) guarantees furthermore that for each En~there exists a convergent sequence (ei’) where (ei’) is contained in the spectrum of H1 and e~’ E0~as 1 . Now, among all the convergent sequences which are formed from the numbers e1 ~ e21 , e22, e31 , e32, e33, ... , the sequences (e11), (e12), (e13), ... are obviously those Sequences with the lowest, second lowest, third lowest, ... limit. This together with the fact that xx is indeed an eigenvalue of HT then yieldsx~= E,~for all X, provided that x~~ xx for all A ~ A’. That the latter is valid can easily be proved by contradiction. i.e. supposing xx = xx and therefore Xx> = IXx’> yields

1

=

=

lim (eixle,x) = 0.

3. Application of the algorithm to a quantum-mechanical hamiltonian In order to demonstrate how the Lanczos algorithm may be applied to a quantum-mechanical system we consider the hamiltonian of the anharmonic oscillator in one dimension. In the coordinate representation the eigenproblem of this hamiltonian is given by:

2/2m) d2/dx2 + x4] (xIE> = E(xIE>. (16) [—(h As a trial state one may take, for example, a superposi-

based in essence on the fact that the sequence (H

1)

tion of normalized oscillator-type basis states, i.e.

converges to H in ~T in the norm resolvent sense  *~ Namely ifit is assumed thatxx is not equal to an eigenvalue ofHT then it is guaranteed that [41*5

IT> = II

IIIRxx(Uhi)

where

Rxx(HT)III

0

—~

as I

°o.

-~

1

~

Ie 1xXeix’I,

(17)

and ~ is an arbitrary positive parameter which in part

and

IIRxx(Hi)Ieix)II

c~I~p~>,

2/2). (18) (xl~>= ;x~exp(—~3yx Here ~ are normalization constants, ‘y = (/12/2m)—h/2,

On the other hand,however, we have Rxx(Hl) = ~ x’=i ~

K

determines the asymptotic behaviour of (xI~)and =

xx

1 —

e/x

—~ °°

as I

-+

00,

Cf- also the definition and theorem VIII.19 in ref. [41, ~. 284. ~ lIlA Ill denotes the norm of an operator A and R~(A) (z — A)~, the resolvent.

(xIT), respectively. This trial state possesses the required property that every power of the hamiltonian may be applied to it. Furthermore, this trial state has the advantagesince thatHXIT) I~> may be calculated analytically, for the and present hamiltonian and the aforementioned choice of basis states the expansion coefficients a~and b~of the vectors

431

Volume 81A, number 8

PHYSICS LETTERS

A =

v=0

In order to apply the algorithm to a given self-ad-

A+4

a p >,

HI~>=

joint operator, an appropriate trial vector must be available, i.e. a vector which belongs to the domains

b I~> V

V0

16 February 1981

V

are related through a recursion relation. This simplifies the numerical calculation considerably. Choosing as the trial state IT> = Ip

0> [i.e. setting K Din eq. (17)] and 13 = 3.5 yields after relatively few iterations the energies of the lowest even-parity eigenstates. In the 10th iteration step, for example, the values of the first three even-parity eigenenergies agree to 5 significant fIgures with the results obtained by Hirschfelder . As can be seen from fig. 1 the convergence rate is the same as for the ordinary Lanczos algorithm for hermitian matrices, i.e. it is more rapid for the lower lying eigenvalues. =

4. Summary and discussion. In the present work we have attempted to demonstrate how the method of moments for a self-adjoint operator may be used to determine eigenvalues and eigenvectors of the latter, in particular in the case where the operator is not necessarily bounded. The resulting algorithm is identical to the Lanczos algorithm when it is formally extended to self-adjoint operators. ______________________________________

of all powers of the operator and which can be expanded in terms of eigenvectors only. Though for self-adjoint operators vectors satisfying the first requirernent exist and even form a dense subset in the Hilbert space ,for a realistic haniiltonian such a trial state need not necessarily be readily available. The second requirement may also cause some difficulties, since the trial state is not allowed to contain scattering states. If it does then it is not intuitively ohvious whether the aforementioned convergence properties remain valid. On the other hand, the possibility of being able to iterate with the full (unprojected) harniltonian in a manner which should yield after a relatively small number of iterations the exact lower lying eigenstates is intriguing even though at least at present the applicability of the algorithm appears to be restricted to quantum-mechanical systems with a small number of degrees of freedom. Finally we would like to mention that if one is only interested in the ground state of the full (unprojected) hamiltonian, then one can use instead of the extended Lanczos algorithm the extended (2 X 2) algorithm .

50

40

We wish to thank the Computing Center of the University of Frankfurt for making their facilities available. One of us (K.G.K.) would like to thank Professor R.M. Dreizler for stimulating his interest in the subject

of the present paper.

£ 6

20 —

£4

£7

10

References Ill D. Masson, in: The Padé approximant in theoretical physics, ed G.A. Baker (Academic Press, New York, 1970) p.197. 121 DK. Faddeew and W.N. Faddejewa, Numerische Methoden der linearen Algebra (Oldenburg, Munich, 1964).

131 U.V. Vorobyev, Moments method in applied mathematics 1

I

1

2

3

I

4

5 ITERATIONS

I

6

7

I

9

10

-

Fig. 1. The convergence rate of the sequences (eli) (el5) converging to the first five even-parity eigenenergies 2/2m = i - of the hamiltonian. Units are chosen such that h

432

t~

(Gordon and Breach, New York, 1962).

[41 M. Reed and B. Simon, Methods of modern mathematical physics I (Academic Press, New York, 1972). 151 D. Secrest, K. Cashion and JO. Hirschfelder, J. Chem. Phys. 37 (1962) 830. [61 E. Nelson, Ann. Math. 70 (1959) 573. [71 1. KG. Phys. Kreuzer, Al 3 (1980) l-I.G. Miller, 2645. R.M. Dreizier and W.A. Berger,