Python实现CG共轭梯度算法详解

"CG共轭梯度算法是求解线性方程组 Ax = b 的一种迭代方法，尤其适用于大型稀疏矩阵问题。在Python中，可以使用NumPy库进行实现。本文档提供了两个CG算法的实现版本：CG2 和 CG。这两个函数接受矩阵 A、向量 b 和初始猜测向量 x 作为输入参数，以及可选的迭代次数上限 imax 和收敛阈值 epsilon。CG算法的核心思想是寻找与当前残差正交的搜索方向，并通过阿伦茨公式更新解向量。” 在CG算法中，有以下几个关键步骤： 1. **初始化**: 创建一个存储迭代步长的数组 `steps`，设置迭代计数器 `i` 为0，计算初始残差 `r = b - A * x`，这同时也作为负梯度方向。 2. **计算搜索方向**: 搜索方向 `d` 初始化为残差 `r` 的副本。在CG2版本中，使用 `np.dot(r.T, r)` 计算残差的平方范数 `delta_new` 作为初始值。 3. **主循环**: 在满足迭代次数不超过 `imax` 和残差平方范数大于 `epsilon` 倍的初始残差平方范数 `delta_0` 的条件下，执行以下操作： - 计算对称矩阵乘积 `q = A * d`。 - 使用阿伦茨公式计算步长 `alpha`，即 `alpha = delta_new / (d.T * q)`。CG2版本中，`np.float` 用于确保 `alpha` 是实数。 - 更新解向量 `x = x + alpha * d`。 - 每隔50次迭代，重新计算残差 `r = b - A * x`。 - 若不需重新计算残差，更新残差 `r = r - alpha * q`。 - 计算新的残差平方范数 `delta_new`，并更新搜索方向 `d = r + beta * d`，其中 `beta` 是根据残差变化率计算的。 4. **结束条件**: 循环结束后，返回包含所有迭代步长的数组 `steps`。 这两个函数的主要区别在于： - CG2 版本使用 `np.linalg.norm` 来比较残差范数，而原始CG版本使用的是残差的平方范数。 - CG2 版本在更新残差时考虑了是否需要每隔50次迭代重计算，这可能会影响算法的效率。 CG算法的优点在于其对大型稀疏矩阵的高效处理，因为每次迭代仅需要矩阵向量乘法，而不是完整的矩阵逆或行列式计算，从而大大降低了计算复杂度。在实际应用中，它常被用于求解大规模科学计算和工程问题中的线性系统。