PA=LU 分解

本文的閱讀等級:中級

A 是一個 n\times n 階可逆矩陣。LU 分解 A=LU 是高斯消去法的一種表達形式,下三角矩陣 L=[l_{ij}] 記錄消元過程使用的乘數,上三角矩陣 U=[u_{ij}] 儲存約化結果,其中 LU 的主對角元分別滿足 l_{jj}=1u_{jj}\neq 0 (見“LU 分解”)。不過,並非每一可逆矩陣都存在 LU 分解。在執行高斯消去法的化簡程序中,若 0 出現在軸元 (pivot) 位置,即 (j,j) 元,列取代運算便無法消去軸元底下各元,這時標準 LU 分解不復存在。可逆矩陣 A 的軸元總數 (即 \mathrm{rank}A) 等於 n,透過列交換運算必能獲得一 (非零) 軸元,所以仍可繼續下一階段的化簡步驟。從實際面來看,縱使未發生「零軸元」的情況 (軸元所在位置為零),為了避免因不當使用消去法而引發災難,部分軸元法 (partial pivoting) 也會適時地交換列 (見“特殊矩陣 (12):對角佔優矩陣”)。見下例 (取自[1]):

A=\begin{bmatrix} 0.0001&1\\ 1&1 \end{bmatrix}=\begin{bmatrix} 1&0\\ 10,000&1 \end{bmatrix}\begin{bmatrix} 0.0001&1\\ 0&-9999 \end{bmatrix}=LU

U 同時包含數值很小和很大的主對角元可知 U 不具數值穩定性,也就是說,標準 LU 分解將原本是良置的 (well-conditioned) A 分解出病態的 (ill-conditioned) U[2]。問題發生的根源在於軸元的數值太小,交換 A 的列以取得數值較大的軸元是最簡單的解決方法。所謂部分軸元法是指在每一次消元步驟進行之前,先搜尋軸元所在位置及其底下的所有元,從其中選出數值 (絕對值) 最大的元,並交換兩列將最大元置於軸元位置。回到前例,令 P=\begin{bmatrix} 0&1\\ 1&0 \end{bmatrix}PA 的 LU 分解可有效改善數值穩定性,結果如下:

PA=\begin{bmatrix} 1&1\\ 0.0001&1 \end{bmatrix}=\begin{bmatrix} 1&0\\ 0.0001&1 \end{bmatrix}\begin{bmatrix} 1&1\\ 0&0.9999 \end{bmatrix}=LU

我們可能認為在計算 LU 分解之前,必須先確定 A 的列排序,即排列矩陣 P。這是錯覺,其實無此必要。但如果化簡 A 的過程中混合了列交換與消元運算,還能保證可逆矩陣 A 必可化約成 PA=LU 形式嗎?答案是肯定的,本文將探討這個問題,並給出一個實用的演算法。

 
在高斯消去法將 A 約化成上三角矩陣 U 的過程中,如欲消去 (i,j) 元,i>j,令列 i 減去列 j 通乘一乘數 c,此運算可表示為消元基本矩陣 E_{ij}(-c)=I-c\mathbf{e}_i\mathbf{e}_j^T (見“特殊矩陣 (10):基本矩陣”)。如欲交換列 i 和列 j,則使用排列基本矩陣 P_{ij}=I+(\mathbf{e}_i-\mathbf{e}_j)(\mathbf{e}_j-\mathbf{e}_i)^T。例如,

E_{31}(-2)=\left[\!\!\begin{array}{rccc} 1&0&0&0\\ 0&1&0&0\\ -2&0&1&0\\ 0&0&0&1 \end{array}\!\!\right],~P_{23}=\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1 \end{bmatrix}

下面我們用一個例子來說明容許列交換運算的高斯消去法。為便利說明,我們不使用部分軸元法,僅在遇見「零軸元」時才執行列交換;同時為容易辨識,在消元產生的零元下加註一底線。過程如下:

\begin{aligned} A=&\left[\!\!\begin{array}{crrc} 2&-3&4&2\\ 6&-9&12&5\\ 4&-5&10&5\\ 2&2&11&9 \end{array}\!\!\right]\xrightarrow[]{~E_{41}(-1)E_{31}(-2)E_{21}(-3)~}\left[\!\!\begin{array}{crcr} 2&-3&4&2\\ \underline{0}&0&0&-1\\ \underline{0}&1&2&1\\ \underline{0}&5&7&7 \end{array}\!\!\right]\\ \xrightarrow[]{~P_{23}~}&\left[\!\!\begin{array}{crcr} 2&-3&4&2\\ \underline{0}&1&2&1\\ \underline{0}&0&0&-1\\ \underline{0}&5&7&7 \end{array}\!\!\right]\xrightarrow[]{~E_{42}(-5)E_{32}(0)~}\left[\!\!\begin{array}{crrr} 2&-3&4&2\\ \underline{0}&1&2&1\\ \underline{0}&\underline{0}&0&-1\\ \underline{0}&\underline{0}&-3&2 \end{array}\!\!\right]\\ \xrightarrow[]{~P_{34}~}&\left[\!\!\begin{array}{crrr} 2&-3&4&2\\ \underline{0}&1&2&1\\ \underline{0}&\underline{0}&-3&2\\ \underline{0}&\underline{0}&0&-1 \end{array}\!\!\right]\xrightarrow[]{~E_{43}(0)~}\left[\!\!\begin{array}{crrr} 2&-3&4&2\\ \underline{0}&1&2&1\\ \underline{0}&\underline{0}&-3&2\\ \underline{0}&\underline{0}&\underline{0}&-1 \end{array}\!\!\right]=U,\end{aligned}

