想象你被空投到一个浓雾弥漫的碗形山谷，要寻找最低点。你被蒙着眼睛，只能用登山杖试探脚下坡度，一小步一小步摸索。最快的方法是“哪里最陡往哪里走”——这就是最速下降法，但很快你会发现自己在狭长的山谷里来回打转，走出令人崩溃的“之”字形。更奢侈的做法是动用直升机航拍整片地形，一步直达谷底——这就是牛顿法，可惜对百万维问题，航拍成本高到无法承受。今天的主人公——共轭梯度法（CG），是一位“会记账的徒步者”：它每完成一个方向，就把这个方向从脑海中划掉，下次只朝与所有老方向“互不干扰”的新方向走。因为从不回头，它最多用 n 步就能到达谷底，而每一步只需轻轻掂量一下背包——做一次矩阵向量乘法。这就是解大规模线性方程组的利器。

岗位关联

掌握共轭梯度法，等于握住了科学计算、机器学习与大规模优化的通用钥匙——它是从经典优化走向现代工程求解器的桥梁，也是大模型时代每个 AI 工程师值得拥有的底层数值能力。

目录
  1. 可视化演示
  2. 问题背景：矩阵与方程
  3. 最速下降：之字形困境
  4. 共轭方向：互不干扰的路径
  5. CG算法：动态构造方向
  6. 收敛分析：山谷形状决定速度
  7. 预处理技术：把山谷改造成圆碗
  8. 非线性扩展：翻越连绵群山
  9. 工程实践：从公式到求解器
  10. AI 应用场景
  11. Python 演示
  12. 常见误区
  13. 面试题

可视化演示

展示共轭梯度法在二次优化问题上的共轭方向搜索路径，对比最速下降的锯齿收敛，交互切换模式观察预处理加速效果。

本节配有交互动效（进入视口后自动加载，约 3KB）。

共轭梯度法：共轭方向搜索 靛蓝为CG路径（A-共轭方向，有限步收敛） 橙色为起点（初始解x0） 绿色为最优解（目标x*） p_i^T A p_j = 0 (i!=j), CG converges in n steps 共轭梯度法利用A-共轭方向在Krylov子空间中搜索，对n维系统最多n步收敛

问题背景：矩阵与方程

继续我们的雾谷之旅。在真实工程里，Ax=b 中的 A 往往巨大无比——流体力学仿真有上百万个网格点，金融风控的协方差矩阵横跨几万只股票，但它们大都“稀疏”：每个方程只和相邻少数变量打交道，就像山谷中每寸土地只和旁边的土地直接相连。更重要的是，A 必须是对称正定矩阵（SPD）。它保证了碗形山谷是规则的、四周收紧的，绝对没有马鞍形的假谷底或山脊线，这让“顺坡向下走”的策略永远不会被误导。

把找谷底变成数学语言，我们会惊讶地发现：解方程组和找最低点是同一件事。考虑二次函数 f(x) = ½xᵀAx − bᵀx，对它求梯度，得到 ∇f(x) = Ax − b。令梯度等于零，方程正是 Ax = b。所以谷底 x* 就是方程的解。对称正定保证 f 是严格凸函数，像一只光滑的碗——任何“脚下发滑”的局部最低点，其实就是全局最低点，不存在迷惑人的小坑。

于是，我们的战役从“硬解方程组”巧妙地变成了“在碗里找底”。接下来的所有方法，都是不同风格的徒步策略：最速下降是冲动型选手，共轭梯度则是聪明而节制的探险家。

f(x)=12xTAx−bTx,∇f(x)=Ax−b\displaystyle f(\mathbf{x}) = \frac12 \mathbf{x}^{\mathsf{T}} A \mathbf{x} - \mathbf{b}^{\mathsf{T}} \mathbf{x}, \qquad \nabla f(\mathbf{x}) = A\mathbf{x} - \mathbf{b}
二次函数 f 的地形图完全由 A 决定；让梯度为零，就得到原始线性方程组。
  1. 二次型½xᵀAx 表示“自身与自身经过 A 耦合后的能量”，它决定了碗壁如何弯折。
    12xTAx\displaystyle \frac12 \mathbf{x}^{\mathsf{T}} A \mathbf{x}
  2. 线性项−bᵀx 像一个外来的推力，把整个碗的底部整体搬离原点。
    −bTx\displaystyle -\mathbf{b}^{\mathsf{T}} \mathbf{x}
  3. 梯度为零∇f = 0 意味着脚下完全平坦、不再有下降方向，正是谷底的位置。
    ∇f(x)=Ax−b=0\displaystyle \nabla f(\mathbf{x}) = A\mathbf{x} - \mathbf{b} = \mathbf{0}

