告别梯度下降的震荡：用Python手把手实现共轭梯度法（CG）求解线性方程组

告别梯度下降的震荡：用Python手把手实现共轭梯度法（CG）求解线性方程组

在数值计算领域，求解大型线性方程组是一个永恒的话题。当矩阵规模达到数千甚至数万维度时，传统的直接解法如LU分解会消耗惊人的内存和计算资源。这时，迭代法便成为更明智的选择——它们不需要存储完整的矩阵，只需矩阵向量乘积运算，特别适合稀疏矩阵场景。

共轭梯度法（Conjugate Gradient, CG）作为迭代法中的明星算法，兼具理论优雅与实用价值。它最初由Hestenes和Stiefel在1952年提出，专门用于对称正定矩阵系统。与最速下降法相比，CG通过精心设计的共轭方向搜索策略，能有效避免"之字形"震荡收敛，通常能在远小于矩阵维度的迭代次数内获得满意解。

本文将带您从零实现CG算法，重点解决三个核心问题：

1. 如何构造相互共轭的搜索方向？
2. 怎样确定最优步长保证每次迭代最大程度降低误差？
3. 何时终止迭代才能平衡精度与效率？

我们将用Python和NumPy搭建完整求解框架，通过可视化对比揭示CG的收敛优势，最终封装成可直接调用的工业级求解器。所有代码均提供详细注释，确保即使线性代数基础薄弱的读者也能理解每个步骤的设计思想。

1. 理论基础：共轭梯度法为何高效

1.1 从二次优化看线性方程组

考虑对称正定矩阵A∈ℝⁿˣⁿ和向量b∈ℝⁿ构成的线性方程组Ax=b。这等价于求解如下二次函数的极小值点：

def quadratic_form(A, b, x): return 0.5 * x.T @ A @ x - b.T @ x

该函数的梯度∇ϕ(x)=Ax-b恰好是方程组的残差。当A的条件数κ(A)=λ_max/λ_min较大时，最速下降法会因反复横跨狭窄山谷而收敛缓慢。下图展示了二维情况下的典型震荡路径：

1.2 共轭方向的魔力

一组向量{p₀,p₁,...,pₖ}关于A共轭意味着：

pᵀApⱼ = 0, ∀i≠j

这种精心设计的方向具有超正交性，能保证：

• 每个方向上的优化互不干扰
• 最多n步即可收敛到精确解（理论值）
• 实际中常在O(√κ(A))次迭代内达到满意精度

关键性质证明：

设误差eₖ=xₖ-x∗，则有∥eₖ∥ₐ ≤ 2[(√κ-1)/(√κ+1)]ᵏ∥e₀∥ₐ
其中∥v∥ₐ=√(vᵀAv)是能量范数

2. 算法实现：从公式到代码

2.1 CG算法核心步骤

完整CG算法的伪代码如下：

输入：A, b, x₀, max_iter, tol r₀ = b - A @ x₀ p₀ = r₀ for k in 0..max_iter: αₖ = (rₖᵀ @ rₖ) / (pₖᵀ @ A @ pₖ) xₖ₊₁ = xₖ + αₖ * pₖ rₖ₊₁ = rₖ - αₖ * A @ pₖ if ∥rₖ₊₁∥ < tol: break βₖ = (rₖ₊₁ᵀ @ rₖ₊₁) / (rₖᵀ @ rₖ) pₖ₊₁ = rₖ₊₁ + βₖ * pₖ return xₖ₊₁

2.2 Python实现细节

import numpy as np from scipy.sparse import csr_matrix def conjugate_gradient(A, b, x0=None, max_iter=1000, tol=1e-6): """共轭梯度法求解Ax=b 参数： A : 对称正定矩阵(n,n) b : 右端向量(n,) x0 : 初始猜测(n,)，默认为零向量 max_iter : 最大迭代次数 tol : 残差容限 返回： x : 解向量 res : 残差历史记录 """ n = len(b) x = np.zeros(n) if x0 is None else x0.copy() r = b - (A @ x) p = r.copy() res = [np.linalg.norm(r)] for _ in range(max_iter): Ap = A @ p alpha = (r @ r) / (p @ Ap) x += alpha * p r_new = r - alpha * Ap if np.linalg.norm(r_new) < tol: break beta = (r_new @ r_new) / (r @ r) p = r_new + beta * p r = r_new res.append(np.linalg.norm(r)) return x, np.array(res)

关键优化技巧：

