非线性最小二乘问题的数值方法 —— 狗腿法 Powell‘s Dog Leg Method (I - 原理与算法)

Title: 非线性最小二乘问题的数值方法 —— 狗腿法 Powell’s Dog Leg Method (I - 原理与算法)


文章目录

  • I. 前言
  • II. 线搜索类型和信赖域类型
    • 1. 线搜索类型 —— 最速下降法
    • 2. 信赖域类型
    • 3. 柯西点
  • III. 狗腿法的原理
    • 1. 狗腿法的构建
    • 2. 狗腿法的优化说明
    • 3. 狗腿法的插值权重
  • IV. 狗腿法的流程
    • 1. 狗腿法的信赖域控制
    • 2. 狗腿法的停止条件
      • 条件一. 梯度不再下降
      • 条件二. 迭代点不更新
      • 条件三. 残差足够小
      • 条件四. 达到最大迭代数
    • 3. 狗腿法的实现流程
  • IV. 总结
  • 参考文献


关联博客文章

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (I)

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (II)

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (III)

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (IV)

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (实例篇 V)

非线性最小二乘问题的数值方法 —— 从高斯-牛顿法到列文伯格-马夸尔特法 (I)

非线性最小二乘问题的数值方法 —— 从高斯-牛顿法到列文伯格-马夸尔特法 (II, Python 简单实例)


I. 前言

之前的博文涉及到最速下降法 (Steepest Descent Method, 也被称为梯度下降法)、高斯-牛顿法 (Gauss-Newton Method) 和列文伯格-马夸尔特法 (Levenberg-Marquardt Method).

列文伯格-马夸尔特法比起高斯-牛顿法法比较明显的一项优势不需要初始点靠近极小值点没有 Hessian 矩阵的正定性要求.

但是列文伯格-马夸尔特法也有其缺点, 比如计算过程中大量的 “不合格” 计算步被舍弃而作无用功.

下面要介绍的狗腿法 (Dog Leg Method 或 Powell’s Dog Leg Method) 正是要克服列文伯格-马夸尔特法低效计算问题, 不去弃置其中的复杂计算.

由之前博文中的介绍可以知道:

- 最速下降法可以获得稳定地收敛效果, 对初值不敏感;

- 高斯-牛顿法在极小值附近可以快速收敛到极小值, 即二阶局部收敛.

利用两种方法的特点, 取长补短可能获得更好地计算性能.

事实上, 列文伯格-马夸尔特法和狗腿法都融合了最速下降法和高斯-牛顿法.

列文伯格-马夸尔特法通过调节控制阻尼参数 μ \mu μ, 让该方法表现出最速下降法特性或/和高斯-牛顿法特性.

而下面的狗腿法通过控制信赖域 (Trust Region) 的大小来显式地调节最速下降法和高斯-牛顿法在狗腿法结果中的比重.

列文伯格-马夸尔特法隐性地平衡最速下降法与高斯-牛顿法; 而狗腿法显式地控制最速下降法与高斯-牛顿法.


II. 线搜索类型和信赖域类型

在求解如下无约束的非线性最小二乘问题时,
m i n i m i z e g ( x ) = 1 2 ∥ r ( x ) ∥ 2 2 = 1 2 ∑ i = 1 m r i ( x ) 2 (II-1) {\rm minimize}\quad {g}(\mathbf{x}) = \frac{1}{2}\|\mathbf{r}(\mathbf{x})\|_2^2 = \frac{1}{2}\sum_{i=1}^{m} r_i(\mathbf{x})^2 \tag{II-1} minimizeg(x)=21r(x)22=21i=1mri(x)2(II-1)
可能遇到各种各样的方法、改进和变种, 这些方法一般都可以被归类为线搜索类型和信赖域类型[1].

线搜索类型的算法都需要先确定搜索的方向 (比如梯度等), 再沿着这个搜索方向寻找下一迭代点. 最速下降法、高斯-牛顿法就属于线搜索类型.

信赖域类型的算法在一个给定的区域内, 使用二阶模型近似原问题, 通过不断利用二阶近似模型的最小值来迭代地逼近原问题的极小值点. 列文伯格-马夸尔特法、狗腿法就属于信赖域类型.


1. 线搜索类型 —— 最速下降法

线搜索类型算法中最典型的就是最速下降法. 利用最速下降法寻找代价函数/目标函数 g ( x ) g(\mathbf{x}) g(x) 最小值的过程, 就像是下山的过程. 先确定往哪个方向走可以下山, 即 g ( x ) g(\mathbf{x}) g(x) 下降方向. 再确定这一步走多远可以最快下山, g ( x ) g(\mathbf{x}) g(x) 下降最大的步长. 到了新的点继续寻找下山方向和最优步长, 周而复始直到到达极小值点.

最速下降法中, 当前迭代点记为 x [ i ] \mathbf{x}_{[i]} x[i], 当前步的搜索方向记为 s d h [ i ] ^{sd}\mathbf{h}_{[i]} sdh[i], 当前步的步长记为 α [ i ] ≥ 0 \alpha_{[i]}\geq 0 α[i]0, 则下一迭代点为
x [ i + 1 ] = x [ i ] + α [ i ] s d h [ i ] (II-1-1) \mathbf{x}_{[i+1]} = \mathbf{x}_{[i]} + \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}} \tag{II-1-1} x[i+1]=x[i]+α[i]sdh[i](II-1-1)
α [ i ] s d h [ i ] \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}} α[i]sdh[i] 为最速下降步.

定义 r ( x ) \mathbf{r}(\mathbf{x}) r(x) 相对于 x \mathbf{x} x 在迭代点 x [ i ] \mathbf{x}_{[i]} x[i] 的 Jacobian 矩阵
J r ( x [ i ] ) ≜ ∂ r ( x ) ∂ x ∣ x [ i ] (II-1-2) \mathbf{J}_r(\mathbf{x}_{[i]}) \triangleq \left.\frac{\partial\mathbf{r}(\mathbf{x})}{\partial \mathbf{x}} \right|_{\mathbf{x}_{[i]}} \tag{II-1-2} Jr(x[i])xr(x) x[i](II-1-2)
最速下降法的**搜索方向**是 g ( x ) g(\mathbf{x}) g(x) 的梯度下降最快的方向, 即为梯度向量的反向.
s d h [ i ] = − ∇ g ( x ) ∣ x [ i ] = − J r ( x [ i ] ) T r ( x [ i ] ) (II-1-3) ^{sd}\mathbf{h}_{[i]} =- \left. \nabla {g}(\mathbf{x}) \right|_{\mathbf{x}_{[i]}} = - \mathbf{J}_r(\mathbf{x}_{[i]}) ^{\rm\small T}\, \mathbf{r}(\mathbf{x}_{[i]}) \tag{II-1-3} sdh[i]=g(x)x[i]=Jr(x[i])Tr(x[i])(II-1-3)
上式只是找到了下山的最优方向, 在这个最优方向上设置多长的步长才可以可让目标函数下降最多呢?

下面需要求**最优步长** α [ i ] \alpha_{[i]} α[i], 使得下一迭代步上的代价函数最小, 即下山最快.

