【求助】python实现最速下降算法的时候 同样的for循环两次出来的点竟然不一样

Question

import copy
import numpy as np
from functools import reduce

DELTA = 0.001


def fun(x):
        return x[0]**2+x[1]**2

def partial(f, x, d):
        _x = copy.copy(x)
        x_ = copy.copy(x)
        _x[d] -= DELTA
        x_[d] += DELTA
        return (f(x_) - f(_x))/(2*DELTA)

y=[]
def down(f,x,lr,decay):
        """
        多元函数f在积分区间interval上的定积分
        :param f: 函数
        :param x: 求梯度的位置，值为向量
        :param lr: 学习率
        :param decay: 衰减率
        :return: 最值的点的位置，值为向量
        """
        for i in range(2):
                t=1
                a=partial(f, x, i)
                x[i]-=a*lr/(1+decay*t)
                b=partial(f, x, i)

                while abs(a-b)>0.00000001:
                        t=t+1
                        a=b
                        x[i]-=a*lr/(1+decay*t)
                        b=partial(f, x, i)
                y.append(x[i])

        return y

z=[9,9]
r=down(fun,z,1,0.5)
h=fun(r)
print(r)
print(h)

python实现最速下降算法的时候，同样的for循环两次出来的点竟然不一样，
不知道为什么求助、、

Answer 1

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

抄来的 您看看
https://blog.csdn.net/qq_26479655/article/details/83038685