最速下降：之字形困境

先看最直接、也最符合直觉的走法：每次用登山杖探出当前最陡的下坡方向，然后大步走下去。数学上，这个方向就是负梯度 −∇f(x_k)，即残差 r_k = b − Ax_k。选定方向后，还要决定走多远——步长 α 要让 f 沿着这条直线降得最低。对二次函数，这个一维优化有闭式解：αₖ = (rₖᵀrₖ)/(rₖᵀ A rₖ)。

按此规则前进，新位置的残差 rₖ₊₁ 与刚走过的方向 rₖ 天然正交。听起来很合理，但问题恰恰出在这里：一旦山谷是瘦长形的，相邻两次的搜索方向几乎互相垂直，路径就会像乒乓球一样在谷壁之间来回弹跳。每一步都在“横穿山谷”，真正的推进却很慢。这就是最速下降法著名的锯齿现象。

最速下降教会我们两件事：一是沿某方向走到最优后，继续沿原方向毫无意义；二是把所有希望押在“当下最陡”上，会遗忘已经走过的路。那位聪明的徒步者于是有了新想法：既然重复走老方向是浪费，那就把所有走过的方向记在账本上，只挑真正没走过的方向走。

αk=rkTrkrkTArk,xk+1=xk+αkrk\displaystyle \alpha_k = \frac{\mathbf{r}_k^{\mathsf{T}} \mathbf{r}_k}{\mathbf{r}_k^{\mathsf{T}} A \mathbf{r}_k}, \qquad \mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{r}_k
步长公式保证沿当前方向一次性走到最低点，但下一方向又可能与之前的方向“打架”，导致反复横跳。
  1. 残差即负梯度rₖ = b − Axₖ 告诉我们脚下最陡的下坡方向。
    rk=b−Axk\displaystyle \mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k
  2. 步长搜索αₖ 使 f(xₖ + αrₖ) 沿该直线最小，相当于在此方向上走到谷底边缘。
    αk=rkTrkrkTArk\displaystyle \alpha_k = \frac{\mathbf{r}_k^{\mathsf{T}} \mathbf{r}_k}{\mathbf{r}_k^{\mathsf{T}} A \mathbf{r}_k}
  3. 位置更新新地点 xₖ₊₁ 是旧地点沿最陡方向迈出的最优一步。
    xk+1=xk+αkrk\displaystyle \mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{r}_k

共轭方向：互不干扰的路径

徒步者的账本，数学上叫“共轭方向”。两条方向 dᵢ 和 dⱼ 被称为 A-共轭（A-正交），如果它们满足 dᵢᵀ A dⱼ = 0。这个条件不像普通正交那么直观，但它有一个绝妙的生活对应：想象你在一座山谷中推行独轮车，推完第 i 条坡道后，再推第 j 条坡道时，车轮不会滑回刚才的老路——因为两条坡道在“地形几何”上完全独立。

为什么要坚持这种独立？因为二次函数的最小化可以分解成 n 个一维最小化问题。如果有 n 条两两 A-共轭的方向，那么沿着每一条方向做一次精确一维搜索，结果就是全局最优：沿 d₁ 走到谷底在 d₁ 方向上的投影，再沿 d₂ 走时不会干扰 d₁ 上已经取得的最优。逐次做完 n 个方向，恰好走完整个 n 维空间。这就像把一个大拼图拆成 n 块互不重叠的拼板，每块单独摆好，整体自然完成。

共轭方向法听起来很美，但小字印刷写着：手动构造一整套共轭方向需要事先知道 A 的全部特征信息，这就像提前画完整张山谷地图——成本又回到了牛顿法的水平。好在共轭梯度法用递归方式绕开了这个难题。

diTAdj=0(i≠j)\displaystyle \mathbf{d}_i^{\mathsf{T}} A \mathbf{d}_j = 0 \quad (i \neq j)
A-共轭不是普通垂直，而是在 A 定义的内积下正交；它保证沿新方向优化时，不破坏旧方向的成果。
  1. A-内积把矩阵 A 当作度量，重新定义两个向量的“夹角”。
    ⟨di,dj⟩A=diTAdj\displaystyle \langle \mathbf{d}_i, \mathbf{d}_j \rangle_A = \mathbf{d}_i^{\mathsf{T}} A \mathbf{d}_j
  2. 正交条件当这个内积为零时，两方向在地形几何下互不纠缠。
    diTAdj=0\displaystyle \mathbf{d}_i^{\mathsf{T}} A \mathbf{d}_j = 0
  3. n步分解若取得 n 条两两 A-共轭的方向，则 x* 可表示为这些方向上的最优步长之和。
    x∗=∑i=0n−1αidi\displaystyle \mathbf{x}^* = \sum_{i=0}^{n-1} \alpha_i \mathbf{d}_i

