分块矩阵的预条件共轭残差法（Block Preconditioned Conjugate Residual, Block PCR）求解对称正定多重右端项线性方程组

题目描述

考虑对称正定线性方程组 \(AX = B\)，其中 \(A \in \mathbb{R}^{n \times n}\) 是一个大型、稀疏的对称正定矩阵，而 \(B \in \mathbb{R}^{n \times s}\) 是一个具有 \(s\) 个右端项的矩阵（通常 \(s \ll n\)，例如 \(s\) 为几到几十）。我们的目标是同时求解这 \(s\) 个线性方程组：\(A\mathbf{x}_j = \mathbf{b}_j\)，其中 \(\mathbf{x}_j\) 和 \(\mathbf{b}_j\) 分别是矩阵 \(X\) 和 \(B\) 的第 \(j\) 列。

直接对每个右端项单独应用迭代法（如共轭梯度法，CG）需要 \(s\) 次独立的求解过程。然而，分块算法（如分块共轭梯度法，Block CG）可以利用多个右端项之间的信息，在一次迭代过程中同时更新所有解向量，通常能获得比单独求解更快的收敛速度（以迭代次数或总计算时间衡量）。当系数矩阵 \(A\) 的条件数较大或病态时，引入预条件子（Preconditioner）\(M \approx A^{-1}\) 是加速收敛的关键。预条件共轭残差法（PCR）是 CG 算法的一个变体，它要求残差在每一步都正交（CG 要求搜索方向共轭），在有限精度算术中有时比 CG 更鲁棒。因此，将 PCR 推广到分块情形并加入预条件处理，就形成了分块预条件共轭残差法（Block PCR），用于高效求解对称正定多重右端项问题。

解题过程详解

我们将分步推导并解释分块预条件共轭残差法的算法流程、数学原理和关键细节。

步骤 1：问题形式化与算法目标
我们的目标是求解：

\[A X = B, \quad A = A^T \succ 0, \quad X, B \in \mathbb{R}^{n \times s} \]

分块算法将同时处理所有 \(s\) 个右端项。算法会生成一系列近似解矩阵 \(X^{(k)}\)，使得残差矩阵 \(R^{(k)} = B - A X^{(k)}\) 的范数（例如 Frobenius 范数）逐渐减小至容忍误差以下。

步骤 2：单右端项预条件共轭残差法（PCR）回顾
为了理解分块版本，先回顾标准 PCR（一个右端项，\(s=1\)）。给定初始猜测 \(x^{(0)}\)，计算初始残差 \(r^{(0)} = b - A x^{(0)}\)。预条件步骤：求解 \(M z^{(0)} = r^{(0)}\) 得到预条件残差 \(z^{(0)}\)，并设初始搜索方向 \(p^{(0)} = z^{(0)}\)。
对于 \(k = 0, 1, 2, \dots\)，执行：

1. \(\alpha_k = (r^{(k)}, z^{(k)}) / (A p^{(k)}, p^{(k)})\) （计算步长）
2. \(x^{(k+1)} = x^{(k)} + \alpha_k p^{(k)}\) （更新解）
3. \(r^{(k+1)} = r^{(k)} - \alpha_k A p^{(k)}\) （更新残差）
4. 求解 \(M z^{(k+1)} = r^{(k+1)}\) （预条件步骤）
5. \(\beta_k = (r^{(k+1)}, z^{(k+1)}) / (r^{(k)}, z^{(k)})\) （计算正交化系数）
6. \(p^{(k+1)} = z^{(k+1)} + \beta_k p^{(k)}\) （更新搜索方向）
   其中 \((\cdot, \cdot)\) 表示内积。PCR 的关键性质是残差向量 \(\{r^{(k)}\}\) 相互正交：\((r^{(i)}, r^{(j)}) = 0, i \ne j\)。

步骤 3：扩展到分块情形——思想
分块算法的核心思想是同时对 \(s\) 个向量进行运算，将它们组织成矩阵。我们维持解矩阵 \(X^{(k)}\)、残差矩阵 \(R^{(k)}\)、预条件残差矩阵 \(Z^{(k)}\) 和搜索方向矩阵 \(P^{(k)}\)，每个都是 \(n \times s\) 维。目标是使得每个残差向量相互正交，并且算法能高效利用矩阵-矩阵乘法（Level-3 BLAS 操作），这通常比多次矩阵-向量乘更高效。

步骤 4：分块预条件共轭残差法（Block PCR）算法推导
给定初始猜测 \(X^{(0)}\)，计算初始残差 \(R^{(0)} = B - A X^{(0)}\)。预条件步骤：求解 \(M Z^{(0)} = R^{(0)}\)（这通常是近似求解 \(M \approx A^{-1}\)，例如通过不完全 Cholesky 分解 \(M = \tilde{L}\tilde{L}^T\) 的前代回代过程）。设初始搜索方向 \(P^{(0)} = Z^{(0)}\)。

对于迭代 \(k = 0, 1, 2, \dots\)：

1. 计算步长矩阵：步长不再是一个标量，而是一个 \(s \times s\) 矩阵 \(\alpha_k\)，由以下最小化问题确定：我们希望最小化新残差的某种范数。在分块 PCR 中，我们要求新的残差矩阵 \(R^{(k+1)}\) 与当前的 Krylov 子空间正交。经过推导（类似于分块 CG），步长矩阵 \(\alpha_k\) 满足：

