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
目錄：