數值計算方法（二）·矩陣及其分解

作者：由腦洞的窗發表于歷史時間：2020-09-16

這個系列是數值計算方法的筆記，目錄請見：

這一部分很長，主線是各種“使解線性方程組

$\bm {Ax}=\bm b$

輕鬆一點的辦法”，把基本的矩陣分解法講完了。下面是它的小目錄：

在介紹矩陣分解之前，首先要知道一些特殊的矩陣和性質：

一些特殊矩陣

條件數

奇異值

然後才能講矩陣的常見分解方法，我會在標題裡寫出這種分解方法適用的矩陣：

方陣的Doolittle分解（Gauss消去、LDU分解、Crout分解）

方陣的PLU分解（帶列主元的Gauss消去法）

三對角矩陣的三角分解

對稱正定矩陣的Cholesky分解

滿秩陣（不一定是方陣，特別適用於病態陣）的QR分解

方陣（特別是複數方陣）的Schur分解

方陣（特別是不可對角化矩陣）的Jordan分解

任意矩陣（特別是長方陣）的奇異值分解

一些特殊矩陣

酉陣（實數範圍內即為正交陣）：

$\bm U^H\bm U=\bm I$

，具有F-範數不變性（實際上也具有奇異值不變性、2-範數不變性、2-條件數不變性）

正規矩陣：

$\bm A^H\bm A=\bm{AA}^H$

Hermite矩陣（實數範圍內即為對稱陣）：

$\bm A^H=\bm A$

反（斜）Hermite矩陣：

$\bm A^H=-\bm A$

作為了解

條件數

$\mathrm{cond}(\bm A)=||\bm A||\cdot||\bm A^{-1}||$

是矩陣的條件數

（使用運算元範數）

。

1。條件數越大，

$\bm{Ax}=\bm b$

裡

$\bm A,\bm b$

的一點變動

$\delta \bm A,\delta \bm b$

對解的影響越大，越病態。

2。根據所選擇的範數有1-條件數、2-條件數和∞-條件數，其中

$\mathrm{cond}_2(\bm A)=\sqrt{\dfrac{\lambda_{\max}(\bm A^H\bm A)}{\lambda_\min(\bm A^H\bm A)}}\\$

3。矩陣的條件數有如下性質：

$\begin{array}{rl} 1)&\mathrm{cond}(\bm A)\ge 1\\ 2)&\mathrm{cond}(\alpha\bm A)=\mathrm{cond}(\bm A)=\mathrm{cond}(\bm A^{-1}),\ \alpha\ne0\\ 3)&\mathrm{cond}(\bm {AB})\le\mathrm{cond}(\bm A)\mathrm{cond}(\bm B) \end{array}\\$

4。第3條的3）說明矩陣乘法算出來的結果，條件數最大可能到兩者條件數的積，實際上放大了病態性，但是對酉矩陣

$\bm U$

來說有：

$\mathrm{cond}_2(\bm U)=1,\mathrm{cond}_2(\bm U\bm A)=\mathrm{cond}_2(\bm A\bm U)=\mathrm{cond}_2(\bm A)\\$

說明乘酉矩陣不會增加病態性。

以上定理作為了解

奇異值

方陣有特徵值，外拓到長方陣，還有奇異值：

1。奇異值具有酉變換不變性

2。非零特徵值個數≥非零奇異值個數=秩

3。矩陣的2-範數是其最大奇異值，F-範數是奇異值平方和開根號（用F-範數酉不變性解釋）

4。矩陣的行列式=特徵值之積，矩陣的跡=特徵值之和

方陣的Doolittle分解（Gauss消去、LDU分解、Crout分解）

Doolittle分解裡的

$\bm L$

是對角元為1的下三角陣（單位下三角陣），LDU分解裡

$\bm L,\bm U$

是對角元為1的下、上三角陣，Crout分解裡

$\bm U$

是對角元為1的上三角陣。

手工計算一般是採用Gauss消去的方法：

1。用初等行變換，將

$\bm A$

的第1列的

$a_{11}$

