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.

Answer 2

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.