2.2 特殊方程组的求解

如果系数矩阵具有一定的特殊结构或性质, 则可以充分利用这些特殊结构或性质来设计高效的数值算法. 本节考虑以下特殊方程组的求解:

• 对称正定矩阵

• 对称不定矩阵

• 三对角矩阵

• 带状矩阵

• Toeplitz 矩阵

2.2.1 对称正定线性方程组

考虑线性方程组

其中 $A\in \mathbb{R}^{n\times n}$ 对称正定的

定理 2.7 (Cholesky 分解) 设 $A \in \mathbb{R}^{n \times n}$ 对称正定, 则存在唯一的对角线元素为正的下三角矩阵 $L$ , 使得

该分解称为Cholesky分解.

(板书)

证明. 首先证明存在性, 用数学归纳法

当 $n = 1$ 时, 由 $A$ 的对称正定性可知 $a_{11} > 0$ . 取 $l_{11} = \sqrt{a_{11}}$ 即可.

假定结论对所有不超过 $n - 1$ 阶的对称正定矩阵都成立. 设 $A \in \mathbb{R}^{n \times n}$ 是 $n$ 阶对称正定, 则 $A$ 可分解为

其中 $\tilde{A}_{22} = A_{22} - A_{12}^{\mathrm{T}}A_{12} / a_{11}$ . 由于 $A$ 对称正定, 所以 $\left[ \begin{array}{cc}1 & 0\\ 0 & \tilde{A}_{22} \end{array} \right]$ 也对称正定, 故 $\tilde{A}_{22}$ 是 $n - 1$ 阶对称正定矩阵. 根据归纳假设, 存在唯一的对角线元素为正的下三角矩阵 $\tilde{L}$ , 使得 $\tilde{A}_{22} = \tilde{L}\tilde{L}^{\mathrm{T}}$ . 令

易知, $L$ 是对角线元素均为正的下三角矩阵, 且

由归纳法可知, 对任意对称正定实矩阵 $A$ , 都存在一个对角线元素为正的下三角矩阵 $L$ , 使得

唯一性可以采用反证法, 留做练习.

A 该定理也可以通过 LU 分解的存在唯一性来证明.

Cholesky分解的实现

我们利用待定系数法来计算对称正定矩阵的Cholesky分解.设 $A = LL^{\mathsf{T}}$ ，即

直接比较等式两边的元素可得

先计算 $l_{11} = \sqrt{a_{11}}$ ，然后计算 $L$ 第1列的其他元素：

依此类推, 可逐次计算出 $L$ 的第 2, 3, $\ldots, n$ 列. 具体算法描述如下.

算法2.9.Cholesky分解

1: for $j = 1$ to $n$ do
2: $l_{jj} = \left(a_{jj} - \sum_{k=1}^{j-1} l_{jk}^2\right)^{1/2} \quad \%$ 先计算对角线元素
3: for $i = j + 1$ to $n$ do
4: $l_{ij} = \frac{1}{l_{jj}}\left(a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk}\right) \quad \%$ 计算 $L$ 的第 $j$ 列
5: end for
6: end for

关于Cholesky算法的几点说明

• 与 LU 分解一样, 可以利用 $A$ 的下三角部分来存储 $L$ ;

• Cholesky 分解算法的运算量: 乘法次数 $\frac{1}{6} n^{3} + \frac{1}{6} n = \frac{1}{3} n^{3} + \mathcal{O}(n^{2})$ , 加法次数 $\frac{1}{6} n^{3} + \frac{1}{6} n = \frac{1}{3} n^{3} + \mathcal{O}(n^{2})$ , 除法次数 $\frac{1}{2} n^{2} - \frac{1}{2} n$ , 开方次数 $n$ .

• Cholesky 分解算法是稳定的（稳定性与全主元 Gauss 消去法相当），故不需要选主元。

• 利用Cholesky分解求解对称正定线性方程组的算法也称为平方根法.

例2.3 已知矩阵 $A = \begin{bmatrix} 4 & 2 & 8 & 0\\ 2 & 10 & 10 & 9\\ 8 & 10 & 21 & 6\\ 0 & 9 & 6 & 34 \end{bmatrix}$ ，计算 $A$ 的Cholesky分解. (板书)

解. 设 $A$ 的Cholesky分解为 $A = LL^{\top}$ , 即

比较等式两边的第一行可得

比较等式两边的第二行可得

比较等式两边的第三行可得

比较等式两边的第四行可得

