Quality of eigenvalue approximation in Lanczos method

I try to familiarize myself with iterative eigenvalue solvers such as Lanczos. So I tried rewrite it to python directly according to wiki. But it doesn't seem to work.

The problem:

• it approximates quite well the largest eigenvalue $\epsilon_{max}$ of matrix $A$
• but the rest of eigenvalues obtained by solution of projected tridiagonal matrix $T$ is distributed like uniformly over the interval $( \epsilon_{min}$, $\epsilon_{max} )$ of original matrix $A$
• I woud expect that it will find $m$ larges eigenvalues
• the vectors $v_i$ (matrix $V$) generated by Lancozs iteration seems to be still orthogonal to each other, so I guess the numerical unstability of Lancozs is not the problem

Python code:

def Lanczos( A, v, m=100 ):
n = len(v)
if m>n: m = n;
# from here https://en.wikipedia.org/wiki/Lanczos_algorithm
V = np.zeros( (m,n) )
T = np.zeros( (m,m) )
vo   = np.zeros(n)
beta = 0
for j in range( m-1 ):
w    = np.dot( A, v )
alfa = np.dot( w, v )
w    = w - alfa * v - beta * vo
beta = np.sqrt( np.dot( w, w ) )
vo   = v
v    = w / beta
T[j,j  ] = alfa
T[j,j+1] = beta
T[j+1,j] = beta
V[j,:]   = v
w    = np.dot( A,  v )
alfa = np.dot( w, v )
w    = w - alfa * v - beta * vo
T[m-1,m-1] = np.dot( w, v )
V[m-1]     = w / np.sqrt( np.dot( w, w ) )
return T, V

# ---- generate matrix A
n = 50; m=10
sqrtA = np.random.rand( n,n ) - 0.5
A = np.dot( sqrtA, np.transpose(sqrtA) )

# ---- full solve for eigenvalues for reference
esA, vsA = np.linalg.eig( A )

# ---- approximate solution by Lanczos
v0   = np.random.rand( n ); v0 /= np.sqrt( np.dot( v0, v0 ) )
T, V = Lanczos( A, v0, m=m )
esT, vsT = np.linalg.eig( T )
VV = np.dot( V, np.transpose( V ) ) # check orthogonality

#print "A : "; print A
print "T : "; print T
print "VV :"; print VV
print "esA :"; print np.sort(esA)
print "esT : "; print np.sort(esT)

plt.plot( esA, np.ones(n)*0.2,  '+' )
plt.plot( esT, np.ones(m)*0.1,  '+' )
plt.ylim(0,1)
plt.show( m )


illustration:

• blue - eigenvalues of matrix $A$
• green - eigenvalues of matrix $T$

