[矩阵计算]Lanczos方法：求稀疏矩阵特征值
更新： 29 JUL 2016
由QR方法知，求矩阵$A$的特征值，大多需要先将其三对角化（详细方法见徐树方先生的教材。此处外链一个例子），即
$$ T=Q^TAQ $$
即找到正交矩阵$Q$使得$T$成为三对角矩阵。然而若$A$为大型稀疏矩阵，常用的方法如Householder和Givens变换都无法充分利用$A$的稀疏性，因此考虑直接计算$T$和$Q$的矩阵元以利用$A$的稀疏性加速运算。
一、Lanczos方法基本原理
将以上分解式中的$Q$写成
$$ Q=[q_1,q_2,\cdots,q_n] $$
其中$ q_i $为$Q$的列向量。$T$写成
$$ T=\begin{bmatrix} \alpha_1 & \beta_1 & & & 0 \\ \beta_1 & \alpha_2 &\ddots & & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \alpha_{n-1} & \beta_{n-1} \\ 0 & & & \beta_{n-1} & \alpha_n \end{bmatrix} $$
比较
$$ AQ=QT $$
两边矩阵的每一列，得到
$$ Aq_i = \beta_{i-1}q_{i-1}+a_iq_i+\beta_iq_{i+1}, \quad i=1,2,\cdots, n$$
由于$ \beta_0, \beta_n, q_0, q_{n+1} $未定义，补充定义$ \beta_0q_0 = \beta_nq_{n+1} = 0 $，这样上式两边乘$q_i^T$得到
$$ \alpha_i = q_i^TAq_i, \qquad \beta_i = q_{i+1}^TAq_i=||Aq_i-\beta_{i-1}q_{i-1}-\alpha_iq_i||_2$$
可以看出只要给定任意的$q_1 \in \mathbb{R}^n$且$||q_1||_2=1$，就能够利用递推关系得到全部$ q_i, \alpha_i, \beta_i $，迭代格式为
$ \alpha_1 = q_1^TAq_1, $ //初始值
$ r_i=Aq_i-\alpha_iq_i-\beta_{i-1}q_{i-1}, $
$ \beta_i=||r_i||_2, $ //由上一轮$\alpha$，$q$产生新的值
$ q_{i+1}=r_i/\beta_1 \ (\beta_i \neq 0)$
$ \alpha_{i+1} = q_{i+1}^TAq_{i+1}, \ i=1,2,\cdots, n-1$
此即Lanczos迭代。其中$q_i$称为Lanczos向量。在第$j$步产生的矩阵$T_j$称为j阶Lanczos矩阵，其特征值可能是$A$的部分特征值的很好的近似。详细内容参考Krylov子空间。
注意其中若某一步$\beta_i=0$，则此时得到的$T_j$的特征值将都是$A$的特征值。
二、具体算法
Lanczos算法求大型稀疏矩阵$A$特征值（三对角矩阵）
1. 输入$A \in S\mathbb{R}^{n\times n}, q_1 \in \mathbb{R}^n\ (||q_1||_2=1)$
2. $u_1:=Aq_1,\ j:=1$
3. $a_j:=q_j^Tu_j,\ r_j:=u_j-a_jq_j,\ \beta_j:=||r_j||_2$
4. 若$\beta_j=0$则结束，否则
$q_{j+1}:=r_j/\beta_j,\ u_{j+1}:=Aq_{j+1}-\beta_jq_j$
$j:=j+1$，转步3
这一算法的主要工作量集中在计算矩阵$A$与向量$v$的乘积$Av$上。在实际使用时，应根据$A$的具体特点，设计一个计算$Av$的子程序，使算法的运算量尽可能少。
三、Lanczos现象
如果Lanczos算法是在没有摄入误差的情况下执行的，则所得到的Lanczos向量$q_1,\cdots,q_j$是相互正交的，而且之多$n$步必然终止。但是，在误差出现的情况下，计算得到的Lanczos向量的正交性很快就会失去，有时甚至还是线性相关的。C.C.Paige指出失去正交性恰与近似特征值的精度提高有关。
在Lanczos矩阵$T_j$中，其特征值$\mu_j$只要与其他$\mu_i$不要太靠近，特征向量$||z_j||_2$就和1很接近。而迭代次数$k$越大（可以超过$n$），则$T_k$含有越多的$A$的较好地近似特征值。当$k$充分大时，$T_k$就含有$A$的所有相异的特征值，此即Lanczos现象。
粗略地讲，位于$A$谱区间两端且分离较好地特征值$\lambda$，在$k\ll n$时$T_k$的特征值内就含有$\lambda$的很好的近似值；位于区间内部而又与其他特征值分离得不好的特征值$\lambda$，需$k\gg n$，$T_k$的特征值中才会有$\lambda$的较好地近似值。
因此，当我们只需计算大型稀疏对称矩阵$A$的少数几个两端特征值时，通常只需迭代很少几步（$k\ll n$），$T_k$的两端特征值就是$A$两端特征值的很好的近似值。
一个直观的例子是Parllet的数据，当$n=10^4$时，取$k=300$就可以求出10个$A$的两端特征值和对应的特征向量的很好的近似值。
若求$A$全部的特征值时，一般$k$要远大于$n$。Cullum和Willoughby的经验是对于绝大多数矩阵来讲，只需$k\leqslant 3n$就可以求出其几乎全部的不同特征值达到机器精度的近似值。
显然，对于$\textbf{HC}=E\textbf{C}$问题的精确基态和几个低能态来说，Lanczos方法是有用而且可以进一步发展的。在这个问题上最高能级一般不要求求出，那么得到低能级的精确值的方法还可以进一步加速。具体的方法就是计算化学中常用的Davidson对角化。
