cg迭代法(Conjugate gradient method)求解线性方程组

Author : zbzhen,        Modified : Sat Sep 9 14:25:34 2023

1. 模型问题

矩阵$A$为满秩对称矩阵，求线性方程组

$$Ax=b$$

2. 分析与推导

设$x_0=0$, 如果非零向量组$\{P_0,P_1,P_2,\dots\}$线性无关, 当$k+1$等于向量$x_{k+1}$的维数时，则向量$x_{k+1}$可写成

$$x_{k+1}=\sum_{j=0}^k \alpha_jP_j$$

其中$\alpha_j$为待定标量系数. 或者可写成迭代格式

$$x_{k+1}=x_k+ \alpha_k P_k$$

记残差

$$r_k = b - A x_k$$

联立上面两式可得

$$r_{k+1} = r_k - \alpha_k AP_k$$

或可改写成

$$AP_k = \dfrac{1}{\alpha_k} (r_k-r_{k+1})$$

接下来就是需要想办法找到合适的$\alpha_k$与$P_k$使得计算复杂度尽可能少

取$P_0=r_0$, 同时把$P_{k}$写成$r_0,r_1,r_2,\cdots,r_k$的线性组合， 使得

$$P_j^TAP_k = 0, \quad j\neq k.$$

也可写成递推格式，

$$P_{k+1}=r_{k+1}+ \beta_k P_k$$

可选取合适的$\beta_k$使得

$$r_k^T r_j = 0,\quad k\neq j.$$

于是联立上面四个等式，在上面倒数第二式左乘$P_k^TA$, 然后由正交性可得

$$\begin{aligned} \beta_k &= -\dfrac{P_k^TAr_{k+1}}{P_k^TAP_k} = -\dfrac{r_{k+1}^TAP_k}{P_k^TAP_k} \\&= -\dfrac{r_{k+1}^T(r_{k}-{r}_{k+1})}{\alpha_k P_k^TAP_k} = \dfrac{r_{k+1}^Tr_{k+1}}{\alpha_k P_k^TAP_k} \end{aligned}$$

最后就是$\alpha_k$的推导. 根据

$$\begin{aligned} P_k^T A P_k &= (r_{k}+ \beta_{k-1} P_{k-1})^T A P_k = r_k^TAP_k \\&= \dfrac{1}{\alpha_k}r_k^T (r_k-r_{k+1}) = \dfrac{1}{\alpha_k} r_k^T r_k \end{aligned}$$

即得

$$\alpha_k = \dfrac{r_k^T r_k }{P_k^T A P_k}$$

进而得到

$$\beta_k = \dfrac{r_{k+1}^Tr_{k+1}}{r_k^T r_k}$$

需要注意的是，上面的公式推导过程似乎并没有用到完全正交性， 只是用到了
$P_{k+1}^TAP_k=r_{k+1}^Tr_k=0$
但事实上， 根据两个递推公式($P_0=r_0$)
$P_{k+1}=r_{k+1}+ \beta_k P_k$ 与 $r_{k+1} = r_k - \alpha_k AP_k$
可以得到完全正交公式$P_j^TAP_k = r_k^T r_j = 0, \quad j\neq k$

3. 算法

3.1. 便于理解的粗糙算法

给定$A,b,x_0$, 求出$x=A^{-1}b$

初始化：$r_0=b - Ax_0$, $P_0=r_0$, $k=0$

如果没有收敛：

• $\alpha_k = \dfrac{r_k^T r_k }{P_k^T A P_k}$

• $x_{k+1}=x_k+ \alpha_k P_k$

• $r_{k+1} = r_k - \alpha_k AP_k$

• $\beta_k = \dfrac{r_{k+1}^Tr_{k+1}}{r_k^T r_k}$

• $P_{k+1}=r_{k+1}+ \beta_k P_k$

• $k \leftarrow k+1$

3.2. 优化后的算法

为了使得不重复计算$AP_k$以及$r_k^Tr_k$,因此算法可以优化为:

给定$A,b,x_0$, 求出$x=A^{-1}b$

初始化：$r_0=b - Ax_0$, $P_0=r_0$, $k=0$, $s_0=r_0^Tr_0$

如果没有收敛：

• $d_k = AP_k$
• $\alpha_k = \dfrac{s_k }{P_k^T d_k}$
• $x_{k+1}=x_k+ \alpha_k P_k$
• $r_{k+1} = r_k - \alpha_k d_k$
• $s_{k+1} = r_{k+1}^Tr_{k+1}$
• $P_{k+1}=r_{k+1}+ \dfrac{s_{k+1}}{s_k} P_k$
• $k \leftarrow k+1$

4. MATLAB / GNU Octave 程序实现

function x = conjgrad(A, b, x)
    r = b - A * x;
    p = r;
    rsold = r' * r;

    for i = 1:length(b)
        Ap = A * p;
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        if sqrt(rsnew) < 1e-10
            break
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew;
    end
end

5. python实现

5.1. 基本cg程序

import numpy as np
def cg(A, b, x):
    r = b - A @ x
    p = r
    rsold = r @ r
    for i in range(len(b)):
        Ap = A @ p
        alpha = rsold / (p @ Ap)
        x += alpha * p
        r -= alpha * Ap
        rsnew = r @ r
        if np.sqrt(rsnew) < 1e-10:
            break
        p = r + (rsnew / rsold) * p
        rsold = rsnew
    print(np.sqrt(rsnew), i)
    return x

5.2. 预处理cg程序

def pcg(A, b, x, Minv):
    r = b - A @ x
    z = Minv @ r
    p = z
    rsold = z @ r
    for i in range(len(b)):
        Ap = A @ p
        alpha = rsold / (p @ Ap)
        x += alpha * p
        r -= alpha * Ap
        z = Minv @ r
        rsnew = z @ r
        if np.sqrt(rsnew) < 1e-10:
            break
        p = z + (rsnew / rsold) * p
        rsold = rsnew
    print(np.sqrt(rsnew), i)
    return x