6 矩陣特征值的數值計算

6 矩陣特征值的數值計算

6.1 特征值與特征向量

設\(A\)是\(n\)階矩陣，\(x\)是非零列向量，如果存在數\(\lambda\)滿足

\[Ax=\lambda x \]

那么稱\(\lambda\)是矩陣\(A\)的一個特征值，\(x\)則是屬於\(\lambda\)的一個特征向量。

理論上要求解矩陣\(A\)的所有特征值需要求解特征多項式\(f(\lambda)=\det(\lambda I-A)\)的所有根，而\(f(\lambda)\)是一個\(n\)階多項式，因為5次以上代數多項式是沒有求根公式的，因此直接求解特征多項式的難度過大。於是，就發展出了一系列的數值算法來計算矩陣的特征值。

格施哥林(Gerschgorin)圓盤定理 設\(A=[a_{ij}]\)是一個\(n\)階矩陣，則\(A\)的每一個特征值必定屬於下述某個圓盤之中：

\[|\lambda-a_{ii}|\le\sum_{k=1\\k\ne i}^n|a_{ik}|\quad(i=1,2,\cdots,n) \]

通過上述圓盤定理，我們可以估計矩陣特征值的大致范圍。

6.2 冪法

冪法是通過矩陣特征向量來求出特征值的一種迭代法，主要用來求解矩陣按模最大的特征值（稱為主特征值）及其對應的特征向量。

假設n階矩陣A的特征值為\(\lambda_1,\lambda_2,\cdots,\lambda_n\)滿足

\[|\lambda_1|>|\lambda_2|\ge\cdots\ge|\lambda_n| \]

並且假設A具有n個線性無關的特征向量\(x_1,x_2,\cdots,x_n\)，需要注意的是這里只考慮了主特征值唯一的情況。冪法的目標就是求解\(\lambda_1,x_1\)。

冪法的基本思想是從任意選取的初始向量\(V^{(0)}\)出發，構造如下向量序列

\[V^{(0)}\ne0,V^{(1)}=AV^{(0)},\cdots,V^{(k+1)}=AV^{(k)},\cdots \]

因為A的特征向量構成了\(\mathbb{R}^n\)空間中的一組基，所以\(V^{(0)}\)可以表示為

\[V^{(0)}=a_1x_1+a_2x_2+\cdots+a_nx_n \]

這里假設\(a_1\ne0\)，於是

\[\begin{aligned} V^{(k)}&=AV^{(k-1)}=\cdots=A^{k}V^{(0)}\\ &=a_1A^kx_1+a_2A^kx_2+\cdots+a_nA^kx_n\\ &=a_1\lambda_1^kx_1+a_2\lambda_2^kx_2+\cdots+a_n\lambda_n^kx_n\\ &=\lambda_1^k(a_1x_1+a_2\left(\frac{\lambda_2}{\lambda_1}\right)^kx_2+\cdots+ a_n\left(\frac{\lambda_n}{\lambda_1}\right)^kx_n) \end{aligned} \]

當\(k\)充分大的時候，我們有

\[\left(\frac{\lambda_i}{\lambda_1}\right)^k\approx0,i=2,3,\cdots,n \]

於是，我們得到\(V^{(k+1)}\approx a_1\lambda_1^{k+1}x_1,V^{(k)}\approx a_1\lambda_1^{k}x_1\)。此時這兩個向量的第\(j\)個分量，有如下關系：

\[V^{(k+1)}_j\approx{a_1\lambda_1^{k+1}[x_1]_j}=\lambda_1(a_1\lambda^k_1[x_1]_j) \approx \lambda_1V_j^{(k)} \]

所以

\[\lambda_1\approx\frac{V^{(k+1)}_j}{V^{(k)}_j} \]

並且此時\(V^{(k)}\)就可以作為屬於特征值\(\lambda_1\)的一個近似特征向量，這是因為特征向量乘以任意常數仍然都是當前所屬特征值的一個特征向量。

上述分析已經說明了冪法的本質，但如果直接使用上述方法則會出現嚴重的計算溢出問題。因為，如果\(|\lambda_1|>1\)，就會有\(\lim_{k\rightarrow\infty}\lambda_1^k=\infty\)。如果\(|\lambda_1|<1\)，則會有\(\lim_{k\rightarrow\infty}\lambda_1^k=0\)。於是在計算過程中會發生上溢或下溢。因此實際應用時，應該采用改進過后的歸一化冪法，算法如下：

