Cholesky 分解

Posted on 09/16/2010 by ccjou

本文的閱讀等級：中級

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

Andre-Louis Cholesky (1875-1918) From http://www-history.mcs.st-andrews.ac.uk/BigPictures/Cholesky.jpeg

 
Cholesky 分解

假設 A 是一個 n\times n 階實對稱矩陣，如欲延伸至 Hermitian 複矩陣，僅需將以下轉置運算 (\cdot)^T 替換為共軛轉置 (\cdot)^\ast (見“特殊矩陣 (9)：Hermitian 矩陣”)。在不使用任何列交換的情況下，若 A 可分解為 A=LDU，其中 L 為下三角矩陣，U 為上三角矩陣，L 和 U 的主對角元皆為 1，D=\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 矩陣。

 
蒙地卡羅法

蒙地卡羅法 (Monte Carlo method) 是一種統計模擬方法，透過電腦模擬產生的亂數計算複雜的和或積分，以及某個隨機事件的機率。假設我們要產生多個隨機變數 \mathbf{x}=(x_1,\ldots,x_n)^T 且共變異數矩陣為 \hbox{cov}[\mathbf{x}]=\Sigma。首先產生一個隨機向量 \mathbf{z} 滿足 \hbox{cov}[\mathbf{z}]=I，接著計算 Cholesky 分解 \Sigma=GG^T，並令 \mathbf{x}=G\mathbf{z}，就有 (見“共變異數矩陣的性質”)

\hbox{cov}[\mathbf{x}]=\hbox{cov}[G\mathbf{z}]=G\,\hbox{cov}[\mathbf{z}]G^T=GG^T=\Sigma。

 
計算方法

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 分解。

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

22 Responses to Cholesky 分解

  1. Vahi Chen says:
    10/10/2015 at 10:30 am

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

    Reply
    • ccjou says:
      10/11/2015 at 6:26 am

      我在上文補充了說明。因為 \mathbf{y} 是任意非零向量，可設 \mathbf{y} 為標準單位向量 \mathbf{e}_i，則 \sum_{i=1}^nd_{ii}y_i^2=d_{ii}>0。

      Reply
      • ryan says:
        02/09/2020 at 12:40 am

        不好意思 請教老師 我也有相同的問題
        設 \mathbf{y} 為標準單位向量 \mathbf{e}_i
        那不就是原本只有一個\mathbf{y} (維度:n*1) 然後因為要有n個標準單位向量\mathbf{e}_i
        所以從原本一個向量變成n個單位向量?
        可以這樣令嗎?
        希望老師可以幫我解惑，麻煩老師了

        Reply
        • ryan says:
          02/09/2020 at 12:42 am

          不好意思 請教老師 我也有相同的問題
          設 y為標準單位向量 ei
          那不就是原本只有一個y (維度:n*1) 然後因為要有n個標準單位向量ei
          所以從原本一個向量變成n個單位向量?
          可以這樣令嗎?
          希望老師可以幫我解惑，麻煩老師了

          Reply
  2. Rongwei Sun says:
    12/23/2015 at 7:52 pm

    老師，請問您，複數矩陣可以Cholesky 分解嗎？

    Reply
  3. Murphy Ma says:
    07/17/2016 at 3:14 pm

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

    Reply
  4. ccjou says:
    07/19/2016 at 4:31 am

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

    Reply
    • Murphy Ma says:
      07/19/2016 at 10:10 am

      非常感謝老師的回復！抱歉我的問題沒有說清楚：問題是求解方程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，不明白為什麼？請老師指點迷津，謝謝！

      Reply
  5. 徐君 says:
    01/20/2017 at 11:13 am

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

    Reply
  6. Chloe says:
    02/07/2017 at 5:46 pm

    教授您好，最近在閱讀論文文獻時，發現若是要執行蒙地卡羅模擬轉化隨機數至有相關的標準常態隨機數時，通常會利用 Cholesky 分析方法將相關係數矩陣轉化成下三角矩陣乘以單行隨機變數向量(請見圖示 https://goo.gl/fUvudl 以及 https://goo.gl/KDT1MT )，然而我的疑惑在於，就計算來說採用下三角矩陣或其Transpose的上三角矩陣結果似乎是一樣的(乘上的都是隨機變數)，那教科書中採用下三角矩陣的原因背後有一定的規則嗎，還是其實用上三角矩陣來執行模擬也可以呢?

    Reply
  7. TK says:
    02/14/2017 at 7:10 pm

    很冒昧的想要請問老師，您在文中提及Cholesky的發明動機是為了解決測地計算問題，然而身為初學者的我還是不太清楚使用Cholesky分解出C=AT·A之後的用意以及Cholesky分解在其他實質地方的應用，能否煩請您舉例說明呢？

    Reply

Leave a Reply

Fill in your details below or click an icon to log in:

Gravatar
WordPress.com Logo

You are commenting using your WordPress.com account. ( Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. ( Log Out /  Change )

Cancel

Connecting to %s