常用优化方法及其原理

机器学习的各类算法中,常常需要用到最优化算法来求解参数.
这篇文章将总结一下机器学习中常用到的最优化算法,偏数理.
本文将会提到以下5种最优化算法

1 梯度下降法(gradient descent)
2 牛顿法(Newton’s method)
3 拟牛顿法(Quasi-Newton Methods)
3.1 DFP算法
 3.2 BFGS算法
 3.3 L-BFGS算法

梯度下降法

梯度下降方法的导出

梯度下降法最核心的部分在于:”当 $f(x)$ 在 $P_{0}$ 点可微时, $f(x)$ 在 $P_{0}$ 的梯度方向是 $f$ 的值增长最快的方向,且沿该方向的变化率就是梯度的模.”
那么相应的,减少最快的方向就是 $f(x)$ 在 $P_{0}$ 的负梯度方向了
下面给出其证明:

方向导数衡量函数在某特定方向l上的变化率,其定义为:
$f_{l}(P_{0})=\lim_{\rho \rightarrow 0}\frac{f(P)-f(P_{0})}{\rho }$
其中 $P$ 为方向 $l$ 上一点, $\rho$ 为 $P$ 与 $P_{0}$ 的距离
可以证明(用一阶泰勒展开证),方向导数可表示为梯度与 $l$ 的方向余弦的内积,即:
$\begin{align*} f_{l}(P_{0})&=f_{x_{1}}(P_{0})cos(\alpha )+f_{x_{2}}(P_{0})cos(\beta )+f_{x_{3}}(P_{0})cos(\gamma )+...\\ &=\triangledown f(P_{0})^{T}\cdot l_{0}\\ &=\left \| \triangledown f(P_{0}) \right \|\cdot \left \| l_{0} \right \|\cdot cos(\theta ) \end{align*}$
其中,限制 $l_{0}$ 为 $l$ 方向上的单位向量,它的模为1; $cos(\theta)$ 为梯度向量与 $l_{0}$ 的夹角
这个其实挺容易理解,就是将 $l$ 方向上的梯度分解到各坐标轴方向.要看详细证明的可以看<数学分析>"方向导数与梯度"章节.
注意:要使这个式子成立,有前提: rho趋于0.
因为求方向导数本来就是求 rho趋于0时的极限.这也是后面讲到的迭代步长要小的原因
可以看到,当 $\theta$ 为0,即 $l$ 方向为梯度方向时,方向导数取最大值.且为梯度向量的模.
因此, $f$ 在 $P_{0}$ 的梯度方向为 $f$ 的值增长最快的方向.
显然,这是在欧氏度量意义下才正确.但不在欧氏空间下的梯度一般没什么意义.

下面用两个小例题解释一下什么叫” $f(x)$ 沿着梯度方向改变”:

例一: $f(\overrightarrow{x})=x_{1}^{2}+x_{2}^{2}$ ,求 $f(\overrightarrow{x})$ 在 $P_{0}(1,2)$ 处的梯度.
解: $\triangledown f(x)=(\frac{\partial f(x)}{\partial x_{1}},\frac{\partial f(x)}{\partial x_{2}})=(2x_{1},2x_{2}) \Rightarrow \triangledown f(P_{0})=(2,4)$
即 $P_{0}$ 处梯度方向是(2,4)
例二: $P_{0}(1,2)$ ,从 $P_{0}$ 发出一条射线 $l$ , $l$ 方向为(2,4),求 $l$ 方向上的点.
解: 显然,答案为
$P_{0}+\lambda l= \begin{bmatrix} 1\\ 2 \end{bmatrix} +\lambda \begin{bmatrix} 2\\ 4 \end{bmatrix},\lambda >0$

看到这应该已经非常直白了,f(x)沿着梯度方向改变,就是指其自变量x沿着梯度方向移动.
注意移动的步长要足够小,因为 $\rho \rightarrow 0^{+}$

很自然推出梯度下降算法:

梯度下降法的收敛性

梯度下降法在”f(x)是连续可微实函数”的条件下是收敛的,且收敛于极值点
严格的数学描述及证明看<最优化理论与算法第二版=""> 陈宝林 P285

牛顿法