\[\left\{\begin{aligned} &j=\mathop{\arg\max}_i|V_i^{(k)}|\\ &m_k=V_j^{(k)}\\ &U^{(k)}=V^{(k)}/m_k\\ &V^{(k+1)}=AU^{(k)} \end{aligned}\right.\quad k=0,1,2,\cdots \]

其中\(m_k\)表示第\(k\)次迭代過程中\(V^{(k)}\)的按模最大的分量。當\(|m_k-m_{k-1}|<\varepsilon\)或達到最大迭代上限時，算法停止並返回\(m_k\)作為主特征值，\(U^{(k)}\)作為特征向量。

算法分析：使用歸一化冪法，當\(k\)充分大的時候有

\[U^{(k-1)}\approx\frac{a_1\lambda_1^{k-1}}{M_{k-1}}x_1 \]

其中\(M_{k-1}=m_0m_1\cdots m_{k-1}\)。由於\(U^{(k-1)}\)的每一個分量的絕對值都是小於等於1，並且絕對值最大的分量一定等於1。因為

\[V^{(k)}=AU^{(k-1)}\approx\frac{a_1\lambda_1^{k}}{M_{k-1}}x_1\approx\lambda_1U^{(k-1)} \]

所以\(V^{(k)}\)的按模最大的分量就約等於\(\lambda_1\)，這意味着\(m_k\approx\lambda_1\)。

綜上可知，當\(k\)充分大的時候，可以將\(m_k\)作為主特征值的近似值，\(U^{(k)}\)則是屬於主特征值的近似特征向量。

上述分析對於主特征值是重根的情況也是成立的，例如

\[\lambda_1=\lambda_2=\cdots=\lambda_r \]

並且滿足

\[|\lambda_1|=\cdots=|\lambda_r|>|\lambda_{r+1}|\ge\cdots\ge|\lambda_n| \]

那么此時\(x_1,x_2,\cdots,x_r\)均是是屬於\(\lambda_1\)的特征向量，只要此時\(x_1,x_2,\cdots,x_n\)仍然構成\(\mathbb{R}^n\)空間中的一組基，那么冪法仍然是可用的。

冪法的幾個問題：

1. 在冪法中如果初始向量選取不當，可能會出現\(a_1=0\)的情況，但是考慮到計算過程中的舍入誤差的影響，迭代到某一步驟的時候，在\(x_1\)方向的分量就不為零了。因此，實際初始向量的選擇可以是任意的。
2. 影響冪法的收斂速度主要是\(\lambda_2/\lambda_1\)的值，如果這個值越小，收斂速度越快。

6.3 反冪法

設矩陣\(A\in\mathbb{R}^{n\times n}\)是非奇異的，\(\lambda\)表示\(A\)的特征值，\(x\)則是相應的特征向量，如果有\(\lambda\ne0\)，那么

\[Ax=\lambda x\Leftrightarrow A^{-1}x=\lambda^{-1}x \]

這說明\(\lambda^{-1}\)是\(A^{-1}\)的特征值，\(x\)則是相應的特征向量。

基於上述這個性質，反冪法常用來求矩陣按模最小的特征值及其對應的特征向量。

假設\(A\)的特征值滿足

\[|\lambda_1|\ge|\lambda_2|\ge\cdots\ge|\lambda_{n-1}|>|\lambda_n|>0 \]

那么\(A^{-1}\)的特征值（\(\lambda_1^{-1},\cdots,\lambda_n^{-1}\)）滿足

\[\frac{1}{|\lambda_n|}>\frac{1}{|\lambda_{n-1}|}\ge\cdots\ge \frac{1}{|\lambda_1|} \]

通過冪法求出\(A^{-1}\)按模最大的特征值就是\(A\)按模最小的特征值。算法如下：

\[\left\{\begin{aligned} &j=\mathop{\arg\max}_i|V_i^{(k)}|\\ &m_k=V_j^{(k)}\\ &U^{(k)}=V^{(k)}/m_k\\ &V^{(k+1)}=A^{-1}U^{(k)} \end{aligned}\right.\quad k=0,1,2,\cdots \]

當\(|1/m_k-1/m_{k-1}|<\varepsilon\)或達到最大迭代上限時，算法停止並返回\(1/m_k\)作為按模最小特征值，\(U^{(k)}\)作為特征向量。

直接計算\(A\)的逆矩陣非常耗時，通常通過從方程\(AV^{(k+1)}=U^{(k)}\)中求解出\(V^{(k+1)}\)，來代替算法的最后一步。此外，為了進一步節省工作量可以事先將\(A\)進行LU分解。那么之后迭代過程中，就只需要求解兩個三角形方程組。當矩陣階數很大時，將有效減少計算量。

6.4 原點平移法

設矩陣\(A\in\mathbb{R}^{n\times n}\)的特征值為\(\lambda_1,\cdots,\lambda_n\)，\(x_1,\cdots,x_n\)則是相應的特征向量，那么

\[Ax_i=\lambda_ix_i\Leftrightarrow (A-qI)x_i=(\lambda_i-q)x_i \]

其中\(q\)是一個實數。

通過上述式子可知，矩陣\(A-qI\)的特征值為\(\lambda_1-q,\cdots,\lambda_n-q\)，這相當於是\(A\)特征值進行了步長為\(q\)的平移，並且保持特征向量不變。

利用這一點性質，如果我們想求出矩陣\(A\)距離\(q\)最近的特征值，那么就可以通過原點平移法將這個最接近\(q\)的特征值平移到0附近。然后使用反冪法求解即可，具體做法就是：

1. 令\(B=A-qI\)。
2. 使用反冪法求解\(B\)的按模最小特征值\(1/m_k\)，特征向量\(U^{(k)}\)。
3. 返回\(1/m_k+q\)作為特征值，\(U^{(k)}\)作為特征向量。

6.5 矩陣的正交分解

6.5.1 Householder變換

定義 設\(u\in\mathbb{R}^n\)，且\(||u||_2=1\)，則矩陣

\[H=I-2uu^T \]

稱為Householder矩陣或者Householder變換，也可以簡稱為H矩陣或H變換。

定理1 H矩陣具有如下性質：

1. \(H\)是對稱矩陣：\(H=H^T\)。
2. \(H\)是正交矩陣：\(H^TH=I,H^T=H^{-1}\)。
3. \(H\)變換保持向量長度不變：\(\forall v\in\mathbb{R}^n,||Hv||_2=||v||_2\)。
4. 設\(S\)是經過原點的超平面，且法向量為\(u\)，那么對於任意非零向量\(v\in\mathbb{R}^n\)，都有\(Hv\)與\(v\)關於超平面\(S\)對稱。

證明：

性質1，2可以很容易被驗證。

因為

\[||Hv||_2^2=(Hv)^T(Hv)=v^T(H^TH)v=v^Tv=||v||_2^2 \]

即\(||Hv||_2=||v||_2\)，性質3成立。

為了證明性質4，假設\(v=v_{\parallel}+v_{\perp}\)，其中\(v_{\parallel}\)是與超平面平行的分量，那么\(u^Tv_{\parallel}=0\)；\(v_{\perp}\)是與超平面垂直的分量，因此必然存在實數\(\lambda\)，使得\(v_{\perp}=\lambda u\)。

\[\begin{aligned} Hv&=Hv_{\parallel}+Hv_{\perp}\\ &=(v_{\parallel}-2uu^Tv_{\parallel})+(v_{\perp}-2uu^Tv_{\perp})\\ &=v_{\parallel}+(\lambda{u}-2\lambda{uu^Tu})\\ &=v_{\parallel}-\lambda{u}\\ &=v_{\parallel}-v_{\perp} \end{aligned} \]

\(Hv=v_{\parallel}-v_{\perp}\)表明\(Hv\)和\(v\)是關於超平面\(S\)對稱的。

定理2 任意非零向量\(v\in\mathbb{R}^n\)，都可以構造出一個合適的H矩陣使得\(Hv=ce_1\)。

其中\(c=-sign(v_1)||v||_2,e_1=[1,0,\cdots,0]^T\)。\(v_1\)表示向量\(v\)的第一個分量。

構造方法：

1. 計算\(w=v-ce_1=v+sign(v_1)||v||_2e_1\)，此時\(w\)就是超平面\(S\)的一個法向量。
2. 令\(u=w/||w||_2\)，便可以計算出\(H\)。

\[H=I-2\frac{ww^T}{||w||_2^2} \]

令\(\beta=2/||w||_2^2\)，則

\[\begin{aligned} \beta=\frac{2}{||w||_2^2}&=\frac{2}{(v-ce_1)^T(v-ce_1)}\\ &=\frac{2}{v^Tv-cv^Te_1-ce_1^Tv+c^2e_1^Te_1}\\ &=\frac{1}{||v||_2(||v||_2+|v_1|)} \end{aligned} \]

於是\(H=I-\beta{ww^T}\)。

定理2的推廣 可以找到一個H矩陣，將一個非零向量的某幾個連續分量變為零。

假設\(v\in\mathbb{R}^n,c=-sign(v_i)\sqrt{v_i^2+v_{i+1}^2+\cdots+v_j^2},2\le i,j\le n\)。我們的目標就是通過H變換將\(v\)變換為\(Hv=[v_1,\cdots,c,0,\cdots,0,v_{j+1},\cdots,v_n]\)，即將第\(i+1\)到\(j\)個分量變成0。

構造方法：

1. 計算\(w=[0,\cdots,0,v_i-c,v_{i+1},\cdots,v_j,0,\cdots,0]\)。
2. 令 \(u=w/||w||_2\)即可算出\(H\)。

令\(\bar{v}=[v_i,\cdots,v_j]^T\)，與定理2中的分析相同，我們可以得到

\[\beta=\frac{1}{||\bar{v}||_2(||\bar{v}||_2+|v_i|)} \]

於是\(H=I-\beta{ww^T}\)。

6.5.2 Hessenberg矩陣

定義 設矩陣\(A=[a_{ij}]_{n\times{n}}\)，當\(i>j+1\)時\(a_{ij}=0\)，則稱矩陣\(A\)為上hessenberg矩陣（或擬上三角形矩陣）；當\(j>i+1\)時\(a_{ij}=0\)，則稱矩陣\(A\)為下hessenberg矩陣（或稱擬下三角形矩陣）。

以4階矩陣為例：

\[\left[\begin{matrix} a_{11}&a_{12}&a_{13}&a_{14}\\ a_{21}&a_{22}&a_{23}&a_{24}\\ 0&a_{32}&a_{33}&a_{34}\\ 0&0&a_{43}&a_{44}\\ \end{matrix}\right]\quad \left[\begin{matrix} a_{11}&a_{12}&0&0\\ a_{21}&a_{22}&a_{23}&0\\ a_{31}&a_{32}&a_{33}&a_{34}\\ a_{41}&a_{42}&a_{43}&a_{44}\\ \end{matrix}\right] \]

左側是上hessenberg矩陣，右側是下hessenberg矩陣。

定理3 對於任意n階矩陣\(A\)，總存在正交矩陣\(Q\)，使得\(B=QAQ^{-1}\)為上hessenberg矩陣。即任意n階矩陣總相似於一個上hessenberg矩陣。

證明：證明的思路就是構造householder矩陣，將任意的n階矩陣逐步變成上hessenberg矩陣，如果

\[A=\left[\begin{matrix} a_{11}&\cdots&\cdots&\cdots&a_{1n}\\ a_{21}&\cdots&\cdots&\cdots&a_{2n}\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ a_{n1}&\cdots&\cdots&\cdots&a_{nn} \end{matrix}\right] \]

我們的目標就是找到一個\(\bar{H_1}\in\mathbb{R}^{(n-1)\times(n-1)}\)使得

\[\bar{H_1}\left[\begin{matrix} a_{21}\\a_{31}\\\vdots\\a_{n1} \end{matrix}\right]= \left[\begin{matrix} c_1\\0\\\vdots\\0 \end{matrix}\right] \]

顯然根據定理2我們可以找到這樣的\(\bar{H_1}\)，然后構造一個n階矩陣

\[H_1=\left[\begin{matrix} 1&0\\ 0&\bar{H_1} \end{matrix}\right] \]

可以驗證\(H_1\)是正交且對稱的矩陣（\(H_1=H_1^T=H_1^{-1}\)），於是利用矩陣分塊乘積可以得到

\[\begin{aligned} H_1AH_1&=\left[\begin{matrix} 1&0\\ 0&\bar{H_1} \end{matrix}\right] \left[\begin{matrix} a_{11}&\cdots&\cdots&\cdots&a_{1n}\\ a_{21}&\cdots&\cdots&\cdots&a_{2n}\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ a_{n1}&\cdots&\cdots&\cdots&a_{nn} \end{matrix}\right] \left[\begin{matrix} 1&0\\ 0&\bar{H_1} \end{matrix}\right]\\ &=\left[\begin{matrix} *&\cdots&\cdots&\cdots&*\\ c_1&\cdots&\cdots&\cdots&*\\ 0&\cdots&\cdots&\cdots&\vdots\\ \vdots&\cdots&\cdots&\cdots&\vdots\\ 0&\cdots&\cdots&\cdots&* \end{matrix}\right]=A_1 \end{aligned} \]

其中\(*\)表示一個數。

因此\(A_1=H_1AH_1\)，即矩陣\(A_1\)相似於\(A\)。

基於這個思想，第\(k\)步變換的矩陣為

\[H_k=\left[\begin{matrix} I_{k\times k}&0\\ 0&\bar{H_k} \end{matrix}\right] \]

一共需要\(n-2\)步就能夠得到上hessenberg矩陣，即

\[H_{n-2}\cdots H_2H_1AH_1H_2\cdots H_{n-2}=\left[\begin{matrix} *&\cdots&\cdots&\cdots&*\\ c_1&\cdots&\cdots&\cdots&*\\ 0&c_2&\cdots&\cdots&\vdots\\ \vdots&\ddots&\ddots&\cdots&\vdots\\ 0&\cdots&0&c_{n-1}&* \end{matrix}\right]=A_{n-2} \]

於是\(A_{n-2}=HAH^{-1}\)，其中\(H=H_{n-2}\cdots H_2H_1\)。這說明經過\(n-2\)步變換之后，一定可以得到上hessenberg矩陣，並且\(A_{n-2}\)相似於\(A\)。

因此所求的正交矩陣\(Q=H\)。

6.5.3 QR分解

定理4 對於任意n階矩陣\(A\)，總存在正交矩陣\(Q\)和上三角陣\(R\)，使得\(A=QR\)。

證明：證明方法與定理3中的類似，根據householder變換，構造出一系列的\(H_1,H_2,\cdots,H_{n-1}\)，使得

\[H_{n-1}\cdots H_2H_1A=R\Leftrightarrow A=(H_1H_2\cdots H_{n-1})R \]

其中\(Q=H_1H_2\cdots H_{n-1}\)，\(R\)就是上三角陣。

需要注意的是，任意的n階矩陣\(A\)總可以進行QR分解，但是其分解一般是不唯一的。但是如果\(A\)是可逆矩陣，且上三角陣\(R\)的對角線元素均取正值，則\(A\)的QR分解是唯一的。

QR方法 是冪法的一種推廣和變形，可以用來求解任意矩陣的特征值，能夠有效地解決階數不高的任意實方陣的全部特征值問題。

因為實矩陣\(A\)總可以進行QR分解，於是\(A=QR\)。並且令\(B=RQ\)，那么

\[B=RQ=(Q^{-1}A)Q=Q^{-1}AQ \]

即\(B\)與\(A\)相似，因此它們有相同的特征值。基於這個想法，可以得到基本的QR算法：

1. 令\(A_1=A\)。
2. 對於\(k=1,2,\cdots\)，作\(A_k\)的QR分解：\(A_k=Q_kR_k\)。
3. \(A_{k+1}=R_kQ_k\)。

可以驗證：\(A_{k+1}=Q_k^{-1}A_kQ_k\)，於是\(A_{k+1}\)總是和\(A_k\)相似。因此，在整個迭代過程中總是保持矩陣\(A\)的特征值不變。

實Schur分解定理 對於任意方陣\(A\in\mathbb{R}^{n\times n}\)，QR方法產生的\(\{A_k\}\)序列收斂於形如下式的分塊上三角陣：

\[\left[\begin{matrix} R_{11}&R_{12}&\cdots&R_{1n}\\ 0&R_{22}&\cdots&R_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&R_{nn}\\ \end{matrix}\right] \]

其中\(R_{ii}\)為實數或者\(2\times2\)的矩陣塊。如果\(R_{ii}\)是一個實數，那么\(R_{ii}\)就是矩陣\(A\)的一個特征值；如果\(R_{ii}\)是一個\(2\times2\)的矩陣塊，那么\(R_{ii}\)的特征值就是\(A\)的一對共軛復數特征值。

QR方法需要注意的問題

1. 在實際中直接使用QR方法時，為了加快收斂，可以先將\(A\)轉化為上hessenberg矩陣，然后在使用QR分解求特征值。
2. 迭代停止的條件：（1）設置最大迭代上限。（2）設定一個誤差上限\(\varepsilon\)，只要某一個范數滿足\(||diag(A_{k+1}-A_k)||<\varepsilon\)成立，停止迭代。