METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS（六）

之后的两个方法都是针对无法进行解析求导情况下的处理方法，大多数情况下并不会用到

3.5. L-M 方法的割线版本

本手册中讨论的方法假设向量函数 f\pmb{f}fff 是可微的，即雅各比矩阵
J(x)=[∂fi∂xj]\pmb{J}(\pmb{x}) = \begin{bmatrix} \frac{\partial \pmb{f}_i}{\partial \pmb{x}_j} \end{bmatrix} JJJ(xxx)=[∂xxxj∂fffi]

存在。在许多实际优化问题中，我们无法为 J\pmb{J}JJJ 中的元素给出公式，例如因为 f\pmb{f}fff 是由“黑匣子”给出的。 L-M 方法的割线版本旨在解决此类问题。

最简单的补救方法是将 J(x)\pmb{J}(\pmb{x})JJJ(xxx) 替换为通过数值微分获得的矩阵 B\pmb{B}BBB ：第 (i,j)(i, j)(i,j) 个元素由有限差分逼近近似
∂fi∂xj(x)≈fi(x+δej)−fi(x)δ=bij(3.28)\frac{\partial f_i}{\partial x_j}(\pmb{x}) \approx \frac{f_i(\pmb{x}+\delta \pmb{e}_j) - f_i(\pmb{x})}{\delta} = b_{ij} \tag{3.28} ∂xj∂fi(xxx)≈δfi(xxx+δeeej)−fi(xxx)=bij(3.28)

其中 ej\pmb{e}_jeeej 是第 jjj 个坐标方向上的单位向量，而 δ\deltaδ 是一个适当小的实数。使用这种策略，每次迭代 x\pmb{x}xxx 都需要对 f\pmb{f}fff 进行 n+1n+1n+1 次评估，并且由于 δ\deltaδ 可能远小于距离 ∣∣x−x∗∣∣||\pmb{x}-\pmb{x}^*||∣∣xxx−xxx∗∣∣，因此我们无法从仅对 f(x)\pmb{f}(\pmb{x})fff(xxx) 的评估中获得更多关于 f\pmb{f}fff 全局行为的信息。我们想要更高的效率。

示例 3.14.

令 m=n=1m = n = 1m=n=1 并考虑一个非线性方程
f:R→Rfindx^suchthatf(x^)=0f:\mathbb{R} \to \mathbb{R} \quad find \, \hat{x} \, such \, that \, f(\hat{x}) = 0 f:R→Rfindx^suchthatf(x^)=0

对于这个问题，我们可以将牛顿-拉夫森算法（3.6）写成如下形式
f(x+h)≈l(h)=f(x)+f˙(x)hsolvethelinearprobleml(h)=0xnew:=x+h(3.29)f(x+h) \approx l(h) = f(x) + \dot{f}(x)h \quad \\ solve \, the \, linear \, problem \, l(h) = 0 \\ x_{new}:= x +h \quad \quad \quad \quad \quad \quad \quad \quad \tag{3.29} f(x+h)≈l(h)=f(x)+f˙(x)hsolvethelinearprobleml(h)=0xnew:=x+h(3.29)

如果我们不能实现 f˙(x)\dot{f}(x)f˙(x)，那么我们可以将其近似为
(f(x+δ)−f(x))/δ(f(x+\delta) - f(x))/ \delta (f(x+δ)−f(x))/δ

δ\deltaδ需要选择的适当小。更一般地，我们可以将（3.29）替换为
f(x+h)≈λ(h)=f(x)+bhwithb≈f˙(x)solvethelinearprobleml(h)=0xnew:=x+h(3.30 a)f(x+h) \approx \lambda(h) = f(x) + bh \quad with \quad b \approx \dot{f}(x) \\ solve \, the \, linear \, problem \, l(h) = 0 \quad \quad \quad \quad \quad \\ x_{new}:= x +h \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \tag{3.30 a} f(x+h)≈λ(h)=f(x)+bhwithb≈f˙(x)solvethelinearprobleml(h)=0xnew:=x+h(3.30 a)

