python实现共轭梯度法

 更新时间:2019年07月03日 14:09:53   作者:tomator01   我要评论

这篇文章主要介绍了python实现的共轭梯度法,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧

共轭梯度法是介于最速下降法与牛顿法之间的一个方法,它仅需利用一阶导数信息,但克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存储和计算hesse矩阵并求逆的缺点,共轭梯度法不仅是解决大型线性方程组最有用的方法之一,也是解大型非线性最优化最有效的算法之一。 在各种优化算法中,共轭梯度法是非常重要的一种。其优点是所需存储量小,具有步收敛性,稳定性高,而且不需要任何外来参数。

算法步骤:

import random
import numpy as np
import matplotlib.pyplot as plt
 
def goldsteinsearch(f,df,d,x,alpham,rho,t):
 '''
 线性搜索子函数
 数f,导数df,当前迭代点x和当前搜索方向d,t试探系数>1,
 '''
 flag = 0
 
 a = 0
 b = alpham
 fk = f(x)
 gk = df(x)
 
 phi0 = fk
 dphi0 = np.dot(gk, d)
 alpha=b*random.uniform(0,1)
 
 while(flag==0):
  newfk = f(x + alpha * d)
  phi = newfk
  # print(phi,phi0,rho,alpha ,dphi0)
  if (phi - phi0 )<= (rho * alpha * dphi0):
   if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
    flag = 1
   else:
    a = alpha
    b = b
    if (b < alpham):
     alpha = (a + b) / 2
    else:
     alpha = t * alpha
  else:
   a = a
   b = alpha
   alpha = (a + b) / 2
 return alpha
 
 
def wolfesearch(f,df,d,x,alpham,rho,t):
 '''
 线性搜索子函数
 数f,导数df,当前迭代点x和当前搜索方向d
 σ∈(ρ,1)=0.75
 '''
 sigma=0.75
 
 flag = 0
 
 a = 0
 b = alpham
 fk = f(x)
 gk = df(x)
 
 phi0 = fk
 dphi0 = np.dot(gk, d)
 alpha=b*random.uniform(0,1)
 
 while(flag==0):
  newfk = f(x + alpha * d)
  phi = newfk
  # print(phi,phi0,rho,alpha ,dphi0)
  if (phi - phi0 )<= (rho * alpha * dphi0):
   # if abs(np.dot(df(x + alpha * d),d))<=-sigma*dphi0:
   if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
    flag = 1
   else:
    a = alpha
    b = b
    if (b < alpham):
     alpha = (a + b) / 2
    else:
     alpha = t * alpha
  else:
   a = a
   b = alpha
   alpha = (a + b) / 2
 return alpha
 
def frcg(fun,gfun,x0):
 
 
 # x0是初始点,fun和gfun分别是目标函数和梯度
 # x,val分别是近似最优点和最优值,k是迭代次数
 # dk是搜索方向,gk是梯度方向
 # epsilon是预设精度,np.linalg.norm(gk)求取向量的二范数
 maxk = 5000
 rho = 0.6
 sigma = 0.4
 k = 0
 epsilon = 1e-5
 n = np.shape(x0)[0]
 itern = 0
 w = np.zeros((2, 20000))
 
 f = open("共轭.txt", 'w')
 
 while k < maxk:
   w[:, k] = x0
   gk = gfun(x0)
   itern += 1
   itern %= n
   if itern == 1:
    dk = -gk
   else:
    beta = 1.0 * np.dot(gk, gk) / np.dot(g0, g0)
    dk = -gk + beta * d0
    gd = np.dot(gk, dk)
    if gd >= 0.0:
     dk = -gk
   if np.linalg.norm(gk) < epsilon:
    break
 
   alpha=goldsteinsearch(fun,gfun,dk,x0,1,0.1,2)
   # alpha=wolfesearch(fun,gfun,dk,x0,1,0.1,2)
   x0+=alpha*dk
 
   f.write(str(k)+' '+str(np.linalg.norm(gk))+"\n")
   print(k,alpha)
   g0 = gk
   d0 = dk
   k += 1
 
 w = w[:, 0:k+1] # 记录迭代点
 return [x0, fun(x0), k,w]
 