其中 E_{32}(0)=E_{43}(0)=I 是不需要執行的虛擬運算。將以上化簡過程表示成矩陣乘積:

E_{43}P_{34}E_{42}E_{32}P_{23}E_{41}E_{31}E_{21}A=U

想像吾人有未卜先知能力,在消元步驟前提早將 A 的各列排序妥當,即得 P_{34}P_{23}A。不過如此一來,基本矩陣 E_{ij} 勢必跟著改變,設為 \tilde{E}_{ij},就有

\tilde{E}_{43}\tilde{E}_{42}\tilde{E}_{32}\tilde{E}_{41}\tilde{E}_{31}\tilde{E}_{21}P_{34}P_{23}A=U

P=P_{34}P_{23}L=(\tilde{E}_{43}\tilde{E}_{42}\tilde{E}_{32}\tilde{E}_{41}\tilde{E}_{31}\tilde{E}_{21})^{-1},上式即可表示為

PA=\tilde{E}_{21}^{-1}\tilde{E}_{31}^{-1}\tilde{E}_{41}^{-1}\tilde{E}_{32}^{-1}\tilde{E}_{42}^{-1}\tilde{E}_{43}^{-1}U=LU

接下來只要確定 \tilde{E}_{ij} 也是消元基本矩陣,L 即符合要求:l_{ii}=1l_{ij}=0i<j。因為 P_{ij}^2=I,可得

\begin{aligned}P_{23}E_{41}E_{31}E_{21}&=P_{23}E_{41}P_{23}^2E_{31}P_{23}^2E_{21}P_{23}^2\\ &=(P_{23}E_{41}P_{23})(P_{23}E_{31}P_{23})(P_{23}E_{21}P_{23})P_{23}, \end{aligned}

其中

\begin{aligned}P_{23}E_{21}P_{23}&=\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1 \end{bmatrix}\left[\!\!\begin{array}{rccc} 1&0&0&0\\ -3&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{array}\!\!\right]\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1 \end{bmatrix}=\left[\!\!\begin{array}{rccc} 1&0&0&0\\ 0&1&0&0\\ -3&0&1&0\\ 0&0&0&1 \end{array}\!\!\right]\\ P_{23}E_{31}P_{23}&=\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1 \end{bmatrix}\left[\!\!\begin{array}{rccc} 1&0&0&0\\ 0&1&0&0\\ -2&0&1&0\\ 0&0&0&1 \end{array}\!\!\right]\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1 \end{bmatrix}=\left[\!\!\begin{array}{rccc} 1&0&0&0\\ -2&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{array}\!\!\right]\\ P_{23}E_{41}P_{23}&=\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1 \end{bmatrix}\left[\!\!\begin{array}{rccc} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ -1&0&0&1 \end{array}\!\!\right]\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&1&0&0\\ 0&0&0&1 \end{bmatrix}=\left[\!\!\begin{array}{rccc} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ -1&0&0&1 \end{array}\!\!\right].\end{aligned}

兩件事情值得注意。第一,上例中 P_{23}E_{i1}P_{23}i=2,3,4,皆為消元基本矩陣。一般情況也成立,理由如下:列交換運算總是在消去軸元 (j,j) 底下各元後啟動,令 k>jl>j,則 \mathbf{e}_j^TP_{kl}=\mathbf{e}_j^T,可知

P_{kl}E_{ij}P_{kl}=P_{kl}(I-c\mathbf{e}_i\mathbf{e}_j^T)P_{kl}=P_{kl}^2-cP_{kl}\mathbf{e}_i\mathbf{e}_j^TP_{kl}=I-c(P_{kl}\mathbf{e}_i)\mathbf{e}_j^T

也是一消元基本矩陣。第二,E_{ij} 所含的乘數 -c 隨列交換運算 P_{kl} 移動位置。若 i=k,則 P_{kl}E_{ij}P_{kl}=E_{lj};若 i=l,則 P_{kl}E_{ij}P_{kl}=E_{kj};若 i\neq ki\neq l,則 P_{kl}E_{ij}P_{kl}=E_{ij}

 
根據以上討論結果,我們可以設計一個演算法來記錄高斯消去法的執行過程與結果,包括下三角矩陣 L=[l_{ij}] 所儲存的乘數 l_{ij}i>j,排列矩陣 P 的列排序,以及上三角矩陣 U=[u_{ij}]i\le j。實際作法是引進一增廣矩陣,在給定矩陣 A 的最右側加入一行代表列排序,設初始值為 (1, 2, 3, 4),並於消元過程產生的零元位置 (即上例加註底線者) 填入所使用的乘數,同樣也加註一底線以利辨識。請特別注意:執行消元運算時,我們不計算代表列排序的最右行與儲存的乘數,但列交換運算則不受此限。完整過程如下:

