2.1 LU分解与Gauss消去法

2.1.1 LU分解

Gauss 消去法本质上就是对系数矩阵 $A$ 进行 LU 分解, 即将 $A$ 分解成两个矩阵的乘积

其中 $L$ 是单位下三角矩阵, $U$ 为非奇异上三角矩阵, 然后再对两个三角方程进行求解: 假定矩阵 $A$ 存在 LU 分解 (2.2), 则方程组 (2.1) 就转化为求解下面两个三角方程组

显然，上式中的两个三角方程组都非常容易求解

A 将一个复杂问题分解成若干相对简单的问题, 是矩阵分解的一个主要目的.

基于 LU 分解的 Gauss 消去法描述如下:

算法2.1.Gauss消去法

1: 对 $A$ 进行 LU 分解: $A = LU$ , 其中 $L$ 为单位下三角矩阵, $U$ 为非奇异上三角矩阵;
2: 利用向前回代, 求解 $Ly = b$ , 即得 $y = L^{-1}b$ ;
3. 利用向后回代, 求解 $Ux = y$ , 即得 $x = U^{-1}y = (LU)^{-1}b = A^{-1}b$

我们知道, 当系数矩阵 $A$ 非奇异时, 方程组 (2.1) 总是存在唯一解. 但是, 并不是每个非奇异矩阵都存在 LU 分解.

定理2.1 (LU分解的存在性和唯一性) 矩阵 $A \in \mathbb{R}^{n \times n}$ 存在LU分解（即存在单位下三角矩阵 $L$ 和非奇异上三角矩阵 $U$ 使得 $A = LU$ ）的充要条件是 $A$ 的所有顺序主子矩阵都非奇异。进一步，若 $A$ 存在LU分解，则分解是唯一的。

(板书)

证明. 必要性: 设 $A_{11}$ 是 $A$ 的 $k$ 阶顺序主子矩阵, 其中 $1 \leq k \leq n$ . 将 $A = LU$ 写成分块形式, 即

可得 $A_{11} = L_{11}U_{11}$ 由于 $L_{11}$ 和 $U_{11}$ 均非奇异, 所以 $A_{11}$ 也非奇异.

充分性: 用归纳法.

当 $n = 1$ 时, 结论显然成立.

假设结论对所有 $n - 1$ 阶矩阵都成立, 即对任意 $n - 1$ 阶矩阵, 如果其所有的顺序主子矩阵都非奇异, 则存在 LU 分解.

考虑 $n$ 阶的矩阵 $A$ , 写成分块形式

其中 $A_{11} \in \mathbb{R}^{(n-1) \times (n-1)}$ 是 $A$ 的 $n-1$ 阶顺序主子矩阵. 由归纳假设可知, $A_{11}$ 存在 LU 分解, 即存在单位下三角矩阵 $L_{11}$ 和非奇异上三角矩阵 $U_{11}$ 使得

令

则

因此可得 $A$ 的LU分解 $A = LU$ ，其中 $L \triangleq \begin{bmatrix} L_{11} & 0 \\ L_{21} & 1 \end{bmatrix}$ 为单位下三角矩阵， $U \triangleq \begin{bmatrix} U_{11} & U_{12} \\ 0 & U_{22} \end{bmatrix}$ 为非奇异的上三角矩阵（ $U$ 的非奇异性可由 $A$ 的非奇异性可得）.

由归纳法可知, 结论成立.

下面证明唯一性. 设 $A$ 存在两个不同的 LU 分解:

其中 $L$ 和 $\tilde{L}$ 为单位下三角矩阵, $U$ 和 $\tilde{U}$ 为非奇异上三角矩阵. 则有

该等式左边为下三角矩阵, 右边为上三角矩阵, 所以只能是对角矩阵. 又单位下三角矩阵的逆仍然是单位下三角矩阵, 所以 $L^{-1}\tilde{L}$ 的对角线元素全是 1 , 故

