跳转至

梯度下降法及其 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\),否则会使这个不等式不成立。

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

image-20221007215257754

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()

unnamed-chunk-5-1

绘制迭代轨迹与等高线图

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=.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$')

迭代轨迹与等高线图

评论