机器学习算法最终总是能转化为最优化问题,习惯上会转化为最小化问题。
个人总结迭代优化算法关键就两点:
(1) 找到下降方向
(2) 确定下降步长

梯度下降算法(Gradient Descent)

梯度下降算法是以最优化函数的负梯度方向为下降方向,更新公式如下: x k + 1 = x k − η g k x_{k+1}=x_k-\eta g_k xk+1=xkηgk其中 g k g_k gk为梯度,学习率 η \eta η用来调整步长。学习率可以设置为随着步长的增加而衰减: η = η 0 1 + i t e r ∗ d e c a y \eta = \frac{\eta_0}{1+iter*decay} η=1+iterdecayη0其中, i t e r 是 当 前 迭 代 步 数 , d e c a y 是 衰 减 率 。 iter是当前迭代步数,decay是衰减率。 iterdecay

最速梯度下降

最速梯度下降算法采用一维精确搜索的反式确定下降步长。即通过极小化 ϕ ( λ ) = L ( x k − λ g k ) \phi (\lambda)=L(x_k - \lambda g_k) ϕ(λ)=L(xkλgk)获得最佳的 λ \lambda λ,即 x k + 1 = x k − λ g k x_{k+1}=x_k-\lambda g_k xk+1=xkλgk最速梯度下降有一个最大的缺点就是 锯齿现象 ,因为相邻两次下降方向是正交的。如下图所示, d 0 垂 直 于 d 1 d_0垂直于d_1 d0d1,造成学习速度下降。 ϕ ′ ( λ ) = ∇ L ( x k − λ g k ) ∗ ( − g k ) = − ∇ L ( x k + 1 ) ∗ g k = g k + 1 g k = 0 \phi^{'}(\lambda)=\nabla L(x_k-\lambda g_k)*(-g_k)=-\nabla L(x_{k+1})*g_k=g_{k+1}g_k=0 ϕ(λ)=L(xkλgk)(gk)=L(xk+1)gk=gk+1gk=0

g

批量梯度下降(BGD)

一次迭代过程,所有训练数据均参与训练。误差的计算是所有样本误差的均值,并以该误差的均值参与梯度的计算。

随机梯度下降(SGD)

一次迭代过程,仅仅一个样本参与训练,且该样本是随机挑选的。

小批量梯度下降(MSGD)

BGD与SGD的结合。一次迭代,随机挑选出K个样本组成小样本集合参与训练,后续训练过程与BGD相同。小批量实际应用最为广泛,原因有:

  1. 虽然批量梯度下降算法计算梯度更精确,但回报是小于线性的。
  2. 小批量适合多核运算。
  3. 小批量有正则化的效果。

一般来讲,仅仅基于梯度的优化算法可以设置较小的K,但需要计算海塞矩阵的优化算法需要较大的K,因为矩阵可能是病态的。

动量算法(momentum)

前面讲到最速梯度下降算法的"锯齿现象"影响了收敛速度,即便不采用一维精确搜索,下降方向也是震荡得厉害,即随机梯度方差较大。因此引入动量项调整梯度方向,减少震荡,加入动量项后,下降的方向不再是梯度方向,而是由当前损失函数梯度与上一轮梯度的合成方向,即 g t = ∇ θ i L ( θ i ) v t = α v t − 1 − η g t θ t = θ t − 1 + v t g_t =\nabla _{\theta_i} L(\theta_i) \\ v_t=\alpha v_{t-1} - \eta g_t \\ \theta_t = \theta_{t-1} + v_t gt=θiL(θi)vt=αvt1ηgtθt=θt1+vt 可见 α \alpha α越大,上一轮梯度的影响越大。之所以取名动量,是类比物理学中单位质量的动量。另外可以采用Nesterov 动量,区别仅仅在对 g t g_t gt进行一个矫正,即梯度计算在施加当前速度之后: θ ^ i = θ t − 1 + α v t − 1 g t = ∇ θ ^ i L ( θ ^ i ) \hat \theta_i = \theta_{t-1} + \alpha v_{t-1} \\ g_t =\nabla _{\hat \theta_i} L(\hat \theta_i) θ^i=θt1+αvt1gt=θ^iL(θ^i) 在凸批量梯度的情况下,Nesterov 动量可将额外误差收敛率从 O ( 1 / k ) \mathcal{O}(1/k) O(1/k)改进到 O ( 1 / k 2 ) \mathcal{O}(1/k^2) O(1/k2),如 Nesterov (1983) 所示。但在随机梯度的情况下,收敛率并未得到改进。

自适应学习率算法(AdaGrad)

AdaGrad算法每个参数拥有独立的,自适应的学习率。用梯度分量平方的累积量来调节学习率以某个参数分量为例: g t = ∇ θ i L ( θ i ) r t = r t − 1 + g t 2 θ t = θ t − 1 − η g t r t + ϵ g_t =\nabla _{\theta_i} L(\theta_i) \\ r_t = r_{t-1} + g_t^2 \\ \theta_t = \theta_{t-1} - \frac{\eta g_t}{\sqrt r_t + \epsilon} gt=θiL(θi)rt=rt1+gt2θt=θt1r t+ϵηgt 其中, η \eta η是初始学习率, ϵ \epsilon ϵ意在防止 r t r_t rt为0,通常设置非常小(如1e-10)。
随着迭代次数的逐步累积, r r r越来越大,学习率越来越小,符合对学习过程先快后慢的要求,但缺点是 r r r会过量过早。

均方根反向传播算法(RMSprop)

指数移动平均(EMA)

一个序列 { a 0 , a 1 , a 3 , . . . , a t , . . . } \{a_0, a_1,a_3,...,a_t,...\} {a0,a1,a3,...,at,...},则这个序列在t时刻的指数移动平均值为: g t = ρ g t − 1 + ( 1 − ρ ) a t g_t = \rho g_{t-1} + (1-\rho)a_t gt=ρgt1+(1ρ)at这是一个迭代式,求和得: g t = ∑ i = 1 t ρ t − i ∗ ( 1 − ρ ) a i g_t=\sum_{i=1}^t\rho ^{t-i}*(1-\rho)a_i gt=i=1tρti(1ρ)ai也就是说t时刻的指数平均值仅仅与t时刻前的序列相关。

参数更新

RMSprop应用梯度分量平方的指数移动平均值来取代AdaGrad算法中梯度分量平方。以某个参数分量为例: g t = ∇ θ i L ( θ i ) m t = ρ m t − 1 + ( 1 − ρ ) g t 2 θ t = θ t − 1 − η g t m t + ϵ g_t =\nabla _{\theta_i} L(\theta_i) \\ m_t = \rho m_{t-1} + (1-\rho)g_t^2 \\ \theta_t=\theta_{t-1}-\frac{\eta g_t }{\sqrt{m_t}+\epsilon}\\ gt=θiL(θi)mt=ρmt1+(1ρ)gt2θt=θt1mt +ϵηgt 其中,超参数 ρ \rho ρ是衰减系数。

自适应矩估计优化算法(Adam)

Adam算法可以看作是动量算法与RMSprop算法的结合,RMS算法仅仅应用了梯度分量的平方来调节参数更新量,而Adam应用梯度分量及其平方共同调节参数更新量,以某个参数分量为例: g t = ∇ θ i L ( θ i ) m t = β 1 m t − 1 + ( 1 − β 1 ) g t v t = β 2 v t − 1 + ( 1 − β 2 ) g t 2 b i a s   c o r r e c t : m t ^ = m t 1 − β 1 t , v t ^ = v t 1 − β 2 t θ t = θ t − 1 − η m t ^ v ^ t + ϵ g_t =\nabla _{\theta_i} L(\theta_i) \\ m_t = \beta_1 m_{t-1} + (1-\beta_1)g_t \\ v_t = \beta_2 v_{t-1} + (1-\beta_2)g_t^2 \\ bias\ correct: \hat {m_t} = \frac{m_t}{1-\beta_1^t},\hat {v_t}=\frac{v_t}{1-\beta_2^t} \\ \theta_t=\theta_{t-1}-\frac{\eta \hat{m_t}}{\sqrt{\hat v_t}+\epsilon} \\ gt=θiL(θi)mt=β1mt1+(1β1)gtvt=β2vt1+(1β2)gt2bias correct:mt^=1β1tmt,vt^=1β2tvtθt=θt1v^t +ϵηmt^ 可见与RMSprop相比,Adam对最后一步参数更新公式中 g t g_t gt进行了修正,即加入了动量项 m t − 1 m_{t-1} mt1。另外, m , v m,v m,v通常初始化为0,而 β 1 , β 2 \beta1, \beta2 β1,β2通常设置较大,例如0.9, 0.999,那么在起始阶段偏执将非常小,因此引入 1 1 − β t \frac{1}{1-\beta^t} 1βt1,随着步数t增大,该值逐渐减小。

牛顿法

梯度下降法仅仅利用了目标函数的一阶导数信息,这意味着未考虑目标函数曲面的曲率,牛顿法考虑了目标函数的二阶导数,迭代更新公式如下: x k + 1 = x k − H k − 1 g k x_{k+1}=x_k-H_k^{-1}g_k xk+1=xkHk1gk其中, H k − 1 H_k^{-1} Hk1为海塞矩阵的逆矩阵, g k g_k gk为梯度。
牛顿法本质是使用海塞矩阵获取沿负梯度方向下降的最优步长,例如在曲率较大时,下降越快,当率较小时则使步长变小。牛顿法只适用于海塞矩阵是正定的情况,在鞍点时就需要正则化,保证朝正确的方向下降。海塞矩阵计算量大,因此深度学习领域很少用到牛顿法。

拟牛顿法

牛顿法虽然收敛速度很快,但海塞矩阵的计算量大。故采用一个正定矩阵 G k G_k Gk近似 H k − 1 H_k^{-1} Hk1 B k B_k Bk近似 H k H_k Hk,称之为拟牛顿算法。BFGS是应用最为广泛的拟牛顿算法,该算法用 B k B_k Bk近似 H k H_k Hk B k B_k Bk更新公式如下: B k + 1 = B k + y k y k T y k T δ k − B k δ k δ k T B k δ k T B k δ k B_{k+1}=B_k+\frac{y_ky_k^T}{y_k^T\delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k} Bk+1=Bk+ykTδkykykTδkTBkδkBkδkδkTBk其中, δ k \delta_k δk为当前自变量增量。
BFGS算法通常需要一维搜索合适的步长。

一维不精确搜索

搜索合适的步长可以有效加快算法收敛速度。一维搜索可分为一维精确搜索和不精确搜索,精确搜索就是每一次迭代,都选取特定步长使目标函数一定能沿着指定方向达到极小(不是最小)。若选取的步长仅使得目标函数下降程度达到用户要求即可,就是不精确搜索。实际中很难获取准确的步长,也没必要耗费计算资源,所以不精确搜索更为实用。一维不精确搜索准则通常要保证两点:1. 目标函数值有足够下降 2.自变量有足够变动幅度
有如下几种一维搜索步长的准则:

  1. Armijo 准则:就是保证目标函数值有足够下降。 f ( X k + α d k ) − f ( X k ) ⩽ ρ α g k T d k α = β γ m β > 0 , γ ∈ ( 0 , 1 ) , ρ ∈ ( 0 , 1 2 ) , m 为 非 负 整 数 f(X_k + \alpha d_k) - f(X_k) \leqslant \rho \alpha g_k^Td_k\\\alpha=\beta\gamma^m\\\beta>0, \gamma \in (0,1), \rho\in(0,\frac{1}{2}),m为非负整数 f(Xk+αdk)f(Xk)ραgkTdkα=βγmβ>0,γ(0,1),ρ(0,21),m
  2. Goldstein准则:通常与Armijo联合使用称为Armijo-Goldstein准则。就是在Armijo准则的基础上,再添加一条准则保证自变量有足够变动幅度。 f ( X k + α d k ) − f ( X k ) ⩾ ( 1 − ρ ) α g k T d k ρ ∈ ( 0 , 1 2 ) f(X_k + \alpha d_k) - f(X_k) \geqslant (1-\rho) \alpha g_k^T d_k\\\rho\in(0,\frac{1}{2}) f(Xk+αdk)f(Xk)(1ρ)αgkTdkρ(0,21)
  3. Wolfe-Powell准则:与Armijo准则联用,采用更简单的式子替换Goldstein准则。 g k + 1 T d k ⩾ σ g k T d k σ ∈ ( ρ , 1 ) g_{k+1}^Td_k\geqslant \sigma g_k^Td_k\\\sigma \in(\rho, 1) gk+1TdkσgkTdkσ(ρ,1)

代码

以最常用的平方误差函数为例,实现各个优化算法,未经严格测试,仅做参考学习。

"""
优化算法的核心两点:1.寻找下降方向  2.一维搜索下降的步长
梯度下降算法以负梯度方向为下降方向
牛顿法以海塞矩阵逆矩阵与梯度乘积为下降方向(此时学习率恒为1)
拟牛顿法同牛顿法,仅仅以B代替海塞矩阵

一维不精确搜索:Wolfe和Armijo确定搜索步长(通常应用于拟牛顿法)
Armijo: f(Xk + eta * Dk) - f(Xk) <= rho * eta * Grad.T * Dk, eta=beta*gamma**m(m为非负整数,beta>0)
Goldstein(需要联合Armijo): f(Xk + eta * Dk) - f(Xk) >= (1-rho) * eta * Grad.T * Dk
Wolfe-Powell(需要联合Armijo):  Grad(k+1).T * Dk >= sigma * Grad(k).T * Dk, rho<sigma<1

以线性回归为例:
平方误差损失函数的最小化的优化问题
J(X)=1/2m(f(X)-y)**2, f(X)=wX,m为样本总数, 最优化函数J(X)对w的梯度为X(wX-y)

输入X:n*d, y:n*1
初始化W:d*1
"""

import numpy as np
from numpy.linalg import norm, inv
from itertools import cycle


def grad(X, y, W):
    return X.T @ (X @ W - y) / X.shape[0]


def initial(X, y, batch_size):
    # 小批量梯度下降初始化
    assert batch_size > 1
    W = np.ones((X.shape[1], 1))
    indx = np.random.permutation(np.arange(X.shape[0]))
    X, y = X[indx], y[indx]
    indx_cycle = cycle(np.arange(X.shape[0] // batch_size + 1) * batch_size)
    return W, indx_cycle


def bgd(X, y, epsilon=1e-2, max_iter=1000, eta=0.1):
    # 批量梯度下降算法
    i = 0
    W = np.ones((X.shape[1], 1))  # 初始化权值 d*1
    while i < max_iter:
        Grad = grad(X, y, W)  # 1*d 计算梯度
        if norm(Grad) < epsilon:  # 如果梯度的第二范数小于阈值,则终止训练
            break
        W -= eta * Grad  # 负梯度方向下降,更新自变量,即权值w
        i += 1

    return W


def sgd(X, y, epsilon=1e-2, max_iter=3600, eta=0.1):
    # 随机梯度下降
    i = 0
    W = np.ones((X.shape[1], 1))
    while i < max_iter:
        indx = np.random.randint(X.shape[0], size=1)  # 随机选择一个样本
        X_s, y_s = X[(indx)], y[(indx)]
        Grad = grad(X_s, y_s, W)
        if norm(Grad) < epsilon:
            break
        W -= eta * Grad
        i += 1
    return W


def msgd(X, y, epsilon=1e-2, max_iter=1000, eta=0.1, batch_size=4):
    # 小批量梯度下降
    W, indx_cycle = initial(X, y, batch_size)
    i = 0
    while i < max_iter:
        j = next(indx_cycle)
        X_s, y_s = X[j:j + batch_size], y[j:j + batch_size]
        Grad = grad(X_s, y_s, W)
        if norm(Grad) < epsilon:
            break
        W -= eta * Grad
        i += 1
    return W


def momentum(X, y, epsilon=1e-2, max_iter=3600, eta=0.1, batch_size=4, alpha=0.01, nesterov=False):
    # 随机梯度下降
    W, indx_cycle = initial(X, y, batch_size)
    i = 0
    v = np.zeros((X.shape[1], 1))  # 初始化动量
    while i < max_iter:
        j = next(indx_cycle)
        X_s, y_s = X[j:j + batch_size], y[j:j + batch_size]
        if nesterov:
            W += alpha * v
        Grad = grad(X_s, y_s, W)
        if norm(Grad) < epsilon:
            break
        v = alpha * v - eta * Grad
        W += v
        i += 1
    return W


def adag(X, y, epsilon=1e-2, max_iter=1000, eta=0.1, batch_size=4, eps_station=1e-10):
    # adaptive gradient descent自适应学习率梯度下降
    W, indx_cycle = initial(X, y, batch_size)
    r = np.zeros((X.shape[1], 1))
    i = 0
    while i < max_iter:
        j = next(indx_cycle)
        X_s, y_s = X[j:j + batch_size], y[j:j + batch_size]
        Grad = grad(X_s, y_s, W)
        if norm(Grad) < epsilon:
            break
        r += Grad ** 2
        W -= eta * Grad / (np.sqrt(r) + eps_station)
        i += 1
    return W


def rms_prop(X, y, epsilon=1e-2, max_iter=1000, eta=0.01, rho=0.9, batch_size=4, eps_station=1e-10):
    # 均方根反向传播算法
    W, indx_cycle = initial(X, y, batch_size)
    m = np.zeros((X.shape[1], 1))
    i = 0
    while i < max_iter:
        j = next(indx_cycle)
        X_s, y_s = X[j:j + batch_size], y[j:j + batch_size]
        Grad = grad(X_s, y_s, W)
        if norm(Grad) < epsilon:
            break
        m = rho * m + (1 - rho) * Grad ** 2
        W -= eta * Grad / (np.sqrt(m) + eps_station)
        i += 1
    return W


def adam(X, y, epsilon=1e-2, max_iter=1000, eta=0.01, beta1=0.9, beta2=0.999, batch_size=4, eps_station=1e-10):
    # adam算法
    W, indx_cycle = initial(X, y, batch_size)
    m = np.zeros((X.shape[1], 1))
    v = np.zeros((X.shape[1], 1))
    i = 0
    while i < max_iter:
        j = next(indx_cycle)
        X_s, y_s = X[j:j + batch_size], y[j:j + batch_size]
        Grad = grad(X_s, y_s, W)
        if norm(Grad) < epsilon:
            break
        m = beta1 * m + (1 - beta1) * Grad
        v = beta2 * v + (1 - beta2) * Grad ** 2
        m_bar = m / (1 - beta1 ** (i + 1))
        v_bar = v / (1 - beta2 ** (i + 1))
        W -= eta * m_bar / (np.sqrt(v_bar) + eps_station)
        i += 1
    return W


def newton(X, y, epsilon=1e-2, max_iter=1000):
    # 牛顿法
    i = 0
    W = np.ones((X.shape[1], 1))
    x1 = X[:, 0]  # 变量维度1,n维向量
    x2 = X[:, 1]
    while i < max_iter:
        err = X @ W - y  # n*1
        Grad = X.T @ err / X.shape[0]
        if norm(Grad) < epsilon:  # 如果梯度的第二范数小于阈值,则终止训练
            break
        err = err.reshape(-1)
        H12 = 2 * x1 @ x2
        H11 = 2 * err @ x1
        H22 = 2 * err @ x2
        H = np.array([[H11, H12], [H12, H22]])  # 计算海塞矩阵
        W -= inv(H) @ Grad  # 负梯度方向下降,更新自变量,即权值w
        i += 1
    return W


def bfgs(X, y, epsilon=1e-4, max_iter=1000):
    # 拟牛顿算法
    i = 0
    W = np.ones((X.shape[1], 1))
    N = X.shape[0]
    B = np.eye(X.shape[1])  # 初始化B d*d
    while i < max_iter:
        err = X @ W - y
        fx = err.T @ err / 2 / N
        Grad = X.T @ err / N  # d*1 计算梯度
        if norm(Grad) < epsilon:  # 如果梯度的第二范数小于阈值,则终止训练
            break
        Dk = - inv(B) @ Grad  # d*1, 下降方向
        eta, W = wp_search(W, Dk, fx, Grad, X, y, N)  # 下降步长以及更新w, 注意弱用WP规则,还可以返回新的梯度,减少重复计算
        delta = eta * Dk  # d*1 自变量w的增量
        yk = B @ delta  # d*1, 更新yk
        B = B + yk @ yk.T / (yk.T @ delta)[0] - B @ delta @ (delta.T @ B) / (delta.T @ B @ delta)[0]  # 更新B
    return W


def wp_search(W, Dk, fx, Grad, X, y, N, sigma=0.75, gamma=0.5, rho=1e-4, beta=1, maxm=100):
    # 基于Wolfe-Powell条件的不精确一维搜索
    assert ((rho < 1.0 / 2) and (rho > 0))
    assert ((gamma < 1.0) and (gamma > 0.0))
    assert ((sigma > rho) and (sigma < 1))
    assert (beta > 0)
    m = 0
    W_new = None
    eta = None
    while m < maxm:
        eta = beta * gamma ** m  # 一维搜索合适的m,进而更新eta
        W_new = W + eta * Dk
        err_new = X @ W_new - y
        fx_new = err_new.T @ err_new / 2 / N  # 下降后的函数值
        diff_val = fx_new - fx  # 下降量
        gdk = Grad.T @ Dk
        exp_diff = eta * gdk
        Grad_new = X.T @ err_new / N  # 更新w后的梯度
        if (diff_val[0] <= rho * exp_diff[0]) and (Grad_new.T @ Dk >= sigma * gdk):
            break
        m += 1
    return eta, W_new


def ag_search(W, Dk, fx, Grad, X, y, N, gamma=0.5, rho=1e-4, beta=1, maxm=100):
    # 基于Armijo-Goldstein条件的不精确一维搜索
    assert ((rho < 1.0 / 2) and (rho > 0))
    assert ((gamma < 1.0) and (gamma > 0.0))
    assert (beta > 0)
    m = 0
    eta = None
    W_new = None
    while m < maxm:
        eta = beta * gamma ** m  # 一维搜索合适的m,进而更新eta
        W_new = W + eta * Dk
        err_new = X @ W_new - y
        fx_new = err_new.T @ err_new / 2 / N  # 下降后的函数值
        diff_val = fx_new - fx  # 下降量
        exp_diff = eta * Grad.T @ Dk
        if (diff_val[0] <= rho * exp_diff[0]) and (diff_val[0] >= (1 - rho) * exp_diff[0]):
            break
        m += 1
    return eta, W_new


def predict(X, W):
    return X @ W


if __name__ == '__main__':
    from sklearn.metrics import mean_squared_error
    from pprint import pprint

    train_data = np.array([[1.1, 1.5, 2.5],
                           [1.3, 1.9, 3.2],
                           [1.5, 2.3, 3.9],
                           [1.7, 2.7, 4.6],
                           [1.9, 3.1, 5.3],
                           [2.1, 3.5, 6.0],
                           [2.3, 3.9, 6.7],
                           [2.5, 4.3, 7.4],
                           [2.7, 4.7, 8.1],
                           [2.9, 5.1, 8.8]])
    X_train, y_train = train_data[:, :-1], train_data[:, [-1]]
    test_data = np.array([[3.1, 5.5, 9.5],
                          [3.3, 5.9, 10.2],
                          [3.5, 6.3, 10.9],
                          [3.7, 6.7, 11.6],
                          [3.9, 7.1, 12.3]])
    X_test, y_test = test_data[:, :-1], test_data[:, [-1]]

    W = bfgs(X_train, y_train, epsilon=1e-5, max_iter=3000)
    y_pred = predict(X_test, W)
    pprint(W)
    print(y_pred, '\n', mean_squared_error(y_test, y_pred))

参考资料

  1. https://arxiv.org/pdf/1412.6980.pdf
  2. http://wemedia.ifeng.com/69799959/wemedia.shtml

我的GitHub
注:如有不当之处,请指正。

Logo

CSDN联合极客时间,共同打造面向开发者的精品内容学习社区,助力成长!

更多推荐