41、快速线性求解器：原理、选择与实践

快速线性求解器：原理、选择与实践

1. 回代求解系数

在求解线性方程组时，常常需要进行回代操作以获得系数。以下是一段实现回代过程的代码：

z.Initialize(0.0);
if(i>=(m-1)){
    nr = m;
    z(nr-1) = y(nr-1)/H(nr-1,nr-1);
}
for(k=nr-2;k>=0;k--){
    z(k) = y(k);
    for(ll=k+1;ll<nr;ll++)
        z(k) -= H(k,ll)*z(ll);
    z(k) = z(k)/H(k,k);
}

这段代码的主要作用是通过回代的方式，从已知的向量 y 和矩阵 H 中计算出向量 z 的各个元素。具体步骤如下：
1. 初始化向量 z 的所有元素为 0.0。
2. 如果 i 大于等于 m - 1 ，则将 nr 设为 m ，并计算 z 的最后一个元素。
3. 从 nr - 2 开始，逆向遍历到 0，对于每个 k ，先将 z(k) 初始化为 y(k) ，然后减去 H(k, ll) * z(ll) 的累加和，最后将结果除以 H(k, k) 。

2. Krylov 空间基向量的线性组合

计算出系数向量 z 后，需要进行 Krylov 空间基向量的线性组合，以得到最终的解向量 x 。代码如下：

for(i=0;i<nr;i++){
    V.GetColumn(i,r);
    x = x + z(i)*r;
}

这段代码的步骤如下：
1. 遍历从 0 到 nr - 1 的每个索引 i 。
2. 从矩阵 V 中获取第 i 列向量，存储在向量 r 中。
3. 将向量 r 乘以系数 z(i) ，并累加到向量 x 上。

3. 残差判断与更新

在得到解向量 x 后，需要判断残差是否满足要求，如果不满足则更新残差向量。代码如下：

if(rho<tol)
    break;
r = b - A*x;

步骤如下：
1. 检查残差 rho 是否小于容差 tol ，如果是，则跳出循环。
2. 计算新的残差向量 r ，即 b - A * x 。

4. 注意事项

在使用上述方法时，有两个重要的注意事项：
- 存储问题 ：该函数需要存储矩阵 H 和 V ，它们的大小可能与原始矩阵 A 相同。这是引入重启参数的主要动机之一，以便使用更小的 Krylov 空间，从而减少存储需求。
- Givens 矩阵的处理 ：可以利用矩阵的结构，避免显式地存储旋转矩阵，而是通过循环遍历适当的位置进行更新。

5. 预条件 GMRES

当矩阵 A 条件数不佳时，GMRES 方法通常会收敛缓慢，此时需要使用预条件器。这里使用不完全 LU（ILU）作为预条件器，具体形式为 M = LU ，其中 L 和 U 是对应于矩阵 A 的下三角和上三角矩阵，但在 aij = 0 的位置不进行填充。

预条件器不应显式构造，而应融入 Arnoldi 迭代过程。需要在修改后的 Arnoldi 迭代算法中插入以下代码：

...
My = v j
w = Ay
for i = 1, . . . , j
    hij = (w, vi)
    w = w − hijvi
...

GMRES 使用长向量来实现正交性，与对称系统的三项递推公式不同。对于非对称系统，也可以使用三项递推公式，如双共轭梯度（BiCG）方法，它使用两个相互正交的向量序列。BiCG 的更稳定版本是准最小残差（QMR）方法，它可以避免可能的崩溃，并且收敛速度与 GMRES 相当。

6. 并行 GMRES

尝试并行化串行的 GMRES 算法时，一个直接的问题是并行化修改后的 Arnoldi 组件。暴力实现的可扩展性很差。可以采用以下两种方法来提高效率和并行性：
- 块操作 ：首先生成 Krylov 空间的基向量 v1, Av1, ..., Akv1 ，然后对整个集合进行正交化。与标准 GMRES 方法中每个新向量立即与所有先前向量正交化不同，这种方法显著提高了数据局部性。
- 使用 BLAS2 例程 ：在 GMRES 中尽可能使用 BLAS2 例程，而不是效率最低的 BLAS1 例程。一种算法是用标准 Gram - Schmidt 方法替换修改后的 Gram - Schmidt 方法，并应用两次。双重正交化已被证明可以减少经典 Gram - Schmidt 方法的数值不稳定性，并且可以并行计算所有点积。

7. 求解器的选择

选择合适的求解器和预条件器是一个复杂的问题，需要考虑多个因素：
- 矩阵性质 ：包括矩阵是否对称、是否正定、 A 和 AT 是否可用、矩阵是否稀疏、矩阵是否病态以及是否存在好的预条件器。
- 计算要求和资源 ：如并行或串行计算、向量或标量处理器、多线程、缓存数据的重用、间接寻址和内存大小。

