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

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

拟牛顿法 Quasi-Newton Methods

为什么引入拟牛顿法

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

下面是公式层面的一个推导。
min ⁡ x f ( x ) \min_x f(x) xmin​f(x)

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

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

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

拟牛顿法基础

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

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

DFP形式推导

对上面的式中，我们知道存在
Δ x = − h t + 1 − 1 Δ g t \Delta x = -h_{t+1}^{-1} \Delta g_t Δx=−ht+1−1​Δgt​
这里我们构造一个矩阵 D t D_t Dt​来逼近这个函数，认为存在
Δ x = D t + 1 Δ g t \Delta x = D_{t+1} \Delta g_t Δx=Dt+1​Δgt​
这里注意， D t D_t Dt​是我们构造的一个矩阵，本身是不准确的，我们想要逐步迭代去逼近真实的 D t D_t Dt​，有点自举式算法的味道。因此我们令
D t + 1 = D t + Δ D t D_{t+1} = D_t + \Delta D_t Dt+1​=Dt​+ΔDt​
代入上面的式子，有
Δ x = D t Δ g t + Δ D t Δ g t \Delta x = D_t \Delta g_t + \Delta D_t \Delta g_t Δx=Dt​Δgt​+ΔDt​Δgt​
我们重点关注的是“自举”的过程，因此将上面的式子变式为
Δ D t Δ g t = Δ x − D t Δ g t \Delta D_t \Delta g_t = \Delta x - D_t \Delta g_t ΔDt​Δgt​=Δx−Dt​Δgt​
在这里我们假设存在一个向量 q t q_t qt​, w t w_t wt​，使得下面的式子成立：
Δ D t = Δ x q t T + D t Δ g t w t T \Delta D_t = \Delta x q_t^T + D_t \Delta g_t w_t^T ΔDt​=ΔxqtT​+Dt​Δgt​wtT​
注意hessian是一个对称矩阵，因此我们认为 Δ D t \Delta D_t ΔDt​也应该是对称的，又参照上面的两个式子，可以得到
q t T Δ g t = I q_t^T \Delta g_t = I qtT​Δgt​=I
w t T Δ g t = I w_t^T \Delta g_t = I wtT​Δgt​=I
不妨设
q t = α t Δ x q_t = \alpha_t \Delta x qt​=αt​Δx
w t = β t D t Δ g t w_t = \beta_t D_t \Delta g_t wt​=βt​Dt​Δgt​
带入到上面的式子，可得
α t = 1 Δ g t T Δ x \alpha_t = \frac{1}{\Delta g_t^T \Delta x} αt​=ΔgtT​Δx1​
β t = 1 Δ g t T D t Δ g t \beta_t= \frac{1}{\Delta g_t^T D_t \Delta g_t} βt​=ΔgtT​Dt​Δgt​1​
代入可得
Δ D t = Δ x t Δ x t T Δ g t T Δ x t − D t Δ g t Δ g t T D t T Δ g t T D t Δ 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 \Delta g_t} ΔDt​=ΔgtT​Δxt​Δxt​ΔxtT​​−ΔgtT​Dt​Δgt​Dt​Δgt​ΔgtT​DtT​​

BFGS方法推导

同理可得BFGS方法的推导。
设
Δ x t = B t + 1 − 1 Δ g t \Delta x_t = B_{{t+1}}^{-1} \Delta g_t Δxt​=Bt+1−1​Δgt​
因此也有
B t + 1 Δ x t = Δ g t B_{t+1} \Delta x_t = \Delta g_t Bt+1​Δxt​=Δgt​
同理，我们认为 B t + 1 = B t + Δ B t B_{t+1}=B_t + \Delta B_t Bt+1​=Bt​+ΔBt​
代入上式得到
B t Δ x t + Δ B t Δ x t = Δ g t B_t \Delta x_t + \Delta B_t \Delta x_t = \Delta g_t Bt​Δxt​+ΔBt​Δxt​=Δgt​
Δ B t Δ x t = Δ g t − B t Δ x t \Delta B_t \Delta x_t = \Delta g_t - B_t \Delta x_t ΔBt​Δxt​=Δgt​−Bt​Δxt​
同样的，有
Δ B t = Δ g t q t T − B t Δ x t w t T \Delta B_t = \Delta g_t q_t^T - B_t \Delta x_t w_t^T ΔBt​=Δgt​qtT​−Bt​Δxt​wtT​
观察可知存在
Δ x t q t T = I \Delta x_t q_t^T = I Δxt​qtT​=I
Δ x t w t T = I \Delta x_t w_t^T = I Δxt​wtT​=I
设
q t T = α t Δ g t T q_t^T = \alpha_t \Delta g_t^T qtT​=αt​ΔgtT​
w t T = β t Δ x t T B t T w_t^T = \beta_t \Delta x_t^T B_t^T wtT​=βt​ΔxtT​BtT​
因此有
α t = 1 Δ x t Δ g t T \alpha_t = \frac{1}{\Delta x_t \Delta g_t^T} αt​=Δxt​ΔgtT​1​
β t = 1 Δ x t T Δ B t T Δ x t \beta_t = \frac{1}{\Delta x_t^T \Delta B_t^T \Delta x_t} βt​=ΔxtT​ΔBtT​Δxt​1​
代入得到
Δ B t = Δ g t Δ g t T Δ x t Δ g t T − B t Δ x t Δ x t T B t T Δ x t T Δ B t T Δ x t \Delta B_t = \frac{\Delta g_t \Delta g_t^T}{\Delta x_t \Delta g_t^T} - \frac{B_t \Delta x_t \Delta x_t^T B_t^T}{\Delta x_t^T \Delta B_t^T \Delta x_t} ΔBt​=Δxt​ΔgtT​Δgt​ΔgtT​​−ΔxtT​ΔBtT​Δxt​Bt​Δxt​ΔxtT​BtT​​
B t + 1 − 1 = ( I n − Δ x Δ g T Δ x t T Δ g t ) B t − 1 ( I n − Δ g t Δ x t T Δ x t T Δ g t ) + Δ x t Δ x t T Δ x t T Δ g t B_{t+1}^{-1} = (I_n - \frac{\Delta x \Delta g^T}{\Delta x_t^T \Delta g_t})B_t^{-1}(I_n - \frac{\Delta g_t \Delta x_t^T}{\Delta x_t^T \Delta g_t}) + \frac{\Delta x_t \Delta x_t^T}{\Delta x_t^T \Delta g_t} Bt+1−1​=(In​−ΔxtT​Δgt​ΔxΔgT​)Bt−1​(In​−ΔxtT​Δgt​Δgt​ΔxtT​​)+ΔxtT​Δgt​Δxt​ΔxtT​​