假设我们已经知道 xprevx_{prev}xprev 和 f(xprev)f(x_{prev})f(xprev)。然后我们可以通过构造下式来求解 bbb
f(xprev)=λ(xprev−x)(3.30 b)f(x_{prev}) = \lambda(x_{prev} - x) \tag{3.30 b} f(xprev)=λ(xprev−x)(3.30 b)

注意，这里h=xprev−xh = x_{prev} - xh=xprev−x

这个式子给出 b=(f(x)−f(xprev)/(x−xprev))b = (f(x) - f(x_{prev}) / (x-x_{prev}))b=(f(x)−f(xprev)/(x−xprev))，并且通过选择 bbb 我们将（3.30）看作为割线方法，例如 Elden et al (2004)中的77页。

作为牛顿-拉夫森方法的有限差分近似的替代，割线方法的主要优点是我们每个迭代步骤只需要一个函数评估，而不是两个。

现在，考虑 f:Rn→Rmf:\mathbb{R}^n \to \mathbb{R}^mf:Rn→Rm 的线性模型（3.7a），
f(x+h)≈l(h)=f(x)+J(x)h\pmb{f}(\pmb{x}+\pmb{h}) \approx \pmb{l}(\pmb{h}) = \pmb{f}(\pmb{x})+\pmb{J}(\pmb{x})\pmb{h} fff(xxx+hhh)≈lll(hhh)=fff(xxx)+JJJ(xxx)hhh

我们将其替换为
f(x+h)≈λ(h)=f(x)+Bh\pmb{f}(\pmb{x}+\pmb{h}) \approx \lambda(\pmb{h}) = \pmb{f}(\pmb{x}) + \pmb{B}\pmb{h} fff(xxx+hhh)≈λ(hhh)=fff(xxx)+BBBhhh

其中 B\pmb{B}BBB 是 J(x)\pmb{J}(\pmb{x})JJJ(xxx) 的当前近似值。在下一个迭代步骤中，我们需要 Bnew\pmb{B}_{new}BBBnew 满足

f(xnew+h)≈f(xnew)+Bnewh\pmb{f}(\pmb{x}_{new} + \pmb{h}) \approx \pmb{f}(\pmb{x}_{new}) + \pmb{B}_{new}\pmb{h} fff(xxxnew+hhh)≈fff(xxxnew)+BBBnewhhh

特别是，我们希望这个模型对 h=x−xnew\pmb{h} = \pmb{x}−\pmb{x}_{new}hhh=xxx−xxxnew 成立，即
f(x)=f(xnew)+Bnew(x−xnew)(3.31 a)\pmb{f}(\pmb{x}) = \pmb{f}(\pmb{x}_{new}) + \pmb{B}_{new}(\pmb{x} - \pmb{x}_{new}) \tag{3.31 a} fff(xxx)=fff(xxxnew)+BBBnew(xxx−xxxnew)(3.31 a)

这给了我们关于 Bnew\pmb{B}_{new}BBBnew 中 m⋅nm\cdot nm⋅n 个未知元素中的 mmm 个方程，所以我们需要更多的条件。 Broyden (1965) 建议使用下列方程补充（3.31 a）
Bnewv=Bvforallv⊥(x−xnew)(3.31 b)\pmb{B}_{new}\pmb{v} = \pmb{B}\pmb{v} \quad for \, all \, \pmb{v} \perp (\pmb{x} - \pmb{x}_{new}) \tag{3.31 b} BBBnewvvv=BBBvvvforallvvv⊥(xxx−xxxnew)(3.31 b)

很容易验证条件（3.31a-b）在使用以下更新方程的情况下被满足

定义 3.32. 布罗伊登秩一更新

Bnew=B+uhT(3.32)\pmb{B}_{new} = \pmb{B} + \pmb{u}\pmb{h}^T \tag{3.32} BBBnew=BBB+uuuhhhT(3.32)
其中
h=xnew−x,u=1hTh(f(xnew)−f(x)−Bh)\pmb{h} = \pmb{x}_{new} - \pmb{x}, \quad \pmb{u} = \frac{1}{\pmb{h}^T\pmb{h}}(\pmb{f}(\pmb{x}_{new}) - \pmb{f}(\pmb{x}) - \pmb{B}\pmb{h}) hhh=xxxnew−xxx,uuu=hhhThhh1(fff(xxxnew)−fff(xxx)−BBBhhh)

