LU分解（LU Factorization）计算方法（手算+MATLAB），关于置换矩阵（Permutation Matrix），部分主元消去法（Partial Pivoting）

背景：

  求解一些列具有相同系数矩阵的线性方程，如：$Ax=b_1$，$Ax=b_2$，…$Ax=b_p$等。当矩阵$A$可逆时，可以先求出矩阵$A$的逆$A^{-1}$，再计算$A^{-1}{b}_1$，$A^{-1}{b}_2$等。

  但是，也可以用LU分解法来解这一系列方程：先使用初等行变换化简解出$Ax=b_1$，并同时得到矩阵$A$的LU分解，剩下的方程使用LU分解法求解即可。

方法：

  设矩阵$A$是$m\times n$的矩阵，它可以使用行变换（不进行行的交换）化简为阶梯阵，则矩阵$A$可以写成$A=LU$的形式。$L$是$m\times m$的单位下三角矩阵（即主对角线元素全为1的下三角矩阵）；$U$是$m\times n$的上三角阶梯形矩阵。

用途

  当$A=LU$时，方程$Ax=b$可以写成$L(Ux)=b$，令$Ux=y$，则有：
$$Ly=bUx=y$$

  即先求解$Ly=b$得到$y$，再求解$Ux=y$得到$x$，因为矩阵$L$和$U$都是三角矩阵，所以求解上述两个方程比直接求解$Ax=b$要简单。

例题

（1）对矩阵$A=\left[\begin{array}{ccc}-5& 3& 4\\ 10& -8& -9\\ 15& 1& 2\end{array}\right]$进行LU分解。
解：
$A=\left[\begin{array}{ccc}-5& 3& 4\\ 10& -8& -9\\ 15& 1& 2\end{array}\right]$有3行，所以，$L$矩阵应该为$3\times 3$的矩阵。

先求$U$矩阵：

先把矩阵$A$一步步消成上三角矩阵就能得到$U$矩阵：

$A=\left[\begin{array}{ccc}-5& 3& 4\\ 10& -8& -9\\ 15& 1& 2\end{array}\right]\to \left[\begin{array}{ccc}-5& 3& 4\\ 0& -2& -1\\ 0& 10& 14\end{array}\right]\to \left[\begin{array}{ccc}-5& 3& 4\\ 0& -2& -1\\ 0& 0& 9\end{array}\right]=U$

再求$L$矩阵：

把上面消元过程中的主元列的主元及其下面的元素都放入一个矩阵，然后每列除以此列的主元即得到$L$矩阵：

$\left[\begin{array}{ccc}-5& 0& 0\\ 10& -2& 0\\ 15& 10& 9\end{array}\right]\to \left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ -3& -5& 1\end{array}\right]=L$

(-5, 10, 15这一列每个元素除以主元-5；-2,10这一列每个元素除以主元-2；9这一列每个元素除以主元9，即得到$L$矩阵)

验算：
$LU=\left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ -3& -5& 1\end{array}\right]\left[\begin{array}{ccc}-5& 3& 4\\ 0& -2& -1\\ 0& 0& 9\end{array}\right]=\left[\begin{array}{ccc}-5& 3& 4\\ 10& -8& -9\\ 15& 1& 2\end{array}\right]=A$

所以上述分解是正确的。

（2）对矩阵$A=\left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ -4& -5& 3& -8& 1\\ 2& -5& -4& 1& 8\\ -6& 0& 7& -3& 1\end{array}\right]$进行LU分解。

解：
$A=\left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ -4& -5& 3& -8& 1\\ 2& -5& -4& 1& 8\\ -6& 0& 7& -3& 1\end{array}\right]$有4行，所以，$L$矩阵应该为$4\times 4$的矩阵，$L$矩阵的第一列应该是矩阵$A$的第一列除以其主元元素（2）的结果：

$L=\left[\begin{array}{cccc}1& 0& 0& 0\\ -2& 1& 0& 0\\ 1& & 1& 0\\ -3& & & 1\end{array}\right]$

先求$U$矩阵：

先把矩阵$A$一步步消成上三角矩阵就能得到$U$矩阵：