上面讲的梯度下降法,我们提到它是用一阶泰勒展开来证明的
我们看一下它一阶泰勒展开是什么样子:
$f(\vec{x})=f(\vec{x}^{(k)}+\triangledown f(\vec{x}^{(k)})^{T}(\vec{x}-\vec{x}^{(k)}))+o(\vec{x}-\vec{x}^{(k)})$
那如果我们用二阶泰勒展开,又会得到什么呢?我们来看一下
$f(\vec{x})=f(\vec{x}^{(k)}+\triangledown f(\vec{x}^{(k)})^{T}(\vec{x}-\vec{x}^{(k)}))+\frac{1}{2}(\vec{x}-\vec{x}^{(k)})^{T}\triangledown ^{2}f(\vec{x}^{k})(\vec{x}-\vec{x}^{(k)})+o((\vec{x}-\vec{x}^{(k)})^{2})$
泰勒展开的实质是用多项式去逼近原函数.
对于梯度下降法,它只用了一阶泰勒展开,只有一次多项式的信息,因此要找到极值点只能通过方向导数,或者说斜率来寻找.
而对于二阶泰勒展开,我们有了二次多项.
二次多项式有个特点就是,它包含了极点的信息.由于该二次多项式是原函数的逼近,因此找到二次多项式的极值位置就能大概知道原函数的极值位置!
那我们在每次迭代中找近似函数的极值点就好了.
这样的优化可不是一星半点.我们可以想象,梯度下降就是往周围坡度最陡的方向一步一步走;而牛顿法是根据参照物直接找到了谷底大概在哪,然后一下子就跳了过去.
因此,牛顿法的收敛速度一般比梯度下降法快.
直观的看就是这样:

图中展示了 $f=sin(x)$ 在x=0.7处的泰勒展开.从x=0.7开始迭代,实际极大值点在C.步长设为0.1

梯度上升法一次迭代后到达A;
牛顿法一次迭代后到达B.

显然牛顿法收敛更快.

牛顿法的导出

根据上面的分析,我们进行下面推导

对f(x)进行泰勒展开,保留其二阶项,记为 $\phi (\vec{x})$ :
$\phi (\vec{x})=f(\vec{x}^{(k)}+\triangledown f(\vec{x}^{(k)})^{T}(\vec{x}-\vec{x}^{(k)}))+\frac{1}{2}(\vec{x}-\vec{x}^{(k)})^{T}\triangledown ^{2}f(\vec{x}^{k})(\vec{x}-\vec{x}^{(k)})$
$\phi (\vec{x})$ 对x求偏导并令其为0:
$\triangledown \phi (x)=\triangledown f(x^{(k)})+\triangledown ^{2}f(x^{(k)})(x-x^{(k)})=0\ \ \ \ \ \ \ \ \ \ \ (*)$
$\Rightarrow x=x^{(k)}-\triangledown ^{2}f(x^{(k)})^{-1}\cdot \triangledown f(x^{(k)})$
其中,
$\triangledown ^{2}f(x^{(k)})=H(x^{(k)})$
$H(x^{(k)})_{i,j}=\begin{bmatrix} \frac{\partial^{2}f}{\partial x_{i}\partial x_{j}} \end{bmatrix}_{n\times n}$
$H(X)$ 为海赛矩阵

看,迭代公式出来了,就是 $x^{(k+1)}=x^{(k)}-\triangledown ^{2}f(x^{(k)})^{-1}\cdot \triangledown f(x^{(k)})$

牛顿法流程

牛顿法的收敛性

非常值得注意的是,牛顿法的最后结果与初始值的选取有很大关系.
从算法提出过程来看,牛顿法是每次迭代都是为了到达近似函数的极值点,因此可能到达极大值点,也可能到达极小值点.
从迭代目标公式可以看到,迭代方向由”负梯度”与”海赛矩阵”共同构成,因此若海赛矩阵非正定,目标函数值可能会上升.若初始点远离极值点,也可能不收敛.

举个简单的例子,上面的 $f=sin(x)$
若从x=0开始迭代,一阶导为1,二阶导=0,因此它会不知道往哪走,干脆就原地不动了,因此它不会收敛;
若从x=0.5开始迭代,一阶导>0,二阶导<0,因此它会往右走,最后到达极大值点;
若从x=-0.5开始迭代,一阶导>0,二阶导>0,因此它会往左走,最后到达极小值点;

拟牛顿法

