1 梯度下降
1.1 特点
求解无约束最优化问题的一种最常用方法
实现简单
是一种迭代算法,每一步需要求解目标函数的梯度向量
1.2 思想
假设
f
(
x
)
f(x)
f(x)是
R
n
R^n
Rn上具有一阶连续偏导数的函数。要求解得无约束最优化问题是
(1.1)
min
x
∈
R
n
f
(
x
)
\tag{1.1}\min_{x\in R^n} f(x)
x∈Rnminf(x)(1.1)
选取适当的初始值
x
(
0
)
x^{(0)}
x(0) ,不断迭代,更新
x
x
x的值,进行目标函数的最小化,直到收敛。由于负梯度方向是函数下降最快的方向,因此,在迭代的每一步,以负梯度方向更新
x
x
x的值,从而达到减少函数值的目的。
1.3 数学基础
由于
f
(
x
)
f(x)
f(x)具有一阶连续偏导数,若第
k
k
k次迭代值为
x
(
k
)
x^{(k)}
x(k),则可将
f
(
x
)
f(x)
f(x)在
x
(
k
)
x^{(k)}
x(k)附近进行一阶泰勒展开:
(1.2)
f
(
x
)
=
f
(
x
(
k
)
)
+
g
k
T
(
x
−
x
(
k
)
)
\tag{1.2}f(x)=f(x^{(k)})+g_k^T(x-x^{(k)})
f(x)=f(x(k))+gkT(x−x(k))(1.2)
这里,
g
k
=
g
(
x
(
k
)
)
=
∇
f
(
x
(
k
)
)
g_k=g(x^{(k)})=\nabla f(x^{(k)})
gk=g(x(k))=∇f(x(k)) 为
f
(
x
)
f(x)
f(x)在
x
(
k
+
1
)
x^{(k+1)}
x(k+1)的梯度。
求出第
k
+
1
k+1
k+1次迭代值
x
(
k
+
1
)
x^{(k+1)}
x(k+1):
(1.3)
x
(
k
+
1
)
←
x
(
k
)
+
λ
k
p
k
\tag{1.3}x^{(k+1)} \leftarrow x^{(k)}+\lambda_{k}p_{k}
x(k+1)←x(k)+λkpk(1.3)
其中,
p
k
p_{k}
pk是搜索方向,取负梯度方向
p
k
=
−
∇
f
(
x
(
k
)
)
p_{k}=-\nabla{f(x^{(k)})}
pk=−∇f(x(k)),
λ
k
\lambda_{k}
λk是步长,由一维搜索确定,即
λ
k
\lambda_{k}
λk使得
(1.4) f ( x ( k ) + λ k p k ) = min λ ≥ 0 x ( k ) + λ p k \tag{1.4}f(x^{(k)}+\lambda_{k}p_{k})=\min_{\lambda \geq 0}{x^{(k)}+\lambda p_{k}} f(x(k)+λkpk)=λ≥0minx(k)+λpk(1.4)
1.4 具体算法
输入:目标函数 f ( x ) f(x) f(x), 梯度函数 g ( X ) = ∇ f ( x ) g(X)=\nabla f(x) g(X)=∇f(x),计算精度 ϵ \epsilon ϵ。
输出: f ( x ) f(x) f(x)的极小值点 x ∗ x^* x∗。
- 取初始值 x ( 0 ) ∈ R n x^{(0)}\in R^n x(0)∈Rn,置 k = 0 k=0 k=0
- 计算 f ( x ( k ) ) f(x^{(k)}) f(x(k))
- 计算梯度 g k = g ( x ( k ) ) g_k=g(x^{(k)}) gk=g(x(k)),当 ∣ ∣ g k ∣ ∣ < ϵ ||g_k||<\epsilon ∣∣gk∣∣<ϵ时,停止迭代,令x*=x{(k)};否则,令 p k = − g ( x ( k ) ) p_k=-g(x^{(k)}) pk=−g(x(k)),求 λ k \lambda_k λk,使
f ( x ( k ) + λ k p k ) = min λ ≥ 0 f ( x ( k ) + λ p k ) f(x^{(k)}+\lambda_k p_k)=\min_{\lambda \geq0 }f(x^{(k)}+\lambda p_k) f(x(k)+λkpk)=λ≥0minf(x(k)+λpk)- 置 x ( k + 1 ) = x ( k ) + λ k p k x^{(k+1)}=x^{(k)}+\lambda_k p_k x(k+1)=x(k)+λkpk,计算 f ( x ( k + 1 ) ) f(x^{(k+1)}) f(x(k+1))
当 ∣ ∣ f ( x ( k + 1 ) ) − f ( x ( k ) ) ∣ ∣ < ϵ ||f(x^{(k+1)})-f(x^{(k)})||<\epsilon ∣∣f(x(k+1))−f(x(k))∣∣<ϵ或 ∣ ∣ x ( k + 1 ) − x ( k ) ∣ ∣ < ϵ ||x^(k+1)-x^(k)||<\epsilon ∣∣x(k+1)−x(k)∣∣<ϵ时,停止迭代,令 x ∗ = x ( k + 1 ) x^*=x^(k+1) x∗=x(k+1)- 否则,置 k = k + 1 k=k+1 k=k+1,转 3
当目标函数是凸函数时,梯度下降法的解释全局最优解。一般情况下,其解不保证是全局最优解。梯度下降法的收敛速度也未必是最快的。
2 牛顿法和拟牛顿法
2.1 特点
也是求解无约束最优化问题的常用方法
有收敛速度快的特点
是迭代算法,每一步需要求解目标函数的海赛矩阵的逆矩阵
计算比较复杂
拟牛顿法通过正定矩阵近似海赛矩阵的逆矩阵或海赛矩阵,简化了这一计算过程
2.2 牛顿法
2.2.1 数学基础
考虑无约束最优化问题
(2.1)
min
x
∈
R
n
f
(
x
)
\min_{x \in R^n} f(x) \tag{2.1}
x∈Rnminf(x)(2.1)
其中
x
∗
x^*
x∗为目标函数的最小点
假设
f
(
x
)
f(x)
f(x)具有二阶连续偏导数,若第
k
k
k次迭代值为
x
(
k
)
x^{(k)}
x(k),则可将
f
(
x
)
f(x)
f(x)在
x
(
k
)
x^{(k)}
x(k)附近进行二阶泰勒展开
(2.2)
f
(
x
)
=
f
(
x
(
k
)
)
+
g
k
T
(
x
−
x
(
k
)
)
+
1
2
(
x
−
x
(
k
)
)
T
H
(
x
(
k
)
)
(
x
−
x
(
k
)
)
f(x)=f(x^{(k)})+g_k^T(x-x^{(k)})+ {1\over2}(x-x^{(k)})^T H(x^{(k)})(x-x^{(k)}) \tag{2.2}
f(x)=f(x(k))+gkT(x−x(k))+21(x−x(k))TH(x(k))(x−x(k))(2.2)
这里,
g
k
=
g
(
x
(
k
)
)
=
∇
f
(
x
(
k
)
)
g_k=g(x^{(k)})=\nabla f(x^{(k)})
gk=g(x(k))=∇f(x(k))是
f
(
x
)
f(x)
f(x)的梯度向量在点
x
(
k
)
x^{(k)}
x(k)的值,
H
(
x
(
k
)
)
H(x^{(k)})
H(x(k))是
f
(
x
)
f(x)
f(x)的海赛矩阵(Hesse matrix)
(2.3)
H
(
X
)
=
[
∂
2
f
∂
x
i
∂
x
j
]
n
×
n
H(X)= \left[ \partial^2 f \over \partial x_i \partial x_j \right]_{n \times n} \tag{2.3}
H(X)=[∂xi∂xj∂2f]n×n(2.3)
在点
x
(
k
)
x^{(k)}
x(k)的值。函数
f
(
x
)
f(x)
f(x)有极值的必要条件是在极值点处一阶导数为0,即梯度向量为0。特别是当
H
(
X
(
k
)
)
H(X^{(k)})
H(X(k))是正定矩阵时,函数
f
(
x
)
f(x)
f(x)的极值为极小值。
2.2.2 思想
牛顿法利用极小点的必要条件
(2.4)
∇
f
(
x
)
=
0
\nabla f(x)=0 \tag{2.4}
∇f(x)=0(2.4)
每次迭代中从点
x
(
k
)
x^{(k)}
x(k)开始,求目标函数的极小点,作为第
k
+
1
k+1
k+1次迭代值
x
(
k
+
1
)
x^{(k+1)}
x(k+1)。具体地,假设
x
(
k
+
1
)
x^{(k+1)}
x(k+1)满足:
(2.5)
∇
f
(
x
(
k
+
1
)
)
=
0
\nabla f(x^{(k+1)})=0 \tag{2.5}
∇f(x(k+1))=0(2.5)
由式(2.2)有
(2.6)
∇
f
(
x
)
=
g
k
+
H
k
(
x
−
x
(
k
)
)
\nabla f(x)=g_k + H_k(x-x^{(k)}) \tag{2.6}
∇f(x)=gk+Hk(x−x(k))(2.6)
其中
H
k
=
H
(
x
(
k
)
)
H_k=H(x^{(k)})
Hk=H(x(k))。这样,式(2.5)变为
(2.7)
g
k
+
H
k
(
x
(
k
+
1
)
−
x
(
k
)
)
=
0
g_k + H_k(x^{(k+1)}-x^{(k)})=0 \tag{2.7}
gk+Hk(x(k+1)−x(k))=0(2.7)
因此,
(2.8)
x
(
k
+
1
)
=
x
(
k
)
−
H
k
−
1
g
k
x^{(k+1)}=x^{(k)}-H_k^{-1}g_k \tag{2.8}
x(k+1)=x(k)−Hk−1gk(2.8)
或者
(2.9)
x
(
k
+
1
)
=
x
(
k
)
+
p
k
x^{(k+1)}=x^{(k)}+p_k \tag{2.9}
x(k+1)=x(k)+pk(2.9)
其中,
(2.10)
H
k
p
k
=
−
g
k
H_kp_k=-g_k \tag{2.10}
Hkpk=−gk(2.10)
用式(2.8)作为迭代公式的算法就是牛顿法。
2.2.3 具体算法
输入:目标函数 f ( X ) f(X) f(X),梯度 g ( x ) = ∇ f ( x ) g(x)=\nabla f(x) g(x)=∇f(x),海赛矩阵 H ( x ) H(x) H(x),精度要求 ϵ \epsilon ϵ;
输出: f ( x ) f(x) f(x)的极小点 x ∗ x^* x∗
- 取初始点 x ( 0 ) x^{(0)} x(0),置 k = 0 k=0 k=0
- 计算 g k = g ( x ( k ) ) g_k=g(x^{(k)}) gk=g(x(k))
- 若 ∣ ∣ g k ∣ ∣ < ϵ ||g_k||<\epsilon ∣∣gk∣∣<ϵ,则停止计算,得近似解 x ∗ = x ( k ) x^*=x^{(k)} x∗=x(k)
- 计算 H k = H ( x ( k ) ) H_k=H(x^{(k)}) Hk=H(x(k)),并求 p k p_k pk:
H k p k = − g k H_kp_k=-g_k Hkpk=−gk- 置 x ( k + 1 ) = x ( k ) + p k x^{(k+1)}=x^{(k)}+p_k x(k+1)=x(k)+pk
- 置 k = k + 1 k=k+1 k=k+1,转(2)
步骤(4)求 p k p_k pk, p k = − H k − 1 g k p_k=-H_k^{-1}g_k pk=−Hk−1gk,要求 H k − 1 H_k^{-1} Hk−1,计算比较复杂,所以有其他改进的方法。
2.3 拟牛顿法
2.3.1 数学基础
在牛顿法的迭代中,需要计算海赛矩阵的逆矩阵
H
−
1
H^{-1}
H−1,这一计算比较复杂,考虑用一个
n
n
n阶矩阵
G
k
=
G
(
x
(
k
)
)
G_k=G(x^{(k)})
Gk=G(x(k))来近似替代
H
k
−
1
=
H
−
1
(
x
(
k
)
)
H^{-1}_k=H^{-1}(x^{(k)})
Hk−1=H−1(x(k))。这就是拟牛顿法的基本想法。
先看牛顿法迭代中海赛矩阵
H
k
H_k
Hk满足的条件。首先,
H
k
H_k
Hk满足以下关系。在式(2.6)中取
x
=
x
(
k
+
1
)
x=x^{(k+1)}
x=x(k+1),即得
(2.11)
g
k
+
1
−
g
k
=
H
k
(
x
(
k
+
1
)
−
x
(
k
)
)
g_{k+1}-g_k=H_k({x^{(k+1)}-x^{(k)}}) \tag{2.11}
gk+1−gk=Hk(x(k+1)−x(k))(2.11)
记
y
k
=
g
k
+
1
−
g
k
y_k=g_{k+1}-g_{k}
yk=gk+1−gk,
δ
k
=
x
(
k
+
1
)
−
x
(
k
)
\delta_k=x^{(k+1)}-x^{(k)}
δk=x(k+1)−x(k),则
(2.12)
y
k
=
H
k
δ
k
y_k=H_k\delta_k \tag{2.12}
yk=Hkδk(2.12)
或
(2.13)
H
k
−
1
y
k
=
δ
k
H_k^{-1}y_k=\delta_k \tag{2.13}
Hk−1yk=δk(2.13)
式(2.12)或(2.13)称为拟牛顿条件。
如果
H
k
H_k
Hk是正定的(
H
k
−
1
H_k^{-1}
Hk−1也是正定的),那么可以保证牛顿法搜索方向
p
k
p_k
pk是下降方向,这是因为搜索方向
p
k
=
−
H
k
−
1
g
k
p_k=-H_k^{-1}g_k
pk=−Hk−1gk,由式(2.8)有
(2.14)
x
=
x
k
+
λ
p
k
=
x
(
k
)
=
x
(
k
)
−
λ
H
(
−
1
)
g
k
x=x^{k}+\lambda p_k=x^{(k)}=x^{(k)}-\lambda H^{(-1)}g_k \tag{2.14}
x=xk+λpk=x(k)=x(k)−λH(−1)gk(2.14)
所以
f
(
x
)
f(x)
f(x)在
x
(
k
)
x^{(k)}
x(k)的泰勒展开式(2.2)可以近似写成:
(2.15)
f
(
x
)
=
f
(
x
(
k
)
)
−
λ
g
k
T
H
k
−
1
g
k
f(x)=f(x^{(k)})-\lambda g_k^TH_k^{-1}g_k \tag{2.15}
f(x)=f(x(k))−λgkTHk−1gk(2.15)
因
H
k
−
1
H_k^{-1}
Hk−1正定,固有
g
k
T
H
k
−
1
g
k
>
0
g_k^TH_k^{-1}g_k>0
gkTHk−1gk>0。当
λ
\lambda
λ为一个充分小的正数时,总有
f
(
x
)
<
f
(
x
(
k
)
)
f(x)<f(x^{(k)})
f(x)<f(x(k)),也就是说
p
k
p_k
pk是下降方向。
2.3.2 思想
拟牛顿法将
G
k
G_k
Gk作为
H
k
−
1
H_k^{-1}
Hk−1的近似,要求矩阵
G
k
G_k
Gk满足同样的条件。首先,每次迭代矩阵
G
k
G_k
Gk是正定的。同时,
G
k
G_k
Gk满足下面的拟牛顿条件:
(2.16)
G
k
+
1
y
k
=
δ
k
G_{k+1}y_k=\delta_k \tag{2.16}
Gk+1yk=δk(2.16)
按照拟牛顿条件选择
G
k
G_k
Gk作为
H
k
−
1
H_k^{-1}
Hk−1的近似或选择
B
k
B_k
Bk作为
H
k
H_k
Hk近似的算法称为拟牛顿法。
2.3.4 具体算法
按照拟牛顿条件,在每次迭代中可以选择更新矩阵
G
k
+
1
G_{k+1}
Gk+1:
(2.17)
G
k
+
1
=
G
k
+
Δ
G
k
G_{k+1}=G_k+\Delta G_k \tag{2.17}
Gk+1=Gk+ΔGk(2.17)
这种选择有一定的灵活性,因此有多重具体实现方法。2.6节介绍Broyden类牛顿法
2.4 DFP算法
2.4.1 数学基础&&思想
DFP算法选择
G
k
+
1
G_k+1
Gk+1的方法是,假设每一步迭代中矩阵
G
k
+
1
G_{k+1}
Gk+1是由
G
k
G_k
Gk加上两个附加项构成的,即
(2.18)
G
k
+
1
=
G
k
+
P
k
+
Q
k
G_{k+1}=G_k+P_k+Q_k \tag{2.18}
Gk+1=Gk+Pk+Qk(2.18)
其中
P
k
P_k
Pk,
Q
k
Q_k
Qk是待定矩阵,这时,
(2.19)
G
k
+
1
y
k
=
G
k
y
k
+
P
k
y
k
+
Q
k
y
k
G_{k+1}y_k=G_ky_k+P_ky_k+Q_ky_k \tag{2.19}
Gk+1yk=Gkyk+Pkyk+Qkyk(2.19)
为使
G
k
+
1
G_{k+1}
Gk+1满足拟牛顿条件,可使
P
k
P_k
Pk和
Q
k
Q_k
Qk满足:
(2.20)
P
k
y
k
=
δ
k
P_ky_k=\delta_k \tag{2.20}
Pkyk=δk(2.20)
(2.21)
Q
k
y
k
=
−
G
k
y
k
Q_ky_k=-G_ky_k \tag{2.21}
Qkyk=−Gkyk(2.21)
事实上,不难找出这样的P_k和Q_k,例如取
(2.22)
P
k
=
δ
k
δ
k
T
δ
k
T
y
k
P_k= \frac{\delta_k\delta_k^T}{\delta_k^Ty_k} \tag{2.22}
Pk=δkTykδkδkT(2.22)
(2.23)
Q
k
=
−
G
k
y
k
y
k
T
G
k
y
k
T
G
k
y
k
Q_k=-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} \tag{2.23}
Qk=−ykTGkykGkykykTGk(2.23)
这样就可得到矩阵
G
k
+
1
G_{k+1}
Gk+1的迭代公式:
(2.24)
G
k
+
1
=
G
k
+
δ
k
δ
k
T
δ
k
T
y
k
−
G
k
y
k
y
k
T
G
k
y
k
T
G
k
y
k
G_{k+1}=G_k+\frac{\delta_k\delta_k^T}{\delta_k^Ty_k}-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} \tag{2.24}
Gk+1=Gk+δkTykδkδkT−ykTGkykGkykykTGk(2.24)
称为DFP算法。
可以证明,如果初始矩阵
G
0
G_0
G0是正定的,则迭代过程中的每个矩阵
G
k
G_k
Gk都是正定的。
2.4.2 具体算法
输入:目标函数 f ( x ) f(x) f(x),梯度 g ( x ) = ∇ f ( x ) g(x)=\nabla f(x) g(x)=∇f(x),精度要求 ϵ \epsilon ϵ;
输出: f ( x ) f(x) f(x)的极小点 x ∗ x^* x∗。
- 选定初始点 x ( 0 ) x^{(0)} x(0),取 G 0 G_0 G0为正定对称矩阵,置 k = 0 k=0 k=0
- 计算 g k = g ( x ( k ) ) g_k=g(x^{(k)}) gk=g(x(k))。若 ∣ ∣ g k ∣ ∣ < ϵ ||g_k||<\epsilon ∣∣gk∣∣<ϵ,则停止计算,得近似解 x ∗ = x ( k ) x^*=x^{(k)} x∗=x(k);否则转 3
- 置 p k = − G k g k p_k=-G_kg_k pk=−Gkgk
- 一维搜索:求 λ \lambda λ使得
f ( x ( k ) + λ k p k ) = min λ ≥ 0 f ( x ( k ) + λ p k ) f(x^{(k)}+\lambda_kp_k)=\min_{\lambda\geq0}f(x^{(k)}+\lambda p_k) f(x(k)+λkpk)=λ≥0minf(x(k)+λpk)- 置 x ( k + 1 ) = x ( k ) + λ k p k x^{(k+1)}=x^{(k)}+\lambda_kp_k x(k+1)=x(k)+λkpk
- 计算 g k + 1 = g ( x ( k + 1 ) ) g_{k+1}=g(x^{(k+1)}) gk+1=g(x(k+1)),若 ∣ ∣ g k + 1 ∣ ∣ < ϵ ||g_{k+1}||<\epsilon ∣∣gk+1∣∣<ϵ,则停止计算,得近似解 x ∗ = x ( k + 1 ) x^*=x^{(k+1)} x∗=x(k+1);否则,按式(2.24)算出 G k + 1 G_{k+1} Gk+1
- 置 k = k + 1 k=k+1 k=k+1,转 3
2.5 BFGS算法
2.5.1 特点
BFGS算法是最流行的拟牛顿算法。
2.5.2 数学基础
可以考虑用 G k G_k Gk逼近海赛矩阵的逆矩阵 H − 1 H_{-1} H−1,也可以考虑用 B k B_k Bk逼近海赛矩阵 H H H。
2.5.3 思想
这时,相应的拟牛顿条件是
(2.25)
B
k
+
1
δ
k
=
y
k
B_{k+1}\delta_k=y_k \tag{2.25}
Bk+1δk=yk(2.25)
可以用同样的方法得到另一迭代公式。首先令
(2.26)
B
k
+
1
=
B
k
+
P
k
+
Q
k
B_{k+1}=B_k+P_k+Q_k \tag{2.26}
Bk+1=Bk+Pk+Qk(2.26)
(2.27)
B
k
+
1
δ
k
=
B
k
δ
k
+
P
k
δ
k
+
Q
k
δ
k
B_{k+1}\delta_k=B_k\delta_k+P_k\delta_k+Q_k\delta_k \tag{2.27}
Bk+1δk=Bkδk+Pkδk+Qkδk(2.27)
考虑使
P
k
P_k
Pk和
Q
k
Q_k
Qk满足:
(2.28)
P
k
δ
k
=
y
k
P_k\delta_k=y_k \tag{2.28}
Pkδk=yk(2.28)
(2.29)
Q
k
δ
k
=
−
B
k
δ
k
Q_k\delta_k=-B_k\delta_k \tag{2.29}
Qkδk=−Bkδk(2.29)
找出适合条件的
P
k
P_k
Pk和
Q
k
Q_k
Qk,得到
B
F
G
S
BFGS
BFGS算法矩阵
B
k
+
1
B_{k+1}
Bk+1的迭代公式:
(2.30)
B
k
+
1
=
B
k
+
y
k
y
k
T
y
k
T
δ
k
−
B
k
δ
k
δ
k
T
B
k
δ
k
T
B
k
δ
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} \tag{2.30}
Bk+1=Bk+ykTδkykykT−δkTBkδkBkδkδkTBk(2.30)
可以证明,如果初始矩阵
B
0
B_0
B0是正定的,则迭代过程中的每个矩阵
B
k
B_k
Bk都是正定的。
2.5.4 具体算法
输入:目标函数 f ( x ) f(x) f(x),梯度 g ( x ) = ∇ f ( x ) g(x)=\nabla f(x) g(x)=∇f(x),精度要求 ϵ \epsilon ϵ;
输出: f ( x ) f(x) f(x)的极小点 x ∗ x^* x∗。
- 选定初始点 x ( 0 ) x^{(0)} x(0),取 B 0 B_0 B0为正定对称矩阵,置 k = 0 k=0 k=0
- 计算 g k = g ( x ( k ) ) g_k=g(x^{(k)}) gk=g(x(k))。若 ∣ ∣ g k ∣ ∣ < ϵ ||g_k||<\epsilon ∣∣gk∣∣<ϵ,则停止计算,得近似解 x ∗ = x ( k ) x^*=x^{(k)} x∗=x(k);否则转 3
- 由 B K p k = − g k B_Kp_k=-g_k BKpk=−gk求出 p k p_k pk
- 一维搜索:求 λ k \lambda_k λk使得
f ( x ( k ) + λ k + p k ) = min λ ≥ 0 f ( x ( k ) + λ p k ) f(x^{(k)}+\lambda_k+p_k)=\min_{\lambda\geq0}f(x^{(k)}+\lambda p_k) f(x(k)+λk+pk)=λ≥0minf(x(k)+λpk)- 置 x ( k + 1 ) = x ( k ) + λ k p k x^{(k+1)}=x^{(k)}+\lambda_kp_k x(k+1)=x(k)+λkpk
- 计算 g k + 1 = g ( x ) ( k + 1 ) g_{k+1}=g(x)^{(k+1)} gk+1=g(x)(k+1),若 ∣ ∣ g k + 1 ∣ ∣ < ϵ ||g_{k+1}||<\epsilon ∣∣gk+1∣∣<ϵ,则停止计算,得近似解 x ∗ = x ( k + 1 ) x^*=x^{(k+1)} x∗=x(k+1);否则,按式(2.30)算出 B k + 1 B_{k+1} Bk+1
- 置 k = k + 1 k=k+1 k=k+1,转 3
2.6 Broyden类牛顿法
2.6.1 数学基础——Sherman-Morrison公式
假设
A
A
A是
n
n
n阶可逆矩阵,
u
u
u,
v
v
v是
n
n
n维向量,且
A
+
u
v
T
A+uv^T
A+uvT也是可逆矩阵,则
(
A
+
u
v
)
−
1
=
A
−
1
−
A
−
1
u
v
T
A
−
1
1
+
v
T
A
−
1
u
(A+uv)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}
(A+uv)−1=A−1−1+vTA−1uA−1uvTA−1
2.6.2 思想
我们可以从 B F G S BFGS BFGS算法矩阵 B k B_k Bk的迭代公式(2.30)得到 B F G S BFGS BFGS算法关于 G k G_k Gk的迭代公式。
2.6.3 具体算法
事实上,若记
G
k
=
B
k
−
1
G_k=B_k^{-1}
Gk=Bk−1,
G
k
+
1
=
B
k
+
1
−
1
G_{k+1}=B_{k+1}^{-1}
Gk+1=Bk+1−1,那么对式(2.30)两次应用
S
h
e
r
m
a
n
−
M
o
r
r
i
s
o
n
Sherman-Morrison
Sherman−Morrison公式即得
(2.31)
G
k
+
1
=
(
I
−
δ
k
y
k
T
δ
k
T
y
k
)
G
k
(
I
−
δ
k
y
k
T
δ
k
T
y
k
)
T
+
δ
k
δ
k
T
δ
k
T
y
k
G_{k+1}=\left( I-\frac{\delta_ky_k^T}{\delta_k^Ty_k} \right) G_k \left( I-\frac{\delta_ky_k^T}{\delta_k^Ty_k} \right)^T+\frac{\delta_k\delta_k^T}{\delta_k^Ty_k} \tag{2.31}
Gk+1=(I−δkTykδkykT)Gk(I−δkTykδkykT)T+δkTykδkδkT(2.31)
称为
B
F
G
S
BFGS
BFGS算法关于
G
k
G_k
Gk的迭代公式。
由
D
F
P
DFP
DFP算法
G
k
G_k
Gk的迭代公式(2.23)得到的
G
k
+
1
G_{k+1}
Gk+1记作
G
D
F
P
G_{DFP}
GDFP,由
B
F
G
S
BFGS
BFGS算法
G
k
G_k
Gk的迭代公式(2.31)得到的
G
k
+
1
G_{k+1}
Gk+1记作
G
B
F
G
S
G_{BFGS}
GBFGS,它们都满足方程拟牛顿条件时,所以它们的线性组合
(2.32)
G
k
+
1
=
α
G
D
F
P
+
(
1
−
α
G
B
F
G
S
)
G_{k+1}=\alpha G^{DFP}+(1-\alpha G^{BFGS}) \tag{2.32}
Gk+1=αGDFP+(1−αGBFGS)(2.32)
也满足拟牛顿条件式,而且是正定的。其中
0
≤
α
≤
1
0\leq\alpha\leq1
0≤α≤1。这样就得到了一类拟牛顿法,称为
B
r
o
y
d
e
n
Broyden
Broyden类算法。
3 随机梯度下降
3.1 特点
随机梯度下降( S G D SGD SGD)及其变种很可能是一般机器学习中应用最多的优化算法,特别是在深度学习中。随机梯度下降可以很大程度地加速,沿着随机挑选的小批量数据的梯度下降方向。
3.2 思想
按照数据生成分布抽取 m m m个小批量(独立同分布)样本,通过计算它们梯度均值,我们可以得到梯度的无偏估计。
3.3 具体算法
当前:随机梯度下降( S G D SGD SGD)在第k个训练迭代的更新。
输入:学习率 ϵ k \epsilon_k ϵk, 初始参数 θ \theta θ
\quad while 停止条件不达到 do
\quad \quad 从训练集中采集包含 m m m个样本{ x ( 1 ) , . . . , x ( m ) x^{(1)},...,x^{(m)} x(1),...,x(m)}的小批量,其中 x ( i ) x^{(i)} x(i)对应目标为 y ( i ) y^{(i)} y(i)。
\quad \quad 计算梯度估计: g ← + 1 m ∇ θ Σ i L ( f ( x ( i ) ; θ ) , y ( i ) ) g \leftarrow + \frac{1}{m}\nabla_\theta\Sigma_iL(f(x^{(i)};\theta),y^{(i)}) g←+m1∇θΣiL(f(x(i);θ),y(i))
\quad \quad 应用更新: θ ← θ − ϵ g \theta \leftarrow \theta - \epsilon g θ←θ−ϵg
\quad end while
3.4 讲解
S
G
D
SGD
SGD算法中的一个关键参数是学习率。之前,我们介绍的SGD使用固定的学习率。在实践中,有必要随着时间的推移逐渐降低学习率,因此我们将第
k
k
k步迭代的学习率记作
ϵ
k
\epsilon_k
ϵk。
这是因为
S
G
D
SGD
SGD中梯度估计引入的噪声源(
m
m
m个训练样本的随机采样)并不会在极小点处消失。相比之下,当我们使用批量梯度下降达到极小点时,整个代价函数的真实梯度会变得很小,之后为0,因此批量梯度下降可以使用固定的学习率。保证
S
G
D
SGD
SGD收敛的一个充分条件是
(3.1)
Σ
k
=
1
∞
ϵ
k
=
∞
\Sigma_{k=1}^{\infty}\epsilon_k=\infty \tag{3.1}
Σk=1∞ϵk=∞(3.1)
且
(3.2)
Σ
k
=
1
∞
ϵ
k
2
<
∞
\Sigma_{k=1}^{\infty}\epsilon_k^2<\infty \tag{3.2}
Σk=1∞ϵk2<∞(3.2)
实践中,一般会线性衰减学习率直到第
r
r
r次迭代:
(3.3)
ϵ
k
=
(
1
−
α
)
ϵ
0
+
α
ϵ
r
\epsilon_k=(1-\alpha)\epsilon_0+\alpha\epsilon_r \tag{3.3}
ϵk=(1−α)ϵ0+αϵr(3.3)
其中
α
=
k
r
\alpha=\frac{k}{r}
α=rk。在
r
r
r步迭代之后,一般使
ϵ
\epsilon
ϵ保持常数。
学习率可通过实验和误差来选取,通常最好的选择方法是检测目标函数值随时间变化的学习曲线。与其说是科学,这更像是一门艺术,我们应该谨慎地参考关于这个问题的大部分指导。使用线性策略时,通常选择的参数为
ϵ
0
,
ϵ
r
,
和
r
\epsilon_0,\epsilon_r,和r
ϵ0,ϵr,和r。通常
r
r
r被设为需要反复遍历训练集几百次的迭代次数。通常
ϵ
r
\epsilon_r
ϵr应设为大约
ϵ
0
的
1
%
\epsilon_0的1\%
ϵ0的1%。主要问题是如何设置
ϵ
0
\epsilon_0
ϵ0。若
ϵ
0
\epsilon_0
ϵ0太大,学习曲线将会剧烈振荡,代价函数值通常会明显增加。温和的振荡是良好的,容易在训练随机代价函数(例如使用
D
r
o
p
o
u
t
Dropout
Dropout的代价函数)时出现。如果学习率太小,那么学习过程会很缓慢。如果初始学习率太低,那么学习可能会卡在一个相当高的代价值。通常,就总训练时间和最终代价值而言,最优初始学习会高于大约迭代100次后达到最佳效果的学习率。因此,通常最好是检测最早的几轮迭代,选择一个比在效果上表现最佳的学习率更大的学习率,但又不能太大导致严重的震荡。
S
G
D
SGD
SGD及相关的小批量亦或是更广义的基于梯度优化的在线学习算法,一个重要的性质是每一步更新的计算时间不依赖训练样本数目的多寡。即使训练样本数目非常大时,它们也能收敛。对于足够大的数据集,
S
G
D
SGD
SGD可能会在处理整个训练集之前就收敛到最终测试集误差的某个固定容差范围内。
研究优化算法的收敛率,一般会衡量额外误差(excess error)
J
(
θ
)
−
min
θ
J
(
θ
)
J(\theta)-\min_\theta J(\theta)
J(θ)−minθJ(θ),即当前代价函数超出最低可能代价的量。
S
G
D
SGD
SGD应用于凸问题时,
k
k
k次迭代后的额外误差量级是
O
(
1
k
)
O\left(\frac{1}{\sqrt{k}} \right)
O(k1),在强凸情况下是
O
(
1
k
)
O(\frac{1}{k})
O(k1)。除非假定额外的条件,否则这些界限不能进一步改进。批量梯度下降在理论上比随机梯度下降有更好的收敛率。然而,
C
r
a
m
e
r
−
R
a
o
Cramer-Rao
Cramer−Rao界限指出,泛化误差的下降速度不会快于
O
(
1
k
)
O(\frac{1}{k})
O(k1)。
B
o
t
t
o
u
a
n
d
B
o
u
s
q
u
e
t
Bottou\quad and\quad Bousquet
BottouandBousquet因此认为对于机器学习任务,不值得探寻快于
O
(
1
k
)
O(\frac{1}{k})
O(k1)的优化算法——更快的收敛可能对应着过拟合。因此,渐进分析掩盖了随机梯度下降在少量更新步之后的很多优点。对于大数据集,
S
G
D
SGD
SGD只需非常少量样本计算梯度从而实现初始快速更新,远远超过了其缓慢的渐进收敛。后文讨论的大多数算法在实践中都受益于这种性质,但是损失了常数倍
O
(
1
k
)
O(\frac{1}{k})
O(k1)的渐进分析。我们也可以在学习过程中逐渐增大小批量的大小,以此权衡批量梯度下降和随机梯度下降两者的优点。
4 使用动量的随机梯度下降
4.1 背景
虽然随机梯度下降仍然是非常受欢迎的优化方法,但其学习过程有时会很慢。动量方法旨在加速学习,特别是处理高曲率、小但一致的梯度,或是带噪声的梯度。动量算法累积了之前梯度指数级衰减的移动平均,并且继续沿该方向移动。
4.2 物理学基础
从形式上看,动量算法引入了变量 v v v充当速度角色——它代表参数在参数空间移动的方向和速率。速度被设为负梯度的指数衰减平均。名称动量(momentum)来自物理类比,根据牛顿运动定律,负梯度是移动空间中粒子的力。动量在物理学上定义为质量乘以速度。
4.3 思想
在动量学习算法中,我们假设是单位质量,因此速度向量
v
v
v也可以看作粒子的动量。超参数
α
∈
[
0
,
1
)
\alpha \in [0,1)
α∈[0,1)决定了之前梯度的贡献衰减得有多慢。更新规则如下:
(3.4)
v
←
α
v
−
ϵ
∇
θ
(
1
m
Σ
i
=
1
m
L
(
f
(
x
(
i
)
;
θ
)
,
y
(
i
)
)
)
v \leftarrow \alpha v - \epsilon \nabla_\theta \left(\frac{1}{m} \Sigma_{i=1}^{m}L(f(x^{(i)};\theta),y^{(i)}) \right) \tag{3.4}
v←αv−ϵ∇θ(m1Σi=1mL(f(x(i);θ),y(i)))(3.4)
(3.5)
θ
←
θ
+
v
\theta \leftarrow \theta + v \tag{3.5}
θ←θ+v(3.5)
速度
v
v
v累积了梯度元素
∇
θ
(
1
m
Σ
i
=
1
m
L
(
f
(
x
(
i
)
;
θ
)
,
y
(
i
)
)
)
\nabla_\theta \left(\frac{1}{m} \Sigma_{i=1}^{m}L(f(x^{(i)};\theta),y^{(i)}) \right)
∇θ(m1Σi=1mL(f(x(i);θ),y(i)))。相对于
ϵ
\epsilon
ϵ,
α
\alpha
α越大,之前梯度对现在方向的影响也越大。带动量的
S
G
D
SGD
SGD算法如下所示。
4.4 具体算法
使用动量的随机梯度下降(SGD)
输入:学习率 ϵ \epsilon ϵ,动量参数 α \alpha α
输入:初始参数 θ \theta θ,初始速度 v v v
\quad while 停止准则未满足 do
\quad \quad 从训练集中采包含 m m m个样本 x ( 1 ) , . . . , x ( m ) {x^{(1)},...,x^{(m)}} x(1),...,x(m)的小批量,对应目标为 y ( i ) y^{(i)} y(i)
\quad \quad 计算梯度估计: g ← 1 m Σ i = 1 m L ( f ( x ( i ) ; θ ) , y ( i ) ) g \leftarrow \frac{1}{m} \Sigma_{i=1}^{m}L(f(x^{(i)};\theta),y^{(i)}) g←m1Σi=1mL(f(x(i);θ),y(i))
\quad \quad 计算速度更新: v ← α v − ϵ g v \leftarrow \alpha v - \epsilon g v←αv−ϵg
\quad \quad 应用更新: θ ← + v \theta \leftarrow + v θ←+v
\quad end while
4.5 讲解
之前,步长只是梯度范数乘以学习率。现在,步长取决于梯度序列的大小和排列。当许多连续的梯度指向相同的方向时,步长最大。如果动量算法总是观测到梯度
g
g
g,那么它会在方向
−
g
-g
−g上不停加速,直到达到最终速度,其中步长大小为
(3.6)
ϵ
∣
∣
g
∣
∣
1
−
α
\frac{\epsilon ||g||}{1-\alpha} \tag{3.6}
1−αϵ∣∣g∣∣(3.6)
因此将动量的超参数视为
1
1
−
α
\frac{1}{1-\alpha}
1−α1有助于理解。例如,
α
=
0.9
\alpha=0.9
α=0.9对应着最大速度10倍于梯度下降算法。
在实践中,
α
\alpha
α的一般取值为
0.5
,
0.9
和
0.99
0.5,0.9和0.99
0.5,0.9和0.99。和学习率一样,
α
\alpha
α也会随着时间不断调整。一般初始值是一个较小的值,随后会慢慢变大。随着时间推移调整
α
\alpha
α没有收缩
ϵ
\epsilon
ϵ重要。
5 自适应学习率算法
5.1 背景
神经网络研究员早就意识到学习率肯定是难以设置的超参数之一,因为它对模型的性能有显著的影响。损失通常高度敏感与参数空间中的某些方向,而不敏感与其他。动量算法可以在一定程度缓解这些问题。但这样做的代价是引入了另一个超参数。在这种情况下,自然会问有没有其他方法。如果我们相信方向敏感度是轴对齐的,那么每个参数设置不同的学习率,在整个学习过程中自动适应这些学习率是有道理的。
5.2 AdaGrad
5.2.1 特点
独立地适应所有模型参数的学习率,缩放每个参数反比于其所有梯度历史平方值总和的平方根。具有损失最大偏导的参数相应地有一个快速下降的学习率,而具有小偏导的参数在学习率上有相对较小的下降。净效果是在参数空间中更为平缓的倾斜方向会取得更大的进步。
在凸优化背景中,
A
d
a
G
r
a
d
AdaGrad
AdaGrad具有一些令人满意的理论性质。然而,经验已经发现,对于训练深度神经网络模型而言,从训练开始时累计梯度平方会导致有效学习率过早或过量的减少。
A
d
a
G
r
a
d
AdaGrad
AdaGrad在某些深度学习模型上效果不错,但不是全部。
5.2.2 具体算法
A d a G r a d AdaGrad AdaGrad 算法
输入:全局学习率 ϵ \epsilon ϵ
输入:初始参数 θ \theta θ
输入:小常数 δ \delta δ,为了数值稳定大约设置为 1 0 − 7 10^{-7} 10−7
\quad 初始化梯度累计变量 r = 0 r=0 r=0
\quad while 没有达到停止标准
\quad \quad 从训练集中采集包含 m m m个样本{ x ( 1 ) , . . . , x ( m ) x^{(1)},...,x^{(m)} x(1),...,x(m)}的小批量,其中 x ( i ) x^{(i)} x(i)对应目标为 y ( i ) y^{(i)} y(i)
\quad \quad 计算梯度: g ← 1 m Σ i L ( f ( x ( i ) ; θ ) , y ( i ) ) g \leftarrow \frac{1}{m} \Sigma_iL(f(x^{(i)};\theta),y^{(i)}) g←m1ΣiL(f(x(i);θ),y(i))
\quad \quad 累计平方梯度: r ← r + g ⨀ g r \leftarrow r + g \bigodot g r←r+g⨀g
\quad \quad 计算参数更新: ∇ θ = − ϵ δ + r ⨀ g \nabla\theta=-\frac{\epsilon}{\sqrt{\delta+r}}\bigodot g ∇θ=−δ+rϵ⨀g(逐元素地应用除和求平方根)
\quad \quad 应用更新: θ ← θ + Δ θ \theta \leftarrow \theta + \Delta\theta θ←θ+Δθ
\quad end while
5.3 RMSProp
5.3.1 特点
R M S P r o p RMSProp RMSProp算法修改 A d a G r a d AdaGrad AdaGrad以在非凸设定下效果更好,改变梯度累计为指数加权的移动平均。 A d a G r a d AdaGrad AdaGrad旨在应用于凸问题时快速收敛。当应用于非凸函数训练神经网络时,学习轨迹可能穿过了很多不同的结构,最终达到一个局部是凸碗的区域。 A d a G r a d AdaGrad AdaGrad根据平方梯度的整个历史收缩学习率,可能使得学习率在达到这样的凸结构前就变得发太小了。 R M S P r o p RMSProp RMSProp使用指数衰减平均以丢弃遥远过去的历史,使其能够在找到凸碗状结构后迅速收敛,它就像一个初始化于该碗状结构的 A d a G r a d AdaGrad AdaGrad算法实例。
5.3.2 具体算法
R M S P r o p RMSProp RMSProp的标准形式如下所示。相比于 A d a G r a d AdaGrad AdaGrad,使用移动平均引入一个新的超参数 ρ \rho ρ,用来控制移动平均的长度范围。
输入:全局学习率 ϵ \epsilon ϵ
输入:初始参数 θ \theta θ
输入:小常数 δ \delta δ,为了数值稳定大约设置为 1 0 − 6 10^{-6} 10−6(用于被小数除时的数值稳定)
\quad 初始化梯度累计变量 r = 0 r=0 r=0
\quad while 没有达到停止标准
\quad \quad 从训练集中采集包含 m m m个样本{ x ( 1 ) , . . . , x ( m ) x^{(1)},...,x^{(m)} x(1),...,x(m)}的小批量,其中 x ( i ) x^{(i)} x(i)对应目标为 y ( i ) y^{(i)} y(i)
\quad \quad 计算梯度: g ← 1 m Σ i L ( f ( x ( i ) ; θ ) , y ( i ) ) g \leftarrow \frac{1}{m} \Sigma_iL(f(x^{(i)};\theta),y^{(i)}) g←m1ΣiL(f(x(i);θ),y(i))
\quad \quad 累计平方梯度: r ← ρ r + ( 1 − ρ ) g ⨀ g r \leftarrow \rho r + (1-\rho)g \bigodot g r←ρr+(1−ρ)g⨀g
\quad \quad 计算参数更新: ∇ θ = − ϵ δ + r ⨀ g \nabla\theta=-\frac{\epsilon}{\sqrt{\delta+r}}\bigodot g ∇θ=−δ+rϵ⨀g( ϵ δ + r \frac{\epsilon}{\sqrt{\delta+r}} δ+rϵ逐元素应用)
\quad \quad 应用更新: θ ← θ + Δ θ \theta \leftarrow \theta + \Delta\theta θ←θ+Δθ
\quad end while
5.3.3 讲解
经验上, R M S P r o p RMSProp RMSProp已被证明是一种有效且使用的深度神经网络优化算法。目前它是深度学习从业者经常采用的优化方法之一。
5.4 Adam
5.4.1 特点
A d a m Adam Adam是另一种学习率自适应的优化算法。如下所示:
5.4.2 具体算法
输入:步长 ϵ \epsilon ϵ(建议默认为:0.001)
输入:矩估计的指数衰减速率, ρ 1 和 ρ 2 \rho_1和\rho_2 ρ1和ρ2在区间 [ 0 , 1 ) [0,1) [0,1)内。(建议默认为:分别为 0.9 0.9 0.9和 0.999 0.999 0.999)
输入:用于数值稳定的小常数 δ \delta δ(建议默认为: 1 0 − 8 10^{-8} 10−8)
输入:初始参数 θ \theta θ
\quad 初始化一阶和二阶矩变量 s = 0 , r = 1 s=0,r=1 s=0,r=1
\quad 初始化时间步 t = 0 t=0 t=0
\quad while 没有达到停止准则 do
\quad \quad 从训练集中采集包含 m m m个样本{ x ( 1 ) , . . . , x ( m ) x^{(1)},...,x^{(m)} x(1),...,x(m)}的小批量,其中 x ( i ) x^{(i)} x(i)对应目标为 y ( i ) y^{(i)} y(i)
\quad \quad 计算梯度: g ← 1 m Σ i L ( f ( x ( i ) ; θ ) , y ( i ) ) g \leftarrow \frac{1}{m} \Sigma_iL(f(x^{(i)};\theta),y^{(i)}) g←m1ΣiL(f(x(i);θ),y(i))
\quad \quad t ← t + 1 t \leftarrow t + 1 t←t+1
\quad \quad 更新有偏一阶矩估计: s ← ρ 1 s + ( 1 − ρ 1 ) g s \leftarrow \rho_1 s + (1-\rho_1) g s←ρ1s+(1−ρ1)g
\quad \quad 更新有偏二阶矩估计: r ← ρ 2 r + ( 1 − ρ 2 ) g ⨀ g r \leftarrow \rho_2 r + (1-\rho_2) g \bigodot g r←ρ2r+(1−ρ2)g⨀g
\quad \quad 更新修正一阶矩的偏差: s 2 ← s 1 − p 1 t s_2 \leftarrow \frac{s}{1-p_1^t} s2←1−p1ts
\quad \quad 更新修正二阶矩的偏差: r 2 ← r 1 − p 2 t r_2 \leftarrow \frac{r}{1-p_2^t} r2←1−p2tr
\quad \quad 计算更新: Δ θ = − ϵ s 2 r 2 + δ ⨀ g \Delta_\theta=-\epsilon \frac{s_2}{\sqrt{r_2+\delta}}\bigodot g Δθ=−ϵr2+δs2⨀g(逐元素应用操作)
\quad \quad 应用更新: θ ← θ + Δ θ \theta \leftarrow \theta + \Delta\theta θ←θ+Δθ
\quad end while
5.4.3 讲解
“Adam"这个名字派生自短语"adaptive moments”。早期算法背景下,它也许是最好被看作结合 R M S P r o p RMSProp RMSProp和具有一些重要区别的动量的变种。首先,在Adam中,动量直接并入了梯度一阶矩(指数加权)的估计。将动量加入RMSProp最直观的方法是将动量应用于缩放后的梯度。结合缩放的动量使用没有明确的理论动机。其次,Adam包括偏置修正,修正从原点初始化的一阶矩(动量项)和(非中心的)二阶矩的估计。RMSProp也采用了(非中心的)二阶矩估计,然而缺失了修正因子。因此,不像Adam,RMSProp二阶矩估计可能在训练初期有很高的偏置。Adam通常被认为对超参数的选择相当鲁棒,尽管学习率有时需要从建议的默认修改。
6 总结——选择正确的优化算法
在本文中,我们讨论了一系列算法,通过自适应每个模型参数的学习率以解决优化深度模型中的难题。此时,一个自然的问题是:该选择哪种算法呢?
遗憾的是,目前在这一点上没有达成共识。Schaul et al展示了许多优化算法在大量学习任务上极具价值的比较。虽然结果表明,据欧自适应学习率(以RMSProp和AdaDelta为代表)得算法族表现得相当鲁棒,不分伯仲,但没有那个算法能够脱颖而出。
目前,最流行并且使用很高的优化算法包括SGD、具动量的SGD、RMSProp、具动量的RMSProp、AdaDelta和Adam。此时,选择那一个算法似乎主要取决于使用者对算法的熟悉程度(以便调节超参数)。
参考文献:
1 统计学习方法 清华大学出版社
2 深度学习 人民邮电出版社