求解 Ax=b —— 共轭梯度法(上)

为了进行线性方程组 $Ax=b$ 的数值求解，这里介绍一种数值求解方法 —— 共轭梯度法（Conjugate gradient method），简称 CG。 共轭梯度法适用于矩阵 $A$ 是对称正定的情况，如果不考虑计算机处理浮点数的误差，该方法能够给出线性方程组的精确解。

共轭梯度法的条件

为了求解线性方程组(或者说线性系统) $Ax=b$，其中 $A$ 是一个 $n\times n$ 的矩阵，即 $A \in R^{n \times n}$；$b \in R^n$；$x \in R^n$。 此外，还要求矩阵 $A$ 是对称正定(symmetric position definite，SPD)的，这样才能保证最后共轭梯度法的迭代会收敛到所要求的解。 所谓对称，是指 $A^T = A$，即矩阵 $A$ 的转置等于其本身。 所谓正定，是指 $\forall x \neq 0$，都有 $x^TAx>0$。
可以总结，为了使用共轭梯度法，对矩阵 $A$ 有如下对约束条件：

1. 矩阵$A$必须是对称矩阵；
2. 矩阵$A$必须是正定矩阵；
3. 矩阵$A$必须是$n\times n$ 的方阵（其中这条已经隐含在条件1和2里面了，因为矩阵的对称和正定都要求矩阵是方阵）；

共轭梯度法迭代方法

共轭梯度法迭代方法还是相对比较简单的，下面直接给出其算法：

$r_0 = b -Ax_0$
$p_0 = r_0$
$\blue{\rm{for}} \; i=0,1,2, \cdots$
$~~~~$ $\alpha_{i}=r_{i}^{T} r_{i} / p_{i}^{T} A p_{i}$
$~~~~$ $x_{i+1}=x_{i}+\alpha_{i} p_{i}$
$~~~~$ $r_{i+1}=r_{i}-\alpha_{i} A p_{i}$
$~~~~$ $\beta_{i+1}=r_{i+1}^{T} r_{i+1} / r_{i}^{T} r_{i}$
$~~~~$ $p_{i+1}=r_{i+1}+\beta_{i+1} p_{i}$
$\blue{\rm{endfor}}$

这里，在每一轮 for 循环的迭代中，只需要更新向量 $r$ 和 $q$ 的值即可。而且在计算时也只需要存储几个向量及矩阵$A$本身即可，对内存需求并不大。
这里的 $\alpha$ 是迭代搜索步长，$p$ 是搜索方向；每一轮迭代中，可以通过步长和搜索方向来找到一个更新的解：即 $x_{i+1}=x_{i}+\alpha_{i} p_{i}$。 有了新的解后，需要更新残留 $r$（残留可以理解为当前的解 $x$ 离精确解的误差）；最后计算 $\beta$ 并以此更新搜索方向 $p$。

CG 迭代求解举例

我们以下面这个线性方程组 $Ax=b$ 的求解为例，来更清晰地说明 GC 方法的求解步骤。

$$\begin{pmatrix} 10 & -1 & 2 & 0 \\ -1 & 11 & -1 & 3 \\ 2 & -1 & 10 & -1 \\ 0 & 3 & -1 & 8 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} = \begin{pmatrix} 6 \\ 25 \\ -11 \\ 15 \end{pmatrix}$$

可以发现这个矩阵 $A$ 是对称的，同时，也可以证明这个矩阵是正定的。
给定初始值$x_0 = \left(0,0,0,0\right)^T$，迭代过程中 $x_i$ 的值记录如下：
(其中，$x_{ij}$ 表示第 $i$ 次迭代后，$x_j$ 的值，如 $x_{42} = 2$ 表示第4次迭代结束后$x_2$的值为2)

	$x_{i1}$	$x_{i2}$	$x_{i3}$	$x_{i4}$
$i=0$	0	0	0	0
$i=1$	0.4716	1.9651	-0.8646	1.1791
$i=2$	0.9964	1.9766	-0.9098	1.0976
$i=3$	1.0015	1.9833	-1.0099	1.0197
$i=4$	1.0000	2.0000	-1.0000	1.0000

可以看到，这里迭代4步后就收敛到所期望的解了: $x^*=\left(1,2,-1,1\right)^T$，实际上这个解也是该线性方程组的解析解。 可以证明，在 CG 方法中，对于 $n\times n$ 的对称正定矩阵，经过 $n$ 次迭代，必然会收敛到该线性方程组的解析解，这点将在后面讨论。

针对该例，对应的 CG 算法的 C 和 Go 语言实现如下:

C

main.c

#include<stdio.h>
#include<string.h>

// return inner product of two vector a and b of lenght len: a*b
double inner_product(double *a, double *b, size_t len) {
    double r = 0;
    for (size_t i = 0; i < len; i++) {
        r += a[i] * b[i];
    }
    return r;
}

