一、共轭梯度算法简介

共轭梯度(Conjugate Gradient)方法是一种迭代算法，可用于求解无约束的优化问题，例如能量最小化。常见的优化算法还有梯度下降法，相比于梯度下降法，共轭梯度法具有收敛更快，以及算法稳定性更好的优点。

从上图可以看出来，梯度下降法优化过程中函数是沿着梯度的反方向逐步优化，后一步优化的结果会对前一步的优化结果造成影响，收敛较慢。而共轭梯度方法每一步的优化方向与前一步的优化方向是共轭的，因此并不会对前一步的优化结果造成影响，同时优化过程中保证在每一个方向上函数优化到最小值，从而保证沿着这些共轭方向优化完成后，函数能够达到全局最小值点。具体解释可以参考这篇文章：https://zhuanlan.zhihu.com/p/338838078

二、实现共轭梯度方法的两块重要积木

1.共轭方向的确定

共轭梯度方法中的新的共轭方向仅由上一步的梯度，新的梯度和上一步的优化方向决定，即：
d⃗k+1=−g⃗k+1+βkd⃗k\vec{d}_{k+1}=-\vec{g}_{k+1}+\beta_k\vec{d}_k dk+1=−gk+1+βkdk
其中βk\beta_kβk的定义方式有很多种，常用的有FR和PRP两种：
βkFR=gk+1Tgk+1gkTgk=∣g⃗k+1∣2∣g⃗k∣2\beta_k^{FR}=\frac{g_{k+1}^Tg_{k+1}}{g_k^Tg_k}=\frac{|\vec{g}_{k+1}|^2}{|\vec{g}_k|^2} βkFR=gkTgkgk+1Tgk+1=∣gk∣2∣gk+1∣2
βkPRP=gk+1T(gk+1−gk)gkTgk=g⃗k+1⋅(g⃗k+1−g⃗k)∣g⃗k∣2\beta_k^{PRP}=\frac{g_{k+1}^T(g_{k+1}-g_k)}{g_k^Tg_k}=\frac{\vec{g}_{k+1}\cdot(\vec{g}_{k+1}-\vec{g}_k)}{|\vec{g}_k|^2} βkPRP=gkTgkgk+1T(gk+1−gk)=∣gk∣2gk+1⋅(gk+1−gk)

2.方向优化步长的确定

