用Python实现最速下降法求极值的方法

对于一个多元函数，用最速下降法（又称梯度下降法）求其极小值的迭代格式为

其中为负梯度方向，即最速下降方向，αkαk为搜索步长。

一般情况下，最优步长αkαk的确定要用到线性搜索技术，比如精确线性搜索，但是更常用的是不精确线性搜索，主要是Goldstein不精确线性搜索和Wolfe法线性搜索。

为了调用的方便，编写一个Python文件，里面存放线性搜索的子函数，命名为linesearch.py，这里先只编写了Goldstein线性搜索的函数，关于Goldstein原则，可以参看最优化课本。

线性搜索的代码如下（使用版本为Python3.3）：

'''
线性搜索子函数
'''

import numpy as np
import random

def goldsteinsearch(f,df,d,x,alpham,rho,t):

  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
    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

上述函数的输入参数主要包括一个多元函数f，其导数df，当前迭代点x和当前搜索方向d，返回值是根据Goldstein准则确定的搜索步长。

我们仍以Rosenbrock函数为例，即有

于是可得函数的梯度为

最速下降法的代码如下：

"""
最速下降法
Rosenbrock函数
函数 f(x)=100*(x(2)-x(1).^2).^2+(1-x(1)).^2
梯度 g(x)=(-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1)),200*(x(2)-x(1)^2))^(T)
"""

import numpy as np
import matplotlib.pyplot as plt
import random
import linesearch
from linesearch import goldsteinsearch

def rosenbrock(x):
  return 100*(x[1]-x[0]**2)**2+(1-x[0])**2

def jacobian(x):
  return np.array([-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0]),200*(x[1]-x[0]**2)])