下面的元素消成0（此時左乘

$\bm L_1$

），將第二列的

$a_{22}$

元素下面的元素消成0（此時再左乘

$\bm L_2$

）……最後得到了

$\bm U$

；

2。把剛才的行變換陣取逆再按出現的順序乘起來：

$\bm L=\bm L_1^{-1}\bm L_2^{-1}\cdots\\$

這一步有兩個技巧：

1。行變換陣取逆和單位下三角陣一樣，把1下面的元素全部取負號就行：

2。行變換陣的乘積就是把所有數字寫進來就行：

如果要寫演算法的話，把矩陣寫開，對應未知量設好，比如Doolittle分解：

對應就有：

Doolittle分解解出的兩個矩陣存在且唯一的充分條件是

階方陣

$\bm A$

的前

個順序主子式均不為0。

如果順序主子式有為0的情況，分解可能還存在，但是即使存在也不唯一。比如下面這個矩陣，第二個順序主子式是0，但是仍然可以分解：

$\left(\begin{array}{cccc} 1 & 2 & 3 & 4\\ 0 & 0 & 2 & 3\\ 0 & 0 & 1 & 2\\ 0 & 0 & 0 & 1 \end{array}\right)=\left(\begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{array}\right)\left(\begin{array}{cccc} 1 & 2 & 3 & 4\\ 0 & 0 & 2 & 3\\ 0 & 0 & 1 & 2\\ 0 & 0 & 0 & 1 \end{array}\right)\\$

Doolittle分解的原理雖然是Gauss消去，但是計算量大大減小了：

1。Doolittle分解要使

$\bm A=\bm{LU}$

，Gauss消去法要使

$(\bm A,\bm b)\rightarrow(\bm U,\bm c)$

，以乘法次數計，這兩種計算量都是

$o(\dfrac{n^3}3)$

（思考：為什麼生成

$\bm L$

並不需要計算量？）

2。解上下三角陣×向量=常向量計算量是

$o(\dfrac{n^2}{2}) \ll o(\dfrac{n^3}{3})$

3。Doolittle矩陣好在係數矩陣不變時，不用重新分解

方陣的PLU分解（帶列主元的Gauss消去法）

帶列主元的Gauss消去法和PLU分解在每一步多了個選主元的步驟：第

步左乘

$\bm L_k$

實際上是讓

$k+1,\cdots,n$

行分別除以

$a_{kk}$

，這個數如果太小，誤差就會相當大，所以要把第

行下具有最大

$|a_{ik}|\ (i=k+1,\cdots,n)$

的行交換上來

無論是手工計算還是程式設計，都可以採用先交換主元，再利用新的矩陣

$\bm A_{temp}$

進行LU分解的辦法，比如對矩陣

$\bm A=\left(\begin{array}{ccc} 0 & -6 & -1\\ 1 & 2 & 2\\ 2 & -2 & 1 \end{array}\right)\\$

先將1與3兩行交換，然後新的2與3兩行再交換，這兩步的操作分別對應左乘矩陣

$\bm P_1=\left(\begin{array}{ccc} & & 1\\ & 1 & \\ 1 & & \end{array}\right),\ \bm P_2=\left(\begin{array}{ccc} 1 & & \\ & & 1 \\ & 1 & \end{array}\right) \\$

由此已經產生了

$\bm P=\bm{P}_1\bm{P}_2=\left(\begin{array}{ccc} & & 1\\ 1 & & \\ & 1 & \end{array}\right),\ \bm A_{temp}=\bm {PA}=\left(\begin{array}{ccc} 2 & -2 & 1\\ 0 & -6 & -1\\ 1 & 2 & 2 \end{array}\right) \\$

然後對

$\bm A_{temp}$

進行LU分解：

$\bm L=\left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ \frac{1}{2} & - \frac{1}{2} & 1 \end{array}\right),\ \bm U=\left(\begin{array}{ccc} 2 & -2 & 1\\ 0 & -6 & -1\\ 0 & 0 & 1 \end{array}\right) \\$

