LUP分解方法，矩阵求逆

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

---

LU矩阵分解

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

该分解可以采用高斯消元法来进行计算。若想将给定矩阵 $\mathbf{A}$ 分解为下三角矩阵 $\mathbf{L}$ 和上三角矩阵 $\mathbf{U}$ ，一个思路就是通过一系列的初等变换将 (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 初始存储矩阵 $\mathbf{A}$， 经过计算后存储上三角矩阵 和 除主对角线以外的下三角矩阵 的元素。

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

第一行 $u_{11}:u_{14}$ 可以直接解出。例如在第一轮外循环后，下三角矩阵第一列被解出，上三角第二行被解出。数组 U 中 存储的矩阵元素为：
$$\left[\begin{array}{cccc}{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{array}\right]$$

第二轮外循环后，下三角矩阵第二列被解出，上三角矩阵第三行被解出。数组 U 中 存储的矩阵元素为：
$$\left[\begin{array}{cccc}{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{array}\right]$$

第三轮外循环…

直到解出 $\mathbf{L},\mathbf{U}$ 矩阵。

LUP矩阵分解

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

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

算法如下：

图2. LUP算法

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

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

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

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

利用LUP分解对矩阵求逆

由(2)式，有：

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

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