梯度下降法及其 Python 实现

用数值近似法求函数在某点的梯度，用回溯线搜索法控制步长，应用梯度下降法求函数极值。

梯度下降的思想是：对某一初始值，不断改变这一初始值，且每一步都朝着能使函数值减小的方向改变，最终函数值几乎不再变小，我们就认为达到了极小值点。

如何计算梯度

显式解

对于偏导函数有显示解的函数，当然可以把偏导数直接写出来，再将坐标点代入偏导函数中，得到梯度向量。

\[ \nabla f(x,y) = \left ( {\partial f \over \partial x}, {\partial f \over \partial y} \right ). \]

例子

求函数在\((x,y,z)=(1,2,\pi/3)\)的梯度

\[ f(x,y,z)=x^2y\sin(z) \]

一阶偏导数为：

\[ \nabla f(x,y,x)=(2xy \sin(z), x^2 \sin(z), x^2 y \cos(z)) = (2 \sqrt{3},\sqrt{3},1) \]

数值方法近似

若不方便求出梯度向量的显式解，可以用数值方法近似求出梯度向量。

\[ \nabla f(x) \approx\left[\begin{array}{c} \frac{f\left(x+\varepsilon e_1\right)-f\left(x-\varepsilon e_1\right)}{2 \varepsilon} \\\ \frac{f\left(x+\varepsilon e_2\right)-f\left(x-\varepsilon e_2\right)}{2 \varepsilon} \\\ \vdots \\\ \frac{f\left(x+\varepsilon e_n\right)-f\left(x-\varepsilon e_n\right)}{2 \varepsilon} \end{array}\right] \]

如何确定步长（学习率）

得到函数在某一点的梯度之后，我们就可以沿着梯度的反方向更新坐标点。

固定步长

如果每一步迭代时移动的步长都相同，可能会出现问题：

• 步长太大，会导致迭代过快，甚至有可能错过最优解。
• 步长太小，迭代速度太慢，很长时间算法都不能结束。

一些算法可以选择自适应的步长，例如回溯线搜索。

自适应步长

回溯线搜索

回溯线搜索是一种自适应步长的方法，它的算法为：

取常数\(\alpha, \beta\)满足\(0<\alpha<0.5,0<\beta<1\)., 及\(\epsilon>0\)

1. 令\(k=0\)给定迭代初始值\(x^0\)及初始步长\(t_0=1\).
2. 更新\(t_k=\beta t_k\)

3. 重复 2 直到

   \[ f\left(x^k-t_k \nabla f\left(x^k\right)\right) \leq f\left(x^k\right)-\alpha t_k\left\|\nabla f\left(x^k\right)\right\|^2 \]

4. 计算\(x^{k+1}=x^k-t^k \nabla f\left(x^k\right)\)并更新\(k=k+1\)

5. 转至 (1) 直到\(\nabla f\left(x^k\right) \leq \epsilon\)

其中的第 (3) 步最关键，几何上理解第 (3) 步：

在\(t=0\)点的切线斜率就是梯度，把这个斜率乘以\(\alpha\)，就得到一条更平缓的虚线，即下图中处在上方的虚线。

现在我们要求

\[ f\left(x^k-t_k \nabla f\left(x^k\right)\right) \leq f\left(x^k\right)-\alpha t_k\left\|\nabla f\left(x^k\right)\right\|^2 \]

也就是更新坐标点后的函数值小于等于上方的虚线的函数值。因此，\(t\)只能处在\((0,t_0)\)之间，不能大于\(t_0\)，否则会使这个不等式不成立。

这样就限制了步长不能太大，也就是不能跨越、偏离极值点太多。

Python 实现

问题描述

考虑\(R^2\)上的二元函数：

\[ f\left(x_1, x_2\right)=e^{x_1+3 x_2-0.1}+e^{x_1-3 x_2-0.1}+e^{-x_1-0.1} \]

利用梯度下降法求其极小值与极小值点。注意： 1. 需利用数值算法计算函数的梯度向量； 2. 需利用回溯线搜索 (Backtracking Line Search) 或精确的线搜索 (Exact line search) 选择自适应的步长。

代码

应用梯度下降法，迭代更新极小值点

Python

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

# 支持中文
matplotlib.rcParams["font.sans-serif"] = ["SimSun"]
matplotlib.rcParams["axes.unicode_minus"] = False
# 渲染公式
plt.rcParams["text.usetex"] = True

