2.6 A 的 LU 分解

一、A = LU

线性代数很多关键的概念实际上就是矩阵的分解（factorization）—— 原始矩阵 $A$ 变成两个或三个特殊矩阵的乘积。这节是第一个分解，实际上也是最重要的分解。$A=LU$ 的分解来自于消元法，其中因子 $L$ 和 $U$ 都是三角形矩阵
矩阵 $U$ 是上三角矩阵，它的主元都在对角线上，消元步骤将 $A$ 变为 $U$。现在我们要反向这些步骤（将 $U$ 变为 $A$ ），通过一个下三角矩阵 $L$ 就可以。$L$ 的元素正好就是乘数 $l_{ij}$ —— 从行 $i$ 减去乘数乘主元行 $j$。
以 $2\times 2$ 的矩阵为例，矩阵 $A$ 有四个元素 $2,1,6,8$，要消去的元素是 $6$。从行 $2$ 减去 $3$ 乘行 $1$，这个正向步骤使用消元矩阵 $E_{21}$，乘数 $l_{21}=3$。从 $U$ 到 $A$ 的反向步骤使用 $L=E_{21}^{-1}$（行 $2$ 加上 $3$ 乘行 $1$）。$$\begin{gathered} 正向从A到U：E_{21}A=\left[\begin{array}{cc}1& 0\\ -3& 1\end{array}\right]\left[\begin{array}{cc}2& 1\\ 6& 8\end{array}\right]=\left[\begin{array}{cc}2& 1\\ 0& 5\end{array}\right]=U \\ 反向从U到A：E_{21}^{-1}U=\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right]\left[\begin{array}{cc}2& 1\\ 0& 5\end{array}\right]=\left[\begin{array}{cc}2& 1\\ 6& 8\end{array}\right]=A \end{gathered}$$上面第二行就是分解 $LU=A$，将 $E_{21}^{-1}$ 用 $L$ 代替。更大的矩阵会有很多 $E^{\prime }s$，$L$ 包含它们所有的逆矩阵。
$A$ 到 $U$ 的每一个步骤都要左乘一个矩阵 $E_{ij}$，将位置 $(i,j)$ 的元素变为 $0$。为了清晰起见，假设没有行交换，如果 $A$ 是 $3\times 3$ 的矩阵，需要在左边依次乘 $E_{21}$，$E_{31}$，$E_{32}$，乘数 $l_{ij}$ 将会使得位置 $(2,1)$，$(3,1)$，$(3,2)$ 位置处的元素都变为 $0$，它们都在对角线下方，消元法在得到一个上三角矩阵后结束。
现在将 $E^{\prime }s$ 移到另外一边，将它们的逆矩阵乘上 $U$：

$$(E_{32}{E}_{31}{E}_{21})A=U变为A=(E_{21}^{-1}{E}_{31}^{-1}{E}_{32}^{-1})U就是A=LU(\mathrm{2.6.1})$$

逆矩阵是反序相乘，三个逆矩阵的乘积就是 $L$。我们得到了 $A=LU$。

二、解释与例子

第一点：每一个逆矩阵 $E^{-1}$ 都是下三角矩阵。它的非主对角线元素是 $l_{ij}$，用来恢复 $-l_{ij}$ 产生的减法。$E$ 和 $E^{-1}$ 的主对角线都是 $1$。上面的例子中 $l_{21}=3$，$E=\left[\begin{array}{cc}1& 0\\ -3& 1\end{array}\right]$，$L=E^{-1}=\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right]$。
第二点：式（2.6.1）展示了一个下三角矩阵（$E_{ij}$ 的乘积）乘 $A$，也展示了所有的 $E_{ij}^{-1}$ 乘 $U$ 会得到 $A$。$E_{ij}$ 的逆矩阵的乘积得到的下三角矩阵就是 $L$。
我们处理这些逆矩阵的一个原因是想要分解 $A$，而不是 $U$。它的 “反向形式” 得到了 $A=LU$。另一个原因是我们会得到更多的信息，$L$ 是一个很好的矩阵，这也是第三点。
第三点：每个乘数 $l_{ij}$ 可以直接放入 $L$ 的 $(i,j)$ 位置，不需要改变。通常矩阵的乘法会将这些位置弄乱，但是在 $L$ 里不会。因为逆矩阵的正确顺序，使得 $l$ 没有发生变化。在式（2.6.3）中会给出原因。
第四点：每个 $E^{-1}$ 的对角线都是 $1$，$L$ 也是如此。