即 $\tilde{L} = L, \tilde{U} = U$ . 唯一性得证

记 $D$ 为 $U$ 的对角线部分, 则 $A = L\tilde{D}\tilde{U}$ , 其中 $\tilde{U} = D^{-1}U$ 是单位上三角矩阵. 因此我们就有下面的LDU分解.

推论2.2(LDU分解)设 $A\in \mathbb{R}^{n\times n}$ 的所有顺序主子矩阵都非奇异,则 $A$ 存在LDU分解,即存在单位下三角矩阵 $L$ 和单位上三角矩阵 $U$ ，以及非奇异对角矩阵 $D$ ，使得 $A = LDU$ ,其中 $L,U,D$ 都是唯一的.反之，若 $A$ 存在LDU分解，则 $A$ 的所有顺序主子矩阵都非奇异

对角占优情形

一般的非奇异矩阵不一定存在 LU 分解, 但如果 $A$ 是对角占优的, 则存在 LU 分解. 我们可以证明, 严格对角占优矩阵在 LU 分解中保持严格对角占优性 (见课后练习), 因此通过数学归纳法可以证明严格对角占优矩阵一定存在 LU 分解. 事实上, 只要矩阵列对角占优且非奇异, 则一定存在 LU 分解 [57].

定理2.3 [57] 设 $A \in \mathbb{R}^{n \times n}$ 非奇异且列对角占优, 则 $A$ 存在 LU 分解且 $L$ 中的元素的绝对值都不超过1. (留作练习, 数学归纳法)

需要注意的是, 在实际计算中, 由于舍入误差的影响, 即使是严格对角占优矩阵, 如果对角占优性不是很强的话, 很有可能会失去对角占优性, 从而导致 LU 分解失败.

2.1.2 LU分解的实现

矩阵的 LU 分解可以通过初等行变换来实现. 给定矩阵

• 第一步: 假定 $a_{11} \neq 0$ , 构造矩阵

易知 $L_{1}$ 的逆为

用 $L_{1}^{-1}$ 左乘 $A$ , 并将所得到的矩阵记为 $A^{(1)}$ , 则

即左乘 $L_{1}^{-1}$ 后, $A$ 的第一列中除第一个元素外其它都变为 0.

• 第二步: 类似地, 我们可以将上面的操作作用在 $A^{(1)}$ 的子矩阵 $A^{(1)}(2:n,2:n)$ 上, 将其第一列除第一个元素外都变为 0. 也就是说, 假定 $a_{22}^{(1)} \neq 0$ , 构造矩阵

用 $L_{2}^{-1}$ 左乘 $A^{(1)}$ , 并将所得到的矩阵记为 $A^{(2)}$ , 则

• 依此类推, 假定 $a_{kk}^{(k - 1)} \neq 0 (k = 3,4,\ldots ,n - 1)$ , 则我们可以构造一系列的矩阵 $L_{3}, L_{4}, \ldots, L_{n - 1}$ , 使得

为一个上三角矩阵. 我们将这个上三角矩阵记为 $U$ , 并记

则可得

这就是 $A$ 的LU分解

将上述过程写成算法，描述如下：

算法2.2.LU分解

1: Set $L = I, U = 0$ % 将 $L$ 设为单位矩阵, $U$ 设为零矩阵
2: for $k = 1$ to $n - 1$ do
3: for $j = k$ to $n$ do
4: $u_{kj} = a_{kj} \quad \%$ 更新 $U$ 的第 $k$ 行
5: end for
6: for $i = k + 1$ to $n$ do
7: $l_{ik} = a_{ik} / a_{kk}$ $\%$ 计算 $l_{ik}$
8: for $j = k + 1$ to $n$ do
9: $a_{ij} = a_{ij} - l_{ik}u_{kj}$ $\%$ 更新 $A(i,k + 1:n)$
10: end for

11: end for
12: end for