CG算法：动态构造方向

共轭梯度法的天才之处在于：不需要预先准备 n 条共轭方向，而是边走边生成。第 k+1 步的新方向 dₖ₊₁ 由当前残差 rₖ₊₁ 与上一条方向 dₖ 线性组合而成：dₖ₊₁ = rₖ₊₁ + βₖ dₖ。系数 βₖ 的作用是让新方向与前面所有方向都保持 A-共轭。残差在更新时也不必重新计算 Ax，只需用上一步的矩阵向量乘积递推：rₖ₊₁ = rₖ − αₖ A dₖ。因此，每一步真正的“重活”只有一次矩阵向量乘法 A dₖ。

这正是“会记账的徒步者”的秘密：账本上不需要记录 n 条旧方向，只需要记住刚走过的那一条，因为已经走过的信息全都凝固在残差和方向的关系里。在数学上，CG 在第 k 步是在 Krylov 子空间 Kₖ = span{r₀, Ar₀, …, A^{k−1}r₀} 上做最优搜索，而这个子空间只需反复“乘以 A”就能扩张。于是，存储成本从 O(n²) 降到 O(n)，每次迭代成本也只剩一次稀疏矩阵向量乘。

现在，这位徒步者每次只需记一条路、迈一步新的“棱线”，背上的包越来越轻。理论上，走完 n 步必到谷底；而在更复杂的真实地形中，它通常也能远快于最速下降。但“能走多快”还取决于山谷本身的形状——这就是下一节要谈的收敛分析。

βk=rk+1Trk+1rkTrk,dk+1=rk+1+βkdk,αk=rkTrkdkTAdk\displaystyle \beta_k = \frac{\mathbf{r}_{k+1}^{\mathsf{T}} \mathbf{r}_{k+1}}{\mathbf{r}_k^{\mathsf{T}} \mathbf{r}_k}, \qquad \mathbf{d}_{k+1} = \mathbf{r}_{k+1} + \beta_k \mathbf{d}_k, \qquad \alpha_k = \frac{\mathbf{r}_k^{\mathsf{T}} \mathbf{r}_k}{\mathbf{d}_k^{\mathsf{T}} A \mathbf{d}_k}
CG 通过 β 递推构造共轭方向，通过 α 确定步长；全程只使用矩阵向量乘，不触碰完整矩阵。
  1. 残差递推rₖ₊₁ = rₖ − αₖ A dₖ，避免重新计算 Ax，减少一次大矩阵乘积。
    rk+1=rk−αkAdk\displaystyle \mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k A\mathbf{d}_k
  2. 方向更新ß 系数把残差与旧方向组合，生成与之前所有方向都共轭的新方向。
    dk+1=rk+1+βkdk\displaystyle \mathbf{d}_{k+1} = \mathbf{r}_{k+1} + \beta_k \mathbf{d}_k
  3. 步长确定αₖ 使当前方向上能量误差最小化，保证新残差与旧残差正交。
    αk=rkTrkdkTAdk\displaystyle \alpha_k = \frac{\mathbf{r}_k^{\mathsf{T}} \mathbf{r}_k}{\mathbf{d}_k^{\mathsf{T}} A \mathbf{d}_k}

收敛分析：山谷形状决定速度

同样是会记账的徒步者，在圆形水碗里走一步就能到中心，在一条细长峡谷里却要挣扎很久。决定这个差异的量叫条件数 κ(A)，它等于 A 的最大特征值与最小特征值之比 λmax/λmin，描述的是“山谷最长轴与最短轴的比例”。κ 越接近 1，碗越圆，CG 收敛越快；κ 越大，山谷越瘦长，每步可推进的有效距离就越短。

CG 的收敛速率有一个漂亮的界：在 A-范数意义下，第 k 步的误差至多是初始误差的 2[(√κ−1)/(√κ+1)]ᵏ 倍。这个式子告诉我们两件事：第一，收敛是线性的（在每步缩减率的意义上），且缩减率由 √κ 决定而不是 κ；第二，理论上 n 维问题最多 n 步必达精确解，这是共轭方向法家族的“遗传天赋”。实际计算中由于舍入误差，这个有限步性质会退化，但收敛速率公式仍然是预测行为的最佳指南。