$$\begin{array}{c}A=LU\end{array}\begin{array}{c}消元过程中没有行交换.上三角矩阵U的\\ 主元在它的对角线上.下三角矩阵L的主\\ 元都是1，乘数l_{ij}在L的对角线下方.\end{array}$$

【例1】消元法从行 $2$ 减去 $\frac{1}{2}$ 乘行 $1$，最后一步从行 $3$ 减去 $\frac{2}{3}$ 乘行 $2$。下三角矩阵 $L$ 有 $l_{21}=\frac{1}{2}$，$l_{32}=\frac{2}{3}$。$LU$ 的乘积得到 $A$：$$A=\left[\begin{array}{ccc}2& 1& 0\\ 1& 2& 1\\ 0& 1& 2\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ \frac{1}{2}& 1& 0\\ 0& \frac{2}{3}& 1\end{array}\right]\left[\begin{array}{ccc}2& 1& 0\\ 0& \frac{3}{2}& 1\\ 0& 0& \frac{4}{3}\end{array}\right]=LU$$因为 $A$ 的 $(3,1)$ 元素是 $0$，所以 $(3,1)$ 位置的乘数是 $0$，即不需要进行操作。

【例2】将 $A$ 左上角的元素 $2$ 改为 $1$，变成 $B$。则所有的主元都是 $1$，所有的乘数也是 $1$。保持这种模式，当 $B$ 是 $4\times 4$ 的矩阵时：$$B=\left[\begin{array}{cccc}1& 1& 0& 0\\ 1& 2& 1& 0\\ 0& 1& 2& 1\\ 0& 0& 1& 2\end{array}\right]=\left[\begin{array}{cccc}1& & & \\ 1& 1& & \\ 0& 1& 1& \\ 0& 0& 1& 1\end{array}\right]\left[\begin{array}{cccc}1& 1& 0& 0\\ & 1& 1& 0\\ & & 1& 1\\ & & & 1\end{array}\right]$$假设没有行交换，如何可以得知 $L$ 和 $U$ 中哪些元素为 $0$ 呢？$$\begin{gathered} 当A的某一行从0开始，则L的该行也从0开始 \\ 当A的某一列从0开始，则U的该列也从0开始 \end{gathered}$$如果某一行从 $0$ 开始，那么就不需要消元，$L$ 相应的位置就是 $0$，这将节省电脑的时间。同样的，如果某一列从 $0$ 开始，$U$ 相应的位置也为 $0$。但是，如果 $0$ 在中间，因为消元法是前向消除，这些位置在 $L$ 或 $U$ 的对应的位置大概率就不再是 $0$。那么为什么 $L$ 相应的位置是乘数 $l_{ij}$，而不发生混乱呢？

关键原因是 $A$ 为什么等于 $LU$：当主元行下方的行减去时，它们还是 $A$ 的原始行吗？不是！因为在消元过程中它们被改变了。那么它们是 $U$ 的行吗？是的！因为主元不再改变。当计算 $U$ 的第三行时，我们会减去乘数乘 $U$ 前面的行（不是 $A$ 的行）：$$Row3ofU=(Row3ofA)-l_{31}(Row1ofU)-l_{32}(Row2ofU)(\mathrm{2.6.2})$$改写这个方程，看看行 $\left[\begin{array}{ccc}{l}_{31}& l_{32}& 1\end{array}\right]$ 是如何与 $U$ 相乘的：$$(Row3ofA)=l_{31}(Row1ofU)+l_{32}(Row2ofU)+1(Row3ofU)(\mathrm{2.6.3})$$正好就是 $(A=LU)$ 的第三行。$L$ 的行 $3$ 的分量是 $l_{31}$，$l_{32}$，$1$。无论 $A$ 有多大，所有的行都是这样。如果没有行交换，则有 $A=LU$。

平衡形式 $LDU$：$A=LU$ 是不对称的，因为 $U$ 的对角线上是主元，而 $L$ 的对角线都是 $1$。将 $U$ 分成一个举着 $D$ 和一个新的矩阵 $U$ 相乘，矩阵 $D$ 是对角线矩阵，它的对角线就是 $U$ 的主元，而新的矩阵 $U$ 的对角线就会变成 $1$，其它元素除以该行的主元：

新的上三角矩阵 $U$ 的对角线都是 $1$。与正常的 $LU$ 不同的是，新形式中间有一个 $D$：下三角矩阵 $L$ 乘对角线矩阵 $D$ 乘上三角矩阵 $U$。

$$矩阵的三角分解可以写A=LU或A=LDU$$

