MIT线性代数总结笔记——LU分解

矩阵分解

矩阵分解（Matrix Factorizations）就是将一个矩阵用两个以上的矩阵相乘的等式来表达。而矩阵乘法涉及到数据的合成(即将两个或多个线性变换的效果组合成一个矩阵)，因此可以说，矩阵分解是对数据的一种分析。

矩阵分解	LU分解
分解形式	$A=LU$ ( $L$ 代表下三角矩阵， $U$ 代表上三角矩阵)
目的	提高计算效率
前提	(1)矩阵是方阵（ $LU$ 分解主要是针对方阵）； (2)矩阵是可逆的，也就是该矩阵是满秩矩阵，每一行都是独立向量； (3)消元过程中没有 $0$ 主元出现，也就是消元过程中不能出现行交换的初等变换。

LU分解

分解形式

$A = LU = \left[\begin{matrix} 1 & 0 & 0 & 0 \\ l_{21} & 1 & 0 & 0 \\ l_{31} & l_{32} & 1 & 0 \\ l_{41} & l_{42} & l_{43} & 1\end{matrix}\right]\left[\begin{matrix} u_{11} & u_{12} & u_{13} & u_{14} \\ 0 & u_{22} & u_{23} & u_{24} \\ 0 & 0 & u_{33} & u_{34} \\ 0 & 0 & 0 & u_{44}\end{matrix}\right]$

简述

$LU$ 分解常用于求解工业和商业问题中的序列方程。它是最常见的求解线性系统 $Ax=b$ 的方法，主要思路是：把 $A$ 分解成一个下三角矩阵（Lower Triangular Matrix）和一个上三角矩阵（Upper Triangular Matrix），简称 $LU$ 。分解后，等式 $Ax=b$ 可以写成 $L(Ux)=b$ ，这样我们令 $y=Ux$ ，这样就可以通过分开求解两个等式来得到 $x$ ：

$Ly=b \\ Ux=y$

即可以先通过 $Ly=b$ 求出 $y$ ，最后再用 $Ux=y$ ，求出 $x$ 。

本质上， $LU$ 分解是高斯消元法的一种表达方式，为了更好地理解 $LU$ 分解，我们需要解释一下高斯消元法。

高斯消元法（Gaussian Elimination）

简述

为了求解线性系统，我们基本的策略是用相等的式子去替代要求解的式子，来让求解变得更容易。对于一个线性系统中有若干个未知数，我们需要将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其代入到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的，并最终得到线性系统的解。这一求解算法称之为高斯消元法（Gaussian Elimination）。

示例

下文将展示高斯消元法的求解过程，设需求解的线性方程组如下：

$\begin{cases} x_1-2x_2+x_3=0 \\ 2x_2-8x_3=8\\5x_1-5x_3=10 \end{cases}$

将其表达为矩阵形式为：

$\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]$

为了将多余的未知数 $x_1$ 消除，我们需要将第三行的 $5$ 消除，矩阵的第一行乘 $-5$ 倍再和第三行相加，得到：

$\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10 \end{array}\right]$

接着再进行对多余的未知数 $x_2$ 的消除，这里将第二行乘 $-5$ 倍再与第三行相加，得到：

$\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]$

至此，我们得到了一个新的三角形矩阵，转换为线性方程组的形式为：

$\begin{cases} x_1-2x_2+x_3=0 \\ x_2-4x_3=4\\x_3=-1 \end{cases}$

通过简单的回代，我们最终可以得到线性方程组的解：

$\begin{cases} x_1=1 \\ x_2=0\\x_3=-1 \end{cases} \quad\quad\quad \left[\begin{array}{ccc|c} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & -1 \end{array}\right]$

高斯消元法与LU分解之间的联系

那么为什么说 $LU$ 分解其实就是高斯消元法的一种表达方式呢？这是因为我们可以将上述示例中每一步的操作用矩阵相乘的形式（行变换的本质就是左乘矩阵）表达出来。

例如上文示例中的第一步：“为了将多余的未知数 $x_1$ 消除，我们需要将第三行的 $5$ 消除，矩阵的第一行乘 $-5$ 倍再和第三行相加” ，可以表达为：

$\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & 0 & 1 \end{matrix}\right] \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]=\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10 \end{array}\right]$

同样地，第二步：“继续进行对多余的未知数 $x_2$ 的消除，这里将第二行乘 $-5$ 倍再与第三行相加” 可以表达为：

$\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -5 & 1 \end{matrix}\right] \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10 \end{array}\right] = \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]$

因此，我们可以用一个矩阵相乘的等式来记录整个高斯消元的过程：

$\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -5 & 1 \end{matrix}\right]\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & 0 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]=\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]$

将三个消元矩阵合并，用 $E$ 表示则为：

$EA =\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & -5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right]=\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right] = U$