以下是一些迭代求解器的主要操作和存储需求的表格：
| 方法 | ddot | daxpy | dgemv | 存储 |
| ---- | ---- | ---- | ---- | ---- |
| Jacobi | 1 | 3n | | |
| SOR | 1 | 1 | | 2n |
| CGM | 2 | 3 | 1 | 6n |
| GMRES | j + 1 | j + 1 | 1 | (j + 5)n |
| QMR | 2 | 12 | 2 | 16n |

对于涉及微分方程的问题，所处理的微分方程类型和选择的离散化方法决定了要求解的线性系统。例如，对于 Poisson 和 Helmholtz 方程，会得到对称系统。以下是不同直接和迭代求解器在求解 Poisson 方程时的计算工作量和存储需求的表格：
| 方法 | 直接/迭代 | 工作量 | 存储 |
| ---- | ---- | ---- | ---- |
| 多重网格 | I | n | n |
| SSOR/Chebyshev | I | n5/4 | n |
| SOR | I | n3/2 | n |
| CGM | I | n3/2 | n |
| Gauss - Seidel | I | n2 | n |
| Jacobi | I | n2 | n |
| 高斯消元/稀疏 | D | n3/2 | n · log n |
| 高斯消元/稠密 | D | n3 | n2 |

对于大型问题，由于内存考虑，通常建议使用迭代求解器。以下是一个可能的决策树：

graph TD;
    A{矩阵 A 是否对称?} -->|是| B{矩阵 A 是否正定?};
    A -->|否| C{AT 是否可用?};
    B -->|是| D[预条件共轭梯度法];
    B -->|否| E[尝试 CGNR];
    C -->|是| F[GMRES 或 GMRES(k)];
    C -->|否| G[QMR 算法];

8. 可用的软件

基本的算法相对容易编程，但更复杂的版本，如预条件、重启、正交化和并行实现等，可以在免费软件 www.netlib.org 中找到。
- 直接求解器 ：是 ScaLAPACK 和 LAPACK++ 的一部分。例如， SGETRF 例程基于 BLAS3 例程执行带主元的 LU 分解，非常高效； SPOTRF 和 SPOTRS 例程分别对对称正定矩阵进行 Cholesky 分解和求解线性系统。
- 迭代求解器 ：是 Templates 包的一部分。例如， cg.h 和 cheby.h 例程分别实现了对称正定系统的共轭梯度法和预条件 Chebyshev 方法； cgs.h 、 gmres.h 、 qmr.h 和 bicg.h 适用于非对称系统。

9. 作业问题

