隐式重新启动Lanczos方法

Question

我想为隐式重启Lanczos方法编写简单的玩具代码。 如果没有隐式重新启动，代码可以完美地工作，但是当我打开重启时，我无法得到正确的解决方案

据我所知，新构造的w应该与所有新的Lanczos向量正交。对于第一次重新启动，正交性得到了很好的保留，但从第二次重新启动开始，正交性被严重破坏，程序无法找到合适的特征值。在

我已经花了好几个小时修了。我差点放弃。。。。。。这是我的python代码

"""
Author: Sunghwan Choi
Date Created: June 19, 2017
Python Version: 2.7 or 3.5

Reference for Lanczos algorithm
http://www.netlib.org/utk/people/JackDongarra/etemplates/node104.html
Reference for implicit restart
http://www.netlib.org/utk/people/JackDongarra/etemplates/node118.html
"""

import numpy as np
from scipy.sparse.linalg import eigsh
#from scipy.sparse import eye
from scipy.sparse import coo_matrix
from numpy import eye

def clustering(eigvals,eigvecs,tol=1e-2):
    ret_eigvals=[]
    ret_eigvecs=[]
    for i in range(len(eigvals)):
        for ret_eigval, ret_eigvec in zip (ret_eigvals,ret_eigvecs):
            if (np.abs(eigvals[i]/ret_eigval-1.0)<tol ):
                break
        else:
            ret_eigvals.append(eigvals[i])
            ret_eigvecs.append(eigvecs[:,i])
    ret_eigvals=np.array(ret_eigvals)
    ret_eigvecs=np.array(ret_eigvecs).T
    return ret_eigvals,ret_eigvecs

def check_conv(matrix, cal_eigval, cal_eigvec, tol):
    indices=[]
    for i in range(len(cal_eigval)):
        if(np.linalg.norm(matrix.dot(cal_eigvec[:,i]) - cal_eigval[i]*cal_eigvec[:,i])< tol):
            indices.append(i)
    return indices

################ input
size=1600
max_step=20000
which='SA'
#implicit=False
implicit=True
energy_range=[0.0,6.0]
tol = 1e-5 
n_eig=6
n_tol_check=40 # n_tol_check>n_eig ==0
######################

# generate 1D harmonic oscillator
h=0.1
matrix=-5/2*eye(size)
matrix+=4/3*(eye(size,k=1)+eye(size,k=-1))
matrix+=-1/12*(eye(size,k=2)+eye(size,k=-2))
matrix=-0.5*matrix/(h*h)
distance =lambda index: (index-size/2)*h
matrix+=np.diagflat( list(map( lambda i: 0.5*distance(i)**2, range(size))))

# solve eigenvalue problem to check validity
true_eigval,true_eigvec = eigsh(matrix,k=50,which=which)
indices = np.all([true_eigval>energy_range[0], true_eigval<energy_range[1]],axis=0)
true_eigval = true_eigval[indices]
true_eigvec = true_eigvec[:,indices]

#initialize variables
alpha=[]; beta=[]
index_v=0
restart_interval = n_tol_check+n_eig if implicit is not False else max_step
T = np.zeros((restart_interval,restart_interval))
v = np.zeros((size,restart_interval))
#Q=np.eye(restart_interval)

#generate initial vector

np.random.seed(1)
initial_vec = np.random.random(size)

#initial_vec = np.loadtxt("tmp")
w = v[:,index_v] = initial_vec/np.linalg.norm(initial_vec)
init_beta = np.linalg.norm(w)
# start Lanczos i_step
for i_step in range(max_step):
    if (i_step is 0):
        v[:,index_v] = w/init_beta
    else:
        v[:,index_v] = w/T[index_v,index_v-1]
    w=matrix.dot(v[:,index_v])

    if (i_step is 0):
        w=w-init_beta*v[:,index_v-1]
    else:
        w=w-T[index_v,index_v-1]*v[:,index_v-1]
    T[index_v,index_v]=np.dot(w,v[:,index_v])
    w -=T[index_v,index_v]*v[:,index_v]

    #check convergence
    if ((i_step+1)%n_tol_check==n_eig and i_step>n_eig):
        # calculate eigenval of T matrix
        cal_eigval, cal_eigvec_= np.linalg.eigh(T[:index_v+1,:index_v+1])
        cal_eigvec = np.dot(v[:,:index_v+1],cal_eigvec_)
        #check tolerance
        conv_indices = check_conv(matrix, cal_eigval, cal_eigvec,tol)
        #filter energy_range
        indices = np.all([cal_eigval[conv_indices]>energy_range[0], cal_eigval[conv_indices]<energy_range[1]],axis=0)
        #check clustering
        conv_cal_eigval,conv_cal_eigvec = clustering((cal_eigval[conv_indices])[indices], (cal_eigvec[conv_indices])[indices])
        if (len(conv_cal_eigval)>=n_eig):
            break
        # implicit restarting 
        if (implicit is True):
            Q=np.eye(restart_interval)
            # do shift & QR decomposition
            indices = np.argsort(np.abs(cal_eigval-np.mean(energy_range)))
            for index in indices[n_eig:]:
                new_Q,new_R = np.linalg.qr(T-cal_eigval[index]*np.eye(len(T)))
                T = np.dot(new_Q.T,np.dot(T,new_Q))
                v = np.dot(v,new_Q)
                Q = np.dot(Q,new_Q)

            w=v[:,n_eig]*T[n_eig,n_eig-1]+w*Q[-1,n_eig-1]
            v[:,n_eig:]=0.0

            T[:,n_eig:] = 0.0
            T[n_eig:,:] = 0.0
            #for debug
            #print(np.dot(w.T, v))

            # reset index 
            index_v=n_eig-1
    index_v+=1
    T[index_v,index_v-1]=np.linalg.norm(w)
    T[index_v-1,index_v]=np.linalg.norm(w)
else:
    print("not converged")
    exit(-1)

print ("energy window: (", energy_range[0],",",energy_range[1],")")
print ("true eigenvalue")
print(true_eigval)
print ("eigenvalue from Lanczos w/ implicit restart (",i_step+1,")")
print(conv_cal_eigval)