共轭梯度法的简单直观理解

参考资料

本文是参考以下内容，结合自己的理解做的笔记。请尽量直接访问下述网页。

1. 矩阵与数值计算（11）——共轭梯度法
2. 共轭梯度法（一）：线性共轭梯度
3. 无痛版共轭梯度法介绍(更新到第五章)
4. 共轭梯度法通俗讲义
   第4个资料尤其清晰完备，很多都是参考它的。

What is: 什么是共轭梯度法？(简单直观理解）

共轭梯度法可以看作是梯度下降法（又称最速下降法）的一个改进。

对梯度下降来说
$${\vec{x}}_{i+1}={\vec{x}}_i-\alpha \nabla f$$
其中$\alpha$控制了一步要走多远，因此被称为步长，在机器学习里面又称为学习率。

梯度下降法x移动的方向正是函数f的负梯度方向，这代表了局部上f减小最快的方向。

但是局部上减小最快的方向并不代表全局上指向最终解的方向。所以梯度下降法会出现像醉汉下山一样走出zig-zag的路线。如下图

图1 梯度下降法在2维解空间（也就是解向量只有两个维度）走出的路径示意图。
注：假如A是正定对称阵，其2维解空间必定是椭圆的。

为什么会走出这一Z形线呢？因为梯度下降的方向恰好与f垂直，也就是说和等高线垂直。沿着垂直于等高线的方向，一定能让函数减小，也就是最快地下了一个台阶。但是最快下台阶并不意味着最快到达目标位置（即最优解），因为你最终的目标并不是直对着台阶的。

为了修正这一路线，采用另一个方向：即共轭向量的方向。

我们先暂且给出共轭梯度法最后的形式，方便字母的定义：
$${\vec{x}}_{i+1}={\vec{x}}_i-\alpha {\vec{d}}_i$$
对照梯度下降法，每次向下走的方向不是梯度了，而是专门的一个方向$\vec{d}$。除此之外和梯度下降法几乎一样。

在推进下一步之前，我们来看看什么是向量共轭。

共轭向量

下面先简要介绍共轭向量

所谓共轭向量，在数学上即：
$$p_i^{T}A{p}_j=0$$

其中A是一个对称正定矩阵。
$p_i$和$p_j$是一对共轭的向量。

可见，共轭是正交的推广化，因为向量正交的定义为：
$$p_i^T{p}_j=0$$
共轭比正交中间只多了个矩阵A，而矩阵的几何意义正是对一个向量进行线性变换（可见Gilber Strang的线代公开课）。因此共轭向量的意思就是一个向量经过线性变换（缩放剪切和旋转）之后与另一个向量正交。

共轭方向

言归正传，如何找到一个更好的方向呢？我们首先可以看看最完美的方向是什么样的。

下面这张图展示的就是最完美的方向。图中向量e代表的是误差。向量d就是方向向量。

误差e即当前迭代所得到的解与精确解的差值：
$${\vec{e}}_i={\vec{x}}_i-{\vec{x}}^{\ast }$$

可惜我们并不能找到误差向量e，因为我们不知道精确解。

那么退而求其次，我们就找误差向量的共轭向量。因为图中可以看出，误差向量是与方向向量垂直的，即正交。刚才说了，共轭就是正交的推广。一个向量乘以一个矩阵之后与另一个方向正交，就是共轭。

即找到
$${\vec{d}}^{T}A\vec{e}=0$$

但是这个公式里面仍然含有e，我们必须想办法去掉它，换成一个我们可以计算的量。

在推进下一步之前，我们先来看看误差与残差这两个概念的区别。

误差与残差

前面写道：

误差error 即当前迭代所得到的解与精确解的差值：
$${\vec{e}}_i={\vec{x}}_i-{\vec{x}}^{\ast }$$

但是这种定义显然是没法直接用的，因为我们不知道精确解$x^{\ast }$

那么退而求其次，我们想到，当误差收敛为0的时候，必然满足方程Ax=b，那么由此就可以定义出残差residual：

$${\vec{r}}_i=\vec{b}-A{\vec{x}}_i$$

这个定义的精妙之处在于，它定义了Ax接近b的距离，当距离为0的时候，恰好就是精确解。但是又能避开精确解本身。

在实际的程序中，我们还常常定义相对残差，即上一步迭代和这一步迭代的残差的相对变化率，这里就不再赘述。

显然，误差和残差之间就差了一个矩阵A，他们两者的关系是这样的：

$${\vec{r}}_i=\vec{b}-A({\vec{e}}_i+{\vec{x}}^{\ast })=\vec{b}-A{\vec{x}}^{\ast }-A{\vec{e}}_i=-A{\vec{e}}_i$$

可见除了A，还多了个负号。

搜索方向的确定

言归正传，利用残差，我们终于可以把误差e给替换掉了：
于是前面的式子就变成了
$${\vec{d}}_i^{T}A{\vec{e}}_i=-{\vec{d}}_i^T{\vec{r}}_i=0$$

那么，这告诉我们：方向向量d，正是与残差向量正交的方向！

接下来我们只需要构建一个与残差正交的向量就可以了。这部分内容是由施密特正交化（更严谨一点的说法是共轭格莱姆-施密特过程）完成的。由于只是一个计算的方法，对概念的理解没有帮助，所以我们跳过，直接给出结论。

每一步搜索方向的时候，这一方向与残差以及前一步的方向有关
$${\vec{d}}_{i+1}={\vec{r}}_{i+1}+{\beta }_{i+1}{\vec{d}}_i$$
其中系数$\beta$可以这样计算：
$${\beta }_{i+1}=\frac{{\vec{r}}_{i+1}^T{\vec{r}}_{i+1}}{{\vec{r}}_i^T{\vec{r}}_i}$$

这个系数beta其实很好记，因为分子就是残差的内积（下一步），分母也是残差的内积（这一步）。
或者说分子就是残差长度的平方（下一步），分母也是残差长度的平方（这一步）。（向量自己和自己的内积就是它的长度）

从另一个角度额外补充一点理解：
每次走的方向恰好是与残差正交的，这意味着：
每走一步恰好能消除残差的一个方向！
所以，当消除了残差所有投影方向上的值，那么就消除了整个残差！

步长，或者说系数alpha

实际上，还有一点没有解决，就是系数$\alpha$怎么算？

这点的解释我们以后有机会再说，直接给出结论。
$${\alpha }_i=\frac{{\vec{r}}_{i+1}^T{\vec{r}}_{i+1}}{{\vec{d}}_i^{T}A{\vec{d}}_i}$$
这个alpha的分子和beta一样，就是残差的内积。分母则是方向向量在乘以矩阵A之后的内积。

How to: 怎么用共轭梯度法？(完整算法）

1. 设定初值
   $${\vec{d}}_0={\vec{r}}_0=\vec{b}-A{\vec{x}}_0$$

2. 计算系数alpha
   $${\alpha }_i=\frac{{\vec{r}}_{i+1}^T{\vec{r}}_{i+1}}{{\vec{d}}_i^{T}A{\vec{d}}_i}$$

3. 迭代一步（向下走一步）
   $${\vec{x}}_{i+1}={\vec{x}}_i-{\alpha }_i{\vec{d}}_i$$

4. 计算残差（此处已经被修改，原文没有被50整除那一个公式 2022-05-27）
   如果迭代次数可以被50整除
   $${\vec{r}}_{i+1}={\vec{r}}_i-{\alpha }_{i}A\vec{x}$$
   否则
   $${\vec{r}}_{i+1}={\vec{r}}_i-{\alpha }_{i}Ad$$

5. 计算系数beta
   $${\beta }_{i+1}=\frac{{\vec{r}}_{i+1}^T{\vec{r}}_{i+1}}{{\vec{r}}_i^T{\vec{r}}_i}$$

6. 计算搜索方向$\vec{d}$
   $${\vec{d}}_{i+1}={\vec{r}}_{i+1}+{\beta }_{i+1}{\vec{d}}_i$$

重复2~6，直到残差足够小

Python代码

更正2024-4-29: 按照wiki重写的代码已经可以正常运行

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import cg
import scipy.sparse as sp
from time import time

# judge if A is positive definite
# https://stackoverflow.com/a/44287862/19253199
# if A is symmetric and able to be Cholesky decomposed, then A is positive definite
def is_pos_def(A):
    A=A.toarray()
    if np.array_equal(A, A.T):
        try:
            np.linalg.cholesky(A)
            print("A is positive definite")
            return True
        except np.linalg.LinAlgError:
            print("A is not positive definite")
            return False
    else:
        print("A is not positive definite")
        return False

def generate_A_b_psd(n=1000):
    A = sp.random(n, n, density=0.01, format="csr")
    A = A.T @ A
    b = np.random.rand(A.shape[0])
    print(f"Generated PSD A: {A.shape}, b: {b.shape}")
    A = sp.csr_matrix(A)
    assert is_pos_def(A)
    return A, b

def main():
    A,b = generate_A_b_psd()

    t=time()
    x_sp, exit_code = cg(A, b, atol=1e-5)
    print("scipy_cg time:", time()-t)
    t=time()
    x_my = my_cg(A, b)
    print("my_cg time:", time()-t)
    print("error:", np.linalg.norm(x_sp-x_my))


def my_cg(A, b, atol=1e-5):
    x=np.zeros(A.shape[0])
    r0=b-A@x
    p=r0.copy()
    r=r0.copy()
    k=0
    while True:
        Ap = A@p
        rTr = r.T@r
        alpha = rTr / (p.T@Ap)
        x1 = x + alpha * p
        r1 = r - alpha * Ap
        if np.linalg.norm(r1)<atol:
            break
        beta=r1.T@r1/(rTr)
        p1=r1+beta*p
        x=x1.copy()
        r=r1.copy()
        p=p1.copy()
        k+=1
    return x1

if __name__ == "__main__":
    main()

输出

Generated PSD A: (1000, 1000), b: (1000,)
A is positive definite
scipy_cg time: 0.26083922386169434
my_cg time: 0.2570023536682129
error: 0.0002792725223332572

来自wiki

https://en.wikipedia.org/wiki/Conjugate_gradient_method

PCG

def test_pcg():
    A,b = generate_A_b_psd()

    P = sp.diags(1/A.diagonal())

    t=time()
    x_sp, exit_code = cg(A, b, atol=1e-5, M=P)
    print("scipy_cg time:", time()-t)
    t=time()
    x_my = my_pcg(A, b, atol=1e-5, M=P)
    print("my_pcg time:", time()-t)
    print("error:", np.linalg.norm(x_sp-x_my, ord=np.inf))
    print("x(first 5):\n", x_sp[:5],"\n", x_my[:5])


# preconditioned conjugate gradient
# https://en.wikipedia.org/wiki/Conjugate_gradient_method#The_preconditioned_conjugate_gradient_method
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.cg.html
# Note: Based on the scipy(https://github.com/scipy/scipy/blob/7dcd8c59933524986923cde8e9126f5fc2e6b30b/scipy/sparse/linalg/_isolve/iterative.py#L406), parameter M is actually the inverse of M in the wiki's formula. We adopt the scipy's definition.
def my_pcg(A, b, atol=1e-5, M=None):
    def solvez(r):
        z = M@r if M is not None else r
        return z
    x=np.zeros(A.shape[0])
    r=b-A@x
    r=r.copy()
    z = solvez(r)
    p=z.copy()
    k=0
    while True:
        Ap = A@p
        rTz = r.T@z
        alpha = r.T@z / (p.T@Ap)
        x1 = x + alpha * p
        r1 = r - alpha * Ap
        if np.linalg.norm(r1)<atol:
            break
        z1 = solvez(r1)
        beta=r1.T@z1/(rTz)
        p1=z1+beta*p
        x=x1.copy()
        r=r1.copy()
        p=p1.copy()
        z=z1.copy()
        k+=1
    return x1


if __name__ == "__main__":
    test_pcg()

输出

Generated PSD A: (1000, 1000), b: (1000,)
A is positive definite
scipy_cg time: 0.24599957466125488
my_pcg time: 0.24100065231323242
error: 4.3591826397459954e-05
x(first 5):
 [31666.44122946   618.03821774  1318.11402122 -3403.45443571
 12217.37609203] 
 [31666.44123179   618.0382179   1318.11403268 -3403.45444025
 12217.37608712]

Why: 为什么共轭梯度法能求解Ax=b?

说了这么多，其实有一个关键问题没有讲，那就是：为什么共轭梯度法能求解Ax=b?

按理说，共轭梯度法是函数最优化的方法，怎么就扯上了求解Ax=b了呢？

实际上使用共轭梯度法的两个条件

1. A是对称的
2. A是正定的

也和这个原理有关。

数学家求解问题的思路是：把不会的问题转化成会的问题，再套用会的问题的思路求解问题。

为了说明这一点，我们要从线性代数的二次型入手。我们可以先复习一下二次型，了解一下它是什么。

二次型

二次型就是关于向量的二次函数。

我们高中学过的二次函数通用表达式为
$$f(x)=a{x}^2+bx+c$$

如果把其中的x替换为向量x，并且把a x^2 替换为
x^T A x 就得到了

$$f(x)=x^{T}Ax+bx+c$$

这就是二次型。

二次型求导得到
$$f^{\prime }(x)=\frac{1}{2}(Ax+A^{T}x)+b$$

将Ax=b问题转化为最优化问题

我们本来求解的是
$$A\mathbf{x}=\mathbf{b}$$

这个问题被转化为了求某个函数的导数等于0的问题，即驻值问题。

我们设这个函数为g(x)。我们的问题即：

$$\begin{gathered} g^{\prime }(x)=0 \\ {x}^{\ast }=argmi{n}_{x}g(x) \end{gathered}$$
argmin_x的意思就是我们求取最小值的时候的x，而不是最小值本身。

这个$x^{\ast }$就是最终解。

那么怎么联系到Ax=b呢？

我们只要改造这个函数g，让它的导数恰好就是Ax-b=0就好了！！
而这个函数，恰好就是二次型函数！

即
$$g^{\prime }(x)=Ax-b$$

于是求最小值得问题就能够被转化为求Ax=b的问题！

这里有个小小的瑕疵：
实际上，二次型g(x)的导数是
$$g^{\prime }(x)=1/2(A^T+A)x-b$$

所以我们就要限定$A^T=A$，即限定A是对称的。这就是第一个条件的由来！

to be continued
2022-05-20

拓展：改进——预处理的共轭梯度法