评价算法的一个主要指标是执行时间, 但这依赖于计算机硬件和编程技巧等, 因此直接给出算法执行时间是不太现实的. 所以我们通常是统计算法中算术运算 (加减乘除) 的次数. 在矩阵计算中, 大多仅仅涉及加减乘除和开方运算. 一般情况下, 加减运算次数与乘法运次数具有相同的量级, 而除法运算和开方运算次数具有更低的量级.
为了尽可能地减少运算量, 在实际计算中, 数, 向量和矩阵做乘法运算时的先后执行次序为: 先计算数与向量的乘法, 然后计算矩阵与向量的乘法, 最后才计算矩阵与矩阵的乘法. 比如计算 $\alpha A B x$ , 其中 $\alpha$ 是数, $A, B$ 是矩阵, $x$ 是向量, 如果按照从左往右计算的话, 则运算量为 $\mathcal{O}(n^{3})$ , 但是如果先计算 $\alpha x$ , 然后计算 $B(\alpha x)$ , 最后再计算 $A(B(\alpha x))$ 的话, 运算量则为 $\mathcal{O}(n^{2})$ , 相差一个量级.

矩阵 $L$ 和 $U$ 的存储

当 $A$ 的第 $i$ 列 (严格下三角部分) 被用于计算 $L$ 的第 $i$ 列后, 在后面的计算中不再被使用. 而 $A$ 的第 $i$ 行 (上三角部分) 更新后就是 $U$ 的第 $i$ 行. 因此, 为了节省存储空间, 我们可以在计算过程中将 $L$ 的第 $i$ 列存放在 $A$ 的第 $i$ 列 (严格下三角部分, $L$ 的对角线全部为 1, 不需要存储), 将 $U$ 的第 $i$ 行存放在 $A$ 的第 $i$ 行 (上三角部分), 这样就不需要另外分配空间存储 $L$ 和 $U$ . 计算结束后, $A$ 的上三角部分为 $U$ , 其严格下三角部分为 $L$ 的绝对下三角部分. 此时算法可以描述为:

算法2.3.LU分解（用 $A$ 存储 $L$ 和 $U$ ）

1: for $k = 1$ to $n - 1$ do
2: for $i = k + 1$ to $n$ do
3: $a_{ik} = a_{ik} / a_{kk}$
4: for $j = k + 1$ to $n$ do
5: $a_{ij} = a_{ij} - a_{ik}a_{kj}$
6: end for
7: end for
8: end for

LU分解的运算量

由算法2.2可知，LU分解的运算量为

• 乘法次数: $T_{p} = \sum_{k=1}^{n-1} \sum_{i=k+1}^{n} \sum_{j=k+1}^{n} 1 = \sum_{k=1}^{n-1} (n-k)^{2} = \frac{1}{6} (2n^{3}-3n^{2}+n) = \frac{1}{3} n^{3}+\mathcal{O}(n^{2})$ .

• 加法次数: $T_{a} = \sum_{k=1}^{n-1} \sum_{i=k+1}^{n} \sum_{j=k+1}^{n} 1 = \frac{1}{3} n^{3} - \frac{1}{2} n^{2} + \frac{1}{6} n = \frac{1}{3} n^{3} + \mathcal{O}(n^{2})$ .

• 除法次数: $T_{d} = \sum_{k=1}^{n-1} \sum_{i=k+1}^{n} 1 = \frac{1}{2} n(n-1) = \frac{1}{2} n^{2} - \frac{1}{2} n = \frac{1}{2} n^{2} + \mathcal{O}(n)$ .

因此总运算量 (含乘法、除法与加法) 为 $\frac{2}{3} n^{3} - \frac{1}{2} n^{2} - \frac{1}{6} n = \frac{2}{3} n^{3} + \mathcal{O}(n^{2})$ .

数据更新顺序与计算效率

