最速下降法python_最速下降法python实现

标签：

代码参数说明

其中d为下降方向，这里取grad(梯度)的负数方向，其中armijo搜索旨在寻找合适的步长

代码

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.ticker import LinearLocator, FormatStrFormatter

from mpl_toolkits.mplot3d import Axes3D

def sdm(fun, gfun, x0, rho, sigma, epsilon):

'''

最速下降法

:param fun: 目标函数

:param gfun: 梯度函数

:param x0: 初始点

:param rho: armijo搜索参数

:param sigma: 同上

:param epsilon: 终止调节参数

:return:

'''

max_iter_k = 5000

max_m = 20

k = 0

while k < max_iter_k:

grad = gfun(x0)

d = -grad

if np.linalg.norm(d) < epsilon:

break

m = 0

mk = 0

while m < max_m: # armijo 搜索

print('f(x + rho^m * d) = {}'.format(fun(x0 + pow(rho, m) * d)))

print('f(x) + sigma * rho^m * g * d = {}'.format(fun(x0) + sigma * pow(rho, m) * np.dot(grad.T, d)))

if fun(x0 + pow(rho, m) * d) < fun(x0) + sigma * pow(rho, m) * np.dot(grad.T, d):

mk = m

break

m += 1

x0 = x0 + pow(rho, mk) * d

k += 1

print('iterations : {}'.format(k))

return x0, fun(x0)

def obj(x):

'''

目标函数 课本p31

:param x:

:return:

'''

y = x[1]

x = x[0]

return 100 * pow(x * x - y, 2) + pow(x - 1, 2)

def obj_g(x):

y = x[1]

x = x[0]

arr = [400 * x * (x * x - y) + 2 * (x - 1), -200 * (x * x - y)]

return np.array(arr).T

def test_f(x):

'''

测试二元函数

:param x:

:return:

'''

y = x[1]

x = x[0]

return (x - 1)**2 + (y - 2)**2

def test_f_g(x):

y = x[1]

x = x[0]

arr = [2 * (x - 1), 2 * (y - 2)]

return np.array(arr).T

if __name__ == '__main__':

X = np.linspace(-3, 3, 100)

Y = np.linspace(-3, 3, 100)

X, Y = np.meshgrid(X, Y)

# Z = (X - 1)**2 + (Y - 2)**2

Z = 100 * (X**2 - Y)**2 + (X - 1)**2

fig = plt.figure()

ax = Axes3D(fig)

surf = ax.plot_surface(X, Y, Z, cmap=plt.cm.winter)

ax.set_xlabel("x-label", color='r')

ax.set_ylabel("y-label", color='g')

ax.set_zlabel("z-label", color='b')

plt.savefig('graph_for_function/obj_f.png')

plt.show()

x0 = np.array([0.0, 0.0]).T

# print(sdm(test_f, test_f_g, x0, 0.5, 0.4, 1e-5))

print(sdm(obj, obj_g, x0, 0.5, 0.4, 1e-5))

标签：