Quality of eigenvalue approximation in Lanczos method

Question

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$

Answer 1

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.