请注意，条件（3.31a）对应于 n=1n = 1n=1 的情况下的割线条件（3.30b）。我们说这种方法是一种广义割线法。

带有此修改的算法 3.16 的中心部分的简要概述具有以下形式
solve(BTB+μI)hslm=−BTf(x)xnew:=x+hslmUpdateBby(3.32)UpdateμandxasinAlgorithm3.16solve \quad (\pmb{B}^T\pmb{B} + \mu \pmb{I}) \pmb{h}_{slm} = -\pmb{B}^T\pmb{f}(x) \\ \pmb{x}_{new}:= \pmb{x}+\pmb{h}_{slm} \quad \quad \quad \quad \quad \quad \quad \\ Update \, \pmb{B} \, by \, (3.32) \quad \quad \quad \quad \quad \quad \\ \quad Update \, \mu \, and \, \pmb{x} \, as \, in \, Algorithm \, 3.16 solve(BBBTBBB+μIII)hhhslm=−BBBTfff(x)xxxnew:=xxx+hhhslmUpdateBBBby(3.32)UpdateμandxxxasinAlgorithm3.16

Powell 已经证明如果向量集合 x0,x1,x2,...\pmb{x}_0, \pmb{x}_1, \pmb{x}_2,...xxx0,xxx1,xxx2,... 收敛到 x∗\pmb{x}^∗xxx∗ 并且如果步长集合 {hk=xk−xk−1}\{\pmb{h}_k = \pmb{x}_k−\pmb{x}_{k−1}\}{hhhk=xxxk−xxxk−1} 满足 {hk−n+1,...,hk}\{\pmb{h}_{k−n+1},.. ., \pmb{h}_k\}{hhhk−n+1,...,hhhk} 对于每个 k≥nk \geq nk≥n 是线性独立的（它们跨越整个 Rn\mathbb{R}^nRn），那么无论 B0\pmb{B}_0BBB0 的选择如何，逼近集合 {Bk}\{\pmb{B}_k\}{BBBk} 收敛到 J(x∗)\pmb{J}(\pmb{x}^∗)JJJ(xxx∗)。

然而，在实践中，经常发生前 nnn 个步骤不跨越整个 Rn\mathbb{R}^nRn，并且存在这样的风险，即在一些迭代步骤之后，当前 B\pmb{B}BBB 对真正的雅各比矩阵的逼近非常差，以至 −BTf(x)-\pmb{B}^T\pmb{f} (\pmb{x})−BBBTfff(xxx) 甚至不是下坡方向。在这种情况下 x\pmb{x}xxx 将保持不变，而 μ\muμ 会增加。近似值 B\pmb{B}BBB 发生了变化，但仍然可能是一个很差的近似值，导致 μ\muμ 进一步增加等。最终，由于 hslm\pmb{h}_{slm}hhhslm 太小以至于满足（3.15b）导致算法停止，尽管此时 x\pmb{x}xxx 可能远离 x∗\pmb{x}^∗xxx∗。

已经提出了许多策略来克服这个问题，例如，偶尔通过有限差分重新计算 B\pmb{B}BBB。在下面的算法 3.34 中，我们用循环的、坐标系的更新来补充由迭代过程确定的更新：设 h\pmb{h}hhh 表示当前步长，并且设 jjj 为当前坐标序号。如果 h\pmb{h}hhh 和 ej\pmb{e}_jeeej 之间的角度 θ\thetaθ 很大，那么我们计算 J\pmb{J}JJJ 的第 jjj 列的有限差分近似值。更具体地说，在满足以下的条件后我们这样做
cosθ=∣hTej∣∣∣h∣∣∣∣ej∣∣<γ⇔∣hj∣<γ∣∣h∣∣(3.33)cos \theta = \frac{|\pmb{h}^T\pmb{e}_j|}{||\pmb{h}||||\pmb{e}_j||} < \gamma \Leftrightarrow |\pmb{h}_j| < \gamma ||\pmb{h}|| \tag{3.33} cosθ=∣∣hhh∣∣∣∣eeej∣∣∣hhhTeeej∣<γ⇔∣hhhj∣<γ∣∣hhh∣∣(3.33)

