非线性最小二乘

非线性最小二乘
tags:数值优化

---

最小二乘

请思考一下下面的问题

$$\bbox[15px, border:1px solid black] { \begin{array}{l} 定义: \; 1.1 \; 最小二乘问题 \\ 找到一个向量x^*,使得 \\ \qquad F(x) = \frac{1}{2}\sum_{i=1}^m(f_{i}(x))^2 \\ 最小，其中f_i:R^n \mapsto R, i=1, \cdots,m,m \ge n \qquad \end{array} }$$

最小二乘的一个重要用途是数据拟合，来看一个例子。假设有$m$组样本点数据$(t_1,y_1),(t_2,y_2), \cdots, (t_m, y_m)$,现在要找到一个函数使得该函数能够很好的拟合上面的数据点并能够预测新的数据点。现在设这个函数的模型

$$M(x,t) = x_3e^{x_1t} + x_4e^{x_2t}$$

其中向量$x$是该模型的参数，这里设$x=(x_1,x_2,x_3,x_4)^T$,如果存在向量$x^{+}$使得

$$y_i = M(x^{+},t_i)+\epsilon_i$$

其中$\{ \epsilon_i \}$纵坐标的测量误差，假设这个误差服从白噪声(即服从高斯分布$N(0,1)$)。对于任何一个$x$我们可以计算它的残差

$$\begin{array}{l} f_i(x) & = y_i - M(x,t_i) \\ &=y_i - x_3e^{x1t_i} - x_4e^{x_2t_i} & i=1,\cdots,m. \end{array}$$

这里使用最小二乘的定义来求解最小的$x^*$值。现实中，很多问题都是最小二乘的一个特殊的变体：给一个函数$F:R^n \mapsto R$,找满足映射关系的$F$的最小参数值。其中$F$叫做目标函数或者代价函数。

$$\bbox[15px, border:2px] { \begin{array}{} 定义:\; 1.2 \; 全局最小 & \qquad & \qquad & \qquad\\ 其中 F:R^n \mapsto R, 找 \\ \quad x^{+} = argmin_x\{F(x)\} \end{array} }$$

一般情况下，找到全局最小很困难。这里我们仅仅求解$F$的局部最小解，也就是求在$x$变换很小的情况下求解$F$在这个很小的范围$\delta$的最小解。其中$\delta$是很小的一个正数。

$$\bbox[15px, border:2px] { \begin{array}{l} 定义:1.3局部最小 & \qquad & \qquad \\ 其中F:R^n mapsto R, 找 x^*，使得 \\ \quad F(x^*) \le F(x) \quad , \Vert x-x^*\Vert \lt \delta. \end{array} }$$

我们假设目标函数是光滑可微分的，下面的泰勒的截断式是有效的

$$F(x+h) = F(x)+h^Tg+\frac{1}{2}h^THh+O(\Vert h\Vert ^3), \qquad (1.4a)$$

其中$g$是梯度函数，即：

