共轭梯度法全面解析：原理推导、代码实现与预条件实践

共轭梯度法Conjugate Gradient Method简称CG是数值优化里一个承上启下的经典算法。上一篇文章我们聊过最速下降法的锯齿问题也提过牛顿法的二阶信息优势但最速下降法收敛太慢牛顿法要算Hessian矩阵在大规模问题上根本存不下。如果你正在处理稀疏线性方程组、大规模无约束优化或者做神经网络训练中涉及二次子问题的加速共轭梯度法恰好站在两条路中间它只需要计算梯度或矩阵向量积不需要显式存储Hessian却能在线性二次问题上达到n步精确收敛。这篇文章是我这系列第5篇我会从最基础的几何直觉出发把共轭梯度法的数学原理、算法推导、代码实现、预条件处理到非线性扩展一次讲透希望给有数值计算基础但没深究过CG的同学一个完整抓手。1. 共轭梯度法到底解决了什么问题1.1 最速下降法的尴尬最速下降法的思路很直白每次沿着当前梯度反方向走也就是选择函数下降最快的方向。单看一步没问题但连续做几步就露馅了。以二维二次函数 f(x) 1/2 x^T A x - b^T x 为例A为对称正定矩阵。最速下降法走出来的路径是典型的锯齿状由于每次搜索方向都与切线方向垂直下一步的起点梯度又与上一步方向正交两步之间夹角反而可能很大导致路径反复横跳。尤其是当A的条件数很大时椭圆的等值线被拉得很长锯齿现象更加严重。我早期手动实现最速下降法时也验证过一个条件数为1000的二维问题理论上要走几百步才能收敛到足够精度。每一步都看似在下降但两个相邻方向几乎抵消了前面的成果这种“原地打转”的效率谁优化谁知道。1.2 牛顿法的另一面牛顿法之所以强是因为它利用了二阶信息每次迭代更新 x_{k1} x_k - H^{-1} g_k其中H是Hessian矩阵g_k是梯度。在二次问题上牛顿法一步就能到达最优点因为二次函数的Hessian是常数矩阵泰勒展开只有二阶项。但问题是二阶信息是有代价的。首先Hessian矩阵需要显式构造对n维问题要存 n×n 个元素当n达到百万量级时这个存储就是天文数字。其次每一步都需要求解一个线性方程组 H d -g这本身就要做大规模矩阵分解或者迭代求解。所以牛顿法在小规模问题上很香但一放到现代工程场景里比如求解有限元方程组、训练大规模机器学习模型内存和计算量往往双双爆炸。1.3 共轭梯度法:走第三条路共轭梯度法的核心思想是我能不能不存储Hessian矩阵只利用矩阵向量积MatVec来隐式获取二阶信息对于二次函数 f(x) 1/2 x^T A x - b^T x梯度就是 r(x) A x - b而和A做一次乘法 A p 就隐式包含了Hessian信息因为A就是Hessian。共轭梯度法每一步只需要做一次 A p 的矩阵向量积不需要存储完整的A也不需要做矩阵分解却能在线性二次问题上保证最多n步精确收敛。这个折中非常聪明它比最速下降法多了对“历史搜索方向”的利用比牛顿法省掉了Hessian的存储与求逆。正因如此CG成为大规模稀疏正定线性系统最经典的迭代解法也是非线性优化领域L-BFGS出现之前最重要的无约束优化算法之一。2. 共轭方向的数学原理2.1 从正交到A-共轭在标准欧几里得空间中两个向量正交是指它们的内积为零u^T v 0。共轭梯度法把内积换成了“以A为度量”的内积如果 u^T A v 0就称u和v关于A共轭也叫A-正交。为什么要换度量因为目标函数的等值线是椭圆族椭圆的长短轴方向与A的特征向量方向相关而普通正交系并不匹配这个几何结构。若我们在一组A-共轭方向上做搜索每个搜索方向都能独立消除该方向上的误差分量不会像最速下降法那样反复横跳。打个生活化的比方想象你在一间长方形房间里要从一个墙角走到对角墙角。普通最短路径是直接斜穿但如果你的导航只允许你“水平走一步、竖着走一步”地前进就会走出锯齿路线。现在换一套坐标系让坐标轴本身斜着对齐房间的对角线那么你沿着新坐标轴各走一步两步就到对角了。A-共轭本质上就是把坐标系旋转到适应目标函数椭圆形状的方向。2.2 n步收敛性的直觉证明在n维空间中任意一组非零且两两A-共轭的向量 p_0, p_1, ..., p_{n-1}必然是线性无关的因此它们构成R^n的一组基。设二次函数的最小点是x*当前点为x_0。误差 e_0 x_0 - x* 可以在这组基下展开e_0 Σ_{i0}^{n-1} c_i p_i当我们在 p_j 方向做精确线搜索时极小化 f(x_j α p_j) 等价于在该方向上把误差分量调整为零。由于所有搜索方向之间A-共轭调整 p_j 方向的分量不会破坏前面已经调整好的 p_0...p_{j-1} 方向的分量。于是每步消掉一个方向上的误差n步之后误差彻底清零。这比最速下降法强在哪里最速下降法的搜索方向是依据当前梯度临时决定的相邻两步方向正交但不一定与更早的方向兼容误差分量会被反复引入又消掉。共轭梯度法则保证每一步都在之前所有方向张成的空间外寻找新分量不走回头路。2.3 子空间视角渐进最优共轭梯度法还有一个极其重要的性质叫“逐步最优性”经过k步之后当前迭代点 x_k 在仿射空间 x_0 span{p_0,...,p_{k-1}} 中已经使二次函数取最小值。换句话说迭代点是在Krylov子空间 K_k(A, r_0) span{r_0, Ar_0, ..., A^{k-1} r_0} 中对目标函数的最优近似解。这个性质可以理解为CG不是盲目地在n个方向上一扫而过它每走一步都在“用前k步构造的低维子空间里找当前最优点”。所以在实际计算中即使我们等不到n步只做几十步得到的解也往往比最速下降法几百步得到的解好得多。3. 完整推导CG算法的五条迭代公式3.1 精确线搜索的步长在给定搜索方向 p_k 后我们需要确定步长 α_k使得 f(x_k α p_k) 最小。令φ(α) 1/2 (x_k α p_k)^T A (x_k α p_k) - b^T (x_k α p_k)求导并令其为零φ(α) p_k^T A (x_k α p_k) - b^T p_k p_k^T (A x_k - b) α p_k^T A p_k 0定义残差 r_k A x_k - b g_k梯度则α_k - (r_k^T p_k) / (p_k^T A p_k)注意由于梯度是负方向实际公式写作 α_k (r_k^T p_k) / (p_k^T A p_k)。这里分母 p_k^T A p_k 就是“沿p_k方向的二阶方向导数”对正定矩阵A恒为正。3.2 残差的递推更新一个关键技巧在于残差并不需要每次用 r_k A x_k - b 重新计算。因为x_{k1} x_k α_k p_k直接代入残差定义r_{k1} A x_{k1} - b A(x_k α_k p_k) - b (A x_k - b) α_k A p_k r_k α_k A p_k从梯度角度看r_{k1} r_k α_k A p_k。注意α_k是负数时这就是残差沿方向减去。实际代码中我会用 r_{k1} r_k - α_k A p_k 这种形式配合负号的语义差异小心处理。更标准的写法是定义 r_k b - A x_k那么残差更新就是 r_{k1} r_k - α_k A p_k。两种符号体系在算法实现时容易混淆我下面给的代码统一用 r b - A x 的约定这也是数值线性代数最常见的约定。这个递推的意义很大一次迭代只需要一次矩阵向量积 A p_k而不需要额外计算 A x_{k1}。3.3 搜索方向的更新与β有了步长α之后我们需要构造下一个搜索方向 p_{k1}。共轭梯度法的搜索方向由当前残差加上前一个方向的修正构成p_{k1} r_{k1} β_k p_k系数β_k的选择要保证 p_{k1} 与 p_k 关于A共轭即满足p_{k1}^T A p_k 0代入得(r_{k1} β_k p_k)^T A p_k 0因此β_k - (r_{k1}^T A p_k) / (p_k^T A p_k)再利用 A p_k (r_{k1} - r_k) / (-α_k) 的关系以及残差之间的正交性可以化简为β_k (r_{k1}^T r_{k1}) / (r_k^T r_k)这就是共轭梯度法最经典的FRFletcher-Reeves公式。后面我会在非线性扩展部分讨论PRP、HS等其他更新公式它们在二次正定问题中理论等价但数值表现和非线性场景下有差异。3.4 标准算法流程综合上面推导共轭梯度法解线性方程组 A x b 的完整流程如下输入对称正定矩阵A右端项b初始猜测x0容差tol 1. r0 b - A x0 2. p0 r0 3. 对 k 0, 1, 2, ... 重复 a. 计算 A p_k b. α_k (r_k^T r_k) / (p_k^T A p_k) c. x_{k1} x_k α_k p_k d. r_{k1} r_k - α_k A p_k e. 如果 ||r_{k1}|| / ||r0|| tol退出 f. β_k (r_{k1}^T r_{k1}) / (r_k^T r_k) g. p_{k1} r_{k1} β_k p_k 4. 输出 x值得注意这里的α表达式已经用 r_k^T r_k 替换掉了 r_k^T p_k因为p_k由残差和前向修正而成可以证明 r_k^T p_k r_k^T r_k。如果直接用一般线搜索公式 α r^T p / p^T A p 也是正确的但在CG语境下使用简化形式误差传播更小。4. 动手实现从零写一个线性CG求解器4.1 一个足够实用的Python实现下面这段代码可以直接复制使用。它的核心是“只调用矩阵向量积”所以即使A是隐式定义的比如算子也可以正常工作import numpy as np def cg(A, b, x0None, tol1e-8, max_iterNone): 共轭梯度法求解 Ax b 要求 A 对称正定或者至少能计算 A p 的向量积 b np.asarray(b, dtypenp.float64) n b.shape[0] x np.zeros(n, dtypenp.float64) if x0 is None else np.asarray(x0, dtypenp.float64).copy() r b - A x p r.copy() rs_norm_sq r r rhs_norm np.linalg.norm(b) if rhs_norm 0.0: rhs_norm 1.0 if max_iter is None: max_iter n history [] for k in range(max_iter): Ap A p pAp p Ap # 处理数值奇异 if pAp 0.0: break alpha rs_norm_sq / pAp x x alpha * p r r - alpha * Ap rs_norm_sq_new r r history.append(np.sqrt(rs_norm_sq_new) / (np.linalg.norm(b) 1e-12)) if np.sqrt(rs_norm_sq_new) tol * rhs_norm: break beta rs_norm_sq_new / rs_norm_sq p r beta * p rs_norm_sq rs_norm_sq_new return x, history在这个实现里我用了相对残差作为终止条件也就是 ||r_k|| / ||b||这样可以避免因为b的尺度不同而出现过早或过晚停止。这个细节在实际中很重要很多人直接判断 ||r_k|| tol遇到b本身模长极大或极小的场景就会出问题。4.2 数值行为验证用一个小规模正定矩阵测试一下。比如生成一个100阶的对称正定矩阵构造已知解然后跑CGnp.random.seed(42) n 100 A np.random.randn(n, n) A A A.T n * np.eye(n) # 保证正定 x_true np.ones(n) b A x_true x_cg, hist cg(A, b, tol1e-10) print(误差:, np.linalg.norm(x_cg - x_true)) print(迭代步数:, len(hist))实际跑下来在精确算术的模拟下CG大概迭代几次到几十次就能达到1e-10级别的相对残差具体步数取决于A的条件数。如果拿同一个问题用最速下降法跑迭代几百步可能还没到同样精度。直观经验是当条件数在1e3量级时CG大约需要几十步最速下降法可能需要几千步。这也是前面理论中“k步最优”在实践中的体现。4.3 二维可视化对比为了直观感受共轭梯度法和最速下降法的区别可以构造一个二维二次函数分别从同一起点出发把迭代路径画出来。你会在图上看到最速下降法的路径呈锯齿状相邻两步互相垂直走的步数很多共轭梯度法走两步就到最优点而且第二步的方向和第一步关于A共轭路径干净利落。这种可视化建议每个学优化的人都亲自做一次。因为一旦你亲眼见过锯齿路径和两步直达的对比你对共轭梯度法的信任感会牢牢记在脑子里比看十页公式都管用。5. 预条件共轭梯度法让收敛速度飞起来5.1 为什么需要预条件CG的收敛速度理论上依赖于矩阵A的条件数 κ(A)。具体来说误差按 (√κ - 1)/(√κ 1) 的速率衰减。如果κ很大比如1e6这个衰减率约等于 0.998意味着每次迭代误差只缩小千分之二需要几千步才能收敛。实际问题中有限元网格加密、偏微分方程离散化、正则化项跨尺度等都会导致矩阵条件数暴涨。这时候单纯靠裸CG会非常尴尬看似每一步都在工作几百步下来误差纹丝不动。预条件的基本思想是找到矩阵M ≈ A^{-1}或者M真正接近A而逆容易做把原方程等价转换使得新矩阵的条件数大幅下降。极端情况是M A^{-1}直接一步解出但M^{-1}的计算代价要和CG迭代收益做权衡实际中选中等强度的预条件子最划算。5.2 左预条件与对称性破坏问题最朴素的想法是求解 M^{-1} A x M^{-1} b。但问题来了即使M和A都对称正定M^{-1} A未必对称CG的正交性推导全部依赖对称性直接套CG是错的。解法之一是用“A-内积视角”推导PCG或者用拆分形式。更常见的做法是把预条件矩阵做分解 M L L^T然后用变量代换使新方程恢复对称正定形式。但实际代码里没有必要显式做变量代换PCG算法会把预条件操作内嵌到每一步。5.3 PCG算法核心变化标准的预条件共轭梯度算法比CG多一个步骤每次计算残差 r_k 后先解一个辅助系统 M z_k r_k得到“预条件残差”z_k然后用z_k代替r_k去更新方向和计算系数。r0 b - A x0 z0 M^{-1} r0 p0 z0 重复 Ap_k A p_k α_k (r_k^T z_k) / (p_k^T A p_k) x_{k1} x_k α_k p_k r_{k1} r_k - α_k A p_k 如果收敛则退出 z_{k1} M^{-1} r_{k1} β_k (r_{k1}^T z_{k1}) / (r_k^T z_k) p_{k1} z_{k1} β_k p_k注意α和β的分子都用 r^T z 而不是 r^T r。当M I时PCG退回普通CG。5.4 常见的预条件子选型对角预条件JacobiM diag(A)构造代价为零适合对角线占优的矩阵。很多时候能显著改善条件数而且没有额外存储。不完全CholeskyIC对A做不完全分解保留稀疏模式M L L^T。效果通常比对角好但实现复杂并行化略费劲。SSOR对称逐次超松弛M (D/ω L)(D/ω)^{-1}(D/ω L)^T其中D是A的对角L是严格下三角。效果介于两者之间且调参ω取1.5左右常见。经验上如果矩阵来自PDE离散化IC或IC(0)基本是标配如果只是普通稠密对称正定矩阵传统做法可能就是直接Cholesky分解了谈不上迭代法。6. 非线性优化里的共轭梯度法6.1 从线性到非线性如果目标函数不是二次函数比如 f(x) Σ log(1 exp(-y_i x^T z_i))那么这个函数不能写成 1/2 x^T A x - b^T x 的形式前面推导的精确线搜索和残差递推不再直接适用。但非线性共轭梯度法Nonlinear CG保留了算法的骨架每次迭代先选择一个搜索方向 p_k然后做线搜索确定步长 α_k再用新的梯度 g_k 构造下一个方向。最直观的方式是把上述CG公式里的“残差r”替换成“负梯度 -g”p_k -g_k β_k p_{k-1}这个形式非常简洁。剩下的问题只有一个β_k怎么选。6.2 常见的β更新公式下面列出实践中几类最常用的β公式用 p 表示搜索方向g 表示梯度FR公式β (g_{k1}^T g_{k1}) / (g_k^T g_k)。这是线性CG公式的直接推广理论性质好但有时收敛偏慢。PRP公式β (g_{k1}^T (g_{k1} - g_k)) / (g_k^T g_k)。实践中往往比FR更稳尤其是目标函数非二次时能自动应对梯度的变化。HS公式β (g_{k1}^T (g_{k1} - g_k)) / (p_k^T (g_{k1} - g_k))。这个公式与“拟牛顿条件”关系最紧密数值上也很常用。从实际经验来看PRP型公式在大规模非线性问题里表现通常优于FR因为当迭代点前进缓慢时梯度变化不大PRP的分子会显著变小自动减弱历史方向的权重避免搜索方向走偏。FR在理论收敛性分析上更优美但精确算术假设在非线性数值噪声环境下不那么可靠。6.3 非线性CG的两个关键约束第一线搜索必须满足充分下降条件。CG方法不像牛顿法有天然下降方向保证β选得不好时 p_k^T g_k 可能大于0即搜索方向不再是下降方向。这时候必须强制检查 p_k^T g_k 0否则把p_k直接重置为 -g_k。第二Wolfe条件几乎是标配。也就是说要找的步长既要让函数值明显下降Armijo条件又不能让梯度变化太剧烈曲率条件。只做回溯线搜索也能跑但收敛速度可能大打折扣。很多开源优化器的非线性CG实现里还会加一个“重启策略”比如每n步重置一次把方向恢复为负梯度或者当相邻梯度相关性下降时重启。这本质上是放弃陈旧的子空间信息避免累积数值误差。7. 实际使用中的调参技巧和踩坑7.1 终止条件怎么设才合理就线性CG而言我强烈建议用相对残差作为停止指标。标准是 ||r_k|| / ||b|| tol。如果问题本身右端项b为零向量那基准换为1或直接看 ||r_k||。对于非线性CG终止条件一般看梯度范数的相对下降比如 ||g_k|| / max(1, ||g_0||) 1e-6。还要注意不要只看函数值变化因为很多工程问题的函数值尺度很大下降极小但梯度还远没有归零容易提前停止。7.2 浮点误差和正交性丢失理论上CG的残差天然满足 r_i^T r_j 0 (i≠j)但浮点运算下这个正交性会逐步丢失。步数越多错误累积越明显。因此对于病态问题可能需要每几轮做一次“重新正交化”或者用Modified Gram-Schmidt处理搜索方向。另一个更实用的建议是如果迭代轮数超过了n而且还在慢慢收敛别硬加迭代先考虑是不是预条件不足或者直接换用带重启的CG变体。7.3 不是所有矩阵都能用CGCG只适用于对称正定矩阵。如果A非对称但对称部分正定可以用CG处理对称化后的方程但不保证最优。常见替代解是GMRES、BiCGStab、MINRES等Krylov子空间方法。我见过一个真实场景某开发者把CG套在一个非对称矩阵上跑了半天不收敛最后检查发现A的对称部分还好但反对称项很大换BiCGStab后几十步就收敛了。所以第一步永远先检查A的性质。在优化器选型上如果问题规模小n 1000直接上牛顿法Cholesky别犹豫如果问题规模大但梯度易算且内存有限非线性CG和L-BFGS是两大主力L-BFGS往往更快但需要多点存储CG内存更省。7.4 使用预条件时要小心的事预条件主要是为了降条件数但不是越强越好。一个过强的预条件相当于每步都在解一个近似原方程的子问题代价高昂。我的习惯是先跑几十步裸CG看残差曲线如果残差下降得很快就不需要预条件如果残差曲线出现“平缓长尾”再加预条件优先从对角预条件开始。8. 写在最后的一些个人体会把这五年用过CG的经验浓缩一下共轭梯度法真正的价值不在n步收敛这个理论亮点而在于它是一种“给点历史信息就起飞”的算法。最速下降法只在乎当前梯度牛顿法在乎全部二阶信息CG恰好站在中间——只记住过去几个搜索方向就能在Krylov子空间中渐进最优。这种“用最少记忆换最大信息利用”的思路几乎贯穿了整个数值优化领域。在实际操作中我建议每个做优化算法的人都把CG亲手写一遍甚至写在纸上推一遍公式。因为只有当你自己推导出 α r^T r / p^T A p 和 β r^T_ r_ / r^T r 这两个式子你才会真正明白为什么CG要如此设计。另外实现时务必把残差递推和方向更新的代码分开调试先测一步是否精确下降再测两步是否满足共轭关系最后再跑完整收敛测试。这个小习惯帮我排掉了好几个符号错误。共轭梯度法到现在仍然是求解大规模稀疏正定系统最实用的迭代法之一值得你花一晚上把它彻底搞懂。