所以 $A$ 的Cholesky分解为

LDL

为了避免开方运算, 我们可以将 $A$ 分解为: $A = L D L ^{\top}$ , 即

通过待定系数法可得

因此我们也可以依次计算出 $D$ 和 $L$ 的各个元素, 这就是 $\mathrm{LDL}^{\top}$ 分解.

基于 $\mathbf{LDL}^{\mathsf{T}}$ 分解求解对称正定线性方程组的算法称为改进的平方根法, 其优点是不需要计算平方根.

算法2.10.改进的平方根法

1: % 先计算 LDL+ 分解
2: for $j = 1$ to $n$ do

3: $d_{j} = a_{jj} - \sum_{k = 1}^{j - 1}l_{jk}^{2}d_{k}$
4: for $i = j + 1$ to $n$ do
5: $l_{ij} = (a_{ij} - \sum_{k=1}^{j-1} l_{ik} d_k l_{jk}) / d_j$
6: end for
7: end for
8: $\%$ 解方程组: $Ly = b$ 和 $DL^{\top}x = y$
9: $y_{1} = b_{1}$
10: for $i = 2$ to $n$ do
11: $y_{i} = b_{i} - \sum_{k = 1}^{i - 1}l_{ik}y_{k}$
12: end for
13: $x_{n} = y_{n} / d_{n}$
14: for $i = n - 1$ to 1 do
15: $x_{i} = y_{i} / d_{i} - \sum_{k = i + 1}^{n}l_{ki}x_{k}$
16: end for

2.2.2 对称不定线性方程组

设 $A \in \mathbb{R}^{n \times n}$ 是非奇异的对称不定矩阵. 若 $A$ 存在 LU 分解, 即 $A = LU$ , 则可写成

其中 $D$ 是由 $U$ 的对角线元素构成的对角矩阵. 然而, 当 $A$ 不定时, 其 LU 分解不一定存在. 若采用选主元 LU 分解, 则其对称性将被破坏. 为了保持对称性, 在选主元时必须对行和列进行同样的置换, 即选取置换矩阵 $P$ , 使得

通常称 (2.4) 为对称矩阵的 $\mathbf{LDL}^{\mathsf{T}}$ 分解。但遗憾的是，这样的置换矩阵可能不一定存在，即分解 (2.4) 不一定存在。

例2.4 设对称矩阵

由于 $A$ 的对角线元素都是0, 对任意置换矩阵 $P$ , 矩阵 $PAP^{\mathsf{T}}$ 的对角线元素仍然都是0. 因此, 矩阵 $A$ 不存在 $\mathrm{LDL}^{\mathsf{T}}$ 分解 (2.4).

块 $\mathbf{LDL}^{\top}$ 分解

由于对称不定矩阵不一定存在 $\mathbf{LDL}^{\mathsf{T}}$ 分解, 于是人们提出了块 $\mathbf{LDL}^{\mathsf{T}}$ 分解, 不再要求 $D$ 是对角矩阵, 它可以是拟对角矩阵 (即块对角矩阵, 而且对角块至多是 2 阶的).

人们发现, 对于一个对称非奇异矩阵 $A \in \mathbb{R}^{n \times n}$ , 总存在置换矩阵 $P$ 使得 (见习题 2.9)

其中 $B \in \mathbb{R}$ 或 $B \in \mathbb{R}^{2 \times 2}$ , 且非奇异. 因此可以对 $PAP^{\top}$ 进行块对角化, 即

其中 $C - EB^{-1}E^{\top}$ 是Schur补

不断重复以上过程, 就可以得到 $A$ 的块 $\mathbf{LDL}^{\top}$ 分解:

其中 $\tilde{D}$ 是拟对角矩阵

与选主元 LU 分解类似, 我们需要考虑块 $\mathrm{LDL}^{\top}$ 分解的选主元策略, 即如何选取置换矩阵 $P$ . 一方面使得分解可以进行下去, 另一方面提高算法的稳定性. Kahan 于 1965 年首先考虑了选主元块 $\mathrm{LDL}^{\top}$ 分解. 目前常用的策略有

• 全主元策略 (Bunch-Parlett 策略): Bunch 和 Parlett [22] 于 1971 年提出了全主元策略来选取置换矩阵, 并证明了其稳定性, 指出其向后误差上界与全主元 Gauss 消去法几乎一样 [20]. 但需要进行 $n^{3} / 6$ 次比较运算, 代价比较昂贵.