读到这里，聪明的读者已经发现：如果山谷太瘦长，再好的徒步者也会走得气喘吁吁。与其硬着头皮走，不如先改变地形——把瘦长峡谷改造成圆碗。这就引出了预处理技术。

∥xk−x∗∥A≤2(κ−1κ+1)k∥x0−x∗∥A,κ=λmax⁡λmin⁡\displaystyle \|\mathbf{x}_k - \mathbf{x}^*\|_A \le 2 \left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^k \|\mathbf{x}_0 - \mathbf{x}^*\|_A, \qquad \kappa = \frac{\lambda_{\max}}{\lambda_{\min}}
收敛速度由 κ 的平方根主导；条件数越接近 1，每步误差缩小得越狠。
  1. A-范数‖v‖_A = √(vᵀAv) 以地形能量衡量误差，比普通欧氏长度更贴合优化目标。
    ∥v∥A=vTAv\displaystyle \|\mathbf{v}\|_A = \sqrt{\mathbf{v}^{\mathsf{T}} A \mathbf{v}}
  2. 条件数κ 是最大与最小特征值之比，刻画碗的“长宽比”。
    κ(A)=λmax⁡λmin⁡\displaystyle \kappa(A) = \frac{\lambda_{\max}}{\lambda_{\min}}
  3. 缩减因子r = (√κ−1)/(√κ+1) 是每一步误差的上界缩减率；κ 越小，r 越小。
    r=κ−1κ+1\displaystyle r = \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}

预处理技术：把山谷改造成圆碗

预处理的思想非常直白：既然 κ 越大越难走，那就找一个与 A 接近、又很容易求逆的矩阵 M，把原方程 Ax=b 改写成 M⁻¹Ax = M⁻¹b。新矩阵 M⁻¹A 的特征值被“压缩”在 1 附近，等效条件数大幅下降，CG 的每步收敛速度也随之飙升。这相当于在徒步前先雇一支工程队，把狭长峡谷的陡壁削平，让它变成缓和的圆碗——地形改造费花在刀刃上，比多走几百步划算得多。

最常见的两种预处理是雅可比预处理器和对不完全 Cholesky 分解。雅可比方法最简单：M 只取 A 的对角线，相当于只根据每个方向的“局部陡峭度”做缩放，适合对角占优问题。不完全 Cholesky 则更进一步：对 A 做 Cholesky 分解时，只保留稀疏模式内的元素，丢掉分解过程中产生的填充元，形成一个仍然稀疏、但十分接近 A 的因子。这个 M 质量更好，代价是前期构造稍贵。

完成地形改造后，CG 甚至可以用到非线性世界。把碗形山谷换成连绵起伏的群山，把“残差”换成“梯度”，一套同样聪明的记账逻辑仍然适用——非线性共轭梯度法由此诞生。

M−1Ax=M−1b,M≈A,M−1 易于计算\displaystyle M^{-1} A \mathbf{x} = M^{-1} \mathbf{b}, \qquad M \approx A, \quad M^{-1} \text{ 易于计算}
预处理不改变精确解，只改变地形的“等高线”；M 越接近 A，M⁻¹A 越接近单位矩阵，CG 收敛越快。
  1. 等价变换两边同乘 M⁻¹，解不变，但迭代中的曲率被修正。
    M−1Ax=M−1b\displaystyle M^{-1} A \mathbf{x} = M^{-1} \mathbf{b}
  2. 近似要求M≈A 使 M⁻¹A 的特征值聚集在 1 附近，条件数接近 1。
    M−1A≈I\displaystyle M^{-1} A \approx I
  3. 可逆成本M 必须容易分解，否则预处理本身成为新的瓶颈。
    M=LLT,其中 L 保持稀疏\displaystyle M = L L^{\mathsf{T}}, \quad \text{其中 } L \text{ 保持稀疏}

非线性扩展：翻越连绵群山

线性 CG 的公式建立在对二次函数精确的一维搜索上，但真实优化问题往往是高维、非二次的。非线性共轭梯度法的思路，是把残差 r 换成一般目标函数的梯度 g = ∇f(x)，步长不再有闭式解，改用线搜索（如 Wolfe 条件）确定。这样，我们依然可以“记住走过的方向”，但地形本身随时在变化，共轭性质不再保证严格成立，只能近似维持。