为了更好地理解和掌握上述内容，以下是一些作业问题：
1. 求矩阵 A = [[0.001, 1], [1, 1]] 和 B = [[7, 6.990], [4, 4]] 的条件数。
2. 设 A = L DLT 是对称正定矩阵， D = diag (dii) ，证明 κ2(A) ≥ max(dii) / min(dii) 。
3. 对于矩阵 [[1, ϵ], [ϵ, 1]] ，当 ϵ 取何值时矩阵病态？单精度或双精度计算会对结果产生什么影响？
4. 假设矩阵 A 是严格对角占优的，证明对 A 应用带部分主元的 LU 分解过程时，不会发生行交换。
5. 估计几乎三角矩阵（Hessenberg 矩阵）的 LU 分解和回代求解的运算次数。
6. 考虑矩阵 A 以 2×2 块形式表示，证明 Schur 补 S ≡ A22 - A21A−111 A12 在无主元高斯消元 m 步后会覆盖矩阵 A22 。
7. 手动计算一个 3×3 矩阵（ ri = i + j - 1 ）的 Hager 算法，并计算其在 L1 范数下的精确条件数，比较两个值。
8. 对 Hilbert 矩阵（ hij = 1 / (i + j - 2) ），当 n = 4, 16, 32 时，应用 Hager 算法估计条件数，观察结果。
9. 修改并行 LU 程序以实现行主元。
10. 以并行 LU 程序为指导，实现一个 MPI 程序来完成 Cholesky 分解。
11. 对三对角矩阵 (-1, 4, -1) 应用不完全 Cholesky 分解，并与标准 Cholesky 分解结果进行比较，判断结论是否适用于任何带状矩阵。
12. 修改循环约简代码，使存储三对角矩阵 A 只使用必要的内存。
13. 手动使用循环约简方法求解系统 {-2x1 + x2 = F1, x1 - 2x2 + 2x3 = F2, ..., x6 - 2x7 = F7} ，确定系统终止的级别和方程，判断是否有足够信息获得解，如有则进行回代求解。
14. 使用 Householder 变换证明，如果 He 是 Hessenberg 矩阵且 He = Q R ，则矩阵 H∗e = RQ 也是 Hessenberg 矩阵。
15. 求解超定系统 Ax = b （ A 为 m×n 矩阵）的浮点运算次数选择：
- (a) 正规方程为 mn2/2 + n3/6 ，Householder QR 方法为 mn2/2 - n3/6 。
- (b) 与 (a) 相反。
- (c) 以上都不是。
16. 使用 Householder 方法对 n×n 矩阵进行带列主元的 QR 分解的浮点运算次数选择：
- (a) 5/3n3 。
- (b) 5/4n3 。
- (c) 4/5n2 。
17. Cholesky 分解的加法和乘法次数大约是 LU 分解的一半，判断对错。
18. QR 方法是否总是、有时或从不保留矩阵的带宽，选择正确答案。
19. 矩阵 A 的 QR 分解是否总是唯一、当 A 非奇异时唯一或不唯一，选择正确答案。
20. 在并行机器上分配矩阵最有效的方式是按行、按列还是取决于要解决的问题，选择正确答案。
21. 定义 h = 1 / N ， µ = h2 ， q[i][j] = (8π2 + 1)h2sin(2πhi)sin(2πhj) ，使用以下方法求解矩阵系统 Au = q ：
- (a) 串行 Jacobi 方法。
- (b) 并行 Jacobi 方法。
- (c) 串行共轭梯度法。
- (d) 串行预条件共轭梯度法（使用不完全 Cholesky 预条件）。
- (e) 并行共轭梯度法。
- (f) 并行预条件共轭梯度法（使用不完全 Cholesky 预条件）。
分别对 N = 4 和 N = 20 进行求解，对于串行算法，绘制解 u[i][j] 的图形；对于并行算法，展示不同处理器数量下的并行加速比。
22. 使用二阶有限差分离散化，在三维网格 N×N×N 上求解 Poisson 方程，采用以下串行算法：
- (a) 共轭梯度法。
- (b) 带不完全 Cholesky 预条件的共轭梯度法。
- (c) 带 Chebyshev 加速的 SSOR 方法。
- (d) Jacobi 方法。
估计计算工作量的阶数 O(nα) ，并通过求解 n = N3 = 643 验证估计结果。
23. 根据文本中的 QR 分解示例，从找到的 w 值计算矩阵 Q ，并证明 A = QR 。
24. 修改文本中的 QR 例程，使其计算矩阵 Q 的值。
25. 编写一个串行函数，使用 Givens 旋转对三对角矩阵进行 QR 分解，函数应接受矩阵 A 作为输入，矩阵 Q 和 R 作为输出。
26. 创建一个函数，使用前面编写的 QR 代码求解三对角矩阵系统 Ax = b ，并说明如何利用 Q 和 R 的性质高效求解。
27. 编写一个并行程序，使用 Givens 旋转对三对角矩阵进行 QR 分解，按行将矩阵 A 分配到各个处理器上，确保每个处理器拥有与 A 对应行的 Q 和 R 行。
28. 考虑五个对称 Toeplitz 矩阵 A （ a(1)k = 1/k ， a(2)k = 1/√k ， a(3)k = 1/k2 ， a(4)k = k ， a(5)k = cos k/k ）和任意非零向量 b ，使用 PCGM 方法，采用循环预条件器求解系统 Ax = b ，并尝试不同类型的预条件器，比较无 CG 预条件时的迭代次数。
29. 估计固定 k 时 GMRES(k) 的运算次数，并与 GMRES 进行比较，假设矩阵 A 的阶数 n 很大，判断 GMRES(k) 是否对任何 k 值都收敛。
30. 对于文本中的 GMRES 示例问题，尝试使用 GMRES(2)，添加输出语句报告残差，尝试不同的初始猜测 x0 ，绘制每种情况下残差与迭代次数的关系图，判断收敛速度是否相同。
31. 考虑线性对流 - 扩散方程 ∂u/∂t + ∂u/∂x = ν ∂2u/∂x2 ，在 x ∈ [0, 10] 和 t ∈ [0, 5] 区域内，采用二阶迎风格式处理对流项，中心有限差分格式处理扩散项（时间上采用 Crank - Nicolson 方法）进行离散化，使用 GMRES(k) 求解得到的系统，尝试不同的 k 值以最小化求解时间，使用 n = 128 个点进行离散化，判断最优 k 值是否依赖于时间步长 Δt ，以及该方法是否无条件稳定。

快速线性求解器：原理、选择与实践

10. 问题分析与解答示例

下面针对部分作业问题进行详细分析与解答，帮助大家更好地理解相关概念和方法。

