共轭梯度法

共轭梯度法

小象实战讲义 · 人工智能数学基础

在求解大规模线性方程组和优化问题时，传统的高斯消元或矩阵分解方法因立方级的计算复杂度而变得低效。本节介绍的共轭梯度法，作为一种迭代算法，为解决这类问题提供了强有力的工具。学完本节，你将理解共轭梯度法的核心思想，掌握其从共轭方向法推导而来的过程，并能够将其应用于求解二次函数优化问题，了解其在处理高维数据时的优势。

💡 核心导读

• 从优化到方程：理解求解线性方程组 $A\mathbf{x} = \mathbf{b}$ 等价于求解二次函数 $\phi(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x}$ 的最小值问题。
• 共轭的定义：掌握向量关于正定矩阵 $A$ 共轭的概念，它是正交概念在广义内积下的推广。
• 共轭方向法：学习如何利用一组 $A$-共轭的方向向量，在最多 $n$ 步内精确求解 $n$ 维二次优化问题。
• 共轭梯度法：掌握如何不预先给定共轭方向，而是利用梯度和前一步方向动态生成共轭方向，形成高效的迭代算法。
• 推广与应用：了解共轭梯度法如何推广到一般非线性函数的优化，以及其在实际应用中的注意事项。

从线性方程组到二次优化

我们首先回顾求解线性方程组 $A\mathbf{x} = \mathbf{b}$ 的经典方法。若 $A$ 是一般方阵，可采用 LU 分解；若 $A$ 是正定矩阵，则可采用更高效的 Cholesky 分解。这两种方法的计算复杂度均为 $O(n^3)$，其中 $n$ 是矩阵的维数。

那么，线性方程组与优化问题有何关联？考虑如下二次函数： $$\phi(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x}$$ 其中 $A \in \mathbb{R}^{n \times n}$ 是对称正定矩阵。该函数的最小值点满足一阶导数为零： $$\nabla \phi(\mathbf{x}) = A\mathbf{x} - \mathbf{b} = 0$$ 这正是线性方程组 $A\mathbf{x} = \mathbf{b}$。因此，求解该线性方程组等价于求解二次函数 $\phi(\mathbf{x})$ 的无约束极小化问题。

对于迭代法，我们从初始点 $\mathbf{x}_0$ 出发，生成序列 ${\mathbf{x}$k}{xk​}，希望其收敛到解 $\mathbf{x}$。定义第 $k$ 步的余量（残差）为： $$\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k = -\nabla \phi(\mathbf{x}_k)$$ 显然，当 $\mathbf{x}$k \to \mathbf{x} 时，有 $\mathbf{r}_k \to \mathbf{0}$。迭代法的核心在于设计 $\mathbf{x}_k$ 的更新规则，使余量快速收敛。

共轭方向：广义的正交

共轭梯度法的关键在于“共轭”二字。我们首先定义向量的共轭性。

定义（$A$-共轭）：设 $A$ 是一个 $n \times n$ 的对称正定矩阵。一组非零向量 ${\mathbf{p}_0, \mathbf{p}$1, \dots, \mathbf{p}{\ell}} 被称为关于 $A$ 共轭，如果它们满足： $$\mathbf{p}_i^T A \mathbf{p}_j = 0, \quad \forall i \ne j$$

可以证明，$A$-共轭的向量组必然是线性无关的。

如何理解“共轭”？考虑最简单的情况：$A = I$（单位矩阵）。此时共轭条件变为 $\mathbf{p}_i^T \mathbf{p}_j = 0$，这正是向量两两正交的定义。因此，当 $A=I$ 时，共轭即正交。

对于一般的正定矩阵 $A$，共轭可以视为在由 $A$ 定义的“广义内积” $<\mathbf{u}, \mathbf{v}>_A = \mathbf{u}^T A \mathbf{v}$ 下的正交。另一种直观理解是，$A\mathbf{p}_j$ 相当于对向量 $\mathbf{p}_j$ 进行了一次线性变换（旋转和拉伸），共轭条件 $\mathbf{p}_i^T (A\mathbf{p}_j)=0$ 意味着变换后的 $\mathbf{p}_j$ 与原始的 $\mathbf{p}_i$ 正交。因此，共轭是正交概念在更一般度量下的推广。