三對角矩陣的三角分解

三對角矩陣的三角分解是一種特殊的LU分解，它的展開形式是

它的計算首先要將

直接替換成

，然後按照上圖藍箭頭的順序，依次計算。

對角列佔優，側對角線上沒有零元素是三角分解存在且唯一的充分條件。

三對角矩陣特殊在它能用四個向量

$\bm a,\bm b,\bm c,\bm f$

儲存方程裡

$\bm{Ax}=\bm f$

裡需要的所有元素，由此計算出的

仍然一一存放在對應的位置。

對稱正定矩陣的Cholesky分解

Cholesky分解

$\bm A=\bm{LL}^T$

僅適用於

對稱正定矩陣

，其中

$\bm L$

是一個下三角陣

分解一定存在，如果規定

$\bm L$

的對角元是正數，則分解也是唯一的

Cholesky法計算量是Gauss消去法的一半

滿秩陣（不一定是方陣，特別實用於病態陣）的QR分解

矩陣QR分解

$\bm A_{m\times n}=\bm Q_{m\times m}\bm R_{m\times n}$

適用於

滿秩陣

，其中

$\bm Q$

是一個正交陣（思考：最後一步由

$\bm Q^{-1}\bm A=\bm R$

求

$\bm Q$

可以怎麼求？），

$\bm R$

是一個對角元都大於0的上三角陣

由於分解後是乘正交陣

$\bm Q$

，條件數不變，適合分解病態性比較大的矩陣

我們要計算QR分解，就要利用Householder矩陣，即當

$\bm\omega\ne \bm 0$

時的以下矩陣：

$\bm H(\bm\omega)=\bm I-\frac{2}{\bm\omega^T\bm\omega}\bm\omega\bm\omega^T\\$

Householder矩陣的特性是

對稱正交，左乘向量進行鏡面變換

：

（4）把x鏡面變換到座標軸的一個正方向上了

要計算QR分解，就要用Householder矩陣形成的正交陣不斷左乘

$\bm A$

，利用Householder矩陣的性質（4）讓

$\bm A$

的第一列投影到

$\bm e_1=(\underbrace{1,0,\cdots,0}_{n個})^T$

的方向上，第二列從第二個元素開始再投影到

$\bm e_1=(\underbrace{1,0,\cdots,0}_{n-1個})^T$

的方向上……從而達到這樣的效果：

$\bm Q_1\bm A=\bm H(\bm\omega_1)\bm A=\left(\begin{array}{cc} a_{11} & \bm a_1^T\\ \bm 0 & \bm A_{22} \end{array} \right) \\ \bm Q_2(\bm Q_1\bm A)=\left(\begin{array}{cc} 1 & \bm 0^T\\ \bm 0 & \bm H(\bm\omega_2) \end{array} \right) \left(\begin{array}{cc} a_{11} & \bm a_1^T\\ \bm 0 & \bm A_{22} \end{array} \right) = \left(\begin{array}{cc} a_{11} & \bm a_1^T\\ \bm 0 & \bm H(\bm\omega_2\bm)\bm A_{22} \end{array} \right)\\[22pt]\cdots\\$

比如分解下面這個矩陣：

$\bm A=\left( \begin{array}{ccc} 0 & 4 & 1\\ 1 & 1 & 1\\ 0 & 3 & 2 \end{array} \right)\\$

它的完整解如下圖：

對滿秩陣，這種分解是存在的，但並不唯一，當限制

$\bm R$

具有正對角元時才唯一

方陣（特別是複數方陣）的Schur分解

矩陣的Schur分解

$\bm A=\bm {URU}^H$

對所有復空間方陣

$\bm A$

成立

其中

$\bm R$

為上三角陣（即使

$\bm A$

是實的，

$\bm {U,R}$

也可能是復的），

和#FormatImgID_125#具有相同特徵值

。

$\bm U$

是酉陣

透過Schur分解可以得到以下結論：