实验表明（相当悲观的）选择 γ=0.8\gamma = 0.8γ=0.8 提供了良好的性能。通过这种选择，我们可以预期每个迭代步骤（几乎）需要对向量函数 f\pmb{f}fff 进行两次评估。

现在我们准备介绍算法。阻尼参数 μ\muμ 的监控与算法 3.16 中的一样，为了清楚起见，我们在演示文稿中省略了它。

我们有以下说明：

初始化。 x0\pmb{x}_0xxx0 是输入，B0\pmb{B}_0BBB0 要么是输入，要么由（3.28）计算得出。停止标准（3.15）中的参数和（3.28）中使用的步长 δ\deltaδ 也是输入值。
比较（3.33）。mod(j,n)mod(j, n)mod(j,n) 是除以 nnn 后的余数。
系数η\etaη由下式给出
ifxj=0thenη:=δ2elseη=δ∣xj∣if \quad x_j = 0 \quad then \quad \eta:=\delta^2 \quad else \quad \eta = \delta |x_j| ifxj=0thenη:=δ2elseη=δ∣xj∣
而迭代值 x\pmb{x}xxx 仅在满足下降条件（2.1）时更新，而近似值 B\pmb{B}BBB 在每一步中更新。因此，当 f(x)\pmb{f}(\pmb{x})fff(xxx) 不变时，近似梯度 g\pmb{g}ggg 也可能发生变化。

需要注意，（3.32）里已经隐含了使用有限差分近似进行更新B的公式（满足（3.31a）），所以说明 3 处的公式就是使用差分近似重新计算 B\pmb{B}BBB

示例 3.15.

我们使用算法 3.34 处理示例 3.13 中对应 λ=0\lambda = 0λ=0 的修改后的 Rosenbrock 问题。如果我们使用与该示例相同的起点和停止标准，并在差分近似（3.28）中取 δ=10−7\delta = 10^{−7}δ=10−7，我们在29 次迭代后找到了解，共涉及 f(x)\pmb{f}(\pmb{x})fff(xxx) 的 53 次评估。作为比较，“真正的”L-M 算法只需要 17 步，这意味着对 f(x)\pmb{f}(\pmb{x})fff(xxx) 和 J(x)\pmb{J}(\pmb{x})JJJ(xxx) 总共进行了 18 次评估。

我们还对示例 1.1、3.7 和 3.11 中的数据拟合问题使用了割线算法。在 δ=10−7\delta = 10^{-7}δ=10−7 以及与示例 3.7 中相同的起点和停止标准的情况下，迭代在 94 步后由（3.15a）停止，总共涉及 192 次 f(x)\pmb{f}(\pmb{x})fff(xxx) 的评估。相比之下，算法 3.16 需要 62 个迭代步骤。

这两个问题表明算法 3.34 是鲁棒的，但它们也说明了一般的经验法则：如果梯度信息可用，通常使用它是值得的。

在许多应用程序中，数字 mmm 和 nnn 很大，但每个函数 fi(x)f_i(\pmb{x})fi(xxx) 仅取决于 x\pmb{x}xxx 中的少数元素。这种情况下大部分的 ∂fi∂xj\frac{\partial f_i}{\partial x_j}∂xj∂fi 为零，我们说 J(x)\pmb{J}(\pmb{x})JJJ(xxx) 是一个稀疏矩阵。在列文伯格-马尔夸特方程（3.13）的解中，有一些利用稀疏性的有效方法，例如 Nielsen (1997)。然而，在更新公式 (3.32) 中，通常向量 h\pmb{h}hhh 和 u\pmb{u}uuu 中的所有元素都是非零的，因此 Bnew\pmb{B}_{new}BBBnew 将是一个稠密矩阵。讨论如何处理这个问题超出了本手册的范围；我们参考 Gill et al (1984) 和 Toint (1987)。

3.6. Dog Leg 方法的割线版本