共轭方向法及其有限步收敛性

基于共轭方向，我们可以构造一种强大的迭代算法——共轭方向法。

算法框架（共轭方向法）： 给定初始点 $\mathbf{x}_0$ 和一组关于 $A$ 共轭的非零方向向量 ${\mathbf{p}_0, \mathbf{p}$1, \dots, \mathbf{p}{n-1}}，按如下方式迭代：

1. 计算当前余量：$\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k$。
2. 确定步长 $\alpha_k$，使得沿方向 $\mathbf{p}$kpk​ 前进能最小化目标函数： $$\alpha_k = \arg\min$${\alpha} \phi(\mathbf{x}_k + \alpha \mathbf{p}_k) 通过求解 $\frac{d}{d\alpha} \phi(\mathbf{x}_k + \alpha \mathbf{p}_k) = 0$，可得： $$\alpha_k = \frac{\mathbf{p}_k^T \mathbf{r}_k}{\mathbf{p}_k^T A \mathbf{p}_k}$$
3. 更新迭代点：$\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k$。

关键定理：对于上述共轭方向法，序列 ${\mathbf{x}$k}{xk​} 将在最多 $n$ 步内收敛到线性方程组 $A\mathbf{x} = \mathbf{b}$ 的解 $\mathbf{x}$。即，$\mathbf{x}$n = \mathbf{x}。

这个性质非常美妙：对于一个 $n$ 维问题，算法保证在 $n$ 步内得到精确解（在无舍入误差的理想情况下）。其几何解释是，在由共轭方向张成的空间中，算法每一步都在一个独立的“广义坐标轴”方向上找到该方向上的最优解，从而一步步精确重构出全局最优解。

为了直观理解 $n$ 步收敛，考虑两种情况：

1. $A$ 为对角矩阵：此时目标函数的等高线是轴对齐的椭圆。从任意点 $\mathbf{x}_0$ 出发，先沿 $\mathbf{p}_0$（例如坐标轴方向）走到该方向上的最优点 $\mathbf{x}_1$，再沿 $\mathbf{p}$1p1​（另一坐标轴方向）即可一步到达中心 $\mathbf{x}$*。每一步解决一个维度。
2. $A$ 为一般正定矩阵：正定矩阵 $A$ 可通过正交相似变换对角化，即存在正交矩阵 $S$ 使得 $S^T A S = D$（对角阵）。令 $\hat{\mathbf{x}} = S^T \mathbf{x}$，则原问题 $\min \phi(\mathbf{x})$ 在新变量 $\hat{\mathbf{x}}$ 下变为一个对角矩阵的优化问题，回归到情况1。这个变量替换 $S^T$ 相当于对坐标轴进行了一次旋转，使得新的坐标轴方向（即 $S$ 的列向量）关于 $A$ 共轭。算法正是在这些旋转后的“广义坐标轴”上依次寻优。

下面的代码演示了对于一个二维对角矩阵的二次函数，使用其坐标轴方向作为共轭方向，两步即可收敛到最优解。

import numpy as np

# 定义对角正定矩阵A和向量b
A = np.array([[4.0, 0.0],
              [0.0, 1.0]])
b = np.array([0.0, 0.0])
x_star = np.array([0.0, 0.0])  # 真实解

# 选择共轭方向（这里就是坐标轴方向，因为A是对角阵）
p0 = np.array([1.0, 0.0])  # 关于A共轭的方向
p1 = np.array([0.0, 1.0])  # 关于A共轭的方向

# 初始点
x0 = np.array([2.0, 3.0])
x = x0.copy()

print(f"初始点 x0: {x}")
print("="*40)

