2.5 解的改进*

当矩阵 $A$ 是病态时, 即使残量 $r = b - A\hat{x}$ 很小, 所求得的数值解 $\hat{x}$ 仍可能带有较大的误差. 此时需要通过一些方法来提高解的精度.

2.5.1 高精度运算

在计算中, 尽可能采用高精度的运算. 比如, 原始数据是单精度的, 但在计算时都采用双精度运算, 或者更高精度的运算. 但更高精度的运算会带来更大的开销.

2.5.2 矩阵元素缩放或平衡 (Scaling or equilibration)

在求解线性方程组时, 先对系数矩阵元素进行适当的缩放是一种很常用的技术手段.

一方面, 如果 $A$ 的元素在数量级上相差很大, 则在计算过程中很可能会出现大数与小数的加减运算, 这样就可能会引入更多的舍入误差. 为了避免由于这种情况而导致的舍入误差, 我们可以在求解之前先对矩阵元素进行缩放 (Scaling), 即在矩阵两边同时乘以两个适当的对角矩阵. 在对矩阵进行缩放时, 需要 $\mathcal{O}(n^2)$ 运算量, 通常不会产生较大的舍入误差.

另一方面, 由前面的扰动分析可知, 矩阵条件数对数值计算的误差有着很大的影响, 是衡量问题是否病态的一个重要指标. 在实际计算中, 如果遇到病态问题, 我们希望能够通过一些简单的方法降低其条件数. 其中一类简单有效的方法就是在矩阵两边同时乘以一个对角矩阵, 即寻找对角矩阵 $D_{r}$ 和 $D_{c}$ , 使得 $D_{r}^{-1}AD_{c}^{-1}$ 的条件数达到最小 (或者尽可能小). 显然, 寻找这样的对角矩阵是非常困难的, 目前仍然是一个开放问题 [117].

一种比较可行的实用方案是使得缩放后的矩阵的每行或每列具有相同的 $p$ -范数．比如取 $D_r = \mathrm{diag}(d_1,d_2,\ldots ,d_n)$ ，其中 $d_{i} = \sum_{j = 1}^{n}|a_{ij}|$ ，这样 $D_r^{-1}A$ 的每行的1-范数都一样,但没法使得所有列也具有相同的范数.

如果矩阵

的所有行和所有列都具有相同 (或近似相同) 的范数, 则称为 $A$ 的均衡化 (equilibration) [117]. 有学者提出了一种迭代方法对 $A$ 进行均衡化 [82]. 更多信息可参见相关资料, 如 [38].

2.5.3 迭代改进法

设近似解 $\hat{x}$ , 残量 $r = b - A\hat{x}$ . 当 $\hat{x}$ 没达到精度要求时, 可以考虑方程组 $A z = r$ . 设 $z$ 是该方程组的精确解, 则

因此 $\hat{x} + z$ 就是原方程组的精确解

在实际计算中, 我们可能得到的是近似解 $\hat{z}$ , 但通常 $\|r - A\hat{z}\|$ 应该比较小, 特别地, 比 $\|r\|$ 更小. 因此 $\hat{x} + \hat{z}$ 应该比 $\hat{x}$ 更接近精确解.

如果新的近似解 $\hat{x} +\hat{z}$ 仍不满足精度要求, 则可重复以上过程

这就是通过迭代来提高解的精度.

算法2.15.通过迭代改进解的精度

1: 设 $PA = LU, \hat{x}$ 是 $Ax = b$ 的近似解
2: while 近似解 $\hat{x}$ 不满足精度要求, do
3: 计算 $r = b - A\hat{x}$
4: 求解 $Ly = Pr$ ，即 $y = L^{-1}Pr$
5: 求解 $Uz = y$ ，即 $z = U^{-1}y$
6: 令 $\hat{x} = \hat{x} + z$
7: end while

由于每次迭代只需计算一次残量和求解两个三角线性方程组, 因此运算量为 $\mathcal{O}(n^2)$ . 所以相对来讲还是比较经济的.

为了提高计算精度, 在计算残量 $r$ 时最好使用原始数据 $A$ , 而不是 $P^{\mathsf{T}}LU$ , 因此对 $A$ 做PLU分解时需要保留矩阵 $A$ , 不能被 $L$ 和 $U$ 覆盖.
由于在计算残量时, 可能会存在大量的相近的数相减运算, 因此可以采用更高精度的运算.
实际计算经验表明, 当 $A$ 病态不是很严重时, 即 $\varepsilon_{u} \kappa_{\infty}(A) < 1$ , 迭代法可以有效改进解的精度, 最后达到机器精度. 但 $\varepsilon_{u} \kappa_{\infty}(A) \geq 1$ 时, 一般没什么效果. 这里 $\varepsilon_{u}$ 表示机器精度.