共轭梯度法

共轭梯度法（conjugate gradient method，CG）是求解实对称正定线性系统的 Krylov 子空间迭代法。它沿相互 A 共轭的方向逐次修正，使二次能量在不断扩大的子空间内最小。每轮主要需要一次矩阵向量积，适合大型稀疏系统；非对称或不定矩阵不满足本方法的基本假设。

两步解出一个二阶系统

取 $$A=\left(\begin{array}{cc}4& 1\\ 1& 3\end{array}\right),b=(1,2{)}^T,x_0=0.$$ 初始残差 $r_0=b-A{x}_0=(1,2{)}^T$，首方向 $p_0=r_0$。因为 $A{p}_0=(6,7{)}^T$，步长为 ${\alpha }_0=(r_0^T{r}_0)/(p_0^{T}A{p}_0)=5/20=1/4$。故 $$x_1=(1/4,1/2{)}^T,r_1=(-1/2,1/4{)}^T.$$ 新残差与旧残差的内积为零。取 ${\beta }_0=(r_1^T{r}_1)/(r_0^T{r}_0)=1/16$，新方向为 $p_1=r_1+{\beta }_0{p}_0=(-7/16,3/8{)}^T$。计算 $A{p}_1=(-11/8,11/16{)}^T$，所以 $p_1^{T}A{p}_1=55/64$、${\alpha }_1=4/11$。第二步得到 $$x_2=x_1+{\alpha }_1{p}_1=(1/11,7/11{)}^T.$$ 代回两行方程分别为一和二，残差为零。这里只用了两步，来自精确算术下二维空间的有限终止性质，不意味着所有大问题都只要少量步数。

二次能量与共轭方向

对称正定性使 $\varphi (x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$ 有唯一极小点 $x_{\ast }=A^{-1}b$，因为 $$\varphi (x)-\varphi (x_{\ast })=\frac{1}{2}(x-x_{\ast }{)}^{T}A(x-x_{\ast }).$$ 因此最小化能量等价于最小化 A 范数误差。方向 $p_i,p_j$ 若满足 $p_i^{T}A{p}_j=0$，称为 A 共轭，它们通常不是普通欧氏正交。

沿方向 $p_k$ 令 $x_{k+1}=x_k+\alpha p_k$，对 $\varphi$ 求导得到 $-p_k^T{r}_k+\alpha p_k^{T}A{p}_k=0$。在 CG 的正交关系下 $p_k^T{r}_k=r_k^T{r}_k$，所以 $${\alpha }_k=\frac{r_k^T{r}_k}{p_k^{T}A{p}_k},r_{k+1}=r_k-{\alpha }_{k}A{p}_k,p_{k+1}=r_{k+1}+{\beta }_k{p}_k,{\beta }_k=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}.$$ SciPy 的接口明确要求矩阵为 Hermitian 正定，并以残差容限决定是否收敛。

两个方向在 A 内积下共轭。首步点已经沿第一方向最小化能量，第二方向不会破坏此前的最小化条件。

正交递推与有限终止的证明

可用归纳法证明：到第 $k$ 步，方向 $p_0,\dots ,p_k$ 两两 A 共轭，残差 $r_k$ 与此前方向正交。新残差满足 $$p_k^T{r}_{k+1}=p_k^T{r}_k-{\alpha }_k{p}_k^{T}A{p}_k=0.$$ 对更早的 $p_j$，新残差与其内积等于 $p_j^T{r}_k-{\alpha }_k{p}_j^{T}A{p}_k=0$。由于此前残差由此前方向线性组合而成，新残差也与此前残差正交。

要让新方向与当前方向共轭，应取 ${\beta }_k=-r_{k+1}^{T}A{p}_k/(p_k^{T}A{p}_k)$。利用 $A{p}_k=(r_k-r_{k+1})/{\alpha }_k$ 和残差正交，分子为 $-\Vert r_{k+1}{\Vert }^2/{\alpha }_k$，代入便得到上述平方残差比。对 $j<k$，$A{p}_j=(r_j-r_{j+1})/{\alpha }_j$，新残差与两项均正交，所以新方向也与所有早期方向共轭。归纳闭合。