牛顿法收敛速度很快,但也有明显的缺陷.海赛矩阵还有可能是奇异的,根本就就不了逆.就算可逆,每次迭代都需要计算海赛矩阵的逆矩阵.当特征成千上万维的时候,计算量会非常大,是不可接受的.而且于是人们希望找到n阶矩阵来替代海赛矩阵的逆矩阵.这样得到的方法就叫拟牛顿法.不同的替代矩阵就得到不同的算法.
要想替代海赛矩阵的逆矩阵,必须先知道海赛矩阵有什么性质.

对称.显然海赛矩阵是对称矩阵,那么它的逆矩阵也是对阵矩阵.( $(A^{-1})^{T}=(A^{T})^{-1}=A^{-1}$
正定.要保证牛顿法搜索方向是下降方向,必须有海赛矩阵为正定的.当然,其逆矩阵也会是正定的.
从上面牛顿法的导出中我们有公式(*)
$\triangledown \phi (x)=\triangledown f(x^{(k+1)})+\triangledown ^{2}f(x^{(k+1)})(x-x^{(k+1)})$
由于 $\phi (\vec{x})$ 在 $x^{(k+1)}$ 附近是f(x)的近似,因此令 $x=x^{k}$ ,
有
$\Rightarrow \triangledown f (x^{(k)})\approx \triangledown f(x^{(k+1)})+\triangledown ^{2}f(x^{(k+1)})(x^{(k)}-x^{(k+1)})$
$\Rightarrow \triangledown f (x^{(k)})-\triangledown f(x^{(k+1)})\approx \triangledown ^{2}f(x^{(k+1)})(x^{(k)}-x^{(k+1)})$
令:
$y_{k}=\triangledown f (x^{(k+1)})-\triangledown f(x^{(k)})$
$\delta _{k}=x^{(k+1)}-x^{(k)}$
则有:
$y_{k}\approx H_{k+1}\delta _{k}\ \ \ \ \ \ \ \ \ \ \ (1)$
或
$H_{k+1}^{-1}y_{k}\approx \delta _{k}\ \ \ \ \ \ \ \ \ \ \ (2)$
公式(1)或公式(2)即为第三个,也是最重要的条件,称为”拟牛顿条件”.

如果能找到矩阵 $G_{k}$ 满足上面三个条件,我们就认为可以将其作为海赛矩阵 $H_{k}^{-1}$ 的近似.
假设我们现在已经找到了一个初始的近似矩阵 $G_{k}$ .由于每轮迭代中,海赛矩阵都会更新,因此我们也要考虑该矩阵 $G_{k}$ 的迭代更新如何进行.
不妨假设:
$G_{k+1}=G_{k}+\Delta G_{k}$

接下来我们看一下不同的 $G_{k}$ 的构造.

在这之前,我们先整理一下思路.在计算 $H_{k+1}^{-1}$ 之前,我们已知的有 $y_{k}$ , $\delta _{k}$ , $G_{k}$ .我们希望构造出来的 $G_{k+1}$ 与已知变量这些有关.
$G_{0}$ 取正定对称矩阵.最简单的当然是单位矩阵.因此一般取 $G_{0}$ 为单位矩阵.

DFP算法(DFP algorithm)

DFP算法的导出

DFP算法核心在于寻找矩阵 $G_{k}$ 去逼近海赛矩阵的逆矩阵 $H_{k}^{-1}$
即令其满足公式(2)
$G_{k}$ 的寻找我们可以用待定系数法
由于要求 $G_{k}$ 为对称矩阵,我们不妨假设:
$G_{k+1}=G_{k}+\alpha uu^{T}+\beta vv^{T}\ \ \ \ \ \ \ \ \ \ \ (3)$
由于这里面有四个未知数,我们对它一无所知.因此不可避免地我们起码需要作四个假设才能把它们找出来.
回过头来看待定的式子, $G_{k+1}$ 要满足拟牛顿条件(公式(1)),那我们两边先乘个 $y_{k}$ 看看:
$G_{k+1}\cdot y_{k}=G_{k}\cdot y_{k}+\alpha uu^{T}\cdot y_{k}+\beta vv^{T}\cdot y_{k}$
稍微换一下位置:
$\Rightarrow G_{k+1}\cdot y_{k}=G_{k}\cdot y_{k}+u(\alpha u^{T}\cdot y_{k})+v(\beta v^{T}\cdot y_{k})$
我们会发现 $(\alpha u^{T}\cdot y_{k})$ , $(\beta v^{T}\cdot y_{k})$ 都是常数.
既然是数,我们也不妨添加两个假设:

$\alpha u^{T}\cdot y_{k}=1$
$\beta v^{T}\cdot y_{k}=1$

从这两个假设中我们可以得到:
$\alpha =\frac{1}{u^{T}y_{k}}$
$\beta =\frac{1}{v^{T}y_{k}}$
有了这两个假设,我们的式子变成了这样:
$G_{k+1}\cdot y_{k}=G_{k}\cdot y_{k}+u+v$
若满足拟牛顿条件,则有:
$G_{k+1}\cdot y_{k}=G_{k}\cdot y_{k}+u+v=\delta _{k}$
还需要两个假设,那我们可以大胆地令:

$u=\delta _{k}$
$v=-G_{k}\cdot y_{k}$

好啦,四个假设用光了,我们将四个假设代回公式(3),我们有: $G_{k+1}=G_{k}+\frac{\delta _{k}\delta _{k}^{T}}{\delta _{k}^{T}y_{k}}-\frac{G_{k}y_{k}y_{k}^{T}G_{k}^{T}}{y_{k}^{T}G_{k}^{T}y_{k}}$
由于 $G_{k}$ 为对称阵, $G_{k}^{T}=G_{k}$ ,上面的式子可以稍微化简:

$G_{k+1}=G_{k}+\frac{\delta _{k}\delta _{k}^{T}}{\delta _{k}^{T}y_{k}}-\frac{G_{k}y_{k}y_{k}^{T}G_{k}}{y_{k}^{T}G_{k}y_{k}}$

这就是我们找到的迭代目标.
只要给出一个正定对称阵作为初始矩阵 $G_{0}$ ,在之后的迭代里, $G_{k}$ 都将满足上面的三个条件.认为其可以替代海赛矩阵的逆矩阵.
这就是DFP算法的核心

DFP算法流程

BFGS算法

DFP算法虽然很厉害,但却很快被BFGS算法代替.实践证明BFGS算法性能更好.因此也是目前流行的拟牛顿法.

BFGS算法的导出

BFGS算法核心在于寻找矩阵 $B_{k}$ 去逼近海赛矩阵 $H_{k}$
即令其满足公式(1).
与上面的DFP算法一样,先用待定系数法,假设
$B_{k+1}=B_{k}+\alpha uu^{T}+\beta vv^{T}\ \ \ \ \ \ \ \ \ \ \ (4)$
两边乘 $\delta _{k}$ ,然后作以下四个假设:

$\alpha u^{T}\cdot y_{k}=1$
$\beta v^{T}\cdot y_{k}=1$
$u=y_{k}$
$v=-B_{k}\cdot \delta _{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}}$

