大名鼎鼎的LU分解——Matlab解线性方程组（3）

前文已经完成了铺垫，这里我们开始进行矩阵的LU分解。第一节我们看到了，因为求逆矩阵会导致运算时间增长、运算量增大，所以在一般求解线性方程组时，常用高斯消元的方法。将高斯消元法凝练、规范化，给出矩阵的LU分解法，将系数矩阵分解成L（乘子矩阵）、U（消元后的矩阵）和P（排列矩阵）——自己起的名字，反正就这么个玩意。

植物之怒

7953人浏览 · 2022-01-21 20:26:21

植物之怒 · 2022-01-21 20:26:21 发布

前言

一、伟大的LU分解法

其实之前已经接触到了LU分解的内容，最内核的就是高斯消元法，这也是在大学线性代数课程中教给大家的方法。

总体上看来，高斯消元法包括了两个阶段：前向消去（forward elimination）和回代（back substitution）。

前向消去要进行n-1步，其中第k步的工作是：将第k个方程以下的方程分别减去第k个方程的某个倍数，从而消去第k个变量。

回代有n步，意思是当所有前向消去都进行完成后，第n个方程此时应该只剩下一个变量，这时候就可以通过一步除法得出该变量的值，进而反代回第n-1个方程。此时n-1个方程只剩下两个变量，其中一个变量的值由上一步给出，同理循环即可求出所有变量。

这里需要注意一点，在前向消去时，如果第k行第k个变量的系数（主元）很小，则需要将k以下的某个方程中，第k个变量系数最大的那个换上来。这种方法看上去有点多余，但是在计算机运算中是必须的，下一条会讲到。

1.前向消去的步骤

停！

这个时候，需要回顾一下高斯消元整体的过程，要包括到每个细节：

首先找到所有方程第一个变量系数最大的那个方程是哪个，然后把这个变量所在的方程交换到第一个来。这个过程用矩阵语言描述，就是：找到系数矩阵所有行第一列元素中最大的元素所在行，然后将这一行换到第一行来，形成一个新的系数矩阵。进行这一步，需要进行一次行变换，用矩阵语言表示，就是需要堆系数矩阵A左乘一个排列矩阵P1：

$A_{1}=P_{1}A$

然后，确认第一个方程第一个变量的系数为主元，依次消去剩下所有行的第一个变量的系数。具体做法为，先将主元乘以a1倍，然后与第二个方程第一个变量的系数加减抵消成0；再将主元乘以a2倍，与第三行第一个变量的系数加减抵消成0……，直到把第一个方程以下所有方程第一个变量的系数全部变成0. 用矩阵语言描述为：将第一行第一个元素确定为主元，然后把第一行的某个倍数加给第二行，用以消去第二行第一个元素，以此类推将第一行以下所有行的第一列元素变成0.

这里需要展开讲解：

用第一节那个例子。假设矩阵为3*3的方阵，第一行主元是10，第二行第三行的第一个元素分别为-3和5 。想用第一行的10消去第二行的-3，需要给第一行成0.3倍再加到第二行去。想用第一行的10消去第三行的5，需要给第一行成0.5倍再减到第三行去。用矩阵表示为：

$M_{1}=\begin{pmatrix} 1 & 0 &0 \\ 0.3 & 1& 0\\ -0.5 & 0& 1 \end{pmatrix}$

（第一行的0.3倍加在第二行，第一行的0.5倍减在第二行）

这是一个下三角矩阵。可以说，所有进行这种变换的M矩阵都是下三角矩阵。

那么此时，经过一番变换后，只剩下第一行的第一列有非零元素在，其他行第一列都是0.

也即：进行了两个矩阵变化后，第一行的第一列（或者说第一个变量）被保留了下来。

$A_{1}=M_{1}P_{1}A$

这个式子表示，在消第一列（或者说第一个变量）时，进行了一次排列矩阵左乘变换（行变换）和一次下三角矩阵左乘变换（行变换）。这便是前向消去的第一步。总共有n-1步。

这n-1步消元全部完成后，会得到一个崭新的系数矩阵U。

$U=M_{n-1}P_{n-1}......M_{2}P_{2}M_{1}P_{1}A$

