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}(-c)P_{kl}=E_{lj}(-c);若 i=l,則 P_{kl}E_{ij}(-c)P_{kl}=E_{kj}(-c);若 i\neq ki\neq l,則 P_{kl}E_{ij}(-c)P_{kl}=E_{ij}(-c)

 
根據以上討論結果,我們可以設計一個演算法來記錄高斯消去法的執行過程與結果,包括下三角矩陣 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 是病態的。

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

8 Responses to 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 說道:

    有幫助!謝謝!!!

  4. mlrofcloud 說道:

    教授您好,請教您,若是不藉由擴充矩陣的演算法,而是依循您演算法之前的說明來拆解基本矩陣乘積,卻無法求出正確分解,不知何處犯錯呢? 如下:

    E{_{43}}P{_{34}}E{_{42}}E{_{32}}P{_{23}}E{_{41}}E{_{31}}E{_{21}}A = U

    E{_{43}}P{_{34}}E{_{42}}P{_{34}}^{2}E{_{32}}P{_{34}}^{2}P{_{23}}E{_{41}}P{_{23}}^{2}E{_{31}}P{_{23}}^{2}E{_{21}}P{_{23}}^{2} A = U

    E{_{43}}\left \{P{_{34}}E{_{42}}P{_{34}}\right\}\left \{P{_{34}}E{_{32}}P{_{34}}\right\}P{_{34}}\left \{P{_{23}}E{_{41}}P{_{23}}\right\}\left\{P{_{23}}E{_{31}}P{_{23}}\right\}\left\{P{_{23}}E{_{21}}P{_{23}}\right \}P{_{23}} A = U

    E{_{43}}\left \{E{_{32}}\right \}\left \{E{_{42}}\right \}P{_{34}}\left \{E{_{41}}\right \}\left \{E{_{21}}\right \}\left \{E{_{31}}\right \}P{_{23}} A = U

    E{_{43}}E{_{32}}E{_{42}}P{_{34}}E{_{41}}P{_{34}}^{2}E{_{21}}P{_{34}}^{2}E{_{31}}P{_{34}}^{2}P{_{23}} A = U

    E{_{43}}E{_{32}}E{_{42}}\left \{P{_{34}}E{_{41}}P{_{34}}\right \}\left \{P{_{34}}E{_{21}}P{_{34}}\right \}\left \{P{_{34}}E{_{31}}P{_{34}}\right \}P{_{34}}P{_{23}} A = U

    E{_{43}}E{_{32}}E{_{42}}\left \{E_{31}\right \}\left \{E_{21}\right \}\left \{E_{41}\right \}P{_{34}}P{_{23}} A = U

    按照上述變換過程後,
    似乎應該 P = P{_{34}}P{_{23}}
    L = \left \{E{_{43}}E{_{32}}E{_{42}}E_{31}E_{21}E_{41}\right \}^{-1}
    但代入驗算後卻明顯有誤???
    懇請教授賜教,謝謝

    • ccjou 說道:

      E_{ij} 原來的乘數必須跟著移動,譬如,
      P_{34}E_{42}(-5)P_{34}=E_{32}(-5)
      你最後得到的E_{ij}是正確的,但不能用原來的乘數,例如,
      E_{21}(-3)\Rightarrow E_{31}(-3)\Rightarrow E_{41}(-3)

  5. mlrofcloud 說道:

    謝謝教授指導,我明白哪邊有誤了

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s