• 部分选主元策略 (Bunch-Kaufman 策略): 为了减少比较运算, Bunch 和 Kaufman [21] 于 1977 年提出了部分选主元策略. 这种策略在每次选主元时, 只需搜索两列, 因此将比较运算复杂度降低到了 $\mathcal{O}(n^{2})$ 量级. 而且这种选主元策略具有较满意的向后稳定性 [69, 70], 运算量为

$\frac{1}{3} n^{3} + \mathcal{O}(n^{2})$ ，因此被广泛使用.著名的程序库LAPACK一开始就采用了这种策略.具体实施也可以参见[70],或31.

• Rook 策略: 采用 Bunch-Kaufman 策略得到的矩阵 $L$ 的元素不能得到很好的控制. 为了克服由此带来的不稳定性, Ashcraft, Grimes 和 Lewis [4] 于 1998 年提出采用对称的 Rook 选主元策略,以潜在的高比较运算量, 可以将矩阵 $L$ 的元素控制在 1 左右. 这种对称的 Rook 选主元策略最先用在 Gauss 消去法中 [96]. Rook 策略整体上与 Bunch-Kaufman 策略类似, 但在选主元时加了一层迭代, 从而能够提供更高的精度. Cheng 在其博士论文 [25] 中对这两种策略进行了数值测试, 给出了两种例子, 分别说明 Rook 策略和 Bunch-Kaufman 策略各有优势. LAPACK 从 3.5.0 版本开始增加了 Rook 策略.

目前大部分软件都采用部分选主元块 $\mathbf{LDL}^{\top}$ 分解来求解对称线性方程组. 关于 Aasen 算法和块 $\mathbf{LDL}^{\top}$ 分解比较也可参见 [4, 12].
以上算法是针对稠密的对称不定线性方程组. 对于稀疏的对称不定矩阵, 可以采用Duff和Reid的多波前法(Multifrontal) [39, 40], 或者Liu提出的稀疏阈值算法[93].

Aasen算法

Aasen[1]于1971年提出了下面的分解

其中 $P$ 为置换矩阵, $L$ 为单位下三角矩阵, $T$ 为对称三对角矩阵. 分解 (2.5) 本质上与部分选主元 LU 分解是一样的, 具体实施细节可参见 [70, 150].

2011年,Rozloznik,Shklarski和Toledo[107]给出了基于现代化计算机结构和BLAS-3的分块三对角化算法.该算法本质上是Aasen算法的分块形式.