$$g \equiv F^{'}(x) = \begin{bmatrix} \frac{\partial F(x)}{\partial x_1} \\ \vdots \\ \frac{\partial F(x)}{\partial x_n} \end{bmatrix} \qquad (1.4b)$$

$H$是黑塞矩阵，即：

$$H \equiv F^{''}(x) = \begin{matrix} \frac{\partial^2F}{ \partial x_i \partial x_j}(x) \end{matrix} \qquad (1.4c)$$

如果$x^*$是局部最小解同时$\Vert h\Vert$足够的小，那么我们不可能找的一个点$x^*+h$使得$F$的值更小。由1.4a可知：

$$\bbox[15px, border:2px] { \begin{array}{} 理论: \; 1.5 \; 局部最小化的必要条件 & \qquad & \qquad \\ 如果x^*是局部最小值，则有 \\ g^* \equiv F^{'}(x^*) = 0. \end{array} }$$

对于上面的必要条件我们给出一个具体的名字叫稳定点：

$$\bbox[15px, border:2px] { \begin{array}{} 定义: \; 1.6 \; 稳定点 & \qquad & \qquad \\ g_s \equiv F{'}(x_s) = 0 \\ 也就是说x_s是F的稳定点 \end{array} }$$

因此局部最小点也是一个稳定点，但也可能是一个局部最大点。稳定点也有可能既不是局部最小点也不是局部最大点，这个点也被称作鞍点。为了确定给定的点是不是局部最小点。我们需要把第二个条件$x_s$带入1.4a,于是有

$$F(x_s+h) = F(x_s) + \frac{1}{2}h^TH_sh+O(\Vert h \Vert ^3), \qquad H_s = F^{''}(x_s). \qquad (1.7)$$

从定义(1.4c)可知，黑塞矩阵是对称的矩阵，如果我们让$H_s$是正定的，它的特征值要比某些数字$\delta \gt 0$要大，则有

$$h^TH_sh \gt \delta \Vert h \Vert ^2$$

上式说明对于给定的足够小的$\Vert h \Vert$第二个条件起主要作用，同时第二个条件是正的，这样我们可以得到

$$\bbox[15px, border:2px] { \begin{array}{} 理论: \; 局部最小化的局部条件 & \qquad &\qquad \\ 假设x_s是一个稳定点，F^{''}(x_s)是正定的。 \\ 然后我们有x_s是一个局部最小点。 \end{array} }$$

如果$H_s$是负定的，则$x_s$是局部最大解。如果$H_s$是非定的(也就是它有正的和负的特征值)，则$x_s$是一个鞍点。

2. 下降方法

所有非线性优化方法都是迭代方法。从迭代开始点$x_0,x_1,\cdots,$，我们希望它能够收敛到被给函数的局部最小点$x^*$，看定义1.3。大多数的方法使用下面的约束强制函数是下降的

$$F(x_{k+1}) \lt F(x_k) \qquad (2.1)$$

这就防止了函数收敛到局部最大点同时也使函数更少的可能收敛到鞍点。如果被给的函数的局部最小点不止一个，那么收敛的结果值依赖于它的初始解$x_0$。因为我们不知道哪一个最小解将要被发现。所以最小解不有必要与$x_0$很接近。

在很多情况下，该方法可能要经历两个阶段才能收敛到最小解。当$x_0$远离最小解时，我们想让该方法能够稳定的朝着$x^*收敛。因此在整个的迭代过程中我们应保证误差不增加(除了在第一步时)，也就是

$$\Vert e_{k+1} \Vert \lt \Vert e_k \Vert \quad k \gt K$$

其中$e_k$表示当前的误差，

$$e_k = e_k - x^*. \qquad (2.2)$$

当$x_k$接近$x^*$时，我们想让函数更快的收敛，看看下面关于收敛速度的一些定义

$$\bbox[15px, border:2px] { \begin{array}{} 线性收敛： \qquad \qquad \\ & \Vert e_{k+1} \le \alpha \Vert e_k \Vert \\ 二次收敛： \\ & \Vert e_{k+1} \Vert = O(\Vert e_k \Vert ^2) \quad \Vert e_k \Vert 很小时 \\ 超线性收敛： \\ & \Vert e_{k+1} / \Vert e_k \Vert \quad 其中 k \to \infty \\ \end{array} }$$

下面给出下降方法的算法步骤，其中迭代的每一步我们需要考虑下面两个问题

• 找到下降的方向$h_d$
• 找到一个好的迭代步长

算法大概步骤

$$\bbox[15px, border:2px] { \begin{array}{} 算法 \, 2.4 下降方法 \\ begin \\ \quad k:= 0; x:= x_0; found := false; & \qquad \{开始点\} \\ \quad while (not found) \; and \; (k \lt k_{max}) \\ \qquad h_d:= search\_direction(x) & \qquad \{搜索x的下山方向\}\\ \qquad if \; (no such h exiists) \\ \qquad \quad found:= true; & \qquad \{x是稳定点\} \\ \qquad else \\ \qquad \quad \alpha := step_length(x, h_d) & \qquad \{x在h_d上的步长\} \\ \qquad \quad x:= x + \alpha h_d; \; k:= k+1; & \qquad\{下一次迭代\} \\ end \end{array} }$$

思考一下下面的情况，关于以$x$起点，$h$为方向的&F&的值。从泰勒阶段时我们看到

$$F(x+\alpha h) = F(x)+\alpha h^TF^{'}(x)+O(\alpha^2) \simeq F(x)+\alpha h^T F^{'}(x) \qquad 其中\alpha 充分的小 \qquad (2.5)$$

这里我们把$F(x+\alpha)$看成关于$\alpha$的函数

$$\bbox[15px, border:2px] { \begin{array}{} \phi (\alpha) = F(x+\alpha h) ,\; & \phi^{'}(\alpha) = h^TF^{'}(x+\alpha h) \\ \phi(0) = F(x) \;,& \phi^{'}(0) = h^TF^{'}(x) & \qquad (2.5) \end{array} }$$

可以看出，如果$\alpha = 0， \phi^{'}(0) = h^tF^{'}(x) \lt 0 \,$时，$F(x)$为单调减函数。所有我们可以得到下降方向$h$满足的条件

$$\bbox[15px, border:2px] { \begin{array}{} 定义: 2.6 下降方向 \\ 对于函数$F(x)$，给一个$x$,如果有h^TF^{'}(x) \lt 0 ,则h是下降方向 \end{array} }$$

如果没有那样的$h$存在，而此时又有$F^{'}(x) = 0$,则说明$x$是稳定点。否则，我们就需要选择一个步长$\alpha$,也就是选择一个沿着$h$方向下降的长度。以便于我们可以得到$F(x+\alpha h)$的值。其中一种确定$\alpha$的近似方法是找

$$\alpha_e = argmin_{a \gt 0} \{F(x+\alpha) \} \qquad (2.7)$$

这个过程叫做线性搜索。下面我们先介绍两种计算下降方向的方法。

2.1 最速下降法

从(2.5)我们可以知道完成一个$\alpha h, \alpha \gt 0$，函数的相对增益值满足

$$\lim_{\alpha \to 0} = \frac{F(x) - F(x+ \alpha h)}{\alpha \Vert h \Vert} = - \frac{1}{\Vert h \Vert} h^T F^{'}(x) = - \Vert F^{'}(x) \Vert \cos\theta$$

其中$\theta$是$h$和$F^{'}(x)$之间的夹角。所以当$\theta = \pi$时，函数增益最大，最大为$\Vert F^{'}(x)$。因此最速下降方向$h_{sd}$

$$h_{sd} = -F^{'}(x) \qquad (2.8)$$

以(2.8)为基础的方法叫做最速下降方法或者梯度方法。它是下降方向最好的方法，我们可以把它应用于(2.7)。该方法收敛的速度是线性的，比较缓慢。然而最速下降法用于精确线性搜索，在有限计算精度下有可能找不到最小点。但大多数情况下，该方法在最初的迭代过程中能够表现出很好的性能。

思考一下有没有这么一种混合方法，该方法在最初的表现的很好，如梯度方法，在接近真实解时它表现的像牛顿法。接下来我们来看看什么是牛顿法。使用混合方法的一个关键性问题是什么时候切换这两种方法最恰当。

2.2 牛顿法

我们从$x^*$是一个稳点时推导这个方法。根据(1.6)的定义，它满足$F^{'}(x) = 0$，这是一个非线性方程组，由泰勒表达式

$$F^{'}(x+h) = F^{'} + F^{''}(x)h +O(\Vert h \Vert ^2) \simeq F^{'}(x) +F^{''}(x)h \qquad \Vert h \Vert 是充分小的$$

我们推导牛顿方法：找$h_n$的解

$$Hh_n = -F^{'}(x) \qquad 其中 H = F^{''}(x) \qquad (2.9a)$$

计算下一次的迭代

$$x := x + h_n. \qquad (2.9b)$$

若假设$H$是正定非奇异的，对于所有$\mu$有$\mu^T H \mu \gt 0$，因此在(2.9a)的两边都乘以$h^T_{n}$，我们可以得到

$$0 \lt h^T_nHh_n = -h^T_n F^{'}(x) \qquad (2.10)$$

这说明了$h_n$是下降方向，它满足2.6定义的条件。

牛顿方法在迭代的最后阶段(即接近最小解)能够表现的非常的好。它也说明了如果黑塞矩阵是正定的(1.8定义的)，在真实解附近$F^{''}(x)$也是正定的，然后就可以得到二次收敛的速度。另一方面，如果$F^{''}(x)$是负定的，我们将得到二次收敛速度的最大解。我们可以使用下降的条件避免收敛到最大点。

下面是依据梯度方法和牛顿方法的混合方法，根据2.10可知，如果$F^{''}$是正定的，牛顿方法能够保证是下山方向。

$$\bbox[] { \begin{array}{} if \; F^{''}(x) 是正定的 \\ \quad h:= h_n; \\ else & \qquad (2.11)\\ \quad h:= h_{sd}; x:=x+\alpha h \end{array} }$$

这里$h_sd$是最速下降方向，$\alpha$是通过线性搜索找到。检查一个矩阵是不是正定的一个好的工具是使用乔里斯基方法(Cholesk’s method)，如果分解成功还可以用来求解线性方程组，因此对正定的检查是很方面的。

在2.4节我们将介绍同时求解搜索方向$h_n$和步长$\alpha$的方法。混合方法是非常有效的但它是很难被使用的。因为它需要计算$F^{''}(x)$，对于于一些复杂的应用是很困难的。相反我们使用拟牛顿方法(Quasi-Newton method)，它使用一系列的逐渐逼近$H^* = F^{''}(x)$的方法。

2.3 线性搜索

给定一个点$x$和下降方向$h$，现在要求解在$X$的$h$方向移动的步长。我们可以设

$$\bbox[] { \begin{array}{} \phi(\alpha) = F(x+ \alpha h), \qquad 其中x和h被固定，\alpha \ge 0. \qquad (2.12) \\ \phi^{'}(\alpha) = h^TF^{'}(x+ \alpha h ) \end{array} }$$

要保证$h$是下降的方向，则有

$$\phi^{'}(0) = h^TF^{'}(x) \lt 0$$

对于牛顿法，通常情况下，我们让$\alpha = 1$，对于 $\alpha$ 的选择可能引起一下几种情况

• $\alpha$ 很小，使得目标函数增益很小，$\alpha$ 应该增加。
• $\alpha$ 很大，大到$\phi(\alpha) \ge \phi(0)$，这时应该减少 $\alpha$ 使得满足下降方向。
• $\alpha$接近$\phi(\alpha)$的最小值，接受这个$\alpha$值。

精确地线性搜索是求解一系列的$\alpha_1, \alpha_2,\cdots$的迭代过程，我们的目标是寻找最小的$\alpha_e$(2.7)。当迭代步长$\alpha_s$满足下面的条件时停止

$$\vert \phi^{'}(\alpha_s) \le \tau \vert \phi^{'}(0)$$

其中$\tau$是非常小的正数。其中

$$\phi(\alpha_k) = F(x+\alpha_k h), \qquad \phi^{'}(\alpha_k) = h^TF^{'}(x+\alpha h)$$

精确的线性搜索可能会更多的计算时间，那是因为当$x$远离$x^*$时，搜索方向$h$可能偏离$x^* - x$更远。因此没有必要精确寻找$\phi$的最小。如果$\alpha$的值不是前两种情况我们就接受它，这种方法也就叫做松散的线性搜索(soft line search)。我们给出严格的版本关于(2.1)描述的下降条件

$$\phi(\alpha_s) \le \phi(0) + \gamma_1 \cdot \phi^{'}(0) \cdot \alpha \qquad 其中 0 \lt \gamma_1 \lt 1. \qquad (2.13a)$$

这就保证了不会出现第二种情况，第一种情况中$(\alpha, \phi(\alpha))$很接近起始点的切向，我们给出下面的条件

$$\phi^{'}(\alpha_s) \ge \gamma_2 \cdot \phi^{'}(0), \qquad 其中\gamma_1 \lt \gamma_2 < 1. \qquad (2.13b)$$

如果开始的猜想的$\alpha$满足这些约束，我就就接受它，否则我们必须以精确线性搜索为基础迭代搜索。

2.4 信頼域和阻尼方法

假设$L$有如下的模型

$$F(x+h) \simeq L(h) \equiv F(x) + h^Tc + \frac{1}{2}h^TBh \qquad (2.14)$$

其中$c \subseteq R^n$, 矩阵$B \subseteq R^{n \times n}$ 是对称的。$L(h)$是对$F(x+h)$的近似。在$h$充分小的情况下这个近似是正确的。

在信頼域中我们假设知道一个正的数字$\triangle$,在以$x$为中心，$\triangle$为半径的球内$L$模型是充分精确地，迭代步长被确定

$$h = h_{tr} \equiv argmin_{\Vert h \Vert \le \triangle} \{L(h)\} \qquad (2.15)$$

阻尼模型里迭代步长被确定

$$h = h_{dm} \equiv argmin_h\{ L(h) + \frac{1}{2}u h^Th\} \qquad (2.16)$$

其中阻尼参数$u \ge 0$。条件$\frac{1}{2}uh^Th = \frac{1}{2}u \Vert h \Vert ^2$被看做惩罚大步骤。

算法2.4的中心部分若以上述方法为基础则有

$$\bbox[] { if(F(x+h) \lt F(x)) \\ \quad x:= x + h; \\ update \; \triangle \; or \; u }$$

如果$h$满足下降的方向(2.1)，则相应的$\alpha = 1$。否则，$\alpha = 0$，也就是说我们不移动$x$。然而，我们并不是停在$x$(除非$x = x^*$),我们修改$\triangle or u$以期待下次迭代的时候有个好的$h$。

因为$L(h)$是在$h$充分小的情况下对$F(x+h)$的很好的近似。下一步失败的原因是因为$h$步长的太大，此时我们应该减小它。更进一步说，如果我们接受了这个步长，那么在下次迭代的时候我们可以使用更大的步长,也就是说在到达最小点$x^*$之前根据需要减小步长。

计算步长模型的质量可以使用增益率来评估

$$\rho = \frac{F(x) - F(x+h)}{L(0) - L(h)} \qquad (2.18)$$

也是说计算函数值真实的减少与预测的减少之间的比率。在这个式子中，分母是正的，如果分子是负的说明这一步不是下山方向(步长太大了应该减小步长)。

在信頼域中我们一直监测$\triangle$的步长，下面的更新策略被广泛使用的

$$\bbox[] { \begin{array}{} if \; \rho \lt 0.25 \\ \quad \triangle := \triangle / 2 ; \\ elseif \; \rho \gt 0.75 & \qquad \qquad (2.19)\\ \quad \triangle := max\{ \triangle , 3 * \Vert h \Vert \} \end{array} }$$

因此，如果$\rho \lt \frac{1}{4}$，我们使用更小的步长，而$\rho \gt \frac{3}{4}$ 说明我们可以使用更大的步长。信頼域算法对阈值0.25,0.75，除数 $p_1 = 2$或者因子$p_2 = 3$的微小变化不敏感。但是$p_1 = 2, p_2 = 3$的选择是非常重要的，它使得$\triangle$值不震荡。

在阻尼方法中，如果$\rho$很小表明应该加大阻尼因子也就是增加惩罚步长。如果$\rho$是一个大的值表明$L(h)$很好的近似了$F(x+h)$,此时阻尼因子应该减小。下面的策略与(2.19)一样都是被广泛使用的，它最初由马夸尔克(Marquardt)提出

$$\bbox[] { \begin{array}{} if \; \rho \lt 0.25 \\ \quad u:= u * 2; $ \qquad & \qquad & (2.20)\\ elseif \; \rho \gt 0.75 \\ \quad u:= u/3; \end{array} }$$

这个方法与上一个方法一样对阈值0.25, 0.75, $p_1=2, p_2=3$的微小变化不敏感，$p_1,p_2$的选择很重要，它使得$u$值不发生震荡。经验表明通过阈值0.25,0.75的不连续变换可能引起”抖动”的减缓收敛。下面是由尼克森(Nielsen)提出的总体性能较优越的策略

$$\bbox[] {\begin{array}{} if \; \rho \gt 0 \\ \quad u:= u * max\{ \frac{1}{3} , 1-(2 \rho -1)^3 \}; \; v:=2; \qquad (2.21)\\ else \\ \quad u:= u * v; \; v: = 2 * v; \end{array} }$$

2.4.1 计算步长

通过函数的稳定点求解阻尼因子的步长。

$$\varphi_u(h) = L(h) + \frac{1}{2}uh^Th = F(x) +h^Tc+\frac{1}{2}h^TBh+ \frac{1}{2}uh^Th$$

要求 $\varphi_u(h)$最小的解，则意味着$h_{dm}$是

$$\varphi^{'}_u(h) = L^{'}(h) + u h = 0$$
的解( 参考2.16),从 (2.14)的定义可以看到它等价于
$$(B+u I )h_{dm} = -c \qquad (2.22)$$

其中$I$是单位阵，如果$u$足够大，则对称矩阵$B+uI$是正定的。根据1.8的理论可知$h_{dm}$是$L$的局部最小。

例子 2.1 阻尼牛顿法

假设模型为上述的$L(h)$,其中$c = F^{'}(x)$， $B=F^{''}(x)$，根据(2.22)我们有

$$(F^{''}(x) + u I)h_{dm} = -F^{'}(x).$$
这里 $h_{dm}$就是所说的阻尼牛顿步，如果 $u$很大，有
$$h_{dm} \simeq - \frac{1}{u} F^{'}(x)$$

也即是在最速方向上的很小的步。另一方面，如果$u$很小，$h_dm$很接近牛顿的步$h_n$。因此，我们可以认为阻尼牛顿法就是最速下降大和牛顿法的组合。

信頼域方法的步$h_{tr}$通过下面的约束问题求解

$$\bbox[] { \begin{array}{} minimize \quad L(h) \\ subject \; to \quad h^Th \le \triangle ^2 . \qquad (2.23) \end{array} }$$

3. 非线性最小二乘

首先考虑下面的问题

给一个向量函数$f: R^n \mapsto R^m , m \gt n \;$现在要使$\Vert f(x) \Vert \;$最小。即：

$$x^*=argmin_x\{F(x)\} \qquad (3.1a) \;$$

或者等价于时$F(x)$最小,即：

$$F(x) = \frac{1}{2}\sum_{i=1}^m (f_i(x))^2 = \frac{1}{2} {\Vert f(x) \Vert}^2 = \frac{1}{2}f(x)^Tf(x)\,. \qquad (3.1b)$$

最小二乘问题可以通过一般的优化方法，我们将要提出的LMA方法更加有效率，很多情况下LMA要比线性方法要好，有些时候甚至是二次收敛的。

• 当初始解远离真实解时，表现的像是一个梯度下降法(线性搜索)，速度很慢，但是是可收敛的。

• 当初始解很接近真实解时，LMA表现的像是高斯牛顿方法，它的收敛速度二次的。

对于求解$F(x)$的最小值，首先需要对$F(X)$求导，这里假设$f$有连续的二街偏导数。我们可以利用泰勒的阶段二项式：

$$f(x+h)=f(x)+Jh+O\Vert h \Vert^2 \, , \, J \subseteq R^{m \times n} \qquad (3.2a);$$ 是雅克比矩阵，$J$由$f$的一阶偏导构成，即：

$$(J)_{ij}={\frac{\partial f_i}{\partial x_j}}(x) \qquad (3.2b) \quad ,$$
$$J(x) = \begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n } \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n } \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \frac{\partial f_m}{\partial x_2} & \cdots & \frac{\partial f_m}{\partial x_n} \end{pmatrix}_{m \times n} \qquad (3.2c)$$

对$F: R^n \mapsto R \, ,$求一阶导，即：

$$\frac{ \partial F }{\partial x} = \sum_{i=1}^m(f_i \frac{\partial f_i}{\partial x_j}(x)) \qquad (3.3)$$

因此$F(x)$的梯度为：

$$F^{'}(x) = J^T f(x) \qquad (3.4a)$$
$$F^{''}(x)= J^TJ+\sum_{i =1}^m(f_if_i^{''}(x)) \qquad (3.4b)$$

---

例3.1 线性最小二乘问题

若$f(x) = b - Ax \;$，其中$b \subseteq R^m , A \subseteq R^{m \times n} \;$，这是一个简单的线性最小二乘问题。其中：$f^{'}(x) = -A \;$，从1.4b可知

$$F^{'}(x) = -A^Tf(x) = -A^T(b-Ax)$$
假设 $F^{'}(x^*)=0 \;$，则有： $$A^TAx^*=A^Tb \qquad (1.5)$$
上式又叫正规方程组，为了解这个问题可以写成： $$Ax^* \simeq b \; ,$$
我们可以通过正交矩阵 $Q$求解$x^*$,像下面这样： $$Q^TA=\begin{pmatrix} R \\ 0 \\ \end{pmatrix}$$ 其中 $R$是上三角矩阵，可以通过回带法求解$x^*$, $$Rx^*=(Q^Tb)_{1:n}$$

*Note: 用正交矩阵求得的解要比解正规方程组得到的解精确的多。

例3.2 求正规方程的解

$f(x^*) = 0 \;$, where $f:R^n \mapsto R^n$

我们可以使用牛顿迭代法，牛顿迭代法的条件是首先有一个初始解的猜想$x_0$($x_0$很接近真实解)，然后依次计算$x_1,x_2,\ldots$，所有这些$x$必须满足$f(x+h)=0, O(\Vert h \Vert) \;$,是可以忽略。求解过程如下：

$$J(x_k)h_k=-f(x_k) \qquad (3.6)$$ $$x_{k+1}=x_k+h_k$$

这里如果$J(x^*)$是非奇异的，牛顿迭代法的收敛速度是二次的。也就是如果$d_k=\Vert x_k-x^* \Vert$是很小的，$\Vert x_{k+1} - x^* \Vert=O(d^{2}_k)$,如果$x$s是远离$x^*$牛顿法可能不收敛或者收敛到一个错误的值。

---

3.1 高斯牛顿法

当初始解很接近真实解，高斯牛顿法是一种十分有效的方法，它的收敛速度是二次的，高斯牛顿法的实现以向量函数的一阶导为基础。同时高斯牛顿法在求解非线性问题时，也是把它局部化线性处理。由$3.2a$知,当$\Vert h \Vert^2$很小时我们可以得到：

$$f(x+h) \simeq l(h)=f(x)+Jf(x) \qquad (3.7a)$$ 将 3.7a带入 3.1b有： $$F(x+h)\simeq L(h)\equiv \frac{1}{2}l(h)^Tl(h)= \frac{1}{2}f^Tf+h^TJ^Tf+\frac{1}{2}h^TJ^TJh =F(x)+h^TJ^Tf+\frac{1}{2}h^TJ^TJh \qquad (3.7b)$$

最小化$L(h)$得到高斯牛顿步长h_gn，

$$h_gn=argmin_h\{L(h)\}$$ 我们很容易看到 $L$的梯度和黑塞矩阵(Hession martrix)，即：
$$L^{'}(h)=J^Tf+J^TJh \qquad L^{''}(h)=J^TJ \qquad 3.8$$

与3.4a比较可知$L^{'}(0)=F^{'}(x)$。从上面的方程$L^{''}(h)$与$h$无关的并且是对称的。这里如果$J$是满秩，也就是说$J$的列是无关的，那么$L^{''}(h)$是正定的，这也说明了$L(h)$是唯一的最小值，我们可以通过下面的方程求解：

$$(J^TJ)h_{gn}=-J^Tf \qquad (3.9)$$

因此对于$F$来讲$h_{gn}$是最速下降的方向，因为：

$$h_{gn}F^{'}(x)=h_{gn}(J^Tf)=-h_{gn}(J^TJ)h_{gn} < 0 \qquad (3.10)$$

我们可以用下面的方程迭代求解:

$$(J^TJ)h_{gn}=-Jf$$ $$x := x+\alpha h_{gn}$$ 其中 $\alpha$被用来表示线性搜索，典型的高斯牛顿方法在所有的迭代步骤中使用 $\alpha = 1$，如果线性高斯牛顿法满足下面的条件则是二次收敛。

• $\{x \subseteq | F(x)\le F(x_0)\}$

• 在所有的迭代步骤中$J$是满秩。

Note:基于优化的牛顿法通常是二次收敛的，但是高斯牛顿法通常不是这种情况。

下面我们比较这两种方法的搜索方向：

$$F^{''}(x)h_n=-F^{'}x \quad and \quad L^{''}(h)h_{gn}=-L^{'}(0)$$ 根据 3.8a 我们可以看出方程的右边是相等的。根据3.4b 和3.8 可以看出来他们的系数矩阵的不同：

$$F^{''}(x) = L^{''}(x)+\sum_{i=1}^mf_i(x)f^{''}_i(x)$$ 因此，当$f(x^*)=0$，$x$很接近$x^*$时有$F^{''}(x)=L^{''}(x)$，这个时候高斯牛顿法二次收敛。如果函数$f_i$曲率很小或者$\vert f_i \vert$的值很小，则是超线性收敛，其他情况，我们希望它是线性收敛。

Note:$F(x^*)$的值决定了收敛的速度。

例子3.3 考虑下面的问题，其中$n=1,m=2$

$$f(x)=\begin{bmatrix} x+1 \\ \gamma x^2+x-1 \end{bmatrix} \qquad F(x)=\frac{1}{2}(x+1)^2+\frac{1}{2}(\gamma x^2+x-1)^2$$

接着：

$$F^{'}(x)=2 \gamma ^2x^3+3 \gamma x^2-2(\gamma -1)x$$ 因此$x=0$是$F$的稳定点，接着： $$F^{''}(x)=6 \gamma ^2x^2+6 \gamma x - 2(\gamma -1)$$
这就表明如果$\gamma \lt 1$,那么$F^{''}(0)\gt 0$,因此$x=0$是一个局部最小解(更准确的说它是一个全局最小解)。下面在看看雅克比矩阵： $$J(x)= \begin{bmatrix} 1 \\ 2 \gamma x +1 \end{bmatrix}$$

高斯牛顿方法的迭代步骤：

$$x_{k+1} = x_{k} - \frac{2 \gamma ^2 x_{k}^3 + 3 \gamma x_{k}^2 - 2(\gamma -1)x_k}{2+4 \gamma x_k + 4 \gamma ^2 x_{k}^2}$$
如果$\gamma \neq 0$且$x_k$趋近于0则有： $$x_{k+1} = x_k+(\gamma -1)x_k+O(x^{2}_k) =\gamma x_k+O(x^{2}_k)$$ 因此，如果$\vert \gamma \vert \lt 1$则是线性收敛，如果$\gamma \lt -1$则高斯牛顿法不能找到最小值。例如$\gamma = -2 \; and \; x_0=0.1$其迭代过程很混乱。如下表所示：

$$\begin{array}{c|l} k & x_k \\ \hline 0 & 0.1000 \\ 1 & -0.3029 \\ 2 & 0.1368 \\ \vdots & \quad \vdots \end{array}$$

最后，如果$\gamma = 0$则有$x_{k+1}=x_k-x_k=0$,这里一步就找到了解，因为这时$f(x)$退化成线性函数。

例子3.4牛顿迭代法和高斯牛顿法

回想例子1.2,$f(x^*)=0$,其中$f:R^n \mapsto R^n$。如果使用牛顿迭代法去解决这个问题，迭代步骤如下：

$$J(x)h_{nr} = -f(x) \quad x := x + h_{nr}$$

高斯牛顿法最小化$F(x)=\frac{1}{2}f^T(x)f(x)$有如下的步骤：

$$(J^T(x)J(x))h_gn = -J^T(x)f(x); \quad x:=x+h_gn$$
这里$J(x)$是方阵，若假设$J(x)$是非奇异的且$(J^T(x))^{-1}$存在，高斯牛顿方法和牛顿迭代法具有相同的问题。也就是初始解的问题。

3.2 The Levernberg-Marquardt Method

Leverberg-Marquardt Method(LMA)中文称作列文伯格-马夸尔特法，也常被叫做阻尼最小二乘法(damped least-squares method)。常被用来解决非线性最小二乘问题，特别是在最小二乘曲线拟合问题(least squares curve fitting)。LMA算法其实就是下降法与高斯牛顿的结合，当初始解远离真实解时，LMA算法表现的像下降法，当初始解接近真实解时LMA算法表现的像是高斯牛顿法。下面定义了$h_{lm}$：

$$(J^TJ+\mu I)h_{lm}=-g \; , \; with \quad g=J^Tf \; , \mu \ge 0 \qquad (3.13)$$

这里，$J=J(x) \; , f = f(x)$。其中阻尼参数$\mu$的作用如下：

• 如果$\mu \ge 0$,则系数矩阵是正定的，这就保证了$h_{lm}$是下降方向

• 如果$\mu$的值较大我们有
  

  $$h_{lm} \simeq -\frac{1}{u}F^{'}(x)$$

也就是在最速下降的方向减少步长，这也说明了如果当前的迭代解如果远离真实解时这是一个很好的确保正确的方法。

• 如果$\mu$是很小的。则有$h_{lm}= h_{gn}$,这对最终的真实解来说是一个很好地迭代步长。这时候如果$F^{'}(x^*) = 0$(或者$F^{'}{x^*}$非长的小)，我们可能得到二次的收敛速度。

这就是阻尼擦数对方向和步长的影响。这也给了我们一个求解的方法，即使不是线性搜索的时候。$\mu$值得选择应该与$A_0=J^{T}(x)J(x_0)$的元素有关，例如让

$$u_0 = \tau \cdot max_i{a^{0}_{ii}}$$

这里$\tau$是用户自定义的(此算法对于$\tau$的选取不敏感，但是一般的经验都是使用很小的值，例如如果$x_0$很接近真实解时$\tau=10^{-6}$)，否则，$\tau = 10^{-3}$或者$\tau = 1$.)。控制$\tau$更新的条件是$\rho$，

$$\rho = \frac{F(x) - F(x+h_{lm})}{L(0)-L(h_{lm})}$$

其中分母是对线性模型(1.7b)增益的预测

$$L(0) - L(h_{lm}) = -h^{T}_{lm}- \frac{1}{2}h^{T}_{lm}J^TJh_{lm} =- \frac{1}{2}h^T_{lm}\left(2g+\left(J^TJ+\mu J- \mu J\right)h_{lm}\right) = \frac{1}{2}h^T_{lm}(\mu h_{lm}-g)$$

因为$h^{T}_{lm} \; , -h^{T}_{lm}g$都是正的，所以保证了$L(0)-L(h_{lm})$是正的。

如果$\rho$的值很大，说明$L(h_{hm})$能够很好的近似$F(x+h_{lm})$，此时我们可以减少阻尼因子$\mu$,使得下一次LM步长更加接近高斯牛顿步长。如果$\rho$很小(也有可能是负的),那么$L(h_{lm})$不能够近似$F(x+h_{lm})$,此时我们应该增加$\mu$,也就是使得$\mu = 2 \mu$以此让$\mu$更加接近最速下降的方向和减少步长。

算法的停止迭代条件应该反映到全剧最小，我们有$F^{'}(x^*) = g(x^*) = 0$，因此我们可以使用

$$\Vert g \Vert _{\infty} \le \epsilon_1 \qquad (3.15a)$$

其中$\epsilon$是用户设定的非常小的正数，另一个判断停止迭代的条件是利用$x$非常小的变换。

$$\Vert x_{new} - x\Vert \le \epsilon_2 (\Vert x\Vert + \epsilon_2) \qquad (3.15b)$$

这里假设如果$x$趋近于0，当$\Vert x |Vert$大到时$\epsilon_2^{2}$时，这个表达式表示了$x$缓慢的变化。最后在一个有限循环中迭代终止的条件

$$k \gt k_{max} \qquad (3.15c)$$

$\epsilon_2 \; , k_{max}$ 都是用户设定的值，$\epsilon_2 \; , k_{max}的作用在于如果

$\epsilon_1$设置的比较小以至于对舍入误差有很大的影响停止迭代，这也表明了$F$的增益也$L$的增益不一致的情况，也就是说可能出现这么个情况：$\mu$增加的很快同时$\Vert h_{lm} \Vert$很小，这时(3.15b)将要起作用。该LM算法总结如下

$Algorithm Levenberg-Marquart method$
$begin$
$\quad k := 0; \; \nu := 2; \; x := x_0;$
$\quad A := J^{T}(x)J(x); \; g:=J^{T}(x)f(x);$
$\quad found := (\Vert g \Vert_{\infty} \le \epsilon_1); \; \mu = \tau * max\{a_{ii}\};$
$\quad while \; (not found) \; and \; ( k \lt k_{max})$
$\qquad k:=k+1; \; Solve (A+\mu I)h_{lm}=-g;$
$\qquad if \; \Vert h_{lm} \Vert \le \epsilon_2(\Vert x \Vert + \epsilon_2)$
$\qquad \quad found := true;$
$\qquad else$
$\qquad \quad x_{new} := x + h_{lm};$
$\qquad \quad \rho:=(F(x) - F(x_{new})) /(L(0)-L(h_{lm}));$
$\qquad \quad if \; \rho \gt 0$
$\qquad \qquad x:= x_{new};$
$\qquad \qquad A:=J^T(x)J(x); \; g:=J^T(x)f(x);$
$\qquad \qquad found := (\Vert g \Vert_\infty \le \epsilon_1);$
$\qquad \qquad \mu =: \mu * max\{\frac{1}{3},1-(2\rho-1)^3\}; \; \nu :=2;$
$\qquad \quad else$
$\qquad \qquad \mu:=\mu * \nu; \; \nu:= 2 * \nu ;$
$end$

参考文献

[1] Method for non-linear least squares problems. 2nd Edition, April 2004. K.Madsen, H.BNielsen, O.Tingleff.
[2] A Brif Description of the Levenberg-Marquardt Algorithm Implemented by levmar.
[3] 科学计算导论 .