lanczos 算法及C++实现（⼀）框架及简单实现
1. lanczos ⽅法的⼤致思路
为了求m 阶⽅阵X 最⼤的r 个特征值和特征向量： X m ×m ≈U m ×r S r ×r U T m ×r ,
其中U 是列正交矩阵，即 U T U =I ，每⼀列为⼀个特征向量，S 是对⾓阵，对⾓线上每个元素为特征值。r 为分解的秩
lanczos 算法分三步求解:
1） 对X 进⾏正交变换得到⼀个三对⾓阵T ：X =PTP T ，其中P 为正交阵
T =α0β10000…β1α1β2000
…0β2α3β300…00β3α4β40
…000β4α5β5
……………………    2） 对T 进⾏奇异值分解：T =W m ×r SW T m ×r
3） 最后 X ≈PWS (PW )T 即为所求
2. lanczos ⽅法的优点
1） T 矩阵是⼀个三对⾓阵，很稀疏，因此对它的特征值分解会⾮常快。
2） 如果r 很⼩，可以不⽤求整个T ，⽽是求其左上1.5r 阶的⽅阵即可得到很好的近似。这样对T 的特征值分解会更快。() 另外⼀种⽅法动态决定求多少lanczos 向量， 这⾥简单总结⼀下，初始选取n ≤m 个lanczos 向量，然后求T 的最⼤特征值和最⼩特征值，然后逐渐增加n ，并更新这两个特征值，直到这种更新⾮常⼩为⽌。如果⽤在求r 个主特征向量，那么可以更新第1个和第r 个特征值。
3. lanczos ⽅法的导出
考虑幂迭代产⽣的⼀系列向量P =[v ,Xv ,X 2v ,…X m −1v ]，其中v 是任意向量。假设X 的特征值分解为X =∑i λi ξi ξT i ，并且特征值按绝对值从⼤
到⼩排序 |λ1|≥|λ2|≥…。则 Xv =∑i λi ξi ξT i v =∑i λi v ′i ξi ，其中v ′i =ξT i v 是v 在正交基[ξ1,ξ2…]下的坐标。 类似的，我们可以得到
X k v =∑i λk i v ′i ξi 。也就是说特征值⼤的特征向量在X k v 中⽐重越⼤(这也是幂迭代收敛到主特征值的原因) 。换句话说，P 中的每列X k v 都是成特征向量ξ1,ξ2,…的线性组合，这些特征向量构成了 P  的⼀组正交基，并且特征值越⼤的特征向量所占的⽐例越⾼。⼀个⾃然的想法就是，对P 的前l 个向量所构成
的矩阵进⾏特征值分解，应该能得到X 最⼤的⼏个特征向量的近似解，l 越⼤，越近似。
然⽽，这种原始的⽅法并不稳定。因为较⼩的特征值被快速的抛弃，即使是第⼆⼤特征值，也是随l 增⼤，其⽐重指数速度缩⼩(|λl 1|>>|λl 2|)。因此，为了能够保留⼩特征向量，需要作正交化，即P 是列正交的。由此便产⽣了 使⽤施密特正交化：每次迭代所得向量都要去掉之前所有向量的投影：p k +1=Xp k ⊖p k ⊖p k −1⊖p k −2…
其中a ⊖b  表⽰ a 在b 的正交空间的投影：
正交化以后，虽然X 特征向量仍旧是P 的⼀组正交基，但⼤特征向量在P 中的⽐重就不⼀定⾼了。因此P 的最⼤特征向量
并不是X 特征向量的近似。 注意到如果对P 每列归⼀化: 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}\|}论⽂解释了该⽅法实际上就是⽤迭代后的向量和前两次向量做正交，但是并没有证明这样得到的向量和其他向量⼀定正交和为什么得到的是三对⾓阵，本⽂在附录A中，给出了lanczos⽅法的正交以及三对⾓性的证明。
4. 三对⾓阵T的特征值分解
三对⾓化后的矩阵T有三种常⽤的⽅法进⾏分解：
1） QR法，即对T进⾏QR分解，T=QR，其中Q为正交阵，R为上三⾓阵，然后令T'=RQ，如此反复进
⾏，直⾄收敛。其中产⽣的Q矩阵的连乘就得到了T的特征向量。具体算法及实现请参考
⽂章第⼆页给出了伪代码。在2.2.2节中，还提到了针对上海森伯格阵的快速QR算法，正好适⽤于这⾥的T矩阵。这⾥不再赘述。相⽐于另两种⽅法，QR⽅法并不⾼效，但是在求最⼤的少数特征向量时，T通常很⼩，因此为实现简单，以及稳定性考虑，QR⽅法常被使⽤。
2）分治法。这个⽅法将T划均匀分成两块三对⾓阵和⼀个秩为1的2*2块状⽅阵的和。分别对两个三对⾓阵分解，然后合并。合并时对特征向量作秩1的修正。详细讨论见下⼀篇⽂章。参考⽂献以及wiki 中的介绍。还有对其中secular equation解法的⽂献
3） MRRR⽅法。这个⽅法是理论上最快的⽅法，O(l^2)复杂度，主要采⽤逆迭代来加快收敛。详细讨论见后续⽂章。该⽅法作者的博⼠论⽂并行计算框架
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向量采⽤的是
中的伪代码实现。由于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
最后⼀个等号由正交性得出。证毕
本⽂属原创，转载请注明出处
⽬录：
Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。