Aasen算法不如的块 $\mathsf{LDL}^{\top}$ 分解使用广泛.2016年12月，LAPACK发行了更新版本3.7.0,由Yamazaki加入了Aasen算法(http://www.netlib.org/lapack/lapack-3.7.0.html).

2.2.3 三对角线性方程组

考虑三对角线性方程组 $Ax = f$ ，其中 $A$ 是三对角矩阵：

我们假定

且

即 $A$ 是不可约弱对角占优的. 此时, 我们可以得到下面的三角分解

由待定系数法, 我们可以得到递推公式:

为了使得算法能够顺利进行下去，我们需要证明 $\alpha_{i}\neq 0$

定理2.8 设三对角矩阵 $A$ 满足条件(2.6)和(2.7).则 $A$ 非奇异，且

(1) $|\alpha_{1}| = |b_{1}| > 0$
(2) $0 < |\beta_{i}| < 1, i = 1,2,\ldots ,n - 1;$
(3) $0 < |c_{i}|\leq |b_{i}| - |a_{i - 1}| < |\alpha_{i}| < |b_{i}| + |a_{i - 1}|,\quad i = 2,3,\ldots ,n.$ （板书）

证明. 由于 $A$ 是不可约且弱对角占优, 所以 $A$ 非奇异. (见定理 1.61)

结论 (1) 是显然的.

下面我们证明结论 (2) 和 (3).

由于 $0 < |c_{1}| < |b_{1}|$ , 且 $\beta_{1} = c_{1} / b_{1}$ , 所以 $0 < |\beta_{1}| < 1$ . 又 $\alpha_{2} = b_{2} - a_{1}\beta_{1}$ , 所以

再由结论 $(\ref{eq:2})$ 和 $\beta_{2}$ 的计算公式可知 $0 < |\beta_2| < 1$ . 类似于 $(\ref{eq:3})$ 和 $(\ref{eq:4})$ , 我们可以得到

依此类推, 我们就可以证明结论 (2) 和 (3).

由定理2.8可知,分解(2.8)是存在的.因此，原方程就转化为求解 $Ly = f$ 和 $Ux = y.$ 由此便可得求解三对角线性方程组的追赶法也称为Thomas算法(1949)，其运算量大约为 $8n - 6$

算法2.11.追赶法

2: $\beta_{1} = c_{1} / b_{1}$
3: $y_{1} = f_{1} / b_{1}$
4: for $i = 2$ to $n - 1$ do
5: $\alpha_{i} = b_{i} - a_{i - 1}\beta_{i - 1}$
6: $\beta_{i} = c_{i} / \alpha_{i}$
7: $y_{i} = (f_{i} - a_{i - 1}y_{i - 1}) / \alpha_{i}$
8: end for
9: $\alpha_{n} = b_{n} - a_{n - 1}\beta_{n - 1}$
10: $y_{n} = (f_{n} - a_{n - 1}y_{n - 1}) / \alpha_{n}$
11: $x_{n} = y_{n}$
12: for $i = n - 1$ to 1 do
13: $x_{i} = y_{i} - \beta_{i}x_{i + 1}$
14: end for

具体计算时, 由于求解 $Ly = f$ 与矩阵 LU 分解是同时进行的, 因此, $\alpha_{i}$ 可以不用存储.
由于 $|\beta_i| < 1$ ，因此在回代求解 $x_{i}$ 时，误差可以得到有效控制

需要指出的是, 我们也可以考虑下面的分解

但此时 $|\gamma_i|$ 可能大于1.比如 $\gamma_1 = a_1 / b_1$ ，因此当 $|b_{1}| < |a_{1}|$ 时， $|\gamma_1| > 1.$ 所以在回代求解时，误差可能得不到有效控制.另外一方面，计算 $\gamma_{i}$ 时也可能会产生较大的舍入误差(大数除以小数).但如果 $A$ 是列对角占优，则可以保证 $|\gamma_i| < 1$

如果 $A$ 是 (行) 对角占优, 则采用分解 (2.8); 如果 $A$ 是列对角占优, 则采用分解 (2.9).

2.2.4 带状线性方程组

设系数矩阵 $A \in \mathbb{R}^{n \times n}$ 是带状矩阵, 其下带宽为 $b_{l}$ , 上带宽为 $b_{u}$ , 即

当 $i > j + b_{l}$ 或 $i < j - b_u$ 时，有 $a_{ij} = 0$

其形状如右图所示

对于带状矩阵, 不带选主元的 LU 分解算法如下.

算法2.12.带状矩阵的LU分解

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

如果 $A$ 存在LU分解, 则可以验证, $L$ 和 $U$ 也是带状矩阵.

定理2.9 设 $A \in \mathbb{R}^{n \times n}$ 是带状矩阵, 其下带宽为 $b_{l}$ , 上带宽为 $b_{u}$ . 若 $A$ 存在不带选主元的 LU 分解 $A = LU$ , 则 $L$ 为下带宽 $b_{l}$ 的带状矩阵, $U$ 为上带宽 $b_{u}$ 的带状矩阵. (留作课外自习)

若采用部分选主元的 LU 分解, 则 $L$ 和 $U$ 会有所变化, 但仍具有一定的特殊结构.

定理2.10 设 $A \in \mathbb{R}^{n \times n}$ 是带状矩阵, 其下带宽为 $b_{l}$ , 上带宽为 $b_{u}$ . 设 $PA = LU$ 是 $A$ 的部分选主元LU分解, 则 $U$ 为上带宽不超过 $b_{l} + b_{u}$ 的带状矩阵, $L$ 为下带宽为 $b_{l}$ 的“基本带状矩阵”, 即 $L$ 每列的非零元素不超过 $b_{l} + 1$ 个. (留作课外自习)

思考：分别统计带状矩阵 LU 分解和部分选主元 LU 分解的运算量.

2.2.5 Toeplitz线性方程组

设 $T_{n}\in \mathbb{R}^{n\times n}$ 是Toeplitz矩阵，即

可知 Toeplitz 矩阵是反向对称 (persymmetric) 矩阵, 即关于东北-西南对角线对称. 记 $J_{n}$ 为 $n$ 阶反向单位矩阵, 即

易知 $J_{n}^{\top} = J_{n}^{-1} = J_{n},J_{n}^{2} = I_{n}$

引理2.11 矩阵 $A \in \mathbb{R}^{n \times n}$ 是反向对称矩阵当且仅当

若 $A$ 可逆，则可得

即反向对称矩阵的逆也是反向对称矩阵

Toeplitz 矩阵的逆是反向对称矩阵, 但不一定是 Toeplitz 矩阵.

一般情况下，Toeplitz 矩阵的乘积不再是一个 Toeplitz 矩阵，但如果是两个上三角或下三角 Toeplitz 矩阵相乘，则乘积仍然是 Toeplitz 矩阵。

定理2.12 设 $T, S \in \mathbb{R}^{n \times n}$ 是上三角 Toeplitz 矩阵, 则 $TS$ 也是上三角 Toeplitz 矩阵. 进一步, 若 $T$ 非奇异, 则 $T^{-1}$ 也是上三角 Toeplitz 矩阵.

Yule-Walker方程组

我们先考虑一类特殊右端项的线性方程组. 设 $T_{n}$ 对称正定, 考虑线性方程组

其中 $t^{(n)} = [t_1, t_2, \ldots, t_{n-1}, t_n]^\top$ , 即右端项的前 $n-1$ 个分量为 $T_n$ 第一列的后 $n-1$ 个元素, $t_n$ 任意. 这类线性方程组称为 Yule-Walker 方程组.

由于 $T_{n}$ 对称正定, 所以 $t_0 > 0$ . 因此我们可以对 $T_{n}$ 的对角线元素进行单位化处理. 不失一般性, 我们假定 $T_{n}$ 的对角线元素为1, 即

由于方程组右端项的特殊性, 我们可以不进行矩阵分解, 而是通过递推方法来求解.

记 $T_{k}$ 为 $T_{n}$ 的 $k$ 阶顺序主子矩阵, $t^{(k)}$ 为 $t^{(n)}$ 的前 $k$ 个分量组成的向量. 设 $y^{(k)}$ 是 $T_{k}y = -t^{(k)}$ 的解, 下面通过递推方法来计算 $T_{k + 1}y = -t^{(k + 1)}$ 的解 $y^{(k + 1)}$ . 记

则 $T_{k + 1}y^{(k + 1)} = -t^{(k + 1)}$ 可写为

即

由于 $T_{k}$ 既是反向对称矩阵, 也是对称矩阵, 因此有 $T_{k}^{-1}J_{k} = J_{k}T_{k}^{-1}$ . 所以根据 (2.11) 可得

代入 (2.12) 可得

又

由 $T_{k + 1}$ 的对称正定性可知 $1 + \left(t^{(k)}\right)^{\top}y^{(k)} > 0.$ 于是可得 $y^{(k + 1)}$ 的计算公式

运算量为 $\mathcal{O}(k)$ . 因此, 我们就可以从一阶 Yule-Walker 方程出发, 利用递推公式 (2.13) 逐步计算出 $T_{n}y = -r_{n}$ 的解. 总的运算量大约为 $3n^{2}$ .

为了进一步减少运算量, 我们引入一个辅助变量. 记 $\tau_k \triangleq 1 + (t^{(k)})^\top y^{(k)}$ , 则

于是可得求解 Yule-Walker 方程组的 Levinson-Durbin 算法 [37, 89] (由 Levinson 和 Durbin 分别于 1946 年和 1960 年独立提出 [81]), 总运算量大约为 $2n^{2}$ .

算法2.13.求解Yule-Walker方程组的Levinson-Durbin算法

1: 输入数据: $t = \left\lbrack {{t}_{1},{t}_{2},\ldots ,{t}_{n}}\right\rbrack \%$ 注: 这里假定 ${t}_{0} = 1$
2: $y(1) = -t_1, \tau = 1, \alpha = -t_1$
3: for $k = 1$ to $n - 1$ do
4: $\tau = (1 - \alpha^2)\tau$
5: $\alpha = -\frac{1}{\tau}\left(t_{k + 1} + \sum_{i = 1}^{k}t_{i}y(k + 1 - i)\right)$
6: $y(1:k) = y(1:k) + \alpha y(k: -1:1)$
7: $y(k + 1) = \alpha$
8: end for

由(2.13)可知，只要 $T_{k + 1}$ 非奇异，就能确保 $1 + \left(t^{(k)}\right)^{\top}y^{(k)}\neq 0$ ，算法就能顺利进行下去所以Levinson-Durbin算法有意义的充要条件是 $T_{n}$ 的所有顺序主子式非零，此时我们称 $T_{n}$ 是强正则(strongly regular)的.