数学系列：Lanczos算法

Lanczos算法的大致思路

为了求\(m\)阶对称方阵\(A\)最大的\(r\)个特征值和特征向量：

\[A_{m\times m}\approx U_{m_\times r}T_{r\times r}U^T_{m\times r} \]

其中\(U\)是列正交矩阵，即 \(U^TU=I\)，每一列为一个特征向量，\(T\)是三对角阵。\(r\)为分解的秩。对矩阵\(T\)进行奇异值分解，从而得到矩阵\(A\)最大的\(r\)个特征值和特征向量。

Lanczos算法求解步骤:

1）取定Lanczos基向量的数量\(k\)，构造Krylov队列\(\{u_1,Au_1,A^2u_1,\ldots,A^ku_1\}\) ，利用施密特正交化将Krylov队列转化为标准正交向量组构成矩阵\(U_{m\times k}\)。注意：\(u_1\)可以取任意非零向量，可取\(u_1=1_m\)或者随机向量。

2）将矩阵\(A\)转换为三对角阵\(T\)：\(T=U_{m\times k}^TA_{m\times m}U_{m\times k}\)，其数学形式为

\[T=\begin{bmatrix} \alpha_0 & \beta_1 & 0 & 0 & 0 & 0 & \ldots\\ \beta_1 & \alpha_1 & \beta_2 & 0 & 0 & 0 & \ldots\\ 0 & \beta_2 & \alpha_3 & \beta_3 & 0 & 0 & \ldots\\ 0 & 0 & \beta_3 & \alpha_4 & \beta_4 & 0 & \ldots\\ 0 & 0 & 0 & \beta_4 & \alpha_5 & \beta_5 & \ldots\\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \end{bmatrix} \]

3）对矩阵\(T\)进行svd分解：\(T_{k\times k} = U*S*V'\)，取对角矩阵\(S\)的最大的前\(r\)个特征值，即为估计的矩阵\(A\)的最大的前\(r\)个特征值。

仿真例子

​ 随机生成一个\(50\times 50\)的对称阵\(A\)，求取其前3个最大特征值。令\(k=10\)，\(u_1=1_{50}\)，根据步骤1）计算得到\(U_{50\times 10}\)，这里设定一个阈值\(1e-5\)，计算\(T=U^T*A*U\)，当其中元素小于阈值则置为0。对\(T\)进行svd分解，取前三个最大特征值。对\(A\)进行svd分解，取前三个特征值进行对比。得到结果

矩阵\(T\)的前十个特征值	矩阵\(A\)的前十个特征值
2.1000	2.1000
1.7706	1.7729
1.5678	1.5773
1.1686	1.5468
0.9834	1.4200
0.7027	1.2849
0.5092	1.1673
0.2193	1.1377
0.0744	1.0842
0.0127	0.9785

可以看出，从第四个特征值开始就变得很不准确了，增大\(k\)可以得到后面更多准确的特征值，但是却很难保证Krylov队列满足列满秩。总结一下有以下几点结论。

• Lanczos算法通过构造矩阵\(T\)求解前\(r\)个特征值，比用传统的svd分解矩阵\(A\)要快的多，很大原因在于矩阵\(T\)不仅维数小，而且是个三对角阵(稀疏矩阵)。
• 在选取\(k\)的值时，\(k\)应该尽可能的远大于\(r\)，通常求取前三个最大的特征值是比较准确的。如果\(r\)接近于\(m\)，那么Lanczos算法就不太准确了。
• 通过选择较大的Lanczos基向量的数量\(k\)，可以获得更多准确的最大的特征值，但是构造的Krylov队列越来越难满足列满秩的要求，因此无法得到正确的标准正交向量组\(U\)。其原因在于，\(k\)越大，\(A^k\)对应的更多的特征值\(\lambda_i^k\)会趋近于0，当\(0<\lambda_i<1\)。
• Lanczos算法只能估计最大的数个特征值，但是要想估计最小的数个特征值就做不到了，其适用性是比较局限的。

Lanczos算法求取最小的特征值的一些思考

​ 上一节介绍了利用Lanczos算法估计数个最大的特征值。如果要估计最小的数个特征值呢，不妨对矩阵\(A\)的逆进行操作。为了避免由于矩阵\(A\)是奇异的，可以进行正则化，即\(A+\delta I\)。仍然采用上述的Lanczos算法，希望得到

\[(A+\delta I)^{-1}\approx U*T*U^T。 \]

首先，需要构造Krylov队列\(\{u_1,(A+\delta I)^{-1}u_1,(A+\delta I)^{-2}u_1,\ldots,(A+\delta I)^{-k}u_1\}\)，显然求解\(A+\delta I\)的逆是非常消耗计算量的，可以采用LU分解，将高维矩阵\(A+\delta I\)分解成一个下三角矩阵与一个上三角矩阵的乘积。因此，可得

\[(A+\delta I)^{-1} = U^{-1}L^{-1} \]

下面的计算方法参照Lanczos算法求解步骤。

仿真例子

​ 随机生成一个\(50\times 50\)的对称阵\(A\)，求取其前3个最小特征值。令\(k=5\)，\(u_1=1_{50}\)，根据步骤1）计算得到\(U_{50\times 5}\)，这里设定\(\delta=1e^{-5}\)，计算\(T=U^T*(A+\delta I)^{-1}*U\)。对\(T\)进行svd分解，取前三个最大特征值。利用\(\hat{\lambda}_{\min} = \hat{\lambda}^{-1}-\delta\)得到估计的最小的特征。

估计的前五个最小的特征值	矩阵\(A\)的前五个最小的特征值
0.0001	0.0001
0.0011	0.0011
0.0032	0.0032
0.0159	0.0059
0.1728	0.0120

理论上似乎是可行的，但是实验结果并不理想。虽然提供仿真例子表明可以正确估计最小的数个特征值，但是Lanczos算法很不稳定，需要不断重启Krylov空间，来保证满足列满秩（这也是为什么\(k\)只能取5，取不到10的原因）。总结多次仿真结果，可以得到以下结论。

• \(\delta\)取值不能偏大，例如\(\delta=0.5\)，第一最小特征值为\(\lambda_1=1e^{-4}\)，第二最小特征值为\(\lambda_2=1e^{-3}\)，那么\((A+\delta I)^{-1}\)的前两个最大特征值分别为1.9996和1.996，他们之间差异性被大大缩小了，\(\lambda_2/\lambda_1 = 10\)。
• \(\delta\)取值不能过于偏小，例如\(\delta=1e{-8}\)，当\(A\)的最小特征值为0时，那么\((A+\delta I)^{-1}\)的最大特征值就是\(1e^8\)，\(A^5\)的最大特征值就是\(1e^{40}\)，因此\(k\)的取值就非常局限了，\(k\)过小会导致最小的数个特征值估计不准确。
• 当然为了防止计算溢出，可以进行比例缩放，例如\((A+\delta I)^{-1}/1e^{4}\)，计算出特征值后再等比例放大。这种方法测试下来是比较有效的，\(k\)值可以扩大到8，前三个最小的特征值估计是准确的。

参考资料：

lanczos算法及C++实现（一）框架及简单实现
Stewart, G.W. "A Krylov-Schur Algorithm for Large Eigenproblems." SIAM Journal of Matrix Analysis and Applications. Vol. 23, Issue 3, 2001, pp. 601–614.