分块矩阵的块不完全Cholesky分解算法
题目描述
本算法讲解块不完全Cholesky分解。这是一种用于大型稀疏对称正定线性方程组求解的预处理技术。其核心思想是：将矩阵按块结构划分，然后对矩阵进行一种不完全的Cholesky分解，即分解过程中只保留与原始矩阵稀疏结构相似的非零块（或忽略某些块内的元素），得到一个稀疏的近似下三角块矩阵及其转置的乘积。这个分解结果可作为预处理子，与迭代法（如共轭梯度法）结合，大幅加速求解过程。本专题将详解算法的推导、步骤、存储与计算策略。
循序渐进讲解
1. 问题背景与动机
- 考虑线性方程组：
\(A x = b\)
其中 \(A\) 是 \(n \times n\) 的大型稀疏对称正定矩阵，通常来自偏微分方程（PDE）的离散化（如有限元、有限差分）。 - 当 \(n\) 极大时，直接法（如完整Cholesky分解）会因填充（fill-in）而产生大量非零元，内存消耗过高。
- 迭代法（如共轭梯度法）适合稀疏矩阵，但对条件数较大的矩阵收敛缓慢。
- 预处理的目标是构造一个矩阵 \(M \approx A\)，使得 \(M^{-1} A\) 的条件数接近1，从而加速迭代收敛。
理想预处理子 \(M\) 应满足：- 逼近 \(A\) 的效果好；
- 构造和求解 \(M z = r\) （其中 \(r\) 为残差）的计算成本低。
块不完全Cholesky分解通过利用矩阵的块稀疏结构，在近似精度和计算效率之间取得平衡。
2. 块矩阵结构与符号约定
- 将矩阵 \(A\) 划分为 \(m \times m\) 的块，每个块大小为 \(p \times p\)（为简化，假设所有块大小相同，实际可不同）：
\[ A = \begin{pmatrix} A_{11} & A_{12} & \cdots & A_{1m} \\ A_{21} & A_{22} & \cdots & A_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ A_{m1} & A_{m2} & \cdots & A_{mm} \end{pmatrix} \]
其中 \(A_{ij}\) 是 \(p \times p\) 的子矩阵，且因 \(A\) 对称，有 \(A_{ij}^T = A_{ji}\)。
- 分解的目标是得到近似分解：
\(A \approx L_B L_B^T\)
其中 \(L_B\) 是块下三角矩阵：
\[ L_B = \begin{pmatrix} L_{11} & & & \\ L_{21} & L_{22} & & \\ \vdots & \vdots & \ddots & \\ L_{m1} & L_{m2} & \cdots & L_{mm} \end{pmatrix} \]
每个 \(L_{ii}\) 是下三角矩阵，\(L_{ij}\)（\(i > j\)）是一般块矩阵。
3. 块不完全分解的推导
从完整块Cholesky分解出发，步骤为：
- 初始：设 \(A^{(1)} = A\)。
- 对 \(k = 1\) 到 \(m\) 执行：
a. 计算当前主块 \(A_{kk}^{(k)}\) 的Cholesky分解：
\[ A_{kk}^{(k)} = L_{kk} L_{kk}^T \]
b. 更新同一列下方的块：对 \(i = k+1, \dots, m\)，
\[ L_{ik} = A_{ik}^{(k)} (L_{kk}^T)^{-1} \]
c. 对右下子矩阵进行更新：对 \(i, j = k+1, \dots, m\)，
\[ A_{ij}^{(k+1)} = A_{ij}^{(k)} - L_{ik} L_{jk}^T \]
关键的不完全性策略：
- 在步骤 (c) 中，只更新那些在原始矩阵 \(A\) 中非零的块 \(A_{ij}\)（即按块稀疏结构保留更新）。
- 对于原本为零的块位置，即使更新产生非零值，也将其丢弃（不存储），从而控制填充。
- 最终得到的 \(L_B\) 具有与 \(A\) 的下三角部分相同的块稀疏结构。
4. 算法步骤详解
输入：对称正定块矩阵 \(A\)，其块稀疏结构已知（通常用块行列索引表示）。
输出：块不完全Cholesky因子 \(L_B\)（块下三角）。
- 将 \(A\) 划分为 \(m \times m\) 块，记录非零块位置集合：
\(\mathcal{S} = \{(i,j) \mid A_{ij} \neq 0, \; i \ge j\}\) （下三角部分）。 - 初始化 \(L_B\) 为与 \(A\) 下三角部分相同块稀疏结构的零矩阵。
- 对于 \(k = 1\) 到 \(m\)：
a. 从 \(A\) 中提取 \(A_{kk}\)。
b. 对已计算的 \(L_{k,1:k-1}\) 块，按公式更新 \(A_{kk}\)：
\[ S_{kk} = A_{kk} - \sum_{j=1}^{k-1} L_{kj} L_{kj}^T \]
这里求和只对 $ (k,j) \in \mathcal{S} $ 的块进行。
c. 对 \(S_{kk}\) 进行完全Cholesky分解（因为它是稠密小矩阵）：
\[ S_{kk} = L_{kk} L_{kk}^T \]
存储 $ L_{kk} $。
d. 对于 \(i > k\) 且 \((i,k) \in \mathcal{S}\)：
- 提取 \(A_{ik}\)。
- 更新：
\[ T_{ik} = A_{ik} - \sum_{j=1}^{k-1} L_{ij} L_{kj}^T \]
求和只对 $ (i,j), (k,j) \in \mathcal{S} $ 的块进行。
- 计算：
\[ L_{ik} = T_{ik} (L_{kk}^T)^{-1} \]
这里 $ (L_{kk}^T)^{-1} $ 可通过解 $ L_{kk}^T X = T_{ik}^T $ 得到（利用前向/回代）。
- 存储 $ L_{ik} $。
- 完成循环，输出 \(L_B\)。
注意：第3.b和3.d中的求和，如果某些 \(L_{ij}\) 在稀疏结构中不存在（即被忽略），则相应项被跳过。这就是“不完全”的核心。
5. 预处理步骤
得到 \(L_B\) 后，预处理子为 \(M = L_B L_B^T\)。
在共轭梯度法（PCG）中，每一步需解 \(M z = r\)：
- 前代：解 \(L_B y = r\) （块下三角方程组）。
- 回代：解 \(L_B^T z = y\) （块上三角方程组）。
因为 \(L_B\) 是块稀疏的，这两步可通过块向前/向后替换高效完成，每个小块方程用直接法（如小规模Cholesky）求解。
6. 实际计算与存储优化
- 存储：只需存储 \(L_B\) 中对应于原始稀疏结构的非零块，可采用块压缩稀疏行（Block-CSR） 格式。
- 并行性：块分解中，对每个 \(k\)，步骤3.d中对不同 \(i\) 的计算是独立的，可并行。
- 数值稳定性：由于丢弃了一些更新，分解可能不总是保持正定性。实践中可采用对角线补偿（在对角块 \(S_{kk}\) 上增加小正数）或采用 \(L_B D_B L_B^T\) 形式（块不完全LDL^T）。
- 复杂度：约 \(O(\text{nnz}_b \cdot p^3)\)，其中 \(\text{nnz}_b\) 是非零块数，\(p\) 是块大小。当 \(p\) 小时很高效。
7. 示例
考虑一个 \(6 \times 6\) 矩阵（分3块，每块 \(2 \times 2\)），块稀疏结构如下（* 表示非零块）：
\[A = \begin{pmatrix} * & * & 0 \\ * & * & * \\ 0 & * & * \end{pmatrix} \]
（按块视图）。
分解步骤：
- \(k=1\)：分解 \(A_{11} = L_{11} L_{11}^T\)；更新 \(A_{21} \rightarrow L_{21}\)；更新 \(A_{22}\)（因为 \((2,2)\) 是非零块）。
- \(k=2\)：用已更新的 \(A_{22}\) 减去 \(L_{21} L_{21}^T\)，然后分解得到 \(L_{22}\)；更新 \(A_{32} \rightarrow L_{32}\)；更新 \(A_{33}\)。
- \(k=3\)：分解更新后的 \(A_{33}\) 得到 \(L_{33}\)。
最终 \(L_B\) 具有与 \(A\) 下三角相同的块稀疏模式。
总结
块不完全Cholesky分解结合了块结构的计算效率与不完全分解的内存节省，是求解大规模稀疏对称正定系统的有效预处理技术。其关键在于：
- 利用块稀疏性控制填充。
- 每个小块使用完全分解保证局部精度。
- 预处理系统可通过块三角求解快速计算。
此算法在计算流体力学、结构分析等领域有广泛应用，是许多迭代求解器（如PCG）的高效预处理子构造方法。