对二次连续可微的函数 r ( x [ i ] + α [ i ] s d h [ i ] ) \mathbf{r}(\mathbf{x}_{[i]} + \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) r(x[i]+α[i]sdh[i]) 作一阶泰勒近似为
r ( x [ i ] + α [ i ] s d h [ i ] ) ≈ r ( x [ i ] ) + α [ i ] J r ( x [ i ] ) s d h [ i ] (II-1-4) \mathbf{r}(\mathbf{x}_{[i]} + \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) \approx \mathbf{r}(\mathbf{x}_{[i]}) + \alpha_{[i]}\, \mathbf{J}_{r}(\mathbf{x}_{[i]})\, {^{sd}\mathbf{h}_{[i]}} \tag{II-1-4} r(x[i]+α[i]sdh[i])r(x[i])+α[i]Jr(x[i])sdh[i](II-1-4)
其对应的目标函数的二阶近似为
g ( x [ i ] + α [ i ] s d h [ i ] ) ≈ 1 2 ∥ r ( x [ i ] + α [ i ] s d h [ i ] ) ∥ 2 = 1 2 [ r ( x [ i ] ) + α [ i ] J r ( x [ i ] ) s d h [ i ] ] T [ r ( x [ i ] ) + α [ i ] J r ( x [ i ] ) s d h [ i ] ] = g ( x [ i ] ) + α [ i ] s d h [ i ] T J r ( x [ i ] ) T r ( x [ i ] ) + 1 2 α [ i ] 2 ∥ J r ( x [ i ] ) s d h [ i ] ∥ 2 ≜ L ( α [ i ] s d h [ i ] ) (II-1-5) \begin{aligned} g(\mathbf{x}_{[i]} + \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) &\approx \frac{1}{2}\|\mathbf{r}(\mathbf{x}_{[i]} + \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}})\|^2\\ &= \frac{1}{2} \left[ \mathbf{r}(\mathbf{x}_{[i]}) + \alpha_{[i]}\, \mathbf{J}_{r}(\mathbf{x}_{[i]})\, {^{sd}\mathbf{h}_{[i]}} \right]^{\small\rm T} \left[ \mathbf{r}(\mathbf{x}_{[i]}) + \alpha_{[i]}\, \mathbf{J}_{r}(\mathbf{x}_{[i]})\, {^{sd}\mathbf{h}_{[i]}} \right]\\ &= g(\mathbf{x}_{[i]}) + \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}^{\small\rm T}} \mathbf{J}_r(\mathbf{x}_{[i]})^{\small\rm T} \mathbf{r}(\mathbf{x}_{[i]}) + \frac{1}{2} \alpha_{[i]}^2 \left\|{\mathbf{J}_r(\mathbf{x}_{[i]}){^{sd}\mathbf{h}_{[i]}}}\right\|^2 \\ &\triangleq L(\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) \end{aligned}\tag{II-1-5} g(x[i]+α[i]sdh[i])21r(x[i]+α[i]sdh[i])2=21[r(x[i])+α[i]Jr(x[i])sdh[i]]T[r(x[i])+α[i]Jr(x[i])sdh[i]]=g(x[i])+α[i]sdh[i]TJr(x[i])Tr(x[i])+21α[i]2 Jr(x[i])sdh[i] 2L(α[i]sdh[i])(II-1-5)
上式看作是以 α [ i ] \alpha_{[i]} α[i] 为自变量的一元二次函数, 该函数一阶导数为零时可求得最优的步长.
α [ i ] = − s d h [ i ] T J r ( x [ i ] ) T r ( x [ i ] ) ∥ J r ( x [ i ] ) s d h [ i ] ∥ 2 = ∥ s d h [ i ] ∥ 2 ∥ J r ( x [ i ] ) s d h [ i ] ∥ 2 (II-1-6) \begin{aligned} \alpha_{[i]} & = - \frac{{^{sd}\mathbf{h}_{[i]}^{\small\rm T}} \mathbf{J}_r(\mathbf{x}_{[i]})^{\small\rm T} \mathbf{r}(\mathbf{x}_{[i]}) } {\left\|{\mathbf{J}_r(\mathbf{x}_{[i]}) {^{sd}\mathbf{h}_{[i]}}}\right\|^2 }\\ &=\frac{\left\| {^{sd}\mathbf{h}_{[i]}} \right\|^2}{\left\|{\mathbf{J}_r(\mathbf{x}_{[i]}) \,{^{sd}\mathbf{h}_{[i]}}}\right\|^2} \end{aligned} \tag{II-1-6} α[i]= Jr(x[i])sdh[i] 2sdh[i]TJr(x[i])Tr(x[i])= Jr(x[i])sdh[i] 2 sdh[i] 2(II-1-6)
式 (II-1-3) 和式 (II-1-6) 就构成了最速下降法的搜索方向和步长, 可通过式 (II-1-1) 得到下一迭代点[2].

注意此处的 α [ i ] \alpha_{[i]} α[i] 可称为无约束 α [ i ] \alpha_{[i]} α[i], 区别于后面将要介绍的信赖域约束下的对应值.


2. 信赖域类型

线搜索类型方法中将每一步迭代计算分成两部分, 搜索方向 s d h [ i ] ^{sd}\mathbf{h}_{[i]} sdh[i] 的计算和最优步长 α [ i ] \alpha_{[i]} α[i] 的计算. 而信赖域类型方法中, 直接在给定的有界区域 (信赖区域) Δ [ i ] \Delta_{[i]} Δ[i] 内优化计算迭代步 t r h [ i ] ^{tr}\mathbf{h}_{[i]} trh[i]. 有对应关系
α [ i ] s d h [ i ] ⇔ t r h [ i ] (II-2-1) \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}} \;\;\Leftrightarrow \;\; {^{tr}\mathbf{h}_{[i]}} \tag{II-2-1} α[i]sdh[i]trh[i](II-2-1)
也就是说 t r h [ i ] ^{tr}\mathbf{h}_{[i]} trh[i] 既包含迭代步的方向又包含迭代步的步长.

类比式 (II-1-5) 可知, 信赖域类型的目标函数的二阶近似为
g ( x [ i ] + t r h [ i ] ) ≈ 1 2 ∥ r ( x [ i ] + t r h [ i ] ) ∥ 2 = 1 2 [ r ( x [ i ] ) + J r ( x [ i ] ) t r h [ i ] ] T [ r ( x [ i ] ) + J r ( x [ i ] ) t r h [ i ] ] = g ( x [ i ] ) + r ( x [ i ] ) T J r ( x [ i ] ) t r h [ i ] + 1 2 t r h [ i ] T J r ( x [ i ] ) T J r ( x [ i ] ) t r h [ i ] ≜ L ( t r h [ i ] ) (II-2-2) \begin{aligned} g(\mathbf{x}_{[i]} + {^{tr}\mathbf{h}_{[i]}}) &\approx \frac{1}{2}\|\mathbf{r}(\mathbf{x}_{[i]} + {^{tr}\mathbf{h}_{[i]}})\|^2\\ &= \frac{1}{2} \left[ \mathbf{r}(\mathbf{x}_{[i]}) + \mathbf{J}_{r}(\mathbf{x}_{[i]})\, {^{tr}\mathbf{h}_{[i]}} \right]^{\small\rm T} \left[ \mathbf{r}(\mathbf{x}_{[i]}) + \mathbf{J}_{r}(\mathbf{x}_{[i]})\, {^{tr}\mathbf{h}_{[i]}} \right]\\ &= g(\mathbf{x}_{[i]}) + \mathbf{r}(\mathbf{x}_{[i]})^{\small\rm T} \, \mathbf{J}_r(\mathbf{x}_{[i]})\, {^{tr}\mathbf{h}_{[i]}} + \frac{1}{2} {^{tr}\mathbf{h}_{[i]}^{\small\rm T}} \,{\mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]}) \,{^{tr}\mathbf{h}_{[i]}}}\\ & \triangleq L(^{tr}\mathbf{h}_{[i]}) \end{aligned}\tag{II-2-2} g(x[i]+trh[i])21r(x[i]+trh[i])2=21[r(x[i])+Jr(x[i])trh[i]]T[r(x[i])+Jr(x[i])trh[i]]=g(x[i])+r(x[i])TJr(x[i])trh[i]+21trh[i]TJr(x[i])TJr(x[i])trh[i]L(trh[i])(II-2-2)
上式作为目标函数 g ( x [ i ] + t r h [ i ] ) g(\mathbf{x}_{[i]} + {^{tr}\mathbf{h}_{[i]}}) g(x[i]+trh[i]) 的在迭代点 x [ i ] \mathbf{x}_{[i]} x[i] 附近的二阶近似模型, 其实受到近似区域的限制. 只有较小的范围该二阶近似模型 L ( h [ i ] ) L(\mathbf{h}_{[i]}) L(h[i]) 才相对准确地描述 g ( x [ i ] + t r h [ i ] ) g(\mathbf{x}_{[i]} + {^{tr}\mathbf{h}_{[i]}}) g(x[i]+trh[i]). 即对近似模型增加的约束条件为
∥ t r h [ i ] ∥ ≤ Δ [ i ] (II-2-3) \left\| {^{tr}\mathbf{h}_{[i]}} \right\| \leq \Delta_{[i]} \tag{II-2-3} trh[i] Δ[i](II-2-3)
其中 Δ [ i ] \Delta_{[i]} Δ[i] 就被称为信赖域半径. 这样信赖域类型的方法就转变为每一步求解如下子问题[1]
min ⁡ t r h [ i ] L ( t r h [ i ] ) s.t. ∥ t r h [ i ] ∥ ≤ Δ [ i ] (II-2-4) \begin{aligned} {\min_{^{tr}\mathbf{h}_{[i]}}} &\quad L(^{tr}\mathbf{h}_{[i]})\\ \text{s.t.} &\quad \left\| {^{tr}\mathbf{h}_{[i]}} \right\| \leq \Delta_{[i]} \end{aligned}\tag{II-2-4} trh[i]mins.t.L(trh[i]) trh[i] Δ[i](II-2-4)

