欢迎查看我的博客文章合集:我的Blog文章索引::机器学习方法系列,深度学习方法系列,三十分钟理解系列等
拟牛顿法
拟牛顿法可以克服牛顿法计算量大的缺点,不在计算目标函数的 Hesse 矩阵,而是构造一个近似 Hesse 矩阵的对称正定矩阵,根据近似矩阵来优化目标函数,不同的近似构造 Hesse 的方法决定了不同的拟牛顿法,构造 Hesse 矩阵是需要满足拟牛顿条件的,拟牛顿条件是这样求得的,首先将 f(x) 在 x k + 1 x_{k+1} xk+1处做二阶泰勒展开(忽略高阶项):
f ( x ) = f ( x k + 1 ) + ∇ f ( x k + 1 ) ( x − x k + 1 ) + 1 2 ( x – x k + 1 ) T ∇ 2 f ( x k + 1 ) ( x − x k + 1 ) f(x) = f(x_{k+1}) + \nabla f(x_{k+1}) (x-x_{k+1})+ \frac{1}{2} (x –x_{k+1})^T\nabla^2 f(x_{k+1} )(x-x_{k+1}) f(x)=f(xk+1)+∇f(xk+1)(x−xk+1)+21(x–xk+1)T∇2f(xk+1)(x−xk+1)
注意在这个式子中,
x
x
x是变量,而
x
k
+
1
x_{k+1}
xk+1是一个值,对
x
x
x求导得到:
∇
f
(
x
)
=
0
+
∇
f
(
x
k
+
1
)
+
H
k
+
1
(
x
−
x
k
+
1
)
\nabla f(x) = 0 + \nabla f(x_{k+1}) + H_{k+1}(x -x_{k+1})
∇f(x)=0+∇f(xk+1)+Hk+1(x−xk+1)整理得到:
g
=
g
k
+
1
+
H
k
+
1
(
x
−
x
k
+
1
)
g = g_{k+1}+ H_{k+1}(x -x_{k+1})
g=gk+1+Hk+1(x−xk+1)
令
x
=
x
k
x=x_k
x=xk ,整理可得
g
k
+
1
–
g
k
=
H
k
+
1
(
x
k
+
1
–
x
k
)
g_{k+1} – g_k = H_{k+1} (x_{k+1} – x_k)
gk+1–gk=Hk+1(xk+1–xk)
这个便是拟牛顿条件了,迭代过程中对
H
k
+
1
H_{k+1}
Hk+1 做出约束,根据约束构造一个近似矩阵
B
k
+
1
B_{k+1}
Bk+1 ,来模拟 Hesse 矩阵就可以了,为了简便起见,引入记号
s
k
s_k
sk 与
y
k
y_k
yk ,令
s
k
=
x
k
+
1
–
x
k
,
y
k
=
g
k
+
1
–
g
k
s_k = x_{k+1} –x_k , y_k = g_{k+1} –g _k
sk=xk+1–xk,yk=gk+1–gk
y
k
=
B
k
+
1
⋅
s
k
y_k = B_{k+1} \cdot s_k
yk=Bk+1⋅sk
因为牛顿法中的迭代方向为
−
H
−
1
⋅
g
-H^{-1} \cdot g
−H−1⋅g,所以令
D
k
+
1
=
H
k
+
1
−
1
D_{k+1} = H_{k+1}^{-1}
Dk+1=Hk+1−1,拟牛顿条件还可以写作:
s
k
=
D
k
+
1
⋅
y
k
s_k = D_{k+1} \cdot y_k
sk=Dk+1⋅yk
拟牛顿法本身是一类算法,下面介绍一下BFGS,算是比较著名的方法了:
BFGS算法(Broyden–Fletcher–Goldfarb–Shanno)[3]
BFGS 是一种拟牛顿方法,通过迭代构建近似 Hesse 矩阵,省去了求解 Hesse 的复杂的步骤,而且 BFGS 构造出来的近似 Hesse 矩阵一定是正定的,这完全克服了牛顿法的缺陷,虽然搜索方向不一定最优,但始终朝着最优的方向前进的。首先初始化 Hesse 矩阵
B
0
=
I
B_0=I
B0=I ,接下来每次迭代对矩阵
B
k
B_k
Bk 进行更新即可:
B
k
+
1
=
B
k
+
Δ
B
k
,
k
=
1
,
2
,
…
B_{k+1} = B_k+ \Delta B_k , \ k = 1,2,…
Bk+1=Bk+ΔBk, k=1,2,…
迭代构建近似矩阵的关键是矩阵 ΔBk 的构造了,将其写作:
Δ
B
k
=
α
u
u
T
+
β
v
v
T
\Delta B_k = \alpha uu^T + \beta vv^T
ΔBk=αuuT+βvvT
这里的向量 u 和 v 是待定的,知道了这两个向量,就可以构造构造 ΔBk 了,且这样构造出的矩阵是对称的,根据拟牛顿条件:
y
k
=
B
k
+
1
s
k
=
(
B
k
+
Δ
B
k
)
s
k
=
(
B
k
+
α
u
u
T
+
β
v
v
T
)
s
k
=
B
k
s
k
+
(
α
u
T
s
k
)
⋅
u
+
(
β
v
T
s
k
)
⋅
v
\begin{aligned} y_k &= B_{k+1} s_k \\ &= (B_k + \Delta B_k)s_k \\ &= (B_k + \alpha uu^T + \beta vv^T)s_k \\ &= B_k s_k + (\alpha u^Ts_k) \cdot u+ (\beta v^Ts_k) \cdot v \end{aligned}
yk=Bk+1sk=(Bk+ΔBk)sk=(Bk+αuuT+βvvT)sk=Bksk+(αuTsk)⋅u+(βvTsk)⋅v
这里
α
u
T
s
k
αu^Ts_k
αuTsk 与
β
v
T
s
k
βv^Tsk
βvTsk 均为实数,代表了 在 u 与 v 方向的拉伸程度,为了计算简单,做如下赋值运算:
α
u
T
s
k
=
1
,
β
v
T
s
k
=
–
1
\alpha u^Ts_k = 1 , \ \beta v^Ts_k = –1
αuTsk=1, βvTsk=–1
代入上式便可得:
u
−
v
=
y
k
–
B
k
s
k
u - v = y_k – B_ks_k
u−v=yk–Bksk
这就得到得到了 u 与 v 的一个近似:
u
=
y
k
,
v
=
B
k
s
k
u = y_k , \ v = B_k s_k
u=yk, v=Bksk
继而求 α 与 β 的值
α
=
1
y
T
s
k
,
β
=
−
1
(
B
k
s
k
)
T
s
k
=
−
1
s
k
T
B
k
s
k
\alpha = \frac{1}{y^Ts_k}, \beta= -\frac{1}{(B_ks_k)^Ts_k} = -\frac{1}{s_k^TB_ks_k}
α=yTsk1,β=−(Bksk)Tsk1=−skTBksk1
α 、 β 、 u 与 v都求得后,便得到了 ΔBk 的更新公式:
Δ
B
k
=
y
k
y
k
T
y
k
T
s
k
–
B
k
s
k
s
k
T
B
k
s
k
T
B
k
s
k
\Delta B_k = \frac{y_ky_k^T}{y_k^Ts_k} – \frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}
ΔBk=ykTskykykT–skTBkskBkskskTBk
因此
B
k
B_k
Bk的迭代公式是:
B
k
+
1
=
B
k
+
y
k
y
k
T
y
k
T
s
k
–
B
k
s
k
s
k
T
B
k
s
k
T
B
k
s
k
B_{k+1} = B_k +\frac{y_ky_k^T}{y_k^Ts_k} – \frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}
Bk+1=Bk+ykTskykykT–skTBkskBkskskTBk
由与牛顿法的方向是 – H k − 1 g k –H^{−1}_kg_k –Hk−1gk 的,所以最好可以直接计算出 B k − 1 B^{−1}_k Bk−1 ,这样就不用再进行求逆运算了,直接根据Sherman-Morrison 公式:可得关于矩阵B 的逆的更新方式:
B k + 1 − 1 = B k − 1 + ( 1 s k T y k + y k T B k − 1 y k ( s k T y k ) 2 ) s k s k T − 1 s k T y k ( B k − 1 y k s k T + s k y k T B k − 1 ) B^{-1}_{k+1} = B^{-1}_k + \left (\frac{1}{s_k^Ty_k}+\frac{y_k^TB_k^{-1}y_k}{(s_k^Ty_k)^2} \right )s_ks_k^T - \frac{1}{s_k^Ty_k} \left (B_k^{-1}y_ks_k^T + s_ky_k^TB^{-1}_k \right) Bk+1−1=Bk−1+(skTyk1+(skTyk)2ykTBk−1yk)skskT−skTyk1(Bk−1ykskT+skykTBk−1)
B
k
−
1
B^{−1}_k
Bk−1 这里用
D
k
D_k
Dk 来表示,给出最终的 BFGS 算法[3]:
停止条件为人为设定,可设定为两次迭代目标函数差的阈值或者梯度差的阈值,或者梯度本身(的模)小于阈值。
其中,步骤2.2搜索步长的方法采用[7]:
比较好理解,就是在搜索方向p上,找到步长 α \alpha α满足Armijo条件,初始步长 α 0 = 1 \alpha_0=1 α0=1是一种常用的设定。
DFP算法
DFP算法也是类似的思想,可以参考[4],写的很详细,我这里简单贴一个图以备查阅:
稍微看下步3,选用的方法就是上面介绍过的Backtracking line search算法,只是选用的符号不一样而已,内容是一样的。非负整数m就是代表了迭代次数。
L-BFGS [3]
工业中实用的拟牛顿法的便是 L-BFGS (Limited-memory BFGS)了,对于近似 Hesse 矩阵
D
k
D_k
Dk :
D
k
+
1
=
D
k
+
(
1
s
k
T
y
k
+
y
k
T
D
k
y
k
(
s
k
T
y
k
)
2
)
s
k
s
k
T
−
1
s
k
T
y
k
(
D
k
y
k
s
k
T
+
s
k
y
k
T
D
k
)
D_{k+1} = D_k + \left (\frac{1}{s_k^Ty_k}+\frac{y_k^T D_ky_k}{(s_k^Ty_k)^2} \right )s_ks_k^T - \frac{1}{s_k^Ty_k}(D_ky_ks_k^T + s_ky_k^T D_k )
Dk+1=Dk+(skTyk1+(skTyk)2ykTDkyk)skskT−skTyk1(DkykskT+skykTDk)
而是存储向量序
s
k
s_k
sk,
y
k
y_k
yk,而且向量序列也不是都存,而是存最近的 m 次的, m 为人工指定,计算
D
k
D_k
Dk 时,只用最新的 m 个向量模拟计算即可。在第 k 次迭代,算法求得了
x
k
x_k
xk ,并且保存的曲率信息为
(
s
i
,
y
i
)
k
−
1
k
−
m
(si,yi)_{k−1}^{k−m}
(si,yi)k−1k−m。为了得到
H
k
H_k
Hk,算法每次迭代均需选择一个初始的矩阵
H
0
K
H_0^K
H0K,这是不同于 BFGS 算法的一个地方,接下来只用最近的 m 个向量对该初始矩阵进行修正,实践中
H
0
K
H_0^K
H0K 的设定通常如下:
H
k
0
=
r
k
I
r
k
=
s
k
−
1
T
y
k
−
1
y
k
−
1
T
y
k
−
1
\begin{aligned} H_k^0 &=r_kI \\ r_k &=\frac{s{k-1}^Ty_{k-1}}{y_{k-1}^Ty_{k-1}} \end{aligned}
Hk0rk=rkI=yk−1Tyk−1sk−1Tyk−1
其中 r k r_k rk 表示比例系数,它利用最近一次的曲率信息来估计真实海森矩阵的大小,这就使得当前步的搜索方向较为理想,而不至于跑得“太偏”,这样就省去了步长搜索的步骤,节省了时间。在L-BFGS算法中,通过保存最近 m 次的曲率信息来更新近似矩阵的这种方法在实践中是很有效的,虽然 L-BFGS 算法是线性收敛,但是每次迭代的开销非常小,因此 L-BFGS 算法执行速度还是很快的,而且由于每一步迭代都能保证近似矩阵的正定,因此算法的鲁棒性还是很强的。
总结下 BFGS 与 L-BFGS 的: BFGS算法在运行的时候,每一步迭代都需要保存一个 n×n 的矩阵,现在很多机器学习问题都是高维的,当 n 很大的时候,这个矩阵占用的内存是非常惊人的,并且所需的计算量也是很大的,这使得传统的 BFGS 算法变得非常不适用。而 L-BFGS 则是很对这个问题的改进版,从上面所说可知,BFGS 算法是通过曲率信息$ (s_k,y_k)$ 来修正 H k H_k Hk 从而得到 H k + 1 H_{k+1} Hk+1 ,L-BFGS 算法的主要思路是:算法仅仅保存最近 m 次迭代的曲率信息来计算 H k + 1 H_{k+1} Hk+1 。这样,我们所需的存储空间就从 n×n 变成了 2m×n 而通常情况下 m << n。
其他拟牛顿算法[6]
共轭梯度法
共轭梯度法是介于梯度下降法和牛顿法,拟牛顿法之间的算法[6]。
待补充…
参考资料
[1]https://blog.csdn.net/batuwuhanpei/article/details/51979831
[2]https://blog.csdn.net/u011722133/article/details/53518134
[3]无约束优化方法(梯度法-牛顿法-BFGS- L-BFGS)
[4]优化算法——拟牛顿法之DFP算法
[5]牛顿法与拟牛顿法
[6]牛顿法,拟牛顿法, 共轭梯度法
[7]【原创】回溯线搜索 Backtracking line search
[8]【原创】牛顿法和拟牛顿法 – BFGS, L-BFGS, OWL-QN
[9]无约束最优化方法——牛顿法、拟牛顿法、BFGS、LBFGS
[10] 无约束最优化的常用方法
[11]机器学习与运筹优化(三)从牛顿法到L-BFGS
[12]ECE236C - Optimization Methods for Large-Scale Systems (Spring 2019),UCLA的优化课程,课件很好。