当看到 $LDU$ 的形式时，我们就可以知道 $U$ 的对角线都是 $1$，每一行都除以其第一个非零元素 —— 主元（ $LU$中 $U$ 的主元 ）：$$\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right]\left[\begin{array}{cc}2& 8\\ 0& 5\end{array}\right]进一步分解为\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right]\left[\begin{array}{cc}2& \\ & 5\end{array}\right]\left[\begin{array}{cc}1& 4\\ 0& 1\end{array}\right](\mathrm{2.6.4})$$主元 $2$ 和 $5$ 进入了 $D$，行分别除以 $2$ 和 $5$，得到新的 $U$，它的对角线都是 $1$。乘数 $3$ 仍然在 $L$ 内。

三、一个方形系统 = 两个三角形系统

矩阵 $L$ 包含了高斯消元法的记忆，它保存了每次进行消元时的乘数。我们可以使用这些求解 $A\mathbf{x}=\mathbf{b}$。
当存在右侧的 $\mathbf{b}$ 时，则需要 $L$。因子 $L$ 和 $U$ 完全取决于左侧（矩阵 $A$）。在 $A\mathbf{x}=\mathbf{b}$ 的右侧，我们先使用 $L^{-1}$ 再使用 $U^{-1}$。该求解步骤会处理两个三角形矩阵。

1、分解（通过对左侧的矩阵 $A$ 进行消元，得到 $L$ 和 $U$）
2、求解（使用 $L$ 对 $\mathbf{b}$ 进行前向消元，然后使用 $U$ 进行回代求得 $\mathbf{x}$）

以前，我们使用增广矩阵 $\left[\begin{array}{cc}A& \mathbf{b}\end{array}\right]$ 同时处理 $A$ 和 $\mathbf{b}$。但是电脑大多数会将两侧分开，$L$ 和 $U$ 都保存有消元的记忆，无论何时我们都可以处理 $\mathbf{b}$。因为这样求解一个单一的系统只需要一个子程序。
那么如何使用 $\mathbf{b}$ 呢？首先对右侧使用前向消元法（乘数存储在 $L$ 中，现在可以使用），它会将 $\mathbf{b}$ 变成一个新的右侧 $\mathbf{c}$，我们现在求解的是 $L\mathbf{c}=\mathbf{b}$，然后使用回代求解 $U\mathbf{x}=\mathbf{c}$，原始系统 $A\mathbf{x}=\mathbf{b}$ 就被分解成了两个三角系统：

$$前向与反向求解L\mathbf{c}=\mathbf{b}，然后求解U\mathbf{x}=\mathbf{c}(\mathrm{2.6.5})$$

要验证 $\mathbf{x}$ 就是要求的解，$U\mathbf{x}=\mathbf{c}$ 两侧同时左乘 $L$，得到 $LU\mathbf{x}=L\mathbf{c}$ 就是 $A\mathbf{x}=\mathbf{b}$。
强调： 这些步骤并没有新的知识。我们使用前向消元求解三角系统 $L\mathbf{c}=\mathbf{b}$，然后回代求解 $U\mathbf{x}=\mathbf{c}$。

【例3】以 $A\mathbf{x}=\mathbf{b}$ 前向消元开始，在 $U\mathbf{x}=\mathbf{c}$ 结束：$$A\mathbf{x}=\mathbf{b}\begin{array}{c}u+2v=5\\ 4u+9v=21\end{array}变为\begin{array}{c}u+2v=5\\ v=1\end{array}U\mathbf{x}=\mathbf{c}$$乘数 $4$ 保存在 $L$ 中，右侧使用 $4$ 将 $21$ 变成了 $1$：$$L\mathbf{c}=\mathbf{b}下三角系统\left[\begin{array}{cc}1& 0\\ 4& 1\end{array}\right][\mathbf{c}]=\left[\begin{array}{c}5\\ 21\end{array}\right]求得\mathbf{c}=\left[\begin{array}{c}5\\ 1\end{array}\right]$$$$U\mathbf{x}=\mathbf{c}上三角系统\left[\begin{array}{cc}1& 2\\ 0& 1\end{array}\right][\mathbf{x}]=\left[\begin{array}{c}5\\ 1\end{array}\right]求得\mathbf{x}=\left[\begin{array}{c}3\\ 1\end{array}\right]$$$L$ 和 $U$ 所使用的也就是以前 $A$ 所使用的 $n^2$ 的存储空间。

四、消元法的成本