那么我们可以得到关于 $A$ 的表达式：

$A = \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right] = E^{-1}U = \left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right] = LU$

这样一来我们就得到了 $A$ 的 $LU$ 分解，即将矩阵A分解为一个下三角矩阵（ $L$ ）和一个上三角矩阵（ $U$ ）的乘积。

伪代码

U = A, L = I
for j = 1 : n - 1 do
    for i = j + 1 : n do
        l_ij = u_ij / u_jj
        for k = j : n do
            u_ik = u_ik - l_ij * u_jk
        end for
    end for
end for

此处是矩阵index从1开始计数。
Credit to http://iacs-courses.seas.harvard.edu/courses/am205/slides/am205_lec07.pdf

算法效率

接下来我们来思考一下算法效率的问题，假如我们有一个 $100*100$ 的矩阵，并且各项均非 $0$ ，那我们需要计算多少次才能得到矩阵 $U$ 呢？

对于这个问题我们先从列的角度来考虑，第一列消元运算结束以后，矩阵变为：

$\left[\begin{matrix} 1 & \cdots & \cdots & \cdots \\ 0 & \ddots & \cdots & \cdots \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & \cdots & \ddots\end{matrix}\right]$

这一步中，第一列的元素运算了 $100$ 次，而第一行共有 $100$ 个元素，于是仅第一行与第一列消元结束后，我们就计算了 $100^2$ 次。之后我们要研究剩下的 $99*99$ 的矩阵，以此类推可知，最后的计算量为 $\sum^{n}_{k=1}{k^2}$ 。

接着，我们可以利用微积分计算得到， $A$ 的操作次数为 $\frac{1}{3}n^3$ ， $b$ 的操作次数为 $n^2$ 次。

结合上文伪代码，因为有三重循环嵌套，因此 $LU$ 分解的复杂度为 $O(n^3)$ ，其中加法操作为 $\frac{1}{3}n^3$ ，乘法操作为 $\frac{1}{3}n^3$ ，总共为 $\frac{2}{3}n^3$ 次，求解 $b$ 和 $y$ 各需要 $n^2$ 次

前提条件

(1)矩阵是方阵（LU分解主要是针对方阵）；
(2)矩阵是可逆的，也就是该矩阵是满秩矩阵，每一行都是独立向量；
(3)消元过程中没有0主元出现，也就是消元过程中不能出现行交换的初等变换。

LUD分解（LDU Decomposition）

上文中我们已经得到了 $A=LU$ ，我们还可以继续分解下去得到 $A=LDU$ ， $D$ 为对角矩阵(Diagonal Matrix)，矩阵的 $LDU$ 分解是在 $LU$ 分解之后，把 $U$ 再次分解，目的是把 $U$ 的对角线元素都化为 $1$ ，上文示例中我们得到的 $A=LU$ 如下

$A = \left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 &10 \end{array}\right] = LU = \left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right]$

接下来，让 $U$ 矩阵中位于对角线上的元素都为 $1$ ，就可以提出一个对角矩阵 $D$

$LU = \left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30 \end{array}\right] =\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1 \end{matrix}\right]\left[\begin{matrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 30 \end{matrix}\right]\left[\begin{array}{ccc|c} 1 & -2 & 1 & 0 \\ 0 & 1 & -4 & 4 \\ 0 & 0 & 1 & -1 \end{array}\right]=LDU$

LUP分解（LU Decomposition with Partial Pivoting）

我们来考虑一个矩阵 $A=\left[\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right]$ ，尽管矩阵 $A$ 非奇异，并且很简单，但是使用 $LU$ 分解的算法会失败，因为第一个主元 $A_{0,0} = 0$ 。在高斯消元法中，我们常常通过人为的行变换来消除主元为 $0$ 时的影响。例如我们需要把矩阵A的第一行和第二行进行交换得到 $A=\left[\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right]$ ，这样我们就能将消元进行下去。这样的操作也可以用矩阵相乘的方法来表达，那么在该示例中，即为 $PA=\left[\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}\right]\left[\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right]=\left[\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right]$ ，这里的矩阵 $P$ 称为置换矩阵（Permutation Matrix）。那么就有了LUP分解（LU decomposition with partial pivoting）： $PA = LU$ 。这一算法的出现显著改善了 $LU$ 分解算法的稳定性。

伪代码

U = A, L = I, P = I
for j = 1 : n - 1 do
    Select i (>= j) that maximizes |u_ij|
    Interchange rows of U: u_(j,j:n) <-> u_(i,j:n)
    Interchange rows of L: l_(j,1:j-1) <-> l_(i,1:j-1)
    Interchange rows of P: p_(j,:) <-> p_(i,:)
    for i = j + 1 : n do
        l_ij = u_ij / u_jj
        for k = j : n do
            u_ik = u_ik - l_ij * u_jk
        end for
    end for
end for