什么是LU分解
将一个矩阵 $A$ 表示成 $A = LU$ 的形式，其中 $L$ 是下三角矩阵， $U$ 是上三角矩阵。这就是所谓的LU分解。例如：
\[A = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 13 \end{pmatrix} \\ L = \begin{pmatrix} 1 & 0 & 0 \\ 4 & 1 & 0 \\ 7 & 2 & 1 \end{pmatrix} \quad U = \begin{pmatrix} 1 & 2 & 3 \\ 0 & -3 & -6 \\ 0 & 0 & 4 \end{pmatrix}\]上面的例子中 $A$ 是方阵。如果 $A$ 不是方阵呢？这种情况下，我们希望 $L$ 和 $U$ 中有一个是方阵，并且是三角矩阵；而另一个尽量做成三角矩阵。比如这样的形式：
\[B = \begin{pmatrix} 1 & 0 & 0 & 0 \\ * & 1 & 0 & 0 \\ * & * & 1 & 0 \\ * & * & * & 1 \\ * & * & * & * \\ * & * & * & * \\ \end{pmatrix} \begin{pmatrix} * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \\ 0 & 0 & 0 & * \end{pmatrix}\]为什么要做LU分解
LU分解其实就是一个运算优化。无论是求解行列式，还是线性方程组等，完全可以直接使用 $A$ 本身来求解。但使用 $LU$ 矩阵的计算量非常小，小到相对来说是可以忽略的。更何况实际编程中，我们还可以将 $LU$ 矩阵保存，即可以用来求行列式，也可以用来求解线性方程组，可谓“一次计算，到处方便”。下面就分别举例看一下。
利用LU分解求行列式
行列式有一个性质，即如果 $A=LU$ ，则 $detA = det(LU) = (detL)(detU)$ 。
行列式还有一个性质，即对于对角矩阵 $B$ ，其行列式的值对角元素的乘积。
因此对于文章开头的例子中的矩阵 $A$ ，其行列式 $detA = (1 * 1 * 1)(1 * (-3) * 4) = -12$
利用LU分解求解线性方程组
如果要求 $Ax = y$ 中的向量 $x$ ，将 $A$ 进行 $LU$ 分解以后，则 $LUx = y$ 。因此求解 $x$ 的问题变成了先求 $Lz = y$ 中的向量 $z$ ，再求 $Ux = z$ 中的向量 $x$ 。乍一看LU分解没有简化问题，反而使问题复杂化了。其实不然，因为 $L$ 和 $U$ 都是三角矩阵，而求解三角矩阵构成的方程组是非常简单的。仍然以文章开头的例子举例，要求
\[Ax = \begin{pmatrix} 1\\ 1 \\ 1 \end{pmatrix}\]的解 $x$ 。
因为 $A=LU$ ，令 $Ux=z$ ，所以 $Lz = y$ 。因此我们先求向量 $z$ ：
\[z = (z_1, z_2, z_3)^T \\ Lz = y \quad => \\ \begin{pmatrix} 1 & 0 & 0 \\ 4 & 1 & 0 \\ 7 & 2 & 1 \end{pmatrix} \begin{pmatrix} z_1\\z_2\\z_3\end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} \quad =>\]由此可以得到一个很简单的线性方程组：
\[z_1 = 1 \\ 4z_1 + z_2 = 1 \\ 7z_1 + 2z_2 + z_3 = 1\]所以 $z = (1, -3, 0)^T$
有了 $z$ 向量以后，接下来我们求解 $Ux = z$ ，即
\[\begin{pmatrix} 1 & 2 & 3 \\ 0 & -3 & -6 \\ 0 & 0 & 4 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 1 \\ -3 \\ 0 \end{pmatrix}\]我们也可以得到一个很简单的线性方程组：
\[x_1 + 2x_2 + 3x_3 = 1 \\ -3x_2 - 6x_3 = -3 \\ 4x_3 = 0\]由此可得到线性方程组的解 $x = (-1, 1, 0)^T$。
从上面的过程中可以看出来，进行LU分解以后，整个求解过程都非常简单。
使用python库进行LU分解
虽然这篇文章主要介绍LU分解，但我目前并不打算学习和介绍使用笔算进行LU分解。取而代之的是使用python库求解。
import sympy
A = sympy.Matrix([[2, 4, 2], [1, 5, 2], [4, -1, 9]])
L, U, P = A.LUdecomposition()
#output:
#Matrix([
#[ 1, 0, 0],
#[1/2, 1, 0],
#[ 2, -3, 1]])
#U
#Matrix([
#[2, 4, 2],
#[0, 3, 1],
#[0, 0, 8]])
#P
#[]
import scipy.linalg
A = scipy.matrix([[2, 4, 2], [1, 5, 2], [4, -1, 9]])
P, L, U = scipy.linalg.lu(A)
#output:
#P:
#array([[ 0., 0., 1.],
# [ 0., 1., 0.],
# [ 1., 0., 0.]])
#L:
#array([[ 1. , 0. , 0. ],
# [ 0.25 , 1. , 0. ],
# [ 0.5 , 0.85714286, 1. ]])
#U:
#array([[ 4. , -1. , 9. ],
# [ 0. , 5.25 , -0.25 ],
# [ 0. , 0. , -2.28571429]])
上面代码中，$P$ 称为置换矩阵。
总结
这篇文章介绍了什么是LU分解，更主要的是，为什么要进行LU分解：进行分解后，可以利用分解结果多次求解，从而增加程序效率。
参考资料 ：
- 《程序员的数学 - 线性代数》