企业官网建设流程全解析

用Python从零实现共轭梯度法：代码拆解与可视化实战

在机器学习与科学计算的领域中，线性方程组求解是一个基础但至关重要的问题。当面对大规模稀疏矩阵时，传统直接解法往往效率低下，而迭代法中的共轭梯度法（Conjugate Gradient, CG）因其优异性能成为工程师的首选工具。本文将以Python为实践工具，从零开始构建CG算法，通过可视化手段揭示其内在数学之美。

1. 环境准备与问题建模

共轭梯度法要求系数矩阵A对称正定。我们首先生成一个满足条件的测试矩阵：

import numpy as np from scipy.sparse import diags # 生成对称正定矩阵 def generate_spd_matrix(n, kappa=10): """生成n×n对称正定矩阵，条件数为kappa""" A = diags([1, -0.25, -0.25], [0, 1, -1], shape=(n, n)).toarray() A = (A + A.T)/2 # 确保对称 # 调整条件数 eigvals = np.linspace(1, kappa, n) Q, _ = np.linalg.qr(np.random.randn(n,n)) return Q @ np.diag(eigvals) @ Q.T n = 50 A = generate_spd_matrix(n) b = np.random.randn(n)

关键参数说明：

• n: 矩阵维度
• kappa: 条件数控制
• A: 对称正定矩阵
• b: 随机生成的右侧向量

2. 算法核心实现

CG算法的精妙之处在于其共轭方向的构造方式。下面我们分步骤实现算法主体：

def conjugate_gradient(A, b, x0=None, max_iter=1000, tol=1e-6): if x0 is None: x0 = np.zeros_like(b) x = x0.copy() r = b - A @ x p = r.copy() rs_old = r.dot(r) residuals = [] solutions = [x.copy()] for i in range(max_iter): Ap = A @ p alpha = rs_old / p.dot(Ap) x += alpha * p r -= alpha * Ap rs_new = r.dot(r) residuals.append(np.sqrt(rs_new)) solutions.append(x.copy()) if np.sqrt(rs_new) < tol: break beta = rs_new / rs_old p = r + beta * p rs_old = rs_new return x, residuals, solutions

算法关键变量说明：

3. 收敛分析与可视化

理解CG的收敛行为对于实际应用至关重要。我们通过可视化手段展示其特性：

import matplotlib.pyplot as plt # 运行CG算法 x_cg, residuals, solutions = conjugate_gradient(A, b) # 绘制残差下降曲线 plt.figure(figsize=(10,6)) plt.semilogy(residuals, 'o-', label='CG Residual') plt.xlabel('Iteration') plt.ylabel('Residual Norm') plt.title('Convergence History') plt.grid(True) # 添加理论收敛边界 kappa = np.linalg.cond(A) theoretical_bound = [2*( (np.sqrt(kappa)-1)/(np.sqrt(kappa)+1) )**k for k in range(len(residuals))] plt.semilogy(theoretical_bound, '--', label='Theoretical Bound') plt.legend()

典型输出结果分析：

1. 实际收敛曲线通常位于理论边界下方
2. 前几次迭代往往能大幅降低残差
3. 条件数越大，收敛速度越慢

注意：实际应用中，残差可能不会严格单调下降，这与浮点运算精度有关

4. 性能对比与工程实践

我们将自实现的CG与标准库函数进行对比：

from scipy.sparse.linalg import cg as scipy_cg # 使用相同初始值比较 x0 = np.random.randn(n) x_sp, info = scipy_cg(A, b, x0=x0, tol=1e-6) x_our, _, _ = conjugate_gradient(A, b, x0=x0, tol=1e-6) # 验证解的正确性 exact_solution = np.linalg.solve(A, b) print(f"Scipy CG error: {np.linalg.norm(x_sp - exact_solution):.2e}") print(f"Our CG error: {np.linalg.norm(x_our - exact_solution):.2e}")

工程实践建议：

1. 预处理技术能显著改善条件数
2. 对于非对称系统，考虑使用GMRES等算法
3. 实际应用中应监控残差而非依赖固定迭代次数

5. 算法变体与扩展应用

CG算法有多种变体适用于不同场景：

以非线性CG为例，其Python实现核心差异在于：

def nonlinear_conjugate_gradient(f, grad_f, x0, max_iter=100): x = x0.copy() grad = grad_f(x) p = -grad for _ in range(max_iter): # 使用线搜索确定步长 alpha = line_search(f, grad_f, x, p) x_new = x + alpha * p grad_new = grad_f(x_new) # PR+方法计算beta beta = max( (grad_new.dot(grad_new - grad)) / grad.dot(grad), 0 ) p = -grad_new + beta * p x, grad = x_new, grad_new return x

在机器学习中的应用场景：

• 逻辑回归参数优化
• 神经网络权重训练
• 支持向量机求解