8
$\begingroup$

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$

enter image description here

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Christian Clason Mar 6 '17 at 19:47
13
$\begingroup$

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.

enter image description here

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.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – GoHokies Mar 6 '17 at 18:34
  • $\begingroup$ @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. $\endgroup$ – Bill Greene Mar 6 '17 at 23:57
  • $\begingroup$ 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 ? $\endgroup$ – Prokop Hapala Mar 7 '17 at 6:48
  • 1
    $\begingroup$ @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... $\endgroup$ – GoHokies Mar 7 '17 at 21:01
  • 1
    $\begingroup$ ... 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. $\endgroup$ – GoHokies Mar 7 '17 at 21:01
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Prokop Hapala Jul 13 at 8:28
  • $\begingroup$ 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. $\endgroup$ – Luke Jul 13 at 18:05
  • $\begingroup$ 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. $\endgroup$ – Luke Jul 13 at 18:14
  • $\begingroup$ Sure, Lancozs is not good for getting full spectrum. I meant speed up getting good approximation of few (e.g. 10) highets eigenvalues. $\endgroup$ – Prokop Hapala Jul 14 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.