β 的计算也有多种版本。最经典的 Fletcher-Reeves（FR）公式直接用相邻两步梯度的范数平方之比：βₖ^FR = (gₖ₊₁ᵀgₖ₊₁)/(gₖᵀgₖ)。另一种流行的 Polak-Ribière-Polyak（PRP）公式引入相邻梯度之差，在目标函数高度非线性时通常更稳健：βₖ^PRP = (gₖ₊₁ᵀ(gₖ₊₁−gₖ))/(gₖᵀgₖ)。这两种公式都只需梯度信息，不需要海森矩阵，因此特别适合大规模无约束优化。

非线性 CG 是线性 CG 的“远房亲戚”：血统相同，但更灵活。在实际工程问题中，它常作为大规模优化算法的核心引擎。终于，我们走到旅程的终点——看看这位徒步者如何在真实世界中大显身手。

βkFR=gk+1Tgk+1gkTgk\displaystyle \beta_k^{\mathrm{FR}} = \frac{\mathbf{g}_{k+1}^{\mathsf{T}} \mathbf{g}_{k+1}}{\mathbf{g}_k^{\mathsf{T}} \mathbf{g}_k}
FR 公式用梯度模长平方之比更新方向；它把线性 CG 的递推逻辑平移到一般光滑优化。
  1. 梯度替代残差在非线性问题中，负梯度仍是最陡下降方向。
    gk=∇f(xk)\displaystyle \mathbf{g}_k = \nabla f(\mathbf{x}_k)
  2. 梯度的平方范数分母与分子分别记录新旧梯度的“力度”。
    gkTgk=∥gk∥22\displaystyle \mathbf{g}_k^{\mathsf{T}} \mathbf{g}_k = \|\mathbf{g}_k\|_2^2
  3. 方向组合β 控制旧方向对新方向的“记忆程度”，避免锯齿。
    dk+1=−gk+1+βkFRdk\displaystyle \mathbf{d}_{k+1} = -\mathbf{g}_{k+1} + \beta_k^{\mathrm{FR}} \mathbf{d}_k

工程实践：从公式到求解器

当我们真正面对百万变量的工程系统时，理论上的“n 步收敛”要打折扣，因为浮点误差会破坏共轭关系，实际迭代次数往往超过 n。这时工程实践的关键是：用合适的终止准则、预处理策略和重启机制。最常用的终止判据是相对残差范数 ‖rₖ‖₂ ≤ ε‖b‖₂，比如 ε=1e−6，相当于徒步者走到“足够平”的地方就停下，而不是强求绝对谷底。

更深层的视角是 Krylov 子空间：CG 的每次迭代都在扩张子空间 Kₖ = span{r₀, Ar₀, …, A^{k−1}r₀}，并在其中寻找最优近似。这个子空间结构解释了为什么 CG 适合与预处理、重启配合使用。实际求解器中，CG 常用于有限元结构分析、CFD 压力泊松方程、机器学习中的正则化最小二乘。选择求解器时，工程师会先判断：矩阵是否对称正定？规模多大？稀疏模式怎样？然后决定用 CG 还是双共轭梯度稳定法（BiCGSTAB）、广义最小残差法（GMRES）等变体。

从雾谷里蒙眼的徒步者，到工程求解器中百万变量的幕后英雄，共轭梯度法始终贯彻一个朴素智慧：不重复走旧路，每次只做最便宜的计算。这把“有限步收敛”的数学之美与“一次矩阵向量乘”的工程之轻，完美地融合在了一起。

∥rk∥2∥b∥2≤ε,Kk(A,r0)=span{r0,Ar0,…,Ak−1r0}\displaystyle \frac{\|\mathbf{r}_k\|_2}{\|\mathbf{b}\|_2} \le \varepsilon, \qquad \mathcal{K}_k(A, \mathbf{r}_0) = \mathrm{span}\{\mathbf{r}_0, A\mathbf{r}_0, \ldots, A^{k-1}\mathbf{r}_0\}
工程中用相对残差判据控制精度；CG 的每一步都在 Krylov 子空间内寻找最优近似，预处理则让子空间更快“长满”。
  1. 相对残差分子是当前残差，分母是初始右端项，用来归一化精度。
    ∥rk∥2=∥b−Axk∥2\displaystyle \|\mathbf{r}_k\|_2 = \|\mathbf{b} - A\mathbf{x}_k\|_2
  2. 容差判据当相对残差小于 ε 时迭代停止，平衡精度与成本。
    ∥rk∥2≤ε∥b∥2\displaystyle \|\mathbf{r}_k\|_2 \le \varepsilon \|\mathbf{b}\|_2
  3. Krylov 子空间反复用 A 乘初始残差，生成的一系列向量张成搜索空间。
    Kk(A,r0)=span{r0,Ar0,…,Ak−1r0}\displaystyle \mathcal{K}_k(A, \mathbf{r}_0) = \mathrm{span}\{\mathbf{r}_0, A\mathbf{r}_0, \ldots, A^{k-1}\mathbf{r}_0\}