$\bm A是正規矩陣\leftrightarrow \bm A=\bm {UDU}^H，\bm{D}是對角陣\\ \bm A是Hermite陣\leftrightarrow \bm A=\bm {UDU}^H，\bm{D}是實的對角陣\\ \bm A是反Hermite陣\leftrightarrow \bm A=\bm {UDU}^H，\bm{D}是對角元為0或純虛數的對角陣\\ \bm A是酉陣\leftrightarrow \bm A=\bm {UDU}^H，\bm{D}是對角元模均為1的對角陣\\$

這個主要是理論分析用的，計算方法略。

方陣（特別是不可對角化矩陣）的Jordan分解

Jordan分解

$\bm A=\bm{TJT}^{-1}$

適用於任何方陣，其中

$\bm J$

是Jordan陣，

$\bm T$

是一個可逆陣

Jordan分解首先需要知道Jordan塊和Jordan陣，其中Jordan塊在階數為

，對角元始終為

$\lambda$

的對角陣的基礎上，上對角線上填了1，通常寫作

$\bm J_n(\lambda)$

，如：

Jordan陣

$\bm J$

是對角分塊矩陣，對角位置填入Jordan塊。對角矩陣也是Jordan陣，因為它的階數為1。

求解Jordan分解需要先求

$\bm J$

，然後求正交陣

$\bm T$

。求

$\bm J$

的時候需要求解每一個特徵值

$\lambda_i$

的代數重數

和幾何重數

$\alpha_i$

1。代數重數

=特徵方程

$\mathrm{det}(\lambda\bm I-\bm A)=0$

中根

$\lambda_i$

的重數=Jordan標準型對角元中

$\lambda_i$

的個數

2。幾何重數

$\alpha_i=n-\mathrm{rank}(\lambda_i\bm I-\bm A)$

$\lambda_i$

對應的線性無關特徵向量的個數=

$\lambda_i$

對應的Jordan塊的個數

3。具體階數為

的Jordan塊有

$r_{l+1}+r_{l-1}-2r_l$

個，其中

$r_0=n,r_l=\mathrm{rank}(\lambda_i\bm I-\bm A)^l$

4。矩陣可以進行對角化的充要條件是每一個特徵值

$\lambda_i$

都有

$m_i=\alpha_i$

求解

$\bm T$

要將

$\bm A=\bm{TJT}^{-1}$

寫成

$\bm{AT}=\bm{TJ}$

，再按Jordan塊分解開：

Jordan鏈需要一個鏈首

$\bm t_1^i$

，它要在特徵值

對應的線性無關特徵向量裡找，需要注意，直接求出來的線性無關特徵向量可能都不能做鏈首，要能根據以下式子，取

求出非零的向量：

然後再根據該式，逐步求解整個鏈中的向量，直至

達到該Jordan塊的階數

。這個過程需要懂解齊次方程組和非齊次方程，見

［1］

，剛才的例子繼續求解如下：

不計Jordan塊排列順序，這種分解存在且

$\bm J$

是唯一的（即使

$\bm J$

唯一，

$\bm T$

仍有可能不同）

任意矩陣（特別是長方陣）的奇異值分解

$\bm A_{m\times n}=\bm U_{m\times m} \left(\begin{array}{cc} \bm\Sigma_r & \bm O\\ \bm O & \bm O \end{array}\right)_{m\times n}\bm V^H_{n\times n} 或\bm A=\bm{U_1\Sigma V_1}^H$

，前者是滿的奇異值分解，後者是約化的奇異值分解

式中

$\bm{U,V}$

都是酉陣（思考：為什麼求其分量的時候必須

標準化？

），

$\bm\Sigma$

是含有

$\bm A$

的所有非0奇異值的對角陣

求解奇異值分解要根據以下步驟：

計算

$\bm A_{m\times n}$

的奇異值，組成矩陣

$\bm \Sigma_{r\times r}$

已知

$\bm V^H(\bm A^H\bm A)\bm V=\left( \begin{array}{cc} \bm\Sigma^2_{r\times r} &\bm O\\ \bm O & \bm O_{(n-r)\times (n-r)} \end{array} \right)_{n\times n}\\$