即使是线搜索类型中的最速下降法, 在最优计算过程中也利用了原目标函数的一阶或者二阶近似, 故也存在范围过大而使得近似程度降低, 最终影响最优计算结果的情况. 所以说信赖域的引入, 保证了近似模型对原目标函数的近似程度, 是有必要的.


3. 柯西点

根据信赖域的思想, 在无约束的最速下降法中增加截断步长限制, 即
∥ α ˉ [ i ] s d h [ i ] ∥ ≤ Δ [ i ] (II-3-1) \left \| {{\bar\alpha}_{[i]} {^{sd}\mathbf{h}_{[i]}}} \right \| \leq \Delta_{[i]} \tag{II-3-1} αˉ[i]sdh[i] Δ[i](II-3-1)
就可以得到柯西点的定义.

柯西点[1]

L ( α ˉ [ i ] s d h [ i ] ) L({\bar\alpha}_{[i]} {^{sd}\mathbf{h}_{[i]}}) L(αˉ[i]sdh[i]) g ( x ) g(\mathbf{x}) g(x) 在点 x [ i ] \mathbf{x}_{[i]} x[i] 处的二阶近似. 常数 α ˉ [ i ] \bar\alpha_{[i]} αˉ[i] 为如下优化问题的解
min ⁡ L ( α ˉ [ i ] s d h [ i ] ) s.t. ∥ α ˉ [ i ] s d h [ i ] ∥ ≤ Δ [ i ] α ˉ [ i ] ≥ 0 \begin{aligned} \min\quad& L({\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}})\\ \text{s.t.}\quad & \left \| {\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} \right \| \leq \Delta_{[i]}\\ &{\bar\alpha_{[i]}\geq 0} \end{aligned} mins.t.L(αˉ[i]sdh[i]) αˉ[i]sdh[i] Δ[i]αˉ[i]0
则称 c x [ i ] = x [ i ] + α ˉ [ i ] s d h [ i ] ^{c}\mathbf{x}_{[i]} = \mathbf{x}_{[i]} + {\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} cx[i]=x[i]+αˉ[i]sdh[i] 为柯西点. (其中 α ˉ [ i ] s d h [ i ] {\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} αˉ[i]sdh[i] 为柯西点对应的迭代步 —— 柯西步)

由式 (II-1-5) 可知
L ( α ˉ [ i ] s d h [ i ] ) = g ( x [ i ] ) ⏟ 常数项 + [ − α ˉ [ i ] s d h [ i ] T s d h [ i ] ] ⏟ 第二项 ≤ 0 + 1 2 α ˉ [ i ] 2 s d h [ i ] T J r ( x [ i ] ) T J r ( x [ i ] ) s d h [ i ] ⏟ 第三项 (II-3-2) L(\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) = \underbrace{g(\mathbf{x}_{[i]})}_{常数项} + \underbrace{ \left[-\bar \alpha_{[i]} {^{sd}\mathbf{h}_{[i]}^{\small\rm T}} {^{sd}\mathbf{h}_{[i]}}\right]}_{\color{green}第二项 \,\leq\, 0} + \underbrace{\frac{1}{2} {\bar\alpha}_{[i]}^2 {^{sd}\mathbf{h}_{[i]}^{\small\rm T}} \mathbf{J}_r(\mathbf{x}_{[i]})^{\small\rm T} {\mathbf{J}_r(\mathbf{x}_{[i]}){^{sd}\mathbf{h}_{[i]}}}}_{\color{blue}第三项} \tag{II-3-2} L(αˉ[i]sdh[i])=常数项 g(x[i])+第二项0 [αˉ[i]sdh[i]Tsdh[i]]+第三项 21αˉ[i]2sdh[i]TJr(x[i])TJr(x[i])sdh[i](II-3-2)
如果 J r ( x [ i ] ) T J r ( x [ i ] ) \mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]}) Jr(x[i])TJr(x[i]) 半负定的情况:

首先其中第二项 ≤ 0 \leq 0 0. 然后因为负半定, 故第三项 ≤ 0 \leq 0 0. 则可知 ∥ α ˉ [ i ] s d h [ i ] ∥ \|{\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}}\| αˉ[i]sdh[i] 越大, L ( α ˉ [ i ] s d h [ i ] ) L(\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) L(αˉ[i]sdh[i]) 越小. 故此时约束情况下的极小值在信赖域边界上取得, 即 ∥ α ˉ [ i ] s d h [ i ] ∥ = Δ [ i ] \left \| {\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} \right \| = \Delta_{[i]} αˉ[i]sdh[i] =Δ[i], 所以
α ˉ [ i ] = Δ [ i ] ∥ s d h [ i ] ∥ (II-3-3) \bar \alpha_{[i]} = \frac{\Delta_{[i]}}{\|^{sd}\mathbf{h}_{[i]}\|} \tag{II-3-3} αˉ[i]=sdh[i]Δ[i](II-3-3)
需要说明形如 J r ( x [ i ] ) T J r ( x [ i ] ) \mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]}) Jr(x[i])TJr(x[i]) 的对称矩阵, 不可能负定. 当 J r ( x [ i ] ) T J r ( x [ i ] ) = 0 \mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]}) = \mathbf{0} Jr(x[i])TJr(x[i])=0 时, 式 (II-3-3) 结论仍然成立.

如果 J r ( x [ i ] ) T J r ( x [ i ] ) \mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]}) Jr(x[i])TJr(x[i]) 正定的情况:

进一步假设, 如果 L ( α ˉ [ i ] s d h [ i ] ) L(\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) L(αˉ[i]sdh[i]) 的极小值点在信赖域内, 此时没有触发约束作用, 则有约束 α ˉ [ i ] \bar\alpha_{[i]} αˉ[i] 值和无约束 α [ i ] \alpha_{[i]} α[i] 值一致, 即式 (II-1-6) 所示.

进一步假设, 如果 L ( α ˉ [ i ] s d h [ i ] ) L(\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) L(αˉ[i]sdh[i]) 的极值点在信赖域外, 因为 L ( α ˉ [ i ] s d h [ i ] ) L(\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}) L(αˉ[i]sdh[i]) 为凸函数的缘故, 此时约束域 (信赖域) 内的极小值也在信赖域边界上取得, 如式 (II-3-3).

也就是说, 如果 J r ( x [ i ] ) T J r ( x [ i ] ) \mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]}) Jr(x[i])TJr(x[i]) 正定, 有约束 α ˉ [ i ] \bar\alpha_{[i]} αˉ[i] 值为
α ˉ [ i ] = min ⁡ { ∥ s d h [ i ] ∥ 2 ∥ J r ( x [ i ] ) s d h [ i ] ∥ 2 , Δ [ i ] ∥ s d h [ i ] ∥ } (II-3-4) \bar\alpha_{[i]} = \min \left\{ \frac{\left\| {^{sd}\mathbf{h}_{[i]}} \right\|^2}{\left\|{\mathbf{J}_r(\mathbf{x}_{[i]}) \,{^{sd}\mathbf{h}_{[i]}}}\right\|^2} ,\frac{\Delta_{[i]}}{\|^{sd}\mathbf{h}_{[i]}\|} \right\} \tag{II-3-4} αˉ[i]=min{ Jr(x[i])sdh[i] 2 sdh[i] 2,sdh[i]Δ[i]}(II-3-4)
以上计算柯西点或对应迭代步的示意图如下.

目标函数的极小值点在信赖域 J r ( x [ i ] ) T J r ( x [ i ] ) = 0 \mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]})=\mathbf{0} Jr(x[i])TJr(x[i])=0 或目标函数的极小值点在信赖域
cauchy_point_1cauchy_point_2

III. 狗腿法的原理

1. 狗腿法的构建

有了上面的铺垫, 我们可以正式引入狗腿法了.

狗腿法是高斯-牛顿法和最速下降法的融合, 重点在于信赖域对融合结果的影响, 基本原理可用如下三张图说明.

