什么是梯度下降法

梯度下降法（英语：Gradient descent）是一个一阶最优化算法，通常也称为最速下降法。 要使用梯度下降法找到一个函数的局部极小值，必须向函数上当前点对应梯度（或者是近似梯度）的反方向的规定步长距离点进行迭代搜索。

一元梯度下降法

为了便于理解，笔者先从一元梯度下降算法开始讲起，如图所示的损失函数为一元二次方程，为了找到函数极小值，我们从随机初始化的点向梯度反方向不断搜索，最终找到近似的极小值点。

代码实现

import numpy as np
import matplotlib.pyplot as plt

# 原函数
def J(x):
    return x**2 - x + 4

# 导函数
def dJ(x):
    return 2*x - 1

# 梯度搜索
def gradient_descent(theta, eta, epsilon=1e-8):
    # 保存训练轨迹
    theta_history = []
    while True:
        gradient = dJ(theta)
        last_theta = theta
        theta = last_theta - eta * gradient
        theta_history.append(theta)
        if abs(J(last_theta) - J(theta)) < epsilon:
            break

    return np.array(theta_history)

if __name__ == '__main__':

    x = np.linspace(-10, 10, 1000)
    y = J(x)

    history = gradient_descent(-10, 0.1)
    print(history)

    # 绘制曲线和梯度下降轨迹
    plt.plot(x, y)
    plt.scatter(history, J(history), color='red')
    plt.show()

多元梯度下降法

以上实现的一元梯度下降法只是为了方便学习和理解，接下来笔者将实现多元梯度下降法，并用一个实例直观的感受下梯度下降法的强大之处。

如图，左图是方程z=2x+3y+4的函数图像，我们将采用梯度下降法从已知数据中拟合出这个二元一次方程。右图即是拟合的过程及拟合结果。

代码实现

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 损失函数
def J(theta, x, y):
    return np.sum((y - np.dot(x, theta))**2) / len(x)

# 导函数
def dJ(theta, x, y):
    res = np.empty(len(theta))
    res[0] = np.sum(np.dot(x, theta) - y)
    for i in range(1, len(theta)):
        res[i] = np.sum(np.dot(np.dot(x, theta) - y, x[:, i]))
    return res * 2 / len(x)

# 梯度搜索
def gradient_descent(theta, x, y, eta, epsilon=1e-8):
    theta_history = []
    while True:
        gradient = dJ(theta, x, y)
        last_theta = theta
        theta = theta - gradient * eta
        esp = abs(J(last_theta, x, y) - J(theta, x, y))
        if esp < epsilon:
            break
        if esp > 0.5:
            theta_history.append(theta)

    theta_history.append(theta)

    return np.array(theta_history)

if __name__ == '__main__':
    # 准备数据
    size = 100
    np.random.seed(999)
    x = np.linspace(-1, 1, size)
    y = np.linspace(-1, 1, size)
    x_b = np.hstack([np.ones((len(x), 1)), x.reshape(-1, 1), y.reshape(-1, 1)])
    z = 2*x+3*y+4+np.random.uniform(-0.2, 0.2, size)

    # 通过梯度下降搜索theta
    eta = 0.05
    init_theta = np.zeros(x_b.shape[1])
    theta_history = gradient_descent(init_theta, x_b, z, eta)

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

    fig = plt.figure()
    # 绘制原始图像 z=2*x+3*y+4
    ax1 = fig.add_subplot(1, 2, 1, projection='3d')
    Z = 2*X + 3*Y + 4 # 原函数
    ax1.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')
    plt.title('z=2*x+3*y+4')

    # 绘制预测图像
    ax2 = fig.add_subplot(1, 2, 2, projection='3d')
    for his in theta_history:
        Z = his[1] * X + his[2] * Y + his[0]  # 预测函数
        ax2.plot_surface(X, Y, Z, rstride=1, cstride=1)

    plt.show()

⬆️