共轭梯度法的搜索方向加上负梯度方向是为了什么？

显示全部楼层 · 2022-12-2 13:20:04

模拟不错 · 2022-12-2 14:22:03

在数值代数的语境中，对于对称正定矩阵A，要求解

就可以考虑共轭梯度法。共轭梯度法（与其切比雪夫加速变种）是求解对称正定线性系统的标准算法，但共轭梯度法也经常出现在非线性优化的语境里，这是由于上面线性系统的求解与下面的极小化问题等价：

从优化的角度看有诸多便利，比如可以求梯度：

用梯度下降可以求解线性方程，这个效率不太高，于是有共轭梯度法来救场。
若干假设

接下来将在以下几个假设下推导出共轭梯度法：

首先CG作为一种梯度法，其核心要素在于下降方向与步长：
. d/ G4 Z, P* k( W8 } 其中是方向，是步长。1 ]- E+ V0 o: X2 F+ C
构成矩阵A导出内积空间下的正交基,4 y. a5 f7 j3 V+ W7 W$ c
由此，可以写成5 S! [ a8 [) o# C8 }! @# \( W* v
5 L8 o1 n3 Q, s" {$ Y
对于迭代点, 还希望其具有某种最优性，比如：
{) f M( p% I3 e" l2 ` 其中,随着空间的扩张，会逐渐逼近.! j$ R8 J) @9 m/ B8 ^$ @+ b
这个最优性是共轭梯度法的本质，如果变更最优准则（比如二范数），还可以导出其他算法，不提。
: Q+ Q( z- i, i
也不能随意取，需要充分利用已知信息，取$ o: @3 z+ N$ k* d
其中, 也挺合理，即为之前点所在位置的梯度所构成的子空间。9 e' F/ g$ d8 Y4 C5 M y5 u
另外，给定初始值后，初始方向也有了，
9 E( P X% {1 d- x9 _' z

共轭梯度法的性质

有了上面这些对于算法的假设或者限制，就可以得到更多关于算法的性质了：

步长的实质与的表示如。果没有假设3的话，是可以随意取值的，然而有了性质3，就有：. T3 G- I8 {+ S7 p

由于正交性，立马得到

.这说明了

，即有限步一定收敛，在收敛性分析中，虽然有限步收敛听起来非常强，但通常N非常大，最后分析线性收敛的界更具有参考性。

由假设4，由于, 那么有6 Z% N1 R" o/ C
除了上面这条，还有一条非常类似的,由+ F( k0 f8 f* F) W$ ]6 f% f
推出：. Z2 N! O, o) I! j, q( \% X1 W
当时，由正交性，就有% f& ]$ j/ y; ?+ g5 `* W1 ~+ W
即与.! I2 z. \0 I P' F* R
共轭梯度法已经呼之欲出了，利用正交化，先求出方向：
7 S6 M2 W h* {! B1 c$ P7 Q/ p. c
5 q* U& d2 t- Q
再考虑下具体步长，由于- d$ v8 u5 e" U7 k, p
那么有
5 g7 D! h. S. t6 y. ^/ o 也即
0 J* B' G* Y5 d" I( F
: V: K4 M6 n6 A4 X8 h " M `& V4 }8 B6 g# W: A2 f( I% S! [0 M8 }
此时，就得到共轭梯度法的第一个形式
* F. W% z" Z0 L$ z. M: ]+ ?2 p
! x( d0 ]: I9 {9 { 当然，这个形式并不太像教科书上的共轭梯度法，但理论上已经一致了，接下来需要化简。/ ~- ^/ }$ L( K$ C+ T* {8 V
由于, 那么有,性质 4中的式子可以大幅简化，结合步长公式，变成：8 f8 k" k, F/ j* C
; w0 Y/ N; k: Z" \4 M0 g7 s
上面的推导需要计算内积，这个还可以被优化掉，注意到
! K) X5 t0 u# T- \" k- x9 Q 那么有.' i3 h% l* h9 A c: u' C
最后的共轭梯度法：/ ^7 H' }4 I% u% ^7 S! }4 C
) o5 X* E5 W+ k5 d# G0 Y1 i

算例

CG的实现非常简单，考虑求解一个Poisson边值问题，上次用Deepritz方法求解过：

离散Poisson方程也讲过很多了，这里不讲如何进行方程的理算与构建了，可参考之前的文章
派大西：泊松方程求解与卫星组件热布局快速预测.py网格数量为

, 那么对应的线性系统矩阵

，初始化

：
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

Nx = 201
Ny = 201
dx = 2 / (Nx - 1)
dy = 2 / (Ny - 1)
A = sp.sparse.lil_matrix((Nx * Ny, Nx * Ny))
b = np.zeros(Nx * Ny)
q = lambda x, y: 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)
for i in range(Ny):
for j in range(Nx):
      xp = j * 2 / (Nx - 1) - 1
      yp = i * 2 / (Ny - 1) - 1
      if i == 0 or i == Ny - 1 or j == 0 or j == Ny - 1:
         A[i * Nx + j, i * Nx + j] = 1
         b[i * Nx + j] = 0.
         continue
      A[i * Nx + j, i * Nx + j] = 1 / dx + 1 / dy
      A[i * Nx + j, i * Nx + j + 1] = -1 / 2 / dx
      A[i * Nx + j, i * Nx + j - 1] = -1 / 2 / dx
      A[i * Nx + j, (i - 1) * Nx + j] = -1 / 2 / dy
      A[i * Nx + j, (i + 1) * Nx + j] = -1 / 2 / dy
      b[i * Nx + j] = q(xp, yp)
定义共轭梯度法：
def cg(A, b, x0=None, tor=1e-10):
N = len(b)
if x0 is None:
      x = np.random.randn(N)
else:
      x = x0.copy()
res = np.Inf
r = b - A @ x
d = r
r_inner_old = np.sum(r * r)
iteration = 0
while res > tor:
      Ad = A @ d
      alpha = r_inner_old / np.sum(Ad * d)
      x = x + alpha * d
      if iteration % 50 != 0:
         r = r - alpha * Ad
      else:
         r = b - A @ x
      r_inner_new = np.sum(r * r)
      beta = r_inner_new / r_inner_old
      d = r + beta * d
      r_inner_old = r_inner_new
      res = np.linalg.norm(r_inner_new)
      print(iteration, res)
      iteration += 1
return x

最后求解得到：

共轭梯度法的搜索方向加上负梯度方向是为了什么？-51.jpg

看起来没啥问题，然后是时间与规模的关系：

共轭梯度法的搜索方向加上负梯度方向是为了什么？-52.jpg

时间的增长可能来源于矩阵向量乘法的耗时增加与迭代步数的增加，看起来会比迭代步数增长更快，如果追求更高的规模或精度，就需要引入一些预条件了，下次再说吧。

苏麒麟麟v · 2022-12-2 22:42:37

被邀请了，那就班门弄斧，抛砖引玉一下。（术语我知道中文是啥的，我会用中文，不知道怎么翻译的我就保留英文了）
Reference: Conjugate Gradient Descent Notes (我基本就是照着这个来的，总结和翻译了一下）
一句话总结：为什么要加上负梯度方向？因为CG的原理导致下降方向正好和负梯度方向有关，不是为了什么原因特地加入负梯度。
（详细解释）
共轭梯度法本来是为了求解线性系统

，这等价于最小化

。当然这个方法也被拓展到用于优化一些别的东西。不过以下都假设我们要最小化

。
下面介绍几个必要的概念：
1） Krylov subspaces：
一系列nested subspaces:

.显然

Fact 1:

Proof: From Cayley-Hamiltion theorem. (Refer to the reference page 3)

可以看出来线性系统的最优解就在

里，那么如何找到它？
2）Krylov sequences:
定义

，我们将

称为Krylov sequences。显然

。这说明了如果我们只要有这么

个Krylov sequences，我们就能找到

的最优解。本质上共轭梯度法就是生成这些Krylov sequences，当你生成到

时就结束了。

Krylov sequences 好棒啊，那么怎么找到它们？又要介绍一通概念了。
3）Optimal condition on

:
有两个要求:

4) Residual

我们接下来定义residual

, 可以观察到

。同时我们也看出来

。这里出现负梯度啦！！！
根据residual的性质，我们得到

and

。
5）Conjugate direction

定义：

,有如下性质：

和

. 证明在reference第七页。
我们就得到Krylov的一个&#39;conjugate&#39; basis:

可以看出来，Krylov sequence的之间变化的方向和负梯度有关。这就是为什么看CG算法里会有加入了负梯度的感觉。但是其实这个思路是反过来的，不是为了什么特意加入负梯度，而是在Krylov sequence求解的时候自然而然的出现了负梯度。

共轭梯度法的搜索方向加上负梯度方向是为了什么？

浏览过的版块