1.牛顿法
统计学习方法有了具体形式后就转换为最优化问题。有时最优化问题存在解析解,可以由公式计算,多数情况下没有解析解,需要用数值计算的方法求解,牛顿法和拟牛顿法是求解无约束最优化问题的常用方法,收敛速度快。
牛顿法是迭代算法,每一步需要求解目标函数的海塞矩阵的逆矩阵。
1.1 算法推导
无约束最优化问题如下:
min
x
∈
R
n
f
(
x
)
.
.
.
.
.
.
.
.
.
.
(
1
)
\min_{x \in R^n} f(x)..........(1)
x∈Rnminf(x)..........(1)
一元函数在
x
k
x_k
xk处的泰勒展开式为
f
(
x
)
=
f
(
x
k
)
+
(
x
−
x
k
)
f
′
(
x
k
)
+
1
2
!
(
x
−
x
k
)
2
f
′
′
(
x
k
)
+
o
(
n
)
f(x)=f(x_k)+(x-x_k)f'(x_k)+\frac{1}{2!}(x-x_k)^2f''(x_k)+o(n)
f(x)=f(xk)+(x−xk)f′(xk)+2!1(x−xk)2f′′(xk)+o(n)
二元函数在
(
x
k
,
y
k
)
(x_k,y_k)
(xk,yk)处的泰勒展开式为
f
(
x
,
y
)
=
f
(
x
k
,
y
k
)
+
(
x
−
x
k
)
f
x
′
(
x
k
,
y
k
)
+
(
y
−
y
k
)
f
y
′
(
x
k
,
y
k
)
+
1
2
!
(
x
−
x
k
)
2
f
x
x
′
′
(
x
k
,
y
k
)
+
1
2
!
(
x
−
x
k
)
(
y
−
y
k
)
f
x
y
′
′
(
x
k
,
y
k
)
+
1
2
!
(
x
−
x
k
)
(
y
−
y
k
)
f
y
x
′
′
(
x
k
,
y
k
)
+
1
2
!
(
y
−
y
k
)
2
f
y
y
′
′
(
x
k
,
y
k
)
+
o
(
n
)
f(x,y)=f(x_k,y_k)+(x-x_k)f'_x(x_k,y_k)+(y-y_k)f'_y(x_k,y_k)\\+\frac{1}{2!}(x-x_k)^2f''_{xx}(x_k,y_k)+\frac{1}{2!}(x-x_k)(y-y_k)f''_{xy}(x_k,y_k)\\ +\frac{1}{2!}(x-x_k)(y-y_k)f''_{yx}(x_k,y_k)+\frac{1}{2!}(y-y_k)^2f''_{yy}(x_k,y_k)+o(n)
f(x,y)=f(xk,yk)+(x−xk)fx′(xk,yk)+(y−yk)fy′(xk,yk)+2!1(x−xk)2fxx′′(xk,yk)+2!1(x−xk)(y−yk)fxy′′(xk,yk)+2!1(x−xk)(y−yk)fyx′′(xk,yk)+2!1(y−yk)2fyy′′(xk,yk)+o(n)
多元函数在
x
⃗
k
\vec x_k
xk处的泰勒展开式为
f
(
x
(
1
)
,
x
(
2
)
,
.
.
.
,
x
(
n
)
)
=
f
(
x
k
(
1
)
,
x
k
(
2
)
,
.
.
.
x
k
(
n
)
)
+
∑
i
=
1
n
(
x
(
i
)
−
x
k
(
i
)
)
f
x
(
i
)
′
(
x
k
(
1
)
,
x
k
(
2
)
,
.
.
.
x
k
(
n
)
)
+
1
2
!
∑
i
,
j
=
1
n
(
x
i
−
x
k
(
i
i
)
)
(
x
(
j
)
−
x
k
(
j
)
)
f
i
j
′
′
(
x
k
(
1
)
,
x
k
(
2
)
,
.
.
.
x
k
(
n
)
)
+
o
(
n
)
f(x^{(1)},x^{(2)},...,x^{(n)})=f(x_k^{(1)},x_k^{(2)},...x_k^{(n)})\\ +\sum_{i=1}^n(x^{(i)}-x_k^{(i)})f'_{x^{(i)}}(x_k^{(1)},x_k^{(2)},...x_k^{(n)})+ \\ \frac{1}{2!}\sum_{i,j=1}^n(x^{i}-x_k^{(ii)})(x^{(j)}-x_k^{(j)})f_{ij}''(x_k^{(1)},x_k^{(2)},...x_k^{(n)})+o(n)
f(x(1),x(2),...,x(n))=f(xk(1),xk(2),...xk(n))+i=1∑n(x(i)−xk(i))fx(i)′(xk(1),xk(2),...xk(n))+2!1i,j=1∑n(xi−xk(ii))(x(j)−xk(j))fij′′(xk(1),xk(2),...xk(n))+o(n)
假设f(x)有二阶连续偏导数,若第k次迭代值为
x
k
x_k
xk,则可以将f(x)在
x
k
x_k
xk附近二阶泰勒展开,并将泰勒展开写成矩阵的形式,有
f
(
x
)
=
f
(
x
k
)
+
[
∇
f
(
x
k
)
]
T
(
x
−
x
k
)
+
1
2
!
(
x
−
x
k
)
T
H
(
x
k
)
(
x
−
x
k
)
+
o
(
n
)
.
.
.
.
.
.
.
.
.
.
(
2
)
f(x)=f(x_k)+[\nabla f(x_k)]^T(x-x_k)+\frac{1}{2!}(x-x_k)^TH(x_k)(x-x_k)+o(n)..........(2)
f(x)=f(xk)+[∇f(xk)]T(x−xk)+2!1(x−xk)TH(xk)(x−xk)+o(n)..........(2)
其中x为列向量,
[
∇
f
(
x
k
)
]
[\nabla f(x_k)]
[∇f(xk)]为f(x)的梯度向量在
x
k
x_k
xk的值,
H
(
x
k
)
H(x_k)
H(xk)是f(x)的海塞矩阵在点
x
k
x_k
xk的值,其形式如下
H
(
x
)
=
[
∂
2
f
∂
x
i
∂
x
j
]
n
×
n
.
.
.
.
.
.
.
.
.
.
(
3
)
H(x)=\left[\frac{\partial^2f}{\partial x_i \partial x_j}\right]_{n\times n}..........(3)
H(x)=[∂xi∂xj∂2f]n×n..........(3)
我们知道一元函数取得极值点的条件是导数等于0,二元函数取得极值点的条件是一阶导数等于0,并且二阶导数要大于0,对应到上面的泰勒展开,只要
(
x
−
x
k
)
T
H
(
x
k
)
(
x
−
x
k
)
>
0
(x-x_k)^TH(x_k)(x-x_k)>0
(x−xk)TH(xk)(x−xk)>0即可,而这个满足线性代数中二次型的特征,当H矩阵为正定矩阵时,该不等式成立。即
- 当海塞矩阵H为正定矩阵时,临界点 x k x_k xk为局部极小值。
- 当海塞矩阵H为负定矩阵时,临界点 x k x_k xk为局部极大值。
- 当海塞矩阵H为不定矩阵时,临界点 x k x_k xk不是极值。
牛顿法利用这一点,每次迭代从点
x
k
x_k
xk开始,求目标函数的极小值点,作为第k+1次迭代值
x
k
+
1
x_{k+1}
xk+1,则
x
k
+
1
x_{k+1}
xk+1满足
∇
f
(
x
k
+
1
)
=
0..........
(
4
)
\nabla f(x_{k+1})=0..........(4)
∇f(xk+1)=0..........(4)
式子2求导得,
∇
f
(
x
)
=
∇
f
(
x
k
)
+
H
(
x
k
)
(
x
−
x
k
)
.
.
.
.
.
.
.
.
.
.
(
5
)
\nabla f(x)=\nabla f(x_k)+H(x_k)(x-x_k)..........(5)
∇f(x)=∇f(xk)+H(xk)(x−xk)..........(5)
求导公式可以参考【机器学习总结】向量、矩阵求导公式
则式子4等价于
∇
f
(
x
k
)
+
H
(
x
k
)
(
x
k
+
1
−
x
k
)
=
0.........
(
6
)
\nabla f(x_k)+H(x_k)(x_{k+1}-x_k)=0.........(6)
∇f(xk)+H(xk)(xk+1−xk)=0.........(6)
则
x
k
+
1
=
x
k
−
H
(
x
k
)
−
1
∇
f
(
x
k
)
.
.
.
.
.
.
.
.
.
(
7
)
x_{k+1}=x_k-H(x_k)^{-1}\nabla f(x_k).........(7)
xk+1=xk−H(xk)−1∇f(xk).........(7)
牛顿法以式子7进行迭代求得极小值,其流程总结如下:
- 取初始点 x 0 x_0 x0,令k=0
- 计算 ∇ f ( x k ) \nabla f(x_k) ∇f(xk),若 ∣ ∣ ∇ f ( x k ) ∣ ∣ < ϵ ||\nabla f(x_k)|| \lt \epsilon ∣∣∇f(xk)∣∣<ϵ,则停止计算,得到近似解 x ∗ = x k x^*=x_k x∗=xk
- 计算 H ( x k ) H(x_k) H(xk)并根据 x k + 1 = x k − H ( x k ) − 1 ∇ f ( x k ) x_{k+1}=x_k-H(x_k)^{-1}\nabla f(x_k) xk+1=xk−H(xk)−1∇f(xk)求得 x k + 1 x_{k+1} xk+1
- 令k=k+1,转2
2. 拟牛顿法
牛顿法中每轮迭代都需要计算海塞矩阵的逆矩阵,计算比较复杂,考虑用一个n阶矩阵
G
k
=
G
(
x
k
)
G_k=G(x_k)
Gk=G(xk)来近似代替
H
k
−
1
=
H
−
1
(
x
k
)
H_k^{-1}=H^{-1}(x_k)
Hk−1=H−1(xk)。
式子5中,将
x
=
x
k
+
1
x=x_{k+1}
x=xk+1代入,有
∇
f
(
x
k
+
1
)
−
∇
f
(
x
k
)
=
H
k
(
x
k
+
1
−
x
k
)
.
.
.
.
.
.
.
.
.
(
8
)
\nabla f(x_{k+1})-\nabla f(x_k)=H_k(x_{k+1}-x_k).........(8)
∇f(xk+1)−∇f(xk)=Hk(xk+1−xk).........(8)
令
y
k
=
∇
f
(
x
k
+
1
)
−
∇
f
(
x
k
)
y_k=\nabla f(x_{k+1})-\nabla f(x_k)
yk=∇f(xk+1)−∇f(xk),
δ
k
=
x
k
+
1
−
x
k
\delta_k=x_{k+1}-x_k
δk=xk+1−xk,则
y
k
=
H
k
δ
k
.
.
.
.
.
.
.
.
.
.
(
9
)
y_k=H_k\delta_k..........(9)
yk=Hkδk..........(9)
式子(9)称为拟牛顿条件。
如果
H
k
H_k
Hk是正定矩阵,则可以保证牛顿法搜索方向
p
k
=
−
H
k
−
1
∇
f
(
x
k
)
p_k=-H_k^{-1}\nabla f(x_k)
pk=−Hk−1∇f(xk)是下降方向。由式子7有
x
=
x
k
+
λ
p
k
=
x
k
−
λ
H
k
−
1
∇
f
(
x
k
)
.
.
.
.
.
.
.
.
.
.
(
10
)
x=x_k+\lambda p_k=x_k-\lambda H_k^{-1}\nabla f(x_k)..........(10)
x=xk+λpk=xk−λHk−1∇f(xk)..........(10)
则f(x)在
x
k
x_k
xk的泰勒展开式可以近似为
f
(
x
)
=
f
(
x
k
)
−
λ
[
∇
f
(
x
k
)
]
T
H
k
−
1
∇
f
(
x
k
)
.
.
.
.
.
.
.
.
.
.
(
11
)
f(x)=f(x_k)-\lambda [\nabla f(x_k)]^T H_k^{-1}\nabla f(x_k)..........(11)
f(x)=f(xk)−λ[∇f(xk)]THk−1∇f(xk)..........(11)
由于H的逆矩阵为正定矩阵,故有
[
∇
f
(
x
k
)
]
T
H
k
−
1
∇
f
(
x
k
)
>
0
[\nabla f(x_k)]^T H_k^{-1}\nabla f(x_k)\gt 0
[∇f(xk)]THk−1∇f(xk)>0,当
λ
\lambda
λ是充分小的正数时,总有
f
(
x
)
<
f
(
x
k
)
f(x) \lt f(x_k)
f(x)<f(xk) ,则
p
k
=
[
∇
f
(
x
k
)
]
T
H
k
−
1
∇
f
(
x
k
)
p_k= [\nabla f(x_k)]^T H_k^{-1}\nabla f(x_k)
pk=[∇f(xk)]THk−1∇f(xk)是下降方向。
拟牛顿法寻找
G
k
G_k
Gk来近似
H
−
1
H^{-1}
H−1,
G
k
G_k
Gk满足以下拟牛顿条件
G
k
+
1
y
k
=
δ
k
.
.
.
.
.
.
.
.
.
.
(
12
)
G_{k+1}y_k=\delta_k..........(12)
Gk+1yk=δk..........(12)
2.1 DFP算法
DFP算法假设每步迭代中,有
G
k
+
1
=
G
k
+
P
k
+
Q
k
.
.
.
.
.
.
.
.
.
.
(
13
)
G_{k+1}=G_k+P_k+Q_k..........(13)
Gk+1=Gk+Pk+Qk..........(13)
则
G
k
+
1
y
k
=
G
k
y
k
+
P
k
y
k
+
Q
k
y
k
.
.
.
.
.
.
.
.
.
.
(
14
)
G_{k+1}y_k=G_ky_k+P_ky_k+Q_ky_k..........(14)
Gk+1yk=Gkyk+Pkyk+Qkyk..........(14)
为了满足拟牛顿条件(12),可以让
P
k
y
k
=
δ
k
Q
k
y
k
=
−
G
k
y
k
P_ky_k=\delta_k\\ Q_ky_k=-G_ky_k
Pkyk=δkQkyk=−Gkyk
事实上可以取
P
k
=
δ
k
δ
k
T
δ
T
y
k
.
.
.
.
.
.
.
.
.
(
15
)
Q
k
=
−
G
k
y
k
y
k
T
G
k
y
k
T
G
k
y
k
.
.
.
.
.
.
.
.
.
.
(
16
)
P_k=\frac{\delta_k\delta_k^T}{\delta^Ty_k}.........(15)\\ Q_k=-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k}..........(16)
Pk=δTykδkδkT.........(15)Qk=−ykTGkykGkykykTGk..........(16)
得到
G
k
+
1
G_{k+1}
Gk+1的迭代公式
G
k
+
1
=
G
k
+
δ
k
δ
k
T
δ
T
y
k
−
−
G
k
y
k
y
k
T
G
k
y
k
T
G
k
y
k
.
.
.
.
.
.
.
.
.
.
(
17
)
G_{k+1}=G_k+\frac{\delta_k\delta_k^T}{\delta^Ty_k}--\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k}..........(17)
Gk+1=Gk+δTykδkδkT−−ykTGkykGkykykTGk..........(17)
可以证明如果
G
0
G_0
G0是正定的,则迭代过程中每个矩阵都是正定的。
DFP算法总结如下:
- 选定初始点 x 0 x_0 x0,取 G 0 G_0 G0为正定对称矩阵,置k=0
- 计算 ∇ f ( x k ) \nabla f(x_k) ∇f(xk),若 ∣ ∣ ∇ f ( x k ) ∣ ∣ < ϵ ||\nabla f(x_k)|| \lt \epsilon ∣∣∇f(xk)∣∣<ϵ,则停止计算,得到近似解 x ∗ = x k x^*=x_k x∗=xk
- 置 p k = − G k ∇ f ( x k ) p_k=-G_k\nabla f(x_k) pk=−Gk∇f(xk)
- 一维搜索,求 λ k \lambda _k λk使得 f ( x k + λ p k ) f(x_k+\lambda p_k) f(xk+λpk)最小
- 置 x k + 1 = x k + λ k p k x_{k+1}=x_k+\lambda_kp_k xk+1=xk+λkpk
- 计算 ∇ f ( x k + 1 ) \nabla f(x_{k+1}) ∇f(xk+1),若 ∣ ∣ ∇ f ( x k + 1 ) ∣ ∣ < ϵ ||\nabla f(x_{k+1})|| \lt \epsilon ∣∣∇f(xk+1)∣∣<ϵ,则停止计算,得到近似解 x ∗ = x k + 1 x^*=x_{k+1} x∗=xk+1,否则按式子 G k + 1 = G k + δ k δ k T δ 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^Ty_k}--\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} Gk+1=Gk+δTykδkδkT−−ykTGkykGkykykTGk计算得到 G k + 1 G_{k+1} Gk+1
- 置k=k+1,转3
2.2 BFGS(Broyden-Fletcher-Goldfarb-Shanno)算法
除了可以用
G
k
G_k
Gk逼近
H
−
1
H^{-1}
H−1之外,我们可以用
B
k
B_k
Bk来逼近
H
H
H,此时拟牛顿条件为
B
k
+
1
δ
k
=
y
k
.
.
.
.
.
.
.
.
.
.
(
18
)
B_{k+1}\delta_k=y_k..........(18)
Bk+1δk=yk..........(18)
跟DFP算法类似,
B
k
+
1
B_{k+1}
Bk+1的迭代公式为
B
k
+
1
=
B
k
+
y
k
y
k
T
y
T
δ
k
−
−
B
k
δ
k
δ
k
T
B
k
δ
k
T
B
k
δ
k
.
.
.
.
.
.
.
.
.
(
19
)
B_{k+1}=B_k+\frac{y_ky_k^T}{y^T\delta_k}--\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k}.........(19)
Bk+1=Bk+yTδkykykT−−δkTBkδkBkδkδkTBk.........(19)
BFGS算法总结如下:
- 选定初始点 x 0 x_0 x0,取 G 0 G_0 G0为正定对称矩阵,置k=0
- 计算 ∇ f ( x k ) \nabla f(x_k) ∇f(xk),若 ∣ ∣ ∇ f ( x k ) ∣ ∣ < ϵ ||\nabla f(x_k)|| \lt \epsilon ∣∣∇f(xk)∣∣<ϵ,则停止计算,得到近似解 x ∗ = x k x^*=x_k x∗=xk
- 由 B k p k = − ∇ f ( x k ) B_kp_k=-\nabla f(x_k) Bkpk=−∇f(xk)求出 p k p_k pk
- 一维搜索,求 λ k \lambda _k λk使得 f ( x k + λ p k ) f(x_k+\lambda p_k) f(xk+λpk)最小
- 置 x k + 1 = x k + λ k p k x_{k+1}=x_k+\lambda_kp_k xk+1=xk+λkpk
- 计算 ∇ f ( x k + 1 ) \nabla f(x_{k+1}) ∇f(xk+1),若 ∣ ∣ ∇ f ( x k + 1 ) ∣ ∣ < ϵ ||\nabla f(x_{k+1})|| \lt \epsilon ∣∣∇f(xk+1)∣∣<ϵ,则停止计算,得到近似解 x ∗ = x k + 1 x^*=x_{k+1} x∗=xk+1,否则按式子 B k + 1 = B k + y k y k T y 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^T\delta_k}--\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^TB_k\delta_k} Bk+1=Bk+yTδkykykT−−δkTBkδkBkδkδkTBk计算得到 B k + 1 B_{k+1} Bk+1
- 置k=k+1,转3
2.3 Broyden类算法
迭代公式(19)经过两次应用Sherman-Morrison公式可以得到
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
y
k
T
δ
k
T
y
k
.
.
.
.
.
.
.
.
.
(
20
)
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_ky_k^T}{\delta_k^Ty_k}.........(20)
Gk+1=(I−δkTykδkykT)Gk(I−δkTykδkykT)T+δkTykδkykT.........(20)
称为BFGS算法关于
G
k
G_k
Gk的迭代公式,将式子(20)得到的
G
k
+
1
G_{k+1}
Gk+1记作
G
B
F
G
S
G^{BFGS}
GBFGS,将DFP算法得到的G_{k+1}记作
G
D
F
P
G^{DFP}
GDFP,它们都满足拟牛顿条件,则它们的线性组合也同样满足拟牛顿条件,并且也是正定的。
G
k
+
1
=
a
G
D
F
P
+
(
1
−
a
)
G
B
F
G
S
.
.
.
.
.
.
.
.
.
.
(
21
)
G_{k+1}=aG^{DFP}+(1-a)G^{BFGS}..........(21)
Gk+1=aGDFP+(1−a)GBFGS..........(21)
其中
0
≤
a
≤
1
0\le a \le 1
0≤a≤1,这样就得到了一类拟牛顿算法称为Broyden类算法。