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

大家好，欢迎来到IT知识分享网。

参考资料

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

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

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

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

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

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

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

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

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

共轭向量

下面先简要介绍共轭向量

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

共轭方向

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

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

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

但是这个公式里面仍然含有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，还多了个负号。

搜索方向的确定

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

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

步长，或者说系数alpha

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

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// # 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. my_cg time: 0.82129 error: 0.00032572 

来自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/7dcd8c923cde8e9126f5fc2e6b30b/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. my_pcg time: 0. error: 4.59954e-05 x(first 5): [31666. 618.0 1318. -3403. 12217.] [31666. 618.0 1318. -3403. 12217.] 

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

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

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

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

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

也和这个原理有关。

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

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

二次型

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

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

这就是二次型。

将Ax=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就好了！！
而这个函数，恰好就是二次型函数！

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

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

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