参数设置
A 为对称三对角 4×4 矩阵。固定元素: A[1][2]=A[2][1]=1, A[2][3]=A[3][2]=1, A[0][1]=A[1][0]=1, A[2][2]=2, A[3][3]=1。b 的剩余分量: b[1]=2, b[2]=3, b[3]=4。初值 x_0=0。
A =
b =
残差范数收敛历史 log10(‖r_k‖_2)
蓝色实线＝CG 法／红色实线＝最速下降法／横轴＝迭代步 k／纵轴＝log10(‖r_k‖_2)
理论与主要公式
CG 法在更新时保持残差 $r_k = b - A x_k$ 与搜索方向 $p_k$ 在 A-意义下两两正交（共轭）。
初始化: $r_0 = b - A x_0$, $p_0 = r_0$。迭代步:
$$\alpha_k = \frac{r_k^\top r_k}{p_k^\top A p_k}, \quad x_{k+1} = x_k + \alpha_k p_k$$
$$r_{k+1} = r_k - \alpha_k A p_k, \quad \beta_k = \frac{r_{k+1}^\top r_{k+1}}{r_k^\top r_k}$$
$$p_{k+1} = r_{k+1} + \beta_k p_k$$
在理想算术下，n 维对称正定方程组最多 n 步即可精确求解。误差以 $\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^k$ 的速率衰减，其中 $\kappa = \lambda_\max/\lambda_\min$ 为条件数。
共轭梯度法模拟器是什么
🙋
我在 FEM 教材里反复看到「用 CG 法求解刚度方程 Kd=f」，CG 法到底是什么？
🎓
简单说，就是对对称正定矩阵 A 迭代求解 Ax=b 的方法。每一步先算残差 $r_k = b - A x_k$，再构造一个与过去所有方向都 A-正交（共轭）的搜索方向 $p_k$，最后沿这个方向走最优步长 $\alpha_k$。把模拟器里的「迭代次数」从 1 拖到 10，看蓝色曲线在对数坐标下急剧下降。
🎓
那是最速下降法，每次只看「当前残差方向」往下走的原始版本。把 A[0][0] 调到 10、A[1][1] 调到 1 来增大条件数，你就会看到最速下降法在山谷里来回摆动几乎不收敛。CG 法绝不会重复探索同一个方向，所以对这个 4 维问题只用 4 步就能精确命中真解。
🙋
真的！第 4 步时蓝点突然降到残差 1e-15 左右，这是巧合吗？
🎓
不是巧合，是定理。对 $n$ 维对称正定系统，CG 法在理想算术下最多 $n$ 步就能精确求解，因为 Krylov 子空间 $\mathrm{span}\{r_0, A r_0, A^2 r_0, \ldots\}$ 维数饱和到 n 时一定包含真解。浮点算术下会略多几步，但相对最速下降法仍是数量级的差距。
🙋
条件数 cond(A) 的卡片显示 11 左右，这个数字代表什么？
🎓
条件数是最大特征值与最小特征值的比 $\kappa = \lambda_\max / \lambda_\min$。$\kappa$ 越大，等高线就越像细长的山谷，最速下降法会被困住。CG 法的误差按 $\left(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^k$ 衰减——注意里面是平方根，因此即使 $\kappa=100$，CG 也比最速下降法快约十倍。实际工程里几乎都用「预条件 CG (PCG)」来进一步降低有效条件数。
常见问题
A 必须对称（A^T = A）且正定（任意非零向量 x 都有 x^T A x > 0）。实际中 FEM 刚度矩阵、拉普拉斯矩阵、协方差矩阵等天然满足。对非对称矩阵可用 Bi-CGSTAB、GMRES，对对称不定矩阵可用 MINRES、SYMMLQ，它们都是 CG 思想的近亲。
CG 只需要矩阵向量乘 matVec(A, v)，因此能完整保持稀疏结构，永远不显式构造分解矩阵。直接法分解后会出现「填入 (fill-in)」，对数百万自由度的 FEM 内存就爆掉，所以预条件 CG (PCG) 才是事实标准。但当多个右端项共用一个 A 时，复用 LU 分解可能比迭代法更快。
M 应当近似 A 而 M^{-1}v 计算便宜。Jacobi（对角缩放）最简单但效果有限；不完全 Cholesky 分解 (IC(0), ICT) 是结构分析的主力；流体或复杂物理推荐代数多重网格 (AMG)，能取得与条件数无关的 O(n) 收敛。Trilinos、PETSc 等库提供了大量现成的预条件子。
最常用的是相对残差 ‖r_k‖_2 / ‖b‖_2 < tol，tol 取 1e-6 到 1e-10。条件数大时残差小不代表真误差 ‖x_k - x_*‖ 小，因此还需结合 A-范数误差或有效位数估计。务必再加一个最大迭代上限（例如 2n 或 500-2000），若收敛停滞则重新设计预条件子。
现实应用
有限元结构分析 (FEA): 汽车、航空、建筑的应力分析会得到节点位移为未知数的大型刚度方程 Kd=f。K 在施加约束后是对称正定且高度稀疏，自由度可达数百万。预条件 CG（IC-PCG、AMG-PCG）是 Abaqus、Ansys、Nastran 等商用求解器以及 FreeFEM、FEniCS 等开源框架的主力解法。
偏微分方程的离散化: 用有限差分、有限体积法离散热传导、Poisson 方程、椭圆型 PDE，会得到对称正定方程组。多重网格与 CG 的组合可在地球流体力学、气象模拟中实现与条件数无关的 O(n) 可扩展性。
机器学习与优化: 二次型最小化 min (1/2)x^T A x - b^T x 等价于求解 Ax=b，因此 CG 本质上是一种优化算法。在「Hessian-free 优化」中，CG 不显式构造 Hessian 即可求解 Newton 步长，被广泛用于深度学习与强化学习。
图像处理与断层成像: CT 重建、图像复原最终归结为求解正规方程 (A^T A) x = A^T b。由于 A^T A 是对称半正定（加正则化后正定），CG 的衍生算法 (CGLS、LSQR) 是该领域内存占用小、速度快的标准工具。
常见误区与注意事项
最常见的误区是无条件相信「CG 一定在 n 步达到精确解」。这是「理想算术」下的结论，在浮点运算中舍入误差会逐步破坏共轭性 $p_i^\top A p_j = 0$，n 步之后仍残留一小段非零残差。实际中以相对残差容差作停止准则，必要时使用重启 (restart) 或再正交化。模拟器是 4 维小问题几乎看不到舍入噪声，但在百万维实际问题中表现高度依赖条件数。
第二个陷阱是条件数大就放弃 CG 法。诚然条件数 $\kappa$ 越大收敛越慢，但选个好的预条件子 M 就能彻底改写局面。生产环境几乎不用裸 CG——不完全 Cholesky (IC)、对称 Gauss-Seidel (SSOR)、代数多重网格 (AMG) 都能把实际条件数压到 $\sqrt{\kappa}$ 甚至 $O(1)$。大多数「CG 慢」的情况其实是「预条件子选错了」。在模拟器里把 A[0][0]=10、A[1][1]=1 拉开差距，已经能看到 CG 与最速下降法的差距，好的预条件子还会把差距进一步拉大到只需 1-2 步。
最后一点，切勿对非对称或不定系统使用 CG。CG 是为对称正定矩阵设计的，对非对称矩阵（如含对流项的流体系统）盲目使用，$p_k^\top A p_k$ 可能变为负值或零导致发散。非对称系统应当使用 Bi-CGSTAB、GMRES、IDR(s)，对称不定则用 MINRES、SYMMLQ。把 CG 严格视为「SPD 专用最优解法」，是安全实现的第一步。