Jacobi 特徵值算法

本文的閱讀等級:中級

給定一 n\times n 階實矩陣 A,使用 Householder 變換可求得一正交矩陣 (orthogonal matrix) QQ^T=Q^{-1},使得 B=[b_{ij}]=Q^TAQ 為三對角 (tridiagonal) 矩陣,其中 b_{ij}=0\vert i-j\vert>1 (見“特殊矩陣 (19):Hessenberg 矩陣”)。此外,如果 A 是實對稱矩陣,利用旋轉矩陣可求出一正交矩陣 Q 使得 D=Q^TAQ 為對角矩陣。因為 D 正交相似於 AD 的主對角元即為 A 的特徵值。德國數學家雅可比 (Carl Gustav Jacob Jacobi) 於1846年公開這個對角化實對稱矩陣的計算方法,後人稱之為 Jacobi 特徵值算法。

 
A 為一 n\times n 階實對稱矩陣。為使特徵值不變並維持數值計算的穩定性,雅可比考慮正交相似變換

A\to Q^TAQ

其中 Q 是一正交矩陣。經由特殊設計的正交矩陣 Q,我們可以消滅 A 的任何一個非主對角元。考慮 n=2 的情況。二階正交矩陣可區分為兩種,旋轉矩陣 R 與鏡射矩陣 H (見“旋轉與鏡射”),如下:

R=\left[\!\!\begin{array}{cr}  \cos\theta&-\sin\theta\\  \sin\theta&\cos\theta\end{array}\!\!\right],~~  H=\left[\!\!\begin{array}{cr}  \cos\theta&\sin\theta\\  \sin\theta&-\cos\theta\end{array}\!\!\right]

正交矩陣 Q 滿足 \vert\det Q\vert=1。因為 \det R=1\det H=-1,旋轉矩陣稱為適當的 (proper) 正交矩陣,鏡射矩陣則稱為不適當的正交矩陣。考慮 B=R^TAR,為便利推廣至高階矩陣,我們用指標 ij 取代數字 12,表示如下:

\displaystyle  B=\begin{bmatrix}  b_{ii}&b_{ij}\\  b_{ji}&b_{jj}  \end{bmatrix}=\left[\!\!\begin{array}{rc}  \cos\theta&\sin\theta\\  -\sin\theta&\cos\theta\end{array}\!\!\right]  \begin{bmatrix}  a_{ii}&a_{ij}\\  a_{ij}&a_{jj}  \end{bmatrix}  \left[\!\!\begin{array}{cr}  \cos\theta&-\sin\theta\\  \sin\theta&\cos\theta\end{array}\!\!\right]=R^TAR

乘開可得 b_{ij}=\cos\theta\sin\theta(a_{jj}-a_{ii})+(\cos\theta^2-\sin\theta^2)a_{ij}。若 2\theta=\arctan\frac{2a_{ij}}{a_{ii}-a_{jj}},由倍角公式 \cos 2\theta=\cos^2\theta-\sin^2\theta\sin 2\theta=2\cos\theta\sin\theta 可推得 b_{ij}=b_{ji}=0,證明二階實對稱矩陣可正交對角化。

 
這裡我們要特別強調正交相似變換 B=R^TAR 不改變矩陣的 Frobenius 範數。對於 n\times n 階實矩陣 M=[m_{pq}],Frobenius 範數定義為 (見“矩陣範數”)

\displaystyle\Vert M\Vert_F=\sqrt{\hbox{trace}(M^TM)}=\sqrt{\sum_{p=1}^n\sum_{q=1}^n m^2_{pq}}

使用定義與跡數循環不變性,

\displaystyle\begin{aligned}  \Vert B\Vert_F^2&=\hbox{trace}(B^TB)=\hbox{trace}(R^TA^TRR^TAR)=\hbox{trace}(R^TA^TAR)\\  &=\hbox{trace}(A^TARR^T)=\hbox{trace}(A^TA)=\Vert A\Vert^2_F.\end{aligned}

因此,前述二階正交相似變換滿足 b_{ii}^2+b_{jj}^2=a_{ii}^2+a_{jj}^2+2a^2_{ij},形象化的說法是旋轉矩陣實現的正交相似變換把非主對角元成分「分配」到主對角元。

 
對於 n>2,假設我們想要消滅 A 的最大 (絕對值) 非主對角元 a_{ij}\neq 0,稱為軸元 (pivot)。考慮在 \mathbb{R}^n 中第 i 軸與第 j 軸所張開的平面上旋轉,稱為 Givens 旋轉[1](見“Givens 旋轉於 QR 分解的應用”),表達式如下:

