PA=LU 分解

本文的閱讀等級：中級

令 $A$ 是一個 $n\times n$ 階可逆矩陣。LU 分解 $A=LU$ 是高斯消去法的一種表達形式，下三角矩陣 $L=[l_{ij}]$ 記錄消元過程使用的乘數，上三角矩陣 $U=[u_{ij}]$ 儲存約化結果，其中 $L$ 和 $U$ 的主對角元分別滿足 $l_{jj}=1$ 和 $u_{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}=1$ 且 $l_{ij}=0$ ， $i<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>j$ ， $l>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 k$ ， $i\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$ 是病態的。

9 Responses to PA=LU 分解

老羅 says:

11/17/2014 at 3:25 pm

請問周老師:
為什麼

$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}$ 這個運算要用那個關鍵字搜尋站內的文章呢?

- ccjou says:
  
  11/17/2014 at 4:47 pm
  
  如果沒有特別指明，線性代數社群習慣使用行向量 (column)，故
  $\mathbf{e}_3\mathbf{e}_1^T=\begin{bmatrix} 0\\ 0\\ 1\\ 0 \end{bmatrix}\begin{bmatrix} 1&0&0&0 \end{bmatrix}=\begin{bmatrix} 0&0&0&0\\ 0&0&0&0\\ 1&0&0&0\\ 0&0&0&0 \end{bmatrix}$
  此運算稱為outer product，但稱「外積」容易與cross product混淆。維基百科介紹：
  http://en.wikipedia.org/wiki/Outer_product
  
老羅 says:

11/17/2014 at 5:15 pm

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

- 老羅 says:
  
  11/18/2014 at 11:49 am
  
  參考，“特殊矩陣 (10)：基本矩陣” ，所以 $P_{ij}$ 排列矩陣是
  基本矩陣 $E_{1}$ 型。
  
cheng23 says:

01/15/2017 at 10:34 am

有幫助!謝謝!!!

mlrofcloud says:

04/21/2017 at 6:49 pm

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

$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 says:
  
  04/22/2017 at 10:04 am
  
  $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)$ 。
  
mlrofcloud says:

04/22/2017 at 11:36 am

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

Haodong Zhang says:

05/01/2018 at 11:29 pm

Thanks a lot for your explanation. It helps a lot!

However, I am still wondering the solution of Permutation matrix.
Especially the algorithm for calculating P. Looking forward to your reply.

Apropo, I’m from mainland. Therefore, I am really grateful if your answer is in traditional Chinese.
Cuz I could only read them rather than type and write.

Thanks!

	jianglong on Strassen 演算法──分治矩陣乘法
	jianglong on Strassen 演算法──分治矩陣乘法
	xmj on 內積的定義
	Ning ChingSan on 線性代數的第一堂課──矩陣乘法的定義
	momo on 兩岸線性代數用詞參照
	訪客 on 克拉瑪公式的簡易幾何證明

PA=LU 分解

9 Responses to PA=LU 分解

Leave a comment Cancel reply

搜尋(繁體中文或英文)

訊息看板

近期文章

線性代數專欄

其他主題專欄

每週問題

數據充分性問題

其他分類

Recent Comments

近期最多人點閱

分類

Archives

標籤雲

線代線上影音課程

線代學習網站

線代電子書

矩陣計算器

LaTeX

Blogroll

訂閱

閱讀導引

學習資源

專題探究

急救查詢

其他分頁

Meta

網路狀態

Blog Stats

PA=LU 分解

Share this:

9 Responses to PA=LU 分解

Leave a comment Cancel reply

搜尋(繁體中文或英文)

訊息看板

近期文章

線性代數專欄

其他主題專欄

每週問題

數據充分性問題

其他分類

Recent Comments

近期最多人點閱

分類

Archives

標籤雲

線代線上影音課程

線代學習網站

線代電子書

矩陣計算器

LaTeX

Blogroll

訂閱