# 第一步，沿p0方向
r = b - A @ x
alpha = (p0 @ r) / (p0 @ A @ p0)
x = x + alpha * p0
print(f"第一步后 x1: {x}, 余量 r1: {b - A @ x}")

# 第二步，沿p1方向
r = b - A @ x
alpha = (p1 @ r) / (p1 @ A @ p1)
x = x + alpha * p1
print(f"第二步后 x2: {x}, 余量 r2: {b - A @ x}")
print(f"是否等于精确解 x*? {np.allclose(x, x_star)}")

共轭梯度法：动态生成共轭方向

共轭方向法需要一个预先给定的共轭向量组。共轭梯度法的巧妙之处在于，它能在迭代过程中动态地生成这些共轭方向，而无需预先计算和存储。

其核心思想是：利用当前点的负梯度方向（即余量 $\mathbf{r}$krk​）和前一步的方向 $\mathbf{p}${k-1}，通过一个线性组合来构造新的方向 $\mathbf{p}_k$，并要求 $\mathbf{p}$kpk​ 与 $\mathbf{p}${k-1} 关于 $A$ 共轭。

设方向更新公式为： $$\mathbf{p}$$k = -\mathbf{r}k + \beta_k \mathbf{p}{k-1}pk​=−rk​+βk​pk−1​ 其中 $\beta_k$ 是待定系数。将共轭条件 $\mathbf{p}$k^T A \mathbf{p}{k-1} = 0 代入，可以解出： $$\beta_k = \frac{\mathbf{r}$$k^T A \mathbf{p}{k-1}}{\mathbf{p}{k-1}^T A \mathbf{p}_{k-1}} 可以证明，按此方式生成的方向向量组 ${\mathbf{p}_0, \mathbf{p}_1, \dots}$ 确实是关于 $A$ 共轭的。

结合步长公式 $\alpha_k = \frac{\mathbf{p}_k^T \mathbf{r}_k}{\mathbf{p}_k^T A \mathbf{p}_k}$，我们得到标准的线性共轭梯度法算法。

算法（线性共轭梯度法）：

1. 初始化：给定 $\mathbf{x}_0$，计算 $\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0$，令 $\mathbf{p}_0 = \mathbf{r}_0$。
2. 对 $k=0,1,2,\dots$ 直到收敛，执行： a. 计算步长：$\alpha_k = \frac{\mathbf{r}_k^T \mathbf{r}_k}{\mathbf{p}$k^T A \mathbf{p}k}αk​=pkT​Apk​rkT​rk​​。 b. 更新解：$\mathbf{x}${k+1} = \mathbf{x}k + \alpha_k \mathbf{p}k。 c. 更新余量：$\mathbf{r}${k+1} = \mathbf{r}k - \alpha_k A \mathbf{p}k。 d. 计算系数：$\beta${k+1} = \frac{\mathbf{r}{k+1}^T \mathbf{r}{k+1}}{\mathbf{r}k^T \mathbf{r}k}。 e. 更新方向：$\mathbf{p}${k+1} = \mathbf{r}{k+1} + \beta{k+1} \mathbf{p}_k。

注意：上述第 a 步和第 d 步的公式是经过数学化简后的等价形式，它们比原始定义式在计算上更高效，且避免了显式存储方向向量 $\mathbf{p}_k$ 进行某些内积运算，节省了存储空间。这正是算法的“实用版本”。

共轭梯度法继承了共轭方向法的有限步收敛性质：对于 $n$ 维问题，理论上最多 $n$ 步即可得到精确解。此外，还有一个更强的结论：若 $A$ 有 $r$ 个互不相同的特征值，则算法最多 $r$ 步收敛。

推广至一般非线性优化

共轭梯度法最初用于求解线性方程组（二次优化），但可以自然地推广到求解一般非线性函数 $f(\mathbf{x})$ 的无约束优化问题。