Q=\begin{bmatrix}  1&\cdots&0&\cdots&0&\cdots&0\\    \vdots&\ddots&\vdots&~&\vdots&~&\vdots\\    0&\cdots&c&\cdots&-s&\cdots&0\\    \vdots&~&\vdots&\ddots&\vdots&~&\vdots\\    0&\cdots&s&\cdots&c&\cdots&0\\    \vdots&~&\vdots&~&\vdots&\ddots&\vdots\\    0&\cdots&0&\cdots&0&\cdots&1    \end{bmatrix}

其中 c=\cos\thetas=\sin\theta。Givens 旋轉矩陣 Q=[q_{ij}] 完全複製單位矩陣 I,除了四個元例外:q_{ii}=q_{jj}=cq_{ij}=-sq_{ji}=s。令 B=Q^TAQ。同前,選擇適宜的 \theta 可消滅 a_{ij},也就是說 b_{ij}=b_{ji}=0。對於 k\neq i,j,不難確認 b_{kk}=a_{kk}。為度量 AB 的非主對角元成分的消長,令 n\times n 階矩陣 A^\circ 的主對角元為零,非主對角元等於 A 的非主對角元。同樣地,B^\circ 也按照這個方式定義。使用 \Vert B\Vert_F=\Vert A\Vert_Fb^2_{ii}+b^2_{jj}=a^2_{ii}+a^2_{jj}+2a^2_{ij},可得

\displaystyle  \begin{aligned}  \Vert B^\circ\Vert_F^2&=\Vert B\Vert^2_F-\sum_{k=1}^nb_{kk}^2=\Vert B\Vert^2_F-(b^2_{ii}+b^2_{jj})-\sum_{k\neq i,j}b^2_{kk}\\  &=\Vert A\Vert_F^2-(a^2_{ii}+a^2_{jj}+2a^2_{ij})-\sum_{k\neq i,j}a_{kk}^2\\  &=\Vert A^\circ\Vert_F^2-2a_{ij}^2.  \end{aligned}

a_{ij}\neq 0\Vert B^\circ\Vert_F<\Vert A^\circ\Vert_F 表示經過 Givens 旋轉變換,A 的整體非主對角元成分隨之減少。另一方面,對於所有的 1\le p,q\le na^2_{pq}\le a^2_{ij},就有

\displaystyle  \Vert A^\circ\Vert^2_F=\sum_{p\neq q}a^2_{pq}\le\sum_{p\neq q}a_{ij}^2=(n^2-n)a^2_{ij}

其中 n^2-n 是所有的非主對角元數。合併上面兩式,

\displaystyle  \Vert B^\circ\Vert^2_F=\Vert A^\circ\Vert_F-2a^2_{ij}\le\Vert A^\circ\Vert_F^2-2\frac{\Vert A^\circ\Vert_F^2}{n^2-n}=\left(1-\frac{2}{n^2-n}\right)\Vert A^\circ\Vert_F^2

 
Jacobi 特徵值算法不過就是重複執行上述消 (軸) 元運算。給定一 n\times n 階實對稱矩陣 A,設 A_0=A,計算矩陣序列 A_k=Q_{k}^TA_{k-1}Q_kk=1,2,\ldots。正交矩陣 Q_{k} 的設計方法如下:找出 A_{k-1} 的 (非零) 軸元,設為 a_{ij}。設 \tau=\frac{a_{ii}-a_{jj}}{2a_{ij}}。令 Q_{k}(i,j) 平面上的 Givens 旋轉矩陣,參數值如下:

\displaystyle  c=\frac{1}{\sqrt{1+\rho^2}},~~s=\frac{\rho}{\sqrt{1+\rho^2}}=\rho c

其中 \rho=\hbox{sign}(\tau)(-\vert\tau\vert+\sqrt{1+\tau^2}) (推導見[2])。沿用前述記號。對於 n>2,當 k\to\infty

\displaystyle  \Vert A^\circ_k\Vert_F^2\le\left(1-\frac{2}{n^2-n}\right)^k\Vert A^\circ\Vert_F^2\to 0

令正交矩陣 U_k=Q_1Q_2\cdots Q_kk\ge 1。上式蘊含

\displaystyle  \lim_{k\to\infty}U_k^TAU_k=\lim_{k\to\infty}A_k=D

其中 D 是一對角矩陣。古典 Jacobi 特徵值算法為尋得最大非主對角元必須耗費 \frac{n^2-n}{2} 次比較。為免除這個運算成本,一種 Jacobi 算法的變形,稱為逐列循環 (cyclic-by-row),按照列序逐一遍歷 (上三角) 非主對角元,計算 \frac{n^2-n}{2} 次 (旋轉) 相似變換,稱為一次席捲 (sweep)。若以席捲為計算單位,可以證明逐列循環同樣具有古典 Jacobi 算法的平方收斂性能。

 
最後我介紹美國賓州大學教授維爾夫 (Herbert Wilf) 提出的 Jacobi 特徵值算法所啟發的實對稱矩陣可正交對角化證明[3] (另一個證明見[4])。令 O(n)O(n,\mathbb{R}) 為所有的 n\times n 階實正交矩陣 Q=[q_{ij}] 形成的集合 (稱為正交群)。因為每一 \vert q_{ij}\vert\le 1O(n) 是有界的 (bounded)。另外,Q 滿足 QQ^T=I=[\delta_{ij}],即