def fun(x):
 return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2
def gfun(x):
 return np.array([-400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0]), 200 * (x[1] - x[0] ** 2)])
 
 
if __name__=="__main__":
 x1 = np.arange(-1.5, 1.5 + 0.05, 0.05)
 x2 = np.arange(-3.5, 4 + 0.05, 0.05)
 [x1, x2] = np.meshgrid(x1, x2)
 f = 100 * (x2 - x1 ** 2) ** 2 + (1 - x1) ** 2 # 给定的函数
 plt.contour(x1, x2, f, 20) # 画出函数的20条轮廓线
 
 x0 = np.array([-1.2, 1])
 x=frcg(fun,gfun,x0)
 print(x[0],x[2])
 # [1.00318532 1.00639618]
 w=x[3]
 # print(w[:, :])
 plt.plot(w[0, :], w[1, :], 'g*-') # 画出迭代点收敛的轨迹
 plt.show()

代码中求最优步长用得是goldsteinsearch方法,另外的wolfesearch是试验的部分,在本段程序中不起作用。

迭代轨迹:

三种最优化方法的迭代次数对比:

最优化方法

最速下降法

共轭梯度法

牛顿法

迭代次数

1702

240

5

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持脚本之家。

相关文章

  • 利用aardio给python编写图形界面

    利用aardio给python编写图形界面

    这篇文章主要为大家详细介绍了利用aardio给python编写图形界面的方法,具有一定的参考价值,感兴趣的小伙伴们可以参考一下
    2017-08-08
  • python cookie 读取和保存方法

    python cookie 读取和保存方法

    今天小编就为大家分享一篇python cookie 读取和保存方法,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧
    2018-12-12
  • python通过tcp发送xml报文的方法

    python通过tcp发送xml报文的方法

    今天小编就为大家分享一篇python通过tcp发送xml报文的方法,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧
    2018-12-12
  • python实现在tkinter中使用matplotlib绘制图形的方法示例

    python实现在tkinter中使用matplotlib绘制图形的方法示例

    这篇文章主要介绍了python实现在tkinter中使用matplotlib绘制图形的方法,结合实例形式分析了python使用tkinter与matplotlib进行正弦曲线图形绘制的相关操作技巧,需要的朋友可以参考下
    2018-01-01
  • python实现floyd算法

    python实现floyd算法

    这篇文章主要为大家详细介绍了python实现floyd算法,具有一定的参考价值,感兴趣的小伙伴们可以参考一下
    2018-01-01
  • python输入整条数据分割存入数组的方法

    python输入整条数据分割存入数组的方法

    今天小编就为大家分享一篇python输入整条数据分割存入数组的方法,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧
    2018-11-11
  • 剖析django中模版标签的解析与参数传递

    剖析django中模版标签的解析与参数传递

    这篇文章主要介绍了剖析django中模版标签的解析与参数传递,django是重多高人气python框架中最为著名的一个,需要的朋友可以参考下
    2015-07-07
  • 利用python和ffmpeg 批量将其他图片转换为.yuv格式的方法

    利用python和ffmpeg 批量将其他图片转换为.yuv格式的方法

    今天小编就为大家分享一篇利用python和ffmpeg 批量将其他图片转换为.yuv格式的方法,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧
    2019-01-01
  • 3个用于数据科学的顶级python库

    3个用于数据科学的顶级python库

    今天小编就为大家分享一篇关于3个用于数据科学的顶级python库,小编觉得内容挺不错的,现在分享给大家,具有很好的参考价值,需要的朋友一起跟随小编来看看吧
    2018-09-09
  • python实现过滤单个androids程序日志脚本分享

    python实现过滤单个androids程序日志脚本分享

    这篇文章主要介绍了python实现过滤单个androids程序日志脚本分享,本文讲解了原理、实现代码、使用方法、最新代码等内容,需要的朋友可以参考下
    2015-01-01

最新评论