• 使用稀疏矩阵存储（如CSR格式）提升大矩阵运算效率
• 避免重复计算矩阵乘积，预先存储Ap
• 采用相对残差∥rₖ∥/∥b∥作为终止条件更鲁棒

3. 数值实验：CG vs 最速下降法

3.1 测试案例设计

构造不同条件数的希尔伯特矩阵作为测试用例：

def generate_test_case(n, kappa): """生成条件数为kappa的测试矩阵""" U = np.random.randn(n, n) Q, _ = np.linalg.qr(U) # 随机正交矩阵 lam = np.linspace(1, kappa, n) # 等间距特征值 A = Q @ np.diag(lam) @ Q.T x_true = np.random.randn(n) b = A @ x_true return A, b, x_true

3.2 收敛性对比

设置n=100, κ=1e4，比较两种方法的残差下降速度：

可视化结果显示，CG法（蓝色）呈近似线性收敛，而最速下降法（红色）呈现典型震荡：

3.3 条件数敏感性分析

固定n=200，变化条件数κ观察迭代次数：

kappa_list = [10**i for i in range(1, 6)] iters = [] for kappa in kappa_list: A, b, _ = generate_test_case(200, kappa) _, res = conjugate_gradient(A, b) iters.append(len(res))

结果验证理论预测——迭代次数与√κ成正比：

4. 工程实践：工业级优化技巧

4.1 预处理技术

当κ(A)极大时，可通过预处理矩阵M≈A⁻¹加速收敛。修改后的PCG算法：

def preconditioned_cg(A, b, M, x0=None, max_iter=1000, tol=1e-6): x = np.zeros_like(b) if x0 is None else x0.copy() r = b - A @ x z = M @ r # 预处理步骤 p = z.copy() res = [np.linalg.norm(r)] for _ in range(max_iter): Ap = A @ p alpha = (r @ z) / (p @ Ap) x += alpha * p r_new = r - alpha * Ap if np.linalg.norm(r_new) < tol: break z_new = M @ r_new beta = (r_new @ z_new) / (r @ z) p = z_new + beta * p r, z = r_new, z_new res.append(np.linalg.norm(r)) return x, res

常用预处理子选择：

• Jacobi预处理：M=diag(A)⁻¹
• 不完全Cholesky分解
• 多项式预处理：M≈(diag(A))⁻¹

4.2 实用终止策略

除残差范数外，推荐组合使用以下判据：

1. 相对残差：∥rₖ∥/∥b∥ < ε
2. 迭代停滞：|∥rₖ∥-∥rₖ₋₁∥| < δ∥rₖ∥
3. 最大迭代次数：k ≥ max_iter

实现示例：

def adaptive_stop(r, r_prev, b, k, tol, max_iter): rel_res = np.linalg.norm(r) / np.linalg.norm(b) stagnate = abs(np.linalg.norm(r)-np.linalg.norm(r_prev)) < 1e-3*np.linalg.norm(r) return rel_res < tol or stagnate or k >= max_iter

4.3 稀疏矩阵优化

对于稀疏矩阵，使用专门的存储格式可大幅提升性能：

from scipy.sparse import csr_matrix def sparse_cg(A_csr, b, x0=None): """针对稀疏矩阵优化的CG实现""" x = np.zeros_like(b) if x0 is None else x0 r = b - A_csr.dot(x) p = r.copy() # 其余部分与普通CG相同，注意用dot代替@

内存占用对比（n=10,000的带状矩阵）：

5. 扩展应用：非线性优化问题

虽然CG专为线性系统设计，但其思想可推广到非线性优化。以Fletcher-Reeves方法为例：

def nonlinear_cg(f, grad_f, x0, max_iter=100): x = x0.copy() grad = grad_f(x) p = -grad res = [np.linalg.norm(grad)] for _ in range(max_iter): # 线搜索求步长（需满足Wolfe条件） alpha = line_search(f, grad_f, x, p) x_new = x + alpha * p grad_new = grad_f(x_new) # Fletcher-Reeves公式 beta = (grad_new @ grad_new) / (grad @ grad) p = -grad_new + beta * p x, grad = x_new, grad_new res.append(np.linalg.norm(grad)) return x, res

实际应用中建议使用更鲁棒的Polak-Ribière变体：

beta_pr = max(0, (grad_new @ (grad_new - grad)) / (grad @ grad))

典型应用场景：

• 机器学习模型训练
• 物理系统能量最小化
• 金融资产组合优化