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

背景：

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

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

方法：

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

用途

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

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

例题

（1）对矩阵 A = [ − 5 3 4 10 − 8 − 9 15 1 2 ] A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix} A=⎣⎡​−51015​3−81​4−92​⎦⎤​进行LU分解。
解：
A = [ − 5 3 4 10 − 8 − 9 15 1 2 ] A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix} A=⎣⎡​−51015​3−81​4−92​⎦⎤​有3行，所以， L L L矩阵应该为 3 × 3 3\times 3 3×3的矩阵。

先求 U U U矩阵：

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

A = [ − 5 3 4 10 − 8 − 9 15 1 2 ] → [ − 5 3 4 0 − 2 − 1 0 10 14 ] → [ − 5 3 4 0 − 2 − 1 0 0 9 ] = U A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}\rightarrow \begin{bmatrix}-5&3&4\\0&-2&-1\\0&10&14\end{bmatrix}\rightarrow \begin{bmatrix}-5&3&4\\0&-2&-1\\0&0&9\end{bmatrix}=U A=⎣⎡​−51015​3−81​4−92​⎦⎤​→⎣⎡​−500​3−210​4−114​⎦⎤​→⎣⎡​−500​3−20​4−19​⎦⎤​=U

再求 L L L矩阵：

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

[ − 5 0 0 10 − 2 0 15 10 9 ] → [ 1 0 0 − 2 1 0 − 3 − 5 1 ] = L \begin{bmatrix}-5&0&0\\10&-2&0\\15&10&9\end{bmatrix}\rightarrow \begin{bmatrix}1&0&0\\-2&1&0\\-3&-5&1\end{bmatrix}=L ⎣⎡​−51015​0−210​009​⎦⎤​→⎣⎡​1−2−3​01−5​001​⎦⎤​=L

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

验算：
L U = [ 1 0 0 − 2 1 0 − 3 − 5 1 ] [ − 5 3 4 0 − 2 − 1 0 0 9 ] = [ − 5 3 4 10 − 8 − 9 15 1 2 ] = A LU=\begin{bmatrix}1&0&0\\-2&1&0\\-3&-5&1\end{bmatrix}\begin{bmatrix}-5&3&4\\0&-2&-1\\0&0&9\end{bmatrix}=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}=A LU=⎣⎡​1−2−3​01−5​001​⎦⎤​⎣⎡​−500​3−20​4−19​⎦⎤​=⎣⎡​−51015​3−81​4−92​⎦⎤​=A

所以上述分解是正确的。

（2）对矩阵 A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix} A=⎣⎢⎢⎡​2−42−6​4−5−50​−13−47​5−81−3​−2181​⎦⎥⎥⎤​进行LU分解。

解：
A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix} A=⎣⎢⎢⎡​2−42−6​4−5−50​−13−47​5−81−3​−2181​⎦⎥⎥⎤​有4行，所以， L L L矩阵应该为 4 × 4 4\times 4 4×4的矩阵， L L L矩阵的第一列应该是矩阵 A A A的第一列除以其主元元素（2）的结果：

L = [ 1 0 0 0 − 2 1 0 0 1   1 0 − 3     1 ] L=\begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&\space&1&0\\-3&\space&\space&1\end{bmatrix} L=⎣⎢⎢⎡​1−21−3​01  ​001 ​0001​⎦⎥⎥⎤​

先求 U U U矩阵：

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

A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 − 9 − 3 − 4 10 0 12 4 12 − 5 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 4 7 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 0 5 ] = U A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&-9&-3&-4&10\\0&12&4&12&-5\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&4&7\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&0&5\end{bmatrix}=U A=⎣⎢⎢⎡​2−42−6​4−5−50​−13−47​5−81−3​−2181​⎦⎥⎥⎤​→⎣⎢⎢⎡​2000​43−912​−11−34​52−412​−2−310−5​⎦⎥⎥⎤​→⎣⎢⎢⎡​2000​4300​−1100​5224​−2−317​⎦⎥⎥⎤​→⎣⎢⎢⎡​2000​4300​−1100​5220​−2−315​⎦⎥⎥⎤​=U

再求 L L L矩阵：

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

[ 2 0 0 0 − 4 3 0 0 2 − 9 2 0 − 6 12 4 5 ] → [ 1 0 0 0 − 2 1 0 0 1 − 3 1 0 − 3 4 2 1 ] = L \begin{bmatrix}2&0&0&0\\-4&3&0&0\\2&-9&2&0\\-6&12&4&5\end{bmatrix}\rightarrow \begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&-3&1&0\\-3&4&2&1\end{bmatrix}=L ⎣⎢⎢⎡​2−42−6​03−912​0024​0005​⎦⎥⎥⎤​→⎣⎢⎢⎡​1−21−3​01−34​0012​0001​⎦⎥⎥⎤​=L

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

验算：
L U = [ 1 0 0 0 − 2 1 0 0 1 − 3 1 0 − 3 4 2 1 ] [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 0 5 ] = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] = A LU=\begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&-3&1&0\\-3&4&2&1\end{bmatrix}\begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&0&5\end{bmatrix}=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix}=A LU=⎣⎢⎢⎡​1−21−3​01−34​0012​0001​⎦⎥⎥⎤​⎣⎢⎢⎡​2000​4300​−1100​5220​−2−315​⎦⎥⎥⎤​=⎣⎢⎢⎡​2−42−6​4−5−50​−13−47​5−81−3​−2181​⎦⎥⎥⎤​=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 = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡​141​342​083​⎦⎤​，手动LU分解为：