AI 应用场景

大规模图神经网络的隐式传播

阶段：模型训练与推理中的数值核心

解决的问题：图神经网络（GNN）常需计算归一化拉普拉斯矩阵 L 的隐式表示，如求解 (I + L)x = h。传统直接求逆在千万节点图上内存爆炸；AI 工程师希望只依赖矩阵向量乘，让内存和算力随边数线性增长。

具体案例：在商品推荐系统中，用户-物品二部图有 5000 万节点、2 亿条边。用预处理共轭梯度法求解拉普拉斯系统，设置相对残差 1e-5，结合雅可比预处理，单个 GPU 上数百轮迭代即可完成传播，而直接分解法根本无法装入显存。

物理仿真中有限元刚度方程

阶段：仿真数据生成（物理 AI 训练语料）

解决的问题：结构力学、流体仿真的核心是求解 K u = F，K 是千万自由度的稀疏对称正定刚度矩阵。传统 LU 分解需要 O(n²) 存储，在 HPC 上也不可扩展；CG 配合多重网格预处理成为主流。

具体案例：汽车碰撞仿真中，每帧需求解约 2000 万自由度的位移场。工程师采用共轭梯度法搭配不完全 Cholesky 预处理，在分布式集群上数百次迭代收敛，使仿真数据流水线能以近实时速度生成训练样本，供神经网络学习动力学。

大规模岭回归与高斯过程

阶段：模型训练中的线性系统求解

解决的问题：岭回归 (XᵀX + λI)w = Xᵀy 与高斯过程回归均涉及大矩阵求逆。直接求逆复杂度 O(n³)，而 CG 只需计算矩阵向量乘，配合预处理可把训练时间从数小时压缩到数分钟。

具体案例：广告点击率预估中，特征维度达 800 万、训练样本数亿。算法团队用 CG 迭代求解正则化方程，每次迭代用分块矩阵运算完成 Xᵀ(Xv) 与 λv 的累加，整个训练管线在 CPU 集群上稳定运行。

Python 演示

python68 行
import numpy as np  def conjugate_gradient(A, b, x0=None, tol=1e-8, max_iter=1000):    n = len(b)    if x0 is None:        x0 = np.zeros(n)    x = x0.copy()    r = b - A @ x    d = r.copy()    r_sq = r @ r    history = []    for k in range(max_iter):        history.append(x.copy())        if np.sqrt(r_sq) <= tol:            break        Ad = A @ d        alpha = r_sq / (d @ Ad)        x = x + alpha * d        r_new = r - alpha * Ad        r_sq_new = r_new @ r_new        beta = r_sq_new / r_sq        d = r_new + beta * d        r = r_new        r_sq = r_sq_new    return x, history  def steepest_descent(A, b, x0=None, tol=1e-8, max_iter=5000):    n = len(b)    if x0 is None:        x0 = np.zeros(n)    x = x0.copy()    r = b - A @ x    history = []    for k in range(max_iter):        history.append(x.copy())        r_sq = r @ r        if np.sqrt(r_sq) <= tol:            break        Ar = A @ r        alpha = r_sq / (r @ Ar)        x = x + alpha * r        r = b - A @ x    return x, history  if __name__ == "__main__":    A = np.array([[4.0, 1.0], [1.0, 3.0]])    b = np.array([1.0, 2.0])    x_cg, _ = conjugate_gradient(A, b)    x_sd, _ = steepest_descent(A, b)    x_true = np.linalg.solve(A, b)    print("CG 解:", x_cg)    print("SD 解:", x_sd)    print("真值:", x_true)     n = 200    A = (np.eye(n) * 2.0         + np.eye(n, k=1) * 0.5         + np.eye(n, k=-1) * 0.5)    bvec = np.ones(n)    x_cg, hist_cg = conjugate_gradient(A, bvec, tol=1e-6)    x_sd, hist_sd = steepest_descent(A, bvec, tol=1e-6)    res_cg = np.linalg.norm(A @ x_cg - bvec)    res_sd = np.linalg.norm(A @ x_sd - bvec)    print("n=200 CG 迭代数:", len(hist_cg), "残差:", res_cg)    print("n=200 SD 迭代数:", len(hist_sd), "残差:", res_sd)
目的：用纯 NumPy 实现共轭梯度法与最速下降法，对比两者在小型密集与较大稀疏系统上的收敛表现。
  1. r = b - A @ x计算当前残差向量，它是二次函数的负梯度，相当于登山杖探出的“最陡方向”。
  2. alpha = r_sq / (d @ Ad)沿搜索方向做精确一维搜索，确定这一步走多远。
  3. r_new = r - alpha * Ad用递推方式更新残差，避免每一步都重新计算 Ax，省下一次大矩阵乘法。
  4. beta = r_sq_new / r_sq计算新旧残差范数平方之比，用来构造与历史方向共轭的新方向。
  5. d = r_new + beta * d将当前残差与旧方向组合，形成不重复的搜索方向。
  6. x = x + alpha * d沿新方向迈出最优一步，更新近似解。
  7. history.append(x.copy())记录每一步的解，便于事后统计迭代次数和绘制收敛轨迹。