X1=np.arange(-1.5,1.5+0.05,0.05)
X2=np.arange(-3.5,2+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条轮廓线

def steepest(x0):

  print('初始点为:')
  print(x0,'\n')
  imax = 20000
  W=np.zeros((2,imax))
  W[:,0] = x0
  i = 1
  x = x0
  grad = jacobian(x)
  delta = sum(grad**2) # 初始误差

  while i<imax and delta>10**(-5):
    p = -jacobian(x)
    x0=x
    alpha = goldsteinsearch(rosenbrock,jacobian,p,x,1,0.1,2)
    x = x + alpha*p
    W[:,i] = x
    grad = jacobian(x)
    delta = sum(grad**2)
    i=i+1

  print("迭代次数为:",i)
  print("近似最优解为:")
  print(x,'\n')
  W=W[:,0:i] # 记录迭代点
  return W

x0 = np.array([-1.2,1])
W=steepest(x0)

plt.plot(W[0,:],W[1,:],'g*',W[0,:],W[1,:]) # 画出迭代点收敛的轨迹
plt.show()

为了实现不同文件中函数的调用，我们先用import函数导入了线性搜索的子函数，也就是下面的2行代码

import linesearch
from linesearch import goldsteinsearch

当然，如果把定义goldsteinsearch函数的代码直接放到程序里面，就不需要这么麻烦了，但是那样的话，不仅会使程序显得很长，而且不便于goldsteinsearch函数的重用。

此外，Python对函数式编程也支持的很好，在定义goldsteinsearch函数时，可以允许抽象的函数f，df作为其输入参数，只要在调用时实例化就可以了。与Matlab不同的是，传递函数作为参数时，Python是不需要使用@将其变为函数句柄的。

运行结果为

初始点为:

[-1.2 1. ] 

迭代次数为: 1504

近似最优解为:

[ 1.00318532 1.00639618]

迭代点的轨迹为

由于在线性搜索子程序中使用了随机函数，初始搜索点是随机产生的，因此每次运行的结果不太相同，比如再运行一次程序，得到

初始点为:
[-1.2 1. ] 

迭代次数为: 1994

近似最优解为:
[ 0.99735222 0.99469882]

所得图像为

以上这篇用Python实现最速下降法求极值的方法就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持我们。

赞 (0)

相关推荐

  • python 梯度法求解函数极值的实例

    如下所示: #coding utf-8 a=0.001 #定义收敛步长 xd=1 #定义寻找步长 x=0 #定义一个种子x0 i=0 #循环迭代次数 y=0 dic={} import math def f(x): y=math.sin(x) #定义函数f(X)=sinx return y def fd(x): y=math.cos(x) #函数f(x)导数fd(X)=cosx return y while y>=0 and y<3.14*4: y=y+xd x=y while abs(fd(

  • python 寻找离散序列极值点的方法

    使用 scipy.signal 的 argrelextrema 函数(API),简单方便 import numpy as np import pylab as pl import matplotlib.pyplot as plt import scipy.signal as signal x=np.array([ 0, 6, 25, 20, 15, 8, 15, 6, 0, 6, 0, -5, -15, -3, 4, 10, 8, 13, 8, 10, 3, 1, 20, 7, 3, 0 ])

  • 梯度下降法介绍及利用Python实现的方法示例

    梯度下降法介绍及利用Python实现的方法示例

    本文主要给大家介绍了梯度下降法及利用Python实现的相关内容,分享出来供大家参考学习,下面话不多说,来一起看看详细的介绍吧. 梯度下降法介绍 梯度下降法(gradient descent),又名最速下降法(steepest descent)是求解无约束最优化问题最常用的方法,它是一种迭代方法,每一步主要的操作是求解目标函数的梯度向量,将当前位置的负梯度方向作为搜索方向(因为在该方向上目标函数下降最快,这也是最速下降法名称的由来). 梯度下降法特点:越接近目标值,步长越小,下降速度越慢. 直观上

  • python梯度下降法的简单示例

    梯度下降法的原理和公式这里不讲,就是一个直观的.易于理解的简单例子. 1.最简单的情况,样本只有一个变量,即简单的(x,y).多变量的则可为使用体重或身高判断男女(这是假设,并不严谨),则变量有两个,一个是体重,一个是身高,则可表示为(x1,x2,y),即一个目标值有两个属性. 2.单个变量的情况最简单的就是,函数hk(x)=k*x这条直线(注意:这里k也是变化的,我们的目的就是求一个最优的   k).而深度学习中,我们是不知道函数的,也就是不知道上述的k.   这里讨论单变量的情况: 在不知道

  • 用Python实现最速下降法求极值的方法

    对于一个多元函数,用最速下降法(又称梯度下降法)求其极小值的迭代格式为 其中为负梯度方向,即最速下降方向,αkαk为搜索步长. 一般情况下,最优步长αkαk的确定要用到线性搜索技术,比如精确线性搜索,但是更常用的是不精确线性搜索,主要是Goldstein不精确线性搜索和Wolfe法线性搜索. 为了调用的方便,编写一个Python文件,里面存放线性搜索的子函数,命名为linesearch.py,这里先只编写了Goldstein线性搜索的函数,关于Goldstein原则,可以参看最优化课本. 线性搜

  • 使用Python实现牛顿法求极值

    对于一个多元函数 用牛顿法求其极小值的迭代格式为 其中 为函数 的梯度向量, 为函数 的Hesse(Hessian)矩阵. 上述牛顿法不是全局收敛的.为此可以引入阻尼牛顿法(又称带步长的牛顿法). 我们知道,求极值的一般迭代格式为 其中 为搜索步长, 为搜索方向(注意所有的迭代格式都是先计算搜索方向,再计算搜索步长,如同瞎子下山一样,先找到哪个方向可行下降,再决定下几步). 取下降方向 即得阻尼牛顿法,只不过搜索步长 不确定,需要用线性搜索技术确定一个较优的值,比如精确线性搜索或者Goldste

  • Python求导数的方法

    本文实例讲述了Python求导数的方法.分享给大家供大家参考.具体实现方法如下: def func(coeff): sum='' for key in coeff: sum=sum+'+'+str(key)+'*'+'x'+'**'+str(coeff[key]) return sum[1:] from sympy import * from sympy.core.sympify import SympifyError expr = func({2:0,3:1,4:2,5:7}) x = Sym

  • python求pi的方法

    本文实例讲述了python求pi的方法,是一篇翻译自国外网站的文章,分享给大家供大家参考. 具体实现方法如下: #_*_ coding=utf-8 *_* ## {{{ http://code.activestate.com/recipes/578130/ (r5) def pi(places=10): """Computes pi to given number of decimal places 参数places表示要返回的pi的小数点后位数 方法:先整体扩大10**8(

  • Python编程使用*解包和itertools.product()求笛卡尔积的方法

    本文实例讲述了Python编程使用*解包和itertools.product()求笛卡尔积的方法.分享给大家供大家参考,具体如下: [问题] 目前有一字符串s = "['a', 'b'],['c', 'd']",想把它分开成为两个列表: list1 = ['a', 'b'] list2 = ['c', 'd'] 之后使用itertools.product()求笛卡尔积,应该写成: for i in itertools.product(list1, list2): print i 结果为

  • Python Dataframe 指定多列去重、求差集的方法

    1)去重 指定多列去重,这是在dataframe没有独一无二的字段作为PK(主键)时,需要指定多个字段一起作为该行的PK,在这种情况下对整体数据进行去重. Attention:主要用到了drop_duplicates方法,并设置参数subset为多个字段名构成的数组. 具体代码如下: >>>import pandas as pd >>>data={'state':[1,1,2,2,1,2,2],'pop':['a','b','c','d','b','c','d']} &

  • python分治法求二维数组局部峰值方法

    题目的意思大致是在一个n*m的二维数组中,找到一个局部峰值.峰值要求大于相邻的四个元素(数组边界以外视为负无穷),比如最后我们找到峰值A[j][i],则有A[j][i] > A[j+1][i] && A[j][i] > A[j-1][i] && A[j][i] > A[j][i+1] && A[j][i] > A[j][i-1].返回该峰值的坐标和值. 当然,最简单直接的方法就是遍历所有数组元素,判断是否为峰值,时间复杂度为O(n^2

  • python计算波峰波谷值的方法（极值点）

    python求极值点主要用到scipy库. 1. 首先可先选择一个函数或者拟合一个函数,这里选择拟合数据:np.polyfit import pandas as pd import matplotlib.pyplot as plt import numpy as np from scipy import signal #滤波等 xxx = np.arange(0, 1000) yyy = np.sin(xxx*np.pi/180) z1 = np.polyfit(xxx, yyy, 7) # 用

  • Python实现约瑟夫环问题的方法

    本文实例讲述了Python实现约瑟夫环问题的方法.分享给大家供大家参考,具体如下: 题目:0,1,...,n-1这n个数字排成一个圆圈,从数字0开始每次从这个圆圈里删除第m个数字.求出这个圆圈里剩下的最后一个数字. 定义函数f(n,m),表示每次在n个数字(0,1,...,n-1)中每次删除第m个数字后最后剩下的数字. 在n个数字中,假设第一个被删除的数字为k,那么删除k之后剩下的n-1个数字为0~k-1,k 1~n-1,并且下一次删除从数字k 1开始计数.第二个序列最后剩下的数字也就是我们要求

  • python安装与使用redis的方法

    本文实例讲述了python安装与使用redis的方法.分享给大家供大家参考,具体如下: 1.安装 好吧,我承认我只会最简单的安装: sudo apt-get install redis-server python 支持包: (其实就一个文件,搞过来就能用) sudo apt-get install python-redis 2.配置 配置一下吧,默认配置文件在: "/etc/redis/redis.conf" 绑定ip: "bind 127.0.0.1″ -> &quo

随机推荐