根据指标的循环次序, 算法 2.3 也称为 KIJ 型 LU 分解. 在实际计算中, 我们一般不建议使用这个算法. 因为对于指标 $k$ 的每次循环, 都需要更新 $A$ 的第 $k + 1$ 至第 $n$ 行. 这种反复读取数据的做法会使得计算效率大大降低.

如果数据是按行存储的, 如 $C / C++$ , 为了提高计算效率, 我们可以采用下面的IKJ型LU分解算法.

算法2.4.LU分解（IKJ型）

1: for $i = 2$ to $n$ do
2: for $k = 1$ to $i - 1$ do
3: $a_{ik} = a_{ik} / a_{kk}$
4: for $j = k + 1$ to $n$ do
5: $a_{ij} = a_{ij} - a_{ik}a_{kj}$
6: end for
7: end for
8: end for

IKJ型LU分解算法可以用下图来描述. ([113,p.291])

思考：如果数据是按列存储的，如FORTRAN或MATLAB，则怎样设计算法比较好？

2.1.3 Gauss消去法

假定矩阵 $A$ 存在LU分解 $A = LU$ ，则方程组 $Ax = b$ 就转化为求解下面两个三角方程组

这两个方程组都非常容易求解, 于是 Gauss 消去法描述如下:

算法2.5.Gauss消去法

1: 将 $A$ 进行 LU 分解: $A = LU$ , 其中 $L$ 为单位下三角矩阵, $U$ 为非奇异上三角矩阵;
2: 利用向前回代, 求解 $Ly = b$ , 即得 $y = L^{-1}b$ ;
3. 利用向后回代, 求解 $Ux = y$ , 即得 $x = U^{-1}y = (LU)^{-1}b = A^{-1}b$

得到 $A$ 的LU分解后, 我们最后需要用回代法求解两个三角方程组, 计算过程描述如下.

算法2.6.回代求解 $Ly = b$ 和 $Ux = y$

1: $y_{1} = b_{1} / l_{11}$ % 向前回代求解 $Ly = b$
2: for $i = 2:n$ do
3: for $j = 1 : i - 1$ do
4: $b_{i} = b_{i} - l_{ij}y_{j}$
5: end for
6: $y_{i} = b_{i} / l_{ii}$
7: end for
8: $x_{n} = y_{n} / u_{nn}$ $\%$ 向后回代求解 $U x = y$
9: for $i = n - 1 : -1 : 1$ do
10: for $j = n: -1: i + 1$ do
11: $y_{i} = y_{i} - u_{ij}x_{j}$
12: end for
13: $x_{i} = y_{i} / u_{ii}$
14: end for

如果数据是按列存储的, 则采用列存储方式效率会高一些. 下面是以 $Ux = y$ 为例, 描述按列存储时的回代求解过程.

算法2.7. 向后回代求解 $Ux = y$ (列存储方式)

1: for $k = n : -1 : 1$ do
2: $x_{k} = y_{k} / u_{kk}$
3: for $i = k - 1: -1:1$ do
4: $y_{i} = y_{i} - x_{k}u_{ik}$
5: end for
6: end for

计算复杂度

这两个算法的运算量: 乘法和加法均 $n(n - 1)$ 次, 除法 $2n$ 次. 加上 LU 分解的运算量, Gauss 消去法的总运算量为 $\frac{2}{3} n^3 + \mathcal{O}(n^2)$ .

可以证明, 以上两个算法都是向后稳定的 (componentwise backward stable) [70].

2.1.4 选主元LU分解

在 LU 分解过程中, 我们称 $a_{kk}^{(k-1)}$ 为主元. 如果 $a_{kk}^{(k-1)} = 0$ , 则算法就无法进行下去. 即使 $a_{kk}^{(k-1)}$ 不为零, 但如果 $|a_{kk}^{(k-1)}|$ 的值很小, 由于舍入误差的原因, 也可能会给计算结果带来很大的误差. 此时我们就需要通过选主元来解决这个问题.

