Cholesky 分解

本文的閱讀等級:中級

任何一個實對稱正定矩陣 A 都可唯一分解為 A=GG^T;若 A 為 Hermitian 正定矩陣,則有 A=GG^{\ast},其中 G 是一個下三角矩陣且主對角元皆為正數。這個分解式稱為 Cholesky 分解,由法國砲兵軍官喬列斯基 (André-Louis Cholesky) 所提出,當初的發明動機是為了解決測地計算問題。喬列斯基本人從未公開發表他的演算法,他於1918年的一場戰事中身亡,之後經口耳相傳,這個卓越的矩陣分解法才逐漸為時人所知。本文從 LU 分解推導 Cholesky 分解,並介紹一個演算法。

 
Cholesky 分解

A 為一個 n\times n 階實對稱矩陣,如欲延伸至 Hermitian 複矩陣,僅需將以下轉置運算 (\cdot)^T 替換為共軛轉置 (\cdot)^\ast (見“特殊矩陣 (9):Hermitian 矩陣”)。在不使用任何列交換的情況下,若 A 可分解為 A=LDU,其中 L 為下三角矩陣,U 為上三角矩陣,LU 的主對角元皆為 1D=\mathrm{diag}(d_{11},\ldots,d_{nn}),則 A 的 LDU 分解式可以進一步化約為 A=LDL^T。理由如下:將 A=LDU 代入 A^T=A,推得 U^TDL^T=LDU。如果 D 的主對角不含零元,LDU 分解是唯一的 (見“LU 分解”),因此推論 U=L^T。更進一步,如果 A 是正定矩陣 (見“正定矩陣的性質與判別方法”),則對於所有的非零實向量 \mathbf{x}

\displaystyle \mathbf{x}^TA\mathbf{x}=\mathbf{x}^TLDL^T\mathbf{x}=\mathbf{y}^TD\mathbf{y}=\sum_{i=1}^nd_{ii}y_i^2>0

上面令 \mathbf{y}=(y_1,\ldots,y_n)^T=L^T\mathbf{x}。因為 L 是可逆矩陣,任何一個非零向量 \mathbf{y} 都滿足上述不等式。設 \mathbf{y} 為標準單位向量 \mathbf{e}_i (第 i 元是 1,其餘元是 0),推知每一 d_{ii} 為正數。令 D^{1/2}=\mathrm{diag}(\sqrt{d_{11}},\ldots,\sqrt{d_{nn}})。因此,

A=LDL^T=LD^{1/2}D^{1/2}L^T=(LD^{1/2})(LD^{1/2})^T

實對稱正定矩陣 A 的 Cholesky 分解即為 A=GG^T,其中 G=LD^{1/2} 是包含正主對角元的下三角矩陣。由於可逆矩陣的 LDU 分解是唯一的,正定矩陣的 Cholesky 分解也是唯一的。下為一例,

\left[\!\!\begin{array}{rr}    3&-3\\    -3&5    \end{array}\!\!\right]=\left[\!\!\begin{array}{rc}    1&0\\    -1&1    \end{array}\!\!\right]\begin{bmatrix}    3&0\\    0&2    \end{bmatrix}\left[\!\!\begin{array}{cr}    1&-1\\    0&1    \end{array}\!\!\right]=\left[\!\!\begin{array}{rr}    \sqrt{3}&0\\    -\sqrt{3}&\sqrt{2}    \end{array}\!\!\right]\left[\!\!\begin{array}{rr}    \sqrt{3}&-\sqrt{3}\\    0&\sqrt{2}    \end{array}\!\!\right]

 
正規方程

Cholesky 分解常見於求解最小平方法的正規方程 (normal equation) A^TA\hat{\mathbf{x}}=A^T\mathbf{b},此處 A 為一 m\times n 階矩陣,且 \mathrm{rank}A=n。Cholesky 分解得到 A^TA=GG^T,正規方程可表示為 GG^{T}\hat{\mathbf{x}}=A^T\mathbf{b}=\mathbf{d}。先以正向代入法解 G\mathbf{y}=\mathbf{d},跟著再用反向代入法解 G^T\hat{\mathbf{x}}=\mathbf{y},應用原理與 LU 分解相同。

 
QR 分解

Cholesky 分解和 QR 分解也有些關係。設 A 的 QR 分解為 A=QR (見“Gram-Schmidt 正交化與 QR 分解”),其中 Q 是正交矩陣 (orthogonal matrix),Q^T=Q^{-1}R 是上三角矩陣。將 A=QR 代入計算交互乘積 A^TA,可得

A^TA=(QR)^T(QR)=R^TQ^TQR=R^TR

結果顯示 A^TA 的 Cholesky 下三角矩陣 G 其實就是 A 的 QR 分解的 R^T 矩陣。

 
計算方法

Cholesky 分解有許多種計算方法,n 階方陣的計算量大都接近n^3/3,下面我們用三階實對稱正定矩陣說明一種名為 Cholesky-Crout 演算法。考慮