BFGS算法流程

L-BFGS算法

L-BFGS算法看名字就知道是BFGS算法的优化.BFGS算法又有什么地方需要优化呢?
显然BFGS算法把计算量降了下来,但仍然需要储存矩阵 $G_{k}$ 或者 $B_{k}$ .
这个矩阵是非常大的,我们来计算一下. 假设特征有十万维,每个数字占用4字节.那么这个矩阵占用的内存是:
$\frac{4*10^{5}*10^{5}}{2^{10}*2^{10}*2^{10}}=37.25(GB)$
而十万维在NLP中根本就是稀松平常的事情,甚至说远不止十万维.因此这就是需要优化的地方.
L-BFGS算法的思想是不储存 $G_{k}$ 或 $B_{k}$ ,而是用的时候再通过递推公式计算.因此需要储存的只有 $\begin{Bmatrix} y_{1},...,y_{k} \end{Bmatrix}$ , $\begin{Bmatrix} \delta _{1},...,\delta _{k} \end{Bmatrix}$ 和 $G_{0}$ .
甚至仅储存最新的前m个:
$\begin{Bmatrix} y_{k-m+1},...,y_{k} \end{Bmatrix}$ 和 $\begin{Bmatrix} \delta _{k-m+1},...,\delta _{k} \end{Bmatrix}$
具体的计算推导可以参考这篇文章传送门.