例2.1 用LU分解求解线性方程组 $Ax = b$ ，其中 $A = \begin{bmatrix} 0.02 & 61.3 \\ 3.43 & -8.5 \end{bmatrix}, b = \begin{bmatrix} 61.5 \\ 25.8 \end{bmatrix}$ ，要求在运算过程中保留3位有效数字. (板书)

解．根据LU分解待定系数法，设

直接比较等式两边可得

解方程组 $Ly = b$ 可得

解方程组 $Ux = y$ 可得

□

事实上, 精确解为 $x_{1} = 10.0$ 和 $x_{2} = 1.00$ . 我们发现 $x_{1}$ 的数值解误差非常大. 出现这个问题的原因就是 $\left|a_{11}\right|$ 太小, 用它做主元时会放大舍入误差.

选主元是通过利用置换矩阵(也称排列矩阵)对行进行交换来实现的. 首先介绍置换矩阵的一些基本性质.

引理2.4设 $P\in \mathbb{R}^{n\times n}$ 为置换矩阵， $A\in \mathbb{R}^{n\times n}$ 为任意矩阵，则

(1) $PA$ 相当于将 $A$ 的行进行置换; $AP$ 相当于将 $A$ 的列进行置换;
(2) $P^{-1} = P^{\mathsf{T}}$ ，即 $P$ 是正交矩阵；
(3) $\operatorname{det}(P) = \pm 1$ ;
(4) 置换矩阵的乘积仍然是置换矩阵.

定理2.5 (选主元LU分解的存在性) 设 $A \in \mathbb{R}^{n \times n}$ 非奇异, 则存在置换矩阵 $P_L, P_R$ , 以及单位下三角矩阵 $L$ 和非奇异上三角矩阵 $U$ , 使得

(板书)

证明. 用归纳法

当 $n = 1$ 时, 取 $P_L = P_R = L = 1, U = A$ 即可.

假设结论对所有 $n - 1$ 阶矩阵都成立

考虑 $n$ 阶非奇异矩阵 $A \in \mathbb{R}^{n \times n}$ , 易知 $A$ 至少存在一个非零元, 取置换矩阵 $P_{1}$ 和 $P_{2}$ 使得

其中 $a_{11} \neq 0, A_{22} \in \mathbb{R}^{(n-1) \times (n-1)}$ . 令

则 $u_{11} \neq 0$ ，且有

两边取行列式可得

所以 $\operatorname{det}(U_{22}) \neq 0$ ，即 $U_{22} \in \mathbb{R}^{(n-1) \times (n-1)}$ 非奇异. 由归纳假设可知, 存在置换矩阵 $\tilde{P}_L \in \mathbb{R}^{(n-1) \times (n-1)}$ 和 $\tilde{P}_R \in \mathbb{R}^{(n-1) \times (n-1)}$ , 使得

其中 $\tilde{L}_{22}$ 为单位下三角矩阵, $\tilde{U}_{22}$ 为非奇异上三角矩阵. 令

则有

其中 $L$ 为单位下三角矩阵, $U$ 为非奇异上三角矩阵, 即 $A$ 存在 LU 分解.

因此, 由数学归纳法可知, 结论成立.

定理2.5也可以利用定理2.1的结论来证明，见习题2.3
事实上, 定理2.5中的 $P_L$ 和 $P_R$ 只有一个是必需的

置换矩阵的选取

问题: 第 $k$ 步时, 如何选取置换矩阵 $P_L^{(k)}$ 和 $P_R^{(k)}$ ?

选法一.选取 $P_L^{(k)}$ 和 $P_R^{(k)}$ 使得主元为剩下的矩阵中绝对值最大,这种选取方法称为“全主元Gauss消去法”,简称GECP(Gaussianeliminationwithcomplete pivoting);

