共轭梯度法

o( ￣⌣￣ )o included in Math

2025-12-25 2039 words 5 minutes

Contents

参考资料

简介

Conjugate Gradient Method, 简称 CG 算法，常用于求解线性方程组

共轭梯度法的思想

选择一个优化方向后, 便将该方向优化到极致，后面优化的过程不再影响之前的方向的极小值
这就要求每一次的优化方向都共轭正交

什么是正交？

假设有一组向量 $S_{1}S_{2}…S_{m}$, 以及一个正定矩阵 $A$
若有

$$ \left{ \begin{matrix} S_{i}^{T} A S_{j} = 0 \ S_{i}^{T} A S_{i} \neq 0 \ \end{matrix} \right. $$

则称 $S_{1}S_{2}…S_{m}$ 是 $A$ 的 共轭向量，它们之间线性无关
若 $A$ 还是单位矩阵，则 $S_{1}S_{2}…S_{m}$ 是 $A$ 的 正交向量

共轭是正交的推广，正交的要求更高

什么是共轭？

假设有矩阵 $X、Y$，矩阵 $A$ 是 $n$ 阶正定阵
若有 $X^{T} A Y = 0$，则称 $X$ 与 $Y$ 关于 $A$ 是共轭的

共轭梯度法的优缺点

优

克服了梯度下降收敛慢
避免了牛顿法所需的二阶导数信息
对于 $n$ 维优化问题，最多 $n$ 次迭代就能找到最优解
缺

极端情况下，矩阵 $A$ 的条件数很大，会影响收敛速度
需要注意设置迭代次数，并不是越多越好

算法过程

假设有方程：$Ax = b$
其中 $A$ 是 $n$ 阶的对称正定矩阵，$b$ 已知，求解 $x$

构造

$$ f(x) = \frac{1}{2}x^{T}Ax - b^{T}x + c $$

则

$$ \frac{\partial f(x)}{\partial x} = Ax - b $$

由 $Ax = b$ 得

$$ \frac{\partial f(x)}{\partial x} = 0 $$

即求解方程转换成求能使得 $f(x)$ 最小的向量

$$ x^{*} = \mathop{\arg\min} \limits_{x} (\frac{1}{2}x^{T}Ax - b^{T}x + c) $$

向量 $x$ 任取初值 $x_{0}$，函数 $f(x)$ 在 $x_{0}$ 处的 负梯度方向 为 $r_{0} = b − Ax_{0}$

最速下降（Gradient Descent）法的基本思想是每次迭代都沿着 $x$ 的梯度的方向下降，最终会收敛到最优值上
共轭梯度法则是在由前一步的下山方向 $p_{k-1}$ 与当前步的梯度方向 $r_{k}$ 所张成的「二维平面」上找到新的最优下山方向

如何推导

$x$、$b$、$r$ 是向量，$c$ 是常数，$\alpha$、$\beta$、$\xi$、$\eta$ 是标量
给定初始向量 $x_{0}$，第一步选取的下山方向 $p_{0}$ 仍是 $r_{0} = b − Ax_{0}$
之后第 $k+1$ 步 $（k \geq 1）$，就在前一步的下山方向 $p_{k−1}$ 和当前步的梯度方向 $r_{k}$ 所张成的平面上找最优的下山方向，即第 $k+1$ 步的下山方向 $p_{k}$ 为：

$$ p_{k} = x_{k} + \xi r_{k} + \eta p_{k-1}，其中 \xi, \eta \in \mathbb{R} $$

当然 $p_{k}$ 也不能任意取，需要考虑到二次型 $f(x)$ 的限制？？？？？
将 $p_{k}$ 带入 $f(x)$，构造一个新函数 $\psi(\xi, \eta)$：

$$ \begin{aligned} \psi(\xi, \eta) = \ & f(x_{k} + \xi r_{k} + \eta p_{k-1})\ = \ & \frac{1}{2}(x_{k} + \xi r_{k} + \eta p_{k-1})^\mathrm{T}A(x_{k} + \xi r_{k} + \eta p_{k-1})\ \ & -b^\mathrm{T}(x_{k} + \xi r_{k} + \eta p_{k-1}) + c \end{aligned} \tag{1} $$

式 1 分别对 $\xi$ 和 $\eta$ 求偏导得

$$ \begin{aligned} \frac{\partial\psi}{\partial\xi} = \ & r_{k}^\mathrm{T}A(x_{k} + \xi r_{k} + \eta p_{k-1}) - b^\mathrm{T}r_{k} \ = \ & r_{k}^\mathrm{T}Ax_{k} + \xi r_{k}^\mathrm{T}Ar_{k} + \eta r_{k}^\mathrm{T}Ap_{k-1} - b^\mathrm{T}r_{k}\ \ & 又因为 r_{k} = b - Ax_{k}\ = \ & r_{k}^\mathrm{T}(b - r_k) + \xi r_{k}^\mathrm{T}Ar_{k} + \eta r_{k}^\mathrm{T}Ap_{k-1} - b^\mathrm{T}r_{k}\ \ & 又因为 r_{k}^\mathrm{T}b = b^\mathrm{T}r_{k}\ = \ & -r_{k}^\mathrm{T}r_k + \xi r_{k}^\mathrm{T}Ar_{k} + \eta r_{k}^\mathrm{T}Ap_{k-1}\ \frac{\partial\psi}{\partial\eta} =\ &p_{k-1}^\mathrm{T}A(x_{k} + \xi r_{k} + \eta p_{k-1}) - b^\mathrm{T}p_{k-1}\ = \ & p_{k-1}^\mathrm{T}Ax_{k} + \xi p_{k-1}^\mathrm{T}Ar_{k} + \eta p_{k-1}^\mathrm{T}Ap_{k-1} - b^\mathrm{T}p_{k-1}\ \ & 又因为 r_{k} = b - Ax_{k}\ = \ & p_{k-1}^\mathrm{T}(b - r_k) + \xi p_{k-1}^\mathrm{T}Ar_{k} + \eta p_{k-1}^\mathrm{T}Ap_{k-1} - b^\mathrm{T}p_{k-1}\ = \ & -p_{k-1}^\mathrm{T}r_{k} + \xi p_{k-1}^\mathrm{T}Ar_{k} + \eta p_{k-1}^\mathrm{T}Ap_{k-1}\ \ & 又因为 p_{k-1} 与 r_{k} 共轭\ = \ & \xi p_{k-1}^\mathrm{T}Ar_{k} + \eta p_{k-1}^\mathrm{T}Ap_{k-1}\ \end{aligned} $$

令 $\frac{\partial\psi}{\partial\xi} = \frac{\partial\psi}{\partial\eta} = 0$，即：

$$ \begin{cases} \tilde\xi r_{k}^\mathrm{T}Ar_{k} + \tilde\eta r_{k}^\mathrm{T}Ap_{k-1} - r_{k}^\mathrm{T}r_k = 0\ \tilde\xi p_{k-1}^\mathrm{T}Ar_{k} + \tilde\eta p_{k-1}^\mathrm{T}Ap_{k-1} = 0 \end{cases} \tag{2} $$

可解得极小值点 $\tilde{x} = x_{k} + \tilde \xi r_{k} + \tilde \eta p_{k-1}$，故

$$ \begin{aligned} \tilde{x} - x_{k} = \tilde \xi r_{k} + \tilde \eta p_{k-1} \end{aligned} $$

由于 $r_k \neq 0$ （因为 $r_k$ 是残差，若为 $0$ 则应该停止迭代了），所以 $\tilde{\xi} \neq 0$（否则更新值只与前一次的下山方向有关，这样无法收敛）可作为分母，因此可取：

$$ \begin{aligned} p_k = \frac{1}{\tilde\xi}(\tilde{x} - x_k) = r_k + \frac{\tilde \eta}{\tilde \xi}p_{k-1} \end{aligned} $$

作为新的下山方向。令 $\beta_{k-1} = \frac{\tilde{\eta}}{\tilde{\xi}}$，则可由式 2 得：

$$ \begin{aligned} \beta_{k-1} = \ & \frac{\tilde{\eta}}{\tilde{\xi}} \ \ & 因为 \frac{\partial\psi}{\partial\eta} = 0 得 \xi p_{k-1}^{\mathrm{T}}Ar_{k} = - \eta p_{k-1}^{\mathrm{T}}Ap_{k-1} \ \ & 左右同除以 -\xi 和 p_{k-1}^{\mathrm{T}}Ap_{k-1} \ = \ & -\frac{p_{k-1}^{\mathrm{T}}Ar_{k}}{p_{k-1}^{\mathrm{T}}Ap_{k-1}} \ \ & 因为p、r是向量，A对称，结果是标量，所以可以交换顺序 \ = \ & -\frac{r_k^{\mathrm{T}}Ap_{k-1}}{p_{k-1}^{\mathrm{T}}Ap_{k-1}} \tag{3} \end{aligned} $$

至此就找到了新的下山方向，向量 $x$ 将每次沿着下山方向，以 $\alpha$ 为步长下降，即 $x_{k+1} = x_{k} + \alpha_{k}p_{k}$，将其带入 $f(x)$ 可得：

$$ \begin{aligned} \phi(\alpha_k) = \ & f(x_k + \alpha_kp_k)\ = \ & \frac{1}{2}(x_k + \alpha_kp_k)^{\mathrm{T}}A(x_k + \alpha_kp_k) - b^{\mathrm{T}}(x_k + \alpha_kp_k) + c\ = \ & f(x_k) + \frac{1}{2}\alpha_{k}^2p_k^{\mathrm{T}}Ap_k + \alpha_{k}x_{k}^{\mathrm{T}}Ap_{k} - \alpha_{k}b^{\mathrm{T}}p_{k}\ \ & 又因为 r_{k}^{\mathrm{T}} = (b - Ax_{k})^{\mathrm{T}} = x_{k}^{\mathrm{T}}A - b^{\mathrm{T}}\ = \ & f(x_k) + \frac{1}{2}\alpha_{k}^2p_k^{\mathrm{T}}Ap_k - \alpha_{k}r_k^{\mathrm{T}}p_k \end{aligned} $$

令 ${\phi}’(\alpha_{k}) = 0$ 可解得：

$$ \begin{aligned} \alpha_k = \ & \frac{r_k^{\mathrm{T}}p_k}{p_k^{\mathrm{T}}Ap_k} \ = \ & \frac{r_k^{\mathrm{T}}(r_k + \beta_{k-1}p_{k-1})}{p_k^{\mathrm{T}}Ap_k} \ \ & 因为 p_{k-1} 与 r_{k} 相互共轭 \ = \ & \frac{r_k^{\mathrm{T}}r_k}{p_k^{\mathrm{T}}Ap_k} \end{aligned} \tag{4} $$

根据 $x_{k+1} = x_{k} + \alpha_{k}p_{k}$ 重新整理 $r_{k+1}$ 可得：

$$ \begin{aligned} r_{k+1} = \ & b - Ax_{k+1}\ = \ & b - A(x_k + \alpha_k p_k)\ = \ & r_k - \alpha_k A p_k\ \end{aligned} $$

移项

$$ \begin{aligned} Ap_k = \frac{1}{\alpha}(r_k - r_{k+1}) \tag{5} \end{aligned} $$

将式 7 带入式 3 得：

$$ \begin{aligned} \beta_{k} = , & -\frac{r_{k+1}^{\mathrm{T}}Ap_{k}}{p_{k}^{\mathrm{T}}Ap_{k}} \ = , & -\frac{\frac{1}{\alpha}r_{k+1}^{\mathrm{T}}(r_k - r_{k+1})}{\frac{1}{\alpha}p_{k}^{\mathrm{T}}(r_k - r_{k+1})} \ , & 因为 p_{k-1} 与 r_{k} 相互共轭，r_{k+1} 与 r_{k} 相互共轭 \ = , & \frac{r_{k+1}^{\mathrm{T}}r_{k+1}}{p_{k}^{\mathrm{T}}r_k} \ = , & \frac{r_{k+1}^{\mathrm{T}}r_{k+1}}{r_{k}^{\mathrm{T}}p_k} \ = , & \frac{r_{k+1}^{\mathrm{T}}r_{k+1}}{r_{k}^{\mathrm{T}}(r_k + \beta_{k-1}p_{k-1})} \ = , & \frac{r_{k+1}^{\mathrm{T}}r_{k+1}}{r_{k}^{\mathrm{T}}r_k} \end{aligned} $$

由此可得到共轭梯度法的一般过程：首先向量 $x$ 任意取一初值 $x_{0}$，由 $r_{0} = b - Ax_{0}$ 得到残差 $r_{0}$，第一步的下山方向 $p_{0} = r_{0}$，根据式 6 更新向量 $x$：

$$ \begin{aligned} x_1 = , & x_0 + \alpha_0p_0\ = , & x_0 + \frac{r_0^{\mathrm{T}}r_0}{p_0^{\mathrm{T}}Ap_0}p_0 \end{aligned} $$

之后的第 $k+1$ 步 $（k \geq 1）$，其上一步的残差 $r_{k}$、新下山方向 $p_{k}$、新的向量 $x_{k+1}$ 分别为：

$$ \begin{aligned} r_{k} = , & r_{k-1} - \frac{r_{k-1}^{\mathrm{T}}r_{k-1}}{p_{k-1}^{\mathrm{T}}Ap_{k-1}}Ap_{k-1}\ p_{k} = , & r_{k} - \frac{r_{k}^{\mathrm{T}}r_{k}}{r_{k-1}^{\mathrm{T}}r_{k-1}}p_{k-1}\ x_{k+1} = , & x_k + \frac{r_k^{\mathrm{T}}r_k}{p_k^{\mathrm{T}}Ap_k}p_{k} \end{aligned} $$

其他注意点

$p_{k-1}$ 与 $r_{k}$ 相互共轭 $r_{k+1}$ 与 $r_{k}$ 相互共轭在共轭梯度法中，相邻的两次搜索得到的梯度相互正交共轭梯度法依赖于线搜索，非精确线搜索不能保证共轭性和下降性子空间扩展定理：

$$ r_{k}^\mathrm{T}p_{j} = 0（j=0,1,…k-1） $$

预共轭梯度法

对矩阵 $A$ 进行预处理，减小 $A$ 的条件数，即 PCG 算法包括 Cholesky 分解预处理、Jacobi 预处理等等首先把矩阵 A 分为严格下三角，对角，以及严格上三角矩阵

$$ A = L + D + L^\mathrm{T} $$

那么 jacobi 预处理矩阵就是 $A$ 的对角矩阵

$$ M = D $$

另外两种比较简单的预处理矩阵分别是 GaussSeidel 以及 SOR

$$ M = L + D \ M = \frac{1}{\omega}(D+\omega L) $$

cholesky 分解

$$ A = LL^\mathrm{T}，其中 A 是下三角矩阵 $$