思路是在当前迭代点 $\mathbf{x}_k$ 处，用二阶泰勒展开来局部近似目标函数： $$f(\mathbf{x}) \approx f(\mathbf{x}_k) + \nabla f(\mathbf{x}_k)^T (\mathbf{x} - \mathbf{x}_k) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_k)^T H_f(\mathbf{x}_k) (\mathbf{x} - \mathbf{x}_k)$$ 其中 $H_f(\mathbf{x}_k)$ 是 $f$ 在 $\mathbf{x}_k$ 处的 Hessian 矩阵。这个近似函数是一个二次函数。如果我们将其视为新的目标，那么搜索方向 $\mathbf{p}_k$ 的构造可以沿用线性共轭梯度法的思想，即： $$\mathbf{p}_k = -\nabla f(\mathbf{x}$$k) + \beta_k \mathbf{p}{k-1} 系数 $\beta_k$ 有不同的计算公式，常见的有 Fletcher-Reeves (FR) 公式和 Polak-Ribière (PR) 公式等，例如 FR 公式为： $$\beta_k^{\text{FR}} = \frac{\nabla f(\mathbf{x}$$k)^T \nabla f(\mathbf{x}k)}{\nabla f(\mathbf{x}{k-1})^T \nabla f(\mathbf{x}{k-1})}

与线性情形的关键区别在于步长 $\alpha_k$ 的选取。对于非线性函数，我们无法像二次函数那样通过解析式求出精确最小化步长。因此，需要使用非精确线搜索方法来确定 $\alpha_k$，例如满足 Wolfe 条件的搜索或回溯法。

算法框架（非线性共轭梯度法）：

1. 初始化 $\mathbf{x}_0$，计算 $\mathbf{g}_0 = \nabla f(\mathbf{x}_0)$，令 $\mathbf{p}_0 = -\mathbf{g}_0$。
2. 对 $k=0,1,2,\dots$ 直到收敛，执行： a. 使用非精确线搜索求步长 $\alpha_k$，使得 $f(\mathbf{x}$k + \alpha_k \mathbf{p}k)f(xk​+αk​pk​) 充分下降。 b. 更新解：$\mathbf{x}${k+1} = \mathbf{x}k + \alpha_k \mathbf{p}k。 c. 计算新梯度：$\mathbf{g}${k+1} = \nabla f(\mathbf{x}{k+1})。 d. 计算系数 $\beta${k+1}（如用 FR 公式）。 e. 更新方向：$\mathbf{p}${k+1} = -\mathbf{g}{k+1} + \beta_{k+1} \mathbf{p}kpk+1​=−gk+1​+βk+1​pk​。 f. （可选）如果满足重置条件（如 $k$ 是 $n$ 的倍数或 $|\mathbf{g}${k+1}^T \mathbf{g}k| 很大），则令 $\mathbf{p}${k+1} = -\mathbf{g}_{k+1}（即重置为最速下降方向）。

📝 动手练一练

1. 共轭方向验证：给定正定矩阵 $A = \begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix}$，验证向量 $\mathbf{p}_0 = \begin{bmatrix} 1 \ 0 \end{bmatrix}$ 和 $\mathbf{p}_1 = \begin{bmatrix} -1/2 \ 1 \end{bmatrix}$ 是否关于 $A$ 共轭。并尝试求解线性方程组 $A\mathbf{x} = \mathbf{b}$，其中 $\mathbf{b} = \begin{bmatrix} 3 \ 3 \end{bmatrix}$，使用 $\mathbf{x}_0 = \begin{bmatrix} 0 \ 0 \end{bmatrix}$ 和这两个共轭方向，验证是否两步即可得到精确解。

   参考答案：

   import numpy as np
   A = np.array([[2.0, 1.0], [1.0, 2.0]])
   b = np.array([3.0, 3.0])
   p0 = np.array([1.0, 0.0])
   p1 = np.array([-0.5, 1.0])
   
   # 验证共轭性
   print("p0^T A p1 =", p0 @ A @ p1)  # 应为0
   
   # 共轭方向法求解
   x = np.array([0.0, 0.0])
   # 第一步
   r = b - A @ x
   alpha0 = (p0 @ r) / (p0 @ A @ p0)
   x = x + alpha0 * p0
   print(f"x1 = {x}")
   # 第二步
   r = b - A @ x
   alpha1 = (p1 @ r) / (p1 @ A @ p1)
   x = x + alpha1 * p1
   print(f"x2 = {x}")
   print(f"验证 A*x2 = b? {np.allclose(A @ x, b)}")