L U = [ 1 0 0 4 1 0 1 1 8 1 ] [ 1 3 0 0 − 8 8 0 0 2 ] = [ 1 3 0 4 4 8 1 2 3 ] = F LU=\begin{bmatrix}1&0&0\\4&1&0\\1&\frac{1}{8}&1\end{bmatrix}\begin{bmatrix}1&3&0\\0&-8&8\\0&0&2\end{bmatrix}=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=F LU=⎣⎡​141​0181​​001​⎦⎤​⎣⎡​100​3−80​082​⎦⎤​=⎣⎡​141​342​083​⎦⎤​=F

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

L 1 U 1 = [ 1 0 0 1 4 1 0 1 4 1 2 1 ] [ 4 4 8 0 2 − 2 0 0 2 ] = [ 4 4 8 1 3 0 1 2 3 ] ≠ F L_1U_1=\begin{bmatrix}1&0&0\\\frac{1}{4}&1&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}\begin{bmatrix}4&4&8\\0&2&-2\\0&0&2\end{bmatrix}=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix}\neq F L1​U1​=⎣⎡​141​41​​0121​​001​⎦⎤​⎣⎡​400​420​8−22​⎦⎤​=⎣⎡​411​432​803​⎦⎤​​=F

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

L 2 U 2 = [ 1 4 1 0 1 0 0 1 4 1 2 1 ] [ 4 4 8 0 2 − 2 0 0 2 ] = [ 1 3 0 4 4 8 1 2 3 ] = F L_2U_2=\begin{bmatrix}\frac{1}{4}&1&0\\1&0&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}\begin{bmatrix}4&4&8\\0&2&-2\\0&0&2\end{bmatrix}=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=F L2​U2​=⎣⎡​41​141​​1021​​001​⎦⎤​⎣⎡​400​420​8−22​⎦⎤​=⎣⎡​141​342​083​⎦⎤​=F

  可见输出两个矩阵的方法得到LU矩阵是正确的，第一种方法得到的 L 1 U 1 ≠ F L_1U_1\neq F L1​U1​​=F。两种方法得到的LU矩阵都和手动算的结果不一样。

原因：

  第一种方法：输入[L1 U1 P] = lu(F)命令的时候，MATLAB首先生成 P F = [ 4 4 8 1 3 0 1 2 3 ] PF=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix} PF=⎣⎡​411​432​803​⎦⎤​，其中 P P P为置换矩阵（Permutation Matrix）， P = [ 0 1 0 1 0 0 0 0 1 ] P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix} P=⎣⎡​010​100​001​⎦⎤​，所以，这种方法求得的分解是 L 1 U 1 = P F L_1U_1=PF L1​U1​=PF而不是 L 1 U 1 = F L_1U_1=F L1​U1​=F。

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

L 2 = [ 1 4 1 0 1 0 0 1 4 1 2 1 ] = P T L 1 = [ 0 1 0 1 0 0 0 0 1 ] [ 1 0 0 1 4 1 0 1 4 1 2 1 ] L_2=\begin{bmatrix}\frac{1}{4}&1&0\\1&0&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}=P^TL_1=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}1&0&0\\\frac{1}{4}&1&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix} L2​=⎣⎡​41​141​​1021​​001​⎦⎤​=PTL1​=⎣⎡​010​100​001​⎦⎤​⎣⎡​141​41​​0121​​001​⎦⎤​

第二种方法算出来的 L 2 U 2 = F L_2U_2=F L2​U2​=F是成立的。

关于置换矩阵（Permutation Matrix）

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

  例如： P 21 = [ 0 1 0 1 0 0 0 0 1 ] P_{21}=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix} P21​=⎣⎡​010​100​001​⎦⎤​对于上述矩阵 F = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡​141​342​083​⎦⎤​，让 P 21 P_{21} P21​乘以 F F F就是把 F F F的第二行和第一行交换：
P 21 F = [ 0 1 0 1 0 0 0 0 1 ] [ 1 3 0 4 4 8 1 2 3 ] = [ 4 4 8 1 3 0 1 2 3 ] P_{21}F=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix} P21​F=⎣⎡​010​100​001​⎦⎤​⎣⎡​141​342​083​⎦⎤​=⎣⎡​411​432​803​⎦⎤​

  类似的，还有：

P 32 = [ 1 0 0 0 0 1 0 1 0 ] P_{32}=\begin{bmatrix}1&0&0\\0&0&1\\0&1&0\end{bmatrix} P32​=⎣⎡​100​001​010​⎦⎤​

P 31 = [ 0 0 1 0 1 0 1 0 0 ] P_{31}=\begin{bmatrix}0&0&1\\0&1&0\\1&0&0\end{bmatrix} P31​=⎣⎡​001​010​100​⎦⎤​

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

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

  例如，上述第一种方法求LU分解时，对于矩阵 F = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡​141​342​083​⎦⎤​，MATLAB首先生成 P F = [ 4 4 8 1 3 0 1 2 3 ] PF=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix} PF=⎣⎡​411​432​803​⎦⎤​，把 F F F的第一第二行进行交换，让第一列最大的元素4放在主元位置。

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

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