分类情况与说明示意图
Case I: 高斯-牛顿法的目标函数的极小值在信赖域
如果高斯-牛顿法获得的下一迭代步在信赖域内部, 则狗腿法的当前步直接采用高斯-牛顿法的结果 g n h [ i ] ^{gn}\mathbf{h}_{[i]} gnh[i], 以取得更快的收敛效果.
DC-motor-1
Case II: 无约束最速下降法的目标函数的极小值在信赖域
排除高斯-牛顿法获得在信赖域内的迭代步的情况下 (即排除 Case I 之后), 如果无约束最速下降法的迭代步在信赖域之外, 那么柯西点只能在梯度相反方向与信赖域边界的交点上取得.
这种情况下, 使用柯西点结果 α ˉ [ i ] s d h [ i ] {\bar\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} αˉ[i]sdh[i] 作为狗腿法的当前步.
DC-motor-1
Case III: 其他 (高斯-牛顿法迭代步目标函数的极小值在信赖域, 而无约束最速下降法迭代步目标函数的极小值在信赖域)
这种情况下的柯西点就是无约束最速下降法的迭代步. 此时狗腿法的路径是先沿着梯度反方向到达柯西点, 然后折向在信赖域外的高斯-牛顿步, 交于信赖域边缘, 这个交点就是该情况下的狗腿步.
也就是说, 该情况下的狗腿法结果是最速下降法和高斯-牛顿法之间的线性插值, 这个线性插值的权重分配决定于信赖域,
α [ i ] s d h [ i ] + β [ i ] ( g n h [ i ] − α [ i ] s d h [ i ] ) {\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} + \beta_{[i]} ( ^{gn}\mathbf{h}_{[i]} -{\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} ) α[i]sdh[i]+β[i](gnh[i]α[i]sdh[i])
其中 β [ i ] \beta_{[i]} β[i] 即为两类结果线性插值的权重.
信赖域较大的话, 高斯-牛顿步权重较大;
信赖域大到包含高斯-牛顿步的话, 就是 Case I;
信赖域较小的话, 最速下降步权重较大;
信赖域小到无法包含无约束最速下降步的话, 就是 Case II.
DC-motor-1

