共轭梯度法的简单直观理解
参考资料
本文是参考以下内容,结合自己的理解做的笔记。请尽量直接访问下述网页。
- 矩阵与数值计算(11)——共轭梯度法
- 共轭梯度法(一):线性共轭梯度
- 无痛版共轭梯度法介绍(更新到第五章)
- 共轭梯度法通俗讲义
第4个资料尤其清晰完备,很多都是参考它的。
What is: 什么是共轭梯度法?(简单直观理解)
共轭梯度法可以看作是梯度下降法(又称最速下降法)的一个改进。
对梯度下降来说
x
⃗
i
+
1
=
x
⃗
i
−
α
∇
f
\vec x_{i+1}=\vec x_i - \alpha\nabla f
xi+1=xi−α∇f
其中
α
\alpha
α控制了一步要走多远,因此被称为步长,在机器学习里面又称为学习率。
梯度下降法x移动的方向正是函数f的负梯度方向,这代表了局部上f减小最快的方向。
但是局部上减小最快的方向并不代表全局上指向最终解的方向。所以梯度下降法会出现像醉汉下山一样走出zig-zag的路线。如下图
图1 梯度下降法在2维解空间(也就是解向量只有两个维度)走出的路径示意图。
注:假如A是正定对称阵,其2维解空间必定是椭圆的。
为什么会走出这一Z形线呢?因为梯度下降的方向恰好与f垂直,也就是说和等高线垂直。沿着垂直于等高线的方向,一定能让函数减小,也就是最快地下了一个台阶。但是最快下台阶并不意味着最快到达目标位置(即最优解),因为你最终的目标并不是直对着台阶的。
为了修正这一路线,采用另一个方向:即共轭向量的方向。
我们先暂且给出共轭梯度法最后的形式,方便字母的定义:
x
⃗
i
+
1
=
x
⃗
i
−
α
d
⃗
i
\vec x_{i+1}=\vec x_i - \alpha \vec d_i
xi+1=xi−αdi
对照梯度下降法,每次向下走的方向不是梯度了,而是专门的一个方向
d
⃗
\vec d
d。除此之外和梯度下降法几乎一样。
在推进下一步之前,我们来看看什么是向量共轭。
共轭向量
下面先简要介绍共轭向量
所谓共轭向量,在数学上即:
p
i
T
A
p
j
=
0
p_i^TAp_j=0
piTApj=0
其中A是一个对称正定矩阵。
p
i
p_i
pi和
p
j
p_j
pj是一对共轭的向量。
可见,共轭是正交的推广化,因为向量正交的定义为:
p
i
T
p
j
=
0
p_i^Tp_j=0
piTpj=0
共轭比正交中间只多了个矩阵A,而矩阵的几何意义正是对一个向量进行线性变换(可见Gilber Strang的线代公开课)。因此共轭向量的意思就是一个向量经过线性变换(缩放剪切和旋转)之后与另一个向量正交。
共轭方向
言归正传,如何找到一个更好的方向呢?我们首先可以看看最完美的方向是什么样的。
下面这张图展示的就是最完美的方向。图中向量e代表的是误差。向量d就是方向向量。
误差e即当前迭代所得到的解与精确解的差值:
e
⃗
i
=
x
⃗
i
−
x
⃗
∗
\vec e_i=\vec x_i- \vec x^*
ei=xi−x∗
可惜我们并不能找到误差向量e,因为我们不知道精确解。
那么退而求其次,我们就找误差向量的共轭向量。因为图中可以看出,误差向量是与方向向量垂直的,即正交。刚才说了,共轭就是正交的推广。一个向量乘以一个矩阵之后与另一个方向正交,就是共轭。
即找到
d
⃗
T
A
e
⃗
=
0
\vec d ^T A \vec e =0
dTAe=0
但是这个公式里面仍然含有e,我们必须想办法去掉它,换成一个我们可以计算的量。
在推进下一步之前,我们先来看看误差与残差这两个概念的区别。
误差与残差
前面写道:
误差error 即当前迭代所得到的解与精确解的差值:
e
⃗
i
=
x
⃗
i
−
x
⃗
∗
\vec e_i=\vec x_i- \vec x^*
ei=xi−x∗
但是这种定义显然是没法直接用的,因为我们不知道精确解 x ∗ x^* x∗
那么退而求其次,我们想到,当误差收敛为0的时候,必然满足方程Ax=b,那么由此就可以定义出残差residual:
r ⃗ i = b ⃗ − A x ⃗ i \vec r_i=\vec b-A\vec x_i ri=b−Axi
这个定义的精妙之处在于,它定义了Ax接近b的距离,当距离为0的时候,恰好就是精确解。但是又能避开精确解本身。
在实际的程序中,我们还常常定义相对残差,即上一步迭代和这一步迭代的残差的相对变化率,这里就不再赘述。
显然,误差和残差之间就差了一个矩阵A,他们两者的关系是这样的:
r ⃗ i = b ⃗ − A ( e ⃗ i + x ⃗ ∗ ) = b ⃗ − A x ⃗ ∗ − A e ⃗ i = − A e ⃗ i \vec r_i=\vec b - A(\vec e_i+\vec x^*)=\vec b - A \vec x^* -A\vec e_i = -A\vec e_i ri=b−A(ei+x∗)=b−Ax∗−Aei=−Aei
可见除了A,还多了个负号。
搜索方向的确定
言归正传,利用残差,我们终于可以把误差e给替换掉了:
于是前面的式子就变成了
d
⃗
i
T
A
e
⃗
i
=
−
d
⃗
i
T
r
⃗
i
=
0
\vec d_i ^T A \vec e_i =-\vec d_i ^T \vec r_i=0
diTAei=−diTri=0
那么,这告诉我们:方向向量d,正是与残差向量正交的方向!
接下来我们只需要构建一个与残差正交的向量就可以了。这部分内容是由施密特正交化(更严谨一点的说法是共轭格莱姆-施密特过程)完成的。由于只是一个计算的方法,对概念的理解没有帮助,所以我们跳过,直接给出结论。
每一步搜索方向的时候,这一方向与残差以及前一步的方向有关
d
⃗
i
+
1
=
r
⃗
i
+
1
+
β
i
+
1
d
⃗
i
\vec d_{i+1} = \vec r_{i+1} +\beta_{i+1} \vec d_i
di+1=ri+1+βi+1di
其中系数
β
\beta
β可以这样计算:
β
i
+
1
=
r
⃗
i
+
1
T
r
⃗
i
+
1
r
⃗
i
T
r
⃗
i
\beta_{i+1} = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec r_{i}^T \vec r_{i} }
βi+1=riTriri+1Tri+1
这个系数beta其实很好记,因为分子就是残差的内积(下一步),分母也是残差的内积(这一步)。
或者说分子就是残差长度的平方(下一步),分母也是残差长度的平方(这一步)。(向量自己和自己的内积就是它的长度)
从另一个角度额外补充一点理解:
每次走的方向恰好是与残差正交的,这意味着:
每走一步恰好能消除残差的一个方向!
所以,当消除了残差所有投影方向上的值,那么就消除了整个残差!
步长,或者说系数alpha
实际上,还有一点没有解决,就是系数 α \alpha α怎么算?
这点的解释我们以后有机会再说,直接给出结论。
α
i
=
r
⃗
i
+
1
T
r
⃗
i
+
1
d
⃗
i
T
A
d
⃗
i
\alpha_i = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec d_{i}^T A\vec d_{i} }
αi=diTAdiri+1Tri+1
这个alpha的分子和beta一样,就是残差的内积。分母则是方向向量在乘以矩阵A之后的内积。
How to: 怎么用共轭梯度法?(完整算法)
-
设定初值
d ⃗ 0 = r ⃗ 0 = b ⃗ − A x ⃗ 0 \vec d_0=\vec r_0 = \vec b - A \vec x_0 \\ d0=r0=b−Ax0 -
计算系数alpha
α i = r ⃗ i + 1 T r ⃗ i + 1 d ⃗ i T A d ⃗ i \alpha_i = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec d_{i}^T A\vec d_{i} } αi=diTAdiri+1Tri+1 -
迭代一步(向下走一步)
x ⃗ i + 1 = x ⃗ i − α i d ⃗ i \vec x_{i+1}=\vec x_i - \alpha_i \vec d_i xi+1=xi−αidi -
计算残差(此处已经被修改,原文没有被50整除那一个公式 2022-05-27)
如果迭代次数可以被50整除
r ⃗ i + 1 = r ⃗ i − α i A x ⃗ \vec r_{i+1}=\vec r_i - \alpha_i A\vec x ri+1=ri−αiAx
否则
r ⃗ i + 1 = r ⃗ i − α i A d \vec r_{i+1}=\vec r_i - \alpha_i A d ri+1=ri−αiAd -
计算系数beta
β i + 1 = r ⃗ i + 1 T r ⃗ i + 1 r ⃗ i T r ⃗ i \beta_{i+1} = \frac{ \vec r_{i+1}^T \vec r_{i+1} } {\vec r_{i}^T \vec r_{i} } βi+1=riTriri+1Tri+1 -
计算搜索方向 d ⃗ \vec d d
d ⃗ i + 1 = r ⃗ i + 1 + β i + 1 d ⃗ i \vec d_{i+1} = \vec r_{i+1} +\beta_{i+1} \vec d_i di+1=ri+1+βi+1di
重复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了呢?
实际上使用共轭梯度法的两个条件
- A是对称的
- A是正定的
也和这个原理有关。
数学家求解问题的思路是:把不会的问题转化成会的问题,再套用会的问题的思路求解问题。
为了说明这一点,我们要从线性代数的二次型入手。我们可以先复习一下二次型,了解一下它是什么。
二次型
二次型就是关于向量的二次函数。
我们高中学过的二次函数通用表达式为
f
(
x
)
=
a
x
2
+
b
x
+
c
f(x) = a x^2 +bx+c
f(x)=ax2+bx+c
如果把其中的x替换为向量x,并且把a x^2 替换为
x^T A x 就得到了
f ( x ) = x T A x + b x + c f(x) = x^T A x +bx+c f(x)=xTAx+bx+c
这就是二次型。
二次型求导得到
f
′
(
x
)
=
1
2
(
A
x
+
A
T
x
)
+
b
f'(x) = \frac{1}{2}( A x + A^T x)+b
f′(x)=21(Ax+ATx)+b
将Ax=b问题转化为最优化问题
我们本来求解的是
A
x
=
b
A\mathbf x=\bf b
Ax=b
这个问题被转化为了求某个函数的导数等于0的问题,即驻值问题。
我们设这个函数为g(x)。我们的问题即:
g
′
(
x
)
=
0
x
∗
=
a
r
g
m
i
n
x
g
(
x
)
g'(x)=0\\ x^*=argmin_x g(x)
g′(x)=0x∗=argminxg(x)
argmin_x的意思就是我们求取最小值的时候的x,而不是最小值本身。
这个 x ∗ x^* x∗就是最终解。
那么怎么联系到Ax=b呢?
我们只要改造这个函数g,让它的导数恰好就是Ax-b=0就好了!!
而这个函数,恰好就是二次型函数!
即
g
′
(
x
)
=
A
x
−
b
g'(x)=Ax-b
g′(x)=Ax−b
于是求最小值得问题就能够被转化为求Ax=b的问题!
这里有个小小的瑕疵:
实际上,二次型g(x)的导数是
g
′
(
x
)
=
1
/
2
(
A
T
+
A
)
x
−
b
g'(x)=1/2 (A^T+A)x-b
g′(x)=1/2(AT+A)x−b
所以我们就要限定 A T = A A^T=A AT=A,即限定A是对称的。这就是第一个条件的由来!
to be continued
2022-05-20