lanczos算法及C++实现（二）实对称阵奇异值分解的QR算法
1. 基本的QR算法
我们先讨论一般对阵矩阵的QR算法，再讨论对称三对角阵的QR算法
给定一个实对称阵X,假设其特征值分解为X=PSP'，其中P对正交阵，S是对角阵。求P,S的QR算法如下，其中 $Q_k$为正交阵，$R_k$为上三角阵：
$X_1=X$
for k=1,2, ...
$X_k=Q_k*R_k$ （QR分解）
$X_{k+1}=R_k*Q_k$
endfor
$P=Q_1*Q_2* \dots *Q_{\infty}$
将$R_{\infty}$的非对角线元素置0即得S
用matlab代码表示，即
function [P,S]=qreigen
A=rand(10);
A=A+A' %Symmetry real matrix, otherwise P'AP is a block upper triangle matrix
m=size(A,1);
P=eye(m);
for i=1:1000
[Q,R]=qr(A); %QR factorization
A=R*Q;
P=P*Q;
end
S=diag(diag(R));
end
这个算法很简单，每次迭代做一个QR分解和一个矩阵乘法。最基本的，我们用施密特正交化进行QR分解，这在上一篇lanczos算法及C++实现（一）框架及简单实现中使用。A讨论了QR方法和幂迭代方法的关系。
值得注意的是，如果矩阵X不是对称阵，这个方法往往不能收敛到一个上三角阵，而是一个块状上三角阵。
function [P,A]=qreigen A=rand(10); %not symmetry m=size(A,1); P=eye(m); for i=1:1000 [Q,R]=qr(A); %QR factorization A=R*Q; P=P*Q; end % A is block upper triangle matrix, not upper triangle A end
运行结果如下：
A =4.7120 -0.6211 0.4352 -0.0450 -0.1297 0.5700 -0.3769 0.4790 -0.0754 0.2798 0 0.0273 0.8587 -0.2188 -0.0095 -0.1913 -0.3785 -0.2142 -0.0565 -0.1236 0 -0.9150 0.3761 -0.0914 0.0217 -0.1527 0.0317 0.4235 0.4612 -0.3284 0 -0.0000 0.0000 -0.7400 -0.0686 0.2954 0.1561 -0.1755 0.1879 -0.2237 0 -0.0000 0.0000 -0.0000 0.4096 -0.5238 -0.0330 0.1636 0.3090 -0.4849 0 0.0000 0.0000 -0.0000 0.2002 0.7429 0.2064 -0.4820 -0.2602 -0.2005 0 0 0 0.0000 -0.0000 -0.0000 -0.3812 -0.1792 -0.0188 -0.4469 0 0 0 0.0000 -0.0000 0.0000 0.4117 -0.1513 0.1134 -0.2247 0 0 0 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.3465 0.0887 0 0 0 0 0 0 -0.0000 -0.0000 -0.0000 0.1833
2. 快速QR算法
施密特正交法每次需要O(n^3)的计算，很慢。所以通常不用这个方法。一般来说，我们采用如下QR算法，该算法分成两步：
1、利用Householder变换 (如果你对Householder变换不了解，请见附录B)，将X正交变换为上海森伯格阵(Upper Hessenburg Matrix)，复杂度O(n^3)。
2、在每次迭代中，利用Givens变换 (如果你对Givens变换不了解，请见附录C)，求得上海森伯格阵的QR分解，复杂度O(n^2)。
这样只有在迭代之前需要O(n^3)的预处理，而每次迭代只要平方复杂度，因此比施密特正交法快很多。
2.1 矩阵的上海森伯格化 Hessenburg Reduction
在这一步，对一个方阵X(不要求对称)，我们要求一个正交变换P，使得PXP'为上海森伯格矩阵，即低于对角线下一行的所有元素为0：
$PXP^T=\left(\begin{matrix}* & * & *& *& \dots \\ * & * & *& *& \dots \\ & * & *& *& \dots \\ & & *& *& \dots \\ & & & *& \dots \\ & & & & \dots \end{matrix}\right)$
如何求得这样的P呢？基本思想如下，将矩阵X分成4块：
$\left(\begin{matrix}x_{11} & X_{1,2:n}\\ X_{2:n,1} & X_{2:n,2:n} \end{matrix}\right)$
然后用Householder变换，将左下的一列变为只有一个非零元素：
$P_1 X_{2:n,1} = v_1 \mathbf e_1=\left(\begin{matrix}v_1\\ 0 \\ 0 \\ \dots \end{matrix}\right)$
其中$P_1=I-2u*u^T$是一个正交矩阵，
$u=\frac{X_{2:n,1}-\|X_{2:n,1}\|\mathbf e_1}{\left\|X_{2:n,1}-\|X_{2:n,1}\|\mathbf e_1\right\|}$ ，
$ \mathbf e_i$是一个除第i个分量为1以外，其他分量为0的向量。P1的具体求法，参见附录B
这样我们有
$\left(\begin{matrix}1 & 0 \\ 0 & P_1 \end{matrix}\right) X \left(\begin{matrix}1 & 0 \\ 0 & P^T_1 \end{matrix}\right) = \left(\begin{matrix}x_{11} & X^*_{1,2:n} \\ v_1\mathbf e_1 & X^*_{2:n,2:n} \end{matrix}\right) $
类似的，我们可以构造$\left(\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & P_2\end{matrix}\right)$等共n-1个矩阵，使得X正交相似于一个上海森伯格矩阵。
值得一提的是，如果X是一个对称阵，那么根据对称性，X正交相似于一个三对角阵。
2.2 上海森伯矩阵的QR分解
首先，如果H是一个上海森伯格阵，其QR分解为H=QR，那么RQ也是一个上海森伯格阵。证明省略。
这就是说，如果我们有一个对海森伯格阵的快速QR算法，那么在每次迭代中我们可以反复使用该算法。
我们采用Givens变换来进行QR分解。对于n阶上海森伯格阵H,从上到下对每相邻的两行使用Given变换，每次变换都消去对角线下方的一个元素。共n-1次变换。由于Givens变换是正交变换，那么这n-1次变换的乘积也是正交变换，变换后，
\(G^T_{n-1}*\dots *G^T_1 H\) 变成了上三角阵。其中\(G^T_i\)是第i次Givens行变换的矩阵表示。
根据QR算法，我们需要求RQ。因此在对H进行n-1次行变换后，再从左到右对每相邻的两列使用对应的Given变换的逆变换，同样一共也是n-1次变换：
\(RQ=G^T_{n-1}*\dots *G^T_1 H G_1*\dots *G_{n-1}\)，这样我们就完成了一次迭代。
这个过程可以用下图表示：
QR分解过程,$G_{n-1}^T*\dots*G_1^TH=R$，其中红框代表当前需要考虑的H中的元素，这些元素将被更新。注意每次Givens变换将一个对角线下的元素变成0.
而R*Q的过程如下：
其中
$c_i=\frac{H(i)_{i,i}}{\sqrt{H(i)^2_{i,i}+H(i)^2_{i+1,i}}}$
$s_i=\frac{H(i)_{i+1,i}}{\sqrt{H(i)^2_{i,i}+H(i)^2_{i+1,i}}}$
$H(i)=G^T_{i-1}*\dots*G_1^TH$
将这个过程中所有产生的G乘起来就得到了H的特征向量。我们可以看到，每次迭代只需要O(n^2)的计算量，比施密特正交化快多了。
2.3 三对角阵的特征值分解
在lanczos方法中，我们得到的是一个对称三对角阵，是上海森伯格阵的特例。因此可以用上面所述方法计算特征值。值得注意的是，这里可以利用三对角阵的特性加速计算，在求第一个行变换G'1H的时候，可以只考虑H左上的2*3的小阵，而不是整个前两行。而接下来的每次行变换也只需考虑H的2*3的子阵。在列变换中，每次只需考虑H的3*2的子阵即可。这样一次QR迭代的复杂度只有O(n)了。
2.4 通过平移(shifts)和deflation加速上海森伯格矩阵的特征值分解
由于QR算法的收敛速度取决于特征值之间的差异，差别越大，QR收敛越快。这和幂迭代算法是一致的。那么怎么加快特征值之间的差异呢？通常的做法就是平移。
假设H有两个特征值$\lambda_1>\lambda_2\ge0$，那么QR迭代的速度为 $\left(\frac{\lambda_2}{\lambda_1}\right)^k$。假如$\delta$满足$\lambda_1>\lambda_2\ge\delta\ge0$，那么矩阵$H^*=H-\delta I$和H有相同的特征向量，其特征值为$\lambda_1-\delta>\lambda_2-\delta\ge0$，而求H^*的特征向量会很快，其收敛速度为$\left(\frac{\lambda_2-\delta}{\lambda_1-\delta}\right)^k$
基于这样的观察，在求海森伯格矩阵的特征值时，每次迭代我们进行如下QR分解
$H-\delta I=QR$ （QR分解）
$H=RQ+\delta I$
显然，这样的QR分解不改变H的特征向量，但是如果选取合适的
$\delta$收敛速度要快很多。特别地，如果$\delta$恰好取到某个的特征值时，那么一次迭代即可就可确定该特征值。见文献https://www.math.kth.se/na/SF2524/matber15/qrmethod.pdf中的lemma 2.3.2。事实上我们没那么幸运，但可以用每次迭代中，H对角线上的元素可以看作是特征值的估计。所以每次计算完RQ后，取对角线上固定某个位置的元素（比如始终取第1个）作为$\delta$。但如果矩阵存在一对特征值互为相反数，那么这种方法不能收敛。比如矩阵 $\left(\begin{matrix}0 & 1\\1 & 0\end{matrix}\right)$，因此更好的做法是Wilkinson shift，取H对角线上的2*2的方阵的特征值作为$\delta$。此外对复矩阵还有双位移方法(Double shift)，这里不作详细的讨论，有兴趣的读者可以参考http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf
那么这样只能加速计算某个特征值，比如始终取对角线上的最后一个则加快了求解最小特征值的速度。如何加速求解其他特征值的速度呢？这里就用到了deflation。
抱歉我不知道怎么翻译这个词deflation，因此就用英文了。
Deflation是指在迭代过程中，如果对角线上某个位置的下一行几乎为0的时候，说明该位置上的值很接近于特征值了，此时将改位置所对应的行和列从矩阵中刨去，对剩下的矩阵进行计算。
一般的步骤如下，取对角线最右下的元素做deflation，直到该元素左边的元素几乎为0时，这时最小特征值已经求出，将最后一行一列去掉，对剩下的矩阵继续计算。如此反复，deflate n-1次之后，所有特征值全部求出。
2.5 QR算法的C++实现
最后，附上我的QR算法的C++实现。这个算法实现了这篇文章中的所有内容，包括实对称阵到上海森伯格阵的转化（hessenburgReduction函数），上海森伯格阵特征值分解的QR算法(QRHessenberg函数)，三对角阵的特征值分解的QR算法(QRTridiagonal函数)，每个QR算法中都用了shift和deflation加速，采用的是Wilkinson shift。在Lanczos算法中，只有QRTridiagonal函数被调用。当然如果不用lanczos方法三对角化，也可以通过hessenburgReduction函数和QRHessenberg函数来进行QR迭代求解特征值分解。具体代码参见lanczos算法及C++实现（〇）C++代码
3. 参考文献
https://www.math.kth.se/na/SF2524/matber15/qrmethod.pdf
http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf
附录
附录 A QR算法等价于并行幂迭代算法(Simultaneous Iteration)
QR算法实际上等价于并行幂迭代算法 (Simultaneous Iteration)。在幂迭代中，如果我们用一组线性独立的向量而不是一个向量与X反复相乘，则最终得到的将是一组线性独立的向量，对其正交化后，可以 得到X的特征向量。但实际上由于幂迭代的性质，这些向量快速的收敛到主特征向量，而其他的特征向量被快速丢弃。为了保留较小特征向量的成分，每次与X相乘 后，都做一次正交化，那么这样最终这些向量将收敛到X的不同的特征向量。如果一开始选取的一组向量是单位阵，那么这种并行的幂迭代算法，实际上就是QR算 法。为什么这么说呢？我们可以看到，第一迭代，这组向量相乘X，得到X，然后正交化。而正交化实际上就是QR分解：$X=Q_1*R_1$，接着第二次迭代，我们用X乘$Q_1$并正交化：$X*Q_1=Q_2*R_2$，如此我们得到
\(X=Q_1R_1\)
\(XQ_1=Q_2R_2\)
\(XQ_2=Q_3R_3\)
...
比较一下QR算法
\(X=Q_1R_1\)
\(R_1Q_1=Q_2R_2\)
\(R_2Q_2=Q_3R_3\)
...
用数学归纳法可以证明QR算法中的
\(R_k\)
和并行幂迭代中的
\(R_k\)
相同。因为当k=1时，显然成立，当k=2时
并行幂迭代中
\(Q_1R_1Q_1=Q_2R_2\)
而QR算法中
\(R_1Q_1=Q_2R_2\)
由于
\(Q_1\)
是正交阵，一个矩阵在进行行正交变换（左乘正交阵）后，其QR分解中的R不变。因此两者的
\(R_2\)相同
以此类推。可以证明两个算法中所有的
\(R_k\)相同
然后，我们可以得到在并行幂迭代中，
\(Q_k=X^k*R_k^{-1}*R_{k-1}^{-1}*\dots*R_1^{-1}\)
而QR算法中
\(Q_k=R_{k-1}*\dots*R_1*X*R_k^{-1}*R_{k-1}^{-1}*\dots*R_1^{-1}\)
于是，可以得到两个算法中
\(Q_k\)的联系：
\(Q^{\mathrm{Simultaneous Iteration}}_k=Q_1^{QR}*\dots*Q_k^{QR}\)
附录 B Householder变换
又名Householder Reflections，顾名思义，是一种镜像变换。如图，假设给定一个点x和一个过原点的平面（虚线）,平面用一个垂直于它的单位向量v表示，那么对x的Householder变换就是求x关于平面v的镜像点x'。现在我们来看怎么计算x'
首先$0$指向$x$的向量在向量$v$上的投影为 $v*v^Tx$ ，在平面v上的投影为$x-v*v^Tx$。所以其镜像点x'的坐标为$x-2v*v^Tx = (I-2vv^T)x$
这里矩阵$I-2vv^T$就是一个householder矩阵，是一个对称的正交矩阵，有两个特征值，1和-1
现在我们考虑如何利用householder变换将一个向量变x成只有一个非0元素的向量,如下图所示
关键就是求那个平面，即v向量。由图可得出x-x'平行于v，因此我们有
$v=\frac{u}{\|u\|}$
其中
$u=x-\left(\begin{matrix} \|x\| \\ 0 \\0 \\0 \\ \dots\end{matrix}\right)$
附录 C Givens变换
又名Givens旋转，顾名思义，是对一个向量进行旋转变换，不改变向量的长度，因此是一个正交变换，其逆变换就是反方向旋转回去。
如图所示，我们考虑一个二维向量$v=\left(\begin{matrix}a\\b\end{matrix}\right)$，将它逆时针旋转$\theta$角度后，求所得到的坐标。
这里我们可以将向量重新表示为\(v=|v|\left(\begin{matrix}cos(\alpha)\\sin(\alpha)\end{matrix}\right)\)，其中\(\alpha\)是向量v与x坐标轴的夹角。
这样变换后的坐标为
$w=|v|\left(\begin{matrix}cos(\alpha+\theta)\\sin(\alpha+\theta)\end{matrix}\right)$
$=|v|\left(\begin{matrix}cos(\alpha)cos(\theta)-sin(\alpha)sin(\theta)\\sin(\alpha)cos(\theta) + sin(\theta)cos(\alpha)\end{matrix}\right)$
$= \left(\begin{matrix}cos(\theta) & -sin(\theta)\\sin(\theta) & cos(\theta)\end{matrix}\right)*|v|\left(\begin{matrix}cos(\alpha)\\sin(\alpha)\end{matrix}\right) $
$= \left(\begin{matrix}cos(\theta) & -sin(\theta)\\sin(\theta) & cos(\theta)\end{matrix}\right)v$
其中$G= \left(\begin{matrix}cos(\theta) & -sin(\theta)\\sin(\theta) & cos(\theta)\end{matrix}\right)$为Givens变换矩阵，是一个正交阵，其逆阵为
$G^T= \left(\begin{matrix}cos(\theta) & sin(\theta)\\-sin(\theta) & cos(\theta)\end{matrix}\right)$
本文属作者原创，转载请注明出处：
http://www.cnblogs.com/qxred/p/qralgorithm.html
本系列目录：
lanczos算法及C++实现（二）实对称阵奇异值分解的QR算法
lanczos算法及C++实现（三）实对称三对角阵特征值分解的分治算法