先把狗腿法原理介绍中基于信赖域的分类条件写成数学形式
Case I ⇔ ∥ g n h [ i ] ∥ ≤ Δ [ i ] Case II ⇔ ∥ α [ i ] s d h [ i ] ∥ ≥ Δ [ i ] Case III ⇔ ∥ g n h [ i ] ∥ > Δ [ i ] & ∥ α [ i ] s d h [ i ] ∥ < Δ [ i ] (III-1-1) \begin{aligned} \text{Case I} \quad &\Leftrightarrow \quad \|{^{gn}\mathbf{h}_{[i]}}\| \leq \Delta_{[i]}\\ \text{Case II} \quad &\Leftrightarrow \quad \| {\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}}\| \geq \Delta_{[i]}\\ \text{Case III} \quad &\Leftrightarrow \quad \|{^{gn}\mathbf{h}_{[i]}}\| > \Delta_{[i]} \quad \& \quad \| {\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}}\| < \Delta_{[i]} \end{aligned} \tag{III-1-1} Case ICase IICase IIIgnh[i]Δ[i]α[i]sdh[i]Δ[i]gnh[i]>Δ[i]&α[i]sdh[i]<Δ[i](III-1-1)
再由上面的原理介绍, 我们可知狗腿法的迭代步
d l h [ i ] = { g n h [ i ] , Case I Δ [ i ] ∥ s d h [ i ] ∥ s d h [ i ] , Case II α [ i ] s d h [ i ] + β [ i ] ( g n h [ i ] − α [ i ] s d h [ i ] ) , Case III (III-1-2) {^{dl}\mathbf{h}_{[i]}} = \left\{ {\begin{array}{ll} {^{gn}\mathbf{h}_{[i]}}, & \text{Case I}\\ {\frac{\Delta_{[i]}}{\|^{sd}\mathbf{h}_{[i]}\|} {^{sd}\mathbf{h}_{[i]}}}, & \text{Case II}\\ {\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} + \beta_{[i]} ( ^{gn}\mathbf{h}_{[i]} -{\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}} ), & \text{Case III} \end{array}} \right.\tag{III-1-2} dlh[i]= gnh[i],sdh[i]Δ[i]sdh[i],α[i]sdh[i]+β[i](gnh[i]α[i]sdh[i]),Case ICase IICase III(III-1-2)
其中高斯-牛顿法的迭代步 g n h [ i ] {^{gn}\mathbf{h}_{[i]}} gnh[i] 的计算
g n h [ i ] = − ( H ~ ( x [ i ] ) ) − 1 ∇ g ( x [ i ] ) = − ( [ ∂ r ( x ) ∂ x ] T ∂ r ( x ) ∂ x ∣ x [ i ] ) − 1 ∇ g ( x [ i ] ) = − [ J r ( x [ i ] ) T J r ( x [ i ] ) ] − 1 ∇ g ( x [ i ] ) (III-1-3) \begin{aligned} {^{gn}\mathbf{h}_{[i]}} &= - \left( \widetilde{\mathbf{H}}(\mathbf{x}_{[i]}) \right)^{-1}\, \nabla g(\mathbf{x}_{[i]})\\ &= - \left(\left.\left[ \frac{\partial\mathbf{r}(\mathbf{x})}{\partial \mathbf{x}} \right]^{\rm\small T} \frac{\partial\mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}\right|_{\mathbf{x}_{[i]}}\right)^{-1} \, \nabla g(\mathbf{x}_{[i]})\\ &= - \left[ \mathbf{J}_r(\mathbf{x}_{[i]})^{\small\rm T} \mathbf{J}_r(\mathbf{x}_{[i]}) \right]^{-1} \, \nabla g(\mathbf{x}_{[i]})\\ \end{aligned}\tag{III-1-3} gnh[i]=(H (x[i]))1g(x[i])= [xr(x)]Txr(x) x[i] 1g(x[i])=[Jr(x[i])TJr(x[i])]1g(x[i])(III-1-3)
因为需要计算高斯-牛顿法的迭代步, 由该方法的使用条件可知, 在狗腿法中也要求 J r ( x [ i ] ) T J r ( x [ i ] ) \mathbf{J}_r(\mathbf{x}_{[i]})^{\rm\small T} \mathbf{J}_r(\mathbf{x}_{[i]}) Jr(x[i])TJr(x[i])正定.


2. 狗腿法的优化说明

[1] Case I 和 Case II 的优化特性

狗腿法中 Case I 和 Case II 各自的最优化特性我们已经在 “高斯-牛顿法的优化观点” 和上面的 “柯西点” 中进行了说明.

再考虑到下面的结论 (结论 1)

- 高斯-牛顿步长度总是大于等于柯西步长度;

比较自然可以得到 Case I 和 Case II.

[2] Case III 的优化特性

但是 Case III 呢?

Case III “线性插值” 结构应该是构建性的, 基于这些历史留名的数学家的直觉.

那为什么也要到信赖域的边缘上才取得最优性能呢?

如果考虑到下面的引理 (本篇博文不展开讨论和证明)

- Case III 中, 迭代步长是插值权重参数 β \beta β 单调递增函数; (可以推出前面的 “结论 1”)

- Case III 中, 关于迭代步的近似模型函数是插值权重参数 β \beta β 的单调递减函数.

我们自然希望 Case III 中步长越长越好, 当然最大就是到达信赖域边缘.

以上作为对狗腿法优化特性的不严格说明. 下面分析插值权重参数 β \beta β 的具体计算.


3. 狗腿法的插值权重

下面需要进一步确定 Case III 中的插值权重系数 β [ i ] \beta_{[i]} β[i]​. 为了简化书写, 定义如下符号
a ≜ g n h [ i ] b ≜ α [ i ] s d h [ i ] c ≜ b T ( a − b ) β ≜ β [ i ] Δ ≜ Δ [ i ] (III-3-1) \begin{aligned} \mathbf{a} &\triangleq {^{gn}\mathbf{h}_{[i]}}\\ \mathbf{b} &\triangleq {\alpha_{[i]} {^{sd}\mathbf{h}_{[i]}}}\\ c & \triangleq \mathbf{b}^{\small\rm T} \left(\mathbf{a}- \mathbf{b}\right)\\ \beta & \triangleq \beta_{[i]} \\ \Delta & \triangleq \Delta_{[i]} \end{aligned} \tag{III-3-1} abcβΔgnh[i]α[i]sdh[i]bT(ab)β[i]Δ[i](III-3-1)
(参考了文献 [2], 但此处与文献中的 a \mathbf{a} a b \mathbf{b} b 的记号反过来了, 因为觉得高斯-牛顿法作用更大)

根据 Case III 中狗腿步在信赖域边缘上, 得到约束条件
∥ b + β ( a − b ) ∥ 2 = Δ 2 ⇒ ∥ a − b ∥ 2 β 2 + 2 c β + ∥ b ∥ 2 − Δ 2 = 0 (III-3-2) \| \mathbf{b} + \beta ( \mathbf{a} -\mathbf{b})\|^2 = \Delta^{2} \\ \Rightarrow \qquad \|\mathbf{a} - \mathbf{b}\|^2 {\color{blue}{\beta}^2} + 2 c {\color{blue}\beta} +\|\mathbf{b}\|^2 - \Delta^2 = 0 \tag{III-3-2} b+β(ab)2=Δ2ab2β2+2cβ+b2Δ2=0(III-3-2)
计算一元二次方程式 (III-3-2) 就能获得Case III 中的插值权重系数 β \beta β.

定义一元二次方程对应的一元二次函数
ψ ( β ) = ∥ a − b ∥ 2 β 2 + 2 c β + ∥ b ∥ 2 − Δ 2 (III-3-3) \psi(\beta) = \|\mathbf{a} - \mathbf{b}\|^2 {{\beta}^2} + 2 c {\beta} +\|\mathbf{b}\|^2 - \Delta^2 \tag{III-3-3} ψ(β)=ab2β2+2cβ+b2Δ2(III-3-3)
因为函数 ψ ( β ) \psi(\beta) ψ(β) 的二次项系数大于零 ( ∥ a − b ∥ 2 > 0 \|\mathbf{a} - \mathbf{b}\|^2 > 0 ab2>0), 故函数开口向上.

因为在 Case III 中 ∥ b ∥ 2 − Δ 2 < 0 \|\mathbf{b}\|^2-\Delta^2 < 0 b2Δ2<0, 故函数有两个根.

在 Case III 中,

β → − ∞ \beta \rightarrow -\infty β 时, ψ ( β ) → + ∞ \psi(\beta) \rightarrow +\infty ψ(β)+;

β = 0 \beta=0 β=0 时, ψ ( 0 ) < 0 \psi(0) < 0 ψ(0)<0;

β = 1 \beta = 1 β=1 时, ψ ( β ) = ∥ a ∥ 2 − Δ 2 > 0 \psi(\beta) = \|\mathbf{a}\|^2 - \Delta^2 > 0 ψ(β)=a2Δ2>0 (由式 (III-3-2) 中原始约束条件推到得到)

结合零点定理可知, ψ ( β ) \psi(\beta) ψ(β) 的两个根分布于 ( − ∞ , 0 ) (-\infty, 0) (,0) ( 0 , 1 ) (0,1) (0,1) 这两个区间内[2], 如下图所示.

函数示意
DC-motor-1

首先需要说明 0 < β < 1 0 < \beta < 1 0<β<1, 不然的话 Case III 中的狗腿步就不在柯西点和高斯-牛顿步的中间进行插值了, 即 0 < β < 1 0 < \beta < 1 0<β<1 这是由融合本身的要求决定的. 那么两个根中只要排除 ( − ∞ , 0 ) (-\infty,0) (,0) 区间内的根, 留下的另一就是可行解.

由一元二次方程求根公式
β = − c ± c 2 − ∥ a − b ∥ 2 ( ∥ b ∥ 2 − Δ 2 ) ∥ a − b ∥ 2 (III-3-4) \beta = \frac{-c \,{\color{red}\pm}\, \sqrt{c^2 - \|\mathbf{a}-\mathbf{b}\|^2 (\|\mathbf{b}\|^2-\Delta^2)} }{\|\mathbf{a}-\mathbf{b}\|^2} \tag{III-3-4} β=ab2c±c2ab2(b2Δ2) (III-3-4)
在 Case III 中 ∥ b ∥ 2 − Δ 2 < 0 \|\mathbf{b}\|^2-\Delta^2 < 0 b2Δ2<0, 故 c 2 − ∥ a − b ∥ 2 ( ∥ b ∥ 2 − Δ 2 ) > c 2 c^2 - \|\mathbf{a}-\mathbf{b}\|^2 (\|\mathbf{b}\|^2-\Delta^2) > c^2 c2ab2(b2Δ2)>c2, 则有 ∣ c ∣ < c 2 − ∥ a − b ∥ 2 ( ∥ b ∥ 2 − Δ 2 ) |c| < \sqrt{c^2 - \|\mathbf{a}-\mathbf{b}\|^2 (\|\mathbf{b}\|^2-\Delta^2)} c<c2ab2(b2Δ2) .

如果式 (III-3-4) 中的 “ ± \color{red}\pm ±” 符号取 “ − \color{blue}- ” 号, 此时必为负根, 需要排除.

所以式 (III-3-4) 中的 “ ± \color{red}\pm ±” 符号取 “ + \color{green}+ +” 号, 得到可行解. 即
β = − c + c 2 + ∥ a − b ∥ 2 ( Δ 2 − ∥ b ∥ 2 ) ∥ a − b ∥ 2 (III-3-5) \beta = \frac{-c \,{+}\, \sqrt{c^2 + \|\mathbf{a}-\mathbf{b}\|^2 (\Delta^2-\|\mathbf{b}\|^2)} }{\|\mathbf{a}-\mathbf{b}\|^2} \tag{III-3-5} β=ab2c+c2+ab2(Δ2b2) (III-3-5)

c < 0 c < 0 c<0 时, 分子部分是两部分正值相加, 即 − c > 0 -c> 0 c>0 c 2 + ∥ a − b ∥ 2 ( Δ 2 − ∥ b ∥ 2 ) > 0 \sqrt{c^2 + \|\mathbf{a}-\mathbf{b}\|^2 (\Delta^2-\|\mathbf{b}\|^2)} > 0 c2+ab2(Δ2b2) >0, 数值计算结果会比较稳定, 无需特殊处理.

c > 0 c > 0 c>0 时, 分子部分是一正一负相加, 即 − c < 0 -c< 0 c<0 c 2 + ∥ a − b ∥ 2 ( Δ 2 − ∥ b ∥ 2 ) > 0 \sqrt{c^2 + \|\mathbf{a}-\mathbf{b}\|^2 (\Delta^2- \|\mathbf{b}\|^2)} > 0 c2+ab2(Δ2b2) >0, 分子部分可能较接近于 0, 可能会使得数字计算的误差较大. 此时做一下分子有理化得到
β = Δ 2 − ∥ b ∥ 2 c + c 2 + ∥ a − b ∥ 2 ( Δ 2 − ∥ b ∥ 2 ) (III-3-6) \beta = \frac{\Delta^2 -\|\mathbf{b}\|^2}{c \,{+}\, \sqrt{c^2 + \|\mathbf{a}-\mathbf{b}\|^2 (\Delta^2-\|\mathbf{b}\|^2)}} \tag{III-3-6} β=c+c2+ab2(Δ2b2) Δ2b2(III-3-6)
使得分子和分母都是相对大一点的值, 以消除可能的数值计算病态问题.

最后, 插值权重参数 β \beta β 的具体计算整理为
β = { − c + c 2 + ∥ a − b ∥ 2 ( Δ 2 − ∥ b ∥ 2 ) ∥ a − b ∥ 2 , c ≤ 0 Δ 2 − ∥ b ∥ 2 c + c 2 + ∥ a − b ∥ 2 ( Δ 2 − ∥ b ∥ 2 ) , c > 0 (III-3-7) \beta = \left\{ \begin{array}{ll} \frac{-c \,{+}\, \sqrt{c^2 + \|\mathbf{a}-\mathbf{b}\|^2 (\Delta^2-\|\mathbf{b}\|^2)} }{\|\mathbf{a}-\mathbf{b}\|^2} , & c \leq 0\\ \frac{\Delta^2 -\|\mathbf{b}\|^2}{c \,{+}\, \sqrt{c^2 + \|\mathbf{a}-\mathbf{b}\|^2 (\Delta^2-\|\mathbf{b}\|^2)}}, & c> 0 \end{array} \right. \tag{III-3-7} β= ab2c+c2+ab2(Δ2b2) ,c+c2+ab2(Δ2b2) Δ2b2,c0c>0(III-3-7)


IV. 狗腿法的流程

1. 狗腿法的信赖域控制

在列文伯格-马夸尔特法中, 定义了增益比率来控制该方法中阻尼参数, 从而获得不同的迭代方法特性. 狗腿法如法炮制, 定义增益比率 (Gain Ratio),
ϱ [ i ] = g ( x [ i ] ) − g ( x [ i ] + d l h [ i ] ) L ( 0 ) − L ( d l h [ i ] ) (IV-1-1) \varrho_{[i]} = \frac{g(\mathbf{x}_{[i]})-g(\mathbf{x}_{[i]} + {^{dl}\mathbf{h}_{[i]})}}{L(\mathbf{0}) - L({^{dl}\mathbf{h}_{[i]}})} \tag{IV-1-1} ϱ[i]=L(0)L(dlh[i])g(x[i])g(x[i]+dlh[i])(IV-1-1)
通过增益比率 ϱ [ i ] \varrho_{[i]} ϱ[i] 来评估近似模型 L L L 对目标函数 g g g 的近似程度, 进而控制信赖域半径大小. 信赖域半径大小决定了狗腿法迭代步的优化性能.

如果信赖域半径太大, 近似模型和实际目标函数相差很远, 利用近似模型预测的最小值点 (迭代步到达点) 与实际目标函数最小值点相差很远, 迭代步无法以最优的方式到达更小目标值.

如果信赖域半径太小, 虽然在这个小小的信赖域内近似模型对实际目标函数拟合得很好, 但是通过狗腿法获得的是仅限于这个小小收敛域的优化解, 那么每一步的更新步长都很小, 无法实现快速迭代收敛.

增益比率 ϱ [ i ] \varrho_{[i]} ϱ[i] 的出现就是要解决上述问题.

如果增益比率 ϱ [ i ] \varrho_{[i]} ϱ[i] 为负数或者虽为正数但接近于 0 (数值比较小), 说明在当前的信赖域内近似模型 L L L 无法较好地拟合目标函数 g g g. 此时说明信赖域半径过大, 需要减小信赖域半径.

如果增益比率 ϱ [ i ] \varrho_{[i]} ϱ[i] 接近于 1 (数值比较大), 说明在当前信赖域内近似模型 L L L 非常好地拟合了目标函数 g g g. 此时为了更快地达到收敛, 可以在迭代步长上更加激进, 即可以增加信赖域半径.

得到如下信赖域半径控制算法[2]:
if ϱ > 0.75 Δ : = max ⁡ { Δ , 3 ∗ ∥ d l h ∥ } elseif ϱ < 0.25 Δ : = Δ / 2 \begin{aligned} \textbf{if}\;\; &\varrho > 0.75\\ &\Delta := \max\{\Delta, \,3 *\|^{dl}\mathbf{h}\| \} \\ \textbf{elseif} \;\; &\varrho < 0.25\\ &\Delta := \Delta/2 \\ \end{aligned} ifelseifϱ>0.75Δ:=max{Δ,3dlh}ϱ<0.25Δ:=Δ/2


2. 狗腿法的停止条件

狗腿法的算法终止条件和列文伯格-马夸尔特法中的基本一致.

条件一. 梯度不再下降

∥ ∇ g ( x [ i ] ) ∥ ∞ ≤ ε 1 (IV-2-1) \| \nabla g(\mathbf{x}_{[i]}) \|_{\infin} \leq \varepsilon_1 \tag{IV-2-1} ∥∇g(x[i])ε1(IV-2-1)

其中 ε 1 \varepsilon_1 ε1 为一小量.

条件二. 迭代点不更新

∥ x [ i + 1 ] − x [ i ] ∥ ≤ ε 2 ( ∥ x [ i ] ∥ + ε 2 ) (IV-2-2) \|\mathbf{x}_{[i+1]} - \mathbf{x}_{[i]}\| \leq \varepsilon_2(\|\mathbf{x}_{[i]}\|+\varepsilon_2) \tag{IV-2-2} x[i+1]x[i]ε2(x[i]+ε2)(IV-2-2)

其中 ε 2 \varepsilon_2 ε2 为一小量. 当 ∥ x [ i ] ∥ = 0 \|\mathbf{x}_{[i]}\|=0 x[i]=0 时, 上式右手边为 ε 2 2 \varepsilon_2^2 ε22.

条件三. 残差足够小

∥ r ( x [ i ] ) ∥ ∞ ≤ ε 3 (IV-2-3) \|\mathbf{r}(\mathbf{x}_{[i]})\|_{\infty} \leq \varepsilon_3 \tag{IV-2-3} r(x[i])ε3(IV-2-3)

其中 ε 3 \varepsilon_3 ε3 为一小量. 其实条件三能够被条件一涵盖, 参见式 (II-1-3).

条件四. 达到最大迭代数

i ≥ i max ⁡ (IV-2-4) i \geq i_{\max} \tag{IV-2-4} iimax(IV-2-4)

其中 i max ⁡ i_{\max} imax 为最大迭代步数, 防止程序无限循环.

以上四个停止条件任意一个得到满足, 算法程序就终止运行并返回运算结果.


3. 狗腿法的实现流程

狗腿法的算法流程[2]如下:

Powell’s Dog Leg Method begin i : = 0 x : = x 0 Δ : = Δ 0 ∇ g ( x ) : = J r ( x ) T r ( x ) f o u n d : = ( ∥ r ( x ) ∥ ∞ ≤ ε 3 ) or ( ∥ ∇ g ( x ) ∥ ∞ ≤ ε 1 ) while ( not f o u n d ) and ( i < i max ⁡ ) i : = i + 1 s d h : = − ∇ g ( x ) α : = ∥ s d h ∥ 2 ∥ J r ( x ) s d h ∥ 2 c p h : = α s d h {Steepest descent step} g n h : = − [ J r ( x ) T J r ( x ) ] − 1 ∇ g ( x ) {Gauss-Newton step} ϱ : = − 1 while ( ϱ < 0 ) and ( not f o u n d ) {Until step acceptable} if ∥ g n h ∥ ≤ Δ d l h : = g n h elseif ∥ c p h ∥ ≥ Δ d l h : = Δ [ i ] ∥ s d h [ i ] ∥ s d h [ i ] else c : = c p h T ( g n h − c p h ) if c ≤ 0 β : = − c + c 2 + ∥ g n h − c p h ∥ 2 ( Δ 2 − ∥ c p h ∥ 2 ) ∥ g n h − c p h ∥ 2 else β : = Δ 2 − ∥ c p h ∥ 2 c + c 2 + ∥ g n h − c p h ∥ 2 ( Δ 2 − ∥ c p h ∥ 2 ) d l h : = c p h + β ( g n h − c p h ) if ∥ d l h ∥ ≤ ε 2 ( ∥ x ∥ + ε 2 ) f o u n d : = true else x new : = x + d l h ϱ : = g ( x ) − g ( x new ) L ( 0 ) − L ( d l h ) if ϱ > 0 {Step acceptable} x : = x new ∇ g ( x ) : = J r ( x ) T r ( x ) f o u n d : = ( ∥ r ( x ) ∥ ∞ ≤ ε 3 ) or ( ∥ ∇ g ( x ) ∥ ∞ ≤ ε 1 ) if ϱ > 0.75 {Expanding trust region} Δ : = max ⁡ { Δ , 3 ∗ ∥ d l h ∥ } elseif ϱ < 0.25 {Shrinking trust region} Δ : = Δ / 2 f o u n d : = ( Δ ≤ ε 2 ( ∥ x ∥ + ε 2 ) ) end \begin{array}{ll} \textbf{Powell's Dog Leg Method}\\ \textbf{begin}\\ \qquad i:=0\\ \qquad \mathbf{x}:=\mathbf{x}_0\\ \qquad \Delta:=\Delta_0\\ \qquad \nabla {g}(\mathbf{x}) := \mathbf{J}_r(\mathbf{x}) ^{\rm\small T}\, \mathbf{r}(\mathbf{x})\\ \qquad found:= (\|\mathbf{r}(\mathbf{x})\|_{\infty} \leq \varepsilon_3) \; \textbf{or}\; \left(\| \nabla {g}(\mathbf{x}) \|_{\infin} \leq \varepsilon_1\right)\\ \textbf{while} \;\;(\textbf{not} \;found) \; \textbf{and}\; (i<i_{\max})\\ \qquad i:= i+1\\ \qquad {^{sd}\mathbf{h}} := - \nabla {g}(\mathbf{x}) \\ \qquad \alpha :=\frac{\left\| {^{sd}\mathbf{h}} \right\|^2}{\left\|{\mathbf{J}_r(\mathbf{x}) \,{^{sd}\mathbf{h}}}\right\|^2}\\ \qquad {^{cp}\mathbf{h}} := {\alpha\, {^{sd}\mathbf{h}}} \;\quad\qquad\qquad\qquad\qquad\qquad\qquad\quad \text{\color{grey}\{Steepest descent step\}}\\ \qquad {^{gn}\mathbf{h}} := - \left[ \mathbf{J}_r(\mathbf{x})^{\small\rm T} \mathbf{J}_r(\mathbf{x}) \right]^{-1} \, \nabla g(\mathbf{x}) \qquad\qquad\quad \text{\color{grey}\{Gauss-Newton step\}}\\ \qquad \varrho := -1\\ \qquad \textbf{while} \;\; (\varrho<0)\;\textbf{and}\; (\textbf{not}\; found)\,\quad \qquad\qquad\; \text{\color{grey}\{Until step acceptable\}}\\ \qquad\qquad \textbf{if}\;\; \|{^{gn}\mathbf{h}}\| \leq \Delta\\ \qquad\qquad\qquad {^{dl}\mathbf{h}}:= {^{gn}\mathbf{h}}\\ \qquad\qquad \textbf{elseif}\;\; \| {^{cp}\mathbf{h}} \| \geq \Delta\\ \qquad\qquad\qquad {^{dl}\mathbf{h}}:= {\frac{\Delta_{[i]}}{\|^{sd}\mathbf{h}_{[i]}\|} {^{sd}\mathbf{h}_{[i]}}}\\ \qquad\qquad \textbf{else}\;\; \\ \qquad\qquad\qquad c:= {^{cp}\mathbf{h}}^{\small\rm T} \left({^{gn}\mathbf{h}} - {^{cp}\mathbf{h}}\right)\\ \qquad\qquad\qquad \textbf{if} \;\; c \leq 0\\ \qquad\qquad\qquad\qquad \beta:=\frac{-c \,{+}\, \sqrt{c^2 + \|{^{gn}\mathbf{h}}-{^{cp}\mathbf{h}}\|^2 (\Delta^2-\|{^{cp}\mathbf{h}}\|^2)} }{\|{^{gn}\mathbf{h}} -{^{cp}\mathbf{h}}\|^2} \\ \qquad\qquad\qquad \textbf{else}\\ \qquad\qquad\qquad\qquad \beta:=\frac{\Delta^2 -\|{^{cp}\mathbf{h}}\|^2}{c \,{+}\, \sqrt{c^2 + \|{^{gn}\mathbf{h}}-{^{cp}\mathbf{h}}\|^2 (\Delta^2-\|{^{cp}\mathbf{h}}\|^2)}}\\ \qquad\qquad\qquad {^{dl}\mathbf{h}}:= {^{cp}\mathbf{h}} + \beta ( ^{gn}\mathbf{h} - {^{cp}\mathbf{h}} )\\ \qquad\qquad \textbf{if} \;\; \|{^{dl}\mathbf{h}}\| \leq \varepsilon_{2} (\|\mathbf{x}\|+\varepsilon_2)\\ \qquad\qquad\qquad found:= \textbf{true}\\ \qquad\qquad \textbf{else} \\ \qquad\qquad\qquad \mathbf{x}_{\text{new}} := \mathbf{x} + {^{dl}\mathbf{h}}\\ \qquad\qquad\qquad \varrho := \frac{g(\mathbf{x})-g(\mathbf{x}_{\text{new}})}{L(\mathbf{0}) - L({^{dl}\mathbf{h}})}\\ \qquad\qquad\qquad \textbf{if} \;\; \varrho>0 \qquad\qquad \qquad\qquad \qquad\qquad \text{\color{grey}\{Step acceptable\}}\\ \qquad\qquad\qquad\qquad \mathbf{x} := \mathbf{x}_{\text{new}} \\ \qquad\qquad\qquad\qquad \nabla {g}(\mathbf{x}) := \mathbf{J}_r(\mathbf{x}) ^{\rm\small T}\, \mathbf{r}(\mathbf{x})\\ \qquad\qquad\qquad \qquad found:= (\|\mathbf{r}(\mathbf{x})\|_{\infty} \leq \varepsilon_3) \; \textbf{or}\; \left(\| \nabla {g}(\mathbf{x}) \|_{\infin} \leq \varepsilon_1\right)\\ \qquad\qquad\qquad \textbf{if}\;\; \varrho > 0.75 \;\; \qquad\qquad\qquad \qquad\qquad \text{\color{grey}\{Expanding trust region\}} \\ \qquad\qquad\qquad\qquad \Delta := \max\{\Delta, \,3 *\|^{dl}\mathbf{h}\| \}\\ \qquad\qquad\qquad \textbf{elseif}\;\; \varrho < 0.25 \;\;\qquad\qquad \qquad\qquad \text{\color{grey}\{Shrinking trust region\}} \\ \qquad\qquad\qquad\qquad \Delta := \Delta/2 \\ \qquad\qquad\qquad\qquad found:=(\Delta \leq \varepsilon_{2} (\|\mathbf{x}\|+\varepsilon_2))\\ \textbf{end} \end{array} Powell’s Dog Leg Methodbegini:=0x:=x0Δ:=Δ0g(x):=Jr(x)Tr(x)found:=(r(x)ε3)or(∥∇g(x)ε1)while(notfound)and(i<imax)i:=i+1sdh:=g(x)α:=Jr(x)sdh2sdh2cph:=αsdh{Steepest descent step}gnh:=[Jr(x)TJr(x)]1g(x){Gauss-Newton step}ϱ:=1while(ϱ<0)and(notfound){Until step acceptable}ifgnhΔdlh:=gnhelseifcphΔdlh:=sdh[i]Δ[i]sdh[i]elsec:=cphT(gnhcph)ifc0β:=gnhcph2c+c2+gnhcph2(Δ2cph2) elseβ:=c+c2+gnhcph2(Δ2cph2) Δ2cph2dlh:=cph+β(gnhcph)ifdlhε2(x+ε2)found:=trueelsexnew:=x+dlhϱ:=L(0)L(dlh)g(x)g(xnew)ifϱ>0{Step acceptable}x:=xnewg(x):=Jr(x)Tr(x)found:=(r(x)ε3)or(∥∇g(x)ε1)ifϱ>0.75{Expanding trust region}Δ:=max{Δ,3dlh}elseifϱ<0.25{Shrinking trust region}Δ:=Δ/2found:=(Δε2(x+ε2))end


IV. 总结

本篇博文中, 我们先分析了列文伯格-马夸尔特法的特点, 引出我们要分析的狗腿法.

然后介绍了狗腿法涉及的两个基础算法类别, 线搜索类型和信赖域类型.

接着结合线搜索和信赖域介绍了柯西点.

有了以上的理论基础作铺垫, 我们开始介绍狗腿法的基本原理, 包括该方法的构建、优化性质的说明、插值权重参数的计算.

最后介绍狗腿法的算法流程, 分别是信赖域的控制、算法停止条件, 以及完整的狗腿算法流程.


在完整的狗腿算法流程中, 我们可以发现

- 每一迭代步只要完成一次高斯-牛顿步这类大计算工作即可;

- 即使遇到迭代步不获认可的情况, 也只需更新信赖域并重新组合迭代步, 而不需要重复计算如高斯-牛顿步等大工作.

在算法计算量/效率上明显区别于列文伯格-马夸尔特法的地方.

在列文伯格-马夸尔特法中, 如果遇到迭代步不获认可, 需要针对修改后的阻尼参数重新计算整个阻尼高斯-牛顿步.

所以说狗腿法的计算效率更高, 节省了大量计算消耗.


针对狗腿法的变化与改进也有很多, 比如狗腿法的二维子空间算法、双重狗腿法等, 本篇博文不再展开.


(如有问题, 请指出, 谢谢!)


参考文献

[1] 刘浩洋, 户将, 李勇锋, 文再文, “最优化: 建模、算法与理论”, 高教出版社, 2020

[2] K. Madsen, H.B. Nielsen, O. Tingleff, METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS, 2nd Edition, Informatics and Mathematical Modelling Technical University of Denmark, http://www2.imm.dtu.dk/pubdb/edoc/imm3215.pdf, 2004

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/635628.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Java 全栈知识点问题汇总(上)

Java 全栈知识点问题汇总&#xff08;上&#xff09; 1 Java 基础 1.1 语法基础 面向对象特性&#xff1f;a a b 与 a b 的区别3*0.1 0.3 将会返回什么? true 还是 false?能在 Switch 中使用 String 吗?对equals()和hashCode()的理解?final、finalize 和 finally 的不同…

Git 配置与理解

简述 Git 在 Windows 和 Ubuntu 中的配置&#xff0c;以及对 Git 工作区域划分和 Git 中对于文件状态划分的理解。 git 基础安装与配置 基于 WSL 的 Ubuntu 下的 git 打开或关闭Windows功能 -> Hyper-V、Virtual Machine Platform、Windows Subsystem for Linux # 1.必须…

STM32407用汇顶的GT911触摸芯片调试实盘

这个配置很关键 代码 #include "stm32f4xx.h" #include "GT9147.h" #include "Touch.h" #include "C_Touch_I2C.h" #include "usart.h" #include "delay.h" #include "LCD.h" #incl…

基于SSM的图书馆管理系统(有报告)。Javaee项目。ssm项目。

演示视频&#xff1a; 基于SSM的图书馆管理系统&#xff08;有报告&#xff09;。Javaee项目。ssm项目。 项目介绍&#xff1a; 采用M&#xff08;model&#xff09;V&#xff08;view&#xff09;C&#xff08;controller&#xff09;三层体系结构&#xff0c;通过Spring Sp…

win系统环境搭建(十四)——Windows系统下使用docker安装mysql8和mysql5.7

windows环境搭建专栏&#x1f517;点击跳转 win系统环境搭建&#xff08;十四&#xff09;——Windows系统下使用docker安装mysql8和mysql5.7 文章目录 win系统环境搭建&#xff08;十四&#xff09;——Windows系统下使用docker安装mysql8和mysql5.7MySQL81.新建文件夹2.创建…

设计模式篇章(4)——十一种行为型模式

这个设计模式主要思考的是如何分配对象的职责和将对象之间相互协作完成单个对象无法完成的任务&#xff0c;这个与结构型模式有点像&#xff0c;结构型可以理解为静态的组合&#xff0c;例如将不同的组件拼起来成为一个更大的组件&#xff1b;而行为型更是一种动态或者具有某个…

three.js从入门到精通系列教程016 - three.js通过OrbitControls对立方体实现旋转和缩放

<!DOCTYPE html> <html><head><meta charset"UTF-8"><title>three.js从入门到精通系列教程016 - three.js通过OrbitControls对立方体实现旋转和缩放</title><script src"ThreeJS/three.js"></script><…

EasyRecovery2024免费电脑数据恢复软件下载

easyrecovery是一款功能强大、易于使用的硬盘数据恢复软件。这款软件可以帮助用户非常方便地恢复丢失的数据。软件非常容易使用和高效的数据恢复。感兴趣的朋友们赶快来下载吧。 无论是因为意外删除、格式化、病毒感染、系统崩溃还是其他原因&#xff0c;该软件可以帮助您恢复…

你还在找PDF转Word的工具?一个好软件推荐,赶紧查收!

前言 前段时间朋友跟小白吐槽&#xff1a;为啥PDF文件转Word文档的工具都要收费&#xff1f; 因为他们都收费啊。 小白之前找了很多类似有这种功能的工具&#xff0c;都发现&#xff1a;不但收费&#xff0c;可能还附带全家桶&#xff0c;而且还有……广告&#xff01; 在一次…

Pytest插件“pytest-selenium” - 让自动化测试更简洁

在现代Web应用的开发中,自动化测试成为确保网站质量的重要手段之一。而Pytest插件 pytest-selenium 则为开发者提供了简单而强大的工具,以便于使用Python进行Web应用的自动化测试。本文将深入介绍 pytest-selenium 插件的基本用法和实际案例,助你轻松进入无忧的Web应用测试之…

中文词向量训练-案例分析

1 数据预处理&#xff0c;解析XML文件并分词 #!/usr/bin/env python # -*- coding: utf-8 -*- # process_wiki_data.py 用于解析XML&#xff0c;将XML的wiki数据转换为text格式 import logging import os.path import sys from gensim.corpora import WikiCorpus import jieba…

phpStorm 设置终端为git bash

环境&#xff1a; windows , PhpStorm 2022 为自己的终端配置git样式的使用&#xff0c; 默认终端样式 一、打开设置&#xff0c;选择git bin 二、重新打开终端 不加--login -i 的终端 加了--login -i 的终端 最重要的一点是什么&#xff0c;他可以像mac一样支持 ctrlv 复…

扎克伯格宣布将购买35万个GPU

Meta公司马克.扎克伯格1月18日在Instagram上发表文章称&#xff0c;该公司正在加强人工智能研究团队的力量&#xff0c;并在充实AI基础设施“弹药库“&#xff0c;计划在今年年底前向芯片设计商英伟达购买35万个H100 GPU芯片&#xff0c;从而使该公司的GPU总量达到约60万个&…

利用预训练模型SKEP进行情感分析

项目地址&#xff1a;文本情感分析 - 飞桨AI Studio星河社区 (baidu.com) baidu/Senta: Baidus open-source Sentiment Analysis System. (github.com) 本项目将详细全面介绍情感分析任务的两种子任务&#xff0c;句子级情感分析和目标级情感分析。 同时演示如何使用情感分析…

线性规划案例分享

今天想写一个最优传输的简单实现&#xff0c;结果学歪了&#xff0c;学到线性规划去了&#xff0c;这里我发现了一个宝藏网站 虽然是讲计量经济的&#xff0c;但是里面提供的公式和代码我很喜欢&#xff0c;有时间可以好好读一下 https://python.quantecon.org/lp_intro.html …

如何一键部署本地Java项目到服务器上

一、背景 我开发了一个Java代码&#xff0c;现在想部署到服务器上&#xff0c;当然可以使用Jenkins部署&#xff0c;但是Jenkins配置和维护成本比较高&#xff0c;所以我今天分享的是轻量级的一键部署脚本 演示&#xff1a;本地Window的Java代码 -> Vmware虚拟机Centos7上…

面试题:RabbitMQ 有哪几种消息模式?

文章目录 前言核心组成Rabbitmq 消息模式3.1 Simple 模式ProductorCustomer 3.2 Fanout 模式ProductorCustomer 3.3 Direct 模式Productor 3.4 Topic 模式Productor 3.5 Work 模式3.5.1 轮询分发ProductorWorker1 3.5.2 公平分发Worker1 防止消息丢失机制4.1 消息确认4.2 持久化…

在WIN从零开始在QMUE上添加一块自己的开发板(一)

文章目录 一、前言二、源码编译&#xff08;一&#xff09;安装Msys2&#xff08;二&#xff09;配置GCC工具链&#xff08;三&#xff09;安装QEMU构建依赖&#xff08;四&#xff09;下载编译QEMU源码 二、QUME编程基础&#xff08;一&#xff09;QOM机制&#xff08;二&…

LabVIEW振动筛螺栓松动故障诊断

LabVIEW振动筛螺栓松动故障诊断 概述&#xff1a;利用LabVIEW解决振动筛螺栓松动的故障诊断问题。通过集成的方法&#xff0c;不仅提高了故障检测的准确性&#xff0c;还优化了维护流程&#xff0c;为类似的机械设备故障提供了可靠的解决方案。 由于工作条件复杂&#xff0c;…

Linux系统安装NFS服务器

NFS是一种网络文件系统&#xff0c;英文全称Network File System&#xff0c;通过NFS可以让不同的主机系统之间共享文件或目录。通过NFS&#xff0c;用户可以直接在本地NFS客户端读写NFS服务端上的文件&#xff0c;是非常好的共享存储工具。本篇文章将介绍如何在CentOS7上安装N…