python实现牛顿法_牛顿法和最速下降法的Python实现

1 牛顿法

1.1 牛顿法的Python程序

from sympy import *

import numpy as np

# 假设多元函数是二维形式

# x_init为二维向量(x1, x2)

def newton_dou(step, x_init, obj):

i = 1 # 记录迭代次数的变量

while i <= step:

if i == 1:

grandient_obj = np.array([diff(obj, x1).subs(x1, x_init[0]).subs(x2, x_init[1]), diff(obj, x2).subs(x1, x_init[0]).subs(x2, x_init[1])], dtype=float) # 初始点的梯度值

hessian_obj = np.array([[diff(obj, x1, 2), diff(diff(obj, x1), x2)], [diff(diff(obj, x2), x1), diff(obj, x2, 2)]], dtype=float) # 初始点的hessian矩阵

inverse = np.linalg.inv(hessian_obj) # hessian矩阵求逆

x_new = x_init - np.matmul(inverse, grandient_obj) # 第一次迭代公式

print(x_new)

# print('迭代第%d次：%.5f' %(i, x_new))

i = i + 1

else:

grandient_obj = np.array([diff(obj, x1).subs(x1, x_new[0]).subs(x2, x_new[1]), diff(obj, x2).subs(x1, x_new[0]).subs(x2, x_new[1])], dtype=float) # 当前点的梯度值

hessian_obj = np.array([[diff(obj, x1, 2), diff(diff(obj, x1), x2)], [diff(diff(obj, x2), x1), diff(obj, x2, 2)]], dtype=float) # 当前点的hessian矩阵

inverse = np.linalg.inv(hessian_obj) # hessian矩阵求逆

x_new = x_new - np.matmul(inverse, grandient_obj) # 迭代公式

print(x_new)

# print('迭代第%d次：%.5f' % (i, x_new))

i = i + 1

return x_new

x0 = np.array([0, 0], dtype=float)

x1 = symbols("x1")

x2 = symbols("x2")

newton_dou(5, x0, x1**2+2*x2**2-2*x1*x2-2*x2)

1.2 牛顿法的结果分析

程序执行的结果如下：

[1. 1.]

[1. 1.]

[1. 1.]

[1. 1.]

[1. 1.]

Process finished with exit code 0

经过实际计算函数

的极值点为

，只需一次迭代就能收敛到极值点。

2 最速下降法

2.1 最速下降法的Python程序

# !/usr/bin/python

# -*- coding:utf-8 -*-

from sympy import *

import numpy as np

import matplotlib.pyplot as plt

def sinvar(fun):

s_p = solve(diff(fun)) #stationary point

return s_p

def value_enter(fun, x, i):

value = fun[i].subs(x1, x[0]).subs(x2, x[1])

return value

def grandient_l2(grand, x_now):

grand_l2 = sqrt(pow(value_enter(grand, x_now, 0), 2)+pow(value_enter(grand, x_now, 1), 2))

return grand_l2

def msd(x_init, obj):

i = 1

grandient_obj = np.array([diff(obj, x1), diff(obj, x2)])

error = grandient_l2(grandient_obj, x_init)

plt.plot(x_init[0], x_init[1], 'ro')

while(error>0.001):

if i == 1:

grandient_obj_x = np.array([value_enter(grandient_obj, x_init, 0), value_enter(grandient_obj, x_init, 1)])

t = symbols('t')

x_eta = x_init - t * grandient_obj_x

eta = sinvar(obj.subs(x1, x_eta[0]).subs(x2, x_eta[1]))

x_new = x_init - eta*grandient_obj_x

print(x_new)

i = i + 1

error = grandient_l2(grandient_obj, x_new)

plt.plot(x_new[0], x_new[1], 'ro')

else:

grandient_obj_x = np.array([value_enter(grandient_obj, x_new, 0), value_enter(grandient_obj, x_new, 1)])

t = symbols('t')

x_eta = x_new - t * grandient_obj_x

eta = sinvar(obj.subs(x1, x_eta[0]).subs(x2, x_eta[1]))

x_new = x_new - eta*grandient_obj_x

print(x_new)

i = i + 1

error = grandient_l2(grandient_obj, x_new)

plt.plot(x_new[0], x_new[1], 'ro')

plt.show()

x_0 = np.array([0, 0], dtype=float)

x1 = symbols("x1")

x2 = symbols("x2")

result = msd(x_0, x1**2+2*x2**2-2*x1*x2-2*x2)

print(result)

2.2 最速下降法结果分析

程序执行的结果如下：

[0 0.500000000000000]

[0.500000000000000 0.500000000000000]

[0.500000000000000 0.750000000000000]

[0.750000000000000 0.750000000000000]

[0.750000000000000 0.875000000000000]

[0.875000000000000 0.875000000000000]

[0.875000000000000 0.937500000000000]

[0.937500000000000 0.937500000000000]

[0.937500000000000 0.968750000000000]

[0.968750000000000 0.968750000000000]

[0.968750000000000 0.984375000000000]

[0.984375000000000 0.984375000000000]

[0.984375000000000 0.992187500000000]

[0.992187500000000 0.992187500000000]

[0.992187500000000 0.996093750000000]

[0.996093750000000 0.996093750000000]

[0.996093750000000 0.998046875000000]

[0.998046875000000 0.998046875000000]

[0.998046875000000 0.999023437500000]

[0.999023437500000 0.999023437500000]

[0.999023437500000 0.999511718750000]

None

Process finished with exit code 0

最速下降法选择的停止条件为迭代的梯度值如果小于0.001，则停止迭代，通过结果可看出设定停止条件的情况下，极值点并没有完全收敛到

点，在Python中画出了迭代的极值点（如图1所示）。如果继续缩小停止条件，则程序会延长时间增加迭代次数，最后极值点达到

。

图1 最速下降法迭代的极值点