\displaystyle  \sum_{k=1}^nq_{ik}q_{jk}=\delta_{ij},~~1\le i,j\le n

可知 O(n) 是閉合的 (closed)。對於 \mathbb{R}^{n\times n} (即 n\times n 階實矩陣構成的集合) 的子集 S,Heine-Borel 定理宣稱下列兩個陳述是等價的:(1) S 是有界的,而且是閉合的,(2) S 是緊緻的 (compact)。所以,O(n) 是緊緻的。考慮 n>2,推論步驟如下:

  1. 考慮正交相似變換 T_A(Q)=Q^TAQ,其中 Q\in O(n)。因為 T_A 是在緊集 (compact set) O(n) 的一個連續映射,值域 T_A(O(n)) 也是緊緻的。
  2. 定義於緊集 T_A(O(n)) 的函數 f(B)=\Vert B^\circ\Vert_F\in[0,\infty) 是連續的,因此有最小值。設最小值位於 D=U^TAU\in T_A(O(n))
  3. f(D)=\Vert D^\circ\Vert_F 必定為零,也就是說 D 是一個對角矩陣。使用反證法。假設最小值 f(D)>0,即 D 至少有一個非主對角元不為零。設計 Givens 旋轉矩陣 Q 消滅 D 的軸元可使 f(Q^TDQ)<f(D),其中 Q^TDQ=(UQ)^TA(UQ)\in T_A(O(n)),得到一個矛盾。

 
這個匠心獨具的證法有別於傳統的線性代數證明 (見“實對稱矩陣可正交對角化的證明”),其中最為獨特之處是既沒有闡述實對稱矩陣的特徵值為實數,也沒有演繹特徵向量構成一組單範正交集 (orthonormal set)。論證直面正交對角化,思路簡單扼要:通過正交相似變換消滅實對稱矩陣的 (非主對角) 軸元,使得整體非主對角元成分減少,再套用緊緻性斷言所有的非主對角元定可被悉數殲滅。

 
註解
[1] 有些文獻稱之為 Jacobi 旋轉以便與應用於 QR 分解的 Givens 旋轉有所區別,後者由美國數學家吉文斯 (Wallace Givens) 於1950年代引入數值分析。
[2] 使用倍角公式 \tan 2\theta=\frac{2\tan\theta}{1-\tan^2\theta}。設 \tau=\frac{1}{\tan 2\theta},則 \tan^2\theta+2\tau\tan\theta-1=0,解得 \rho= \tan\theta=-\tau\pm\sqrt{1+\tau^2}。選擇絕對值較小的解,即 \rho=\hbox{sign}(\tau)(-\vert\tau\vert+\sqrt{1+\tau^2})
[3] Herbert S. Wilf, An algorithm-inspired proof of the spectral theorem in E^n, Amer. Math. Monthly, Vol. 88, 1981, pp 49-50.
[4] 因為我們已經得到 \lim_{k\to\infty}U_k^TAU_k=D,其中 U_k\in O(n)D 是對角矩陣,在此條件下,另一個快捷證明是直接引用 Bolzano-Weierstrass 定理:\mathbb{R}^{n\times n} 中每一個有界序列有一個收斂子列 (subsequence)。因此,從有界序列 \{U_k\} 取一個收斂子列,令 U 為極限,即證明 U^TAU=D

Advertisements
本篇發表於 線性代數專欄, 數值線性代數 並標籤為 , , , , 。將永久鏈結加入書籤。

2 則回應給 Jacobi 特徵值算法

  1. Meiyue Shao 說道:

    我加些注释吧
    1. Jacobi 算法被认为是数值线性代数这个分支的开端, 用于将 2 阶实对称阵对角化的旋转矩阵也称为 Jacobi 旋转, Givens 旋转这个名称则常用于 QR 分解的场合, 以示区别.
    2. 按 \rho=\mathrm{sign}(\tau)/(|\tau|+\sqrt{1+|\tau|^2}) 来写比较好, 一来这就是用于实际计算的公式, 二来避免潜在的由“xxx或xxx”产生的误解.
    3. 实际计算时 Jacobi 算法的主元的选取并不是选模最大的非对角元, 而是对非对角元逐一遍历, 当然相应的理论分析要麻烦一些.
    4. 在介绍完 Jacobi 算法的收敛性时其实已经得到了一个对角化的证明, 由 Bolzano–Weierstrass 定理对 U_k 取一个收敛子列即可.

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s