大多数机器学习算法都涉及某种形式的优化。优化指的是改变 x\mathbf x 以最小化或最大化 f(x)f(\mathbf x) 的任务。我们通常以最小化 f(x)f(\mathbf x) 指代大多数最优化问题,最大化可经由最小化 f(x)-f(\mathbf x) 来实现。

函数极值

考虑无约束优化问题

minxf(x)\min\limits_{\mathbf x} f(\mathbf x)

假定目标函数 f(x)f(\mathbf x) 二阶连续可微,自变量 xRn\mathbf x\in\R^n。一元函数对极值的推导可以推广到多元函数,函数取得极值的必要条件是导数为0。无约束的优化问题中

f(x)=f(x)x=0\nabla f(\mathbf x)=\frac{\partial f(\mathbf x)}{\partial\mathbf x}=0

梯度为0的点称为临界点(critical point)或驻点(stationary point)。下图为三种驻点:极大、极小和鞍点

极大或极小值的判定,取决于该函数的二阶导数,Hessian矩阵是由多元函数的二阶偏导数组成的矩阵

H(x)=(2fx122fx1x22fx1xn2fx2x12fx222fx2xn2fxnx12fxnx22fxn2)\mathbf H(\mathbf x)=\begin{pmatrix} \cfrac{\partial^2 f}{\partial x_1^2}&\cfrac{\partial^2 f}{\partial x_1\partial x_2}&\cdots&\cfrac{\partial^2 f}{\partial x_1\partial x_n} \\ \cfrac{\partial^2 f}{\partial x_2\partial x_1}&\cfrac{\partial^2 f}{\partial x_2^2}&\cdots&\cfrac{\partial^2 f}{\partial x_2\partial x_n} \\ \vdots &\vdots &\ddots &\vdots \\ \cfrac{\partial^2 f}{\partial x_n\partial x_1}&\cfrac{\partial^2 f}{\partial x_n\partial x_2}&\cdots&\cfrac{\partial^2 f}{\partial x_n^2} \\ \end{pmatrix}

一般情况下,多元函数的混合二阶偏导数与求导次序无关

2fxixj=2fxjxi\cfrac{\partial^2 f}{\partial x_i\partial x_j}=\cfrac{\partial^2 f}{\partial x_j\partial x_i}

因此Hessian矩阵是一个对称矩阵。多元函数极值的判别如下

  • 如果Hessian矩阵在该点是正定的,则函数在该点有极小值;
  • 如果Hessian矩阵在该点是负定的,则函数在该点有极大值;
  • 如果Hessian矩阵在该点是不定的,则该点为鞍点,它在一个方向上具有极小值,在另一个方向上具有极大值。

对于 nn 阶矩阵 A\mathbf A ,若对任意的nn维非0向量 x\mathbf x 都有
(1) xTAx>0\mathbf x^T\mathbf A\mathbf x>0 ,则称矩阵A\mathbf A为正定矩阵
(2) xTAx<0\mathbf x^T\mathbf A\mathbf x<0 ,则称矩阵A\mathbf A为负定矩阵
xTAx0\mathbf x^T\mathbf A\mathbf x\geqslant 0 ,则称矩阵A\mathbf A为半正定矩阵

在许多情况下,找解析解是一个很困难的问题,这就迫使我们使用数值方法找近似解,最常用的求解算法包括梯度下降法,牛顿法,拟牛顿法等。

拉格朗日乘数法

等式约束

拉格朗日乘数法(Lagrange multipliers)是一种寻找多元函数在一组约束下的极值的方法。通过引入拉格朗日乘数,可将有 nn 个变量与 ll 个约束条件的最优化问题转化为具有 n+ln+l 个变量的无约束优化问题求解。

考虑约束优化问题

minxf(x)s.t.hi(x)=0,i=1,2,l\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf x) \\ \text{s.t.} &\quad h_i(\mathbf x)= 0 ,\quad i=1,2,\cdots l \end{aligned}

几何意义:对于一个具有等式约束的优化问题

minxf(x)s.t.h(x)=0\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf x) \\ \text{s.t.} &\quad h(\mathbf x)= 0 \end{aligned}

从几何角度看,该问题的目标是在由方程 h(x)=0h(\mathbf x)=0 确定的 n1n-1 维曲面上寻找目标函数的极小点。

如图所示, f(x)=cmf(\mathbf x)=c_m 为目标函数的等值线 (蓝色),约束曲线 h(x)=0h(\mathbf x)=0 (红色)。对于约束曲面上的任意点 x\mathbf x ,该点的梯度 h(x)\nabla h(\mathbf x) 正交于约束曲面。假设 x0\mathbf x_0 是约束曲线和等高线的交点,我们沿着 h(x)=0h(\mathbf x)=0 移动,总能走到更高或更低的等高线上,直至约束曲线和等高线相切。也就是说,极值点必是目标函数 f(x)f(\mathbf x) 等值线族中与约束曲线 h(x)h(\mathbf x) 能相切的那个切点 。

假设切点为 x\mathbf x^* ,因为两曲线在切点处必有公法线,则两曲线在切点处的梯度向量共线,即存在 λ0\lambda\neq0 ,使下式成立

f(x)+λh(x)=0s.t.h(x)=0\begin{aligned} &\nabla f(\mathbf x^*)+\lambda\nabla h(\mathbf x^*)=0 \\ &\text{s.t.}\quad h(\mathbf x^*)= 0 \end{aligned}

参数 λ\lambda 称为拉格朗日乘数(Lagrange multiplier)。定义拉格朗日函数(Lagrangian function)

L(x,λ)=f(x)+λh(x)L(\mathbf x,λ)=f(\mathbf x)+λh(\mathbf x)

不难发现原问题等价于