定义函数，用于求原函数、偏导函数、梯度和向量模。

Python

def f(X):
    x1 = X[0]
    x2 = X[1]
    item1 = np.exp(x1 + 3 * x2 - 0.1)
    item2 = np.exp(x1 - 3 * x2 - 0.1)
    item3 = np.exp(-x1 - 0.1)
    return item1 + item2 + item3


def derivative_x1(x1, x2, delta=1e-6):
    X_delta1 = [x1 + delta, x2]
    X_delta2 = [x1 - delta, x2]
    return (f(X_delta1) - f(X_delta2)) / (2 * delta)


def derivative_x2(x1, x2, delta=1e-6):
    X_delta1 = [x1, x2 + delta]
    X_delta2 = [x1, x2 - delta]
    return (f(X_delta1) - f(X_delta2)) / (2 * delta)


def gradient(X, delta=1e-6):
    x1 = X[0]
    x2 = X[1]
    return np.array([derivative_x1(x1, x2, delta), derivative_x2(x1, x2, delta)])


def norm(vector):
    return np.sqrt(sum([i**2 for i in vector]))

设置参数：

• 初始值为np.array([1, 1])
• 回溯线搜索参数\(\alpha=0.4, \beta=0.5\)
• 梯度的模的精度为norm(gradient(X)) <\(1e-6\)
• 初始步长\(t=1\)

Python

# 回溯线性搜索实现梯度下降
alpha = 0.4
beta = 0.5
epsilon = 1e-6
k = 0

X = np.array([1, 1])
t = 1

k_list = [k]
X_all = X

应用梯度下降法，迭代更新极小值点。

Python

# 当梯度的模仍然高于精度要求时，需要一直迭代
while norm(gradient(X)) > epsilon:
    # 回溯线性搜索，检查目前的步长是否满足要求，若步长太大，则需要缩小步长
    while f(np.subtract(X, t * gradient(X))) > f(X) - alpha * t * gradient(
        X
    ).T @ gradient(X):
        t = beta * t
    # 记录此次的自变量值
    k = k + 1
    k_list.append(k)
    X = X - t * gradient(X)
    X_all = np.vstack((X_all, X))
    # print("k = %d, X = %s, f(X) = %f" % (k, X, f(X)))
# 退出循环，说明已经得到了满足精度要求的梯度。此时的自变量值就是极小值点
print("极小值点为：X* = {}\n 极小值为：f(X*) = {}".format(X, f(X)))

Text Only

## 极小值点为：X* = [-3.4657320e-01 -1.0000889e-12]
## 极小值为：f(X*) = 2.559266696658411

绘制绘制迭代次数与自变量取值

Python

# 绘制迭代次数与自变量取值
fig = plt.figure(figsize=(5, 4), dpi=200)
plt.plot(X_all[:, 0], "r-o", markerfacecolor="none", label=r"$x_1$")
plt.plot(X_all[:, 1], "b-o", markerfacecolor="none", label=r"$x_2$")
plt.xlabel("迭代次数", usetex=False)  # 中文不用 tex 渲染，否则会报错
plt.ylabel("自变量的取值", usetex=False)
plt.legend()

绘制迭代轨迹与等高线图

Python

# 绘制迭代轨迹与等高线图
x1 = np.linspace(-2, 1.2, 1000)
x2 = np.linspace(-1, 1.2, 1000)
X1, X2 = np.meshgrid(x1, x2)
Y = np.exp(X1 + 3 * X2 - 0.1) + np.exp(X1 - 3 * X2 - 0.1) + np.exp(-X1 - 0.1)
fig = plt.figure(figsize=(5, 4), dpi=200)
plt.contourf(
    X1, X2, Y, levels=[3, 4, 5, 6, 7, 8, 9, 10], alpha=0.75, cmap="coolwarm"
)  # contourf 可以绘制等高线并为区域涂色

Python

C = plt.contour(
    X1, X2, Y, levels=[3, 4, 5, 6, 7, 8, 9, 10], colors="black"
)  # contour 可以绘制等高线，但不涂色
plt.clabel(C, inline=True)  # 在等高线上显示函数值
plt.plot(
    [i[0] for i in X_all], [i[1] for i in X_all], "k-o", markerfacecolor="none"
)  # 绘制迭代轨迹
plt.xlabel(r"$x_1$")
plt.ylabel(r"$x_2$")

评论