线搜索之用backtracking确定步长（牛顿法、最速下降法）

用某数学软件得Rosenbrock函数：
$$f(x,y)=100\ast (y-x^2)\ast 2+(1-x{)}^2$$的最小值为0.其函数图像为：

但两种线搜索方法的结果差距还挺大的。

import numpy as np

#目标函数
def function(x):
    x1,x2 = x[0], x[1]
    f = 100*pow((x2 - pow(x1,2)),2)+pow((1-x1),2)
    return f

#最速下降的梯度方向
def steepest_decent_gradient(x):
    h = 1e-4
    grad = np.zeros_like(x)
    for i in range(x.size):
        temp = x[i]
        x[i] = temp + h
        fh1 = function(x)
        x[i] = temp - h
        fh2 = function(x)
        grad[i] = (fh1-fh2)/(2*h)
        x[i] = temp
    return grad

#求偏导
def partical_derivative(target, x):
    h = 1e-4
    x_add, x_subtract = x.copy(), x.copy()
    x_add[target] = x_add[target] + h
    x_subtract[target] = x_subtract[target] - h
    result = (function(x_add)-function(x_subtract))/(2*h)
    return result

#牛顿法的梯度方向
def newton_gradient(x):
    A = np.zeros((2,2))
    b = steepest_decent_gradient(x)
    m, n = A.shape
    h = 1e-4
    #写一阶导的函数，然后调用
    for i in range(m):
        for j in range(n):
            x_add, x_subtract = x.copy(), x.copy()
            x_add[j] = x_add[j] + h
            x_subtract[j] = x_subtract[j] - h
            A[i][j] = (partical_derivative(i,x_add) - partical_derivative(i,x_subtract))/(2*h)
    result = np.linalg.solve(A, -1*b)
    return result

#backtracking
def backtracking(x,alpha,row,c,pk,gradient):
    '''
    d = newton_gradient(function(x),x)
    b = np.dot(d,pk)
    con = function(x) + c*alpha*b
    '''
    while function(x+alpha*pk) > function(x) + c*alpha*np.dot(gradient(x),x):
        alpha = row * alpha
    return alpha

#用backtracking的最速下降
def steepest_decent(x0,num):
    steps = []
    for i in range(num):
        pk = steepest_decent_gradient(x0)
        alpha = backtracking(x0,1,0.5,0.8,pk,steepest_decent_gradient)
        steps.append(alpha)
        x1 = x0 - alpha*pk
        if abs(function(x1)-function(x0)) <= 1e-4:
            print(steps)
            return print("找到极小值：{}".format(function(x1)))
        x0 = x1
    print(steps)
    print("迭代{}次内未找到函数极小值。".format(num))
    return 0

#用backtracking的牛顿方法
def newton(x0,num):
    steps = []
    for i in range(num):
        pk = newton_gradient(x0)
        alpha = backtracking(x0,1,0.5,0.8,pk,newton_gradient)
        steps.append(alpha)
        x1 = x0 + alpha*pk
        if abs(function(x1) - function(x0)) <= 1e-4:
            print(steps)
            return print("找到极小值：{}".format(function(x1)))
        x0 = x1
    print(steps)
    print("迭代{}次内未找到函数极小值。".format(num))
    return 0

if __name__ == '__main__':
    x0 = np.array((-1.2,1))
    x1 = np.array((1.2,1.2))
    print('最速下降法，初值为{}'.format(x0))
    steepest_decent(x0,20)
    print('牛顿法，初值为{}'.format(x0))
    newton(x0,20)
    print('最速下降法，初值为{}'.format(x1))
    steepest_decent(x1, 20)
    print('牛顿法，初值为{}'.format(x1))
    newton(x1, 20)

结果为：
最速下降法，初值为[-1.2 1. ]
[4.336808689942018e-19]
找到极小值：24.199999999999996
牛顿法，初值为[-1.2 1. ]
[1, 0.0625, 0.25, 0.5, 1, 1, 1, 1, 0.5, 1, 0.5, 1, 0.5, 1, 1, 1, 1, 1, 1]
找到极小值：2.738165819482374e-10
最速下降法，初值为[1.2 1.2]
[8.673617379884035e-19]
找到极小值：5.8
牛顿法，初值为[1.2 1.2]
[1, 3.469446951953614e-18]
找到极小值：0.03838401538193726
从结果来看，在这个目标函数下，backtracking不适合最速下降。