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++實現（三）實對稱三對角陣特征值分解的分治算法