$A=\left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ -4& -5& 3& -8& 1\\ 2& -5& -4& 1& 8\\ -6& 0& 7& -3& 1\end{array}\right]\to \left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ 0& 3& 1& 2& -3\\ 0& -9& -3& -4& 10\\ 0& 12& 4& 12& -5\end{array}\right]\to \left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ 0& 3& 1& 2& -3\\ 0& 0& 0& 2& 1\\ 0& 0& 0& 4& 7\end{array}\right]\to \left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ 0& 3& 1& 2& -3\\ 0& 0& 0& 2& 1\\ 0& 0& 0& 0& 5\end{array}\right]=U$

再求$L$矩阵：

把上面消元过程中的主元列的主元及其下面的元素都放入一个矩阵，然后每列除以此列的主元就得到$L$矩阵：

$\left[\begin{array}{cccc}2& 0& 0& 0\\ -4& 3& 0& 0\\ 2& -9& 2& 0\\ -6& 12& 4& 5\end{array}\right]\to \left[\begin{array}{cccc}1& 0& 0& 0\\ -2& 1& 0& 0\\ 1& -3& 1& 0\\ -3& 4& 2& 1\end{array}\right]=L$

(2，-,4, 2, -6这一列每个元素除以主元2；3，-9,12这一列每个元素除以主元3；2,4这一列每个元素除以主元2，5这一列除以主元5，即得到$L$矩阵)

验算：
$LU=\left[\begin{array}{cccc}1& 0& 0& 0\\ -2& 1& 0& 0\\ 1& -3& 1& 0\\ -3& 4& 2& 1\end{array}\right]\left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ 0& 3& 1& 2& -3\\ 0& 0& 0& 2& 1\\ 0& 0& 0& 0& 5\end{array}\right]=\left[\begin{array}{ccccc}2& 4& -1& 5& -2\\ -4& -5& 3& -8& 1\\ 2& -5& -4& 1& 8\\ -6& 0& 7& -3& 1\end{array}\right]=A$

所以上述分解是正确的。

关于MATLAB中的LU分解函数

  MATLAB中提供了LU分解的函数lu( )，有两种用法，一种是返回三个矩阵（L,U和置换矩阵P），一种是只返回L和U矩阵。例如：

F =
     1     3     0
     4     4     8
     1     2     3

>> [L1 U1 P] = lu(F)  //使用MATLAB内置的lu( )函数求矩阵F的LU分解，返回三个矩阵L,U和P的用法
L1 =
    1.0000         0         0
    0.2500    1.0000         0
    0.2500    0.5000    1.0000

U1 =
     4     4     8
     0     2    -2
     0     0     2
P =
     0     1     0
     1     0     0
     0     0     1

>> [L2 U2] = lu(F) //使用MATLAB内置的lu( )函数求矩阵F的LU分解，返回2个矩阵L和U的用法
L2 =
    0.2500    1.0000         0
    1.0000         0         0
    0.2500    0.5000    1.0000
U2 =
     4     4     8
     0     2    -2
     0     0     2

对于上述矩阵$F=\left[\begin{array}{ccc}1& 3& 0\\ 4& 4& 8\\ 1& 2& 3\end{array}\right]$，手动LU分解为：

$LU=\left[\begin{array}{ccc}1& 0& 0\\ 4& 1& 0\\ 1& \frac{1}{8}& 1\end{array}\right]\left[\begin{array}{ccc}1& 3& 0\\ 0& -8& 8\\ 0& 0& 2\end{array}\right]=\left[\begin{array}{ccc}1& 3& 0\\ 4& 4& 8\\ 1& 2& 3\end{array}\right]=F$

上面MATLAB第一种方法(L1，U1和P三个矩阵)算出来的结果：

$L_1{U}_1=\left[\begin{array}{ccc}1& 0& 0\\ \frac{1}{4}& 1& 0\\ \frac{1}{4}& \frac{1}{2}& 1\end{array}\right]\left[\begin{array}{ccc}4& 4& 8\\ 0& 2& -2\\ 0& 0& 2\end{array}\right]=\left[\begin{array}{ccc}4& 4& 8\\ 1& 3& 0\\ 1& 2& 3\end{array}\right]\ne F$

