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

机器学习算法最终总是能转化为最优化问题，习惯上会转化为最小化问题。个人总结迭代优化算法关键就两点：(1) 找到下降方向(2) 确定下降步长最速梯度下降算法梯度下降算法是以最优化函数的梯度为下降方向，学习率ηη\eta乘以梯度的模即为下降步长。更新公式如下：xk+1=xk−η∗gkxk+1=xk−η∗gkx_{k+1}=x_k-\eta*g_k其中gkgkg_k为梯度。...

slx_share

11742人浏览 · 2018-05-08 20:17:20

slx_share · 2018-05-08 20:17:20 发布

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

梯度下降算法(Gradient Descent)

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

最速梯度下降

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

批量梯度下降(BGD)

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

随机梯度下降(SGD)

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

小批量梯度下降(MSGD)

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

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

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

动量算法(momentum)

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

自适应学习率算法(AdaGrad)

AdaGrad算法每个参数拥有独立的，自适应的学习率。用梯度分量平方的累积量来调节学习率以某个参数分量为例： $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}$ 其中， $\eta$ 是初始学习率， $\epsilon$ 意在防止 $r_t$ 为0，通常设置非常小(如1e-10)。
随着迭代次数的逐步累积， $r$ 越来越大，学习率越来越小，符合对学习过程先快后慢的要求，但缺点是 $r$ 会过量过早。

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

指数移动平均(EMA)

一个序列 ${a_0, a_1,a_3,...,a_t,...\}$ ,则这个序列在t时刻的指数移动平均值为： $g_t = \rho g_{t-1} + (1-\rho)a_t$ 这是一个迭代式，求和得： $g_t=\sum_{i=1}^t\rho ^{t-i}*(1-\rho)a_i$ 也就是说t时刻的指数平均值仅仅与t时刻前的序列相关。

参数更新

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

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

Adam算法可以看作是动量算法与RMSprop算法的结合，RMS算法仅仅应用了梯度分量的平方来调节参数更新量，而Adam应用梯度分量及其平方共同调节参数更新量，以某个参数分量为例： $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} \\$ 可见与RMSprop相比，Adam对最后一步参数更新公式中 $g_t$ 进行了修正，即加入了动量项 $m_{t-1}$ 。另外, $m, v$ 通常初始化为0，而 $\beta1, \beta2$ 通常设置较大，例如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_ky_k^T}{y_k^T\delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k}$ 其中， $\delta_k$ 为当前自变量增量。
BFGS算法通常需要一维搜索合适的步长。

一维不精确搜索

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

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

参考资料

我的GitHub
注：如有不当之处，请指正。

CSDN学习社区

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

更多推荐

嵌入式作业（七）：基于Ardunio的STM32串口通信

嵌入式作业（七）0作业要求1Ardunio 完成STM32的串口通信（1）安装Ardunio IDE（2）stm32串口通信2关于 stduino IDE0作业要求安装 Ardunio IDE 和相关软件支持库，在Ardunio 完成STM32板子的串口通信程序：（1）持续向串口输出“Hello world！”；（2）当接收到“stop!”时，停止输出。网上有一个国人版的MCU集成开发平台， st

CSDN学习社区

JDBC详解

JDBC文章目录JDBC什么是JDBC?JDBC驱动程序:Java使用JDBC访问数据库的步骤:设置classpath:Oracle连接字符串的书写格式:简单的例子:常用数据库的驱动程序及JDBC URL:Oracle数据库:SQL Server数据库MySQL数据库Access数据库PreparedStatement接口:JNDI-数据源（Data Source）与连接池（Connection

CSDN学习社区

“模式识别与机器学习”学习笔记no2.再谈感知机

接**上篇：上篇主要进行了PLA，Pocket算法的理论过程分析和在给定数据集上利用pocket算法对数据集进行分类学习，得到错分数量最少的分类面。上篇中pocket算法的过程已经进行了编程和测试，框架已经建立了起来，这一篇主要上篇中没有提到或涉及不深的几个问题。1.数据集的构造。上篇是直接使用了题目给的向量，这次来根据正态分布来产生数据集。np.random.normal函数可以根据均值和方差生