当优化方向确定之后，需要利用线性搜索技术确定优化的步长αmin\alpha_{min}αmin，即是寻找α>0使得
f(x⃗+αmind⃗)=minf(x⃗+αd⃗)f(\vec{x}+\alpha_{min}\vec{d})=minf(\vec{x}+\alpha\vec{d}) f(x+αmind)=minf(x+αd)
线搜索技术包含两大类，即精确线搜索技术和非精确线搜索技术，详见这篇文章：https://www.longzf.com/optimization/2/line_search/
精确线搜索技术包括牛顿法和二分法。由于之后的代码使用的线搜索技术是牛顿法，所以下面简单介绍一下牛顿法。
牛顿法被用于求解函数的极小极大值问题，函数在极值点处应有f′(α)=0f'(\alpha)=0f′(α)=0，将函数展开到二阶泰勒展开之后，可以得到α的迭代公式：
αk+1=αk−f′(αk)f′′(αk)\alpha_{k+1} = {\alpha}_k - \frac{f'(\alpha_k)}{f''(\alpha_k)} αk+1=αk−f′′(αk)f′(αk)

note

由于数值计算过程中要用到差分的方法，所以这里简单列出以下使用一阶微分和二阶微分所使用的差分表达式：
f′(x)≈f(x+δ)−f(x−δ)2δf'(x)\approx\frac{f(x+\delta )-f(x-\delta)}{2\delta} f′(x)≈2δf(x+δ)−f(x−δ)
f′′(x)≈f(x+δ)+f(x−δ)−2f(x)δ2f''(x)\approx\frac{f(x+\delta)+f(x-\delta)-2f(x)}{\delta^2} f′′(x)≈δ2f(x+δ)+f(x−δ)−2f(x)
其中δ\deltaδ是小量。

三、共轭梯度算法优化过程

计算初始初始梯度值g⃗0\vec{g}_0g0和优化方向d⃗0\vec{d}_0d0
g⃗0=∇⃗f(x⃗0)\vec{g}_0=\vec{\nabla}f(\vec{x}_0) g0=∇f(x0)
d⃗0=−g⃗0\vec{d}_0 = -\vec{g}_0 d0=−g0
如果g⃗k<ϵ⃗\vec{g}_k<\vec{\epsilon}gk<ϵ，退出迭代过程，否则执行以下步骤
用线性搜索算法(牛顿法)求出使得minf(x⃗k+αd⃗k)minf(\vec{x}_k+\alpha\vec{d}_k)minf(xk+αdk)的步长α，并更新x⃗k+1\vec{x}_{k+1}xk+1:
x⃗k+1=x⃗k+αd⃗k\vec{x}_{k+1}=\vec{x}_k+\alpha\vec{d}_k xk+1=xk+αdk
计算新的梯度
g⃗k+1=∇⃗f(x⃗k+1)\vec{g}_{k+1}=\vec{\nabla}f(\vec{x}_{k+1}) gk+1=∇f(xk+1)
根据前面提到的FR公式或PRP公式计算新的组合系数
计算新的共轭方向：
d⃗k+1=−g⃗k+1+βkd⃗k\vec{d}_{k+1}=-\vec{g}_{k+1}+\beta_k\vec{d}_k dk+1=−gk+1+βkdk
重复执行第2步

四、python实现共轭梯度算法

这里使用的测试函数形式是：
f(x,y)=(3x−2y)2+(x−1)4f(x, y)=(3x-2y)^2+(x-1)^4 f(x,y)=(3x−2y)2+(x−1)4
可以看到函数存在最小值点(1, 1.5)，以下是实现的python代码。

1.FR-CG

############共轭梯度算法#############
import numpy as np
import matplotlib.pyplot as plt
def testFun(x, y):t = 3.0*x - 2.0*yt1 = x - 1.0z = np.power(t, 2) + np.power(t1, 4)return zdef gradTestFun(x, y):'''求函数的梯度'''delta_x = 1e-6      #x方向差分小量delta_y = 1e-6      #y方向差分小量grad_x = (testFun(x+delta_x, y)-testFun(x-delta_x, y))/(2.0*delta_x)grad_y = (testFun(x, y+delta_y)-testFun(x, y-delta_y))/(2.0*delta_y)grad_xy = np.array([grad_x, grad_y])return grad_xydef getStepLengthByNewton(array_xy, array_d):'''采用牛顿法，精确线性搜索确定移动步长'''a0 = 1.0           #初始猜测值e0 = 1e-6          #退出搜索循环的条件delta_a = 1e-6     #对a作差分的小量while(1):new_a = array_xy + a0*array_dnew_a_l = array_xy + (a0-delta_a)*array_dnew_a_h = array_xy + (a0+delta_a)*array_ddiff_a0 = (testFun(new_a_h[0], new_a_h[1]) - testFun(new_a_l[0], new_a_l[1]))/(2.0*delta_a)if np.abs(diff_a0) < e0:breakddiff_a0 = (testFun(new_a_h[0], new_a_h[1]) + testFun(new_a_l[0], new_a_l[1]) - 2.0*testFun(new_a[0], new_a[1]))/(delta_a*delta_a)a0 = a0 - diff_a0/ddiff_a0return a0def plotResult(array_xy_history):x = np.linspace(-1.0, 4.0, 100)y = np.linspace(-4.0, 8.0, 100)X, Y = np.meshgrid(x, y)Z = testFun(X, Y)plt.figure(dpi=300)plt.xlim(-1.0, 4.0)plt.ylim(-4.0, 8.0)plt.xlabel("x")plt.ylabel("y")plt.contour(X, Y, Z, 40)plt.plot(array_xy_history[:,0], array_xy_history[:,1], marker='.', ms=10)xy_count = array_xy_history.shape[0]for i in range(xy_count):if i == xy_count-1:breakdx = (array_xy_history[i+1][0] - array_xy_history[i][0])*0.6dy = (array_xy_history[i+1][1] - array_xy_history[i][1])*0.6plt.arrow(array_xy_history[i][0], array_xy_history[i][1], dx, dy, width=0.1)def mainFRCG():'''使用CG算法优化，用FR公式计算组合系数'''ls_xy_history = []                           #存储初始坐标的迭代结果xy0 = np.array([4.0, -2.0])                   #初始点grad_xy = gradTestFun(xy0[0], xy0[1])d = -1.0*grad_xy                             #初始搜索方向e0 = 1e-6                                    #迭代退出条件xy = xy0while(1):ls_xy_history.append(xy)grad_xy = gradTestFun(xy[0], xy[1])tag_reach = np.abs(grad_xy) < e0if tag_reach.all():breakstep_length = getStepLengthByNewton(xy, d)xy_new = xy + step_length*dgrad_xy_new = gradTestFun(xy_new[0], xy_new[1])b = np.dot(grad_xy_new, grad_xy_new)/np.dot(grad_xy, grad_xy)       #根据FR公式计算组合系数d = b*d - grad_xy_newxy = xy_newarray_xy_history = np.array(ls_xy_history)plotResult(array_xy_history)return array_xy_history

以下是运行结果

最终得到的最小值的坐标为[ 1.00113664, 1.50170496]
运行时间：

2.PRP-CG

换用PRP公式计算只需要将上面求解组合系数部分的代码作少量修改即可，代码为：

def mainPRPCG():'''使用CG算法优化，用PRP公式计算组合系数'''ls_xy_history = []xy0 = np.array([4.0, -2.0])grad_xy = gradTestFun(xy0[0], xy0[1])d = -1.0*grad_xye0 = 1e-6xy = xy0while(1):ls_xy_history.append(xy)grad_xy = gradTestFun(xy[0], xy[1])tag_reach = np.abs(grad_xy) < e0if tag_reach.all():breakstep_length = getStepLengthByNewton(xy, d)xy_new = xy + step_length*dgrad_xy_new = gradTestFun(xy_new[0], xy_new[1])b = np.dot(grad_xy_new, (grad_xy_new - grad_xy))/np.dot(grad_xy, grad_xy)    #根据PRP公式计算组合系数d = b*d - grad_xy_newxy = xy_newarray_xy_history = np.array(ls_xy_history)plotResult(array_xy_history)return array_xy_history

运行结果为：

得到的最小值的坐标是[ 1.00318321, 1.50477481]
运行时间：

3.GD

最后还一起写了利用梯度下降法优化的结果以作对比，以下是实现代码：

def mainGD():'''使用梯度下降法计优化函数'''ls_xy_history = []xy0 = np.array([4.0, -2.0])ls_xy_history.append(xy0)grad_xy = gradTestFun(xy0[0], xy0[1])alpha = 1e-3      e0 = 1e-3xy = xy0while(1):tag_reach = np.abs(grad_xy)<e0if tag_reach.all():breakxy = xy - alpha*grad_xygrad_xy = gradTestFun(xy[0], xy[1])ls_xy_history.append(xy)array_xy_history = np.array(ls_xy_history)plotResult(array_xy_history)return array_xy_history

运行结果为：

得到的最小值坐标是[ 0.9185095 , 1.37763925]
运行时间是：

五、总结

对比可以发现，CG算法在计算的速度和准确度上都相较于GD算法有一定的优势。

python实现共轭梯度算法相关推荐

python实现共轭梯度算法（含误差与运算次数的折线图）
西安交通大学<数值分析>第三章课后题3.2 import numpy as np import math from math import log import matplotlib.py ...
共轭梯度算法求最小值-scipy
1 # coding=utf-8 2 3 #共轭梯度算法求最小值 4 import numpy as np 5 6 from scipy import optimize 7 8 9 10 11 def ...
共轭梯度算法之FR算法
共轭梯度算法之FR算法引理 FR算法算法步骤引理若f(x)=12xTGx+δTx+γf(x)=\frac{1}{2}x^TGx+\delta^T x+\gammaf(x)=21xTGx+δT ...
最优化方法----powell共轭梯度算法
Fletcher_Reeves共轭梯度算法要先计算梯度,然后对梯度进行修正得到一组共轭梯度.而计算梯度和计算函数值一样耗时.powell算法使用一维搜索技术寻找一组共轭梯度,从而不需要计算梯度. #i ...
机器学习：共轭梯度算法（PCG）
今天介绍数值计算和优化方法中非常有效的一种数值解法,共轭梯度法.我们知道,在解大型线性方程组的时候,很少会有一步到位的精确解析解,一般都需要通过迭代来进行逼近,而 PCG 就是这样一种迭代逼近算法. ...
MATLAB从入门到精通-最速下降算法、牛顿算法、BFGS拟牛顿算法、共轭梯度算法无约束极值问题
前言以下是我为大家准备的几个精品专栏,喜欢的小伙伴可自行订阅,你的支持就是我不断更新的动力哟! MATLAB-30天带你从入门到精通 MATLAB深入理解高级教程(附源码) tableau可视化数据 ...
机器学习基础（五十九）—— 高级优化算法（梯度下降、L-BFGS、共轭梯度）
优化算法两大核心,一曰:方向,比如由负梯度方向给出:二曰:步长. 在机器学习领域,不管是基础的梯度下降,还是更为高级的 L-BFGS.共轭梯度,都对应的是参数学习,用来学习相关模型的参数. 1. 梯度 ...
matlab ncg,2步非线性共轭梯度(NCG)法
2步非线性共轭梯度(NCG)法提出一种求解非线性方程组的2步非线性共轭梯度法(2步NCG法),并在一定条件下证明 (本文共9页) 阅读全文>> 作者提出的自适应共轭梯度方法是对共轭梯度方 ...
Opencv（python）图像梯度和边缘检测算法
1.图像梯度图像梯度计算的是图像的边缘信息 ,图像梯度计算的是图像变化的速度.对于图像的边缘部分,其灰度值变化较大,梯度值也较大,对于图像中比较平滑的部分,其灰度值变化较小,相应的梯度值也较小.图像 ...

python实现共轭梯度算法

python实现共轭梯度优化算法