MATLAB第二种方法(L2和U2 两个矩阵)算出来的结果：

$L_2{U}_2=\left[\begin{array}{ccc}\frac{1}{4}& 1& 0\\ 1& 0& 0\\ \frac{1}{4}& \frac{1}{2}& 1\end{array}\right]\left[\begin{array}{ccc}4& 4& 8\\ 0& 2& -2\\ 0& 0& 2\end{array}\right]=\left[\begin{array}{ccc}1& 3& 0\\ 4& 4& 8\\ 1& 2& 3\end{array}\right]=F$

  可见输出两个矩阵的方法得到LU矩阵是正确的，第一种方法得到的$L_1{U}_1\ne F$。两种方法得到的LU矩阵都和手动算的结果不一样。

原因：

  第一种方法：输入[L1 U1 P] = lu(F)命令的时候，MATLAB首先生成$PF=\left[\begin{array}{ccc}4& 4& 8\\ 1& 3& 0\\ 1& 2& 3\end{array}\right]$，其中$P$为置换矩阵（Permutation Matrix），$P=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]$，所以，这种方法求得的分解是$L_1{U}_1=PF$而不是$L_1{U}_1=F$。

  第二种方法：输入[L2 U2] = lu(F)命令的时候，算出来的$U_2$矩阵和第一种方法一样（但是和手算的不一样）；算出来的$L_2$矩阵其实是$L_2=P^T{L}_1$，即：

$L_2=\left[\begin{array}{ccc}\frac{1}{4}& 1& 0\\ 1& 0& 0\\ \frac{1}{4}& \frac{1}{2}& 1\end{array}\right]=P^T{L}_1=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ \frac{1}{4}& 1& 0\\ \frac{1}{4}& \frac{1}{2}& 1\end{array}\right]$

第二种方法算出来的$L_2{U}_2=F$是成立的。

关于置换矩阵（Permutation Matrix）

  上MATLAB计算LU分解时，两种方法都用到了置换矩阵$P$。
  置换矩阵是(0,1)矩阵，即只有0和1两种元素的矩阵，其作用是让矩阵的行进交换的矩阵。置换矩阵$P$是改变相应阶数的单位矩阵得到的，例如对单位矩阵$I$交换第一第二行得到$P_{21}=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]$，那么得到的置换矩阵$P_{21}$的功能就是，使得和它相乘的矩阵也交换第一第二行。

  例如：$P_{21}=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]$对于上述矩阵$F=\left[\begin{array}{ccc}1& 3& 0\\ 4& 4& 8\\ 1& 2& 3\end{array}\right]$，让$P_{21}$乘以$F$就是把$F$的第二行和第一行交换：
$P_{21}F=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 3& 0\\ 4& 4& 8\\ 1& 2& 3\end{array}\right]=\left[\begin{array}{ccc}4& 4& 8\\ 1& 3& 0\\ 1& 2& 3\end{array}\right]$

  类似的，还有：

$P_{32}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right]$

$P_{31}=\left[\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right]$

关于部分主元消去法（Partial Pivoting）

  部分主元消去法在求解线性系统时常被用到（如LU分解和使用A\b求$Ax=b$等），在计算机的浮点数计算中能提高精确度。

  例如，上述第一种方法求LU分解时，对于矩阵$F=\left[\begin{array}{ccc}1& 3& 0\\ 4& 4& 8\\ 1& 2& 3\end{array}\right]$，MATLAB首先生成$PF=\left[\begin{array}{ccc}4& 4& 8\\ 1& 3& 0\\ 1& 2& 3\end{array}\right]$，把$F$的第一第二行进行交换，让第一列最大的元素4放在主元位置。

  类似的，使用部分主元消去法时，每一行主元位置的元素与其下面的元素中最大的元素放在主元位置（使用行交换将最大的元素换到主元位置），后面使用高斯消元法，这就是部分主元消去法。所以，上面求解中MATLAB使用部分主元消去法时，先使用置换矩阵$P$来调整主元使之最大，然后再求解。

  部分主元消去法就是多了一个把最大的值调到主元位置的操作，却能提高计算精度。