KaaPexei 发表于 2023-2-27 13:44

机器学习优化算法(一)

在机器学习领域,数理统计、概率论和优化理论等知识是必不可少的,尤其是以优化理论最为重要。机器学习经常会把问题抽象为一个优化问题,比如最大化回报奖赏、最小化分类错误带来的损失或者最大化似然等等,而如何求解这些问题就必须有相应的优化算法。
本文将主要介绍几个常见的优化算法:

[*]随机梯度下降(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等优化算法,主要是一阶、二阶梯度优化,并给出了详细的公式推导。

unityloverz 发表于 2023-2-27 13:49

为什么可以用k-m+1~k的参数来近似H而不是0~k的所有参数呢
页: [1]
查看完整版本: 机器学习优化算法(一)