lanczos算法及C++实现（一）框架及简单实现
1. lanczos方法的大致思路
为了求$m$阶方阵$X$最大的$r$个特征值和特征向量： $X_{m\times m}\approx U_{m\times r} S_{r\times r} U^T_{m\times r}$,
其中$U$是列正交矩阵，即 $U^TU=I$，每一列为一个特征向量，$S$是对角阵，对角线上每个元素为特征值。$r$为分解的秩
lanczos算法分三步求解:
1） 对$X$进行正交变换得到一个三对角阵$T$：$X=PTP^T$，其中$P$为正交阵
$T= \left[ \begin{matrix} \alpha_0 & \beta_1 & 0 & 0 & 0& 0& \dots \\ \beta_1 & \alpha_1 & \beta_2 & 0 & 0& 0& \dots \\ 0 & \beta_2 &\alpha_3 &\beta_3 & 0& 0 & \dots \\ 0 & 0 & \beta_3 &\alpha_4 &\beta_4 & 0& \dots \\ 0 & 0 & 0 & \beta_4 &\alpha_5 &\beta_5 & \dots \\ \dots & \dots &\dots &\dots &\dots &\dots &\dots\\ \end{matrix} \right]$
2） 对$T$进行奇异值分解：$T=W_{m\times r}SW^T_{m \times r}$
3） 最后 $X\approx PW S (PW)^T$即为所求
2. lanczos方法的优点
1） $T$矩阵是一个三对角阵，很稀疏，因此对它的特征值分解会非常快。
2） 如果$r$很小，可以不用求整个T，而是求其左上$1.5r$阶的方阵即可得到很好的近似。这样对T的特征值分解会更快。(https://en.wikipedia.org/wiki/Lanczos_algorithm) 另外一种方法动态决定求多少lanczos向量，http://peghoty.blog.163.com/blog/static/4934640920071179393238/ 这里简单总结一下，初始选取$n \le m$个lanczos向量，然后求T的最大特征值和最小特征值，然后逐渐增加$n$，并更新这两个特征值，直到这种更新非常小为止。如果用在求$r$个主特征向量，那么可以更新第$1$个和第$r$个特征值。
3. lanczos方法的导出
考虑幂迭代产生的一系列向量$P = [v, Xv, X^2v, \dots X^{m-1}v]$，其中$v$是任意向量。假设$X$的特征值分解为$X = \sum_i \lambda_i \xi_i\xi^T_i$，并且特征值按绝对值从大到小排序 $|\lambda_1| \ge |\lambda_2| \ge \dots$。则 $Xv = \sum_i \lambda_i \xi_i \xi_i^Tv = \sum_i \lambda_i v'_i \xi_i$，其中$v'_i = \xi_i^Tv$是$v$在正交基$[\xi_1, \xi_2\dots ]$下的坐标。 类似的，我们可以得到$X^kv = \sum_i \lambda^k_i v'_i \xi_i$。也就是说特征值大的特征向量在$X^kv$中比重越大(这也是幂迭代收敛到主特征值的原因) 。换句话说，$P$中的每列$X^kv$都是成特征向量$\xi_1, \xi_2, \dots$的线性组合，这些特征向量构成了 $P$ 的一组正交基，并且特征值越大的特征向量所占的比例越高。一个自然的想法就是，对$P$的前$l$个向量所构成的矩阵进行特征值分解，应该能得到$X$最大的几个特征向量的近似解，$l$越大，越近似。
然而，这种原始的方法并不稳定。因为较小的特征值被快速的抛弃，即使是第二大特征值，也是随$l$增大，其比重指数速度缩小($|\lambda_1^l| >> |\lambda_2^l|$)。因此，为了能够保留小特征向量，需要作正交化，即$P$是列正交的。由此便产生了Arnoldi's algorithm 使用施密特正交化：每次迭代所得向量都要去掉之前所有向量的投影：$p_{k+1}=Xp_k \ominus p_k \ominus p_{k-1} \ominus p_{k-2} \dots$
其中$a \ominus b$ 表示 a在b的正交空间的投影：
正交化以后，虽然$X$特征向量仍旧是$P$的一组正交基，但大特征向量在$P$中的比重就不一定高了。因此$P$的最大特征向量并不是$X$特征向量的近似。
注意到如果对$P$每列归一化: $p_k = \frac{p_k}{\|p_k\|}$，则$P$就是一个正交阵，因此$P^TXP$和$X$有着相同的特征值，这意味着如果能求得$P^TXP$的主特征向量$w_1, w_2 \dots$，则$Pw_1, Pw_2, \dots$就是$X$的主特征向量。很重要的是$P^TXP$是三对角阵! 证明可参见后文。而求解三对角阵的特征向量会容易很多！
不过Arnoldi's algorithm使用的是施密特正交化，要求$O(l^2)$次向量内积，相减。随着$l$增大，该方法速度堪忧。 lanczos提出了更快的$O(l)$复杂度的正交化方法。由于$X$是实对称阵，可以证明，每次迭代得到的向量只需要和前两次得到的向量正交后，即可保证该向量和其他所有向量正交。在列归一化后，满足$X=PTP^T$, $T$为三对角阵。 $p_{k+1}=\frac{Xp_k \ominus p_k \ominus p_{k-1}}{\|Xp_k \ominus p_k \ominus p_{k-1}\|}$ 论文 http://www.cs.cornell.edu/~bindel/class/cs6210-f09/lec32.pdf 解释了该方法实际上就是用迭代后的向量和前两次向量做正交，但是并没有证明这样得到的向量和其他向量一定正交和为什么得到的是三对角阵，本文在附录A中，给出了lanczos方法的正交以及三对角性的证明。
4. 三对角阵$T$的特征值分解
三对角化后的矩阵T有三种常用的方法进行分解：
1） QR法，即对T进行QR分解，$T=QR$，其中$Q$为正交阵，$R$为上三角阵，然后令$T'=RQ$，如此反复进行，直至收敛。其中产生的$Q$矩阵的连乘就得到了$T$的特征向量。具体算法及实现请参考http://www.cnblogs.com/qxred/p/qralgorithm.html
文章https://www.math.kth.se/na/SF2524/matber15/qrmethod.pdf第二页给出了伪代码。在2.2.2节中，还提到了针对上海森伯格阵的快速QR算法，正好适用于这里的$T$矩阵。这里不再赘述。相比于另两种方法，QR方法并不高效，但是在求最大的少数特征向量时，$T$通常很小，因此为实现简单，以及稳定性考虑，QR方法常被使用。
2） 分治法。这个方法将$T$划均匀分成两块三对角阵和一个秩为1的$2*2$块状方阵的和。分别对两个三对角阵分解，然后合并。合并时对特征向量作秩$1$的修正。详细讨论见下一篇文章。参考文献 A Serial Implementation of Cuppen's Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem 以及wiki 中的介绍 Divide-and-conquer eigenvalue algorithm 。还有对其中secular equation解法的文献 Chapter 12 Solving secular equations
3） MRRR方法。这个方法是理论上最快的方法，$O(l^2)$复杂度，主要采用逆迭代来加快收敛。详细讨论见后续文章。该方法作者的博士论文 A New O(n2) Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector Problem
5. lanczos方法的实现
这里讨论大型、稀疏、非方阵A的奇异值分解的实现，并求其少数主奇异向量。主要从如下角度考虑：
1）对于非方阵的SVD可以转化为对称方阵，本文采用了通过对$X=A*A^T$的特征值分解从而得到$A$的奇异值分解的方法。$X$的特征向量即为$A$的左奇异向量，而$A$的右奇异向量可以通过左奇异向量左乘$A$得到。
2）由于矩阵非常大，不能完全读入内存，因此是存在硬盘上的。在计算lanczos向量时，每次计算$A*A^T*p$需要扫两遍磁盘。$A$在磁盘上顺序存储，顺序读取，这样可以利用磁盘的读缓存。
3）考虑到多线程并行计算$A*A^T*p$，需要开辟读缓存，每次装载一批数据读入内存。每个线程维护自己的一个$p$向量，迭代完之后合并。
4）由于只求少数主奇异向量，因此lanczos向量也不会很多，因此可以假设lanczos向量可以完全放在内存中。
本文中求lanczos向量采用的是
https://en.wikipedia.org/wiki/Lanczos_algorithm
中的伪代码实现。由于lanczos方法并不稳定，迭代很多次之后，所得向量可能不正交甚至线性相关，关于lanczos的重启即实现会在后续文章中讨论。
附录
附录A：由lanczos方法得到的$P$矩阵是正交阵，$P^TXP$是三对角阵。
首先引入一个符号，对于两个向量$a,b$, $a \oplus b=\lambda a + \theta b$，表示由$a$和$b$张成的空间中的某一个向量
证明：采用数学归纳法。需要证明两个命题
（1）正交性： $\forall i<k, p_i^Tp_k=0$
（2） 三对角：$\forall i<k-1, p_i^TXp_k=0$
先证明正交性。当$k\le2$时，结论显然成立。当$k>2$时，假设结论成立，即对任意的$i < j \le k, p_i^Tp_j=0$，现在要证明对所有$i \le k, p_i^Tp_{k+1}=0$
分两种情况，如果$i < k - 1$，则有$\begin{align}\nonumber & p_{k+1}^Tp_i \\\nonumber =& (Xp_k \ominus p_k \ominus p_{k-1})^Tp_i \\\nonumber =& (Xp_k)^Tp_i \\\nonumber =& p_k^TXp_i \\\nonumber =& p_k^T(p_{i+1}\oplus p_i \oplus p_{i-1})\\\nonumber =& p_k^Tp_{i+1} \\\nonumber = & 0\end{align}$
否则如果$i = k - 1$或者$i = k$,根据$p_{k+1} =Xp_k \ominus p_{k} \ominus p_{k-1}$，直接可以得到$p_i^Tp_{k+1}=0$
再证明三对角性。当$k\le 2$时，结论显然成立，当$k> 2$时，假设结论成立，则当$k+1$时，对于$i+1 < k+1$的任意$i$有
$ p_{k+1}^TXp_i = p_{k+1}(p_{i+1}\oplus p_i \oplus p_{i-1}) =0$
最后一个等号由正交性得出。证毕
本文属原创，转载请注明出处
http://www.cnblogs.com/qxred/p/lanczos.html
目录：