10.1 矩阵条件数的计算

对于矩阵 (A = \begin{bmatrix}0.001 & 1\1 & 1\end{bmatrix}) 和 (B = \begin{bmatrix}7 & 6.990\4 & 4\end{bmatrix})，矩阵的条件数可以通过计算矩阵的最大奇异值与最小奇异值的比值来得到。在实际计算中，可以使用数值计算库（如 Python 的 numpy 库）来完成。以下是 Python 代码示例：

import numpy as np

A = np.array([[0.001, 1], [1, 1]])
B = np.array([[7, 6.990], [4, 4]])

cond_A = np.linalg.cond(A)
cond_B = np.linalg.cond(B)

print(f"矩阵 A 的条件数: {cond_A}")
print(f"矩阵 B 的条件数: {cond_B}")

通过上述代码，可以方便地计算出矩阵 (A) 和 (B) 的条件数。

10.2 对称正定矩阵条件数的证明

设 (A = L D L^T) 是对称正定矩阵，(D = diag (d_{ii}))。根据条件数的定义，(\kappa_2(A) = \frac{\lambda_{max}(A)}{\lambda_{min}(A)})，其中 (\lambda_{max}(A)) 和 (\lambda_{min}(A)) 分别是矩阵 (A) 的最大和最小特征值。

由于 (A) 是对称正定矩阵，其特征值都是正的。又因为 (A = L D L^T)，可以证明 (\lambda_{max}(A) \geq \max(d_{ii})) 且 (\lambda_{min}(A) \leq \min(d_{ii}))，所以 (\kappa_2(A) \geq \frac{\max(d_{ii})}{\min(d_{ii})})。

10.3 矩阵病态性的判断

对于矩阵 (C = \begin{bmatrix}1 & \epsilon\\epsilon & 1\end{bmatrix})，其条件数为 (\kappa_2(C) = \frac{1 + |\epsilon|}{1 - |\epsilon|})。当 (\epsilon) 接近 (\pm1) 时，条件数会变得很大，矩阵就会变得病态。

在单精度或双精度计算中，由于浮点数的精度限制，当矩阵病态时，计算结果的误差会显著增大。双精度计算的精度更高，相对来说能在一定程度上减少误差，但对于病态矩阵，仍然可能会出现较大的误差。

11. 不同求解方法的比较与应用

在实际应用中，需要根据具体问题的特点选择合适的求解方法。以下是几种常见求解方法的比较：

方法	适用矩阵类型	收敛速度	存储需求	并行性
Jacobi	对角占优矩阵	较慢	较小	较好
SOR	对角占优矩阵	比 Jacobi 快	较小	较好
CGM	对称正定矩阵	较快	适中	较好
GMRES	非对称矩阵	适中	较大	较难并行化
QMR	非对称矩阵	较快	较大	较难并行化

对于不同类型的问题，可以参考上述表格选择合适的求解方法。例如，对于对称正定矩阵的问题，CGM 是一个不错的选择；对于非对称矩阵的问题，可以考虑 GMRES 或 QMR 方法。

12. 并行计算的优化策略

在并行计算中，为了提高计算效率，需要采用一些优化策略。以下是一些常见的优化策略：

• 数据划分 ：根据问题的特点，合理地划分矩阵和向量，使得每个处理器能够独立地进行计算。常见的数据划分方式有按行划分、按列划分和分块划分。
• 通信优化 ：减少处理器之间的通信量，避免频繁的数据传输。可以采用局部计算和全局同步的方式，减少通信开销。
• 算法选择 ：选择适合并行计算的算法，如并行 GMRES 中采用的块操作和使用 BLAS2 例程的方法。

13. 总结

快速线性求解器在科学计算和工程应用中具有重要的作用。通过选择合适的求解器和预条件器，可以显著提高计算效率和精度。在实际应用中，需要综合考虑矩阵的性质、计算要求和资源等因素，选择最适合的求解方法。

同时，并行计算是提高计算效率的重要手段，但并行化过程中也会遇到一些挑战，如数据通信和负载均衡等问题。需要采用合适的优化策略，提高并行计算的性能。

通过完成作业问题，可以加深对快速线性求解器的理解和掌握，提高解决实际问题的能力。希望大家在学习和实践中不断探索，更好地应用这些方法解决各种问题。

graph LR;
    A[问题分析] --> B[选择求解方法];
    B --> C[串行或并行计算];
    C --> D[结果验证与优化];
    D --> E[应用到实际问题];

以上流程图展示了使用快速线性求解器解决问题的一般流程，从问题分析开始，选择合适的求解方法，进行计算，然后验证和优化结果，最后将结果应用到实际问题中。这个流程可以帮助大家系统地解决各种线性方程组问题。