机器学习优化算法（一）
在机器学习领域，数理统计、概率论和优化理论等知识是必不可少的，尤其是以优化理论最为重要。机器学习经常会把问题抽象为一个优化问题，比如最大化回报奖赏、最小化分类错误带来的损失或者最大化似然等等，而如何求解这些问题就必须有相应的优化算法。本文将主要介绍几个常见的优化算法：
[*]随机梯度下降（Stochastic Gradient Descent, SGD）
[*]改进的迭代尺度算法（Improved Iterated Scaling, IIS）
[*]牛顿法（Newton） & 拟牛顿法（Quasi-Newton）
[*]DFP
[*]BFGS & L-BFGS
文章主要参考了李航老师的《统计学习方法》里面相关的介绍。
A. SGD
本文首先介绍SGD，这是最简单的优化方法，也是最为常用的方法。在深度学习里面，一般都会使用SGD或其变种进行优化，也是反向传播算法的核心。
那么对于SGD，假设优化目标为：
\min_w f(w) \\ 根据一阶泰勒展开式（假设目标函数有连续的一阶偏导数），那么有下式，其中 w_k 是第 k 次迭代的的变量值， g_k = \nabla f(w_k) 是函数 f 在 w_k 处的梯度：
f(w) = f(w_k) + g_k^T (w - w_k) \\ 为了使得目标函数的值下降的最快，那么选择 w - w_k = -\lambda g_k 最为合适，即有以下更新：
w_{k+1} = w_k - \lambda g_k \\
使用以上权重更新之后的目标值更新的结果为：
f(w_{k+1}) = f(w_k) + g_k^T (w_{k+1} - w_k) = f(w_k) - \lambda g_k^Tg_k \leq f(w_k) \\
可以看出目标值是下降的。其中 \lambda 是学习率，可以使用固定值，也可以使用线性搜索求解，即：
\lambda^{*} = \arg \min_{\lambda \geq 0} f(w_k - \lambda g_k) \\
一般来说，梯度下降法最为简单实用，但收敛速率有点慢。随机梯度下降则是基于批处理（Mini-Batch）的梯度下降，即函数的梯度信息是通过一个Mini-Batch的数据计算得到的，而不是所有的数据的梯度。
在深度学习领域，SGD是最为常用的优化方法，因为深度网络所拟合的函数是非线性、非凸的，很难求解最优解，与其耗费很多功夫设计一个优化算法，倒不如使用最为简单、通用的SGD。SGD会有很多变种，包括Momentum、AdaGrad、RMSProp、Adam等等。
经常会有这么一个问题：为什么深度学习里面的学习率都是固定的或者固定步长衰减？而不用这种线性搜索寻找最优的学习率呢？这是因为在深度网络里面，主要要求的是速度，正如深度网络不使用二阶优化（二阶优化需要大量的存储和计算时间）的目的一样，使用线性搜索带来的代价比较大。此外，由于深度网络里面都是使用的Mini-Batch的形式，计算的梯度只是真实梯度的一个近似，使用精确的线性搜索的结果不一定是最优的，所以倒不如直接固定学习率了。
如果将待优化的参数记为 w ，其梯度信息记作 g = \nabla_w J(w) ，第 t 次迭代对应的参数和梯度分别为 w_t 和 g_t ，那么有下面的结果：
[*]SGD采用MiniBatch数据计算出来的梯度作梯度下降。
w_{t+1} = w_t - \eta g_t\\
[*]Momentum计算梯度时综合考虑了当前步梯度和上一步梯度的信息，一方面，由于动量 m_{t-1} 的影响，梯度方向不会大幅度震荡，尽量和上一次更新的方向保持一致；另一方面，当陷入平缓区域时， g_t 近乎为0，但此时由于有动量在，可以跳出平缓区域。
m_t = \lambda m_{t-1} + g_t \\ w_{t+1} = w_t - \eta m_t
[*]AdaGrad对于不同的参数使用不同的学习率（根据之前轮数梯度大小的平方和进行调整），后期很可能会由于学习率调整的分母太大导致学习率为零。
w_{t+1, i} = w_{t, i} - \frac{\eta}{\sqrt{G_{t, ii}+ \epsilon}}g_{t, i} \\ G_{t,ii} = \sum_{k=1}^{t-1} g_{k,i}^2
[*]RMSProp是使用梯度大小的平方值的平均值（而不是AdaGrad的求和）来约束不同参数的学习率。
E = \alpha E + (1-\alpha)g_t^2\\ w_{t+1, i} = w_{t, i} - \frac{\eta}{\sqrt{E+ \epsilon}}g_{t, i} \\
[*]Adam是Adaptive Moment estimation的简称，将Momentum和自适应学习率两个加速的办法综合了起来，也是最为常用的优化算法之一。
m_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t \\ n_t = \beta_2 n_t + (1-\beta_2) g_t^2 \\ \hat{m}_t = \frac{m_t}{1-\beta_1^t} \\ \hat{n}_t = \frac{n_t}{1-\beta_2^t} \\ w_{t+1} = w_t - \frac{\eta}{\sqrt{\hat{n}_t}+\epsilon}\hat{m}^t
B. IIS
改进的迭代尺度算法（IIS）的核心思想是寻找一组参数 w + \delta ，计算函数值增加（这里假设优化目标是最大化目标函数）的幅度 L(w + \delta ) - L(w) 。一般来说，会寻找函数值增加幅度的下界 A(\delta | w) ，最大化这个下界，从而可以寻找到一组参数的更新 \delta 。具体优化策略为：
\delta^{*} = \arg \max_{\delta} A(\delta | w) \\ s.t. \quad L(w + \delta) - L(w) \geq A(\delta | w)
下面以最大熵模型为例进行介绍，最大熵模型（Maximum Entropy）对条件分布建模为：
p_w(y|x) = \frac{1}{Z_w(x)}\exp \left( \sum_{i=1}^n w_i f_i(x,y) \right) \\ Z_w(x) = \sum_y \exp \left( \sum_{i=1}^n w_i f_i(x,y) \right)
其中 f_i(x,y) 为第 i 个特征函数对应的特征值。那么最大熵模型一般是在经验分布 \tilde{p}(x, y) 上利用最大似然求解，经验分布可以取 \tilde{p}(x, y)=\frac{1}{N} 的均匀分布。优化目标为：
L(w) = \sum_{x,y} \tilde{p}(x, y) \log p_w(y|x) \\
计算 L(w + \delta) - L(w) 得到：
L(w + \delta) - L(w) = \sum_{x,y} \tilde{p}(x, y) \log p_{w+\delta}(y|x) - \sum_{x,y} \tilde{p}(x, y) \log p_w(y|x) \\ = \sum_{x,y} \tilde{p}(x, y) \sum_{i=1}^n \delta_i f_i(x, y) - \sum_x \tilde{p}(x)\log \frac{Z_{w+\delta}(x)}{Z_w(x)} \\ \geq \sum_{x,y} \tilde{p}(x, y) \sum_{i=1}^n \delta_i f_i(x, y) + 1 - \sum_x \tilde{p}(x)\frac{Z_{w+\delta}(x)}{Z_w(x)}
下面推导两个归一化因子的除法：
\frac{Z_{w+\delta}(x)}{Z_w(x)} = \frac{\sum_y \left( \exp \left( \sum_{i=1}^n w_i f_i(x,y) \right) \exp\left(\sum_{i=1}^n \delta_i f_i(x, y)\right)\right)}{\sum_y \exp \left( \sum_{i=1}^n w_i f_i(x,y) \right)} \\ = \sum_y \left( \frac{\exp \left( \sum_{i=1}^n w_i f_i(x,y) \right)}{\sum_y \exp \left( \sum_{i=1}^n w_i f_i(x,y) \right)}\exp \left( \sum_{i=1}^n \delta_i f_i(x,y) \right) \right) \\ = \sum_y p_w(y|x)\exp \left( \sum_{i=1}^n \delta_i f_i(x,y) \right)
那么有：
L(w + \delta) - L(w) \geq \sum_{x,y} \tilde{p}(x, y) \sum_{i=1}^n \delta_i f_i(x, y) + 1 - \sum_x \tilde{p}(x)\sum_y p_w(y|x)\exp \left( \sum_{i=1}^n \delta_i f_i(x,y) \right)
按道理，做到上面的程度就是找到了函数值提升的一个下界了，此时最大化这个下界求解 \delta 即可，那么就需要求偏导，发现求偏导之后 \delta_i 和 \delta_j 会出现在同一个式子里面（主要是最后面的exp那一项决定的），因此想办法将 \delta_i 独立出来。
引入 f^{\#}(x,y) = \sum_i f_i(x, y) ，那么有：
\exp \left( \sum_{i=1}^n \delta_i f_i(x,y) \right) = \exp \left( \sum_{i=1}^n \frac{f_i(x,y)}{f^{\#}(x,y)}f^{\#}(x,y)\delta_i \right) \leq \sum_{i=1}^n \frac{f_i(x,y)}{f^{\#}(x,y)} \exp\left( f^{\#}(x,y)\delta_i \right)
那么代入上面的下界公式得到 L(w + \delta) - L(w) \geq A(\delta|w) ，有：
A(\delta | w) = \sum_{x,y} \tilde{p}(x, y) \sum_{i=1}^n \delta_i f_i(x, y) + 1 - \sum_x \tilde{p}(x)\sum_y p_w(y|x)\sum_{i=1}^n \frac{f_i(x,y)}{f^{\#}(x,y)} \exp\left( f^{\#}(x,y)\delta_i \right)
对其中的 \delta_i 求偏导有：
\sum_{x,y} \tilde{p}(x, y) f_i(x, y) - \sum_x \tilde{p}(x)\sum_y p_w(y|x) f_i(x,y) \exp\left( f^{\#}(x,y)\delta_i \right) = 0 \\
上述公式中可以看出 \delta_i 的求解和其余的 \delta_j 无关，大大简化了计算，所以可以求出 \delta 。至此，IIS用来求解最大熵模型的公式推导完毕。
C. 牛顿法 & 拟牛顿法
牛顿法相比较于梯度下降法，主要是利用了二阶的梯度信息，即海塞阵（Hessian Matrix）信息。Hessian矩阵的定义为：
H = \left[ \frac{\partial^2 f}{\partial{w_i}\partial{w_j}} \right]_{n \times n} \\
类似梯度下降，使用二阶泰勒展开式得到：
f(w) = f(w_k) + g_k^T(w - w_k) + \frac{1}{2}(w-w_k)^TH_k(w-w_k) \\
对 w 求导得到：
\nabla f(w)= g_k + H_k(w - w_k) \\
为了求得极值，选择 \nabla f(w) = 0 ，即：
w_{k+1} = w_k - H_k^{-1}g_k \\
假设 H_k 是正定的，那么选择了这一步权重更新操作之后的目标值也会下降，得到下降的目标值为：
f(w_{k+1}) = f(w_k) + g_k^T(w_{k+1} - w_k) + \frac{1}{2}(w_{k+1}-w_k)^TH_k(w_{k+1}-w_k) \\ = f(w_k) - \frac{1}{2}g_k^TH_k^{-1} g_k < f(w_k)\\
以上便是牛顿法，使用了二阶的梯度信息，要求函数二阶梯度连续。然而，每次对Hessian阵求逆比较复杂，因此便引入了拟牛顿法。
拟牛顿法主要是为了近似 H_k^{-1} 或者 H_k 而提出的算法，主要包括DFP、BFGS和L-BFGS等。
D. DFP
DFP是Davidon-Fletcher-Powell的缩写，主要思想是对 H_k^{-1} 进行近似。通过牛顿法，我们得到：
\nabla f(w) = \nabla f(w_k) + H_k (w - w_k) \\
取 w = w_{k+1}, g_k = \nabla f(w_k), y_k = g_{k+1} - g_k, \delta_k = w_{k+1}-w_k 得到拟牛顿条件：
H_k^{-1}y_k = \delta_k \\
下面我们要通过 w_k, \delta_k, g_k, y_k 来近似地将 H_k^{-1} 近似出来，这样就不用每次迭代都去求一次Hessian矩阵及其逆矩阵了，我们引入近似 G_{k+1} = H_k^{-1} ，并且假设有以下迭代公式：
G_{k+1} = G_k + P_k + Q_k \\
其中 G_0 是初始化为正定矩阵的矩阵，那么由 H_k^{-1}y_k = \delta_k ，我们得到下面的式子：
G_k y_k + P_k y_k + Q_k y_k = \delta_k \\
为了求出 P_k, Q_k （我们只要求出一组合理、简单的解满足上面的式子即可），我们可以令：
G_k y_k + Q_k y_k = 0 \\ P_k y_k = \delta_k
那么可以通过观察寻找到一组合适的解（读者可以自行验证上面的式子是否成立）：
Q_k = -\frac{G_k y_k y_k^T G_k}{y_k^TG_ky_k} \\ P_k = \frac{\delta_k\delta_k^T}{\delta_k^Ty_k}
由于上面得到矩阵 P_k, Q_k 的秩都为1，所以更新公式也可以写为，这就是秩2-校正公式：
G_{k+1} = G_k + \alpha uu^T + \beta vv^T \\
那么 H_k^{-1} = G_{k+1} 即可通过下面的式子迭代计算得到：
G_{k+1}= G_k -\frac{G_k y_k y_k^T G_k}{y_k^TG_ky_k} + \frac{\delta_k\delta_k^T}{\delta_k^Ty_k} \\
假如初始化得到的矩阵 G_0 是正定对称的，那么 G_k 也是正定对称的。可以用数学归纳法证明，假设 G_k 是对称正定的，那么有：
G_{k+1}^T = G_k^T -\frac{G_k^T y_k y_k^T G_k^T}{y_k^TG_ky_k} + \frac{\delta_k\delta_k^T}{\delta_k^Ty_k} = G_{k+1}\\
对称性得到证明。
假设矩阵 G_k 对称正定，那么 存在G_k = LL^T 成立。根据柯西不等式有：
(a^Tb)^2 \leq (a^Ta)(b^Tb) \\ a = L^Tx, \quad b = L^Ty \quad \Rightarrow \quad (x^TG_k y)^2 \leq (x^TG_k x) (y^TG_k y)
下面的式子证明正定性，对于任意的 x \neq 0 ：
x^TG_{k+1}x = x^TG_k x -\frac{x^TG_k y_k y_k^T G_k x}{y_k^TG_ky_k} + \frac{x^T \delta_k\delta_k^T x}{\delta_k^Ty_k} \\ > x^TG_k x - \frac{x^TG_kx y_k^T G_k y_k}{y_k^TG_k y_k} + 0 = 0
所以正定性也可以得到证明。
那么DFP的算法总结一下就是下面的流程：
[*]初始化矩阵 G_0 ，待优化权重 w_0
[*]for k = 0 to maxiter
[*] 计算梯度 g_k = \nabla f(w_k) ，计算更新方向 d_k = - G_k g_k
[*] 线性搜索最优的步长 \eta_k = \arg \min_{\eta} f(w_k + \eta d_k)
[*] 权重更新 w_{k+1} = w_k + \eta d_k
[*] 计算梯度 g_{k+1} = \nabla f(w_{k+1})
[*] 计算 y_k = g_{k+1} - g_k, \delta_k = w_{k+1}-w_k
[*] 更新 G_{k+1}= G_k -\frac{G_k y_k y_k^T G_k}{y_k^TG_ky_k} + \frac{\delta_k\delta_k^T}{\delta_k^Ty_k}
[*] 若 ||g_{k+1}||_2 \leq \epsilon 退出循环，近似得到最优解 w^{*} = w_{k+1}
[*] 令k = k + 1，继续循环迭代
D. BFGS
BFGS同DFP一样，也是拟牛顿法的一种，全称为Broyden–Fletcher–Goldfarb–Shanno。它的主要思想是使用 G_{k+1} 近似 H_k ，而不是近似 H_k^{-1} 。并且一般来说，BFGS比DFP更为常用，是因为它具有自校正的性质，可以在对Hessian矩阵估计出现偏差之后的几轮内校正，而DFP则没有这个性质。
对比一下DFP的核心公式，使用 G_{k+1} \approx H_k^{-1} ：
H_k^{-1}y_k = \delta_k \\ G_{k+1}= G_k -\frac{G_k y_k y_k^T G_k}{y_k^TG_ky_k} + \frac{\delta_k\delta_k^T}{\delta_k^Ty_k} \\
在BFGS里面是使用 G_{k+1} \approx H_k :
y_k = H_k \delta_k \\ G_{k+1}= G_k -\frac{G_k \delta_k \delta_k^T G_k}{\delta_k^TG_k \delta_k} + \frac{y_k y_k^T}{y_k^T \delta_k} \\
所以BFGS的算法流程为：
[*]初始化矩阵 G_0 ，待优化权重 w_0
[*]for k = 0 to maxiter
[*] 计算梯度 g_k = \nabla f(w_k) ，计算更新方向 d_k ，满足 G_k d_k = -g_k
[*] 线性搜索最优的步长 \eta_k = \arg \min_{\eta} f(w_k + \eta d_k)
[*] 权重更新 w_{k+1} = w_k + \eta d_k
[*] 计算梯度 g_{k+1} = \nabla f(w_{k+1})
[*] 计算 y_k = g_{k+1} - g_k, \delta_k = w_{k+1}-w_k
[*] 更新 G_{k+1}= G_k -\frac{G_k \delta_k \delta_k^T G_k}{\delta_k^TG_k \delta_k} + \frac{y_k y_k^T}{y_k^T \delta_k}
[*] 若 ||g_{k+1}||_2 \leq \epsilon 退出循环，近似得到最优解 w^{*} = w_{k+1}
[*] 令k = k + 1，继续循环迭代
下面介绍L-BFGS，L指的是Limited-memory，是指的在存储空间上对BFGS进行优化的变种。
牛顿法需要求Hessian矩阵的逆矩阵，时间和空间复杂度为 O(n^3) ，拟牛顿法使用中间矩阵迭代近似Hessian矩阵或者其逆矩阵，时空复杂度为 O(n^2) ，而LBFGS则不需显式地构造和存储Hessian矩阵或其逆矩阵的近似，而是通过前 m 个历史梯度和目标值进行拟合，时间和空间复杂度为 O(mn) 。
一般来说，对于参数数量 n = 10^4 以上的优化问题，只能选用L-BFGS了，此时的 m 可以选择30左右。
下面简单介绍一下思路，回顾BFGS的更新步骤， G_{k+1} \approx H_k ：
G_{k+1}= G_k -\frac{G_k \delta_k \delta_k^T G_k}{\delta_k^TG_k \delta_k} + \frac{y_k y_k^T}{y_k^T \delta_k}\\
与DFP直接近似 H_k^{-1} 不同，这里直接对 G_{k+1} \approx H_k 求逆得到 H_{k}^{-1} 的近似，主要是利用了Sherman-Morrison公式，当 A 可逆时，有：
(A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{1 + v^TA^{-1}u} \\
利用Sherman-Morrison公式，分别利用两次可以得到（令 B_k = G_k +\frac{y_k y_k^T}{y_k^T \delta_k} ）：
B_k^{-1} = \left( G_k +\frac{y_k y_k^T}{y_k^T \delta_k} \right)^{-1} = G_k^{-1} - \frac{G_k^{-1}y_k y_k^T G_k^{-1}}{y_k^T\delta_k + y_k^TG_k^{-1}y_k} \\
那么有（省略号部分是将 B_k^{-1} 代入即可）：
H_{k+1}^{-1}= (G_k)^{-1} \\ = \left( G_k -\frac{G_k \delta_k \delta_k^T G_k}{\delta_k^TG_k \delta_k} + \frac{y_k y_k^T}{y_k^T \delta_k} \right)^{-1} \\= \left( B_k - \frac{G_k \delta_k \delta_k^T G_k}{\delta_k^TG_k \delta_k}\right)^{-1} \\ = B_k^{-1} - \frac{B_k^{-1}G_k\delta_k \delta_k^T G_k B_k^{-1}}{\delta_k^TG_k\delta_k - \delta_k^T G_k B_k^{-1}G_k\delta_k} \\ = \cdots \\ =\left( I - \frac{y_k \delta_k^T}{y_k^T\delta_k} \right) H_k^{-1} \left( I - \frac{\delta_k y_k^T}{y_k^T \delta_k} \right) + \frac{y_k y_k^T}{y_k^T \delta_k}\\
最后可以得到的式子是：
H_{k+1}^{-1} = \left( I - \frac{y_k \delta_k^T}{y_k^T\delta_k} \right) H_k^{-1} \left( I - \frac{\delta_k y_k^T}{y_k^T \delta_k} \right) + \frac{y_k y_k^T}{y_k^T \delta_k} \\
再重复一遍，这个是使用 G_{k+1} 近似 H_k ，然后对 G_{k+1} 求逆得到的结果，而不是类似于DFP直接近似 H_k^{-1} 得到的结果。
那么下面L-BFGS则容易推导了，对上述公式里面引入一些符号：
V_k =I - \frac{y_k \delta_k^T}{y_k^T\delta_k} \\ \rho_k = \frac{1}{y_k^T\delta_k}
那么更新公式可以写为：
H_{k+1}^{-1} = V_k^T H_k^{-1}V_k + \rho_k y_k y_k^T \\
那么迭代m + 1次得到：
H_{k+1}^{-1} = V_k^TH_k^{-1}V_k + \rho_k y_ky_k^T \\ = V_k^T\left( V_{k-1}^TH_{k-1}^{-1}V_{k-1} + \rho_{k-1}y_{k-1}y_{k-1}^T \right) + \rho_k y_k y_k^T \\ = \cdots \\ = V_k^T\cdots V_{k-m}^T H_{k-m}^{-1}V_{k-m}\cdots V_k \\ + \rho_{k-m} V_k^T \cdots V_{k-m+1}^T y_{k-m} y_{k-m}^T V_{k-m+1} \cdots V_k \\ + \cdots \\ + \rho_k y_k y_k^T
最后，将 H_{k-m}^{-1} 使用固定的 H_0^{-1} 来替代，即得到：
H_{k+1}^{-1} \\ =V_k^T\cdots V_{k-m}^T H_{0}^{-1}V_{k-m}\cdots V_k \\ + \rho_{k-m} V_k^T \cdots V_{k-m+1}^T y_{k-m} y_{k-m}^T V_{k-m+1} \cdots V_k \\ + \cdots\\ + \rho_k y_k y_k^T
一般来说， H_0^{-1} 可以取对角矩阵，然后再把历史过程中的 \{y_j, \delta_j \}_{j=k-m}^k 存下来即可，这样所需要的存储空间只有 O(mn) ，而不用将整个 H_k^{-1} 存储下来了。当然计算的过程中多涉及到对角矩阵乘法，都是稀疏矩阵，可以大大优化运算速度。
总之，L-BFGS就是使用历史的 m 条梯度信息来近似计算Hessian矩阵的逆，并没有显示构造或者存储Hessian矩阵或者其逆矩阵的近似矩阵。
总结一下本文，本文主要介绍了SGD、牛顿法、拟牛顿法和IIS等优化算法，主要是一阶、二阶梯度优化，并给出了详细的公式推导。 为什么可以用k-m+1~k的参数来近似H而不是0~k的所有参数呢
页:
[1]