\[ \alpha_k = (P^{(k)T} A P^{(k)})^{-1} (R^{(k)T} Z^{(k)}) \]

这里，\(P^{(k)T} A P^{(k)}\) 和 \(R^{(k)T} Z^{(k)}\) 都是 \(s \times s\) 矩阵。注意 \(R^{(k)T} Z^{(k)}\) 是对称正定的，因为 \(M\) 对称正定且 \(Z^{(k)} = M^{-1} R^{(k)}\)，所以 \(R^{(k)T} Z^{(k)} = R^{(k)T} M^{-1} R^{(k)} \succ 0\)。

2. 更新解和残差：

\[ X^{(k+1)} = X^{(k)} + P^{(k)} \alpha_k \]

\[ R^{(k+1)} = R^{(k)} - A P^{(k)} \alpha_k \]

注意这里 \(\alpha_k\) 是矩阵，乘法是矩阵乘法，\(P^{(k)} \alpha_k\) 是 \(n \times s\) 矩阵。

3. 预条件步骤：求解 \(M Z^{(k+1)} = R^{(k+1)}\) 得到预条件残差。

4. 计算正交化系数矩阵：为了保持搜索方向矩阵的“\(A\)-共轭性”和残差的正交性，我们需要一个 \(s \times s\) 矩阵 \(\beta_k\)。在分块 PCR 中，我们要求新的搜索方向矩阵 \(P^{(k+1)}\) 满足：

\[ R^{(k+1)T} Z^{(k+1)} \text{ 是对角矩阵（或至少 } R^{(k+1)T} Z^{(j)} = 0 \text{ for } j \le k \text{）} \]

通过构造，我们令：

\[ \beta_k = (R^{(k)T} Z^{(k)})^{-1} (R^{(k+1)T} Z^{(k+1)}) \]

由于 \(R^{(k)T} Z^{(k)}\) 对称正定，其逆存在。

5. 更新搜索方向矩阵：

\[ P^{(k+1)} = Z^{(k+1)} + P^{(k)} \beta_k \]

步骤 5：算法伪代码

输入：对称正定矩阵 A (n×n), 右端项矩阵 B (n×s), 预条件子 M, 初始猜测 X0 (n×s), 容忍误差 tol, 最大迭代次数 maxiter.
输出：近似解矩阵 X.

设置 X = X0.
计算 R = B - A X.
求解 M Z = R 得到 Z.
设置 P = Z.
计算 G = R^T Z.  // G 是 s×s 矩阵

for k = 0 to maxiter-1 do
    计算 AP = A * P.  // 主要的矩阵-矩阵乘，一次计算涉及 s 个矩阵-向量乘
    计算 H = P^T * AP. // H 是 s×s
    求解小规模线性系统 H * α = G 得到 α. // α 是 s×s
    更新 X = X + P * α.
    更新 R = R - AP * α.
    如果 ||R||_F < tol，则终止循环。
    求解 M Z_new = R 得到新的预条件残差。
    计算 G_new = R^T * Z_new.
    求解 G * β = G_new 得到 β. // 注意这里用前一个 G 的逆左乘 G_new
    更新 P = Z_new + P * β.
    更新 Z = Z_new, G = G_new.
end for
返回 X.

注意：实际实现中，求解小系统 \(H \alpha = G\) 和 \(G \beta = G_{\text{new}}\) 通常通过 Cholesky 分解完成，因为 \(H\) 和 \(G\) 都是对称正定的。

步骤 6：关键性质与讨论

1. 正交性：在精确算术下，算法保证残差矩阵满足 \(R^{(k)T} R^{(j)} = 0\) 对于 \(k \ne j\)（对于预条件版本是 \(R^{(k)T} Z^{(j)} = 0\)）。这源自 PCR 的性质推广。
2. 收敛性：分块 PCR 的收敛速度通常受 \(A\) 的极端特征值影响。预条件子的引入改善了条件数，从而加速收敛。
3. 计算优势：主要的计算成本在于矩阵-矩阵乘积 \(A P\)（每次迭代一次）和预条件子求解 \(M Z = R\)。前者可以利用 Level-3 BLAS 优化，后者需要高效求解 \(s\) 个线性系统（如果 \(M\) 是显式分解，如不完全 Cholesky，则可通过前代回代同时处理多右端项）。
4. 数值稳定性：当 \(s\) 较大时，搜索方向矩阵 \(P^{(k)}\) 的列可能变得线性相关（称为“列退化”），导致算法崩溃。实践中，可采用“选择性重新正交化”或“显式块大小缩减”策略来维持数值稳定性。
5. 预条件子选择：常见的预条件子包括不完全 Cholesky（IC）、代数多重网格（AMG）、稀疏近似逆（SPAI）等，选择取决于矩阵 \(A\) 的性质。

总结
分块预条件共轭残差法通过同时处理多个右端项，利用矩阵运算提升计算效率，并结合预条件技术加速收敛，是求解对称正定多重右端项大规模稀疏线性方程组的一种高效迭代算法。理解其推导需要掌握共轭残差法的基本思想和分块算法的推广技巧。