LUP分解方法，矩阵求逆

LU分解
LU、LUP分解
LUP求解线性方程组

LU矩阵分解

LU分解可以将一个矩阵 A\mathbf AA 分解为一个单位下三角矩阵 L\mathbf LL 和一个上三角矩阵 U\mathbf UU 的乘积。单位下三角单位矩阵主对角线上的元素为 1 。有如下关系：
A=LU(1)\mathbf A = \mathbf L \mathbf U \tag1A=LU(1)

该分解可以采用高斯消元法来进行计算。若想将给定矩阵 A\mathbf AA 分解为下三角矩阵 L\mathbf LL 和上三角矩阵 U\mathbf UU ，一个思路就是通过一系列的初等变换将 (1)式化为上三角矩阵，且保证这些变换的乘积是一个下三角。

在程序实现上，（1）式中的矩阵依次对应到下方三个矩阵以及两个三角矩阵相乘后元素之间的关系：

图1. LU分解过程

核心代码
算法如下：

for (int i = 0; i < N - 1; i++)
{double u = U[i * N + i];for (int j = i + 1; j < N; j++){U[j * N + i] /= u;for (int k = i + 1; k < N; k++){U[j * N + k] -= U[i * N + k] * U[j * N + i];}//Of for k}// Of for j
}

数组 U 初始存储矩阵 A\mathbf AA，经过计算后存储上三角矩阵和除主对角线以外的下三角矩阵的元素。

在代码的实现上和图1 的过程有一些区别。代码实现更像是一点一点地把矩阵 A\mathbf AA 给稀释的感觉。

第一行 u11:u14u_{11} : u_{14}u11:u14 可以直接解出。例如在第一轮外循环后，下三角矩阵第一列被解出，上三角第二行被解出。数组 U 中存储的矩阵元素为：
[u11u12u13u14l21u22u23u24l31l32u22l32u23+u33l32u24+u34l41l42u22l42u23+l43u33l42u24+l43u34+u44]\begin{bmatrix} u_{11} & u_{12} & u_{13} & u_{14}\\ l_{21} & u_{22} & u_{23} & u_{24}\\ l_{31} & l_{32}u_{22} & l_{32}u_{23}+u_{33} &l_{32}u_{24}+u_{34}\\ l_{41} & l_{42}u_{22} &l_{42}u_{23}+l_{43}u_{33} & l_{42}u_{24}+l_{43}u_{34}+u_{44} \end{bmatrix}⎣⎢⎢⎡u11l21l31l41u12u22l32u22l42u22u13u23l32u23+u33l42u23+l43u33u14u24l32u24+u34l42u24+l43u34+u44⎦⎥⎥⎤

第二轮外循环后，下三角矩阵第二列被解出，上三角矩阵第三行被解出。数组 U 中存储的矩阵元素为：
[u11u12u13u14l21u22u23u24l31l32u33u34l41l42l43u33l43u34+u44]\begin{bmatrix} u_{11} & u_{12} & u_{13} & u_{14}\\ l_{21} & u_{22} & u_{23} & u_{24}\\ l_{31} & l_{32} & u_{33} & u_{34}\\ l_{41} & l_{42} &l_{43}u_{33} & l_{43}u_{34}+u_{44} \end{bmatrix}⎣⎢⎢⎡u11l21l31l41u12u22l32l42u13u23u33l43u33u14u24u34l43u34+u44⎦⎥⎥⎤

第三轮外循环…

直到解出 L,U\mathbf L, \mathbf UL,U 矩阵。

LUP矩阵分解

LUP分解是在LU分解的基础上增加主元的选取。

相对于LU分解，LUP分解增加了置换矩阵 P\mathbf PP ，有如下关系:
PA=LU(2)\mathbf P \mathbf A = \mathbf L \mathbf U \tag{2}PA=LU(2)

算法如下：

图2. LUP算法

如图2 所示，3~6行增加了对矩阵 L,U,P\mathbf L, \mathbf U, \mathbf PL,U,P 的换行，7~12行算法与LU分解相同。

为什么增加了置换矩阵 P\mathbf PP？其实我也是糊涂的（嘿嘿…

大概有以下几个原因：
1.被分解的矩阵 A\mathbf AA 其实被要求是非奇异矩阵。再有，如果主元直接为0，LU算法的在计算过程中会出现分母为0的情况。
2.如果主元不为 0 但是接近 0 呢？为了避免除数过小引起精度问题，所以会尽量选择大的主元。

经过行变换后的矩阵再经过LU分解就能尽量避免以上问题。

利用LUP分解对矩阵求逆

由(2)式，有：

P−1A−1=L−1U−1\mathbf P^{-1} \mathbf A^{-1} = \mathbf L^{-1} \mathbf U^{-1}P−1A−1=L−1U−1

L,U\mathbf L, \mathbf UL,U 为三角矩阵，P\mathbf PP 为置换矩阵，到这里其实就已经容易求出逆矩阵了。