计算数学中的预处理技术
预处理技术是数值线性代数中的核心概念，旨在加速大型稀疏线性方程组迭代求解的收敛速度。其基本思想是通过一个“预处理矩阵”来改造原系统，使新系统的系数矩阵具有更优的谱性质（例如条件数更小、特征值聚集），从而让迭代法（如共轭梯度法、GMRES法）能用更少的步骤达到收敛。
第一步：问题的提出——为什么需要预处理？
当使用迭代法（如共轭梯度法）求解线性方程组 \(A\mathbf{x} = \mathbf{b}\) 时，其收敛速度严重依赖于系数矩阵 \(A\) 的谱性质。具体来说：
- 条件数过大：如果矩阵 \(A\) 的条件数 \(\kappa(A)\) 很大（即矩阵是“病态”的），迭代法收敛会非常缓慢，甚至停滞不前。
- 特征值分布分散：即使条件数不大，但如果特征值广泛分散在复平面的不同区域，迭代法的收敛性也会很差。
预处理的核心目标就是通过一个线性变换，将原始的“难解”系统转化为一个等价的“易解”系统。
第二步：预处理的基本形式
预处理有两种等价的主要形式：
- 左预处理：我们寻找一个非奇异矩阵 \(M\)，使其在某种意义上近似于 \(A\)（即 \(M \approx A\)），然后用 \(M^{-1}\) 左乘原方程：
\[ M^{-1}A\mathbf{x} = M^{-1}\mathbf{b} \]
这个新系统的系数矩阵是 \(M^{-1}A\)。我们希望 \(M^{-1}A\) 的条件数远小于 \(A\) 的条件数，或者其特征值紧密地聚集在一起。
- 右预处理：我们也可以将预处理矩阵应用到变量上。设 \(\mathbf{x} = M^{-1}\mathbf{y}\)，代入原方程：
\[ A(M^{-1}\mathbf{y}) = \mathbf{b} \]
整理得：
\[ (AM^{-1})\mathbf{y} = \mathbf{b} \]
我们先求解关于 \(\mathbf{y}\) 的方程组，其系数矩阵为 \(AM^{-1}\)，然后再通过 \(\mathbf{x} = M^{-1}\mathbf{y}\) 得到原解。右预处理在算法实现上有时比左预处理更方便。
第三步：预处理矩阵 \(M\) 的构造原则
预处理矩阵 \(M\) 的选择是预处理技术的灵魂。一个“好”的预处理矩阵 \(M\) 需要满足两个看似矛盾的要求：
- 有效性：\(M^{-1}A\)（或 \(AM^{-1}\)）应该具有快速收敛的谱性质。
- 经济性：计算 \(M^{-1}\) 作用于任意向量的操作（即求解形如 \(M\mathbf{z} = \mathbf{r}\) 的方程组）必须计算代价低廉。如果求解 \(M\mathbf{z} = \mathbf{r}\) 和求解 \(A\mathbf{x} = \mathbf{b}\) 一样困难，那么预处理就失去了意义。
因此，预处理矩阵 \(M\) 通常是原矩阵 \(A\) 的一个“廉价且有效”的近似。
第四步：经典的预处理子示例
根据对“廉价近似”的不同理解，产生了多种经典的预处理子：
-
雅可比预处理（对角预处理）：这是最简单的预处理子。取 \(M\) 为 \(A\) 的对角线部分，即 \(M = \text{diag}(A)\)。其逆矩阵极易计算，对于强对角占优的矩阵效果显著。
-
高斯-赛德尔预处理/SOR预处理：取 \(M\) 为 \(A\) 的下三角部分（包括对角线），即 \(M = L + D\)。这对应于高斯-赛德尔迭代法的分裂矩阵。求解 \(M\mathbf{z} = \mathbf{r}\) 是一个前向代入过程，计算代价低。SOR（逐次超松弛）预处理是其推广。
-
不完全LU分解（ILU）预处理：这是最著名且应用最广泛的预处理子之一。其思想是模仿LU分解 \(A = LU\)，但在分解过程中引入一个“填充”策略，只允许在预定的稀疏模式位置（例如，只保留原矩阵 \(A\) 的非零元位置，即ILU(0)）上进行计算和存储。这样得到的 \(L\) 和 \(U\) 是稀疏的下三角和上三角矩阵，且 \(M = LU \approx A\)。求解 \(M\mathbf{z} = \mathbf{r}\) 只需进行稀疏矩阵的前向和后向代入，计算高效。还有更精确的ILU(k)和ILUT等变体。
第五步：预处理迭代法的算法实现
以最常用的预条件共轭梯度法（PCG）为例，展示预处理如何融入迭代法。算法步骤如下：
- 初始化 \(\mathbf{x}_0\)，计算残差 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\)。
- 求解预处理系统：\(\mathbf{z}_0 = M^{-1}\mathbf{r}_0\)。
- 设初始搜索方向 \(\mathbf{p}_0 = \mathbf{z}_0\)。
- 对于 \(k=0,1,2,\dots\)，直到收敛：
a. 计算步长：\(\alpha_k = (\mathbf{r}_k^T \mathbf{z}_k) / (\mathbf{p}_k^T A \mathbf{p}_k)\)。
b. 更新解：\(\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k\)。
c. 更新残差：\(\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k A \mathbf{p}_k\)。
d. 求解预处理系统：\(\mathbf{z}_{k+1} = M^{-1}\mathbf{r}_{k+1}\)。
e. 计算系数：\(\beta_k = (\mathbf{r}_{k+1}^T \mathbf{z}_{k+1}) / (\mathbf{r}_k^T \mathbf{z}_k)\)。
f. 更新搜索方向：\(\mathbf{p}_{k+1} = \mathbf{z}_{k+1} + \beta_k \mathbf{p}_k\)。
注意，算法中我们从不显式地计算 \(M^{-1}\) 或 \(M^{-1}A\)，而是在每一步需要计算 \(M^{-1}\) 作用于当前残差向量 \(\mathbf{r}_k\) 的结果 \(\mathbf{z}_k\)，这通过求解一个以 \(M\) 为系数矩阵的线性方程组实现。