每个非零共轭方向线性独立：若线性组合为零，与其中一个方向作 A 内积，正定性迫使它的系数为零。因此在 $n$ 维精确算术下最多产生 $n$ 个方向。此时残差与整个空间正交，只能为零。相同正交关系还说明 $x_k$ 在 $x_0+{𝒦}_k(A,r_0)$ 上最小化能量，其中 ${𝒦}_k=\mathrm{span}\{r_0,A{r}_0,\dots ,A^{k-1}{r}_0\}$。

浮点运算会破坏精确正交和有限终止。条件数大时，实际步数与残差核验尤其重要。经典误差界 $\Vert e_k{\Vert }_A\le 2[(\sqrt{\kappa }-1)/(\sqrt{\kappa }+1){]}^k\Vert e_0{\Vert }_A$ 使用 $\kappa ={\lambda }_{max}/{\lambda }_{min}$；特征值聚簇时实际表现还可能更好。预条件子正是通过改变有效几何来改善这些性质。

可运行的稠密演示

代码输入实对称正定矩阵、右端、容限和最大步数，返回近似及迭代数。它检查对称性和迭代中出现的正曲率，但不以这些有限检查替代完整正定性证明。演示使用稠密乘法，每步 $O(n^2)$ 时间，工作向量为 $O(n)$；稀疏存储时乘法可降为 $O(\mathrm{nnz}(A))$。不修改输入。

Python 3

import math


def cg(A, b, tol=1e-12, maxiter=None):
    A = [list(map(float, row)) for row in A]
    b = list(map(float, b)); n = len(b)
    if not n or len(A) != n or any(len(row) != n for row in A):
        raise ValueError("square matrix matching nonempty b required")
    if not all(math.isfinite(v) for row in A for v in row) or not all(math.isfinite(v) for v in b):
        raise ValueError("finite data required")
    if any(A[i][j] != A[j][i] for i in range(n) for j in range(i)):
        raise ValueError("symmetric matrix required")
    if not math.isfinite(tol) or tol <= 0:
        raise ValueError("positive finite tolerance required")
    maxiter = n if maxiter is None else maxiter
    if not isinstance(maxiter, int) or maxiter < 0:
        raise ValueError("nonnegative integer maxiter required")
    dot = lambda u, v: sum(a*c for a, c in zip(u, v))
    mv = lambda p: [dot(row, p) for row in A]
    x = [0.0]*n; r = b[:]; p = r[:]; rr = dot(r, r)
    if not math.isfinite(rr):
        raise OverflowError("initial residual square outside floating-point range")
    target = tol*max(1.0, math.sqrt(rr))
    if not math.isfinite(target):
        raise OverflowError("residual threshold outside floating-point range")
    for k in range(maxiter+1):
        actual = [v-w for v, w in zip(b, mv(x))]
        if math.hypot(*actual) <= target:
            return x, k
        if k == maxiter:
            break
        Ap = mv(p); denominator = dot(p, Ap)
        if not math.isfinite(denominator) or denominator <= 0 or not math.isfinite(rr):
            raise ArithmeticError("breakdown or nonpositive curvature")
        alpha = rr/denominator
        x = [v+alpha*d for v, d in zip(x, p)]
        r = [v-alpha*d for v, d in zip(r, Ap)]
        rr_new = dot(r, r)
        if not math.isfinite(rr_new) or not all(math.isfinite(v) for v in x):
            raise OverflowError("nonfinite iteration")
        beta = rr_new/rr
        p = [v+beta*d for v, d in zip(r, p)]; rr = rr_new
    raise RuntimeError("true residual did not reach tolerance")


print(cg([[4, 1], [1, 3]], [1, 2]))  # ([1/11, 7/11], 2), up to rounding

有解变式取 $A=\mathrm{diag}(2,5)$、$b=(2,0{)}^T$。初始残差已经是一个特征向量，首步 $\alpha =1/2$ 就到达 $x=(1,0{)}^T$。若取不定矩阵 $\mathrm{diag}(1,-1)$ 和 $b=(1,1{)}^T$，第一步分母为零，直接暴露正定假设的作用；这时应选与不定系统匹配的方法。