算法：矩阵 LUP 分解

本文着笔于矩阵 LUP 分解算法，以及利用矩阵的 LUP 分解来解线性方程组、求矩阵对应行列式的值、求逆矩阵。

对于矩阵的定义代码如下：

矩阵 A 的 LUP 分解即为找到下三角矩阵 L、上三角矩阵 U、行置换矩阵 P 使得 PA = LU ，实现方法为高斯消元。

E = I = ( 1 0 ⋯ 0 0 1 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 1 ) , e i j = { 1 i = j 0 i ≠ j \begin{matrix} E \end{matrix} = \begin{matrix} I \end{matrix} = \begin{pmatrix} 1 ; 0 ; \cdots ; 0 \\ 0 ; 1 ; \cdots ; 0 \\ \vdots ; \vdots ; \ddots ; \vdots \\ 0 ; 0 ; \cdots ; 1 \end{pmatrix}, e_{ij} = \begin{cases} 1 ; i = j \\0 ; i \neq j \end{cases} E​=I​=⎝⎜⎜⎜⎛​10⋮0​01⋮0​⋯⋯⋱⋯​00⋮1​⎠⎟⎟⎟⎞​,eij​={10​i=ji̸​=j​

初始矩阵 P = E。对于矩阵 A，我们可以将它分解成这样两个矩阵的乘积（这里直接借用了 Introduction to Algorithms 里的描述） A = ( a 11 ω T ν A ′ ) = ( 1 0 ν / a 11 I n − 1 ) ( a 11 ω T 0 A ′ − ν ω T / a 11 ) \begin{matrix} A \end{matrix} = \begin{pmatrix} a_{11} ; \omega^T \\ \nu ; A' \end{pmatrix} = \begin{pmatrix} 1 ; 0 \\ \nu / a_{11} ; I_{n-1} \end{pmatrix} \begin{pmatrix} a_{11} ; \omega^T \\ 0 ; A' - \nu \omega^T / a_{11} \end{pmatrix} A​=(a11​ν​ωTA′​)=(1ν/a11​​0In−1​​)(a11​0​ωTA′−νωT/a11​​)

这样，将这个规模为 n 的问题分解为一个规模为 n - 1 的子问题，迭代之即可求得矩阵 PA 的 LU 分解。最后得到的矩阵 U 相当于高斯消元后的结果，因为求矩阵 U 的实质是用当前规模矩阵的第一行的相应倍数减其它行，将其它行首列元素消为 0。

为了避免除数为零或过小引起精度问题，在迭代到矩阵第 i 行时，我们可以找到 j 使得 ∣ A j , i ∣ \vert A_{j, i} \vert ∣Aj,i​∣ 尽可能大，其中 i ≤ j ≤ n i \leq j \leq n i≤j≤n。当 i 不等于 j 时，交换矩阵 PA 的第 i 行与第 j 行，相当于矩阵 A 不变，矩阵 P 的第 i 行与第 j 行交换。

对于矩阵的 LUP 分解代码如下：

时间复杂度 O ( n 2 ) O(n^2) O(n2)。（若算上矩阵 LUP 分解，总时间复杂度 O ( n 3 ) O(n^3) O(n3)）

由行列式性质可知 ∣ P ∣ ∣ A ∣ = ∣ L ∣ ∣ U ∣ \begin{vmatrix} P \end{vmatrix} \begin{vmatrix} A \end{vmatrix} = \begin{vmatrix} L \end{vmatrix} \begin{vmatrix} U \end{vmatrix} ∣∣​P​∣∣​∣∣​A​∣∣​=∣∣​L​∣∣​∣∣​U​∣∣​

其中，L 为对角线元素均为 1 的下三角矩阵，故其行列式值为 1；U 为上三角矩阵，其行列式值为对角线元素之积；P 初始行列式值为 1，每对 P 进行一次行交换，相当于行列式值取其相反数。

在具体实现过程中，行列式值的初值为 1 或 -1 在 LUP 分解过程中决定，此部分代码见上文。

对于矩阵的行列式值求解代码如下：