\begin{aligned}  A&=\begin{bmatrix}    a_{11}&a_{12}&a_{13}\\    a_{21}&a_{22}&a_{23}\\    a_{31}&a_{32}&a_{33}    \end{bmatrix}=GG^T\\  &=\begin{bmatrix}    g_{11}&0&0\\    g_{21}&g_{22}&0\\    g_{31}&g_{32}&g_{33}    \end{bmatrix}\begin{bmatrix}    g_{11}&g_{21}&g_{31}\\    0&g_{22}&g_{32}\\    0&0&g_{33}    \end{bmatrix}\\    &=\begin{bmatrix}    g_{11}^2&g_{21}g_{11}&g_{31}g_{11}\\    g_{21}g_{11}&g_{21}^2+g_{22}^2&g_{31}g_{21}+g_{32}g_{22}\\    g_{31}g_{11}&g_{31}g_{21}+g_{32}g_{22}&g_{31}^2+g_{32}^2+g_{33}^2    \end{bmatrix}.\end{aligned}

G 矩陣的左上角開始,由左至右逐次解出 G 的各行向量,如下:

\begin{aligned}\displaystyle  g_{11}&=\sqrt{a_{11}}\\    g_{21}&=\frac{a_{21}}{g_{11}}\\    g_{31}&=\frac{a_{31}}{g_{11}}\\    g_{22}&=\sqrt{a_{22}-g_{21}^2}\\    g_{32}&=\frac{1}{g_{22}}\left(a_{32}-g_{31}g_{21}\right)\\    g_{33}&=\sqrt{a_{33}-g_{31}^2-g_{32}^2}.\end{aligned}

由此可以推出一般公式:對於 j=1,\ldots,n

\displaystyle  g_{jj}=\sqrt{a_{jj}-\sum_{k=1}^{j-1}g_{jk}^2}

\displaystyle  g_{ij}=\frac{1}{g_{jj}}\left(a_{ij}-\sum_{k=1}^{j-1}g_{ik}g_{jk}\right),~~i=j+1,j+2,\ldots,n.

注意,正定矩陣 A 的每一個主對角元 a_{jj} 都是正數。在一般情況下,Cholesky 分解具有優異的數值穩定性及精確性,但如果 A 是病態矩陣,以平方根算出的 g_{jj} 會出現極小的數值,捨入誤差 (roundoff error) 甚至可能造成 g_{jj}<0。萬一出現這個狀況,計算程序就必須停止。當 Cholesky 分解不可行時,我們通常轉而使用其他需要較多計算量的分解式,如 QR 分解或 SVD 分解。

This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , , . Bookmark the permalink.

13 則回應給 Cholesky 分解

  1. Vahi Chen 說:

    在Cholesky分解推导过程中,为什么由\sum_{i=1}^n d_{ii}y_i^2>0就能得出每一个d_{ii}都是正数呢?

  2. Rongwei Sun 說:

    老師,請問您,複數矩陣可以Cholesky 分解嗎?

  3. Murphy Ma 說:

    老師,請問您:如果一個full row rank實數矩陣A與其Transpose相乘生成的正定矩陣B=A*A‘,如何利用cholesky分解法求解(A*A‘+ alpha* I)的inverse矩陣(I是單位矩陣)?若A是full column rank實數矩陣,重新計算後兩個的答案是否相同?渴望您的指導,Thanks a lot!

  4. ccjou 說:

    實對稱正定矩陣 C=AA^T+\alpha I 的 Cholesky 分解 C=GG^T 給出 C^{-1}=L^TL,其中 L=G^{-1}。我不明白何謂「重新計算後兩個的答案是否相同」?

    • Murphy Ma 說:

      非常感謝老師的回復!抱歉我的問題沒有說清楚:問題是求解方程Cx=q的解x=C\q,其中q=A’ * b+\alpha * z,C=A’ * A+ \alpha* I(A是full column rank),因C是一個high dimension的超大矩陣,現在想用cholesky分解x=C\q=L’ \ L\ q(A是full column rank,能理解),但是對於A是full row rank時,C=1/ alpha*A*A’+I,看到的solution是x=q / alpha -(A’*(L’ \ L\(A*q)))/ alpha^2,不明白為什麼?請老師指點迷津,謝謝!

  5. 徐君 說:

    是否任意一个正定矩阵可以可以有无穷多个这样的分解:A=GG^T,其中G是维度合适的矩阵。

發表迴響

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

WordPress.com Logo

你正使用 WordPress.com 帳號留言。 登出 / 變更 )

Twitter picture

你正使用 Twitter 帳號留言。 登出 / 變更 )

Facebook照片

你正使用 Facebook 帳號留言。 登出 / 變更 )

Google+ photo

你正使用 Google+ 帳號留言。 登出 / 變更 )

連結到 %s