当然，对雅各比矩阵使用割线近似的想法也可以与第 3.3 节中的 Dog Leg 方法结合使用。在本节中，我们将考虑 m=nm = nm=n 的特殊情况，即求解非线性方程组。 Broyden (1965) 不仅给出了定义 3.32 中的公式，
B=B+(1hhT(y−Bh))hTwhereh=xnew−x,y=f(xnew−f(x))(3.35 a)\pmb{B} = \pmb{B} + (\frac{1}{\pmb{h}\pmb{h}^T}(\pmb{y}-\pmb{B}\pmb{h}))\pmb{h}^T \\ where \quad \pmb{h} = \pmb{x}_{new} - \pmb{x} , \quad \pmb{y} = \pmb{f}(\pmb{x}_{new} - \pmb{f}(\pmb{x})) \tag{3.35 a} BBB=BBB+(hhhhhhT1(yyy−BBBhhh))hhhTwherehhh=xxxnew−xxx,yyy=fff(xxxnew−fff(xxx))(3.35 a)

用于更新近似的雅各比矩阵。他还给出了更新雅雅各比矩阵近似逆 D≈J(x)−1\pmb{D} \approx \pmb{J}(\pmb{x})^{-1}DDD≈JJJ(xxx)−1 的公式。公式是
Dnew=D+(1hTDy(h−Dy))(hTD)(3.35 b )\pmb{D}_{new} = \pmb{D} +(\frac{1}{\pmb{h}^T\pmb{D}\pmb{y}}(\pmb{h} - \pmb{D}\pmb{y}))(\pmb{h}^T\pmb{D}) \tag{3.35 b } DDDnew=DDD+(hhhTDDDyyy1(hhh−DDDyyy))(hhhTDDD)(3.35 b )

其中 h\pmb{h}hhh 和 y\pmb{y}yyy 在（3.35a）中定义。

使用这些矩阵，最速下降方向 hsd\pmb{h}_{sd}hhhsd 和高斯-牛顿步长（在这种情况下与牛顿步长相同，参见示例 3.5）hgn\pmb{h}_{gn}hhhgn（3.18）近似为
hssd=−BTf(x)andhsgn=−Df(x)(3.36)\pmb{h}_{ssd} = -\pmb{B}^T\pmb{f}(\pmb{x}) \quad and \quad \pmb{h}_{sgn} = -\pmb{D}\pmb{f}(\pmb{x}) \tag{3.36} hhhssd=−BBBTfff(xxx)andhhhsgn=−DDDfff(xxx)(3.36)

算法 3.21 很容易修改为使用这些近似值。初始 B=B0\pmb{B} = \pmb{B}_0BBB=BBB0 可由差分逼近（3.28）求得，D0\pmb{D}_0DDD0 计算为 B0−1\pmb{B}_0^{-1}BBB0−1。很容易证明，当前的 B\pmb{B}BBB 和 D\pmb{D}DDD 满足 BD=I\pmb{B}\pmb{D} = \pmb{I}BBBDDD=III。步长参数 α\alphaα 由（3.19）找到，其中 J(x)\pmb{J}(\pmb{x})JJJ(xxx) 由 B\pmb{B}BBB 代替。

与 L-M 方法的割线版本一样，该方法需要额外的更新以保持 B\pmb{B}BBB 和 D\pmb{D}DDD 是当前雅各比矩阵及其逆的良好近似。我们发现，围绕（3.33）讨论的策略在这种情况下也很有效。还应该提到的是，（3.35b）中的分母可能为零或非常小。如果
∣hTDy∣<ϵM∣∣h∣∣|\pmb{h}^T\pmb{D}\pmb{y}| < \sqrt{\epsilon_M}||\pmb{h}|| ∣hhhTDDDyyy∣<ϵM∣∣hhh∣∣

则 D\pmb{D}DDD 不会被更新，而是由D=B−1\pmb{D} = \pmb{B}^{-1}DDD=BBB−1计算。

