共轭梯度法

在数值分析中，共轭梯度法（Conjugate Gradient Method, CG 法）是当今求解对称正定的线性方程组的首选数值方法，从格式上来看它是一种迭代方法，在迭代到$n$步之后，不计机器误差的情况下可以得到精确解。

方法

设我们要解决的线性方程组为

$$Ax = b, \quad A \in \R^{n \times n}, b \in \R^n.$$

设$x_0$为一给定的迭代初始向量，$r_0 = b - Ax_0$是初始残量，取共轭方向$p_0 = r_0$，按照如下迭代格式进行迭代

$$\left\{ \begin{align} \alpha_k & = \dfrac{(r_k, r_k)}{(p_k, p_k)_A}, \\ x_{k+1} & = x_k + \alpha_k p_k, \\ r_{k+1} & = r_k - \alpha_k A p_k, \\ \beta_k & = \dfrac{(r_{k+1}, r_{k+1})}{(r_k, r_k)}, \\ p_{k+1} & = r_{k+1} + \beta_k p_k. \end{align}\right.$$

其中，$(,)$是普通内积，$(,)_A$是$A$内积，即$(\alpha, \beta)_A = (A\alpha, \beta).$

推导过程从略，这里仅指出一些特点：

1. 和完全正交方法一样，CG 法选择的约束空间和搜索空间是一样的，因此是一种正交投影法，但它是考察了对称矩阵的特殊结构，简化了 FOM 方法。
2. CG 方法中我们设定的残量$r_0, r_1, \cdots, r_k, \cdots$相互正交。
3. $p_k$被称为搜索方向，这些方向彼此$A$正交，即$(Ap_i, p_j) = 0, \forall i \ne j.$且当$i > n$时$p_i = 0.$
4. $\alpha_k$是由下式待定的
   $$x_{k+1} = x_k + \alpha_k p_k.$$
   它描述了第$k$次迭代解向精确解靠近的$A$共轭方向上的步长。
5. 未得到精确解时，第$k$步 CG 法得到的向量总是在以下意义下最优：
   $$x_k = \arg \min_{x \in x_0 + K} \| x - x^* \|_A$$
   其中$x^*$是精确解，Krylov 子空间$K = \text{span} \{ r_0, Ar_0, \cdots, A^{k-1} r_0 \}.$

算法

以下程序设计是 CG 法实现的简易 Matlab 代码。

function [z, k] = cg(A, b, x0, eps)
% 共轭梯度法 (CG法) 求解对称正定矩阵的线性方程组 Ax = b 的解的程序.
% A 是系数矩阵, 必须对称正定, b 是常数项矩阵, x0 是迭代解的初始值, eps 是相对误差, 可缺省.
if ~exist('eps', 'var')
    eps = 10e-6;
end
r(:,1) = b - A * x0;
x(:,1) = x0;
q(:,1) = r(:,1);
[~, n] = size(A);
for k=1:n
    r(:,2) = r(:,1); x(:,2) = x(:,1); q(:,2) = q(:,1);
    alpha = dot(r(:,2), r(:,2))/dot(A*q(:,2), q(:,2));
    x(:,1) = x(:,2) + alpha * q(:,2);
    r(:,1) = r(:,2) - alpha * A * q(:,2);
    beta = dot(r(:,1), r(:,1))/dot(r(:,2), r(:,2));
    q(:,1) = r(:,1) + beta * q(:,2);
    if norm(r(:,2))/norm(r(:,1)) < eps
        break;
    end
end
z = x(:,1);
end

参考资料

1. 黄云清, 《数值计算方法》, 科学出版社, 北京, 2012-06, ISBN 978-7-0302-3428-5.