\begin{aligned}  A=&\left[\!\!\begin{array}{crrccc} 2&-3&4&2&\vline&1\\ 6&-9&12&5&\vline&2\\ 4&-5&10&5&\vline&3\\ 2&2&11&9&\vline&4 \end{array}\!\!\right]\xrightarrow[]{E_{41}(-1)E_{31}(-2)E_{21}(-3)}\left[\!\!\begin{array}{crcrcc} 2&-3&4&2&\vline&1\\ \underline{3}&0&0&-1&\vline&2\\ \underline{2}&1&2&1&\vline&3\\ \underline{1}&5&7&7&\vline&4 \end{array}\!\!\right]\\ \xrightarrow[]{P_{23}}&\left[\!\!\begin{array}{crcrcc} 2&-3&4&2&\vline&1\\ \underline{2}&1&2&1&\vline&3\\ \underline{3}&0&0&-1&\vline&2\\ \underline{1}&5&7&7&\vline&4 \end{array}\!\!\right]\xrightarrow[]{E_{42}(-5)E_{32}(0)}\left[\!\!\begin{array}{crrrcc} 2&-3&4&2&\vline&1\\ \underline{2}&1&2&1&\vline&3\\ \underline{3}&\underline{0}&0&-1&\vline&2\\ \underline{1}&\underline{5}&-3&2&\vline&4 \end{array}\!\!\right]\\ \xrightarrow[]{P_{34}}&\left[\!\!\begin{array}{crrrcc} 2&-3&4&2&\vline&1\\ \underline{2}&1&2&1&\vline&3\\ \underline{1}&\underline{5}&-3&2&\vline&4\\ \underline{3}&\underline{0}&0&-1&\vline&2 \end{array}\!\!\right]\xrightarrow[]{E_{43}(0)}\left[\!\!\begin{array}{crrrcc} 2&-3&4&2&\vline&1\\ \underline{2}&1&2&1&\vline&3\\ \underline{1}&\underline{5}&-3&2&\vline&4\\ \underline{3}&\underline{0}&\underline{0}&-1&\vline&2 \end{array}\!\!\right].\end{aligned}

從最後得到的增廣矩陣即可寫出 PA=LU 分解。最末行顯示列排序為 (1,3,4,2),對應排列矩陣

P=\begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ 0&0&0&1\\ 0&1&0&0 \end{bmatrix}

下三角矩陣 L 和上三角矩陣 U 則分別為

L=\begin{bmatrix} 1&0&0&0\\ 2&1&0&0\\ 1&5&1&0\\ 3&0&0&1 \end{bmatrix},~U=\left[\!\!\begin{array}{crrr} 2&-3&4&2\\ 0&1&2&1\\ 0&0&-3&2\\ 0&0&0&-1 \end{array}\!\!\right]

 
註解
[1] Gene H. Golub and Charles F. Van Loan, Matrix Computations, 2nd edition, 1989.
[2] 條件數 (condition number) 是一個判斷矩陣是否為良置或病態的標準方法 (詳見“條件數)。對於一 n\times n 階可逆矩陣 A,令 \sigma_{\max}(A) 為最大奇異值,\sigma_{\min}(A) 為最小奇異值。矩陣 A 的條件數定義為 \kappa(A)=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}。若 \kappa(A) 接近 1,則 A 是良置系統;若 \kappa(A) 很大,則 A 是病態系統。上例中,A=\begin{bmatrix} 0.0001&1\\ 1&1 \end{bmatrix} 的條件數為 \kappa(A)=1.618/0.618\approx 2.618,可知 A 是良置的;U=\begin{bmatrix} 0.0001&1\\ 0&-9999 \end{bmatrix} 的條件數為 \kappa(U)=9999/0.0001\approx 10^8,可知 U 是病態的。

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

5 則回應給 PA=LU 分解

  1. 老羅 說:

    請問周老師:
    為什麼

    e_{3}e_{1}^{T}
    等於
    \begin{bmatrix}  0& 0 & 0 &0 \\   0& 0 & 0 &0 \\   1& 0 & 0 &0 \\   0& 0 & 0 &0  \end{bmatrix}
    , 而
    e_{3}
    \begin{bmatrix} 0 & 0 & 1 & 0  \end{bmatrix}
    的意思嗎?

    \begin{bmatrix}  0& 0 & 0 &0 \\   0& 0 & 0 &0 \\   1& 0 & 0 &0 \\   0& 0 & 0 &0  \end{bmatrix}
    比較易懂,有橫式表達法嗎?
    看到e_{3}e_{1}^{T},腦袋就當機了。
    同理
    (e_{i}-e_{j})(e_{j}-e_{i})^{T} 這個運算要用那個關鍵字搜尋站內的文章呢?

  2. 老羅 說:

    謝謝 ,原來是外積,\bigotimes。之前一直以為是內積。

  3. cheng23 說:

    有幫助!謝謝!!!

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s