4
$\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
11
$\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.

$\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

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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