梯度下降法基础

对于一个一元函数$f(x)$，如果在点$x$附近有微小的变化$\Delta x$，则$f(x)$的变化$f(x+\Delta x)-f(x)$可表示成如下微分形式。
$f(x+\Delta x)-f(x)\backsimeq f^{\prime }(x)\Delta x$
也就是说，在点$x$附近；如果$\Delta x$和$f^{\prime }(x)$的符号相同，那么$f^{\prime }(x)\Delta x$是正数；如果符号相反，那么$f^{\prime }(x)\Delta x$是负数。
例如，取$\Delta x=-\alpha f^{\prime }(x)$（$\alpha$是一个很小的正数），那么$f(x+\Delta x)-f(x)=-\alpha f^{\prime }(x{)}^2$是负数，即$f(x+\Delta x)$的值比$f(x)$小，或者说，$x$沿着$f^{\prime }(x)$的反方向$-f^{\prime }(x)$运动$\Delta x$的距离，其函数值$f(x+\Delta x)$比原来的$f(x)$小。
让$x$沿着其导函数$f^{\prime }(x)$的反方向（$-f^{\prime }(x)$）移动一个微小的距离$-\alpha f^{\prime }(x)$，就能到达$x_{new}=x-\alpha f^{\prime }(x)$。$x_{new}$的函数值$f^{\prime }(x_{new})$肯定小于之前的函数值$f^{\prime }(x)$。随着$x$不断接近最小值点，$f^{\prime }(x)$将不断接近0（因为在函数极值点$x^{\ast }$，$f^{\prime }(x^{\ast })=0$），$x$移动的距离$\Delta x$越来越接近0。
这就是梯度下降法的思想，即从一个初始的$x$出发，不断用以下公式更新$x$的值。
$x=x-\alpha f^{\prime }(x)$
对于当前的$x$，沿着其导数（梯度）的反方向（$-f^{\prime }(x)$）移动，就能使$f(x)$变小。在理想情况下，对于达到最小值$f(x)$的$x$，有$f^{\prime }(x)=0$。此时，如果迭代更新$x$，$x$的值将不会发生变化。
当然，移动的步伐不能太大。这是因为，根据导数的定义，上述近似公式只在点$x$附近适用。如果移动的步伐太大，就有可能跳过最优值所对应的点$x$，导致点$x$的值来回震荡。梯度下降法就是求一个逼近的最优解。为了避免一直迭代下去，可以用下列方法检查是否足够逼近最优解。

• $f^{\prime }(x)$的绝对值足够小。
• 迭代次数已达到预先设定的最大迭代次数。

梯度下降法的示例代码如下。其中，参数df用于计算函数$f(x)$的导数$f^{\prime }(x)$，x是变量的初始值，alpha是学习率，iteration表示迭代次数，epsilon用于检查$f^{\prime }(x)$的值是否接近0。

def gradient_descent(df, x, alpha=0.01, iterations=100, epsilon=1e-8):
    history = [x]
    for i in range(iterations):
        if abs(df(x)) < epsilon:
            print("梯度足够小！")
            break
        x = x - alpha * df(x)
        history.append(x)
    return history

该函数将迭代过程中所有被更新的$x$保存在一个Python的列表对象history中，并返回这个对象。
以函数$f(x)=x^3-3x^2-9x+2$为例，其导函数为$f^{\prime }(x)=3x^2-6x-9$，要想求$x=1$附近函数的极小值，可调用gradient_descent()函数，代码如下。

df = lambda x: 3 * x ** 2 - 6 * x - 9
path = gradient_descent(df, 1., 0.01, 200)
print(path[-1])
###输出###
梯度足够小！
2.999999999256501

执行以上代码，得到了极值点$x=2.999999999256501$。
以下代码可以将迭代过程中$x$所对应的曲线绘制出来，如图所示。

import numpy as np
import matplotlib.pyplot as plt
f = lambda x: np.power(x, 3) - 3 * x ** 2 - 9 * x + 2
x= np.arange(-3, 4, 0.01)
y = f(x)
plt.plot(x, y)

path_x = np.asarray(path)
path_y = f(path_x)
plt.quiver(path_x[:-1], path_y[:-1], path_x[1:]-path_x[:-1], path_y[1:]-path_y[:-1], scale_units='xy', angles='xy', scale=1, color='k')
plt.scatter(path[-1], f(path[-1]))
plt.show()

对于多变量函数，梯度下降法的原理和单变量函数相同，只不过用梯度代替了导数，公式如下。
$f(x+\Delta x)-f(x)\backsimeq \nabla f(x)\Delta x$
以Beale’s函数公式为例，
$f(x,y)=(1.5-x+xy{)}^2+(2.25-x+x{y}^2{)}^2+(2.625-x+x{y}^3{)}^2$
这个函数的全局极小值点是$(3,0.5)$，可以用如下Python代码计算函数值。

f = lambda x, y: (1.5 - x + x * y) ** 2 + (2.25 - x + x * y ** 2) ** 2 + (2.625 - x + x * y ** 3) ** 2

为了绘制这个函数所对应的曲面，可以先取一些均匀分布的坐标值，示例如下。

xmin, xmax, xstep = -4.5, 4.5, .2
ymin, ymax, ystep = -4.5, 4.5, .2
x_list = np.arange(xmin, xmax + xstep, xstep)
y_list = np.arange(ymin, ymax + ystep, ystep)

