超级会员免费看
快速线性求解器：预条件共轭梯度法（PCGM）详解
1. 共轭梯度法（CGM）收敛性
在非对称系统求解器中，虽然某些情况有所不同，但 $A^{-1}$ 残差范数 $(r^TA^{-1}r)$ 会单调递减。由于有限算术运算，CGM 的收敛速率由矩阵 $A$ 的条件数 $\kappa_2(A)$ 控制，具体估计如下：
$\frac{|s - x_k|
2}{|s - x_0|_2} \leq 2\gamma^k(\kappa_2(A))$
其中 $\kappa_2(A) = \frac{\lambda
{max}}{\lambda_{min}}$，$\gamma = \frac{\sqrt{\kappa_2(A)} - 1}{\sqrt{\kappa_2(A)} + 1}$。
当 $\kappa_2$ 值较大时，$\gamma \to 1$，CGM 收敛所需的迭代次数与 $\sqrt{\kappa_2(A)}$ 成正比。例如，对于在 $N \times N$ 网格上使用二阶有限差分离散化的泊松方程，$\kappa_2(A) \propto N^2$，因此迭代次数与 $N$ 成正比。在假设 SOR 具有最优松弛参数的情况下，CGM 的计算成本与 SOR 相似。
2. 预条件器
为了加速 CGM 的收敛，我们使用预条件器。将线性系统 $Ax = b$ 乘以非奇异预条件矩阵 $M$ 转换为 $M^{-1}Ax = M^{-1}b$。预条件器应满足以下性质：
- 谱上接近矩阵 $A$，使得 $\kappa_2(M^{-1}A) \ll \kappa_2(A)$。
- 求逆成本低，因为需要求解 $Mx = b$。
- 由于 $A$ 是对称正定矩阵，$M$ 也必须是对称正定的。
为了“修正”搜索方向且不显著增加计算复杂度，我们需要对预条件系统进行对称化。将 $M$ 用其特征向量和特征值表示为 $M = V\Lambda V^T$，则 $M^{1/2} \equiv V\Lambda^{1/2}V^T$。
将 $M^{-1}Ax = M^{-1}b$ 乘以 $M^{1/2}$ 得到 $(M^{-1/2}AM^{-1/2})(M^{1/2}x) = M^{-1/2}b$，即 $By = f$，其中 $B \equiv M^{-1/2}AM^{-1/2}$，$y \equiv M^{1/2}x$，$f \equiv M^{-1/2}b$。$B$ 和 $M^{-1}A$ 相似，具有相同的特征值。
2.1 预条件共轭梯度算法（PCG）
以下是 PCG 算法的步骤：
1.
初始化
：
- 选择 $x_0$，计算 $r_0 = b - Ax_0$。
- 求解 $M\tilde{r}
0 = r_0$，得到 $p_0 = \tilde{r}_0$。
2.
循环
：对于 $k = 1, \cdots$
- 计算 $\alpha_k = \frac{(\tilde{r}_k, r_k)}{(p_k, Ap_k)}$。
- 更新 $x
{k + 1} = x_k + \alpha_kp_k$。
- 计算 $r_{k + 1} = r_k - \alpha_kAp_k$。
- 求解 $M\tilde{r}
{k + 1} = r
{k + 1}$。
- 如果 $(\tilde{r}
{k + 1}, r
{k + 1}) \leq \epsilon$ 且 $(r_{k + 1}, r_{k + 1}) \leq \epsilon$，则返回结果。
- 计算 $\beta_k = \frac{(\tilde{r}
{k + 1}, r
{k + 1})}{(\tilde{r}
k, r_k)}$。
- 更新 $p
{k + 1} = \tilde{r}_{k + 1} + \beta_kp_k$。
注意 ：停止准则基于实际残差 $r$，但先检查修改后的残差 $\tilde{r}$ 可以节省计算时间。
2.2 不同类型的预条件器
- 对角缩放（Jacobi 预条件） ：对于对角元素大小差异很大的矩阵，使用 $M = diag(a_{11}, a_{22}, \cdots, a_{nn})$ 非常有效，它可以将 $B$ 的条件数降低到最小值的 $n$ 倍。
- 块对角预条件器 ：可以将 $M$ 构建为 $A$ 的块子矩阵组成的块对角矩阵。对称 SOR 可以作为更有效的块预条件器。
-
不完全 Cholesky 分解
：对于有限差分离散化扩散问题或 Toeplitz 型矩阵，当对角元素相等时，对角缩放失效。不完全 Cholesky 分解是一种有效的预条件器。可以通过抑制填充项或修改 Cholesky 算法来实现。例如，对于矩阵
$A =
\begin{bmatrix}
4 & -1 & -1 & 0 \
-1 & 4 & 0 & -1 \
-1 & 0 & 4 & -1 \
0 & -1 & -1 & 4
\end{bmatrix}$
标准 Cholesky 分解得到
$L_c =
\begin{bmatrix}
2.0000 & 0 & 0 & 0 \
-0.5000 & 1.9365 & 0 & 0 \
-0.5000 & -0.1291 & 1.9322 & 0 \
0 & -0.5164 & -0.5521 & 1.8516
\end{bmatrix}$
第一种不完全 Cholesky 分解版本为
$L_{i1} =
\begin{bmatrix}
2.0000 & 0 & 0 & 0 \
-0.5000 & 1.9365 & 0 & 0 \
-0.5000 & 0 & 1.9322 & 0 \
0 & -0.5164 & -0.5521 & 1.8516
\end{bmatrix}$
第二种版本为
$L_{i2} =
\begin{bmatrix}
2.0000 & 0 & 0 & 0 \
-0.5000 & 1.9365 & 0 & 0 \
-0.5000 & 0 & 1.9365 & 0 \
0 & -0.5164 & -0.5164 & 1.8619
\end{bmatrix}$
在实际应用中，可以设置阈值 $\epsilon$，当 $|l_{ij}| \leq \epsilon$ 时，令 $l_{ij} = 0$。不完全 Cholesky 分解只需在迭代循环前进行一次，并存储 $M = LL^T$ 或 $L$。
2.3 预条件器使用注意事项
- 虽然对称正定矩阵的 Cholesky 分解总是可行的，但不完全分解在某些情况下可能无法进行，可能会出现计算 $l_{ii}$ 时需要对负数开平方根的情况。此时，应使用一个大的正数覆盖相应元素，可以任意选择，也可以选择相邻对角元素的和或该行其余元素绝对值的和，以确保对角占优。
- 区域分解是另一种对偏微分方程进行预条件的方法。将感兴趣的区域分解为子区域，这些子区域可以重叠，然后在每个子区域中快速近似求解 PDE。预条件矩阵 $M$ 可以通过组合子问题的解来构建，如果子区域不重叠，则 $M$ 是块对角矩阵；如果子区域重叠，则 $M$ 是块对角子矩阵的乘积。
3. Toeplitz 矩阵和循环预条件器
循环矩阵作为一种特殊的预条件器，对 Toeplitz 矩阵非常有效。Toeplitz 矩阵具有恒定的对角线，常见于信号处理和许多时空不变的问题中。循环矩阵 $C$ 是 Toeplitz 矩阵，且其元素满足 $c_k = c_{k + n}$，其中 $n \times n$ 是矩阵 $C$ 的大小。
一般 Toeplitz 矩阵有 $(2n - 1)$ 个独立元素，由第一行和第一列确定；而一般循环矩阵只有 $n$ 个独立元素。假设对称情况下，Toeplitz 矩阵有 $n$ 个独立元素（第一列），循环矩阵有 $[\frac{n}{2}] + 1$ 个。具体区别如下表所示：
|矩阵类型|独立元素数量|元素表示|
| ---- | ---- | ---- |
|Toeplitz 矩阵|$(2n - 1)$ 个（一般），$n$ 个（对称）|$A =
\begin{bmatrix}
a_0 & a_{-1} & \cdots & \cdots & a_{1 - n} \
a_1 & a_0 & a_{-1} & \cdots & \
\vdots & \vdots & a_0 & \cdots & a_{-1} \
a_{n - 1} & \cdots & \cdots & a_1 & a_0
\end{bmatrix}$|
|循环矩阵|$n$ 个（一般），$[\frac{n}{2}] + 1$ 个（对称）|$C =
\begin{bmatrix}
c_0 & c_{n - 1} & \cdots & \cdots & c_1 \
c_1 & c_0 & c_{n - 1} & \cdots & \
\vdots & \vdots & c_0 & \cdots & c_{n - 1} \
c_{n - 1} & \cdots & \cdots & c_1 & c_0
\end{bmatrix}$|
对于对称 Toeplitz 矩阵 $A$ 构成的系统 $Ax = b$，可以用相应的循环矩阵进行预条件。在大多数应用中，主对角线及其相邻对角线占主导地位，可以用它们构建合适的循环预条件器 $C$。$c_1 = a_1$ 既出现在主对角线旁边，也出现在矩阵的角落。这些相对较大的元素控制着 $C^{-1}A$ 的特征谱，进而控制 PCGM 的收敛速率。其余特征值聚集在 1 附近。通常，我们从 Toeplitz 矩阵中复制一些主对角线到循环矩阵，但不是全部。
4. 并行 PCGM
预条件共轭梯度法的并行化相对简单。假设矩阵 $A$ 按行分布在各个进程中，算法的共轭梯度部分只需四个 MPI 调用：三个 $MPI_Allreduce$ 用于完成点积计算，一个 $MPI_Allgather$。如果使用残差与修改后残差的点积作为停止准则，可以减少一个 $MPI_Allreduce$ 调用。
根据预条件器的选择，可能需要额外的 MPI 调用进行预条件处理。对于对角预条件器，不需要额外的 MPI 调用，因为对角预条件可以在每个进程的所有行上本地完成。而不完全 Cholesky 等其他预条件器可能需要额外的 MPI 调用，在选择预条件器时需要考虑其成本。
以下是并行 PCG 算法迭代部分的示意图：
graph TD;
A[开始迭代] --> B[计算 z = Ap];
B --> C[计算局部点积 local - sum];
C --> D[MPI_Allreduce 得到全局和 sum];
D --> E[计算 alpha = c / sum];
E --> F[更新 x 和 r];
F --> G[预条件处理，计算 mr];
G --> H[计算局部点积 local - sum];
H --> I[MPI_Allreduce 得到全局和 sum];
I --> J[判断是否满足停止准则];
J -- 是 --> K[结束迭代];
J -- 否 --> L[计算 beta];
L --> M[更新 p];
M --> N[MPI_Allgather 更新 p];
N --> O[更新 c];
O --> B;
4.1 MPI 程序示例
我们通过一个 MPI 程序演示使用对角预条件器的 PCGM。示例问题是求解以下方程：
$\frac{d^2u(x)}{dx^2} - c_2e^{c_1(x - 0.5)^2}u(x) = -\sin(2\pi x)(4\pi^2 - c_2e^{c_1(x - 0.5)^2})$
在区间 $x \in [0, 1]$ 上，边界条件为 $u(0) = u(1) = 0$，常数 $c_1 = 20.0$，$c_2 = 1000.0$，精确解为 $u(x) = \sin(2\pi x)$。
为了更好地解释代码，将整个程序分为四个部分：
4.1.1 第一部分：MPI 初始化、初始内存分配和网格及右侧向量生成
#include <iostream.h>
#include <iomanip.h>
#include "SCmathlib.h"
#include<mpi.h>
const int rows-per-proc = 40;
const double c1 = 20.0;
const double c2 = 1000.0;
const double tol = 1.0e-14;
int main(int argc, char *argv[]){
int i,j,k;
int mynode, totalnodes, totalsize, offset;
MPI-Status status;
double sum,local-sum,c,d,alpha,beta;
double ** A, *q, *x, *grid;
double *p,*z,*r,*mr;
MPI-Init(&argc,&argv);
MPI-Comm-size(MPI-COMM-WORLD, &totalnodes);
MPI-Comm-rank(MPI-COMM-WORLD, &mynode);
totalsize = totalnodes*rows-per-proc;
p = new double[totalsize];
z = new double[rows-per-proc];
r = new double[rows-per-proc];
mr = new double[rows-per-proc];
x = new double[rows-per-proc];
q = new double[rows-per-proc];
grid = new double[rows-per-proc];
double dx = 1.0/(totalsize+1);
for(i=0;i<rows-per-proc;i++){
grid[i] = dx*(1+rows-per-proc*mynode+i);
q[i] = -dx*dx*sin(2.0*M_PI*grid[i])*
(-4.0*M_PI*M_PI - c2*exp(c1*(grid[i]-0.5)*
(grid[i]-0.5)));
x[i] = 1.0;
}
}
注意事项
：
- 程序中有四个全局变量：
rows_per_proc
表示每个处理器的行数，
c1
和
c2
是问题特定的常数，
tol
是收敛所需的容差。
- 除了数组
p
外，其他数组只需
rows_per_proc
大小，每个处理器只需维护其部分的残差、修改后的残差和解决方案向量，但每个处理器必须有整个
p
向量的副本。
4.1.2 第二部分：内存分配和矩阵 $A$ 的生成
/* Part 2 */
A = new double*[rows-per-proc];
for(i=0;i<rows-per-proc;i++){
A[i] = new double[totalsize];
for(j=0;j<totalsize;j++)
A[i][j] = 0.0;
}
if(mynode==0){
A[0][0] = 2.0 + dx*dx*c2*exp(c1*(grid[0]-0.5)*
(grid[0]-0.5));
A[0][1] = -1.0;
for(i=1;i<rows-per-proc;i++){
A[i][i] = 2.0 + dx*dx*c2*exp(c1*(grid[i]-0.5)*
(grid[i]-0.5));
A[i][i-1] = -1.0;
A[i][i+1] = -1.0;
}
}
else if(mynode == (totalnodes-1)){
A[rows-per-proc-1][totalsize-1] = 2.0 +
dx*dx*c2*exp(c1*(grid[rows-per-proc-1]-0.5)*
(grid[rows-per-proc-1]-0.5));
A[rows-per-proc-1][totalsize-2] = -1.0;
for(i=0;i<rows-per-proc-1;i++){
offset = rows-per-proc*mynode + i;
A[i][offset] = 2.0 + dx*dx*c2*exp(c1*(grid[i]-0.5)*
(grid[i]-0.5));
A[i][offset-1] = -1.0;
A[i][offset+1] = -1.0;
}
}
else{
for(i=0;i<rows-per-proc;i++){
offset = rows-per-proc*mynode + i;
A[i][offset] = 2.0 + dx*dx*c2*exp(c1*(grid[i]-0.5)*
(grid[i]-0.5));
A[i][offset-1] = -1.0;
A[i][offset+1] = -1.0;
}
}
注意事项 ：矩阵设置分为三种情况，需要仔细处理第一个和最后一个处理器，因为它们分别包含第一行和最后一行。
4.1.3 第三部分：PCGM 初始化
/* Part 3 */
offset = mynode*rows-per-proc;
for(i=0;i<totalsize;i++)
p[i] = 1.0;
for(i=0;i<rows-per-proc;i++){
r[i] = q[i] - dot(totalsize,A[i],p); //计算残差
mr[i] = r[i]/A[i][offset+i]; //计算修改后的残差
}
local-sum = dot(rows-per-proc,mr,r);
MPI-Allreduce(&local-sum,&sum,1,MPI_DOUBLE,MPI_SUM,
MPI_COMM_WORLD);
c = sum;
MPI-Allgather(mr,rows-per-proc,MPI_DOUBLE,p,rows-per-proc,
MPI_DOUBLE,MPI_COMM_WORLD);
注意事项
：初始向量选择为全 1 向量。由于需要计算初始残差，所有处理器都需要完整的初始猜测。这里临时使用
p
数组来实现这一点，而不是分配新的向量，因为在计算修改后的残差之前
p
未被使用，不会产生不良影响。
4.1.4 第四部分：PCGM 主迭代循环
/* Part 4 */
for(k=0;k<totalsize;k++){
for(i=0;i<rows-per-proc;i++)
z[i] = dot(totalsize,A[i],p);
local-sum = dot(rows-per-proc,z,p+offset);
MPI-Allreduce(&local-sum,&sum,1,MPI_DOUBLE,MPI_SUM,
MPI_COMM_WORLD);
alpha = c/sum;
for(i=0;i<rows-per-proc;i++){
x[i] = x[i] + alpha*p[offset+i];
r[i] = r[i] - alpha*z[i];
}
/* 预条件处理阶段 */
for(i=0;i<rows-per-proc;i++)
mr[i] = r[i]/A[i][offset+i];
local-sum = dot(rows-per-proc,mr,r);
MPI-Allreduce(&local-sum,&sum,1,MPI_DOUBLE,MPI_SUM,
MPI_COMM_WORLD);
d = sum; // 包含残差和修改后残差的内积
local-sum = dot(rows-per-proc,r,r);
MPI-Allreduce(&local-sum,&sum,1,MPI_DOUBLE,MPI_SUM,
MPI_COMM_WORLD);
// sum 现在包含残差和残差的内积
if(mynode == 0){
cout << k << "\t" << "dot(mr,r) = " << d << "\t";
cout << "dot(r,r) = " << sum << endl;
}
if(fabs(d) < tol) break;
if(fabs(sum) < tol) break;
beta = d/c;
for(i=0;i<rows-per-proc;i++)
z[i] = mr[i] + beta*p[i+offset];
MPI-Allgather(z,rows-per-proc,MPI_DOUBLE,p,rows-per-proc,
MPI_DOUBLE,MPI_COMM_WORLD);
c = d;
}
delete[] p;
delete[] z;
delete[] r;
delete[] mr;
delete[] x;
delete[] q;
delete[] grid;
for(i=0;i<rows-per-proc;i++)
delete[] A[i];
delete[] A;
MPI_Finalize();
注意事项
：
- 只需要四个 MPI 调用：三个
MPI_Allreduce
和一个
MPI_Allgather
。三个
MPI_Allreduce
用于获取所有处理器的内积，
MPI_Allgather
用于确保所有处理器都有更新后的
p
值。
- 由于使用对角预条件器，预条件器不需要额外的通信。如果使用不完全 Cholesky 作为预条件器，则需要额外的通信。
通过上述示例可以看出，PCGM 比标准 CGM 收敛速度快得多。同时，使用二阶有限差分方法近似导数算子时，误差应该随着网格点数的加倍而减小 4 倍，这与理论分析一致。
4.2 结果分析
4.2.1 误差分析
在这个示例中，我们使用二阶有限差分方法来近似导数算子。根据理论，当我们将网格点数加倍时，误差应该减小 4 倍。下面是 L2 误差与网格点数的对数 - 对数图：
graph LR;
A[网格点数] --> B[L2 误差];
style A fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
style B fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
从图中可以看出，直线的斜率大小约为 2，这与我们使用的二阶近似方法是一致的。这表明随着网格点数的增加，误差按照预期的速度减小。
4.2.2 收敛速度比较
我们还绘制了 CGM 和 PCGM 的残差内积 $(r, r)$ 随迭代次数的变化图：
graph LR;
A[迭代次数] --> B[CGM 残差内积];
A --> C[PCGM 残差内积];
style A fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
style B fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
style C fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
从图中可以明显看出，PCGM 的收敛速度比标准 CGM 快得多。这是因为预条件器有效地改善了矩阵的条件数，使得算法能够更快地收敛到解。
4.3 代码优化建议
- 减少 MPI 通信 ：虽然并行 PCG 算法已经通过合理的设计减少了 MPI 调用次数，但在某些情况下，仍然可以进一步优化。例如，可以考虑使用非阻塞通信来重叠计算和通信，从而提高效率。
- 选择合适的预条件器 ：不同的预条件器在不同的问题上表现不同。在实际应用中，需要根据矩阵的特性和问题的规模选择合适的预条件器。例如，对于稀疏矩阵，不完全 Cholesky 分解可能是一个不错的选择；而对于对角占优的矩阵，对角预条件器可能更合适。
- 内存管理 ：在处理大规模矩阵时，内存管理非常重要。可以考虑使用稀疏矩阵存储格式来减少内存占用，同时避免不必要的内存分配和释放操作。
5. 总结
5.1 关键知识点回顾
- 共轭梯度法（CGM） ：收敛速率受矩阵条件数影响，对于某些问题，迭代次数与矩阵条件数的平方根成正比。
- 预条件共轭梯度法（PCGM） ：通过引入预条件器来改善矩阵的条件数，从而加速收敛。常见的预条件器包括对角缩放、块对角预条件器和不完全 Cholesky 分解等。
- Toeplitz 矩阵和循环预条件器 ：循环矩阵作为预条件器对 Toeplitz 矩阵非常有效，能够控制特征谱，加快收敛速度。
- 并行 PCGM ：可以通过 MPI 实现并行化，根据预条件器的选择，可能需要额外的 MPI 调用。
5.2 实际应用建议
- 问题分析 ：在使用 PCGM 之前，需要对问题进行充分的分析，了解矩阵的特性，如对角元素的差异、稀疏性等，以便选择合适的预条件器。
- 代码实现 ：在实现并行 PCGM 时，要注意 MPI 调用的使用，合理分配内存，避免不必要的通信和计算开销。
- 结果验证 ：通过绘制误差图和收敛曲线，验证算法的正确性和性能，确保达到预期的效果。
5.3 未来研究方向
- 新型预条件器的研究 ：不断探索和开发更有效的预条件器，以适应不同类型的矩阵和问题。
- 并行算法的优化 ：进一步优化并行 PCGM 算法，提高其在大规模并行计算环境下的性能。
- 与其他算法的结合 ：将 PCGM 与其他数值算法相结合，以解决更复杂的问题。
总之，预条件共轭梯度法是一种强大的线性求解器，通过合理选择预条件器和并行化实现，可以在许多实际问题中取得良好的效果。希望本文能够帮助读者更好地理解和应用这一算法。
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
