矩阵A的LU分解 A=LU

矩阵\(\symbfit A\)的LU分解 \(\symbfit A=\symbfit L\symbfit U\)

​ 如果已知\(\symbfit A\)和\(\symbfit B\)的逆矩阵\(\symbfit A^{-1}\)和\(\symbfit B^{-1}\)，那么\(\symbfit {AB}\)的逆是什么？

​ 先给出结论^[1]：\(\symbfit {AB}\)的逆矩阵是\(\symbfit B^{-1}\symbfit {A}^{-1}\)。简单证明一下：\((\mathbf{AB})(\symbfit {B}^{-1}\symbfit {A}^{-1})=\symbfit I\)，应用乘法结合律便可以得到这一结论。Gilbert Strang教授在他的书中做了一个形象的比喻：就好比先脱鞋子，再脱袜子，它的“逆”就是先穿袜子后穿鞋子。对于两步操作（\(\symbfit A\)和\(\symbfit B\)如果右乘一个矩阵的话，自然可以理解为对这个矩阵的操作），如果想要逆运算，那么必然是先逆后一步，再前者。这一结论同样适用于三个及以上的矩阵。

​ 书归正传，那么什么是矩阵的LU分解？其中的“U”可以猜到，在之前的篇章中也有提到，是矩阵经消元之后的结果，不过并没有提到它为什么叫“U"，在本篇随后会有说明。

​ 让我们回到之前讲消元矩阵的时候，以下面这个运算为例：

\[\begin{bmatrix} 1 & 0 \\ -4 & 1 \end{bmatrix} \begin{bmatrix} 2 & 1 \\ 8 & 7 \end{bmatrix} = \begin{bmatrix} 2 & 1 \\ 0 & 3 \end{bmatrix} \]

对矩阵\(\mathbf A\)的\(row2 - 4row1\)，即左乘\(\mathbf E_{2,1}\)，得到\(\mathbf U\)。那如果想把这个等式化成类似\(\mathbf A=\mathbf L\mathbf U\)的形式，那么问题是\(\mathbf L\)是什么，即要得到

\[\begin{bmatrix} 2 & 1 \\ 8 & 7 \end{bmatrix} = \mathbf L \begin{bmatrix} 2 & 1 \\ 0 & 3 \end{bmatrix} \]

不难得出，\(\mathbf L=\begin{bmatrix}1&0\\4&1\end{bmatrix}\)，而且\(\mathbf L=\mathbf E_{2,1}^{-1}\)。重新写一下

\[\begin{bmatrix} 2 & 1 \\ 8 & 7 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 4 & 1 \end{bmatrix} \begin{bmatrix} 2 & 1 \\ 0 & 3 \end{bmatrix} \]

观察到，同时把矩阵\(\mathbf A\)化成了一个下三角矩阵（lower）与一个上三角矩阵（upper）的乘积的形式，这也就是\(\mathbf L\)和\(\mathbf U\)的来历。消元的目的也正是将矩阵化为上三角矩阵的形式。右式也可以表示为

\[\begin{bmatrix} 1 & 0 \\ 4 & 1 \end{bmatrix} \begin{bmatrix} 2 & 0 \\ 0 & 3 \end{bmatrix} \begin{bmatrix} 1 & 1/2 \\ 0 & 1 \end{bmatrix} \]

其中中间的矩阵为对角矩阵，记为\(\mathbf D\)。

​ 再考虑一下\(\mathbf A\)为\(3\times3\)的情况。在假设没有行交换的情况下，有\(\mathbf{E}_{3,2}\mathbf{E}_{3,1}\mathbf{E}_{2,1}\mathbf{A}=\mathbf{U}\)，那么就有

\[\mathbf{A}=\mathbf{E}_{2,1}^{-1}\mathbf{E}_{3,1}^{-1}\mathbf{E}_{3,2}^{-1}\mathbf{U} \]

举个例子，假设\(\mathbf{E}_{3,2}=\begin{bmatrix}1&0&0\\0&1&0\\0&-5&1\end{bmatrix}\)，\(\mathbf{E}_{2,1}=\begin{bmatrix}1&0&0\\-2&1&0\\0&0&1\end{bmatrix}\)，没有用到\(\mathbf{E}_{3,1}\)。那么

\[\mathbf{E}_{3,2}\mathbf{E}_{2,1}= \begin{bmatrix} 1 & 0 & 0\\ -2 & 1 & 0\\ 10 & -5 & 1 \end{bmatrix} =\mathbf{E} \]

如果按照之前对消元矩阵的理解去分析，会发现关于矩阵\(\mathbf A\)的$ row3 \(处理是\) 10row1-5row2+row3\(，这是两个运算综合起来的结果：\) row2-2row1 \(和\) row3-5row2$。如果要计算其逆矩阵，不难得出

\[\begin{bmatrix} 1 & 0 & 0\\ 2 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 5 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0\\ 2 & 1 & 0\\ 0 & 5 & 1 \end{bmatrix} =\mathbf{L} \]

​ 注意到，在\(\mathbf{L}\)中，与\(\mathbf{E}\)不同，矩阵相乘时，\(2\)与\(5\)没有相互干扰得到\(10\)。因此，要想求出\(\mathbf{L}\)，只需把所有消元乘数都写进来即可（实际上对\(\mathbf{E}_{3,2}\)和\(\mathbf{E}_{2,1}\)的逆运算也不过是把该运算减去的数加回去而已，很容易求出来），不过这个的前提是没有行互换。

