矩阵分解算法

接触到通信中的下行链路预编码，其中很多线性预编码方法涉及到矩阵求逆，而矩阵求逆的计算复杂度随着维度增加而剧增，因此出现了很多基于矩阵分解的简化求逆方法。在总结预编码方法时发现当时学线性代数好像只接触了LU分解和SVD分解，故在此对常见的一些矩阵分解算法做个记录以便查询。后续也许会系统学习矩阵计算相关并记录。也许。。。资料来源于网络上各大佬的博客以及潘建瑜老师的《矩阵计算》讲义。侵删~

一、\(LU\)分解

1.1 \(LU\)分解定理

  \(LU\)分解实质就是将系数矩阵\(A\)分解为\(L\)和\(U\)两个矩阵的乘积。其中\(L\)是单位下三角矩阵，\(U\)为非奇异上三角矩阵。

  假定矩阵\(A\)存在\(LU\)分解，则原方程\(Ax=b\)就将转化为求解下面两个三角方程组，而三角方程组的求解非常容易求解。

• \[\left\{\begin{align}&Ly=b,\\&Ux=y.\end{align}\right. \]

  此外，\(LU\)分解也可用于对矩阵求逆，分别对三角矩阵\(L\)和\(U\)求逆比直接对\(A\)求逆的计算复杂度低很多。假设\(PA=LU\)，则

• \[A^{-1}=U^{-1}L^{-1}P \]

  并不是每个非奇异矩阵都存在\(LU\)分解。

定理（\(LU\)分解的存在性和唯一性）：矩阵\(A \in \mathbb{R}^{n \times n}\)存在\(LU\)分解（即存在单位下三角矩阵\(L\)和非奇异上三角矩阵\(U\)使得\(A=LU\)）的充要条件是\(A\)的所有顺序主子式都非奇异（行列式非零）。进一步地，若\(A\)存在 \(LU\)分解, 则分解是唯一的。

1.2 \(LU\)分解的实现方式

1.2.1 Gauss消元法

  Gauss消元法的本质是通过Gauss变换（矩阵初等变换的组合）来构造\(A\)的\(LU\)分解。给定一个矩阵

• \[A = \left[ {\begin{array}{*{20}{c}} a_{11}&a_{12}& \cdots &a_{1n}\\ a_{21}&a_{22}& \cdots &a_{2n}\\ \vdots &{}& \ddots &{}\\ a_{n1}&a_{n2}& \cdots &a_{nn} \end{array}} \right] \in \mathbb{R}{^{n \times n}} \]

• 第一步，假定\(a_{11}\ne 0\)，构造矩阵

  • \[{L_1} = \left[ {\begin{array}{*{20}{c}} 1&0&0& \cdots &0\\ l_{21}&1&0& \cdots &0\\ l_{31}&0&1& \cdots &0\\ \vdots & \vdots &{}& \ddots &{}\\ l_{n1}&0&0& \cdots &1 \end{array}} \right],其中~l_{i1} = \frac{a_{i1}}{a_{11}},i = 2,3, \ldots ,n. \]

  • 易知\(L_1\)的逆为

  • \[{L_1^{-1}} = \left[ {\begin{array}{*{20}{c}} 1&0&0& \cdots &0\\ -l_{21}&1&0& \cdots &0\\ -l_{31}&0&1& \cdots &0\\ \vdots & \vdots &{}& \ddots &{}\\ -l_{n1}&0&0& \cdots &1 \end{array}} \right] \]

  • 用\(L_1^{-1}\)左乘\(A\)，并将所得矩阵记为\(A^{(1)}\)，则

  • \[A^{(1)}=L_{1}^{-1} A\left[\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ 0 & a_{22}^{(1)} & \cdots & a_{2 n}^{(1)} \\ \vdots & \vdots & \ddots & \\ 0 & a_{n 2}^{(1)} & \cdots & a_{n n}^{(1)} \end{array}\right] \]

  • 即左乘\(L_1^{-1}\)后, \(A\)的第一列中除第一个元素外其它都变为0。

• 第二步，类似地，将上面的操作作用在\(A^{(1)}\)的子矩阵\(A^{(1)}(2:n,2:n)\)上，将其第一列除第一个元素外都变为0。也就是说，假定\(a_{22}^{(1)}\ne 0\)，构造矩阵

  • \[{L_1} = \left[ {\begin{array}{*{20}{c}} 1&0&0& \cdots &0\\ 0&1&0& \cdots &0\\ 0&l_{32}&1& \cdots &0\\ \vdots & \vdots &{}& \ddots &{}\\ 0&l_{n2}&0& \cdots &1 \end{array}} \right],其中{l_{i2}} = \frac{a_{i2}^{(1)}}{a_{22}^{(1)}},i = 3,4, \ldots ,n. \]

  • 用\(L_2^{-1}\)左乘\(A^{(1)}\)，所得矩阵记为\(A^{(2)}\)，则

  • \[A^{(2)}=L_{2}^{-1} A^{(1)}=L_{2}^{-1} L_{1}^{-1} A=\left[\begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \cdots & a_{1 n} \\ 0 & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & a_{2 n}^{(1)} \\ 0 & 0 & a_{33}^{(2)} & \cdots & a_{3 n}^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \\ 0 & 0 & a_{n 3}^{(2)} & \cdots & a_{n n}^{(2)} \end{array}\right]. \]

  • 依此类推，假定\(a_{kk}^{(k-1)}\ne 0(k=3,4,\dots,n-1)\)，则可以构造一系列矩阵\(L3,L4,\dots,L_{n-1}\)，使得

  • \[L_{n-1}^{-1} \cdots L_{2}^{-1} L_{1}^{-1} A=\left[\begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \cdots & a_{1 n} \\ 0 & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & a_{2 n}^{(1)} \\ 0 & 0 & a_{33}^{(2)} & \cdots & a_{3 n}^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \\ 0 & 0 & 0 & \cdots & a_{n n}^{(n-1)} \end{array}\right] \]

  • 为一个上三角矩阵，即所需要分解得到的非奇异上三角矩阵\(U\)，并记

  • \[L = {L_1}{L_2} \ldots {L_{n - 1}} = \left[ {\begin{array}{*{20}{c}} 1&0&0& \cdots &0\\ l_{21}&1&0& \cdots &0\\ l_{31}&l_{32}&1& \cdots &0\\ \vdots & \vdots &{}& \ddots &{}\\ l_{n1}&l_{n2}&l_{n3}& \cdots &1 \end{array}} \right] \]

  • 最终得到矩阵\(A\)的\(LU\)分解，即

  • \[A=LU \]

1.2.2 待定系数法

  假定\(A\)存在\(LU\)分解，即\(A=LU\)，或

• \[\left[\begin{array}{ccccc} a_{11} & a_{12} & a_{13} & \cdots & a_{1 n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2 n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3 n} \\ \vdots & & & \ddots & \vdots \\ a_{n 1} & a_{n 2} & a_{n 3} & \cdots & a_{n n} \end{array}\right]=\left[\begin{array}{ccccccccc} 1 & & & & \\ l_{21} & 1 & & & \\ l_{31} & l_{32} & 1 & & \\ \vdots & & & \ddots & \\ l_{n 1} & l_{n 2} & \cdots & l_{n, n-1} & 1 \end{array}\right]\left[\begin{array}{ccccc} u_{11} & u_{12} & u_{13} & \cdots & u_{1 n} \\ & u_{22} & u_{23} & \cdots & u_{2 n} \\ & & u_{33} & \cdots & u_{2 n} \\ & & & \ddots & \vdots \\ & & & & & u_{n n} \end{array}\right] \]

可以通过比较等式两边的元素来计算\(L\)和\(U\)中各元素的值。具体计算过程如下：

1. 比较等式两边的第一行，可得

   \[u_{1j}=a_{1j},\quad j=1,2,\dots,n \]

   再比较等式两边的第一列，可得

   \[a_{i1}=l_{i1}u_{11}\Rightarrow l_{i1}=a_{i1}/u_{11},\quad i=2,3,\dots,n \]

2. 比较等式两边的第二行，可得

   \[a_{2 j}=l_{21} u_{1 j}+u_{2 j} \Rightarrow u_{2 j}=a_{2 j}-l_{21} u_{1 j}, \quad j=2,3, \ldots, n \]

   再比较等式两边的第二列，可得

   \[a_{i2}=l_{i1} u_{12}+l_{i2}u_{22} \Rightarrow l_{i2}=(a_{i2}-l_{i1}u_{12})/u_{22}, \quad j=3,4, \ldots, n \]

3. 以此类推吗，第\(k\)步时，比较等式两边的第\(k\)行，可得

   \[u_{k j}=a_{k j}-\left(l_{k 1} u_{1 j}+\cdots+l_{k, k-1} u_{k-1, j}\right), \quad j=k, k+1, \ldots, n \]

   比较等式两边的第\(k\)列，可得

   \[l_{i k}=\left(a_{i k}-l_{i 1} u_{1 k}-\cdots-l_{i, k-1} u_{k-1, k}\right) / u_{k k}, \quad i=k+1, k+2, \ldots, n \]

1.3 条件更弱的LU分解定理（选主元\(LU\)分解）

  在\(LU\)分解算法中，称\(a_{kk}^{(k-1)}\)为主元，若\(a_{kk}^{(k-1)}=0\)，算法无法进行下去；若\(|a_{kk}^{(k-1)}|\)值非常小，则因为舍入误差导致结果误差非常大。可通过选主元来解决此问题。

  选主元需要引入置换矩阵，其基本性质如下：

设\(P\subset{\mathbb{R}^{n\times n}}\)为置换矩阵，\(X\subset{\mathbb{R}^{n\times n}}\)为任意矩阵，则

（1）\(PX\)相当于将\(X\)的行进行置换；\(XP\)相当于对\(X\)的列进行置换；

（2）\(P^{-1}=P^T\)，即\(P\)是正交矩阵；

（3）\(\det(P)=\pm1\);

（4）置换矩阵的乘积仍然是置换矩阵。

  另外选主元\(LU\)分解也不是一定存在的，其条件如下：

定理（选主元\(LU\)分解的存在性）：设\(A\subset{\mathbb{R}^{n\times n}}\)非奇异，则存在置换矩阵\(P_L,P_R\)，以及单位下三角矩阵\(L\)和非奇异上三角矩阵\(U\)，使得

\[P_LAP_R=LU \]

  注意置换矩阵有两种选择方式：

• “全主元Gauss消元法”（GECP）：使主元为剩下矩阵中绝对值最大。

• “部分选主元Gauss消元法”（GEPP）：使主元为第\(k\)列第\(k\)到第\(n\)个元素中绝对值最大，此时\(P_R^{(k)}=I\)，因此也成为列主元Gauss消元法

  （1）GECP 比 GEPP 更稳定, 但工作量太大, 在实际应用中通常使用 GEPP 算法.

  （2）GEPP 算法能保证 L 所有的元素的绝对值都不超过 1.

二、\(QR\)分解

  \(QR\)分解是将一个矩阵分解为一个正交矩阵（酉矩阵）和一个上三角矩阵的乘积。\(QR\)分解被广泛应用于线性最小二乘问题的求解和矩阵特征值的计算。

定理（\(QR\)分解）：设\(A\in{\mathbb{C}^{m\times n}(m\ge n)}\)，则存在一个单位列正交矩阵\(Q\in \mathbb{C}^{m\times n}(即Q^*Q=I_{n\times n})\)和一个上三角矩阵\(R\in \mathbb{C}^{n\times n}\)，使得

• \[A=QR \]

若\(A\)列满秩，则存在一个具有对角线元素的上三角矩阵\(R\)使得\(A=QR\)成立，且此时\(QR\)分解唯一，即\(Q\)和\(R\)都唯一。

推论：设\(A\in{\mathbb{C}^{m\times n}(m\ge n)}\)，且秩为\(l(0\le l\le n)\)，则存在一个置换矩阵\(P\)，使得

• \[AP=Q\left[\begin{align}&R_{11}~~R_{12}\\&~0~~~~~~~0\end{align}\right]_{n\times n}, \]

其中\(Q\in \mathbb{C}^{m\times n}\)单位列正交，\(R_{11}\in \mathbb{C}^{l\times l}\)是非奇异上三角矩阵。

三、\(SVD\)

  设 \(A\in \mathbb C^{m×n} (m ≥ n)\)，则$ A^∗A\in\mathbb C^{n×n} \(和\)AA^∗\in \mathbb C^{m×m}$都是 Hermite 半正定矩阵，且它们具有相同的非零特征值 (都是正实数)。

定理（\(SVD\)）：设任意矩阵 \(A\in\mathbb C^{m×n} (m ≥ n)\)，则存在酉矩阵\(U\in \mathbb C^{m\times m}\)和\(V\in \mathbb C^{n\times n}\)，使得

• \[U^{*} A V=\left[\begin{array}{l}\Sigma_r&0\\ 0&0\end{array}\right] \quad 或 \quad A=U\left[\begin{array}{l}\Sigma_r&0\\ 0&0\end{array}\right] V^{*} \]

其中\(\sum=\mathrm{diag}(\sigma_1,\sigma_2,\ldots,\sigma_r)\in \mathbb R^{r\times r}\)，且\(\sigma_1\ge\sigma_2\ge\ldots\ge\sigma_r\ge0\)。该分解称为\(A\)的奇异值分解（\(SVD\)）,而 \(\sigma_1,\sigma_2,\ldots,\sigma_r\) 称为\(A\)的奇异值（特征值），\(r\) 是矩阵 \(A\) 的秩。

四、\(Cholesky\)分解

  对称正定矩阵的基本性质：

设\(A\in\mathbb R^{n\times n}\).

• \(A\) 对称正定当且仅当 \(A\) 对称且所有特征值都是正的；
• \(A\) 对称正定当且仅当\(X^TAX\) 对称正定，其中 $X\in R_{n×n} $是一个任意的非奇异矩阵；
• 若 \(A\) 对称正定，则 \(A\) 的任意主子矩阵都对称正定；
• 若 \(A\) 对称正定，则 \(A\) 的所有对角线元素都是正的，且 \(\max_\limits{i\ne j} {|a_{ij}|} < \max_\limits{i} {|a_{ii}|}\)， 即绝对值最大的元素出现在对角线上。

  基于对称正定的\(Cholesky\)分解：

定理（\(Cholesky\)分解）：设\(A\in\mathbb R^{n\times n}\)对称正定，则存在唯一的对角线元素为正的下三角矩阵\(L\)，使得

\[A=LL^T \]

该分解称为\(Cholesky\)分解。

五、\(Jordan\)分解（\(Jordan\)标准型）

  设\(A\in\mathbb C^{n\times n}\)，有\(p\)个不同特征值，则存在非奇异矩阵\(X\in \mathbb C^{n\times n}\)，使得

• \[X^{-1} A X=\left[\begin{array}{cccc}J_{1} & & & \\ & J_{2} & & \\ & & \ddots & \\ & & & J_{q}\end{array}\right] \triangleq J \]

  其中\(J_i\)的维数等于\(\lambda_i\)的代数重数，且具有下面的结构：

• \[J_{i}=\left[\begin{array}{llll}J_{i 1} & & & \\ & J_{i 2} & & \\ & & \ddots & \\ & & & J_{i \nu_{i}}\end{array}\right], \quad J_{i k}=\left[\begin{array}{cccc}\lambda_{i} & 1 & & \\ & \ddots & \ddots & \\ & & \lambda_{i} & 1 \\ & & & \lambda_{i}\end{array}\right] \]

  这里的\(v_i\)为\(\lambda_i\)的几何重数，\(J_{ik}\)为（对应于\(\lambda_i\)的）\(Jordan\)块。

参考链接

• 几种矩阵分解算法： LU分解，Cholesky分解，QR分解，SVD分解，Jordan分解
• 矩阵的LU分解
• 矩阵计算