这里讨论消元法的成本 —— 即计算时间的问题。我们在计算机上解方程，就需要考虑计算成本，在科学计算时我们可能会遇到大型系统，三维空间的问题就很容易有一百万个未知数，如果计算成本太高的话，我们不可能让计算机计算成百上千年。
消元法的第一阶段是将列 $1$ 的第一主元以下的元素全部变为 $0$，第一行以下的元素全部都需要改变，而改变一个元素需要一次乘法和一次减法，所以第一阶段大约需要 $n^2$ 次乘法和 $n^2$ 次减法，实际上会少一些，因为第一行不变，实际上需要 $n(n-1)$ 次乘法和加法，这里计算的是一个大致成本。
第二阶段我们需要将列 $2$ 的第二主元下方的元素变为 $0$，此时我们要考虑的矩阵会小一些，是一个 $(n-1)\times (n-1)$ 的矩阵，所以这一阶段大约是 $(n-1{)}^2$ 次乘法与减法。越往下进行所要考虑的矩阵越小，最终要得到矩阵 $U$ 则粗略估计需要的次数为 $n^2+(n-1{)}^2+\cdots +2^2+1^2$。
上式平方和的公式为 $\frac{1}{3}n(n+\frac{1}{2})(n+1)$，当 $n$ 很大时，就可以忽略里面的 $\frac{1}{2}$ 和 $1$，总和大约就是 $\frac{1}{3}{n}^3$。$x^2$ 从 $0$ 到 $n$ 的积分就是 $\frac{1}{3}{n}^3$，需要注意的是积分是连续的，而这里是离散的。

$$对矩阵A使用消元法大概需要\frac{1}{3}{n}^3次乘法和\frac{1}{3}{n}^3次减法$$

下面考虑右侧的 $\mathbf{b}$，我们要计算 $L\mathbf{c}=\mathbf{b}$，得到 $\mathbf{c}$。首先，我们从 $b_2,\cdots ,b_n$ 减去乘数乘 $b_1$，这里需要 $n-1$ 步，第二阶段就不需要考虑 $b_1$，共需要 $n-2$ 步，最后一阶段需要 $1$ 步。
最后考虑回代，通过 $U\mathbf{x}=\mathbf{c}$ 求解 $\mathbf{x}$。首先，计算 $x_n$ 需要 $1$ 步，仅需要除以最后一个主元；然后计算 $x_{n-1}$ 需要 $2$ 步，这里需要代入 $x_n$，然后除以第 $n-1$ 主元；最后计算 $x_1$ 时需要 $n$ 步，要代入 $(n-1)$ 个未知数，然后除以第一主元。所有计算右侧的 $\mathbf{b}$ 正好需要需要 $n^2$ 步（从前到后再回代）：$$[(n-1)+(n-2)+\cdots +1]+[1+2+\cdots +(n-1)+n]=n^2(\mathrm{2.6.6})$$可以看到，右侧的成本要比左侧小很多。

$$求解右侧需要n^2次乘法和n^2次减法$$

一个带状矩阵 $B$ 只在主对角线的上方和下方有 $w$ 个非零对角线，带状外的其它元素在消元过程中都保持 $0$ 不变（$L$ 和 $U$ 中）。第一列需要 $w^2$ 次乘法和减法（在主元下方产生 $w$ 个 $0$，每个 $0$ 需要使用长度为 $w$ 的主元行），所以要执行完消元过程共需要 $n{w}^2$ 次乘法和减法，得到 $U$。这样会节省很多时间。

$$带状矩阵A到U\frac{1}{3}{n}^3减少到n{w}^2求解n^2减少到2nw$$

一个三对角矩阵（$w=1$）可以计算的很快，不需要存储 $0$。

五、主要内容总结

1. 高斯消元法（没有行交换，且是可逆矩阵）将 $A$ 分解成 $L$ 乘 $U$。
2. 下三角矩阵 $L$ 包含用来乘主元行的乘数 $l_{ij}$，它们使得 $A$ 变成 $U$。乘积 $LU$ 将这些行反向加回去可以恢复成 $A$。
3. 在右侧我们求解 $L\mathbf{c}=\mathbf{b}$（前向），然后求解 $U\mathbf{x}=\mathbf{c}$（反向）。
4. 分解： 左侧共有 $\frac{n^3-n}{3}$ 次乘法和减法（这个结果没有取近似）。
5. 求解： 右侧共有 $n^2$ 次乘法和减法。
6. 对于带状矩阵，需要的步骤 $\frac{1}{3}{n}^3$ 变为 $n{w}^2$，$n^2$ 变为 $2nw$。

六、例题

