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$