讲透一个强大算法模型，共轭梯度法

最新推荐文章于 2025-04-25 21:12:48 发布

我不是小upper

最新推荐文章于 2025-04-25 21:12:48 发布

阅读量988

点赞数 23

分类专栏：数据科学机器学习文章标签：算法机器学习人工智能

本文链接：https://blog.csdn.net/qq_70350287/article/details/147447466

版权

数据科学同时被 2 个专栏收录

29 篇文章

订阅专栏

机器学习

22 篇文章

订阅专栏

哈喽，大家好，我是不是小upper，

今天，咱们来聊聊优化算法里的共轭梯度法。在较为简单的情境下，我们面临这样一个问题：有一个 “碗状” 函数（从某种角度看，这等同于求解线性方程组）

我们的目标是找出能使该函数取值最小（或者让线性方程 $Ax = b$ 组成立）的向量。

1. 为什么最速下降法不够「聪明」？

最速下降法每次沿着当前梯度方向（最陡下坡方向）更新参数，看似高效，却存在致命缺陷：当目标函数呈「扁碗状」（即矩阵 A 的特征值差异大，几何上表现为某些方向极陡、另一些方向极缓），算法会在不同主轴方向间来回震荡，走出「之字形」路径，收敛极慢。比如在狭长的山谷中，最速下降法像一个只盯着脚下最陡下坡的行者，不断在两侧山壁间折返，难以直抵谷底。这种「短视」行为在高维稀疏问题中尤为明显，当变量间存在强耦合（矩阵 A 条件数大）时，收敛速度会随维度增加急剧恶化。

2. 共轭梯度法：让搜索方向「互不干扰」的聪明策略

核心思想：用「共轭方向」打破震荡

共轭梯度法的关键是构造一组 「共轭方向」—— 即每一步的搜索方向与之前所有方向在矩阵 A 的意义下「正交」（数学上称为 A- 共轭）。

数学定义：对于对称正定矩阵 A，向量 $d_i$ 和 $d_j$ 若满足 $d_i^T A d_j = 0$ $(i \neq j)$ ，则称它们是 A- 共轭的。
直观理解：这些方向在「加权空间」中彼此独立，沿其中一个方向优化后，后续方向的优化不会「undo」之前的成果。就像在扁碗中找到一系列互不干扰的下坡路径，每一步都彻底清除当前方向的残差，直指碗底。

优势：

避免震荡：通过共轭性，每次搜索方向都「继承」之前的优化成果，不会在已优化的方向上重复做功，收敛路径更直接。
线性收敛加速：对于 n 维问题，理论上至多 n 步即可精确收敛（针对二次型问题），尤其适合大规模、稀疏的线性系统（如 A 为稀疏矩阵时，计算效率远超最速下降）。

原理详解：从二次型优化到共轭方向构造

1. 问题等价：最小化二次型 = 解线性方程组

目标：最小化二次型函数 $f(x) = \frac{1}{2}x^T A x - b^T x$ （A 对称正定）。

求导得梯度 $g(x) = Ax - b$ ，令梯度为零即得线性方程组 $Ax = b$ 。
因此，找 $f(x)$ 的最小值等价于求解 $Ax = b$ 。

2. 最速下降法的瓶颈：特征值「跨度」的诅咒

最速下降法迭代公式：

$x_{k+1} = x_k - \alpha_k g_k, \quad \alpha_k = \frac{g_k^T g_k}{g_k^T A g_k}$

当 A 的特征值 $\lambda_1 \gg \lambda_n$ （条件数 $\kappa = \lambda_1/\lambda_n$ 大），迭代次数随 $\kappa$ 增长呈线性恶化，高维下几乎不可用。

3. 共轭方向：清扫残差的「正交工具」

共轭方向性质：若方向组 $\{d_0, d_1, \dots, d_{k}\}$ 两两 A- 共轭，则沿每个方向 $d_i$ 做一维精确优化（即找到最优步长）后，残差 $r_{i+1} = b - A x_{i+1}$ 会与所有已用方向正交 $d_j^T r_{i+1} = 0$ , $j \leq i$
有限步收敛：对于 n 维问题，至多 n 步即可找到精确解（忽略计算误差），因为 n 个共轭方向张成整个空间，残差可被完全清除。

4. 迭代步骤（针对二次型问题）：

初始化： $x_0$ 任选，计算初始残差 $r_0 = b - A x_0$ ，取初始方向 $d_0 = r_0$ 。
迭代：对：
- 计算步长 $\alpha_k = \frac{r_k^T r_k}{d_k^T A d_k}$ ，更新 $x_{k+1} = x_k + \alpha_k d_k$ 。
- 计算新残差 $r_{k+1} = r_k - \alpha_k A d_k$ 。
- 构造新方向 $d_{k+1} = r_{k+1} + \beta_k d_k$ ，其中 $\beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$ （Fletcher-Reeves 公式），确保 $d_{k+1}$ 与 $d_k$ 共轭。

5. 核心优势：

无需矩阵求逆：仅通过向量运算和矩阵 - 向量乘法迭代，适合稀疏矩阵。
内存友好：每步仅需存储当前和前一步的残差与方向，空间复杂度低。

