最速下降法（梯度下降法）python实现

大家好，欢迎来到IT知识分享网。

有用请点赞，没用请差评。

欢迎分享本文，转载请保留出处。

最近在写论文，做的是共轭梯度反演方法，所以将论文中的部分内容分享出来吧。

大部分的机器学习算法的本质都是建立优化模型，通过最优化方法对目标函数（或损失函数）进行优化，从而训练出最好的模型。常见的最优化方法有梯度下降法、牛顿法和拟。牛顿法、共轭梯度法等等。在数学领域称为最优化，在地球物理领域也称为反演。

关于最速下降法的原理这里就不赘述了，本文也是实现的最基本的最速下降法，当然还有随机梯度下降等改进后的算法。

当目标函数是凸函数时，梯度下降法的解是全局最优解。一般情况下，其解不保证是全局最优解。梯度下降法的收敛速度也未必是很快的。

缺点：数值实验表明，当目标函数的等值线接近一个圆（球）时，最速下降法下降很快，当目标函数的等值线是一个扁长的椭球时，最速下降法开始几步下降较快，后来就出现锯齿线性，下降就十分缓慢。原因是在相邻两个迭代点上的函数f(x)的两个梯度是互相正交的，即，两个搜索方向互相正交，就产生了锯齿性质。当接近极小点时，步长越小，前进越慢。

梯度下降法算法如下(吐槽：csdn这个文档格式真的是太丑了而且不方便，所以公式还是直接截word图片更清楚)：

以Rosenbrock函数为例，即

于是可得函数的梯度为：

 采用Goldstein原则确定最优步长，同时也有Wolfe法线性搜索，本文采用的是Goldstein原则。

import random
import numpy as np
import matplotlib.pyplot as plt


"""
最速下降法
Rosenbrock函数
函数 f(x)=100*(x(2)-x(1).^2).^2+(1-x(1)).^2
梯度 g(x)=(-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1)),200*(x(2)-x(1)^2))^(T)
"""


def goldsteinsearch(f,df,d,x,alpham,rho,t):
    '''
    线性搜索子函数
    数f，导数df，当前迭代点x和当前搜索方向d
    '''
    flag = 0

    a = 0
    b = alpham
    fk = f(x)
    gk = df(x)

    phi0 = fk
    dphi0 = np.dot(gk, d)
    # print(dphi0)
    alpha=b*random.uniform(0,1)

    while(flag==0):
        newfk = f(x + alpha * d)
        phi = newfk
        # print(phi,phi0,rho,alpha ,dphi0)
        if (phi - phi0 )<= (rho * alpha * dphi0):
            if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
                flag = 1
            else:
                a = alpha
                b = b
                if (b < alpham):
                    alpha = (a + b) / 2
                else:
                    alpha = t * alpha
        else:
            a = a
            b = alpha
            alpha = (a + b) / 2
    return alpha


def rosenbrock(x):
    # 函数:f(x) = 100 * (x(2) - x(1). ^ 2). ^ 2 + (1 - x(1)). ^ 2
    return 100*(x[1]-x[0]**2)**2+(1-x[0])**2


def jacobian(x):
    # 梯度g(x) = (-400 * (x(2) - x(1) ^ 2) * x(1) - 2 * (1 - x(1)), 200 * (x(2) - x(1) ^ 2)) ^ (T)
    return np.array([-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0]),200*(x[1]-x[0]**2)])



def steepest(x0):

    print('初始点为:')
    print(x0,'\n')
    imax = 20000
    W = np.zeros((2, imax))
    epo=np.zeros((2, imax))
    W[:, 0] = x0
    i = 1
    x = x0
    grad = jacobian(x)
    delta = sum(grad ** 2)  # 初始误差

    f=open("最速.txt",'w')

    while i < imax and delta > 10 ** (-5):
        p = -jacobian(x)
        x0 = x
        alpha = goldsteinsearch(rosenbrock, jacobian, p, x, 1, 0.1, 2)
        x = x + alpha * p
        W[:, i] = x
        if i % 5 == 0:

            epo[:,i] =np.array((i,delta))
            f.write(str(i)+"        "+str(delta)+"\n")
            print(i,np.array((i,delta)))
        grad = jacobian(x)
        delta = sum(grad ** 2)
        i = i + 1

    print("迭代次数为:", i)
    print("近似最优解为:")
    print(x, '\n')
    W = W[:, 0:i]  # 记录迭代点

    return [W,epo]

if __name__=="__main__":
    X1 = np.arange(-1.5, 1.5 + 0.05, 0.05)
    X2 = np.arange(-3.5, 4 + 0.05, 0.05)
    [x1, x2] = np.meshgrid(X1, X2)
    f = 100 *(x2 - x1 ** 2) ** 2 + (1 - x1) ** 2  # 给定的函数
    plt.contour(x1, x2, f, 20)  # 画出函数的20条轮廓线

    x0 = np.array([-1.2, 1])
    list_out = steepest(x0)
    W=list_out[0]

    epo=list_out[1]

    plt.plot(W[0, :], W[1, :], 'g*-')  # 画出迭代点收敛的轨迹

    plt.show()

 输出结果（由于在线性搜索子程序中使用了随机函数，初始搜索点是随机产生的，因此每次运行的结果不太相同）：

                                                               图一

                                               图二

迭代误差曲线中蓝色曲线为最速下降法的，橙色线为共轭梯度法的，在下一章会介绍。

最速下降法的优化思想是用当前位置负梯度方向作为搜索方向，因为该方向为当前位置的最快下降方向。从最速下降法迭代路径图中可以看到梯度下降法在接近最优解的区域收敛速度明显变慢，因此利用梯度下降法求解需要很多次的迭代。越接近目标值，步长越小，前进越慢。

注：内容原创，部分文字来源于网络。