【例4】下三角帕斯卡矩阵 $L$ 包含著名的 “帕斯卡三角形”，这里我们分解帕斯卡。
对称帕斯卡矩阵 $P$ 是帕斯卡矩阵 $L$ 和 $U$ 的乘积。对称的 $P$ 矩阵以帕斯卡三角命名，所以它的每个元素都是其上方和左侧元素之和。MATLAB 中，$n\times n$ 的对称 $P$ 矩阵写成 pascal(n)。
问题：建立一个下 - 上三角分解的 $P=LU$。$$\text{pascal(4)}=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 2& 3& 4\\ 1& 3& 6& 10\\ 1& 4& 10& 20\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 1& 0& 0\\ 1& 2& 1& 0\\ 1& 3& 3& 1\end{array}\right]\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 1& 2& 3\\ 0& 0& 1& 3\\ 0& 0& 0& 1\end{array}\right]=LU$$预测并检验 $5\times 5$ 的帕斯卡矩阵的下一行和下一列。
解： 计算 $LU$ 可以得到 $P$。下面从对称 $P$ 矩阵出发，利用消元法得到上三角矩阵 $U$：$$P=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 2& 3& 4\\ 1& 3& 6& 10\\ 1& 4& 10& 20\end{array}\right]\to \left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 1& 2& 3\\ 0& 2& 5& 9\\ 0& 3& 9& 19\end{array}\right]\to \left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 1& 2& 3\\ 0& 0& 1& 3\\ 0& 0& 3& 10\end{array}\right]\to \left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 1& 2& 3\\ 0& 0& 1& 3\\ 0& 0& 0& 1\end{array}\right]=U$$上面步骤所用到的乘数 $l_{ij}$ 会进入到下三角矩阵 $L$，$P=LU$ 是一个特别整洁有序的例子。注意到 $U$ 的在对角线上的主元都是 $1$。
若使用 MATLAB 来计算，指令 lu(pascal(4)) 无法生成上述的 $U$，这是因为 lu 的子程序会在每一列选取最大的主元，这样第二主元就变成了 $3$，而不是 $1$，但是使用乔里斯基（Cholesky）分解不会发生行交换，可以产生上述结果：U = chol(pascal(4))

【例5】问题：求解 $P\mathbf{x}=\mathbf{b}=(1,0,0,0)$。方程的右侧等于 $I$ 的第一列，这表明 $\mathbf{x}$ 会是 $P^{-1}$ 的第一列。这就是高斯 - 若尔当消元法，会匹配 $P{P}^{-1}=I$ 的列。我们已知帕斯卡矩阵 $L$ 和 $U$ 是 $P$ 的两个因子：$$两个三角系统L\mathbf{c}=\mathbf{b}(前向)U\mathbf{x}=\mathbf{c}(后向)$$解： 下三角系统 $L\mathbf{c}=\mathbf{b}$ 由上到下求解$$\begin{array}{c}{c}_1=1\\ c_1+c_2=0\\ c_1+2c_2+c_3=0\\ c_1+3c_2+3c_3+c_4=0\end{array}解得\begin{array}{c}{c}_1=+1\\ c_2=-1\\ c_3=+1\\ c_4=-1\end{array}$$利用 $L^{-1}$ 执行前向消元，得到上三角形系统 $U\mathbf{x}=\mathbf{c}$，使用回代求解 $\mathbf{x}$，上三角系统由下到上求解：$$\begin{array}{c}{x}_1+x_2+x_3+x_4=1\\ x_2+2x_3+3x_4=-1\\ x_3+3x_4=1\\ x_4=-1\end{array}解得\begin{array}{c}{x}_1=+4\\ x_2=-6\\ x_3=+4\\ x_4=-1\end{array}$$使用 inv(pascal(4)) 指令求得 $P$ 的逆矩阵，可以看到解就是 $P^{-1}$ 的第一列。

【例6】对于 $6\times 6$ 的第二差分常对角（constant-diagonal）矩阵 $K$，将主元和乘数放入 $K=LU$。（$L$ 和 $U$ 有两个非零对角线，因为 $K$ 有三个）。找到 $L^{-1}$ 的 $i,j$ 元素的公式。
**解：

format rat;    				 % 切换成分数显示
K = toeplitz([2 -1 0 0 0 0]) % 生成常对角矩阵，特普利兹 toeplitz 矩阵，这里参数是其第一行
[L U] = lu(K);				 % 对 K 进行 LU 分解
inv(L)						 % 求 L 的逆矩阵

从 MATLAB 运行的结果可以得出其 $L^{-1}$ 的 $i,j$ 位置元素的公式是：$$\frac{j}{i},i\ge j$$