前言
本篇继续优化理论的算法学习,牛顿法。
最速下降法
首先回顾上次提到的梯度下降法(其实就是最速下降法):通过求取多元函数在某个点处的梯度,沿着梯度的反方向前进,直到达到迭代停止条件。
对于多元函数(实值向量函数)
f
(
x
)
f(\bf x)
f(x),其在
x
0
\bf x_0
x0处的泰勒展开可表示为:
f
(
x
)
=
f
(
x
0
)
+
∇
f
(
x
0
)
T
(
x
−
x
0
)
+
1
2
(
x
−
x
0
)
T
H
(
x
0
)
(
x
−
x
0
)
f({\bf x}) = f({\bf x_0})+{\bf \nabla} f({\bf x_0})^T({\bf x-x_0}) + \frac {1}{2} ({\bf x-x_0})^T{\bf H(x_0)}({\bf x-x_0})
f(x)=f(x0)+∇f(x0)T(x−x0)+21(x−x0)TH(x0)(x−x0)
其中,
∇
,
H
\nabla, H
∇,H分别是函数梯度,海森矩阵。
梯度下降法直接取一阶梯度的反方向作为优化方向,因此称为最速下降法(每次迭代的方向都是下降最快的方向)。
牛顿法
回到多元泰勒展开:
f
(
x
)
=
f
(
x
0
)
+
∇
f
(
x
0
)
(
x
−
x
0
)
+
1
2
(
x
−
x
0
)
T
H
(
x
0
)
(
x
−
x
0
)
f({\bf x}) = f({\bf x_0})+{\bf \nabla} f({\bf x_0})({\bf x-x_0}) + \frac {1}{2} ({\bf x-x_0})^T{\bf H(x_0)}({\bf x-x_0})
f(x)=f(x0)+∇f(x0)(x−x0)+21(x−x0)TH(x0)(x−x0)
f
(
x
)
f(\bf x)
f(x)的极小值处的导数应当等于
0
\bf 0
0。现在取初始点
x
0
\bf x_0
x0,将
f
(
x
)
f(\bf x)
f(x)在
x
0
\bf x_0
x0处的展开式对
x
\bf x
x进行求导:
∇
f
(
x
0
)
T
+
(
x
−
x
0
)
T
H
(
x
0
)
=
0
∇
f
(
x
0
)
=
−
H
(
x
0
)
(
x
−
x
0
)
x
=
x
0
−
H
(
x
0
)
−
1
∇
f
(
x
0
)
{\bf \nabla} f({\bf x_0})^T+{(\bf x-x_0)}^T{\bf H(x_0)}={\bf 0} \\ {\bf \nabla}f({\bf x_0})=-{\bf H(x_0)}({\bf x-x_0}) \\ \quad \\ {\bf x} = {\bf x_0} - {{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}}
∇f(x0)T+(x−x0)TH(x0)=0∇f(x0)=−H(x0)(x−x0)x=x0−H(x0)−1∇f(x0)
这样就得到了新的
x
\bf x
x,这就是牛顿法。
可以给出更具体的牛顿法求极值过程:
- 确定实值向量函数 f ( x ) f(\bf x) f(x)的初始点 x 0 \bf x_0 x0,梯度 ∇ f {\bf \nabla} f ∇f的表达式,海森矩阵 H ( x ) \bf H(x) H(x)的表达式,终止条件 δ , ϵ \delta, \epsilon δ,ϵ
- 计算梯度 ∇ f ( x 0 ) {\bf \nabla}f({\bf x_0}) ∇f(x0),如果 ∣ ∣ ∇ f ( x 0 ) ∣ ∣ < δ ||{\bf \nabla}f({\bf x_0})||<\delta ∣∣∇f(x0)∣∣<δ,终止迭代,得到解 x = x 0 \bf x=x_0 x=x0;否则进入第3步
- 计算海森矩阵,并迭代 x ′ = x 0 − H ( x 0 ) − 1 ∇ f ( x 0 ) {\bf x'} = {\bf x_0} - {{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}} x′=x0−H(x0)−1∇f(x0)
- 计算 Δ f = f ( x ′ ) − f ( x 0 ) \Delta f= f({\bf x'}) - f(\bf x_0) Δf=f(x′)−f(x0),如果 Δ f < ϵ \Delta f < \epsilon Δf<ϵ,终止迭代,得到解 x = x ′ \bf x=x' x=x′;否则 x 0 = x ′ \bf x_0=x' x0=x′,回到第2步
牛顿法分析
现在考虑牛顿法的优化迭代表达式:
x
′
=
x
0
−
H
(
x
0
)
−
1
∇
f
(
x
0
)
{\bf x'} = {\bf x_0} - {{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}}
x′=x0−H(x0)−1∇f(x0)
将牛顿法的迭代方向与梯度方向做内积:
−
H
(
x
0
)
−
1
∇
f
(
x
0
)
⋅
∇
f
(
x
0
)
=
−
∇
f
(
x
0
)
T
H
(
x
0
)
−
T
∇
f
(
x
0
)
=
−
∇
f
(
x
0
)
T
H
(
x
0
)
−
1
∇
f
(
x
0
)
-{{\bf H(x_0)}}^{-1} {{\bf \nabla} f{\bf(x_0)}} \cdot {{\bf \nabla} f{\bf(x_0)}} \\ = -{{\bf \nabla} f{\bf(x_0)}}^T {{\bf H(x_0)}}^{-T}{{\bf \nabla} f{\bf(x_0)}} \\ = -{{\bf \nabla} f{\bf(x_0)}}^T {{\bf H(x_0)}}^{-1}{{\bf \nabla} f{\bf(x_0)}} \\
−H(x0)−1∇f(x0)⋅∇f(x0)=−∇f(x0)TH(x0)−T∇f(x0)=−∇f(x0)TH(x0)−1∇f(x0)
当海森矩阵是正定矩阵时,迭代方向与梯度方向的内积为负,即函数值在不断减小。这就涉及到了优化理论的重要问题:什么条件能够满足海森矩阵的正定性?
牛顿法的优点:优化收敛速度比梯度下降法更快。原因是牛顿法不仅考虑到函数的一阶梯度,也考虑到了海森矩阵(二阶梯度),在每个优化点处附近都能一次求得最佳迭代值。
牛顿法的缺点:每次迭代的速度较慢。这是因为迭代值的计算中涉及到了海森矩阵的逆,使得单步迭代效率比较低。
代码示例
和上次一样,我写了一个同牛顿法寻找二元函数极小值的代码,二元五次函数,如下所示:
import numpy as np
def partial_derivate_xy(x, y, F):
dx = (F(x + 0.001, y) - F(x, y))/0.001
dy = (F(x, y + 0.001) - F(x, y))/0.001
return dx, dy
def partial_partial_derivate_xy(x, y, F):
dx, dy = partial_derivate_xy(x, y, F)
dxx = (partial_derivate_xy(x + 0.001, y, F)[0] - dx) / 0.001
dyy = (partial_derivate_xy(x, y + 0.001, F)[1] - dy) / 0.001
dxy = (partial_derivate_xy(x, y + 0.001, F)[0] - dx) / 0.001
dyx = (partial_derivate_xy(x + 0.001, y, F)[1] - dy) / 0.001
return dxx, dyy, dxy, dyx
def non_linear_func(x, y):
fxy = 0.5 * (x ** 2 + y ** 2)
return fxy
def non_linear_func_2(x, y):
fxy = x*x + 2*y*y + 2*x*y + 3*x - y - 2
return fxy
def non_linear_func_3(x, y):
fxy = 0.5 * (x ** 2 - y ** 2)
return fxy
def non_linear_func_4(x, y):
fxy = x**4 + 2*y**4 + 3*x**2*y**2 + 4*x*y**2 + x*y + x + 2*y + 0.5
return fxy
def newton_optim(x, y, F):
dx, dy = partial_derivate_xy(x, y, F)
dxx, dyy, dxy, dyx = partial_partial_derivate_xy(x, y, F)
grad = np.array([[dx], [dy]])
hessian = np.array([[dxx, dxy], [dyx, dyy]])
hessian_inv = np.linalg.inv(hessian)
vec_delta = np.matmul(hessian_inv, grad)
vec_opt = np.array([[x], [y]]) - vec_delta
x_opt = vec_opt[0][0]
y_opt = vec_opt[1][0]
return x_opt, y_opt, vec_delta
def optimizer(x0, y0, F, th=0.01):
x = x0
y = y0
counter = 0
while True:
x_opt, y_opt, vec_delta = newton_optim(x, y, F)
if np.linalg.norm(vec_delta) < th:
break
x = x_opt
y = y_opt
counter = counter + 1
print('iter: {}'.format(counter), 'optimized (x, y) = ({}, {})'.format(x, y))
return x, y
if __name__ == '__main__':
x0 = 2.
y0 = 2.
result_x, result_y = optimizer(x0, y0, non_linear_func_4)
print(result_x, result_y)
牛顿法迭代11次后得到优化解。
梯度下降法在lr=0.01时迭代149次得到优化解。
运行1000次时,牛顿法用时0.586s,梯度下降法用时0.285s,可见虽然牛顿法迭代次数少,但速度慢。
后记
下篇将进入高斯牛顿法,后续还会有LM算法、拟牛顿法等。