vSLAMNet（四）-最优化-牛顿法

发表于 | 分类于 ，

1. 概述

  • 是一种在实数域和复数域上近似求解方程的方法
  • 方法使用函数f (x)的泰勒级数的前面几项来寻找方程f (x) = 0的根
  • 最大的特点就在于它的收敛速度很快

2. 算法原理

将函数f(x)在$x_0$处应用泰勒展开得，$f(x) = f(x_0) + \frac{f’(x_0)}{1!}(x - x_0) + \frac{f’’(x_0)}{2!}(x-x_0)^2 + \cdots$，我们取前两项，$f(x) = f(x_0) + f’(x_0)(x-x_0)$。

对等式两边求导数得到，

$f’(x) = f’(x_0) + f’’(x_0)\Delta x$，

当f(x)取最小值时，$f’(x) = 0$，所以我们得到，

$0 = f’(x_0)+ f’’(x_0)\Delta x$，

$\Delta x = -\frac{f’(x_0)}{f’’(x_0)}$，

这样我们就得到了牛顿法的参数迭代更新公式如下，

$x_{n+1} = x_n- \frac{f’(x_n)}{f’’(x_n}$。

例子：

求解最优化问题，函数$\min_{x\in R^n}f(x)$，其中$x^*$为目标函数的极小值。

设f(x)具有二阶连续偏导数，若第k次的迭代值为$x^{(k)}$，则可将f(x)在$x^{(k)}$附近进行二阶泰勒展开，

$f(x) = f(x^{(k)}) + g_k^T(x-x^{(k)}) + \frac{1}{2}(x-x^{(k)})^TH(x^{(k)})(x-x^{(k)})$，

其中，$g_k = g(x^{(k)}) = \bigtriangledown f(x^{(k)})$是f(x)的梯度向量在点$x^{(k)}$，$H(x^{(k)})$是f(x)的海赛矩阵$H(x^{(k)}) = [\frac{\partial^2f}{\partial x_i\partial x_j}]_{nxn}$在点$x^{(k)}$的值。

函数f(x)有极值的必要条件是在极值点处一阶导数为0，即梯度向量为0。特别的当$H(x^{(k)})$是正定矩阵时，函数f(x)的极值为极小值。对上述泰勒公式再进行求导得到，

$\bigtriangledown f(x) = g_k + H_k(x - x^{(k)})$，

其中，$H_k = H(x^{(k)})$，那么，

$g_k + H_k(x^{(k+1)} - x^{(k)}) = 0$，

$x^{(k+1)} = x^{(k)} - H_k^{-1}g_k$。

令$H_kp_k = -g_k$，

那么得到迭代公式，

$x^{(k + 1)} = x^{(k)} + p_k$。
最终可在$\bigtriangledown f(x^*)= 0$处收敛。

算法伪代码：

输入: f(x)，梯度$g(x) = \bigtriangledown f(x)$，海赛矩阵H(x)，精度要求$\epsilon$

输出: f(x)的极小值点$x^*$

  1. 取初始值$x^{(0)}$，置$k = 0$
  2. 计算$g_k = g(x^{(k)})$
  3. 若$||g_k|| \leq \epsilon$，停止得到近似解$x^* = x^{(k)}$
  4. 计算$H_k = H(x^{(k)})$，并求$p_k$，$H_kp_k = -g_k$
  5. 置$x^{(k+1)} = x^{(k)} + p_k$
  6. 置$k = k+1$，转2

优缺点：

  • 二阶收敛，收敛速度快
  • 牛顿法是一种迭代算法，每一步都需要求解目标函数的Hessian矩阵的逆矩阵，计算比较复杂

3. 实例

用牛顿法求解$f(x) = \frac{3}{2}x_1^2 - x_1*x_2 - 2x_1$的极小值。

解：由题可得，

梯度为$\nabla f(x) = (3x_1 - x_2 -2, x_2 - x_1)^T$，

$H(x) = \left[\begin{matrix} 3 & -1 \\ -1 & 1 \end{matrix}\right]$是正定矩阵。

取$x^{(0)} = (0, 0)^T$，则

$\nabla f(x^{(0)}) = (-2, 0)^T$，$H(x^{(0)}) = \left[\begin{matrix} 3 & -1 \\ -1 & 1 \end{matrix}\right]$，$H(x^{(0)})^{-1} = \left[\begin{matrix} \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{3}{2} \end{matrix}\right]$。

$x^{(1)} = x^{(0)} - H(x^{(0)})^{-1} \nabla f(x^{0}) = (1, 1)^T$.

将$x^{(1)}$代入$\nabla f(x)$得$(0, 0)^T$，那么$x^{(1)} = (1, 1)^T$即为极小值点，那么极小值为$-1$。

4. 代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
"""
Newton法
Rosenbrock函数
函数 f(x)=100*(x(2)-x(1).^2).^2+(1-x(1)).^2
梯度 g(x)=(-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1)),200*(x(2)-x(1)^2))^(T)
"""

import numpy as np
import matplotlib.pyplot as plt

def jacobian(x):
    return np.array([-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0]),200*(x[1]-x[0]**2)])

def hessian(x):
    return np.array([[-400*(x[1]-3*x[0]**2)+2,-400*x[0]],[-400*x[0],200]])

X1=np.arange(-1.5,1.5+0.05,0.05)
X2=np.arange(-3.5,2+0.05,0.05)
[x1,x2]=np.meshgrid(X1,X2)
f=100*(x2-x1**2)**2+(1-x1)**2; # 给定的函数
plt.contour(x1,x2,f,20) # 画出函数的20条轮廓线


def newton(x0):

    print('初始点为:')
    print(x0,'\n')
    W=np.zeros((2,10**3))
    i = 1
    imax = 1000
    W[:,0] = x0 
    x = x0
    delta = 1
    alpha = 1

    while i<imax and delta>10**(-5):
        p = -np.dot(np.linalg.inv(hessian(x)),jacobian(x))
        x0 = x
        x = x + alpha*p
        W[:,i] = x
        delta = sum((x-x0)**2)
        print('第',i,'次迭代结果:')
        print(x,'\n')
        i=i+1
    W=W[:,0:i]  # 记录迭代点
    return W

x0 = np.array([-1.2,1])
W=newton(x0)

plt.plot(W[0,:],W[1,:],'g*',W[0,:],W[1,:]) # 画出迭代点收敛的轨迹
plt.show()

5. 参考

-梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现