• As Bill Greene points out, this is exactly what you expect from a projection method such as Lanczos: The extremal Ritz values (i.e., the eigenvalues of the "small" projected matrix monotonically approximate the extremal eigenvalues of the matrix from inside (you can find a proof as Corollary 4.1 in Saad's book on eigenvalue problems). The only thing that I can add to his answer is that for this reason, you usually run Lanczos with $m=2k$ if you want to compute the $k$ largest (or smallest) eigenvalues. Mar 6, 2017 at 19:47

The convergence behavior you are seeing is actually expected. One of things that makes the Lanczos method so interesting is that it does a good job of simultaneously converging eigenvalues at both ends of the spectrum.

I assume your expectation of converging only the largest eigenvalues is based on the fact that, as expected from the Power iteration method, the last Lanczos vector is getting closer and closer to the eigenvector corresponding to the largest eigenvalue. But remember, each computed eigenvector is actually the "best" linear combination of all the computed Lanczos vectors. And, if you start the Lanczos process with a random vector, this vector is as good an approximation to the lowest eigenvector as the largest. That the largest eigenvalues do converge somewhat faster than the smallest, however, is consistent with what you would expect from Power iteration.

Demmel's book, Applied Numerical Linear Algebra, has a nice discussion of the convergence properties of the Lanczos algorithm in chapter 7. The plot below is computed with Demmel's Lanczos code and reproduces figure 7.2 in the book.

He computes approximate eigenvalues to a 1000 x 1000 matrix using increasing numbers of Lanczos vectors. The last column of $+$ signs shows all 1000 eigenvalues from the full matrix. (the last two columns in this plot are essentially the same as your plot above rotated 90 degrees). But this plot also shows that the distribution of the eigenvalues across the spectrum is not accidental; the eigenvalues at both extremes are converging monotonically to the extreme eigenvalues of the full matrix.

• excellent answer. one minor comment: if $A$ is symmetric, the (real) spectrum is symmetrically distributed w.r.t. $( \lambda_\max + \lambda_\min) / 2$ and one uses a random initial vector (just like the OP's and Demmel's setup), Lanczos converges just as fast to the maximum and minimum eigenvalues (as you show in your plot). For complex spectra, or eigenvalues that are not "equally well separated", you start noticing differences in the convergence speed. Mar 6, 2017 at 18:34
• @GoHokies In my experience, the smallest eigenvalues always converge slower than the largest and I thought this was the nature of the algorithm. I'd like to understand better what you're saying. Mar 6, 2017 at 23:57
• aha, so perhaps I had wrong notion that Lanczos is method of choice if I need to find accurate approximation for the first N eigenvalues of system with M eigenvalues (where N<<M) by doing ~N rather than ~M iterations. Can you just correct me what I should search for ? Mar 7, 2017 at 6:48
• @BillGreene what i mean is: Demmel's diagonal matrix elements (the eigenvalues) are sampled from a symmetric Gaussian distribution, hence they are (almost) equi-distributed around $\frac{1}{2}(\lambda_\min + \lambda_\max)$ (see figure 7.1 in his book). Because of this symmetry the Lanczos algorithm does not favor either end of the spectrum, and you notice equal rates of convergence for both the smallest and the largest rates of convergence. Suppose you break this symmetry by modifying Demmel's code to sample the diagonal elements from, say, a lognormal distribution... Mar 7, 2017 at 21:01
• ... then the behavior of the Lanczos routine should be what you expect bc a large fraction of the eigenvalues will be clustered around the lower-end of the spectrum, while the dominant (largest-magnitude) eigenvalues will be well isolated and will converge fastest. Mar 7, 2017 at 21:01

This answer is a bit late, but I was trying to implement the code you wrote and I found some issues.

You have not implemented the Lanczos algorithm correctly. The diagonalized matrix, $$T$$, itself is correct. So when you try and find the eigenvalues you obtain from the diagonalized matrix they are correct, but the matrix $$V$$ is not correct.The matrix $$V$$ is defined by Lanczos algorithm as matrix with the columns $$v_1,v_2,\ldots, v_m$$:

$$V = [v_1,v_2,\ldots, v_m]$$

The way you have it set up the matrix $$V$$ becomes,

$$V = [v_2,v_3,\ldots, v_{m+1}]$$

This is not quite what we want. The modifications you need to make to your code is,

def Lanczos( A, v, m=100 ):
n = len(v)
if m>n: m = n;
# from here https://en.wikipedia.org/wiki/Lanczos_algorithm
V = np.zeros( (m,n) )
T = np.zeros( (m,m) )
V[0, :] = v

# step 2.1 - 2.3 in https://en.wikipedia.org/wiki/Lanczos_algorithm
w = np.dot(A, v[0,:])
alfa = np.dot(w,v[0,:])
w = w - alfa*V[:, 0]
T[0,0] = alfa

# needs to start the iterations from indices 1
for j in range(1, m-1 ):
beta = np.sqrt( np.dot( w, w ) )

V[j,:] = w/beta

# This performs some rediagonalization to make sure all the vectors are orthogonal to eachother
for i in range(j-1):
V[j, :] = V[j,:] - np.dot(np.conj(V[j,:]), V[i, :])*V[i,:]
V[j, :] = V[j, :]/np.linalg.norm(V[j, :])

w = np.dot(A, V[j, :])
alfa = np.dot(w, V[j, :])
w = w - alfa * V[j, :] - beta*V[j-1, :]

T[j,j  ] = alfa
T[j-1 ,j] = beta
T[j,j-1] = beta

return T, V


There may be some small errors in my code as I quickly wrote this up. The reorthogonlization section of the code was added because it was needed for my use of the Lanczos algorithm.

The transition matrix $$V$$ may not be of any importance to you, but it can be useful in several applications. If you are wanting to perform some calculation in the Lanczos basis, you will need the matrix $$V$$ in order to transform any of your matrices from your original basis into the Lanczos basis. I needed to use this for X-ray spectroscopy simulations, in which getting the matrix $$V$$ correct was very important.

• Hi, thanks for correction. Just to be sure - the modification/correction you made does not help with obtaining good approximation of eigenvalues faster? It is only importaint for transfromation of eigenvectors? Jul 13, 2020 at 8:28
• Yes, it does not help with getting a good approximation of the eigenvalues faster. Modifications can be done to get certain values quicker for example, the largest eigenvalues and the smallest eigenvalues. However, I do not know if there is a way to speed up the code if you want to find a good approximation of the entire spectrum of eigenvalues, at least I have not been able to.
– Luke
Jul 13, 2020 at 18:05
• If you wanted to work in the original basis, but you used the Lanczos algorithm to determine your eigenvalues and eigenvectors, then you would need to use the transformation matrix to put the eigenvector back into it's original basis. If you wanted to work in the Lanczos basis then you would need to transform all your operators into the Lanczos basis using the transformation matrix, V. Transforming everything into the Lanczos basis can be very useful as it allows you to work in a smaller basis, which can result in quicker computations.
– Luke
Jul 13, 2020 at 18:14
• Sure, Lancozs is not good for getting full spectrum. I meant speed up getting good approximation of few (e.g. 10) highets eigenvalues. Jul 14, 2020 at 7:31
• Luke, are there matrices like yours on the web, to compare various methods ? How big, how sparse are they ? Thanks Oct 24, 2020 at 10:36