​ 关于这个结论，也不难理解，因为正运算（\(\symbfit{E}\)）从上到下的消元过程，前面的步骤都会影响后续的步骤；而逆运算（\(\symbfit{L}\)）是从下到上的，不会影响到上面的向量（行）。

---

​ 上面正文部分已结束，下面将讨论LU分解的意义，即我们为什么要做LU分解。

​ 首先说一下Gilbert Strang教授的说法。我们对一个\(n\times n\)的矩阵\(\symbfit A\)，需要多少操作，得到上三角矩阵\(\symbfit U\)？假设其中没有\(0\)，因为有\(0\)的话计算次数会减少，考虑最复杂的情况。

​ 如果对第一个主元做处理，让它的下面全为\(0\)，即实现

\[\begin{bmatrix} \boxed{a_{0,0}} & \cdots & a_{1,n}\\ 0 & \cdots & a_{2,n}\\ \vdots & \ddots & \vdots\\ 0 & \cdots & a_{n,n} \end{bmatrix} \]

下面假设先做一次乘法再做一次减法记为一次操作（一步），那么刚才所说的得到第一主元需要\(n\times (n-1)\)步。其下同理，一共需要

\[n\times(n-1)+(n-1)\times(n-2)+\cdots+2\times1=\sum_{i=1}^{n-1}{i(i-1)} \approx\sum_{i=1}^{n-1}{i^2} \approx\dfrac{1}{3}n^3 \]

在\(n\)很大时可以做这个近似。

​ 那如果只对增广矩阵最右侧\(b\)做运算的话，与之同理，需要\(\dfrac{n^2}{2}\)步。

​ 由此我们得出，对形似\(\symbfit A \symbfit x=\symbfit b\)的计算，对\(\symbfit A\)的操作的复杂度是很高的，如果提前将\(\symbfit A\)变为\(\symbfit L\symbfit U\)的形式，那么对于任意的\(\symbfit b\)，计算就容易多了。

​ 最后再说一下在网上找到的矩阵的LU分解的过程以及意义 - 知乎这篇文章，作者的思路也很好。

​ 假设现在有一个系统，系统内涉及到一个运算\(\symbfit A \symbfit x=\symbfit b\)，根据\(\symbfit b\)来求取出信息\(\symbfit x\)，每次都是\(\symbfit b\)不相同，但是\(\symbfit A\)是固定的。这个时候我们就可以通过LU分解，预先将固定的矩阵\(\symbfit A\)分解，然后通过LU来求解\(\symbfit x\)。之所以这么做，是因为LU分解比直接求\(\symbfit x\)的效率要高。

​ 首先，看下述公式

\[\mathbf A\mathbf x=\mathbf L\mathbf U\mathbf x=\mathbf L(\mathbf U\mathbf x)=\mathbf b\\ \Downarrow \\ \mathbf L\mathbf y=\mathbf b,\text{其中}\mathbf y=\mathbf U\mathbf x \]

其第二行可以写成如下形式

\[\begin{bmatrix} l_{1,1}\\ \vdots & \ddots\\ l_{1,n} & \cdots & l_{n,n} \end{bmatrix} \begin{bmatrix} y_1\\ \vdots\\ y_n \end{bmatrix} = \begin{bmatrix} b_1\\ \vdots\\ b_n \end{bmatrix} \]

对于上面的公式，通过逐上自下的方式逐一求取未知数，其运算次数为

\[1+2+\cdots+n\approx \frac{n^2}{2} \]

这样就完成了对\(\symbfit y\)的求解。再代回求解\(\symbfit x\)

\[\begin{bmatrix} u_{1,1} & \cdots & u_{1,n}\\ & \ddots & \vdots\\ & & u_{n,n} \end{bmatrix} \begin{bmatrix} x_1\\ \vdots\\ x_n \end{bmatrix} = \begin{bmatrix} y_1\\ \vdots\\ y_n \end{bmatrix} \]

其运算次数同理，为\(\dfrac{n^2}{2}\)。因此，如果矩阵\(\symbfit A\)已化为\(\symbfit L\symbfit U\)的形式，那么便可以通过以下两步完成对\(\symbfit x\)的求解：

1. 根据\(\symbfit L\symbfit y=\symbfit b\)求解出\(\symbfit y\)；
2. 根据\(\symbfit U\symbfit x=\symbfit y\)将\(\symbfit y\)代入，求解出\(\symbfit x\)。

​ 因此整个过程复杂度为\(\dfrac{n^2}{2}+\dfrac{n^2}{2}=n^2\)。

​ 而直接利用消元法的复杂度在之前已经求解，是\(\dfrac{n^3}{3}+\dfrac{n^2}{2}\)。显然，利用\(\symbfit L\symbfit U\)的形式复杂度更小。

​ 举个不是很恰当的例子，利用\(\symbfit A=\symbfit L\symbfit U\)将\(\symbfit A\)分解的办法就好比代码编译后，根据对其输入的值能很快得到输出一样；而每次直接消元运算就像是将输入参数的值嵌入到程序中，每次运行都需要重新编译。再说一个更不恰当的例子，LU分解就得到了矩阵\(\symbfit A\)的“系统函数”，而不同的输入只需对其做相应的卷积（或乘积）即可。而且，根据正文的讨论，得到\(\symbfit{L}\)和\(\symbfit{U}\)的过程其实并不麻烦，相反，比直接得到\(\symbfit{E}\)要轻松得多。

---

1. 这本应是上一篇给出的，在这一篇里会用到，于是直接在这一篇开头说明了。 ↩︎