2. 非线性共轭梯度法尝试：对于二元函数 $f(x, y) = x^2 + 10y^2$（这是一个简单的二次函数，其 Hessian 矩阵为常数），尝试实现使用 FR 公式和非精确线搜索（简单的回溯法）的非线性共轭梯度法进行优化。设置初始点为 $(5, 5)$，观察其收敛过程。

   参考答案（回溯法线搜索）：

   import numpy as np
   
   def f(x):
       return x[0]**2 + 10 * x[1]**2
   
   def grad_f(x):
       return np.array([2*x[0], 20*x[1]])
   
   def backtracking_line_search(f, grad, x, p, alpha_init=1.0, rho=0.5, c=1e-4):
       """回溯法线搜索"""
       alpha = alpha_init
       while f(x + alpha * p) > f(x) + c * alpha * grad(x) @ p:
           alpha *= rho
       return alpha
   
   # 非线性共轭梯度法 (FR)
   x = np.array([5.0, 5.0])
   g = grad_f(x)
   p = -g
   k = 0
   max_iter = 50
   tol = 1e-6
   
   for k in range(max_iter):
       # 线搜索
       alpha = backtracking_line_search(f, grad_f, x, p)
       # 更新点
       x_new = x + alpha * p
       g_new = grad_f(x_new)
       
       # 检查收敛
       if np.linalg.norm(g_new) < tol:
           print(f"在 {k+1} 步后收敛")
           break
           
       # FR公式计算beta
       beta = (g_new @ g_new) / (g @ g)
       # 更新方向
       p = -g_new + beta * p
       # 为下一次迭代准备
       x, g = x_new, g_new
   
   print(f"最终解: {x}, 最终梯度范数: {np.linalg.norm(g)}")

本章小结

本节深入探讨了共轭梯度法这一高效的最优化算法。

要点回顾：

1. 问题转化：求解对称正定线性方程组 $A\mathbf{x}=\mathbf{b}$ 等价于最小化二次函数 $\phi(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x}$。
2. 共轭性：关于正定矩阵 $A$ 的共轭是正交概念的推广，$A$-共轭的向量组线性无关。
3. 共轭方向法：利用一组预先给定的 $A$-共轭方向，可以在最多 $n$ 步内精确求解 $n$ 维二次优化问题，其几何意义是在旋转后的广义坐标轴上依次寻优。
4. 共轭梯度法：通过将当前负梯度方向与前一步方向进行线性组合，并施加共轭条件，可以动态生成共轭方向，形成高效的迭代算法。其核心是 $\alpha_k$ 和 $\beta_k$ 的更新公式。
5. 推广与应用：通过局部二次近似，共轭梯度法可推广至一般非线性优化，此时需配合非精确线搜索。该算法特别适合大规模、稀疏问题，因为它具有有限步收敛的理论保证和节省存储空间的优点。

行动清单：

• 理论验证：任选一个二维正定矩阵 $A$，手动推导两个关于 $A$ 共轭的向量，并用代码验证其共轭性及共轭方向法的两步收敛性。
• 算法实现：使用 Python 实现线性共轭梯度法，并用于求解一个至少 100 维的随机生成的正定线性方程组，观察其迭代收敛过程。
• 对比思考：思考共轭梯度法与最速下降法、牛顿法在思想、收敛速度、计算开销上的主要区别，总结其各自的适用场景。

— 小象教研组