{xL(x,λ)=0λL(x,λ)=0\begin{cases} \nabla_{\mathbf x} L(\mathbf x,λ)=0 \\ \nabla_{\lambda} L(\mathbf x,λ)=0 \end{cases}

于是,原等式约束优化问题转化为对拉格朗日函数 L(x,λ)L(\mathbf x,λ) 的无约束优化问题。

简易证明(1) 求目标函数 f(x,y,z)f(x,y,z) 在一个约束条件 G(x,y,z)=0G(x,y,z)=0 的极值。

假设约束条件可确定隐函数

z=z(x,y)z=z(x,y)

带入目标函数,原问题就转化为二元函数

φ(x,y)=f(x,y,z(x,y))\varphi(x,y)=f(x,y,z(x,y))

的无条件极值问题。取极值的必要条件为

φx=fx+fzzx=0φy=fy+fzzy=0\frac{\partial\varphi}{\partial x}=f_x+f_z\frac{\partial z}{\partial x}=0 \\ \frac{\partial\varphi}{\partial y}=f_y+f_z\frac{\partial z}{\partial y}=0

另一方面将二元函数 z=z(x,y)z=z(x,y) 带入约束条件得恒等式

G(x,y,z(x,y))=0G(x,y,z(x,y))=0

x,yx,y 分别求导,得

Gx+Gzzx=0Gy+Gzzy=0G_x+G_z\frac{\partial z}{\partial x}=0 \\ G_y+G_z\frac{\partial z}{\partial y}=0

由此,可知

fxGx=fyGy=fzGz\frac{f_x}{G_x}=\frac{f_y}{G_y}=\frac{f_z}{G_z}

方程组形式为

fyGzfzGy=0fxGzfzGx=0fxGyfyGx=0f_yG_z-f_zG_y=0 \\ f_xG_z-f_zG_x=0 \\ f_xG_y-f_yG_x=0

即极值点处,目标函数的梯度 f=(fx,fy,fz)\nabla f=(f_x,f_y,f_z) 和约束条件的梯度 G=(Gx,Gy,Gz)\nabla G=(G_x,G_y,G_z) 外积为零

f×G=ihkfxfyfzGxGyGz=0\nabla f\times \nabla G=\begin{vmatrix} \mathbf i & \mathbf h & \mathbf k \\ f_x & f_y & f_z \\ G_x & G_y & G_z\end{vmatrix}=0

也就是极值点处 f\nabla fG\nabla G 共线。不妨设 f+λG=0\nabla f+\lambda \nabla G=0 ,则上述问题可转化为

{fx(x,y,z)+λGx(x,y,z)=0fy(x,y,z)+λGy(x,y,z)=0fz(x,y,z)+λGz(x,y,z)=0G(x,y,z)=0\begin{cases} f_x(x,y,z)+\lambda G_x(x,y,z)=0 \\ f_y(x,y,z)+\lambda G_y(x,y,z)=0 \\ f_z(x,y,z)+\lambda G_z(x,y,z)=0 \\ G(x,y,z)=0 \end{cases}

引进辅助函数

L(x,y,z,λ)=f(x,y,z)+λG(x,y,z)L(x,y,z,\lambda)=f(x,y,z)+\lambda G(x,y,z)

我们容易看出方程组等价于 L(x,y,z,λ)=0\nabla L(x,y,z,\lambda)=0 。函数 LL 称为拉格朗日函数,参数 λ\lambda 称为拉格朗日乘数。

(2) 求三元函数 f(x,y,z)f(x,y,z) 在两个约束条件 G(x,y,z)=0,H(x,y,z)=0G(x,y,z)=0,H(x,y,z)=0 的极值。

两个约束条件的交集可确定,参数形式可写为

{x=xy=y(x)z=z(x)\begin{cases} x=x \\ y=y(x) \\ z=z(x) \end{cases}

带入目标函数,原问题就转化为一元函数

φ(x)=f(x,y(x),z(x))\varphi(x)=f(x,y(x),z(x))

的无条件极值问题。取极值的必要条件为

φx=fx+fydydx+fzdzdx=0\frac{\partial\varphi}{\partial x}=f_x+f_y\frac{\mathrm dy}{\mathrm dx}+f_z\frac{\mathrm dz}{\mathrm dx}=0

另一方面,将 y=y(x),z=z(x)y=y(x),z=z(x) 带入约束条件得恒等式

G(x,y(x),z(x))=0H(x,y(x),z(x))=0G(x,y(x),z(x))=0 \\ H(x,y(x),z(x))=0

两式对 xx 求导,得

Gx+Gydydx+Gzdzdx=0Hx+Hydydx+Hzdzdx=0G_x+G_y\frac{\mathrm dy}{\mathrm dx}+G_z\frac{\mathrm dz}{\mathrm dx} =0 \\ H_x+H_y\frac{\mathrm dy}{\mathrm dx}+H_z\frac{\mathrm dz}{\mathrm dx} =0

再加上取极值的必要条件,我们知道以下三个梯度

f=(fx,fy,fz)G=(Gx,Gy,Gz)H=(Hx,Hy,Hz)\nabla f=(f_x,f_y,f_z) \\ \nabla G=(G_x,G_y,G_z) \\ \nabla H=(H_x,H_y,H_z)

在极值点处都与向量 (1,dydx,dzdx)(1,\dfrac{\mathrm dy}{\mathrm dx},\dfrac{\mathrm dz}{\mathrm dx}) 垂直。所以这三个梯度在极值点是共面的,即他们的混合积为零

[f G H]=fxfyfzGxGyGzHxHyHz=0[\nabla f\ \nabla G\ \nabla H]=\begin{vmatrix} f_x & f_y & f_z \\ G_x & G_y & G_z \\ H_x & H_y & H_z \end{vmatrix}=0

不妨设 f+λ1G+λ2H=0\nabla f+\lambda_1 \nabla G+\lambda_2 \nabla H=0 ,则原问题可转化为

f+λ1G+λ2H=0G=0,H=0\nabla f+\lambda_1 \nabla G+\lambda_2 \nabla H=0 \\ G=0,H=0

引进拉格朗日函数

L(x,y,z,λ1,λ2)=f(x,y,z)+λ1G(x,y,z)+λ2H(x,y,z)L(x,y,z,\lambda_1,\lambda_2)=f(x,y,z)+\lambda_1G(x,y,z)+\lambda_2H(x,y,z)

则上述方程组等价于 L(x,y,z,λ1,λ2)=0\nabla L(x,y,z,\lambda_1,\lambda_2)=0

(3) 考虑高维优化问题

minxf(x)s.t.hi(x)=0,i=1,2,,l\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf x) \\ \text{s.t.}&\quad h_i(\mathbf x)= 0, \quad i=1,2,\cdots,l \end{aligned}

其中,定义域 DRnD\in\R^n 。现使用行向量

x=(x1,x2,,xn)=(x1,,xr,xr+1,xr+l)=(y,z)\mathbf x=(x_1,x_2,\cdots,x_n)=(x_1,\cdots,x_r,x_{r+1}\cdots,x_{r+l})=(\mathbf y,\mathbf z)

其中 yRr,zRl,n=r+l\mathbf y\in\R^r,\mathbf z\in\R^l,n=r+l 。并使用kk 维向量函数 h:DRl\mathbf h:D\to\R^l 记约束条件 h(x)=0\mathbf h(\mathbf x)= 0

优化问题可记作向量形式

minxf(y,z)s.t.h(y,z)=0\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf{y,z}) \\ \text{s.t.}&\quad \mathbf h(\mathbf{y,z})= 0 \end{aligned}

本节使用分子布局:分子为列向量。分子布局与分母布局只是相差一个转置。

若满足隐函数存在定理,约束条件可确定唯一隐函数

z=z(y)\mathbf z=\mathbf z(\mathbf y)

于是由复合函数的求导法则,知

hy+hzdzdy=0\frac{\partial\mathbf h}{\partial\mathbf y}+\frac{\partial\mathbf h}{\partial\mathbf z}\frac{\mathrm d\mathbf z}{\mathrm d\mathbf y}=0

取极值的必要条件为

fy+fzdzdy=0\frac{\partial f}{\partial\mathbf y}+\frac{\partial f}{\partial\mathbf z}\frac{\mathrm d\mathbf z}{\mathrm d\mathbf y}=0

于是得到极值点处

fy(hy)1fz(hz)1=0\frac{\partial f}{\partial\mathbf y}\left(\frac{\partial\mathbf h}{\partial\mathbf y}\right)^{-1}-\frac{\partial f}{\partial\mathbf z}\left(\frac{\partial\mathbf h}{\partial\mathbf z}\right)^{-1}=0

引入拉格朗日乘数向量 λ=(λ1,λ2,,λl)T\lambda=(\lambda_1,\lambda_2,\cdots,\lambda_l)^T,设

fy(gy)1=fz(gz)1=λT\frac{\partial f}{\partial\mathbf y}\left(\frac{\partial\mathbf g}{\partial\mathbf y}\right)^{-1}=\frac{\partial f}{\partial\mathbf z}\left(\frac{\partial\mathbf g}{\partial\mathbf z}\right)^{-1}=-\lambda^T

可转化为

fy+λThy=0fz+λThz=0\frac{\partial f}{\partial\mathbf y}+\lambda^T\frac{\partial\mathbf h}{\partial\mathbf y}=0 \\ \frac{\partial f}{\partial\mathbf z}+\lambda^T\frac{\partial\mathbf h}{\partial\mathbf z}=0

于是定义拉格朗日函数为

L(y,z,λ)=f(y,z)+λTh(y,z)L(\mathbf y,\mathbf z,λ)=f(\mathbf y,\mathbf z)+\lambda^T\mathbf h(\mathbf{y,z})

原问题等价于

L(x,λ)=0\nabla L(\mathbf x,λ)=0

拉格朗日乘数法的具体步骤如下:

(1) 定义拉格朗日函数

L(x,λ)=f(x)+i=1lλihi(x)L(\mathbf x,λ)=f(\mathbf x)+\sum_{i=1}^{l} λ_ih_i(\mathbf x)

其中 λiλ_i 称作拉格朗日乘数(Lagrange multiplier)。
(2) 令拉格朗日函数关于 x\mathbf x 和拉格朗日乘数 λ\lambda 的梯度等于零

{Lxj=0,j=1,2,,nLλi=0,i=1,2,,l\begin{cases} \dfrac{\partial L}{\partial x_j}=0, &j=1,2,\cdots,n \\ \dfrac{\partial L}{\partial \lambda_i}=0, &i=1,2,\cdots,l \end{cases}

(3) 求解得到的(n+l)(n+l)个方程,即可得到极值点。

拉格朗日乘数法所得的驻点会包含原问题的所有最小解,但并不保证每个驻点都是原问题的最小解。

不等式约束

现在进一步考虑不等式约束优化问题,先考虑一个不等式约束

minxf(x)s.t.g(x)0\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf x) \\ \text{s.t.}&\quad g(\mathbf x)\leqslant 0 \end{aligned}

引入拉格朗日函数

L(x,μ)=f(x)+μg(x)L(\mathbf x,\mu)=f(\mathbf x)+\mu g(\mathbf x)

图中阴影部分代表不等式约束表示的可行域,我们可以根据目标函数 f(x)f(\mathbf x) 的最优解 x\mathbf x^* 是否在可行域内将这类不等式约束优化问题分为两类:

  • g(x)<0g(\mathbf x^*)<0 即最优点在不等式区域内,此时不等式约束不起作用,可直接通过求解 f(x)=0\nabla f(\mathbf x)=0 求解最优点。这等价于设置 μ=0\mu=0 时,求解 L(x,μ)=0\nabla L(\mathbf x^*,\mu)=0 获得最优点。
  • g(x)=0g(\mathbf x^*)=0 即最优点在区域边界上,此时不等式约束退化为等式约束。若要在边界上取得极小值,等值面的梯度必定是指向 g(x)<0g(\mathbf x)<0 区域内,而 g(x)=0g(\mathbf x)=0 的梯度显然是向外。因此,最优解处 f(x)\nabla f(\mathbf x^*) 的方向必须与 g(x)\nabla g(\mathbf x^*)相反,即存在常数 μ>0\mu>0 使得 L(x,μ)=0\nabla L(\mathbf x^*,\mu)=0

整合这两种情形,必满足 μg(x)=0\mu g(\mathbf x)=0。因此,原问题可转化成在如下约束条件下最小化拉格朗日函数

{xL=f(x)+μg(x)=0g(x)0μ0μg(x)=0\begin{cases} \nabla_{\mathbf x}L=\nabla f(\mathbf x)+\mu \nabla g(\mathbf x)=0 \\ g(\mathbf x)\leqslant 0 \\ \mu \geqslant 0 \\ \mu g(\mathbf x)=0 \end{cases}

这些条件合称为 KKT (Karush-Kuhn-Tucker) 条件。不等式约束优化问题中的拉格朗日乘数也称为KKT 乘数。KKT条件将Lagrange乘数法所处理涉及等式的约束优化问题推广至不等式。

注意,拉格朗日乘数在不等式约束中出现,不再是不受限的。

上述做法可推广到多个不等式约束

minxf(x)s.t.gi(x)0i=1,2,,k\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf x) \\ \text{s.t.}&\quad g_i(\mathbf x)\leqslant 0 \quad i=1,2,\cdots,k \end{aligned}

  • 若最优点 x\mathbf x^* 在不等式区域内 gi(x)<0g_i(\mathbf x^*)<0 。此时不等式约束不起作用,可直接通过求解 f(x)=0\nabla f(\mathbf x)=0 求解最优点。这等价于拉格朗日乘子全部为零 μi=0\mu_i=0 ,求解 L(x,μ)=0\nabla L(\mathbf x^*,\mu)=0 获得最优点。
  • 若最优点 x\mathbf x^* 位于第 tt 个边界上 gt(x)=0g_t(\mathbf x^*)=0 。此时不等式约束退化为一个等式约束,f(x)+μtgt(x)=0\nabla f(\mathbf x^*)+\mu_t\nabla g_t(\mathbf x^*)=0,常数 μt>0\mu_t>0 。等价于拉格朗日乘子可全部设置为零,求解 L(x,μ)=0\nabla L(\mathbf x^*,\mu)=0 获得最优点。
  • 若最优点 x\mathbf x^* 位于多个边界的交界面上。此时不等式约束退化为等式约束,交界面对应的拉格朗日乘子大于零,其他拉格朗日乘子可全部设置为零,求解 L(x,μ)=0\nabla L(\mathbf x^*,\mu)=0 获得最优点。

以上约束条件必满足 μigi(x)=0\mu_i g_i(\mathbf x)=0,称为互补松弛性(complementary slackness)。因此,原问题可转化成在如下约束条件下最小化拉格朗日函数

{xL=f(x)+i=1kμigi(x)=0gi(x)0μi0μigi(x)=0\begin{cases} \displaystyle\nabla_{\mathbf x}L=\nabla f(\mathbf x)+\sum_{i=1}^k\mu_i \nabla g_i(\mathbf x)=0 \\ g_i(\mathbf x)\leqslant 0 \\ \mu_i \geqslant 0 \\ \mu_i g_i(\mathbf x)=0 \end{cases}

一般地,考虑具有 kk 个不等式约束和 ll 个等式约束的的优化问题:

minxf(x)s.t.hi(x)=0,i=1,2,lgj(x)0,j=1,2,k\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf x) \\ \text{s.t.} &\quad h_i(\mathbf x)= 0 ,\quad i=1,2,\cdots l \\ &\quad g_j(\mathbf x)\leqslant 0 ,\quad j=1,2,\cdots k \end{aligned}

引入拉格朗日乘子 λ=(λ1,λ2,,λl)T\lambda=(\lambda_1,\lambda_2,\cdots,\lambda_l)^Tμ=(μ1,μ2,,μk)T\mu=(\mu_1,\mu_2,\cdots,\mu_k)^T, 相应的拉格朗日函数为

L(x,λ,μ)=f(x)+i=1lλihi(x)+j=1kμjgj(x)L(\mathbf x,\lambda,\mu)=f(\mathbf x)+\sum_{i=1}^l\lambda_i h_i(\mathbf x)+\sum_{j=1}^k\mu_j g_j(\mathbf x)

其KKT条件为 (i=1,2,,l;j=1,2,,k)(i=1,2,\cdots,l;j=1,2,\cdots,k)

{xf(x)+i=1lλixhi(x)+j=1kμjxgj(x)=0hi(x)=0gj(x)0μj0μjgj(x)=0\begin{cases} \displaystyle\nabla_{\mathbf x}f(\mathbf x)+\sum_{i=1}^l\lambda_i \nabla_{\mathbf x}h_i(\mathbf x)+\sum_{j=1}^k\mu_j \nabla_{\mathbf x}g_j(\mathbf x)=0 \\ h_i(\mathbf x)=0 \\ g_j(\mathbf x)\leqslant 0 \\ \mu_j \geqslant 0 \\ \mu_j g_j(\mathbf x)=0 \end{cases}

KKT 条件是非线性规划(nonlinear programming)最优解的必要条件而不是充分条件。在实际应用上,KKT条件(方程组)一般不存在代数解,许多优化算法可供数值计算选用。

对偶问题

KKT条件虽然从理论上给出了极值的必要条件,但实际上方程解不一定好求。于是,我们又引入了它的对偶问题来辅助求解。对偶问题是利用拉格朗日对偶性将原始问题转换为对偶问题,通过解对偶问题得到原始问题的解。

原始问题(primal problem):考虑约束最优化问题

minxf(x)s.t.hi(x)=0,i=1,2,lgj(x)0,j=1,2,k\begin{aligned} \min\limits_{\mathbf x}& f(\mathbf x) \\ \text{s.t.} &\quad h_i(\mathbf x)= 0 ,\quad i=1,2,\cdots l \\ &\quad g_j(\mathbf x)\leqslant 0 ,\quad j=1,2,\cdots k \end{aligned}

首先引入拉格朗日函数

L(x,λ,μ)=f(x)+i=1lλihi(x)+j=1kμjgj(x)L(\mathbf x,\lambda,\mu)=f(\mathbf x)+\sum_{i=1}^l\lambda_i h_i(\mathbf x)+\sum_{j=1}^k\mu_j g_j(\mathbf x)

其中λi,μj\lambda_i,\mu_j 是拉格朗日乘数,μj0\mu_j \geqslant 0 。我们再定义函数

θp(x)=maxλ,μ;μ0L(x,λ,μ)\theta_p(\mathbf x)=\max_{\lambda,\mu;\mu\geqslant 0} L(\mathbf x,\lambda,\mu)

我们可以通过 x\mathbf x 是否满足原始问题的约束条件来分析函数 θp\theta_p

(1) 如果存在违反原始约束的条件,即有 hi(x)0h_i(\mathbf x)\neq 0gj(x)>0g_j(\mathbf x)>0 。此时可通过调整 λi\lambda_i 使得 λihi(x)+\lambda_i h_i(\mathbf x)\to+\infty ,或者调整 μj+\mu_j\to+\infty 使得 μjgj(x)+\mu_j g_j(\mathbf x)\to+\infty ,而将其余各 λi,μj\lambda_i,\mu_j 均取为零。那么就有

θp(x)=maxλ,μ;μ0[f(x)+i=1lλihi(x)+j=1kμjgj(x)]=+\theta_p(\mathbf x)=\max_{\lambda,\mu;\mu\geqslant 0} \left[f(\mathbf x)+\sum_{i=1}^l\lambda_i h_i(\mathbf x)+\sum_{j=1}^k\mu_j g_j(\mathbf x)\right]=+\infty

(2) 如果满足原始约束条件,则一定是令 gj(x)=0g_j(\mathbf x)=0 ,同时 hi(x)=0h_i(\mathbf x)=0 。因此 θp(x)\theta_p(\mathbf x) 为定值

θp(x)=maxλ,μ;μ0L(x,λ,μ)=f(x)\theta_p(\mathbf x)=\max_{\lambda,\mu;\mu\geqslant 0} L(\mathbf x,\lambda,\mu)=f(\mathbf x)

整合这两种情形

θp(x)={f(x),if hi(x)=0 and gj(x)0+,others\theta_p(\mathbf x)=\begin{cases} f(\mathbf x), & \text{if } h_i(\mathbf x)= 0\text{ and }g_j(\mathbf x)\leqslant 0\\ +\infty, & \text{others} \end{cases}

如果考虑极小化问题

minxθp(x)=minxmaxλ,μ;μ0L(x,λ,μ)=minxmax{+,f(x)}\min_{\mathbf x} \theta_p(\mathbf x)=\min_{\mathbf x}\max_{\lambda,\mu;\mu\geqslant 0}L(\mathbf x,\lambda,\mu)=\min_{\mathbf x}\max\{+\infty,f(\mathbf x)\}

所以无约束 minxθp(x)\min\limits_{\mathbf x} \theta_p(\mathbf x) 等价于在原始约束条件下的 minxf(x)\min\limits_{\mathbf x} f(\mathbf x) ,即它们有相同的解。问题 minxmaxλ,μ;μ0L(x,λ,μ)\min_{\mathbf x}\max_{\lambda,\mu;\mu\geqslant 0}L(\mathbf x,\lambda,\mu) 称为拉格朗日函数的 Min-Max 问题

综上所述,我们可以把原始问题中的约束条件去掉,得到原始问题的等价问题:

minxmaxλ,μL(x,λ,μ)s.t. μj0\min_{\mathbf x}\max_{\lambda,\mu} L(\mathbf x,\lambda,\mu) \\ \text{s.t. }\mu_j\geqslant 0

为了方便,定义原始问题的最优值

p=minxθp(x)p^*=\min_{\mathbf x} \theta_p(\mathbf x)

对偶问题(dual problem):主问题的优化一般比较困难,我们可以通过交换min-max 的顺序来
简化。定义对偶函数

θD(λ,μ)=minxL(x,λ,μ)\theta_D(\lambda,\mu)=\min_{\mathbf x} L(\mathbf x,\lambda,\mu)

对偶函数是以 λ,μ\lambda,\mu 为自变量的函数。对于给定的 λ,μ\lambda,\mu,在可行域内变动 x\mathbf x ,拉格朗日函数 LL 最小值就是 θD\theta_D 的值。对偶函数可令 Lx(x,λ,μ)=0\nabla L_{\mathbf x}(\mathbf x,\lambda,\mu)=0 求得。

再考虑极大化问题

maxλ,μ;μ0θD(λ,μ)=maxλ,μ;μ0minxL(x,λ,μ)\max_{\lambda,\mu;\mu\geqslant 0}\theta_D(\lambda,\mu)=\max_{\lambda,\mu;\mu\geqslant 0}\min_{\mathbf x} L(\mathbf x,\lambda,\mu)

称为拉格朗日函数的 Max-Min 问题,可表示为约束最优化为题

maxλ,μminxL(x,λ,μ)s.t. μj0\max_{\lambda,\mu}\min_{\mathbf x} L(\mathbf x,\lambda,\mu) \\ \text{s.t. }\mu_j\geqslant 0

称为原始问题的对偶问题。定义对偶问题的最优值

d=maxλ,μ;μ0θD(λ,μ)d^*=\max_{\lambda,\mu;\mu\geqslant 0}\theta_D(\lambda,\mu)

原始问题和对偶问题的关系:对任意的 x,λ,μ\mathbf x,\lambda,\mu 都有

θD(λ,μ)=minxL(x,λ,μ)L(x,λ,μ)maxλ,μ;μ0L(x,λ,μ)=θp(x)\theta_D(\lambda,\mu)=\min_{\mathbf x} L(\mathbf x,\lambda,\mu)\leqslant L(\mathbf x,\lambda,\mu)\leqslant \max_{\lambda,\mu;\mu\geqslant 0} L(\mathbf x,\lambda,\mu)=\theta_p(\mathbf x)

即对偶函数满足

θD(λ,μ)θp(x)\theta_D(\lambda,\mu)\leqslant\theta_p(\mathbf x)

则原始问题和对偶问题的最优值满足

d=maxλ,μ;μ0minxL(x,λ,μ)minxmaxλ,μ;μ0L(x,λ,μ)=pd^*=\max_{\lambda,\mu;\mu\geqslant 0}\min_{\mathbf x} L(\mathbf x,\lambda,\mu)\leqslant \min_{\mathbf x}\max_{\lambda,\mu;\mu\geqslant 0} L(\mathbf x,\lambda,\mu)=p^*

在某些条件下,d=pd^*=p^* ,由对偶问题的最优值即可得到原始问题的最优值。

  • 最优值显然满足 dpd^*\leqslant p^* ,这称为弱对偶性(weak duality)。pdp^*-d^* 称为最优对偶间隙(duality gap),是一个非负值。
  • d=pd^*=p^*,则称为强对偶性(strong duality)成立。强对偶成立必须满足以下条件:
    • 原始问题是凸优化问题:f(x)f(\mathbf x)gj(x)g_j(\mathbf x) 均为凸函数,hi(x)h_i(\mathbf x) 为仿射函数;
    • 可行域中至少有一点使不等式约束严格成立。

在强对偶性成立时,将拉格朗日函数分别对原变量和对偶变量求导,再并令导数等于霉,即可得到原变量与对偶变量的数值关系。

为什么要求对偶问题?

  1. 约束减少了,对偶问题只剩 n 个不等式约束;
  2. 主问题不一定是凸优化问题(局部最优不一定是全局最优),对偶问题一定是凸优化问题。

梯度下降法

梯度下降法的一系列算法在众多文章中表述均有差异,本文为了理解方便,以Adam算法公式为基准构造其他梯度下降系公式,符号均和Adam算法保持一致,其本质是一样的,最终结果经过变换也是相同的。统一初始点记为 x0\mathbf x_0,第一次迭代记为 t=1t=1

梯度下降法

梯度下降法(Gradient Descent,GD)也称最速下降法(steepest descent),是一种常用的一阶(first-order)优化方法,是求解无约束优化问题最简单、最经典的迭代方法之一。它被广泛应用于机器学习,是许多算法的基础,比如线性回归、逻辑回归,以及神经网络的早期实现。

考虑无约束优化问题

minxf(x)\min\limits_{\mathbf x} f(\mathbf x)

假定目标函数f(x)f(\mathbf x)连续可微,自变量 xRn\mathbf x\in\R^n

基本思路:梯度下降法希望通过不断执行迭代过程收敛到极小点,即构造一个序列 x0,x1,x2,\mathbf x_0,\mathbf x_1,\mathbf x_2,\cdots ,每次迭代都满足

f(xt)<f(xt1),t=1,2,f(\mathbf x_{t})<f(\mathbf x_{t-1}),\quad t=1,2,\cdots

假设第 tt 次迭代值为 xt\mathbf x_t 。梯度下降法使用泰勒一阶展开式[1]来近似代替目标函数

f(x)f(xt1)+ΔxTgtf(\mathbf x)\approx f(\mathbf x_{t-1})+\Delta\mathbf x^T\mathbf g_t

这里 gt=f(xt1)\mathbf g_t=\nabla f(\mathbf x_{t-1})f(x)f(\mathbf x)xt1\mathbf x_{t-1} 的梯度向量,Δx=xxt1\Delta\mathbf x=\mathbf x-\mathbf x_{t-1}

我们希望找到使 ff 下降的最快的方向,即求解

minΔxf(x)minΔxf(xt1)+ΔxTgt=minΔxΔxTgt=minΔxΔx2gt2cosθ\begin{aligned} \min_{\Delta\mathbf x}f(\mathbf x) &\approx\min_{\Delta\mathbf x}f(\mathbf x_{t-1})+\Delta\mathbf x^T\mathbf g_t \\ &= \min_{\Delta\mathbf x}\Delta\mathbf x^T\mathbf g_t \\ &=\min_{\Delta\mathbf x}\|\Delta\mathbf x\|_2\|\mathbf g_t\|_2\cos\theta \end{aligned}

上式中 θ\thetaΔx\Delta\mathbf x 和梯度向量的夹角。函数值 f(xt1)f(\mathbf x_{t-1}) 为常数,去除后不影响。

当移动方向与梯度方向相反时(即 cosθ=1\cos\theta=-1)取得最小值,也就是沿着负梯度曲线移动能最快到达极值点。

由于梯度曲线是连续曲线,每次迭代步长要无限小才是最短路径。考虑迭代次数有限,加入一个接近于0的正数控制下降的步幅。因此梯度下降法的迭代公式为

xt=xt1λgt\mathbf x_{t}=\mathbf x_{t-1}-\lambda\mathbf g_t

数值 λ>0\lambda>0 称为学习率 (learning rate),作用是控制下降的步幅。梯度向量 gt\mathbf g_t 控制下降的方向。

注意:初始位置的不同可能导致下降到不同的极值点。梯度下降不一定能够找到全局的最优解,有可能是一个局部最优解。当目标函数为凸函数时,局部极小点就对应着函数的全局最小点。

LocalMinimumGradientDescent

学习率:通过选取合适的步长 λ\lambda,就能确保通过梯度下降收敛到局部极小点。

  • 如果 λ\lambda 太小,梯度下降会起作用,但会很慢
  • 如果 λ\lambda 太大,梯度下降可能不断跨过最小值,永不收敛

使用梯度下降法时,通常建议尝试一系列 λ\lambda 值,对于每一个学习率画出少量迭代的代价函数,在尝试了一系列 λ\lambda 后,你可能会选择能快速且持续降低 ffλ\lambda 值。

Values of learning rate to try:0.0010.010.110.0030.030.330.0060.060.66\text{Values of learning rate to try:} \\ \begin{matrix} \cdots & 0.001 & 0.01 & 0.1 & 1 & \cdots \\ \cdots & 0.003 & 0.03 & 0.3 & 3 & \cdots \\ \cdots & 0.006 & 0.06 & 0.6 & 6 & \cdots \\ \end{matrix}

梯度向量gt\mathbf g_t 控制下降的方向,但同时也影响下降的步幅。当我们接近局部最小值时,导数从负值慢慢趋近于0,因此绝对值会自动变小。所以,即使学习率 λ\lambda 保持在某个固定值,更新的步幅也会自动变小。

收敛判断

  • 学习曲线(LearningCurve):横轴是梯度下降的迭代次数,纵轴代表目标函数 f(x)f(\mathbf x)。不同的应用场景中,梯度下降的收敛速度可能有很大差异。事实证明,我们很难事先知道梯度下降要经过多少次迭代才能收敛,所以可以先画个学习曲线后再训练模型。
  • 另一种方法是自动收敛测试 (Automatic convergence test):我们设置一个小容差 ϵ\epsilon (=0.001),如果代价函数在一次迭代中减少的量小于这个值,则可以认为它收敛了。

记住,收敛是指你找到了函数 ff 接近最小值的可能参数。选出正确的 ϵ\epsilon 是相当困难的,所以更倾向于使用学习曲线。

随机梯度下降

给定的数据集

D={(x1,y1),(x2,y2),,(xN,yN)}D=\{(\mathbf x_1,y_1),(\mathbf x_2,y_2),\cdots,(\mathbf x_N,y_N)\}

包含 NN 个样本,pp 个特征。其中,第 ii 个样本的特征向量为 xi=(xi1,xi2,,xip)T\mathbf x_i=(x_{i1},x_{i2},\cdots,x_{ip})^T 。代价函数为

J(θ)=1Ni=1NL(θ;xi,yi)J(\theta)=\frac{1}{N}\sum_{i=1}^N L(\theta;\mathbf x_i,y_i)

即损失函数 L(θ;xi,yi)L(\theta;\mathbf x_i,y_i) 的期望。准确计算这个期望的计算代价非常大,因为我们需要在整个数据集上的每个样本上评估模型。在实践中,我们可以从数据集中随机采样少量的样本,然后计算这些样本上的平均值来估计期望。

根据样本量在参数更新的准确性和执行更新所需的时间之间的权衡,梯度下降法有三种不同的应用方式。

批量梯度下降法(Batch Gradient Descent)是梯度下降法的最原始形式。为了获取准确的梯度,每一 步都在全部数据集上计算梯度,时间花费和内存开销都非常大,无法应用于大数据集、大模型的场景。

θt=θt1λ1Ni=1NθL(θ;xi,yi)\theta_t=\theta_{t-1}-\lambda\frac{1}{N}\sum_{i=1}^N \nabla_\theta L(\theta;\mathbf x_i,y_i)

随机梯度下降法(Stochastic Gradient Descent,SGD)是深度学习最常用的优化方法。SGD 放弃了对梯度准确性的追求,每一步采样单个样本来估计当前的梯度,计算速度快,内存开销小。但由于每步接受的信息量有限,随机梯度下降法对梯度的估计常常出现偏差,造成目标函数曲线收敛得很不稳定,伴有剧烈波动,有时甚至出现不收敛的情况。同时使用单个观测更新也可以在一定程度上增加不确定度,从而减轻陷入局部最小的可能。

小批量梯度下降法(Stochastic Gradient Descent,SGD with mini-batch)是批量梯度下降和随机梯度下降的折中方案,在一定程度上兼顾了以上两种方法的优点,通常缩写 SGD 指的是小批量梯度下降。每次从训练样本集上随机抽取一个小样本集,在抽出来的小样本集上迭代更新权重。小批量大小 bb(batch size)通常在 50 到 256 之间 (b<<Nb<<N),但可能因不同的应用而异。小批量梯度下降通常是训练神经网络时的首选算法。

优点:由于不是在全部训练数据上的损失函数,而是在每轮迭代中,随机优化某一条训练数据上的损失函数,这样每一轮参数的更新速度大大加快。

缺点

(1)准确度下降。由于即使在目标函数为强凸函数的情况下,SGD仍旧无法做到线性收敛。
(2)可能会收敛到局部最优,由于单个样本并不能代表全体样本的趋势。
(3)不易于并行实现。

一般来说SGD步长的选择比 batch GD的步长要小一点,因为梯度下降法使用的是准确梯度,所以它可以朝着全局最优解(当问题为凸问题时)较大幅度的迭代下去,但是随机梯度法不行,因为它使用的是近似梯度,或者对于全局来说有时候它走的也许根本不是梯度下降的方向,故而它走的比较缓,同样这样带来的好处就是相比于梯度下降法,它不是那么容易陷入到局部最优解中去。

Momentum

多维情况下,单个点处梯度的各个维度的分量往往差别很大,这在局部最优附近很常见。在这些情况下,SGD 往往在斜坡上振荡。这时,可以考虑梯度在每个维度上的分量变化,用历史的数据(惯性)修正下降方向,减少这种振荡。

Momentum(动量):一词借鉴了物理中的概念,是一种用来抑制振荡,加速梯度下降的技术。我们定义动量项为梯度方向和历史动量的加权求和

mt=βmt1+gt\mathbf m_t=\beta\mathbf m_{t-1}+\mathbf g_t

其中,gt=f(xt1)\mathbf g_t=\nabla f(\mathbf x_{t-1}) 为梯度向量,0<β<10<\beta<1 是指数衰减系数。加上动量项后的迭代公式为

xt=xt1λmt\mathbf x_{t}=\mathbf x_{t-1}-\lambda\mathbf m_t

动量法相当于每次迭代的时候,都会将之前的梯度考虑进来,每次的移动方向不仅取决于当前的梯度,还取决于过去各个梯度分量在各自方向上是否一致。如果一个梯度分量一直沿着当前方向进行更新,那么每次更新的幅度就越来越大。如果一个梯度分量在一个方向上不断改变,那么其更新幅度就会被衰减。这样我们就可以使用一个较大的学习率,得到了更快的收敛速度以及更小的震荡。

动量项:初始 m0=0\mathbf m_0=0 ,将动量项迭代展开

mt=gt+βgt1+β2gt2++βt1g1\mathbf m_t=\mathbf g_t+\beta\mathbf g_{t-1}+\beta^2\mathbf g_{t-2}+\cdots+\beta^{t-1}\mathbf g_{1}

由于 0<β<10<\beta<1 ,越远的梯度权重也会越小,一般可认为这个影响范围为 1/(1β)1/(1-\beta)。这里使用了指数加权移动平均的数学概念,只是省略了权重系数 (1β)(1-\beta),其本质相同。通常取系数 β=0.9\beta=0.9,相当于取了近10个梯度计算移动平均,因为再往前权重太小了,基本没影响了。

注意,记vt1=λmt1\mathbf v_{t-1}=\lambda\mathbf m_{t-1} ,动量项变为

vt=βvt1+λgt\mathbf v_t=\beta\mathbf v_{t-1}+\lambda\mathbf g_t

迭代公式为

xt=xt1vt\mathbf x_{t}=\mathbf x_{t-1}-\mathbf v_t

此时,会发现和其他文章公式一致。

指数加权移动平均(Exponentially Moving Average,EMA)是一种计算比较相近时刻数据的加权平均值。tt 时刻 EMA 计算公式是:

vt=βvt1+(1β)θtv_t=\beta v_{t-1}+(1-\beta)\theta_t

其中 θt\theta_ttt 时刻的真实值,0<β<10<\beta<1 为权重系数。我们将序列逐一迭代展开得到

v0=0v1=(1β)θ1v2=(1β)(θ2+βθ1)v3=(1β)(θ3+βθ2+β2θ1)vt=(1β)(θt+βθt1+β2θt2++βt1θ1)\begin{aligned} &v_0=0 \\ &v_1=(1-\beta)\theta_1 \\ &v_2=(1-\beta)(\theta_2+\beta\theta_1) \\ &v_3=(1-\beta)(\theta_3+\beta\theta_2+\beta^2\theta_1) \\ &\cdots \\ &v_t=(1-\beta)(\theta_t+\beta\theta_{t-1}+\beta^{2}\theta_{t-2}+\cdots+\beta^{t-1}\theta_1) \end{aligned}

可以看到,本质上就是加权系数呈指数衰减的移动平均,越靠近当前时刻的数值加权系数就越大。

NAG

NAG(Nesterov accelerated gradient)用历史数据和超前梯度修正下降方向,尤其在极值点附近,能有效减少振荡。

定义 Nesterov 动量

vt=βvt1+λxf(xt1βvt1)\mathbf v_{t}=\beta\mathbf v_{t-1}+\lambda\nabla_{\mathbf x} f(\mathbf x_{t-1}-\beta\mathbf v_{t-1})

迭代公式为

xt=xt1vt\mathbf x_{t}=\mathbf x_{t-1}-\mathbf v_{t}

Momentum 算法首先计算当前梯度 f(xt1)\nabla f(\mathbf x_{t-1}),然后沿累积梯度方向 vt1\mathbf v_{t-1} 进行大跳跃,如图中蓝色向量。但 NAG 首先平移到累积梯度方向 vt1\mathbf v_{t-1} 指向的点 (xt1βvt1)(\mathbf x_{t-1}-\beta\mathbf v_{t-1}) ,然后从这个点计算梯度 xf(xtβvt1)\nabla_{\mathbf x} f(\mathbf x_t-\beta\mathbf v_{t-1}),如图中红色向量。

Δxt1=xt1xt2=vt1\Delta\mathbf x_{t-1}=\mathbf x_{t-1}-\mathbf x_{t-2}=-\mathbf v_{t-1} 。则 NAG 迭代公式可以变换为

xt=xt1+βΔxt1λxf(xt1+βΔxt1)\mathbf x_{t}=\mathbf x_{t-1}+\beta\Delta\mathbf x_{t-1}-\lambda\nabla_{\mathbf x} f(\mathbf x_{t-1}+\beta\Delta\mathbf x_{t-1})

这可以理解为,先按上一步的方向移动,然后再按这个点的梯度移动。

AdaGrad

在梯度下降法中,学习率的选取十分重要,因为它关乎到优化过程每一步步幅大小,步幅太大会导致在最优解附近来回振荡,步幅太小容易陷入局部最优解。由于自变量每个维度上收敛速度都不相同,因此根据不同维度的收敛情况分别设置学习率。

AdaGrad(Adaptive Gradient,自适应梯度)是梯度下降法最直接的改进。AdaGrad 累积了到本次迭代为止梯度的历史值信息用于修正不同维度上的学习率

vt=vt1+gt2\mathbf v_{t}=\mathbf v_{t-1}+\mathbf g_t^2

这也是一个向量,计算时分别对向量的每个分量进行。其中 gt2=gtgt\mathbf g_t^2=\mathbf g_t\odot \mathbf g_t 为 Hadamard 积 [2],即梯度向量 gt\mathbf g_t 的各分量的平方组成的向量(square element-wise)。后面所有的计算公式都是对向量的每个分量进行。

迭代公式如下

xt=xt1λvt+ϵgt\mathbf x_{t}=\mathbf x_{t-1}-\frac{\lambda}{\sqrt{\mathbf v_{t}}+\epsilon}\mathbf g_t

其中 λ\lambda 是学习率, gt=f(xt1)\mathbf g_t=\nabla f(\mathbf x_{t-1}) 是第 tt 次迭代时的梯度向量,ϵ\epsilon是一个很小的正数,为了避免分母为零。和标准梯度下降法唯一不同的是多了分母中的这一项。默认取 λ=0.01,ϵ=108\lambda=0.01,\epsilon=10^{-8}

Adagrad的主要优点是它消除了手动调整学习率的需要。采用历史梯度平方和来修正学习率,历史数据变化的越多,学习率减少的越多。本质是解决各方向导数数值量级的不一致而振荡的问题。

直观的来看,在平缓山谷处,更为平缓的方向,会取得更大的进步(因为平缓,所以历史梯度平方和较小,对应学习下降的幅度较小),并且能够使得陡峭的方向变得平缓,从而加快训练速度。

另外,分母中求和的形式实现了退火过程[3],这是很多优化技术中常见的策略,意味着随着时间的推移,学习速率越来越小,从而保证了算法的最终收敛。

AdaGrad 算法的缺点是在经过一定次数的迭代依然没有找到最优点时,由于这时的学习率已经非常小,很难再继续找到最优点。

RMSprop

RMSprop(Root Mean Square Propagation,均方根传播)算法是 Adagrad 的改进,为了降低Adagrad中学习率衰减过快问题,使用梯度平方的指数加权平均限制梯度窗口

vt=βvt1+(1β)gt2\mathbf v_{t}=\beta\mathbf v_{t-1}+(1-\beta)\mathbf g_t^2

迭代公式

xt=xt1λvt+ϵgt\mathbf x_{t}=\mathbf x_{t-1}-\frac{\lambda}{\sqrt{\mathbf v_{t}}+\epsilon}\mathbf g_t

其中 λ\lambda 是学习率, gt=f(xt1)\mathbf g_t=\nabla f(\mathbf x_{t-1}) 是第 tt 次迭代时的梯度向量,ϵ\epsilon是一个很小的正数,为了避免分母为零。默认值 λ=0.001,β=0.9,ϵ=108\lambda=0.001,\beta=0.9,\epsilon=10^{-8} 。经验上,RMSProp被证明有效且实用的深度学习网络优化算法。

优势:能够克服AdaGrad梯度急剧减小的问题,在很多应用中都展示出优秀的学习率自适应能力。尤其在不稳定(Non-Stationary)的目标函数下,比基本的SGD、Momentum、AdaGrad表现更良好。

AdaDelta

AdaDelta 算法和RMSprop类似,也是对 Adagrad 的改进。通过使用梯度平方的指数衰减移动平均,降低学习率衰减过快问题。本节使用统计学符号表示,Adagrad 中的 vt\mathbf v_t 用数学期望 E[g2]t\mathbb E[\mathbf g^2]_t 表示,初始值 E[g2]0=0\mathbb E[\mathbf g^2]_0=0

E[g2]t=βE[g2]t1+(1β)gt2\mathbb E[\mathbf g^2]_t=\beta\mathbb E[\mathbf g^2]_{t-1}+(1-\beta)\mathbf g_t^2

更新量

Δxt=xtxt1=λE[g2]t+ϵgt\Delta\mathbf x_t=\mathbf x_t-\mathbf x_{t-1}=-\frac{\lambda}{\sqrt{\mathbb E[\mathbf g^2]_t}+\epsilon}\mathbf g_t

其中 λ\lambda 是学习率, gt=f(xt1)\mathbf g_t=\nabla f(\mathbf x_{t-1}) 是第 tt 次迭代时的梯度向量,ϵ\epsilon是一个很小的正数,为了避免分母为零。由于分母是梯度的均方根

RMS[g]t=E[g2]t+ϵ\text{RMS}[\mathbf g]_{t}=\sqrt{\mathbb E[\mathbf g^2]_t}+\epsilon

更新量可简写为

Δxt=λRMS[g]tgt\Delta\mathbf x_t=-\frac{\lambda}{\text{RMS}[\mathbf g]_{t}}\mathbf g_t

上式其实还是依赖于全局学习率的,AdaDelta算法采用了牛顿法的思想,寻找最优的全局学习率,并采用Hessian矩阵的对角线近似Hessian矩阵。(原理不理解)

AdaDelta 不需要提前设置全局学习率这一超参数,通过维护更新量 Δx\Delta\mathbf x 平方的指数衰减移动平均来动态确定学习率,初始值 E[Δx2]0=0\mathbb E[\Delta\mathbf x^2]_0=0

E[Δx2]t=βE[Δx2]t1+(1β)Δxt2\mathbb E[\Delta\mathbf x^2]_t=\beta\mathbb E[\Delta\mathbf x^2]_{t-1}+(1-\beta)\Delta\mathbf x_t^2

更新量的均方根为

RMS[Δx]t=E[Δx2]t+ϵ\text{RMS}[\Delta\mathbf x]_{t}=\sqrt{\mathbb E[\Delta\mathbf x^2]_t}+\epsilon

最终的迭代公式为

xt=xt1RMS[Δx]tRMS[g]tgt\mathbf x_t=\mathbf x_{t-1}-\frac{\text{RMS}[\Delta\mathbf x]_{t}}{\text{RMS}[\mathbf g]_{t}}\mathbf g_t

从上式可以看出,AdaDelta 算法将初始全局学习率改为动态计算的 RMS[Δx]t\text{RMS}[\Delta\mathbf x]_{t},在一定程度上平抑了学习率的波动。

**特点:**训练初中期,加速效果不错,很快。训练后期,反复在局部最小值附近抖动

Adam

Adam 是 Adaptive Moment estimation(自适应矩估计)的简称,与传统梯度下降保持同一个学习率不同,Adam通过计算梯度的一阶矩估计(First Moment Estimation)和二阶矩估计(Second Moment Estimation)而为 x\mathbf x 的不同分量设计独立的自适应性学习率。它通常比梯度下降快得多,已经成为实践者训练神经网络的行业标准。

Adam算法可以看做是修正后的动量法和RMSProp的结合。采用指数衰减移动平均(exponential decay average)来对一阶矩和二阶矩进行估计,时间久远的梯度对当前平均值的贡献呈指数衰减。

一阶矩 E[g]\mathbb E[\mathbf g] 表示梯度的均值,这体现了惯性保持,用来控制更新的方向

mt=β1mt1+(1β1)gt\mathbf m_{t}=\beta_1\mathbf m_{t-1}+(1-\beta_1)\mathbf g_t

二阶矩 E[g2]\mathbb E[\mathbf g^2] 表示梯度无中心方差(梯度平方的均值),为不同分量产生自适应的学习速率

vt=β2vt1+(1β2)gt2\mathbf v_{t}=\beta_2\mathbf v_{t-1}+(1-\beta_2)\mathbf g_t^2

其中 gt=f(xt1)\mathbf g_t=\nabla f(\mathbf x_{t-1}) 是第 tt 次迭代时的梯度向量, β1,β2[0,1)\beta_1,\beta_2\in[0,1) 为衰减系数,控制这些移动平均的指数衰减率(exponential decay rates)。 gt2=gtgt\mathbf g_t^2=\mathbf g_t\odot \mathbf g_t 为 Hadamard 积 [2:1],即梯度向量 gt\mathbf g_t 的各分量的平方组成的向量(square element-wise)。后面所有的计算公式都是对向量的每个分量进行。

然而, mt,vt\mathbf m_t,\mathbf v_t 是初始化为0的有偏估计,矩估计值会偏向0,尤其是在训练初期阶段。因此,为了缓解一阶矩和二阶矩初始为0带来的影响,Adam利用了偏差修正(bias-correct)

m^t=mt1β1t,v^t=vt1β2t\hat{\mathbf m}_t=\frac{\mathbf m_t}{1-\beta_1^t},\quad \hat{\mathbf v}_t=\frac{\mathbf v_t}{1-\beta_2^t}

现推导二阶矩 E[gt2]\mathbb E[\mathbf g^2_t]E[vt]\mathbb E[\mathbf v_t] 的关系,得到偏差修正项, vtv_t 代表 vt\mathbf v_t 的任意分量

E[vt]=E[(1β2)i=1tβ2tigi2]=E[gt2](1β2)i=1tβ2ti+ζ=E[gt2](1β2t)+ζ\begin{aligned} \mathbb E[v_t] &= \mathbb E\left[(1-\beta_2)\sum_{i=1}^t\beta_2^{t-i}\cdot g_i^2\right] \\ &=\mathbb E[g_t^2]\cdot(1-\beta_2)\sum_{i=1}^t\beta_2^{t-i}+\zeta \\ &=\mathbb E[g_t^2]\cdot(1-\beta_2^t)+\zeta \end{aligned}

如果真实二阶矩 E[gt2]\mathbb E[g_t^2] 是静态的(stationary),那么 ζ=0\zeta=0。通常可以忽略常数 ζ\zeta ,得到上述修正项。一阶矩估计的推导完全是相似的。

最后,迭代公式为

xt=xt1λv^t+ϵm^t\mathbf x_{t}=\mathbf x_{t-1}-\frac{\lambda}{\sqrt{\hat{\mathbf v}_t}+\epsilon}\hat{\mathbf m}_t

其中 λ\lambda 是学习率,ϵ\epsilon是一个很小的正数,为了避免分母为零。默认λ=0.001,β1=0.9,β2=0.999,ϵ=108\lambda=0.001,\beta_1=0.9,\beta_2=0.999,\epsilon=10^{-8}

物理意义

  • mt\|\mathbf m_t\| 大且 vt\mathbf v_t大时,梯度大且稳定,这表明遇到一个明显的大坡,前进方向明确;
  • mt\|\mathbf m_t\| 趋于零且 vt\mathbf v_t 大时,梯度不稳定,表明可能遇到一个峡谷,容易引起反弹震荡;
  • mt\|\mathbf m_t\| 大且 vt\mathbf v_t 趋于零时,这种情况不太可能出现;
  • mt\|\mathbf m_t\| 趋于零且 vt\mathbf v_t 趋于零时,梯度趋于零,可能到达局部最低点,也可能走到一片坡度极缓的平地,此时要避免陷入平原( plateau)。

优缺点:不需要手动指定学习率;通常收敛速度远大于梯度下降;适用于非稳态(non-stationary)目标;适用于解决包含很高噪声或稀疏梯度的问题;超参数可以很直观地解释,并且基本上只需极少量的调参。

虽然Adam算法是目前主流的优化算法,不过在很多领域里(如计算机视觉的对象识别、NLP中的机器翻译)仍然是使用动量(SGD with Momentum)的效果最佳。

AdaMax

AdaMax 是的Adam变体。在Adam中, vt\mathbf v_t 基于当前梯度平方和过去梯度平方的序列更新,过去梯度平方的权重呈反比例缩放。 vt\mathbf v_t 的任意分量 vtv_t

vt=β2vt1+(1β2)gt2=(1β2)(gt2+β2gt12+β22gt22++β2t1g12)\begin{aligned} v_{t}&=\beta_2v_{t-1}+(1-\beta_2)g_t^2 \\ &=(1-\beta_2)(g_t^2+\beta_2 g_{t-1}^2+\beta_2^2 g_{t-2}^2+\cdots+\beta_2^{t-1}g_1^2) \end{aligned}

其中 gtg_t 是第 tt 次迭代时的梯度向量 f(xt1)\nabla f(\mathbf x_{t-1}) 的分量。可见 vtv_t 基于梯度序列的 2\ell_2 范数[4]更新,我们可以将其推广到基于 LpL_p 范数的更新规则中

vt=β2pvt1+(1β2p)gtpv_{t}=\beta_2^p v_{t-1}+(1-\beta_2^p)|g_t|^p

注意,这里的衰减系数相应的变为 β2p\beta_2^p。然而,pp 值较大的范数会变得不稳定,但当 pp\to\infty 时通常会变得很稳定。为了避免与Adam混淆,我们使用 utu_t 表示无穷范数

ut=limp(β2pvt1+(1β2p)gtp)1/p=max(β2ut1,gt)u_t=\lim\limits_{p\to\infty}\left(\beta_2^p v_{t-1}+(1-\beta_2^p)|g_t|^p\right)^{1/p} =\max(\beta_2u_{t-1},|g_t|)

这里使用的最大值,不需要修正初始化偏差。mt\mathbf m_t 仍然和Adam一样

mt=β1mt1+(1β1)gt\mathbf m_{t}=\beta_1\mathbf m_{t-1}+(1-\beta_1)\mathbf g_t

最终得到简单的迭代公式

xt=xt1λutmt1β1t\mathbf x_t= \mathbf x_{t-1}-\frac{\lambda}{\mathbf u_t}\frac{\mathbf m_t}{1-\beta_1^t}

默认λ=0.002,β1=0.9,β2=0.999\lambda=0.002,\beta_1=0.9,\beta_2=0.999

Nadam

Adam 可以看做是修正后的动量法和RMSProp的结合,我们还可以看到 NAG 优于原版动量,Nadam(Nesterov-accelerated Adaptive Moment Estimation)从而结合了 Adam 和 NAG 。

Adam 动量项更新规则如下

mt=β1mt1+(1β1)gt\mathbf m_t =\beta_1\mathbf m_{t-1}+(1-\beta_1)\mathbf g_t

NAG 的核心在于使用超前点的信息,NAdam 中提出了一种公式变形的思路,即能达到 Nesterov 的效果。在假定连续两次的梯度变化不大的情况下,修改动量项为

mt=β1mt+(1β1)gt\mathbf m_t' =\beta_1\mathbf m_t+(1-\beta_1)\mathbf g_t

更新规则如下

xt=xt1λv^t+ϵ(β1m^t+(1β1)gt1β1t)\mathbf x_{t} =\mathbf x_{t-1}-\frac{\lambda}{\sqrt{\hat{\mathbf v}_t}+\epsilon}\left(\beta_1\hat{\mathbf m}_t+\frac{(1-\beta_1)\mathbf g_t}{1-\beta_1^t}\right)

牛顿法和拟牛顿法

牛顿法和拟牛顿法是求解无约束优化问题常用的方法,其迭代轮数远小于梯度下降法。

牛顿法

考虑无约束优化问题

minxf(x)\min\limits_{\mathbf x} f(\mathbf x)

假定目标函数 f(x)f(\mathbf x) 二阶连续可微,自变量 xRn\mathbf x\in\R^n。下面从两个角度来推导牛顿法(Newton method)迭代公式。

注意,本节牛顿法统一使用第 tt 轮的数据推导,而梯度下降法章节使用 t1t-1 轮推导。

梯度零点迭代:目标函数取极值的必要条件是其梯度为零 f(x)=0\nabla f(\mathbf x)=0 。以一元函数为例,我们可以迭代求解梯度的零值点。图中蓝色曲线为导数曲线 f(x)f'(x)

假设第 tt 次迭代值为 xtx_t ,如图,下一次迭代值使用f(x)f'(x)xtx_t 切线方向得到的零点值,即

f(xt)=f(xt)xtxt+1f''(x_{t})=\frac{f'(x_t)}{x_t-x_{t+1}}

于是,我们就得到

xt+1=xtf(xt)f(xt)x_{t+1}=x_t-\frac{f'(x_{t})}{f''(x_{t})}

二阶泰勒展开:牛顿法可以看作使用更精确的二阶泰勒展开式来近似代替目标函数。

假设第 tt 次迭代值为 xt\mathbf x_t 。根据泰勒二阶展开式有

f(x)f(xt)+gtT(xxt)+12(xxt)THt(xxt)f(\mathbf x)\approx f(\mathbf x_{t})+\mathbf g_t^T(\mathbf x-\mathbf x_{t}) +\frac{1}{2}(\mathbf x-\mathbf x_{t})^T\mathbf H_t(\mathbf x-\mathbf x_{t})

其中,gt\mathbf g_tf(x)f(\mathbf x) 的梯度向量f(x)\nabla f(\mathbf x)x=xt\mathbf x=\mathbf x_{t} 处的值。Ht\mathbf H_tf(x)f(\mathbf x) 的Hessian 矩阵H(x)\mathbf H(\mathbf x)x=xt\mathbf x=\mathbf x_{t} 的值。

f(x)=fx,H(x)=(2fxixj)n×n\nabla f(\mathbf x)=\frac{\partial f}{\partial\mathbf x},\quad \mathbf H(\mathbf x)=\left(\frac{\partial^2 f}{\partial x_i\partial x_j}\right)_{n\times n}

牛顿法使用近似式的极小值点来接近f(x)f(\mathbf x) 的极小点。目标函数取极值的必要条件是其梯度为零。于是,对 f(x)f(\mathbf x) 的二阶泰勒展开求梯度

f(x)gt+Ht(xxt)=0\nabla f(\mathbf x)\approx\mathbf g_t+\mathbf H_t(\mathbf x-\mathbf x_{t})=0

于是得到牛顿法的迭代公式

xt+1=xtHt1gt\mathbf x_{t+1}=\mathbf x_{t}-\mathbf H_t^{-1}\mathbf g_t

牛顿法是一种步长 λ=1\lambda=1 的迭代算法。其中

dt=Ht1gt\mathbf d_t=-\mathbf H_t^{-1}\mathbf g_t

称为牛顿方向。可以证明牛顿法是二次收敛的,但当初始点 x0\mathbf x_0 远离极小值时可能不收敛。

几何意义f(x)f(\mathbf x) 的泰勒二阶展开式为

f(x)f(xt)+(xxt)Tgt+12(xxt)THt(xxt)f(\mathbf x) \approx f(\mathbf x_{t})+(\mathbf x-\mathbf x_{t})^T\mathbf g_t +\frac{1}{2}(\mathbf x-\mathbf x_{t})^T\mathbf H_t(\mathbf x-\mathbf x_{t})

z=xxt\mathbf z=\mathbf x-\mathbf x_{t} ,构造函数 φ(z)\varphi(\mathbf z)

φ(z)=12zTHtz+gtTz+c\varphi(\mathbf z)=\frac{1}{2}\mathbf z^T\mathbf H_t\mathbf z+\mathbf g_t^T\mathbf z+c

f(x)=φ(xxt)f(\mathbf x)=\varphi(\mathbf x-\mathbf x_{t}) 的函数图像可通过平移 xt\mathbf x_{t} 得到。其中,Hessian 矩阵 Ht\mathbf H_t 为对称阵,c=f(xt1)c=f(\mathbf x_{t-1}) 为常数值。易知 φ(z)\varphi(\mathbf z) 的梯度和Hessian 矩阵分别为

φ(z)=Htz+gt,2φ(z)=Ht\nabla \varphi(\mathbf z)=\mathbf H_t\mathbf z+\mathbf g_t,\quad \nabla^2 \varphi(\mathbf z)=\mathbf H_t

可以看出 φ(z)\varphi(\mathbf z) 为二次超曲面,若 Ht\mathbf H_t 为正定矩阵,则φ(z)\varphi(\mathbf z)为凸曲面。所以,从几何上说,牛顿法就是用一个二次曲面去拟合你当前所处位置的局部曲面,而梯度下降法是用一个平面去拟合当前的局部曲面,通常情况下,二次曲面的拟合会比平面更好,所以牛顿法选择的下降路径会更符合真实的最优下降路径。

下图为一元函数和二元函数的示意图:

拟牛顿法

牛顿法使用了二阶导数, 其每轮迭代中涉及到 Hessian 矩阵的求逆,计算复杂度相当高,尤其在高维问题中几乎不可行。另外,如果Hessian矩阵不可逆,则这种方法失效。拟牛顿法(quasi-Newton method)则通过正定矩阵近似 Hessian 矩阵的逆矩阵或Hessian矩阵,可显著降低计算开销。

学习率:牛顿法的搜索方向 dt=Ht1gt\mathbf d_t=-\mathbf H_t^{-1}\mathbf g_t ,拟牛顿法同时加入了步长 λt\lambda_t ,迭代公式变为

xt+1=xt+λtdt\mathbf x_{t+1}=\mathbf x_{t}+\lambda_t\mathbf d_t

最优步长可通过线性搜索(line search)取得

λt=argminλf(xt+λtdt)\lambda_t=\arg\min_\lambda f(\mathbf x_{t}+\lambda_t\mathbf d_t)

拟牛顿条件:先看牛顿法中Hessian 矩阵满足的关系。对 f(x)f(\mathbf x)xt+1\mathbf x_{t+1} 处的泰勒二阶展开式求导

f(x)gt+1+Ht+1(xxt+1)\nabla f(\mathbf x)\approx\mathbf g_{t+1}+\mathbf H_{t+1}(\mathbf x-\mathbf x_{t+1})

其中 gt+1,Ht+1\mathbf g_{t+1},\mathbf H_{t+1} 分别是 f(x)f(\mathbf x)xt+1\mathbf x_{t+1} 处的梯度和 Hessian 矩阵。取 x=xt\mathbf x=\mathbf x_{t} 即得

gtgt+1Ht+1(xtxt+1)\mathbf g_{t}-\mathbf g_{t+1}\approx\mathbf H_{t+1}(\mathbf x_{t}-\mathbf x_{t+1})

这就是拟牛顿条件,他对迭代中的 Hessian 矩阵 Ht+1\mathbf H_{t+1} 做约束。为了简便起见,引入记号

yt=gt+1gtst=xt+1xt\mathbf y_t=\mathbf g_{t+1}-\mathbf g_t \\ \mathbf s_t=\mathbf x_{t+1}-\mathbf x_t

根据约束构造一个 Hessian 矩阵的近似 Bt+1Ht+1\mathbf B_{t+1}\approx \mathbf H_{t+1} ,拟牛顿条件简写为

yt=Bt+1st\mathbf y_t=\mathbf B_{t+1}\mathbf s_t

或者构造一个对 Hessian 矩阵的逆的近似 Dt+1Ht+11\mathbf D_{t+1}\approx \mathbf H_{t+1}^{-1},拟牛顿条件也可写为

st=Dt+1yt\mathbf s_t=\mathbf D_{t+1}\mathbf y_t

下降方向:目标函数 f(x)f(\mathbf x) 的一阶泰勒展开近似为

f(x)f(xt)+gtT(xxt)f(\mathbf x)\approx f(\mathbf x_{t})+\mathbf g_t^T(\mathbf x-\mathbf x_t)

x=xt+1\mathbf x=\mathbf x_{t+1} 点带入上式

f(xt+1)=f(xt)λtgtTHt1gtf(\mathbf x_{t+1})= f(\mathbf x_{t})-\lambda_t\mathbf g_t^T\mathbf H_t^{-1}\mathbf g_t

Ht\mathbf H_t 是正定的(Ht1\mathbf H_t^{-1} 也是正定的),那么二次型 gtTHt1gt>0\mathbf g_t^T\mathbf H_t^{-1}\mathbf g_t>0 。当 λt>0\lambda_t>0 时,总有 f(xt+1)<f(xt)f(\mathbf x_{t+1})<f(\mathbf x_{t}) ,也就是说 dt\mathbf d_t 是下降方向。

因此,拟牛顿法的近似矩阵要同时是正定对称矩阵且满足拟牛顿条件。

DFP

DFP 算法(Davidon-Fletcher-Powell algorithm)是由Davidon,Fletcher,Powell三个人的名字的首字母命名的,是最早的拟牛顿法。算法将 Dt+1\mathbf D_{t+1} 作为 Ht+11\mathbf H_{t+1}^{-1} 的近似,假设 Dt+1\mathbf D_{t+1}Dt\mathbf D_{t} 修正得到。迭代公式为

Dt+1=Dt+ΔDt\mathbf D_{t+1}=\mathbf D_t+\Delta\mathbf D_t

其中,初始化常取单位矩阵 D0=I\mathbf D_0=\mathbf I 。相应的拟牛顿条件是该算法的核心是通过迭代的方法,对 Ht+11\mathbf H_{t+1}^{-1} 做近似。相应的拟牛顿条件是

Dt+1yt=st\mathbf D_{t+1}\mathbf y_t=\mathbf s_t

迭代构建近似矩阵的关键是校正矩阵 ΔDt\Delta\mathbf D_t 的构造了,假定有两个附加项构成

ΔDt=auuT+bvvT\Delta\mathbf D_t=a\mathbf{uu}^T+b\mathbf{vv}^T

其中,a,ba,b为待定系数, u,vRn\mathbf u,\mathbf v\in\R^n 是待定向量。从形式上看,这种待定公式至少保证了矩阵 ΔDt\Delta\mathbf D_t 的对称性。根据拟牛顿条件有

st=Dtyt+u(auTyt)+v(bvTyt)=Dtyt+(auTyt)u+(bvTyt)v\begin{aligned} \mathbf s_t&=\mathbf D_{t}\mathbf y_t+\mathbf u(a\mathbf u^T\mathbf y_t)+\mathbf v(b\mathbf v^T\mathbf y_t) \\ &=\mathbf D_{t}\mathbf y_t+(a\mathbf u^T\mathbf y_t)\mathbf u+(b\mathbf v^T\mathbf y_t)\mathbf v \end{aligned}

这里,括号中的auTyta\mathbf u^T\mathbf y_tbvTytb\mathbf v^T\mathbf y_t均为实数,代表了在 u\mathbf uv\mathbf v 方向的拉伸程度,为了计算简单,做如下赋值

auTyt=1,bvTyt=1a\mathbf u^T\mathbf y_t=1,\quad b\mathbf v^T\mathbf y_t=-1

上式代入待定拟牛顿条件便可得

uv=stDtyt\mathbf u-\mathbf v=\mathbf s_t-\mathbf D_t\mathbf y_t

要使上式成立,不妨直接取

u=st,v=Dtyt\mathbf u=\mathbf s_t,\quad\mathbf v=\mathbf D_t\mathbf y_t

继而求得 a,ba,b 的值

a=1stTyt,b=1(Dtyt)Tyt=1ytTDtyta=\frac{1}{\mathbf s_t^T\mathbf y_t},\quad b=-\frac{1}{(\mathbf D_{t}\mathbf y_{t})^T\mathbf y_{t}}=-\frac{1}{\mathbf y_{t}^T\mathbf D_{t}\mathbf y_{t}}

其中,求解 bb 时用到了 Dt\mathbf D_t 的对称性。至此,我们可得到矩阵 Dt+1\mathbf D_{t+1} 的迭代公式

Dt+1=Dt+ststTstTytDtytytTDtytTDtyt\mathbf D_{t+1}=\mathbf D_{t}+\frac{\mathbf s_t\mathbf s_t^T}{\mathbf s_t^T\mathbf y_t}-\frac{\mathbf D_{t}\mathbf y_{t}\mathbf y_{t}^T\mathbf D_{t}}{\mathbf y_{t}^T\mathbf D_{t}\mathbf y_{t}}

可以证明,如果初始矩阵 D0\mathbf D_0 是正定的,则迭代工程中的每个矩阵 Dt\mathbf D_{t} 都是正定的。

BFGS

BFGS 算法(Broyden–Fletcher–Goldfarb–Shanno)是以发明者Broyden、Fletcher、Goldfarb、Shanno 四个人名字的首字母命名。与DFP算法相比,BFGS算法性能更佳,目前已成为求解无约束非线性优化问题最常用的方法之一。由于BFGS法对一维搜索的精度要求不高,并且由迭代产生的BFGS矩阵不易变为奇异矩阵,因而BFGS法比DFP法在计算中具有更好的数值稳定性。

BFGS 算法中的核心公式推导和DFP完全类似。需要注意的是,BFGS 算法直接逼近Hessian矩阵,即将 Bt+1Ht+1\mathbf B_{t+1}\approx \mathbf H_{t+1} 。仍采用迭代方法,迭代公式为

Bt+1=Bt+ΔBt\mathbf B_{t+1}=\mathbf B_t+\Delta\mathbf B_t

其中初始化常取单位矩阵 B0=I\mathbf B_0=\mathbf I 。相应的拟牛顿条件是

Bt+1st=yt\mathbf B_{t+1}\mathbf s_t=\mathbf y_t

关键是校正矩阵 ΔBt\Delta\mathbf B_t 的构造了,和DFP 算法类似,有两个附加项构成

ΔBt=auuT+bvvT\Delta\mathbf B_t=a\mathbf{uu}^T+b\mathbf{vv}^T

其中,a,ba,b为待定系数, u,vRn\mathbf u,\mathbf v\in\R^n 是待定向量。和DFP一样的求解过程,最后得到矩阵 Bt+1\mathbf B_{t+1} 的迭代公式

Bt+1=Bt+ytytTytTstBtststTBtstTBtst\mathbf B_{t+1}=\mathbf B_{t}+\frac{\mathbf y_t\mathbf y_t^T}{\mathbf y_t^T\mathbf s_t}-\frac{\mathbf B_{t}\mathbf s_t\mathbf s_t^T\mathbf B_{t}}{\mathbf s_t^T\mathbf B_{t}\mathbf s_t}

可以证明,如果初始矩阵 B0\mathbf B_0 是正定的,则迭代工程中的每个矩阵 Bt\mathbf B_{t} 都是正定的。

BFGS 算法中牛顿方向 dt=Bt1gt\mathbf d_t=-\mathbf B_t^{-1}\mathbf g_t 通常是求解线性方程组 Btdt=gt\mathbf B_t\mathbf d_t=-\mathbf g_t 来进行。然而更一般的方法是,通过对 Bt+1\mathbf B_{t+1} 的迭代关系应用 Sherman-Morrison 公式,直接给出 Bt+11\mathbf B_{t+1}^{-1}Bt1\mathbf B_{t}^{-1} 之间的关系式。记 Dt+1=Bt+11,Dt=Bt1\mathbf D_{t+1}=B_{t+1}^{-1},\mathbf D_t=\mathbf B_{t}^{-1},则

Dt+1=(IstytTytTst)Dt(IytstTytTst)+ststTytTst\mathbf D_{t+1}= \left(\mathbf I-\frac{\mathbf s_t\mathbf y_t^T}{\mathbf y_t^T\mathbf s_t}\right) \mathbf D_{t} \left(\mathbf I-\frac{\mathbf y_t\mathbf s_t^T}{\mathbf y_t^T\mathbf s_t}\right)+ \frac{\mathbf s_t\mathbf s_t^T}{\mathbf y_t^T\mathbf s_t}

这样 BFGS 中得到了关于 Dt+1\mathbf D_{t+1} 的迭代公式。

Sherman-Morrison 公式:假设 A\mathbf Ann 阶可逆矩阵,u,v\mathbf u,\mathbf vnn 维向量,且 1+vTA1u01+\mathbf v^T\mathbf A^{-1}\mathbf u\neq 0 则有

(A+uvT)1=A1A1uvTA11+vTA1u(\mathbf A+\mathbf{uv}^T)^{-1}=\mathbf A^{-1}-\frac{\mathbf A^{-1}\mathbf {uv}^T\mathbf A^{-1}}{1+\mathbf v^T\mathbf A^{-1}\mathbf u}

Boyden

由 DFP 算法得到的 Dt+1\mathbf D_{t+1} 记作 Dt+1DFP\mathbf D_{t+1}^{\text{DFP}},由 BFGS 算法得到的 Dt+1\mathbf D_{t+1} 记作 Dt+1BFGS\mathbf D_{t+1}^{\text{BFGS}} 。他们都满足拟牛顿条件,所以他们的线性组合

Dt+1=βDt+1DFP+(1β)Dt+1BFGS\mathbf D_{t+1}=\beta\mathbf D_{t+1}^{\text{DFP}}+(1-\beta)\mathbf D_{t+1}^{\text{BFGS}}

也满足拟牛顿条件,而且是正定的,其中 0β10\leqslant\beta\leqslant1 。这样就得到了一类拟牛顿法,称为Boyden族。

L-BFGS

BFGS的迭代公式为

Dt+1=(IstytTytTst)Dt(IytstTytTst)+ststTytTst\mathbf D_{t+1}= \left(\mathbf I-\frac{\mathbf s_t\mathbf y_t^T}{\mathbf y_t^T\mathbf s_t}\right) \mathbf D_{t} \left(\mathbf I-\frac{\mathbf y_t\mathbf s_t^T}{\mathbf y_t^T\mathbf s_t}\right)+ \frac{\mathbf s_t\mathbf s_t^T}{\mathbf y_t^T\mathbf s_t}

其中

yt=gt+1gtst=xt+1xt\mathbf y_t=\mathbf g_{t+1}-\mathbf g_t \\ \mathbf s_t=\mathbf x_{t+1}-\mathbf x_t

在BFGS中,需要用到一个 nn 阶矩阵 Dt\mathbf D_t ,来计算 Dt+1\mathbf D_{t+1} 。当 nn 很大时,存储这个矩阵会超级占用内存。可以看到,整个过程只与 st\mathbf s_tyt\mathbf y_t有关,每次保存这俩个变量,就能恢复每一轮的 Dt\mathbf D_tst\mathbf s_tyt\mathbf y_t 就是中间保存的变量,其都是 n×1n\times 1 向量。相比巨大的 nn 阶矩阵,其一般小很多。

L-BFGS(Limited-memory BFGS 或 Limited-storage BFGS)对BFGS进行了近似,其基本思想是:不再存储完整的矩阵 Dt\mathbf D_t ,而是存储计算过程中的向量序列 {st,yt}\{\mathbf s_t,\mathbf y_t\}。但也不是所有的都存储,而是保留最新的 mm 个(一般在实践中 3m203\leqslant m \leqslant 20 )。每次计算 Dt\mathbf D_t 时,只利用最新的mm个向量序列 。这样一来,存储空间由 O(n2)O(n^2) 降至 O(mn)O(mn)

接下来,讨论 L-BFGS 算法的具体实现过程。在第 tt 次迭代,已求得 xt+1\mathbf x_{t+1} ,记

ρt1=ytTst,Vt=IρtytstT\rho_t^{-1}=\mathbf y_t^T\mathbf s_t,\quad \mathbf V_t=\mathbf I-\rho_t\mathbf y_t\mathbf s_t^T

则迭代公式写为

Dt+1=VtTDtVt+ρtststT\mathbf D_{t+1}=\mathbf V_t^T\mathbf D_{t}\mathbf V_t+\rho_t\mathbf s_t\mathbf s_t^T

如果给定初始矩阵 D0\mathbf D_0 (通常为正定的对角阵,如 D0=I\mathbf D_0=\mathbf I),递推可得到

Dt+1=(VtTVt1TV1TV0T)D0(V0V1Vt1Vt)+(VtTVt1TV2TV1T)(ρ0s0s0T)(V1V2Vt1Vt)+(VtTVt1TV3TV2T)(ρ1s1s1T)(V2V3Vt1Vt)++(VtTVt1T)(ρt2st2st2T)(Vt1Vt)+VtT(ρt1st1st1T)Vt+ρtststT\begin{aligned} \mathbf D_{t+1} &=(\mathbf V_t^T\mathbf V_{t-1}^T\cdots\mathbf V_1^T\mathbf V_0^T)\mathbf D_{0}(\mathbf V_0\mathbf V_1\cdots\mathbf V_{t-1}\mathbf V_t) \\ &+(\mathbf V_t^T\mathbf V_{t-1}^T\cdots\mathbf V_2^T\mathbf V_1^T)(\rho_0\mathbf s_0\mathbf s_0^T)(\mathbf V_1\mathbf V_2\cdots\mathbf V_{t-1}\mathbf V_t) \\ &+(\mathbf V_t^T\mathbf V_{t-1}^T\cdots\mathbf V_3^T\mathbf V_2^T)(\rho_1\mathbf s_1\mathbf s_1^T)(\mathbf V_2\mathbf V_3\cdots\mathbf V_{t-1}\mathbf V_t) \\ &+\cdots \\ &+(\mathbf V_t^T\mathbf V_{t-1}^T)(\rho_{t-2}\mathbf s_{t-2}\mathbf s_{t-2}^T)(\mathbf V_{t-1}\mathbf V_t) \\ &+\mathbf V_t^T(\rho_{t-1}\mathbf s_{t-1}\mathbf s_{t-1}^T)\mathbf V_t \\ &+\rho_t\mathbf s_t\mathbf s_t^T \end{aligned}

由此可见,计算 Dt+1\mathbf D_{t+1} 仅需要使用向量序列 {si}i=0t\{\mathbf s_i\}_{i=0}^t{yi}i=0t\{\mathbf y_i\}_{i=0}^t 。当然计算的过程中多涉及到对角矩阵乘法,都是稀疏矩阵,可以大大优化运算速度。此时便不再需要记录一个巨大的矩阵 Dt\mathbf D_t ,便可以直接推得 Dt+1\mathbf D_{t+1} 了,然而这样做计算量非常之大,因此可以存储近 mm 次的的向量序列 {si,yi}i=tm+1t\{\mathbf s_i,\mathbf y_i\}_{i=t-m+1}^{t}近似来得到 Dt+1\mathbf D_{t+1} 。算法每次迭代均需选择一个初始的矩阵 D0t\mathbf D_0^t

Dt+1=(VtTVt1TVtm+1T)D0t(Vtm+1Vt1Vt)+(VtTVt1TVtm+2T)(ρtm+1stm+1stm+1T)(Vtm+2Vt1Vt)++(VtTVt1T)(ρt2st2st2T)(Vt1Vt)+VtT(ρt1st1st1T)Vt+ρtststT\begin{aligned} \mathbf D_{t+1} &=(\mathbf V_t^T\mathbf V_{t-1}^T\cdots\mathbf V_{t-m+1}^T)\mathbf D_0^t(\mathbf V_{t-m+1}\cdots\mathbf V_{t-1}\mathbf V_t) \\ &+(\mathbf V_t^T\mathbf V_{t-1}^T\cdots\mathbf V_{t-m+2}^T)(\rho_{t-m+1}\mathbf s_{t-m+1}\mathbf s_{t-m+1}^T)(\mathbf V_{t-m+2}\cdots\mathbf V_{t-1}\mathbf V_t) \\ &+\cdots \\ &+(\mathbf V_t^T\mathbf V_{t-1}^T)(\rho_{t-2}\mathbf s_{t-2}\mathbf s_{t-2}^T)(\mathbf V_{t-1}\mathbf V_t) \\ &+\mathbf V_t^T(\rho_{t-1}\mathbf s_{t-1}\mathbf s_{t-1}^T)\mathbf V_t \\ &+\rho_t\mathbf s_t\mathbf s_t^T \end{aligned}

由 BFGS 算法知道 Dt+1\mathbf D_{t+1} 的作用仅用来获取搜索方向 dt+1=Dt+1gt+1\mathbf d_{t+1}=-\mathbf D_{t+1}\mathbf g_{t+1} 。L-BFGS 给出了一种快速算法,使用双循环迭代(two-loop recursion)计算 Dt+1gt+1\mathbf D_{t+1}\mathbf g_{t+1} 。该算法的计算过程如下

算法L-BFGS的步骤如下所示

坐标下降法

坐标下降法(Coordinate Descent)求解流程为每次选择一个分量 xjx_j 进行优化,将其他分量固定住不动,这样将一个多元函数的极值问题转换为一元函数的极值问题。如果要求解的问题规模很大,这种做法能有效的加快速度。

假设要求解的优化问题为

minxf(x1,x2,,xn)\min\limits_{\mathbf x} f(x_1,x_2,\cdots,x_n)

那么第 tt 轮第 jj 个变量得迭代公式为

xj(t)=argminyf(x1(t),,xj1(t),y,xj+1(t1),xn(t1))x_j^{(t)}=\arg\min_{y} f(x_1^{(t)},\cdots,x_{j-1}^{(t)},y,x_{j+1}^{(t-1)}\cdots,x_n^{(t-1)})

与梯度下降法类似,通过迭代执行该过程,能收敛到所期望的极值。坐标下降法不需计算目标函数的梯度,在每步迭代中仅需求解一维搜索问题,对于某些复杂问题计算较为简便。


  1. 泰勒展开式 f(x+Δx)=f(x)+f(x)Δx+12f(x)(Δx)2+f(x+\Delta x)=f(x)+f'(x)\Delta x+\dfrac{1}{2}f''(x)(\Delta x)^2+\cdots ↩︎

  2. Hadamard 积:设 A=(aij)m×n\mathbf A=(a_{ij})_{m\times n}B=(bij)m×n\mathbf B=(b_{ij})_{m\times n} 是两个同型矩阵,则称各元素乘积组成的同型矩阵 AB=(aijbij)m×n\mathbf A\odot\mathbf B=(a_{ij}b_{ij})_{m\times n}A,B\mathbf A,\mathbf B 的哈达玛积(Hadamard product)。 ↩︎ ↩︎

  3. 从经验上看,学习率在一开始要保持大些来保证收敛速度,在收敛到最优点附近时要小些以避免来回振荡.比较简单的学习率调整可以通过学习率衰减(Learning Rate Decay)的方式来实现,也称为学习率退火(Learning Rate Annealing)。 ↩︎

  4. 范数:向量 xRn\mathbf x\in\R^n 的 Lp 范数定义如下 xp=(x1p+x2p++xnp)1/p\|\mathbf x\|_p=(|x_1|^p+|x_2|^p+\cdots+|x_n|^p)^{1/p}。L2 范数在机器学习中出现地十分频繁,经常简化表示为x\|\mathbf x\|,略去下标2。 LL_\infty 范数也被称为最大范数,这个范数表示向量中具有最大幅值的元素的绝对值 x=max1inxi\|\mathbf x\|_\infty=\max\limits_{1\leqslant i\leqslant n}|x_i| ↩︎