机器学习中优化算法总结以及Python实现

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

梯度下降算法(Gradient Descent)

梯度下降算法是以最优化函数的负梯度方向为下降方向，更新公式如下：$$x_{k+1}=x_k-\eta g_k$$其中$g_k$为梯度，学习率$\eta$用来调整步长。学习率可以设置为随着步长的增加而衰减：$$\eta =\frac{{\eta }_0}{1+iter\ast decay}$$其中，$iter是当前迭代步数，decay是衰减率。$

最速梯度下降

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

批量梯度下降(BGD)

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

随机梯度下降(SGD)

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

小批量梯度下降(MSGD)

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

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

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

动量算法(momentum)

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

自适应学习率算法(AdaGrad)

AdaGrad算法每个参数拥有独立的，自适应的学习率。用梯度分量平方的累积量来调节学习率以某个参数分量为例：$$\begin{gathered} 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+ϵ} \end{gathered}$$ 其中，$\eta$是初始学习率，$ϵ$意在防止$r_t$为0，通常设置非常小(如1e-10)。
随着迭代次数的逐步累积，$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=\rho g_{t-1}+(1-\rho )a_t$$这是一个迭代式，求和得：$$g_t=\sum _{i=1}^t{\rho }^{t-i}\ast (1-\rho )a_i$$也就是说t时刻的指数平均值仅仅与t时刻前的序列相关。

参数更新

RMSprop应用梯度分量平方的指数移动平均值来取代AdaGrad算法中梯度分量平方。以某个参数分量为例：$$\begin{gathered} 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}+ϵ} \end{gathered}$$ 其中，超参数$\rho$是衰减系数。

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

Adam算法可以看作是动量算法与RMSprop算法的结合，RMS算法仅仅应用了梯度分量的平方来调节参数更新量，而Adam应用梯度分量及其平方共同调节参数更新量，以某个参数分量为例：$$\begin{gathered} 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 \\ biascorrect:\widehat{m_t}=\frac{m_t}{1-{\beta }_1^t},\widehat{v_t}=\frac{v_t}{1-{\beta }_2^t} \\ {\theta }_t={\theta }_{t-1}-\frac{\eta \widehat{m_t}}{\sqrt{{\hat{v}}_t}+ϵ} \end{gathered}$$ 可见与RMSprop相比，Adam对最后一步参数更新公式中$g_t$进行了修正，即加入了动量项$m_{t-1}$。另外,$m,v$通常初始化为0，而$\beta 1,\beta 2$通常设置较大，例如0.9, 0.999,那么在起始阶段偏执将非常小，因此引入$\frac{1}{1-{\beta }^t}$,随着步数t增大，该值逐渐减小。

牛顿法

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

拟牛顿法

牛顿法虽然收敛速度很快，但海塞矩阵的计算量大。故采用一个正定矩阵$G_k$近似$H_k^{-1}$或$B_k$近似$H_k$，称之为拟牛顿算法。BFGS是应用最为广泛的拟牛顿算法，该算法用$B_k$近似$H_k$。$B_k$更新公式如下：$$B_{k+1}=B_k+\frac{y_k{y}_k^T}{y_k^T{\delta }_k}-\frac{B_k{\delta }_k{\delta }_k^T{B}_k}{{\delta }_k^T{B}_k{\delta }_k}$$其中，${\delta }_k$为当前自变量增量。
BFGS算法通常需要一维搜索合适的步长。

一维不精确搜索

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

1. Armijo 准则：就是保证目标函数值有足够下降。$$\begin{gathered} f(X_k+\alpha d_k)-f(X_k)⩽\rho \alpha g_k^T{d}_k \\ \alpha =\beta {\gamma }^m \\ \beta >0,\gamma \in (0,1),\rho \in (0,\frac{1}{2}),m为非负整数 \end{gathered}$$
2. Goldstein准则：通常与Armijo联合使用称为Armijo-Goldstein准则。就是在Armijo准则的基础上，再添加一条准则保证自变量有足够变动幅度。$$\begin{gathered} f(X_k+\alpha d_k)-f(X_k)⩾(1-\rho )\alpha g_k^T{d}_k \\ \rho \in (0,\frac{1}{2}) \end{gathered}$$
3. Wolfe-Powell准则:与Armijo准则联用，采用更简单的式子替换Goldstein准则。$$\begin{gathered} g_{k+1}^T{d}_k⩾\sigma g_k^T{d}_k \\ \sigma \in (\rho ,1) \end{gathered}$$

代码

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

"""
优化算法的核心两点：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
注：如有不当之处，请指正。