每次使用（3.35）花费 10n210n^210n2 flops，并且使用（3.36）计算两个步长向量加上使用（3.19）计算 α\alphaα 花费 6n26n^26n2 flops。因此，使用 Dog Leg 方法的无梯度版本的每个迭代步骤花费大约 16n216n^216n2 flops 加上 f(xnew)\pmb{f}(\pmb{x}_{new})fff(xxxnew) 的评估。为了比较，算法 3.21 的每一步花费大约23n3+6n2\frac{2}{3}n^3+6n^232n3+6n2 flops 加上f(xnew)\pmb{f}(\pmb{x}_{new})fff(xxxnew) 和 J(xnew)\pmb{J}(\pmb{x}_{new})JJJ(xxxnew) 的评估。因此，对于较大的 nnn 值，无梯度版本的每次迭代更便宜。然而，迭代步数通常要大得多，如果雅各比矩阵可用，那么梯度版本通常更快。

一个“flop”是两个浮点数之间的简单算术运算。

示例 3.16.

我们在由下式给出的 Rosenbrock 函数（参见示例 3.13）$f:\mathbb{R}^2 \to \mathbb{R}^2 $上使用了算法 3.21 和无梯度的 Dog Leg 方法。
f(x)=[10(x2−x12)1−x1]f(x) = \begin{bmatrix} 10(x_2 - x_1^2) \\ 1- x_1 \end{bmatrix} f(x)=[10(x2−x12)1−x1]

该函数有一个根，x∗=[1,1]T\pmb{x}^* = [1, 1 ]^Txxx∗=[1,1]T，在这两种方法中，我们都使用了起点 x0=[−1.2,1]T\pmb{x}_0 = [ -1.2, 1 ]^Txxx0=[−1.2,1]T 和 ϵ1=ϵ2=10−12,kmax=100\epsilon_1 = \epsilon_2 = 10^{−12}, k_{max} = 100ϵ1=ϵ2=10−12,kmax=100 在停止标准（3.15）中以及 δ=10−7\delta = 10^{−7}δ=10−7 在差分近似（3.28）中。算法 3.21 在 17 次迭代步骤后停止在解上，即在 f\pmb{f}fff 及其雅各比矩阵的 18 次评估之后。割线版本也停在解上；这需要 28 个迭代步骤和对 f\pmb{f}fff 的总共 49 次评估。

METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 翻译（六）相关推荐

METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 翻译（一）
METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS(一) 分章节更新,剩下两部分一周内上传 1. 介绍和定义在本手册中,我们考虑了以下问题定义1.1. 最小 ...
METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 翻译（三）
METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS(三) 3. 非线性最小二乘问题在本手册的其余部分中,我们将讨论求解非线性最小二乘问题的方法.给定一个向量函 ...
METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 翻译（七）
METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS(七) 3.7. 最后的话我们已经讨论了许多解决非线性最小二乘问题的算法.它们都出现在任何好的程序库中,并且 ...
数学知识--Methods for Non-Linear Least Squares Problems（第二章）
Methods for Non-Linear Least Squares Problems 非线性最小二乘问题的方法 2nd Edition, April 2004 K. Madsen, H.B. N ...
数学知识--Methods for Non-Linear Least Squares Problems（第三章）
Methods for Non-Linear Least Squares Problems 非线性最小二乘问题的方法 2nd Edition, April 2004 K. Madsen, H.B. N ...
数学知识--Methods for Non-Linear Least Squares Problems（第一章）
Methods for Non-Linear Least Squares Problems 非线性最小二乘问题的方法 2nd Edition, April 2004 K. Madsen, H.B. N ...
Python3 LinAlgError: SVD did not converge in Linear Least Squares
有用的话记得回过头请给本文点个赞,感谢您的支持! LinAlgError: SVD did not converge in Linear Least Squares 说明在拟合时,y值里存在nan值, ...
numpy.polyfit：LinAlgError: SVD did not converge in Linear Least Squares
原因:有NAN值解决方法:去除NAN值,例如去除a,b,c列中有NAN值的行 data = data.dropna(axis=0,subset=['a', 'b','c']) reference: ...
MULLS: Versatile LiDAR SLAM via Multi-metric Linear Least Square 论文翻译
目录论文题目代码: 摘要 I. 介绍 II. 相关工作 III. 方法 A.Motion compensation 运动补偿 B. Geometric feature points extract ...

METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 翻译（六）