然后，根据以上代码中的x_list和y_list，用np.meshgrid()函数求出它们的交叉位置的坐标点(x, y)，并计算这些坐标点所对应的函数值，示例如下。

x, y = np.meshgrid(x_list, y_list)
z = f(x, y)

最后，调用plot_surface()函数，绘制曲面，示例如下。

ax.plot_surface(x, y, z, norm=LogNorm(), rstride=1, cstride=1, edgecolor='none', alpha=.8, cmap=plt.cm.jet)

完整代码如下所示。

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

f = lambda x, y: (1.5 - x + x * y) ** 2 + (2.25 - x + x * y ** 2) ** 2 + (2.625 - x + x * y ** 3) ** 2
minima = np.array([3.0, 0.5])
minima_ = minima.reshape(-1, 1)
xmin, xmax, xstep = -4.5, 4.5, .2
ymin, ymax, ystep = -4.5, 4.5, .2
x_list = np.arange(xmin, xmax+xstep, xstep)
y_list = np.arange(ymin, ymax+ystep, ystep)
x, y = np.meshgrid(x_list, y_list)
z = f(x, y)

fig = plt.figure(figsize=(8, 5))
ax = plt.axes(projection='3d', elev=50, azim=-50)
ax.plot_surface(x, y, z, norm=LogNorm(), rstride=1, cstride=1, edgecolor='none', alpha=.8, cmap=plt.cm.jet)
ax.plot(*minima_, f(*minima_), 'r*', markersize=10)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
ax.set_xlim((xmin, xmax))
ax.set_ylim((ymin, ymax))
plt.show()

接着计算$f(x,y)$关于$(x,y)$的偏导数，并使用Matplotlib的quiver()函数在二维坐标平面上绘制这些坐标点的梯度的方向，代码如下。

df_x = lambda x, y: 2 * (1.5 - x + x * y) * (y - 1) + 2 * (2.25 - x + x * y ** 2) * (y ** 2 - 1) + 2 * (2.625 - x + x * y ** 3) * (y ** 3 - 1)
df_y = lambda x, y: 2 * (1.5 - x + x * y) * x + 2 * (2.25 - x + x * y ** 2) * (2 * x * y) + 2 * (2.625 - x + x * y ** 3) * (3 * x * y ** 2)
dz_dx = df_x(x, y)
dz_dy = df_y(x, y)

fig, ax = plt.subplots(figsize=(10, 6))
ax.contour(x, y, z, levels=np.logspace(0, 5, 35), norm=LogNorm(), cmap=plt.cm.jet)
ax.quiver(x, y, x-dz_dx, y-dz_dy, alpha=.5)
ax.plot(*minima_, 'r*', markersize=18)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xlim((xmin, xmax))
ax.set_ylim((ymin, ymax))
plt.show()

为了直接利用之前给出的梯度下降法的代码，可用一个Numpy向量表示x，即将

if abs(df(x)) < epsilon:

修改为：

if np.max(np.abs(df(x))) < epsilon:

将分离的x, y坐标点数组组合起来，示例如下。

x_ = np.vstack((x.reshape(1, -1), y.reshape(1, -1)))

可以定义一个针对这个向量化的坐标点x的梯度函数df，示例如下。在以下代码中，还实现了修改后的向量化梯度下降法。

df = lambda x: np.array([2 * (1.5 - x[0] + x[0]*x[1]) * (x[1] - 1) + 2 * (2.25 - x[0]
                         + x[0]*x[1]**2) * (x[1]**2 - 1) + 2 * (2.625 - x[0]
                        + x[0]*x[1]**3) * (x[1]**3 - 1), 2 * (1.5 - x[0]
                        + x[0]*x[1])*x[0] + 2 * (2.25 - x[0]
                        + x[0]*x[1]**2) * (2 * x[0] * x[1]) + 2 * (2.625 - x[0]
                        + x[0]*x[1]**3) * (3 * x[0] * x[1] ** 2)])

def gradient_descent(df, x, alpha=0.01, iterations=100, epsilon=1e-8):
    history = [x]
    for i in range(iterations):
        if np.max(np.abs(df(x))) < epsilon:
            print("梯度足够小！")
            break
        x = x - alpha * df(x)
        history.append(x)
    return history

执行以下代码，从[3, 4]出发，求这个曲面的极值点。

x0 = np.array([3., 4.])
print("初始点", x0, "的梯度"，df(x0))
path = gradient_descent(df, x0, 0.000005, 300000)
print("极值点: ", path[-1])

因为x的初始梯度值很大，所以学习率$\alpha$必须取很小的值，否则会导致震荡或得到的值很大。最后收敛到点$[2.70735828,0.41689171]$，但它不是最优点。
以下代码将绘制迭代过程中x的变化情况。

def plot_path(path, x, y, z, minima_, xmin, xmax, ymin, ymax):
    fig, ax = plt.subplots(figsize=(10,  6))
    ax.contour(x, y, z, levels=np.logspace(0, 5, 35), norm=LogNorm(), cmap=plt.cm.jet)
    ax.quiver(path[:-1, 0], path[:-1, 1], path[1:, 0]-path[:-1, 0], path[1:, 1]-path[:-1, 1], scale_units='xy', angles='xy', scale=1, color='k')
    ax.plot(*minima_, 'r*', markersize=18)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_xlim((xmin, xmax))
    ax.set_xlim((ymin, ymax))
    plt.show()

path = np.asarray(path)
plot_path(path, x, y, z, minima_, xmin, xmax, ymin, ymax)