矩阵A的LU分解 A=LU

news/2024/12/2 22:18:25/文章来源:https://www.cnblogs.com/zzzwwwqqq/p/18582862

矩阵$\mathbf A$的LU分解 $\mathbf A=\mathbf L\mathbf U$

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

先给出结论^[1]：$\mathbf{AB}$的逆矩阵是$\mathbf{B}^{-1}\mathbf{A}^{-1}$。简单证明一下：$(\mathbf{AB})(\mathbf{B}^{-1}\mathbf{A}^{-1})=\mathbf I$，应用乘法结合律便可以得到这一结论。Gilbert Strang教授在他的书中做了一个形象的比喻：就好比先脱鞋子，再脱袜子，它的“逆”就是先穿袜子后穿鞋子。对于两步操作（$\mathbf A$和$\mathbf 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}$的逆运算也不过是把该运算减去的数加回去而已，很容易求出来），不过这个的前提是没有行互换。

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

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

首先说一下Gilbert Strang教授的说法。我们对一个$n\times n$的矩阵$\mathbf A$，需要多少操作，得到上三角矩阵$\mathbf 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}$步。

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

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

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

首先，看下述公式

\[\mathbf A\mathbf x=\mathbf L\mathbf U\mathbf x=\mathbf L(\mathbf U\mathbf x)=\mathbf b\\ \Rightarrow \\ \\ \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} \]

这样就完成了对$\mathbf y$的求解。再代回求解$\mathbf 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}$。因此，如果矩阵$\mathbf A$已化为$\mathbf L\mathbf U$的形式，那么便可以通过以下两步完成对$\mathbf x$的求解：

根据$\mathbf L\mathbf y=\mathbf b$求解出$\mathbf y$；
根据$\mathbf U\mathbf x=\mathbf y$将$\mathbf y$代入，求解出$\mathbf x$。

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

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

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