实现代码

实现代码如下

import numpy as npdef dfp(f, x0, eps=1e-6, max_iter=100):"""DFP法最优化函数Args:f: 目标函数x0: 初始值eps: 精度max_iter: 最大迭代次数Returns:tuple: 最优化的结果"""n = len(x0)x = x0H = np.eye(n)grad = np.ones(n)k = 0while np.linalg.norm(grad) > eps and k < max_iter:grad = np.gradient(f, x)p = -np.dot(H, grad)alpha = 1while f(x + alpha * p) > f(x) + 0.5 * alpha * np.dot(grad, p):alpha = 0.5 * alphas = alpha * px_new = x + sy = np.gradient(f, x_new) - gradrho = 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_newH = H_newk += 1return x, f(x), k

import numpy as npdef bfgs(f, x0, eps=1e-6, max_iter=100):"""BFGS法最优化函数Args:f: 目标函数x0: 初始值eps: 精度max_iter: 最大迭代次数Returns:tuple: 最优化的结果"""n = len(x0)x = x0H = np.eye(n)grad = np.ones(n)k = 0while np.linalg.norm(grad) > eps and k < max_iter:grad = np.gradient(f, x)p = -np.dot(H, grad)alpha = 1while f(x + alpha * p) > f(x) + 0.5 * alpha * np.dot(grad, p):alpha = 0.5 * alphas = alpha * px_new = x + sy = np.gradient(f, x_new) - gradrho = 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_newH = H_newk += 1return x, f(x), k

Sherman-Morrison 公式

对于任意非奇异方阵$A$， u , v ∈ R n u,v \in R^n u,v∈Rn，若 1 + v T A − 1 u ≠ 0 1+v^TA^{-1}u \neq 0 1+vTA−1u=0则
( A + u v T ) − 1 = A − 1 − ( A − 1 u ) ( v T A − 1 ) 1 + v T A − 1 u (A+uv^T)^{-1} = A^{-1} - \frac{(A^{-1}u)(v^TA^{-1})}{1+v^TA^{-1}u} (A+uvT)−1=A−1−1+vTA−1u(A−1u)(vTA−1)​
该公式描述了在矩阵 $A$发生某种变化时，如何利用之前求好的逆，求新的逆。

对迭代公式引入两次 Sherman-Morrison 公式就能得到
B t + 1 − 1 = ( I n − Δ x t Δ g t T Δ x t T Δ g t ) B t − 1 ( I n − Δ g t Δ x t T Δ x t T Δ g t ) + Δ x t Δ x t T Δ x t T Δ g t B_{t+1}^{-1} = (I_n - \frac{\Delta x_t \Delta g_t^T}{\Delta x_t^T \Delta g_t})B_t^{-1}(I_n - \frac{\Delta g_t \Delta x_t^T}{\Delta x_t^T \Delta g_t}) + \frac{\Delta x_t \Delta x_t^T}{\Delta x_t^T \Delta g_t} Bt+1−1​=(In​−ΔxtT​Δgt​Δxt​ΔgtT​​)Bt−1​(In​−ΔxtT​Δgt​Δgt​ΔxtT​​)+ΔxtT​Δgt​Δxt​ΔxtT​​
之后有空会更新下面的一些算法。

• 凸且光滑的函数的BFGS优化算法
• 非凸但平滑的函数BFGS优化算法
• L-BFGS优化算法
• 非凸非平滑函数的BFGS优化算法

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！