常见误区

误区CG 名字里带“梯度”，听起来像通用迭代法，新手常拿来直接求解非对称、不定或非方阵系统，还很自信地觉得能收敛。

正解共轭梯度法只对对称正定矩阵有理论保证：能量范数递减、最多 n 步收敛。对非对称矩阵，残差无法再用简单的三项递推保持共轭性，迭代可能振荡甚至发散。工程中遇到非对称问题应改用 BiCGSTAB、GMRES 等 Krylov 子空间变体。

为何常见：许多数值库里 CG 是默认迭代求解器之一，使用者容易忽略其底层对“对称正定”这一前提的依赖。

误区初学者看到“免存海森矩阵”“无需存储 A”等宣传语，误以为 CG 是零内存算法，连 A 本身都不需要。

正解CG 的核心优点是不显式存储和分解稠密矩阵，但每次迭代仍要计算 A 与向量的乘积。如果 A 是稀疏矩阵，可以只存非零元素或以算子形式隐式给出；如果 A 本身不便于向量乘，CG 的优势就不复存在。

为何常见：把“不存稠密矩阵”过度解读为“不需要矩阵”，忽略了矩阵向量乘法是算法真正的生命线。

误区以为 n 步保证收敛等价于一定很快，直接用百万维矩阵跑裸 CG，发现迭代次数远超 n 或迟迟不收敛，便怀疑算法有 bug。

正解n 步收敛是在精确算术下成立；浮点舍入误差会破坏共轭关系，实际往往需要更多迭代。而且每步的矩阵向量乘在 n 很大时并不便宜。正确做法是配合预处理、重启策略和合理的终止容差。

为何常见：只背诵定理结论，忽略了数值计算中舍入误差与计算复杂度的现实约束。

误区认为共轭梯度法必须先一次性构造出所有共轭方向，再逐条走完，等同于共轭方向法。

正解经典共轭方向法需预先知道整个共轭方向组，代价接近获取 Hessian 信息。CG 的价值在于用β递推逐条生成方向，每次只需上一条方向和当前残差，内存和计算都大幅减少。“动态生成”才是 CG 的关键创新。

为何常见：从名字“共轭梯度”望文生义，混淆了“事先构造基”与“迭代递推”的两种算法。

面试题

Q1. 共轭梯度法求解线性方程组时，矩阵 A 必须满足什么条件？

单选题 · easy

  • A.可逆且行列式非零
  • B.对称正定矩阵← 正确答案
  • C.大型稀疏非对称阵
  • D.主对角线元素全为1

答案：B。CG 的理论基础依赖 A 的对称正定性：二次函数严格凸、所有方向保持 A-共轭、能量范数单调下降。普通可逆矩阵或非对称矩阵都无法保证这些性质，因此 B 正确。

Q2. 两条方向 d₁ 和 d₂ 称为 A-共轭，其数学条件是？

单选题 · easy

  • A.d₁ᵀd₂ = 0
  • B.d₁ᵀ A d₂ = 0← 正确答案
  • C.A d₁ = d₂
  • D.d₁ 与 d₂ 长度相等

答案：B。A-共轭本质是在 A 定义的内积下正交，即 ⟨d₁,d₂⟩_A = d₁ᵀAd₂ = 0。普通正交只是 A 为单位矩阵时的特例，所以选 B。