总结：共轭梯度法的「聪明之处」

共轭梯度法通过构造 A- 共轭的搜索方向，将最速下降法的「短视震荡」转化为「全局协同」，在二次型优化和线性方程组求解中实现高效收敛。尤其在高维稀疏场景（如图像处理、有限元分析）中，其性能远超最速下降法，成为大规模优化问题的核心工具之一。

完整案例

我们以二维 Poisson 方程离散化生成的稀疏矩阵作为测试对象，该矩阵广泛用于模拟有限元 / 有限差分方法中常见的稀疏正定线性系统。具体而言，考虑一个 $n \times n$ 的网格，对应的线性系统矩阵维度为 $N = n^2$ 。以 $n = 100$ 为例，此时矩阵维度为 $10000 \times 10000$ ，具有典型的稀疏结构。

矩阵生成函数 `generate_poisson_2d` 的构造逻辑：

离散化模型：采用有限差分法对二维 Poisson 方程进行离散，内部节点的差分格式对应五点模板：
- 主对角元素：每个节点的中心系数为 4（对应离散化后节点自身的权重）。
- 邻居连接权重：上下左右四个相邻节点的连接权重均为 -1（对应二阶空间导数的中心差分近似）。
边界处理：
- 严格遵循网格的物理边界，避免行尾节点与下一行首节点连接（即不引入跨边界的非物理连接），确保矩阵结构符合有限差分 / 有限元方法的离散规则，维持稀疏矩阵的正确性与对称性。

该构造方法生成的矩阵具有对称正定特性，是验证稀疏矩阵求解算法（如共轭梯度法）的经典测试算例，其稀疏结构（每行约 5 个非零元素）反映了二维空间中节点的局部连接特性。

一、共轭梯度法（CG）与预条件共轭梯度法（PCG）原理

1. 标准共轭梯度法（CG）

核心思想：针对对称正定线性系统 $A\mathbf{x} = \mathbf{b}$ ，通过构造共轭方向序列逐步逼近精确解，避免直接计算矩阵逆。