选法二.选取 $P_L^{(k)}$ 和 $P_R^{(k)}$ 使得主元为第 $k$ 列中第 $k$ 到第 $n$ 个元素中,绝对值最大,这种选取方法称为“部分选主元Gauss消去法”,简称GEPP(Gaussianeliminationwithpartial pivoting),此时 $P_R^{(k)} = I$ ,因此也称为列主元Gauss消去法.

(1) GECP 比 GEPP 更稳定, 但工作量太大, 在实际应用中通常使用 GEPP 算法.
(2) GEPP 算法能保证 $L$ 所有的元素的绝对值都不超过 1.

算法2.8.部分选主元LU分解

1: $p = 1:n$ $\%$ 用于记录置换矩阵
2: for $k = 1$ to $n - 1$ do

3: $[a_{\max}, l] = \max_{k < i < n} |a_{ik}| \quad \%$ 选列主元, 其中 $l$ 表示主元所在的行
4: if $l \neq k$ then

5: for $j = 1$ to $n$ do
6: $a_{tmp} = a_{kj}, a_{kj} = a_{lj}, a_{lj} = a_{tmp} \quad \%$ 交换 $A$ 的第 $k$ 行与第 $l$ 行

7: end for

8: $p_{tmp} = p(k), p(k) = p(l), p(l) = p_{tmp} \quad \%$ 更新置换矩阵
9: end if

10: for $i = k + 1$ to $n$ do
11: $a_{ik} = a_{ik} / a_{kk}$
12: for $j = k + 1$ to $n$ do
13: $a_{ij} = a_{ij} - a_{ik}a_{kj}$
14: end for
15: end for
16: end for

相应的MATLAB程序见LE_PLU.m.

例2.2 用部分选主元LU分解求解线性方程组 $Ax = b$ ，其中 $A = \begin{bmatrix} 0.02 & 61.3 \\ 3.43 & -8.5 \end{bmatrix}$ ， $b = \begin{bmatrix} 61.5 \\ 25.8 \end{bmatrix}$ ，要

求在运算过程中保留3位有效数字

(板书)

解．由于 $|a_{21}| > |a_{11}|$ ，根据部分选主元LU分解算法，我们需要将第一行与第二行交换，即取 $P = \left[ \begin{array}{ll}0 & 1\\ 1 & 0 \end{array} \right],$ 然后计算 $\tilde{A} = PA = \left[ \begin{array}{ll}3.43 & -8.5\\ 0.02 & 61.3 \end{array} \right]$ 的LU分解，即设

直接比较等式两边可得

即

解方程组 $Ly = Pb$ 可得

解方程组 $Ux = y$ 可得

与精确解 $x_{1} = 10, x_{2} = 1$ 相比，数值解具有3位有效数字

2.1.5 矩阵求逆

我们可以通过部分选主元 LU 分解来计算矩阵的逆. 设 $PA = LU$ , 则

等价于求解下面 $2n$ 个三角线性方程组

思考：也可以分别计算 $L^{-1}$ 和 $U^{-1}$ ，然后相乘.哪种方法划算？

2.1.6 分块LU分解

为了提高计算效率, 实际计算中通常采用分块 LU 分解, 即

其中 $A_{ii}$ 是非奇异的方阵

与定理2.1相类似, 对于分块LU分解, 我们有下面的结论

定理2.6(分块LU分解的存在性和唯一性)矩阵 $\mathbf{A}$ 存在唯一分块LU分解的充要条件是 $\mathbf{A}$ 的所有顺序分块主子矩阵都非奇异. (留作练习)

几点注记

分块LU分解不是LU分解，若 $A$ 存在LU分解，则一定存在分块LU分解，反之不成立.
分块 LU 分解与普通 LU 分解的总运算量是一样的, 但分块 LU 分解可以借助 3 级 BLAS 运算, 因此效率更高.
在计算分块 LU 分解过程中需要求解一系列小规模的线性方程组 (以 $A_{ii}$ 为系数矩阵), 可以采用普通 (选主元) LU 分解方法求解.