Python最速下降法实例

最速下降法的实现需要通过符号计算。

首先笔算一步如下，然后通过程序验证：

 

python程序如下，需要pip install sympy：

import numpy as np
from sympy import *
import math
import matplotlib.pyplot as plt# 定义符号
x1, x2, t = symbols('x1, x2, t')

def func():
    # 自定义一个函数
    return (x1-2) * (x1-2) + 4 * (x2-3) * (x2-3)

def grad(data):
    f = func()
    grad_vec = [diff(f, x1), diff(f, x2)]
    grad = []
    for item in grad_vec:
        grad.append(item.subs(x1, data[0]).subs(x2, data[1]))
    return grad

def grad_len(grad):
    vec_len = math.sqrt(pow(grad[0], 2) + pow(grad[1], 2))
    return vec_len

def zhudian(f):
    t_diff = diff(f)
    t_min = solve(t_diff)
    return t_min

绘制目标函数的图形如下：

def func_numeric(x):
    # 自定义一个函数
    return (x[0]-2) * (x[0]-2) + 4 * (x[1]-3) * (x[1]-3)
plt.figure(3)
plt_x = np.arange(0, 4, 0.05)
plt_y = np.arange(0, 4, 0.05)
X, Y = np.meshgrid(plt_x, plt_y)
Z1 = func_numeric([X, Y])
ax = plt.gca(projection='3d')
ax.plot_surface(X, Y, Z1, cmap=plt.get_cmap('rainbow'), linewidth=0.2)
plt.title("steepest descent method")
plt.xlabel("x1")
plt.ylabel("x2")
plt.show()

进行迭代计算：

X0 = [0, 0]
theta = 0.00001

data_x = [X0[0]]
data_y = [X0[1]]
data_f = [func().subs(x1, X0[0]).subs(x2, X0[1])]

f = func()
grad_vec = grad(X0)
print(f"grad: {grad_vec}")
grad_length = grad_len(grad_vec)  # 梯度向量的模长
print('grad_length', grad_length)
k = 1
while grad_length > theta:  # 迭代的终止条件
    k += 1
    p = -np.array(grad_vec)
    # 迭代
    X = np.array(X0) + t*p
    t_func = f.subs(x1, X[0]).subs(x2, X[1])
    t_min = zhudian(t_func)
    print(f"t_min: {t_min}")
    X0 = np.array(X0) + t_min*p
    print(f"X0: {X0}")
    # print(floatat(X0[0]))
    print(f"fx0: {f.subs(x1, float(X0[0])).subs(x2, float(X0[1]))}")
    grad_vec = grad(X0)
    print(f"grad: {grad_vec}")
    grad_length = grad_len(grad_vec)
    print('grad_length', grad_length)
    print('坐标', X0[0], X0[1])
    data_x.append(X0[0])
    data_y.append(X0[1])
    data_f.append(f.subs(x1, float(X0[0])).subs(x2, float(X0[1])))

打印的第一步如下：

grad: [-4, -24]
grad_length 24.331050121192877
t_min: [37/290]
X0: [74/145 444/145]
fx0: 2.2344827586206

和笔算结果一致。

绘制等高线图和迭代过程折线图：

# 绘图
fig = plt.figure(4)
plt.title(r'$Gradient \ method - steepest \ descent \ method$')
plt.plot(data_x, data_y, color='k', label=r'$f(x_1,x_2)=(x_1-2)^2+4*(x_2-3)^2$')

def func_numeric(x):
    # 自定义一个函数
    return (x[0]-2) * (x[0]-2) + 4 * (x[1]-3) * (x[1]-3)
plt_x = np.arange(0, 4, 0.05)
plt_y = np.arange(0, 4, 0.05)
X, Y = np.meshgrid(plt_x, plt_y)
Z1 = func_numeric([X, Y])
# ctf = plt.contourf(X, Y, Z1, 15)
ct = plt.contour(X, Y, Z1)
plt.clabel(ct, inline=True, fontsize=10)
# plt.colorbar(ctf)

plt.legend()
# plt.scatter(1, 1, marker=(5, 1), c=5, s=1000)
# plt.grid()
plt.xlabel(r'$x_1$', fontsize=20)
plt.ylabel(r'$x_2$', fontsize=20)
plt.show()

如图所示，结果正确：