初始条件：取初始猜测 $\mathbf{x}_0$ ，计算初始残差 $\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0$ ，初始搜索方向 $\mathbf{p}_0 = \mathbf{r}_0$ 。
迭代更新： $\alpha_k = \frac{\mathbf{r}_k^\top \mathbf{r}_k}{\mathbf{p}_k^\top A\mathbf{p}_k}, \quad \mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k,\)\(\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k A\mathbf{p}_k, \quad \beta_k = \frac{\mathbf{r}_{k+1}^\top \mathbf{r}_{k+1}}{\mathbf{r}_k^\top \mathbf{r}_k}, \quad \mathbf{p}_{k+1} = \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k,$

至残差范数 $\|\mathbf{r}_k\| \leq \text{tol}$ 终止。

2. Jacobi 预条件共轭梯度法（PCG）

预条件技术：通过预条件器 $M \approx A$ 改善矩阵条件数，压缩特征值分布，加速收敛。

Jacobi 预条件器：取 $M = \text{diag}(A)$ ，即对角矩阵，其逆为 $M^{-1}\mathbf{r} = \mathbf{r} / \text{diag}(A)$ 。
预条件化残差：将残差映射为 $\mathbf{z}_k = M^{-1}\mathbf{r}_k$ ，并在迭代中用 $\mathbf{z}_k$ 替代原残差计算搜索方向：

$\alpha_k = \frac{\mathbf{r}_k^\top \mathbf{z}_k}{\mathbf{p}_k^\top A\mathbf{p}_k}, \quad \beta_k = \frac{\mathbf{r}_{k+1}^\top \mathbf{z}_{k+1}}{\mathbf{r}_k^\top \mathbf{z}_k}.$

二、算法实现与测试问题

1. 2D Poisson 问题建模

生成 $n \times n$ 网格上的稀疏对称正定矩阵 A，对应离散化的 Poisson 方程： $A = \begin{bmatrix} 4I & -I & & \\ -I & 4I & -I & \\ & \ddots & \ddots & -I \\ & & -I & 4I \end{bmatrix} \in \mathbb{R}^{n^2 \times n^2},$

其中 I 为 $n \times n$ 单位矩阵，非零元素描述网格点与上下左右邻居的耦合关系。

2. 算法实现要点

残差历史记录：在迭代中跟踪 $\|\mathbf{r}_k\|$ ，用于分析收敛行为。
预条件器集成：通过函数 $M^{-1}(\mathbf{r}) = \mathbf{r} / \text{diag}(A)$ 实现 Jacobi 预条件，无缝嵌入 CG 框架。

Jacobi 预条件共轭梯度法（PCG）的代码实现

预条件器设计：Jacobi 预条件器利用矩阵对角线元素构造对角矩阵 $M = \text{diag}(A)$，其逆运算为元素级除法，计算复杂度极低，适合稀疏矩阵场景。

# Jacobi 预条件器：M ≈ diag(A)
diag_A = A.diagonal()
M_inv = lambda r: r / diag_A  # 元素级除法实现 M^{-1}r

3. 共轭梯度法核心实现（支持预条件）

以下代码实现了标准 CG 与预条件 CG 的统一框架，通过参数 M_inv 控制预条件器的启用：

def conjugate_gradient(A, b, x0=None, tol=1e-8, maxiter=1000, M_inv=None):
    """
    共轭梯度法（支持预条件）
    参数:
        A: 稀疏矩阵 (scipy.sparse.csr_matrix)
        b: 右端项向量 (numpy.ndarray)
        M_inv: 预条件器函数，None表示无预条件
    返回:
        x: 近似解, residuals: 残差范数历史
    """
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()
    r = b - A.dot(x)  # 初始残差
    z = M_inv(r) if M_inv is not None else r.copy()  # 预条件残差
    p = z.copy()       # 初始搜索方向
    
    residuals = [np.linalg.norm(r)]  # 记录残差历史
    
    for k in range(maxiter):
        Ap = A.dot(p)          # 矩阵-向量乘法
        alpha = (r @ z) / (p @ Ap)  # 步长计算（利用预条件残差）
        x += alpha * p         # 解更新
        r_next = r - alpha * Ap  # 新残差
        
        # 提前终止检查
        res_norm = np.linalg.norm(r_next)
        residuals.append(res_norm)
        if res_norm < tol:
            break
        
        # 预条件残差更新与搜索方向构造
        z_next = M_inv(r_next) if M_inv is not None else r_next
        beta = (r_next @ z_next) / (r @ z)  # 共轭系数计算
        p = z_next + beta * p              # 新搜索方向（共轭方向）
        
        r, z = r_next, z_next  # 状态更新
    
    return x, residuals

4. 测试问题：2D Poisson 方程的离散化矩阵

通过有限差分法生成 $n \times n$ 网格上的稀疏对称正定矩阵，对应离散 Poisson 方程 $-\nabla^2 u = 1$ ，边界条件为 Dirichlet 零边界：

def generate_poisson_2d(n):
    """生成n x n网格的2D Poisson矩阵（CSR格式）"""
    N = n * n
    diag = 4 * np.ones(N)
    off1 = -1 * np.ones(N - 1)       # 水平邻居
    offn = -1 * np.ones(N - n)       # 垂直邻居
    # 构造稀疏矩阵（主对角、左右、上下）
    A = sp.diags([diag, off1, off1, offn, offn], [0, -1, 1, -n, n], shape=(N, N), format='csr')
    # 修复边界连接（避免行末连接到下一行开头）
    for i in range(1, n):
        A[i*n, i*n - 1] = 0  # 移除无效的左边界连接
    return A

# 生成50x50网格的测试矩阵
n = 50
A = generate_poisson_2d(n)
N = n * n
b = np.ones(N)  # 右端项为全1向量

5. 收敛行为对比实验

通过求解相同右端项 $\mathbf{b} = \mathbf{1}$ ，对比标准 CG 与 Jacobi 预条件 CG 的残差下降曲线：

# 标准CG与PCG求解
x_cg, res_cg = conjugate_gradient(A, b, tol=1e-6, maxiter=500)
x_pcg, res_pcg = conjugate_gradient(A, b, tol=1e-6, maxiter=500, M_inv=M_inv)

# 可视化残差收敛曲线（图1）
plt.figure(figsize=(8, 5))
plt.semilogy(res_cg, 'r-', linewidth=2, label='CG')
plt.semilogy(res_pcg, 'b-', linewidth=2, label='PCG (Jacobi)')
plt.xlabel('Iteration', fontsize=12)
plt.ylabel('Residual Norm (log scale)', fontsize=12)
plt.title('Convergence Comparison of CG and PCG', fontsize=14)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.show()

三、收敛性与谱分析

1. 残差收敛对比

CG（红色曲线）：无预条件时，残差呈几何衰减但速度较慢，需约 200 次迭代降至 $10^{-6}$ 。
PCG（蓝色曲线）：引入 Jacobi 预条件后，残差在前 100 次迭代内快速降至 $10^{-8}$ ，收敛速度显著提升。
关键观察：预条件通过改善矩阵条件数，使搜索方向更接近 “最优共轭方向”，减少迭代中的 “锯齿效应”。

2. 矩阵谱特性

特征值分布：前 300 个最小特征值范围为 $[\lambda_{\text{min}}, \lambda_{\text{max}}] \approx [0, 1.4]$ ，条件数 $\kappa(A) = \lambda_{\text{max}}/\lambda_{\text{min}}$ 较大，系统呈现病态。
预条件作用：理想预条件器可将特征值 “压缩” 至更小区间，使 CG 的理论收敛率 $O((\sqrt{\kappa} - 1)/(\sqrt{\kappa} + 1))^k$ 显著提升。

四、算法优化实验：动量 CG 的探索

1. 动量项引入

受深度学习启发，尝试在搜索方向中加入历史位移动量： $\mathbf{x}_{\text{new}} = \mathbf{x} + \alpha \mathbf{p} + \mu (\mathbf{x} - \mathbf{x}_{\text{prev}}),$ 其中 $\mu = 0.05$ 为动量系数，试图利用历史梯度信息加速收敛。