Q3. 最速下降法中，相邻两步搜索方向满足什么关系？

单选题 · mid

  • A.关于 A 共轭
  • B.普通欧氏内积为零← 正确答案
  • C.线性相关且同向
  • D.与 A 的特征向量平行

答案：B。精确一维搜索使 f 沿当前方向达到极小，意味着新梯度与旧搜索方向正交，即 rₖ₊₁ᵀrₖ = 0。因此相邻两步在欧氏空间正交，这恰恰是锯齿现象的来源。

Q4. CG 中 βₖ 的主要作用是？

单选题 · mid

  • A.决定迭代步长大小
  • B.构造与历史方向共轭的新方向← 正确答案
  • C.修正目标函数的常数项
  • D.消除舍入误差影响

答案：B。βₖ 把当前残差与上一条方向组合成新方向 dₖ₊₁，并通过精心设计的公式保证新方向与之前所有方向关于 A 共轭。步长由 αₖ 决定，因此 B 正确。

Q5. 在精确算术下，n 维对称正定系统用 CG 最多迭代多少次收敛？

单选题 · easy

  • A.n 次← 正确答案
  • B.n² 次
  • C.2n 次
  • D.无限次

答案：A。共轭方向法在 n 维空间中最多需要 n 个两两 A-共轭的方向即可张成整个空间，CG 每次迭代产生一个新共轭方向，因此至多 n 步抵达精确解。这是 CG 家族的核心性质。

Q6. 实际大规模计算中，CG 的“n 步收敛”性质为何会失效？

单选题 · extreme

  • A.矩阵 A 的元素太大导致溢出
  • B.浮点舍入误差破坏共轭关系← 正确答案
  • C.预处理矩阵不可逆
  • D.终止容差设置过大

答案：B。n 步终止是精确算术下的结论。浮点运算中，残差与方向的共轭关系会被舍入误差逐渐破坏，导致实际迭代超过 n 步。预处理、重启策略和合理容差都是工程上应对这一现象的手段。

Q7. 相比直接分解法，CG 的工程优势包括哪些？

多选题 · mid

  • A.无需显式存储完整稠密矩阵← 正确答案
  • B.单次迭代只需一次矩阵向量乘← 正确答案
  • C.内存随非零元数量近似线性增长← 正确答案
  • D.可以无条件下处理非对称系统

答案：ABC。CG 不构造稠密分解，只以矩阵向量乘形式接触 A，因此内存和单步成本都低，非常适合大规模稀疏问题。但它要求对称正定，D 把适用范围夸大了，故排除。

Q8. 预处理加速 CG 的机理包括哪些？

多选题 · hard

  • A.选择 M≈A 使 M⁻¹A 特征值聚集← 正确答案
  • B.将原方程改写为 M⁻¹Ax = M⁻¹b← 正确答案
  • C.使等效条件数更接近 1← 正确答案
  • D.将把条件数提升到最大以加速下降

答案：ABC。预处理的目标是让新系统的特征值谱更集中、条件数更小，从而加快线性收敛速率。M⁻¹A 越接近单位矩阵，收敛越快。D 说法与预处理目标完全相反。

Q9. CG 免于存储海森矩阵的关键在于？

单选题 · mid

  • A.只使用矩阵向量乘积 A·d← 正确答案
  • B.对矩阵做 LU 分解
  • C.提前求出 A 的逆矩阵
  • D.用随机采样替代确定性计算

答案：A。CG 的全部迭代只需要计算 A 与方向向量的乘积，算法本身不创建、不分解、不存储完整的 A 或海森矩阵。正因为只依赖 A·d 操作，稀疏矩阵算子或隐式向量乘都能接入。

Q10. 关于非线性共轭梯度法，下列说法正确的有？

多选题 · hard

  • A.用梯度 gₖ 替代线性 CG 中的残差← 正确答案
  • B.目标函数允许是非二次的← 正确答案
  • C.FR 公式用相邻梯度范数平方之比计算 β← 正确答案
  • D.收敛行为与线性 CG 在任何条件下完全相同

答案：ABC。非线性 CG 沿用线性 CG 的记账思想，把残差换为梯度，并用线搜索确定步长，适用于一般光滑优化。FR 公式 β = ‖gₖ₊₁‖²/‖gₖ‖² 是典型实现。但由于目标非二次，共轭性只能近似，收敛行为与线性情形并不完全相同，故 D 错误。