机器人中的数值优化|【三】无约束优化，拟牛顿法，共轭梯度法理论与推导

拟牛顿法 Quasi-Newton Methods

为什么引入拟牛顿法

在前面的章节中，我们学习了牛顿法，牛顿法的核心是先通过将函数泰勒展开，近似为一个二阶项目，对这个二阶项求导，可以得到极值点，则直接找到了在函数展开点附近的最优点。注意，我们这里说的是函数展开点附近的最优点。因为泰勒展开存在截断误差，我们是不能认为该点就是精确解的。

下面是公式层面的一个推导。
$min_x f(x)$

对 $x$ ,我们于 $x_t$ 处（第 $t$ 次迭代的 $x$ ）位置进行二阶泰勒展开，有
$\approx f(x_t)+f^{(1)}(x_t)(x-x_t)+\frac{1}{2}f^{(2)}(x_t)(x-x_t)^2$
令 $x-x_t=\Delta x$
$\approx f(x_t)+f^{(1)}(x_t)\Delta x+\frac{1}{2}f^{(2)}(x_t)\Delta x^2$
$\hat f(x)= f(x_t)+f^{(1)}(x_t)\Delta x+\frac{1}{2}f^{(2)}(x_t)\Delta x^2$
很容易的，我们对上式求一阶导数可以得到极值点
$\hat f'(x)= f^{(1)}(x_t)+f^{(2)}(x_t)\Delta x=0$
即当 $\Delta x=-\frac{ f^{(1)}(x_t)}{f^{(2)}(x_t)}$ 时，有极值点。
此处 $x-x_t=\Delta x$ 可以认为是优下降的方向。令 $g_t=f^{(1)}(x_t)$ , $h_t=f^{(2)}(x_t)$ ，分别代表近似函数的梯度和hessian，那么牛顿法的迭代过程就可以表示为
$\Delta x = -h^{-1} g$
不难看出，牛顿法存在两个问题

hessian不一定存在
hessian的求逆比较复杂

因此，我们引入了拟牛顿法来解决相关问题。

拟牛顿法基础

在概念上，我们使用一个矩阵 $U$ 来近似hessian，在二次条件下，hessian满足如下条件
$x_{t+1}-x_t = h_t^{-1}(f'(x_{t+1})-f'(x_t))$
要求近似矩阵U也应当满足相应条件
$x_{t+1}-x_t = U_t(f'(x_{t+1})-f'(x_t))$
拟牛顿法的基础形式一般有DFP和BFGS法两种。DFP法和BFGS法都是求解无约束优化问题的二次型拟牛顿法，其核心思想是通过构建二次模型来近似原始函数，利用该模型求解最优解的方向和步长，从而迭代地逼近全局最优解。

具体来说，DFP法和BFGS法都通过逐步构建Hessian矩阵的逆矩阵来求解最优解，但它们的不同之处在于构建逆矩阵的方式不同。

DFP形式推导

对上面的式中，我们知道存在
$\Delta x = -h_t^{-1} \Delta g_t$
这里我们构造一个矩阵 $D_t$ 来逼近这个函数，认为存在
$\Delta x = D_t \Delta g_t$
这里注意， $D_t$ 是我们构造的一个矩阵，本身是不准确的，我们想要逐步迭代去逼近真实的 $D_t$ ，有点自举式算法的味道。因此我们令
$D_{t+1} = D_t + \Delta D_t$
代入上面的式子，有
$\Delta x = D_t \Delta g_t + \Delta D_t \Delta g_t$
我们重点关注的是“自举”的过程，因此将上面的式子变式为
$\Delta D_t \Delta g_t = \Delta x - D_t \Delta g_t$
在这里我们假设存在一个向量 $q_t$ , $w_t$ ，使得下面的式子成立：
$\Delta D_t = \Delta x q_t^T + D_t \Delta g_t w_t^T$
注意hessian是一个对称矩阵，因此我们认为 $\Delta D_t$ 也应该是对称的，又参照上面的两个式子，可以得到
$q_t^T \Delta g_t = I$
$w_t^T \Delta g_t = I$
不妨设
$q_t = \alpha_t \Delta x$
$w_t = \beta_t D_t \Delta g_t$
带入到上面的式子，可得
$\alpha_t = \frac{1}{\Delta x^T \Delta g_t}$
$\beta_t= \frac{1}{\Delta g_t^T D_t^T \Delta g_t}$
代入可得
$\Delta D_t = \frac{\Delta x_t \Delta x_t^T}{\Delta g_t^T \Delta x_t} - \frac{D_t \Delta g_t \Delta g_t^T D_t^T}{\Delta g_t^T D_t^T \Delta g_t}$

BFGS方法推导

同理可得BFGS方法的推导。

实现代码如下

import numpy as np

def dfp(f, x0, eps=1e-6, max_iter=100):
    """
    DFP法最优化函数

    Args:
        f: 目标函数
        x0: 初始值
        eps: 精度
        max_iter: 最大迭代次数

    Returns:
        tuple: 最优化的结果
    """
    n = len(x0)
    x = x0
    H = np.eye(n)
    grad = np.ones(n)
    k = 0
    while np.linalg.norm(grad) > eps and k < max_iter:
        grad = np.gradient(f, x)
        p = -np.dot(H, grad)
        alpha = 1
        while f(x + alpha * p) > f(x) + 0.5 * alpha * np.dot(grad, p):
            alpha = 0.5 * alpha
        s = alpha * p
        x_new = x + s
        y = np.gradient(f, x_new) - grad
        rho = 1 / np.dot(y, s)
        A = np.eye(n) - rho * np.outer(s, y)
        H_new = np.dot(A, np.dot(H, A.T)) + rho * np.outer(s, s)
        x = x_new
        H = H_new
        k += 1
    return x, f(x), k

import numpy as np

def bfgs(f, x0, eps=1e-6, max_iter=100):
    """
    BFGS法最优化函数

    Args:
        f: 目标函数
        x0: 初始值
        eps: 精度
        max_iter: 最大迭代次数

    Returns:
        tuple: 最优化的结果
    """
    n = len(x0)
    x = x0
    H = np.eye(n)
    grad = np.ones(n)
    k = 0
    while np.linalg.norm(grad) > eps and k < max_iter:
        grad = np.gradient(f, x)
        p = -np.dot(H, grad)
        alpha = 1
        while f(x + alpha * p) > f(x) + 0.5 * alpha * np.dot(grad, p):
            alpha = 0.5 * alpha
        s = alpha * p
        x_new = x + s
        y = np.gradient(f, x_new) - grad
        rho = 1 / np.dot(y, s)
        A = np.eye(n) - rho * np.outer(s, y)
        H_new = np.dot(A.T, np.dot(H, A)) + rho * np.outer(y, y)
        x = x_new
        H = H_new
        k += 1
    return x, f(x), k

之后有空会更新下面的一些算法。文章来源地址https://uudwc.com/A/Nx6Db

机器人中的数值优化|【三】无约束优化，拟牛顿法，共轭梯度法理论与推导