// return product of matrix A and vector b of size len: A*b
// the result will be saved in vector c.
double *matrix_vec_product(double *A, double *b, double *c, size_t len) {
    for (size_t i = 0; i < len; i++) {
        c[i] = 0.0;
        for (size_t j = 0; j < len; j++) {
            c[i] += A[i * len + j] * b[j];
        }
    }
    return c;
}

// return a + b, where a and b are both vector of both size len.
// results is saved in c.
double *vec_add(double *a, double *b, double *c, size_t len) {
    for (size_t i = 0; i < len; i++) {
        c[i] = a[i] + b[i];
    }
    return c;
}

// return c*a, where a is a vector of size len. and c is a real number
// result is saved in vector b.
double *vec_multiple(double *a, double *b, double c, size_t len) {
    for (size_t i = 0; i < len; i++) {
        b[i] = c * a[i];
    }
    return b;
}

void print_vec(double *a, size_t len) {
    for (size_t i = 0; i < len; i++) {
        printf("%f ", a[i]);
    }
    printf("\n");
}

int main() {
    double A[4][4] = {
            {10, -1, 2,  0},
            {-1, 11, -1, 3},
            {2,  -1, 10, -1},
            {0,  3,  -1, 8}};
    double b[4] = {6, 25, -11, 15};
    int n = 4;

    // init
    double x[] = {0.0, 0.0, 0.0, 0.0};
    double r[n], p[n];
    vec_add(b, vec_multiple(matrix_vec_product(A[0], x, r, n), r, -1.0, n), r, n); // r = b-Ax
    memcpy(p, r, sizeof(double) * n); // p <- r
    // iteration
    double temp[n]; // temp vector
    for (int i = 0; i < n; i++) {
        double alpha = inner_product(r, r, n) / inner_product(p, matrix_vec_product(A[0], p, temp, n), n);
        vec_add(x, vec_multiple(p, temp, alpha, n), x, n);
        vec_add(r, vec_multiple(matrix_vec_product(A[0], p, temp, n), temp, -alpha, n), temp,
                n); // temp = r - alpha*A*P
        double beta = inner_product(temp, temp, n) / inner_product(r, r, n);
        memcpy(r, temp, sizeof(double) * n); // r <- temp
        vec_add(r, vec_multiple(p, p, beta, n), p, n);
        printf("iteration %d: ", i);
        print_vec(x, n);
    }
    return 0;
}

Go

main.go

package main

import (
    "fmt"
)

// a matrix with M rows and N cols.
type Matrix ([][]float64)

// a vector with N data
type Vec ([]float64)

// return nner product of two vector v and v1
func (v Vec) inner_product(v1 Vec) float64 {
    var r float64
    for i := 0; i < len(v); i++ {
        r += v[i] * v1[i]
    }
    return r
}

// add two vectors: v + eff*v1
func (v Vec) add(v1 Vec, eff float64) Vec {
    v2 := make(Vec, len(v))
    for i := 0; i < len(v); i++ {
        v2[i] = v[i] + eff*v1[i]
    }
    return v2
}

// return a vector of eff*m*vec
func (m Matrix) product(vec Vec, eff float64) Vec {
    c := make(Vec, len(vec))
    for i := 0; i < len(m); i++ {
        c[i] = 0
        for j := 0; j < len(m[i]); j++ {
            c[i] = c[i] + eff*m[i][j]*vec[j]
        }
    }
    return c
}

func main() {
    A := Matrix{
        {10, -1, 2, 0},
        {-1, 11, -1, 3},
        {2, -1, 10, -1},
        {0, 3, -1, 8}}
    b := Vec{6, 25, -11, 15}
    x := Vec{0.0, 0.0, 0.0, 0.0}

    // init
    r := b.add(A.product(x, -1.0), 1.0)
    var p = make(Vec, len(r))
    copy(p, r)
    // iteration
    for i := 0; i < len(b); i++ {
        alpha := r.inner_product(r) / p.inner_product(A.product(p, 1.0))
        x = x.add(p, alpha)
        r_temp := r.add(A.product(p, -alpha), 1.0)
        beta := r_temp.inner_product(r_temp) / r.inner_product(r)
        r = r_temp
        p = r.add(p, beta)
        fmt.Printf("iteration %d: ", i+1)
        for k := 0; k < len(x); k++ {
            fmt.Print(x[k], " ")
        }
        fmt.Println()
    }
}

将代码保存为main.c 或者 main.go, 然后编译运行：

C

$ clang main.c -o cg
$ ./cg
iteration 1: 0.471626 1.965108 -0.864648 1.179065 
iteration 2: 0.996432 1.976565 -0.909847 1.097591 
iteration 3: 1.001525 1.983269 -1.009858 1.019696 
iteration 4: 1.000000 2.000000 -1.000000 1.000000

Go

$ go run main.go 
iteration 1: 0.4716259464522676 1.9651081102177816 -0.8646475684958239 1.1790648661306689 
iteration 2: 0.996432359996456  1.9765653145545585 -0.9098469449042644 1.097591134432166 
iteration 3: 1.0015248100222705 1.9832687659087387 -1.009858497868728  1.019695902152845 
iteration 4: 1 2 -1.0000000000000002 1