计算科学 - Lanczos 方法中特征值逼近的质量 - 吾爱随笔录

Lanczos 方法中特征值逼近的质量

计算科学特征值迭代法本征系统

2021-12-09 07:20:14

我尝试熟悉诸如Lanczos之类的迭代特征值求解器。所以我尝试根据wiki直接将其重写为python。但这似乎不起作用。

问题：

它非常接近最大的特征值 $\epsilon_{max}$ 矩阵 $A$
但是通过投影三对角矩阵的解获得的其余特征值在区间上均匀分布,原始矩阵
- 我希望它会找到 $m$ 大特征值
向量 $v_i$ （矩阵 $V$ ）由Lancozs迭代生成的似乎仍然相互正交，所以我猜Lancozs的数值不稳定不是问题

Python代码：

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 )

插图：

蓝色 - 矩阵的特征值 $A$
绿色 - 矩阵的特征值 $T$

2个回答

您看到的收敛行为实际上是预期的。使 Lanczos 方法如此有趣的一件事是它在同时收敛频谱两端的特征值方面做得很好。

我假设您对仅收敛最大特征值的期望是基于这样一个事实，即正如 Power 迭代方法所预期的那样，最后一个 Lanczos 向量越来越接近与最大特征值对应的特征向量。但请记住，每个计算的特征向量实际上是所有计算的 Lanczos 向量的“最佳”线性组合。并且，如果您使用随机向量启动 Lanczos 过程，则该向量与最大特征向量的近似值一样好。然而，最大的特征值确实比最小的特征值收敛得更快，这与您对 Power 迭代的期望是一致的。

Demmel 的书Applied Numerical Linear Algebra在第 7 章中对 Lanczos 算法的收敛特性进行了很好的讨论。下图是使用 Demmel 的 Lanczos 代码计算的，并重现了书中的图 7.2。

他使用越来越多的 Lanczos 向量来计算 1000 x 1000 矩阵的近似特征值。最后一栏 $+$ 符号显示来自完整矩阵的所有 1000 个特征值。（此图中的最后两列与旋转 90 度上方的图基本相同）。但该图也表明，特征值在整个频谱上的分布并非偶然；两个极端的特征值都单调地收敛到整个矩阵的极端特征值。

这个答案有点晚了，但我试图实现你写的代码，我发现了一些问题。

您没有正确实现 Lanczos 算法。对角化矩阵， $T$ ，本身是正确的。因此，当您尝试找到从对角化矩阵中获得的特征值时，它们是正确的，但是矩阵 $V$ 不正确。矩阵 $V$ 由Lanczos 算法定义为具有列的矩阵 $v_1,v_2,\ldots, v_m$ ：

V = [v_{1}, v_{2}, \dots, v_{m}]

$V = [v_1,v_2,\ldots, v_m]$

你设置矩阵的方式 $V$ 变成，

V = [v_{2}, v_{3}, \dots, v_{m + 1}]

$V = [v_2,v_3,\ldots, v_{m+1}]$

这不是我们想要的。您需要对代码进行的修改是，

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

当我快速写下这段代码时，我的代码中可能存在一些小错误。添加了代码的重新正交化部分，因为我使用 Lanczos 算法需要它。

转移矩阵 $V$ 可能对您来说并不重要，但它可以在多个应用程序中有用。如果您想在 Lanczos 基础上执行一些计算，您将需要矩阵 $V$ 为了将您的任何矩阵从您的原始基础转换为 Lanczos 基础。我需要将其用于 X 射线光谱模拟，其中获取矩阵 $V$ 正确非常重要。

其它你可能感兴趣的问题

上一篇高阶泽尼克多项式的数值稳定性下一篇Cholesky 分解的全等级更新