值得注意的是，此时产生的U矩阵，由于第1行以下第1列元素全是0，……第k行以下第k列元素全是0……，第n-1行以下第n-1列元素全是0，所以U矩阵是一个妥妥的上三角矩阵。

此时，将所有的M变换矩阵全都挪到方程右边去，经过这个变换后的M矩阵全变成了L矩阵。

$L_{1}L_{2}......L_{n-1}U=P_{n-1}......P_{2}P_{1}A$

因为M矩阵和P矩阵都是对A进行行变换的矩阵，所以只要按顺序挪过去，不会影响最终结果。值得注意的是，M矩阵是倍加矩阵（通过标准矩阵进行倍加变换得到的矩阵，名字也是我编的），倍加矩阵移到等式另一边，可以直接给非对角线元素添上负号。所以移到等式左边的L矩阵依然是下三角矩阵。

不妨令：

$L=L_{1}L_{2}...L_{n-1}$

$P=P_{1}P_{2}...P_{n-1}$

方程可以简化为：

$LU=PA$

其中：乘子矩阵L包含了所有的乘子，是一个下三角矩阵；U矩阵是消元结束后剩下的系数矩阵，是一个上三角矩阵；P矩阵是初等变换经过排列得到的排列矩阵。

这个过程，称为矩阵A的LU分解（LU factorization）或者三角分解（triangular decomposition）。

2.LU分解有啥用？

表面上看，LU分解进行了大批量的矩阵初等变换，得到了一大堆没有什么用的东西。

但是，如果我们把这东西带回去，会发生什么事呢？

对于方程：

$Ax=b$

由LU分解可知：

$LU=PA$

$A=P^{-1}LU$

原方程变为：

$P^{-1}LUx=b$

$LUx=Pb$

不妨令Ux=y,则又有：

$\left\{\begin{matrix} Ly=Pb\\ Ux=y \end{matrix}\right.$

其中，Pb看作对于已知结果构成的列向量进行行交换变换。通俗点就是把方程结果放一个列向量，起名叫b，对这个b的行与行按照P矩阵规定的那样去交换，得到一个交换过的列向量Pb。

通过上述式子，已知系数矩阵L和结果向量Pb，可以求出未知向量y。再把未知向量y当成已知向量，结合第二个系数矩阵U求出未知数向量x。

所以这有什么用呢？

回想一下上一节的内容——虽然对于一般系数矩阵求解方程有难度，但是对于三角形矩阵当作系数矩阵的线性方程组，有非常简洁快速的算法去求其解。

而求y时候的系数矩阵L和求x时候的系数矩阵U都是三角形矩阵。

这样，就把求一般系数矩阵下的方程，化简为——①进行LU分解，②化成两个三角形矩阵的求解。

二、多出的一步

1.为什么要选主元？

在LU分解中，有一个重要的步骤被我提到——那就是选主元。用人话说，就是在轮到第k行消元的时候，要选择第k行到第n行所有的第k个元素中最大的那个元素作为主元，并把他换到第k行来。

but why？这一步显得有点多余。正常看来，不论第k行的第k个元素是大是小，经过矩阵基本初等变换后求解出的U矩阵都是唯一的，解出来的x也是唯一的。

但是不要忘了，我们目前讨论的这些有的没的，都是在思考计算机的算法。包括未来对于大批量数据处理，得到的数据量绝对不是人类可以用纸笔能处理得了的。而用到计算机算法，就要考虑到编程软件自身的精确度。

考虑上两节讲的那个例子，把这个例子中的系数矩阵简单变一下：

$\begin{pmatrix} 10 & -7 &0 \\ -3 &2.099 &6 \\ 5 & -1& 5 \end{pmatrix}\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3} \end{pmatrix}=\begin{pmatrix} 7\\ 3.901\\ 6 \end{pmatrix}$

易知，这个解仍然是

$\begin{pmatrix} 0\\ -1\\ 1 \end{pmatrix}$

但是如果我们计算用的电脑计算性能不佳，只能进行5位有效数字的浮点数运算，那么会怎样呢？

$\begin{pmatrix} 10 & -7 &0 \\ 0 &-0.001 &6 \\ 0 &2.5& 5 \end{pmatrix}\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3} \end{pmatrix}=\begin{pmatrix} 7\\ 6.001\\ 2.5 \end{pmatrix}$