给定 $n$ 阶方阵 $A$ 与系数向量 $\boldsymbol b \in \mathbb R^n$，求解满足

$$A \boldsymbol x = \boldsymbol b$$

的未知向量 $\boldsymbol x$ 的问题被称为线性方程式问题。
本节讨论计算机中如何求解该类问题

线性代数中，如果 $\det A \neq 0$，则可以写出解

$$\boldsymbol x = A^{-1} \boldsymbol b$$

但实际应用上，由于阶数 $n$ 可能非常巨大，使得逆矩阵计算困难。
此时需要考虑基于 Gauss 消元法的求解

如果系数矩阵是一个上三角或者下三角矩阵，那么求解该方程将极其容易
将系数矩阵 $A$ 改写为上三角矩阵 $U$ 与下三角矩阵 $L$ 的乘积，即 $A = LU$，称为 LU 分解
此时，方程可以改写为

$$L(U \boldsymbol x) = \boldsymbol b$$

我们可以很轻松地求解出 $U \boldsymbol x$，进而求解出 $\boldsymbol x$。

并且，我们可以要求下三角矩阵 $L$ 的对角线元素全为 $1$，以此来保证分解的唯一性

接下来让我们看看如何实现 LU 分解

# 无 Pivoting 的 LU 分解

LU 分解基于 Gauss 消元，而消元法中要求每一步所用的主元都不能为 $0$
如果出现了这种情况，就需要进行行交换，称为 Pivoting
我们先暂时考虑不会遇到所用主元为 $0$ 的情况，也就是说对前 k 阶主子式都要求非零

$$\det(A_1)\neq 0,\quad \det(A_2)\neq 0,\quad \dots,\quad \det(A_n)\neq 0$$

假设我们已经实现了 LU 分解，也就是 $A = LU$
对 $L = (\ell_{ij})$ 进行列向量分解，对 $U = (u_{ij})$ 进行行向量分解，得到

$$A = LU = \begin{pmatrix} \boldsymbol l_1 & \boldsymbol l_2 & \cdots & \boldsymbol l_n \end{pmatrix} \begin{pmatrix} \boldsymbol u_1^T \\ \boldsymbol u_2^T \\ \vdots \\ \boldsymbol u_n^T \end{pmatrix} = \sum_{k=1}^n \boldsymbol l_k \boldsymbol u_k^T$$

我们关注求和项，得到

$$\begin{aligned} \boldsymbol l_k \boldsymbol u_k^T &= \begin{pmatrix} 0 \\ 0 \\ \vdots \\ \ell_{kk} \\ \vdots \\ \ell_{nk} \end{pmatrix} \begin{pmatrix} 0 & 0 & \cdots & u_{kk} & \cdots & u_{kn} \end{pmatrix} \\ &= \begin{pmatrix} 0 & 0 & \cdots & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \ell_{kk} u_{kk} & \cdots & \ell_{kk} u_{kn} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \ell_{nk} u_{kk} & \cdots & \ell_{nk} u_{kn} \end{pmatrix} \end{aligned}$$

观察不难得到，每一个第 $k$ 个求和项实际上的贡献在于右下角开始的一个边长 $\ell(k)$ 的方阵
因此，通过这样的逐级分解，就可以在该种简单情况下得到 LU 分解

具体实现方法如下：
无 Pivoting 的 LU 分解算法

• 关注 $A$ 的第一行第一列，这部分只能由单一的求和项给出，因此可以确定 $\boldsymbol \ell_1$ 与 $\boldsymbol u_1^T$ 的值。只需要将 $\boldsymbol u_1^T$ 指定为 $A$ 的第一行，再相除得到 $\boldsymbol \ell_1$
• 令 $A_1 = A - \boldsymbol \ell_1 \boldsymbol u_1^T$
• 关注 $A_1$ 的第二行第二列，这部分只能由单一的求和项给出，因此可以确定 $\boldsymbol \ell_2$ 与 $\boldsymbol u_2^T$ 的值。只需要将 $\boldsymbol u_2^T$ 指定为 $A_1$ 的第二行，再相除得到 $\boldsymbol \ell_2$
• 以此类推

以上，通过 $n-1$ 次除法，$n(n^2-1)/3$ 次乘法，$n(n^2-1)/3$ 次减法，就可以得到 LU 分解。全体的计算量为 $2n^3/3$，因此时间复杂度为 $O(n^3)$

算法实现

for(k=1; k<=n-1; k++){

    w = 1/a[k][k];

    for(i=k+1; i<=n; i++){

        a[i][k] = a[i][k] * w;

        for(j=k+1; j<=n; j++){

            a[i][j] = a[i][j] - a[i][k] * a[k][j];

        }

    }

}

这样的 LU 分解有一个极其显著的优点：不占用内存

在完成 LU 分解后，就可以快速求解方程了
令 $\boldsymbol y = U \boldsymbol x$，我们先求解 $L \boldsymbol y = \boldsymbol b$
$L$ 是一个下三角矩阵，通过每一次向矩阵的右侧移动消元，可以得到 $\boldsymbol y$ 的值。这个方法叫做 前代 (Forward Substitution)「前進消去」。
算法实现

for(i=1; i<=n; i++){

    s = 0;

    for(j=1; j<=i-1; j++){

        s = s + a[i][j] * y[j];

    }

    y[i] = b[i] - s;

}

接下来求解 $U \boldsymbol x = \boldsymbol y$，$U$ 是一个上三角矩阵，通过每一次向矩阵的左侧移动消元，可以得到 $\boldsymbol x$ 的值。这个方法叫做 回代 (Backward Substitution)「後退代入」。
算法实现

for(i=n; i>=1; i--){

    s = 0;

    for(j=i+1; j<=n; j++){

        s = s + a[i][j] * x[j];

    }

    x[i] = (y[i] - s) / a[i][i];

}

最终我们就可以得到方程的解 $\boldsymbol x$ 了

# 部分 Pivoting 的 LU 分解