從一個特徵值到全部特徵值
上一篇裡，冪法靠著反覆用 A 去乘一個向量、看著它對齊主特徵向量，給了你那個唯一最大的特徵值；而反冪法則讓你能聚焦到任何一個你已經猜得出來的特徵值。兩者都珍貴，但兩者都是「零售」：它們一次只交出一個特徵值，要拿下一個就得對矩陣做降階（deflation）再重來一遍。對一個稠密矩陣、而你想要它全部的特徵值時——這在振動分析、量子化學、穩定性研究裡是家常便飯的需求——那種零售做法就笨拙了。我們要一個「批發」的方法。
而我們絕不能被教科書的定義誘惑。特徵值是特徵多項式 det(A - lambda I) = 0 的根，那為什麼不把那個行列式展開、直接解多項式就好？因為，正如這一級第一篇所論證的，那條路是陷阱：一個多項式的根對它的係數可以極度敏感，所以一旦組出特徵多項式，你就丟掉了再也找不回來的精度。誠實的路是直接從矩陣計算特徵值，絕不經過它的特徵多項式。QR 演算法正是那條路。
那個不講道理的把戲：分解，然後對調
整個想法就在這裡，而它看起來簡單到幾乎不可能管用。從 A_0 = A 開始。把它分解成 A_0 = Q_0 R_0。現在用「相反的順序」把這兩塊乘回去，組出下一個矩陣：A_1 = R_0 Q_0。再把 A_1 = Q_1 R_1 分解、再對調得到 A_2 = R_1 Q_1，如此繼續。就這樣。反覆地：分解成 Q R，再以 R Q 重組。
the bare QR iteration:
A_0 = A
repeat for k = 0, 1, 2, ... :
Q_k, R_k = qr(A_k) # factor: A_k = Q_k R_k
A_{k+1} = R_k Q_k # swap the order back
key fact: A_{k+1} = R_k Q_k = Q_k^T (Q_k R_k) Q_k = Q_k^T A_k Q_k
so each step is an ORTHOGONAL SIMILARITY transform of A_k --
same eigenvalues at every step, just a friendlier basis.
as k grows, A_k drifts toward upper-triangular form,
and its DIAGONAL entries converge to the eigenvalues of A.為什麼這個沉悶的小小對調，會收斂到特徵值？兩個事實扣在一起。第一，因為 Q_k 是正交的，A_{k+1} = Q_k^T A_k Q_k 是一個相似變換——A_{k+1} 跟 A_k 有著一模一樣的特徵值，因此也跟原本的 A 相同。所以沒有任何一步會改變答案；迭代只改變你看這個矩陣所用的基底。第二，也更微妙：這個特定的對調，其實偷偷地在「所有方向上同時跑冪法」：累積起來的那些 Q 的正交行向量，本身正收斂到特徵向量，並按特徵值的大小排好序。結果就是 A_k 慢慢攤平成上三角的形狀，而一個三角矩陣會把它的特徵值大方地擺在對角線上。
一個小小的心像：想像 A 是一個傾斜了某個角度的對稱 2x2 矩陣。每一個 QR 步驟，都把你的座標軸朝著矩陣的「天然軸」——也就是它的特徵向量——稍微轉一點。一直轉下去，非對角的元素就朝零縮小；當它們消失時，對角線就握著特徵值。QR 演算法在精神上，是一台耐心的機器，一次一個分解，去尋找那個對的旋轉。
讓它快到能用：先化成上海森堡形
上面那個赤裸的迴圈很美，卻不能用。一個稠密 n x n 矩陣的每一次 QR 分解要花 O(n^3) 的工作，而你可能需要幾十次迭代，所以天真地跑下去是 O(n^4) 甚至更糟——貴得離譜。把 QR 從一個趣味玩物變成主力的那個修正，是一個一次性的準備步驟：在迭代之前，先把 A 化約成上海森堡形（upper Hessenberg form）——上三角，再加上主對角線正下方那一條額外的非零次對角線。這透過一連串正交反射、一次到位地（差個捨入）完成，而它本身是一個相似變換，所以對特徵值毫無更動。
為什麼這一個準備步驟回報這麼巨大？一個近乎魔法的結構事實：QR 步驟會保持海森堡形。如果 A_k 是海森堡形，那 A_{k+1} 也是，免費。而分解一個海森堡矩陣很便宜——只有那一條次對角線需要被清成零，用一輪吉文斯旋轉（或豪斯霍爾德反射）就能在 O(n^2) 的工作量內辦到，而不是 O(n^3)。一次 O(n^3) 的前置，接著每一步都只要 O(n^2) 的迭代：這就是玩具與工具之間的差別。對於對稱的情形，回報更加銳利，下一節會說明。
- 用正交相似變換把 A 一次化成上海森堡形（O(n^3)，只做一次，差個捨入即為精確）——特徵值不變。
- 對海森堡矩陣施行 QR 步驟；因為海森堡結構被保持，每一步只花 O(n^2)。
- 盯著某個次對角線元素朝零縮小；當它小到可忽略時，角落裡就有一個特徵值收斂了——把它分離出去（降階），再在更小的子塊上繼續。
- 反覆到每一個次對角線元素都消失為止；此時對角線就握著全部的特徵值。
位移與降階：速度真正的來源
海森堡形讓每一步便宜，但純 QR 仍可能需要很多步，因為就像它底層的冪法一樣，它的收斂速率由特徵值的比值掌管——如果兩個特徵值大小相近，相關的那個非對角元素就會龜速地爬向零。解藥跟當初讓反冪法脫胎換骨的那一招相同：位移（shift）。不要分解 A_k，而是分解經過位移的矩陣 A_k - mu_k I = Q_k R_k，然後設 A_{k+1} = R_k Q_k + mu_k I（位移被加回去，好讓特徵值留在原地）。一個選得好、落在某特徵值附近的位移 mu_k，會讓底部角落的次對角線元素猛地下墜——在對稱情形下，收斂典型地變成三次的（cubic），意思是一旦靠近，每一步正確位數大致翻三倍。
在還不知道特徵值的情況下，你怎麼挑一個好位移？右下角那個 2x2 的尾端子塊，本身就是一個小矩陣，它的特徵值你用手就能算；取其中一個（威爾金森位移，Wilkinson shift）就能對「即將收斂的那個特徵值」給出極佳的估計。一旦角落的次對角線元素掉到容忍度以下，右下角那個對角元素就被鎖定為一個已收斂的特徵值。接著你就降階：把那一列和那一行剝掉，在剩下的更小子塊上跑迭代。這正是你在冪法那裡見過的位移與降階節奏，如今被工業化了——讓一個角落收斂、鎖住它、縮小、重複，直到矩陣完全變成三角形。
對稱情形、誠實面，以及它的定位
當 A 是對稱的時，一切都變好，這就是為什麼對稱特徵值問題值得擁有自己的名字。一個對稱矩陣的海森堡化約會產生一個三對角矩陣——只有對角線和它兩側相鄰的兩條對角線上有非零元素——因為對稱性也強迫第一上次對角線之上的元素都消失。QR 迭代於是在三對角的形狀上運作，每步 O(n)，配上威爾金森位移就三次收斂。特徵值保證是實數、特徵向量正交，而整個計算是逆向穩定的。用這種方式算出一個對稱矩陣的全部特徵值，總成本大約 O(n^3)——真正付得起，是那個你會不假思索地呼叫上千次的常式。
現在來談這一級堅持的那份誠實。QR 演算法回傳的每一個特徵值都是近似值——迭代是在某個次對角線元素低於門檻時停下，而不是在它恰好為零時，因為在浮點數裡它基本上永遠到不了恰好的零。好消息是 QR 是逆向穩定的：你得到的特徵值，是某個與 A 非常接近的矩陣的精確特徵值。但逆向穩定只是精度故事的一半。另一半是條件性：如果 A 的特徵值本身就是病態的——對矩陣的擾動高度敏感，例如接近虧損結構的非對稱矩陣就會這樣——那麼即使是一個完美的逆向穩定結果，也可能離真正的特徵值很遠。精度等於條件性乘以穩定性，這裡也跟別處一樣。對稱矩陣是那個幸運的例外：它們的特徵值是完美條件的，所以逆向穩定就替你買到了完整的精度。
最後，談談 QR 所屬的疆界。對於稠密矩陣、邊長到幾千、而你想要全部特徵值的情形，它是對的工具——因為它花 O(n^3) 的工作與 O(n^2) 的記憶體去碰整個矩陣。至於迭代求解器那一級裡那些巨大的稀疏矩陣，你只想要少數幾個極端特徵值、又根本付不起儲存一個稠密海森堡形的代價，QR 就是錯的形狀了，這時換Arnoldi 與 Lanczos 方法上場——它們是 Krylov 子空間的表親，跟冪法一樣，只倚賴那個便宜的乘積 A x。把這個分工記在心裡：稠密又要全部就用 QR，稀疏又只要少數就用 Krylov。如今稠密特徵值已在手，下一篇要轉向奇異值分解——恰如其分地，它正是靠跑一個與這個演算法極近的親戚算出來的。