使用Python实现牛顿法求极值

对于一个多元函数 用牛顿法求其极小值的迭代格式为

其中 为函数 的梯度向量， 为函数 的Hesse（Hessian）矩阵。

上述牛顿法不是全局收敛的。为此可以引入阻尼牛顿法（又称带步长的牛顿法）。

我们知道，求极值的一般迭代格式为

其中 为搜索步长， 为搜索方向（注意所有的迭代格式都是先计算搜索方向，再计算搜索步长，如同瞎子下山一样，先找到哪个方向可行下降，再决定下几步）。

取下降方向 即得阻尼牛顿法，只不过搜索步长 不确定，需要用线性搜索技术确定一个较优的值，比如精确线性搜索或者Goldstein搜索、Wolfe搜索等。特别地，当 一直取为常数1时，就是普通的牛顿法。

以Rosenbrock函数为例，即有

于是可得函数的梯度

函数 的Hesse矩阵为

编写Python代码如下（使用版本为Python3.3）：

"""
Newton法
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

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

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

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 newton(x0):

 print('初始点为:')
 print(x0,'\n')
 W=np.zeros((2,10**3))
 i = 1
 imax = 1000
 W[:,0] = x0
 x = x0
 delta = 1
 alpha = 1

 while i<imax and delta>10**(-5):
  p = -np.dot(np.linalg.inv(hessian(x)),jacobian(x))
  x0 = x
  x = x + alpha*p
  W[:,i] = x
  delta = sum((x-x0)**2)
  print('第',i,'次迭代结果:')
  print(x,'\n')
  i=i+1
 W=W[:,0:i] # 记录迭代点
 return W

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

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

上述代码中jacobian(x)返回函数的梯度，hessian(x)返回函数的Hesse矩阵，用W矩阵记录迭代点的坐标，然后画出点的搜索轨迹。

可得输出结果为

初始点为:
[-1.2 1. ] 

第 1 次迭代结果:
[-1.1752809 1.38067416] 

第 2 次迭代结果:
[ 0.76311487 -3.17503385] 

第 3 次迭代结果:
[ 0.76342968 0.58282478] 

第 4 次迭代结果:
[ 0.99999531 0.94402732] 

第 5 次迭代结果:
[ 0.9999957 0.99999139] 

第 6 次迭代结果:
[ 1. 1.]

即迭代了6次得到了最优解，画出的迭代点的轨迹如下：

由于主要使用了Python的Numpy模块来进行计算，可以看出，代码和最终的图与Matlab是很相像的。

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

赞 (0)

相关推荐

  • TensorFlow 多元函数的极值实例

    flyfish python实现 设函数 的某个邻域内有定义,对于该邻域内异于的点,如果都适合不等式 则称函数在点有极大值. 如果都适合不等式 则称函数在点有极小值． 极大值.极小值统称为极值.使函数取得极值的点称为极值点. 有极小值的例子 函数 在点(0,0)处有极小值.因为对于点 (0,0)的任一邻域内异于(0,0)的点,函数值都为正,而在点(0,0)处的函数值为零.从几何上看这是显然的,因为点(0,0,0)是开口朝上的椭圆抛物面 的顶点. 代码 from matplotlib import

  • Python使用min、max函数查找二维数据矩阵中最小、最大值的方法

    本文实例讲述了Python使用min.max函数查找二维数据矩阵中最小.最大值的方法.分享给大家供大家参考,具体如下: 简单使用min.max函数来得到二维数据矩阵中的最大最小值,很简单,这是因为工作需要用到一个东西所以先简单来写了一下: #!usr/bin/env python #encoding:utf-8 ''''' __Author__:沂水寒城 功能:找出来随机生成矩阵中的最大.最小值 ''' import time import random def random_matrix_ge

  • python求绝对值的三种方法小结

    如下所示: 1.条件判断 2.内置函数abs() 3.内置模块 math.fabs abs() 与fabs()的区别 abs()是一个内置函数,而fabs()在math模块中定义的. fabs()函数只适用于float和integer类型,而abs()也适用于复数. abs()返回是float和int类型,math.fabs()返回是float类型 以上这篇python求绝对值的三种方法小结就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持我们.

  • python求最大值最小值方法总结

    python求最大值最小值方法总结

    方法一(常规): 代码: count = int(input('输入数据个数:\n')) a = 1 while a <= count: num = int(input('请输入第{}个数:'.format(a))) #字符串中的方法 if a == 1: #这句一定会执行,而且只执行一次,目的就是让你输入的第一个数作为根据与之后的数比较 max = min = num #第二个及以后的数都会走else, else: #第一次走else时,比较中的min和max都是你第一次输入的数,以后走els

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

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

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

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

  • Python算法之求n个节点不同二叉树个数

    问题 创建一个二叉树 二叉树有限多个节点的集合,这个集合可能是: 空集 由一个根节点,和两棵互不相交的,分别称作左子树和右子树的二叉树组成 创建二叉树: 创建节点 再创建节点之间的关系 Python代码示例 # !/usr/bin/env python # -*-encoding: utf-8-*- # author:LiYanwei # version:0.1 class TreeNode(object): def __init__ (self, data, left = None, righ

  • 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 Numpy中求向量和矩阵的范数实例

    np.linalg.norm(求范数):linalg=linear(线性)+algebra(代数),norm则表示范数. 函数参数 x_norm=np.linalg.norm(x, ord=None, axis=None, keepdims=False) ①x: 表示矩阵(也可以是一维) ②ord:范数类型 向量的范数: 矩阵的范数: ord=1:列和的最大值 ord=2:|λE-ATA|=0,求特征值,然后求最大特征值得算术平方根 ord=∞:行和的最大值 ③axis:处理类型 axis=1表

  • python、Matlab求定积分的实现

    python求定积分 计算 from sympy import * x = symbols('x') print(integrate(sin(2*x)/(1+x**2), (x, 0, 3))) sympy库中integrate函数 integrate(f, (x, lower_bound, upper_bound)) # f-函数,x-变量,lower_bound-下限,upper_bound-上限 但是发现求不出来,如果是sin(2*x)就可以,为什么? matlab求定积分 syms x

  • Python基于递归算法求最小公倍数和最大公约数示例

    本文实例讲述了Python基于递归算法求最小公倍数和最大公约数.分享给大家供大家参考,具体如下: # 最小公倍数 def lcm(a, b, c=1): if a * c % b != 0: return lcm(a, b, c+1) else: return a*c test_cases = [(4, 8), (35, 42), (5, 7), (20, 10)] for case in test_cases: print('lcm of {} is {}'.format(*case, lcm

  • Python根据欧拉角求旋转矩阵的实例

    利用numpy和scipy,我们可以很容易根据欧拉角求出旋转矩阵,这里的旋转轴我们你理解成四元数里面的旋转轴 import numpy as np import scipy.linalg as linalg import math #参数分别是旋转轴和旋转弧度值 def rotate_mat(self, axis, radian): rot_matrix = linalg.expm(np.cross(np.eye(3), axis / linalg.norm(axis) * radian)) a

  • 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实现迭代法求方程组的根过程解析

    这篇文章主要介绍了python实现迭代法求方程组的根过程解析,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 有方程组如下: 迭代法求解x,python代码如下: import numpy as np import matplotlib.pyplot as plt A = np.array([[8, -3, 2], [4, 11, -1], [6, 3, 12]]) b = np.array([[20, 33, 36]]) # 方法一:消元法求解

随机推荐