Numerical Optimization Ch10. Least-Squares Problems

第十章: 最小二乘问题

在 最小二乘问题(Least-square problems)中, 目标函数 $f$具有如下特殊形式: $$f(x)=\frac{1}{2}\sum _{j=1}^m{r}_j^2(x),$$其中每个 $r_j:{\mathbb{R}}^n\to \mathbb{R}$为光滑函数. 我们称每个 $r_j$为 残差(residual). 本章中假设 $m\ge n$.

最小二乘问题出现在许多应用领域中, 并且(事实上)可能是无约束优化问题最大的源头. 医药、物理、金融等领域的研究者构建参数模型都会考虑上面$f$的形式来表征模型与观测之间的差距. 通过极小化这一函数, 他们就能获得使模型最佳拟合数据的参数. 本章我们将展示如何通过探究$f$及其导数的特殊结构, 设计高效的、鲁棒的极小化算法.

上述优化问题比之一般的无约束极小化问题, 具有最优值非负的特点. 其次, 上述$f$的特殊结构能使得最小二乘问题比一般的无约束极小化问题更易求解. 首先将每个残差组分$r_j$组装成一个残差向量$r:{\mathbb{R}}^n\to {\mathbb{R}}^m$:$$r(x)=(r_1(x),r_2(x),\dots ,r_m(x){)}^T.$$利用这一表示, 我们可以将$f$写作$f(x)=\frac{1}{2}\Vert r(x){\Vert }_2^2$. $f$的导数就可以$m\times n$Jacobi矩阵$J(x)$的形式表示:$$J(x)={\left[\frac{\partial r_i}{\partial x_j}\right]}_{i=1,2,\dots ,m,j=1,2,\dots ,n}=\left[\begin{array}{c}\nabla r_1(x{)}^T\\ \nabla r_2(x{)}^T\\ ⋮\\ \nabla r_m(x{)}^T\end{array}\right],$$其中每个$\nabla r_j(x),j=1,2,\dots ,m$为$r_j$的梯度. 于是$f$的梯度与Hessian矩阵为:$$\begin{array}{rl}\nabla f(x)& =\sum _{j=1}^m{r}_j(x)\nabla r_j(x)=J(x{)}^{T}r(x),\\ {\nabla }^{2}f(x)& =\sum _{j=1}^m\nabla r_j(x)\nabla r_j(x{)}^T+\sum _{j=1}^m{r}_j(x){\nabla }^2{r}_j(x)\\ & =J(x{)}^{T}J(x)+\sum _{j=1}^m{r}_j(x){\nabla }^2{r}_j(x).\end{array}$$一般, 残差的一阶偏导(从而Jacobi矩阵$J(x)$)相对容易(或相对便宜)计算. 从而$\nabla f(x)$的表达式是可使用的. 而有了$J(x)$, 我们也可以计算${\nabla }^{2}f(x)$中的$J(x{)}^{T}J(x)$. 这一过程完全不需要计算$r_j$的任何二阶导数. ${\nabla }^{2}f(x)$的这一部分"免费"可用性是最小二乘问题的关键特征. 进一步地, $J(x{)}^{T}J(x)$往往要比第二项来得更重要. 这可能是由于在解附近

• 残差$r_j$接近于线性(也就是说, ${\nabla }^2{r}_j(x)$相对较小); 或者
• 残差较小(即$r_j(x)$相对较小).

绝大多数非线性最小二乘的算法均会探究和利用Hessian的结构性质.

求解最小二乘问题最广为使用的算法均以之前介绍的线搜索和信赖域为框架. 它们主要基于牛顿和拟牛顿法, 其中会考虑$f$的特殊结构.
本章结构如下: 第1节涵盖最小二乘的一些应用背景; 第2节介绍线性最小二乘问题, 这将启发我们对非线性情形的算法的讨论; 第3节将介绍主要的算法; 第4节概括地介绍最小二乘的一种变体——正交距离回归(也称总体最小二乘). 第五章小谈大规模问题.
本章如不说明, 我们以$\Vert \cdot \Vert$表示欧式范数$\Vert \cdot {\Vert }_2$.

1. 背景介绍

我们引入一个简单的带参模型以展示最小二乘方法将如何帮助我们选择能最佳拟合观测数据的参数.

例1 现在我们想要研究一种特定药物在病人身上的药效如何. 在病人注射药物后, 我们在特定的时间点抽取血液样本、测定样本中药物的凝聚程度, 最终得到时间$t_j$与凝聚度$y_j$构成的数据表.

基于我们过往的经验, 我们发现函数$\varphi (x;t)$对药物在$t$时刻的凝聚度有很好的预测能力, 其中$x$为五维参数向量$x=(x_1,x_2,x_3,x_4,x_5)$:$$\varphi (x;t)=x_1+t{x}_2+t^2{x}_3+x_4{e}^{-x_{5}t}.$$我们待定$x$, 最终需要我们的模型以某种方式最佳匹配拟合我们的观测数据. 一种较好的表示预测模型值域观测值差距的方式就是如下最小二乘函数:$$\frac{1}{2}\sum _{j=1}^m[\varphi (x;t_j)-y_j{]}^2.$$我们定义$r_j(x)=\varphi (x;t_j)-y_j$. 几何上, 每个$|r_j|$均表示点$(t_j,y_j)$与曲线$\varphi (x;t)$(视作$t$的函数, $x$为固定的参数向量)的垂直距离. 可见下图.

最小二乘问题的极小点$x^{\ast }$就使得图中虚线长度平方和极小. 有了$x^{\ast }$, 我们就可以使用$\varphi (x^{\ast };t)$预测$t$时刻病人血液中药物的凝聚程度.

这是固定回归模型(fixed-regressor model)的一个例子: 它假设抽取血液样本的时间点$t_j$具有高精度, 而观测$y_j$则可能(由于仪器或实验人员的限制)或多或少包含随机误差.
在一般的如刚才描述的数据拟合问题中, 模型$\varphi (x;t)$中的坐标$t$还可能是向量. 例如刚才问题中的$t$, 还可以包含病人的其他指标, 例如身高、体重等. 指标涵盖得越全面, 理论上能说明的规律就越丰富.
平方和并不是唯一度量差异的方式. 其他常用的方法还有:

• 最大绝对值. $$\underset{j=1,2,\dots ,m}{\max }|\varphi (x;t_j)-y_j|.$$
• 绝对值和. $$\sum _{j=1}^m|\varphi (x;t_j)-y_j|.$$

利用$l_{\infty }$和$l_1$范数, 我们可以把这两种函数分别写作$$f(x)=\Vert r(x){\Vert }_{\infty },f(x)=\Vert r(x){\Vert }_1.$$我们将在后面的章节说明如何将上述重构为光滑的约束优化问题. 而本章中我们仅讨论$l_2$范数下的问题.

有时, 选取最小二乘标准也有统计上的动机. 我们稍微改变下记号, 令${ϵ}_j$表示模型与观测之间的差距, 即$${ϵ}_j=\varphi (x;t_j)-y_j.$$通常我们假设 { ϵ j } \{\epsilon_j\} {ϵj​}相互独立且同分布(independent and identically distributed, i.i.d.), 它们的概率密度函数为$g_{\sigma }(\cdot )$, 具有一定的方差${\sigma }^2$. 这一假设一般是符合实际的. 例如当模型能够准确反映过程, 且当测定$y_j$的误差不包含系统误差时. 基于这样的假设, 给定参数向量$x$, 特定观测集 { y j } \{y_j\} {yj​}出现的似然为$$p(y;x,\sigma )=\prod _{j=1}^m{g}_{\sigma }({ϵ}_j)=\prod _{j=1}^m{g}_{\sigma }(\varphi (x;t_j)-y_j).$$给定观测$y_1,y_2,\dots ,y_m$, $x$"最可能"的值就是使得$p(y;x,\sigma )$(作为$x$的函数)最大的位置. 此时得到的值称作极大似然估计(maximum likelihood estimate).
当我们假设差异服从正态分布时, 我们有$$g_{\sigma }(ϵ)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{{ϵ}^2}{2{\sigma }^2}\right).$$代入$p(y;x,\sigma )$可得$$p(y;x,\sigma )=(2\pi {\sigma }^2{)}^{-m/2}\exp \left(-\frac{1}{2{\sigma }^2}\sum _{j=1}^m[\varphi (x;t_j)-y_j{]}^2\right).$$对任意固定方差${\sigma }^2$, 取对数可得$p$最大当且仅当残差平方和极小. 这就是说, 当我们假设差距i.i.d.且服从同一个正态分布时, 极大似然估计与极小化残差平方和得到的极小点是一样的. 满足条件的对${ϵ}_j$的假设并不唯一.

进一步我们还有推广形式的目标函数$$r(x{)}^{T}Wr(x),$$其中$W\in {\mathbb{R}}^{m\times m}$对称. 这可以看做是加权最小二乘(weighted least-square problems).

2. 线性最小二乘问题及求解算法

许多数据拟合问题中的模型函数$\varphi (x;t)$是$x$的线性函数. 此时, 残差$r_j(x)$也是线性的, 称相应的问题是线性最小二乘问题. 进一步地, 我们可将残差向量写作$r(x)=Jx-y$, 其中$J$为Jacobi矩阵, $y$为观测(以及可能存在的截距)构成的向量, 二者均独立于$x$存在, 从而目标函数为$$f(x)=\frac{1}{2}\Vert Jx-y{\Vert }^2,$$其中$y=r(0)$. 我们有$$\nabla f(x)=J^T(Jx-y),{\nabla }^{2}f(x)=J^{T}J.$$(此时原本${\nabla }^{2}f(x)$的第二项因${\nabla }^2{r}_j=0,j=1,2,\dots ,m$而不计.) 易知上述定义的$f(x)$是凸的——这对于一般的非线性问题并不必要. 因此, 任意满足$\nabla f(x^{\ast })=0$的$x^{\ast }$均是$f$的全局极小点. 这等价于$x^{\ast }$必须满足如下线性方程组:$$J^{T}J{x}^{\ast }=J^{T}y.$$这称作是对应于$f$的正规方程组.

对于无约束的线性最小二乘问题, 我们介绍三种主要的算法. 在讨论的过程中, 大部分时间我们均假设$m\ge n$以及$J$列满秩.

2.1 基于Cholesky分解的直接法

最直接的求解正规方程组的方法分以下三步:

1. 计算系数矩阵$J^{T}J$以及右端项$J^{T}y$;
2. 计算对称矩阵$J^{T}J$的Cholesky分解;
3. 经前代、回代求得解$x^{\ast }$.

这里, Cholesky分解$$J^{T}J={\bar{R}}^T\bar{R}$$(其中$\bar{R}$为$n\times n$的上三角矩阵, 其对角元为正数)在$m\ge n$且$J$列满秩时必定存在. 此法在实际中被频繁使用且通常比较高效, 只是有一个重要的缺陷: $J^{T}J$的条件数为$J$的条件数的平方. 由于一个问题计算解的相对误差通常与条件数成比例, 因此在求解精度上基于Cholesky分解的方法要逊于不增加条件数的方法. 特别地, 当$J$足够病态, Cholesky分解便无法实施, 这是因为舍入误差可能导致分解过程中对角元出现负值.

2.2 基于QR分解的直接法

第二种方法基于对矩阵$J$的QR分解. 由于正交变换保长, 所以$$\Vert Jx-y\Vert =\Vert Q^T(Jx-y)\Vert ,$$其中$Q$为任意$m\times m$的正交阵. 假设我们对$J$有带列选主元的QR分解:$$J\Pi =Q\left[\begin{array}{c}R\\ O\end{array}\right]=\left[\begin{array}{cc}{Q}_1& Q_2\end{array}\right]\left[\begin{array}{c}R\\ O\end{array}\right]=Q_{1}R,$$其中

• $\Pi$为$n\times n$的排列矩阵(因此正交);
• $Q$为$m\times m$的正交阵;
• $Q_1$为$Q$的前$n$列, $Q_2$为$Q$的后$m-n$列;
• $R$为对角元为正的$n\times n$上三角阵;

因此我们有$$\begin{array}{rl}\Vert Jx-y{\Vert }_2^2& ={\Vert \left[\begin{array}{c}{Q}_1^T\\ Q_2^T\end{array}\right](J\Pi {\Pi }^{T}x-y)\Vert }_2^2\\ & ={\Vert \left[\begin{array}{c}R\\ O\end{array}\right]({\Pi }^{T}x)-\left[\begin{array}{c}{Q}_1^{T}y\\ Q_2^{T}y\end{array}\right]\Vert }^2\\ & =\Vert R({\Pi }^{T}x)-Q_1^{T}y{\Vert }_2^2+\Vert Q_2^{T}y{\Vert }^2.\end{array}$$对于上面的第二项我们无能为力, 不过我们可以把第一项消成0, 即令$$x^{\ast }=\Pi R^{-1}{Q}_1^{T}y.$$(实际操作时, 我们先回代求解$Rz=Q_1^{T}y$, 再排列$z$的元素得到$x^{\ast }=\Pi z$.)

基于QR分解的直接法并不会恶化问题的条件数, 最终计算解$x^{\ast }$上的相对误差通常是与$J$的条件数成比例, 而不是其平方, 从而我们可以说这种方法(相较于第一种方法)是数值稳定的.

2.3 基于奇异值分解(SVD)的直接法

在一些特殊情形下, 我们要求解对数据$(J,y)$上的扰动更加鲁棒. 此时基于$J$的SVD的方法就可派上用场. 矩阵$J$的SVD为:$$J=U\left[\begin{array}{c}S\\ O\end{array}\right]V^T=\left[\begin{array}{cc}{U}_1& U_2\end{array}\right]\left[\begin{array}{c}S\\ O\end{array}\right]V^T=U_{1}S{V}^T,$$其中

• $U$为$m\times m$的正交阵;
• $U_1$为$U$的前$n$列, $U_2$为$U$的后$m-n$列;
• $V$为$n\times n$的正交阵;
• $S$为$n\times n$对角阵, 其中对角元${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_n>0$.

(注意$J^{T}J=V{S}^2{V}^T$, 因此$V$的列就是$J^{T}J$的对应于奇异值${\sigma }_j^2,j=1,2,\dots ,n$的特征向量.) 类似地我们有$$\begin{array}{rl}\Vert Jx-y{\Vert }^2& ={\Vert \left[\begin{array}{c}S\\ O\end{array}\right](V^{T}x)-\left[\begin{array}{c}{U}_1^T\\ U_2^T\end{array}\right]y\Vert }^2\\ & =\Vert S(V^{T}x)-U_1^{T}y{\Vert }^2+\Vert U_2^{T}y{\Vert }^2.\end{array}$$我们令第一项为0, 即有$$x^{\ast }=V{S}^{-1}{U}_1^{T}y.$$进一步地, 以$u_i\in {\mathbb{R}}^m,v_i\in {\mathbb{R}}^n$分别表示$U,V$的第$i$列, 我们有$$x^{\ast }=\sum _{i=1}^n\frac{u_i^{T}y}{{\sigma }_i}{v}_i.$$这一公式给予了我们关于$x^{\ast }$敏度的丰富的信息: 当${\sigma }_i$很小, $x^{\ast }$对于$y$中(以及$J$中)影响$u_i^{T}y$的部分就非常敏感. 这些信息在$J$几近亏秩时尤其有用, 也就是说当${\sigma }_n/{\sigma }_1\ll 1$时. 我们可以适当地采取一些防护措施避免数值上的不稳定. 有时以SVD的高昂计算量换取这些有效信息是值当的.

2.4 三种方法的讨论

以上三种方法均有它们的适用情形.

• 基于Cholesky分解的算法尤其适用于$m\gg n$的情形, 此时储存$J^{T}J$(而不仅仅是$J$)是可以接受的. 当$m\gg n$且$J$稀疏时, 其计算量也并不大. 然而此法在$J$亏秩或病态时必须经过修正, 以允许在$J^{T}J$的对角元上选主元.
• 基于QR分解的算法避免了条件数的爆发, 因此更加数值稳定.
• 基于SVD的算法尽管计算昂贵, 但确实最鲁棒也是最可靠的. 当$J$亏秩时, 一些奇异值${\sigma }_i$为0, 此时任一具有以下形式$$x^{\ast }=\sum _{{\sigma }_i̸=0}\frac{u_i^{T}y}{{\sigma }_i}{v}_i+\sum _{{\sigma }_i=0}{\tau }_i{v}_i$$(系数${\tau }_i$任意)的$x^{\ast }$均是问题的极小点. 通常我们最希望得到具有最小范数的解, 此时就令${\tau }_i=0$即可. 当$J$列满秩但病态时, 最后几个奇异值${\sigma }_n,{\sigma }_{n-1},\dots$相对于${\sigma }_1$较小. 在${\sigma }_i$较小时, 系数$u_i^{T}y/{\sigma }_i$对$u_i^{T}y$中的扰动尤其敏感. 因此我们可以直接忽略对那些敏感项的求和得到更加稳定的近似解.
• 当问题规模较大时, 使用迭代法求解正规方程组将更加高效, 例如共轭梯度法. 最直接的共轭梯度法的一次迭代仅需一次矩阵$(J^{T}J)$-向量乘积. 这一步可通过接连与$J,J^T$相乘得到. 至今已有许多共轭梯度法的修正版本, 它们的单步计算量并无大变, 但却具有更优越的数值性质. 例如Paige和Saunders提出的称为是LSQR的算法.

3. 求解非线性最小二乘问题的算法

3.1 Gauss-Newton法

3.1.1 Gauss-Newton法介绍

下面我们介绍目标函数非线性的情形. 我们将充分挖掘梯度$\nabla f$和Hessian矩阵${\nabla }^{2}f$的结构. 这其中最简单的算法——Gauss-Newton法——可视为线搜索框架下的修正牛顿法. 我们并不求解标准的牛顿方程${\nabla }^{2}f(x_k)p=-\nabla f(x_k)$, 转而求解$$J_k^T{J}_k{p}_k^{\mathrm{G}\mathrm{N}}=-J_k^T{r}_k$$获得相应的搜索方向$p_k^{\mathrm{G}\mathrm{N}}$. 这一简单的改动带来了许多好处:

1. ${\nabla }^2{f}_k\approx J_k^T{J}_k$的近似省去了我们计算每个残差的Hessian${\nabla }^2{r}_j,j=1,2,\dots ,m$的功夫. 事实上, 若我们在计算梯度$\nabla f_k=J_k^T{r}_k$的过程中就已经计算了Jacobi矩阵$J_k$的话, 这一近似是根本不会牵涉到任何的导数计算的. 这在某些应用上可以节省大量的计算时间.
2. 实际上我们经常看到第一项$J^{T}J$(相对于第二项)占主的场景(当然这得离解$x^{\ast }$足够近), 从而近似是得当的, Gauss-Newton法的收敛速度也并不会逊于Newton法太多. 具体说, 比如第二项中每一小项的范数(即$|r_j(x)|\Vert {\nabla }^2{r}_j(x)\Vert$)比$J^{T}J$的特征值要小得多. 在之前我们提到一种情形: 当残差$r_j$较小或者它们近似于线性时. 实际中, 许多最小二乘问题在解附近均有较小的残差, 从而保证了Gauss-Newton法的收敛速度.
3. 只要$J_k$列满秩以及梯度$\nabla f_k$非零, 得到的方向$p_k^{\mathrm{G}\mathrm{N}}$就是下降方向, 从而可用于线搜索中. 事实上, $$(p_k^{\mathrm{G}\mathrm{N}}{)}^T\nabla f_k=(p_k^{\mathrm{G}\mathrm{N}}{)}^T{J}_k^T{r}_k=-(p_k^{\mathrm{G}\mathrm{N}}{)}^T{J}_k^T{J}_k{p}_k^{\mathrm{G}\mathrm{N}}=-\Vert J_k{p}_k^{\mathrm{G}\mathrm{N}}{\Vert }^2\le 0.$$其中最后一个不等式只有在$J_k{p}_k^{\mathrm{G}\mathrm{N}}=0$的时候取等, 此时也应当有$J_k^T{r}_k=\nabla f_k=0$. 这就是说$x_k$已经是个稳定点了.
4. Gauss-Newton法所解的方程与线性情形的正规方程方程有一定的相似度. 具体说来, $p_k^{\mathrm{G}\mathrm{N}}$可以看做是(也实际上就是)以下线性最小二乘问题的解:$$\underset{p}{\min }\frac{1}{2}\Vert J_{k}p+r_k{\Vert }^2.$$因此, 我们可以用求解线性最小二乘问题的算法求解以上子问题. 进一步, 若我们使用基于QR分解或SVD的算法, 我们甚至不需要显式地计算出Hessian的近似$J_k^T{J}_k$. 若使用共轭梯度法也是一样: 我们只需计算矩阵$J_k^T{J}_k$-向量乘积, 而这一步可以通过先后对$J_k,J_k^T$操作得到.
5. 大规模情形. 若残差的数量$m$很大而变量数$n$相对较小, 此时显式地存储$J$似乎就显得不太符合情理. 不过我们可以通过连续计算$r_j,\nabla r_j,j=1,2,\dots ,m$再求和得到:$$J^{T}J=\sum _{j=1}^m(\nabla r_j)(\nabla r_j{)}^T,J^{T}r=\sum _{j=1}^m{r}_j(\nabla r_j).$$
   以上子问题还启发我们给出得到Gauss-Newton法的另外一种途径. 我们利用Taylor展开得到向量值函数的近似$r(x_k+p)\approx r_k+J_{k}p$. 因此$$f(x_k+p)=\frac{1}{2}\Vert r(x_k+p){\Vert }^2\approx \frac{1}{2}\Vert J_{k}p+r_k{\Vert }^2,$$再选取$p_k^{\mathrm{G}\mathrm{N}}$作为这一近似模型的极小点.

Gauss-Newton法的实施通常需要沿着$p_k^{\mathrm{G}\mathrm{N}}$做线搜索, 这其中需要步长参数${\alpha }_k$满足第三章中提到的条件, 例如Armijo条件、Wolfe条件.

3.1.2 Gauss-Newton法的收敛性

第三章中的理论可以用于研究Gauss-Newton法的收敛性质. 我们将利用Zoutendijk定理, 证明在假定Jacobi矩阵$J(x)$(这里的$x$落在需要研究的区域内)的所有奇异值一致远离(uniformly bounded away)0, 即$\exists \gamma >0$使得$$\Vert J(x)z\Vert \ge \gamma \Vert z\Vert ,\forall x\in \mathcal{N}$时(其中$\mathcal{N}$为水平集 L = { x ∣ f ( x ) ≤ f ( x 0 ) } \mathcal{L}=\{x|f(x)\le f(x_0)\} L={x∣f(x)≤f(x0​)}的一个邻域, 这里$x_0$为算法初始点. 我们称之为一致满秩条件), Gauss-Newton法的全局收敛性质. 本章从这里开始假设$\mathcal{L}$是有界的.

定理1 设每个残差函数$r_j$在水平集的一个邻域$\mathcal{N}$内Lipschitz连续可微, Jacobi矩阵$J(x)$在$\mathcal{N}$上满足一致满秩条件. 则由Gauss-Newton法产生的迭代序列 { x k } \{x_k\} {xk​}(其中步长参数${\alpha }_k$满足Wolfe条件)成立$$\underset{k\to \infty }{\lim }{J}_k^T{r}_k=0.$$
证明: 首先注意有界水平集$\mathcal{L}$的邻域可以选取得足够小使得对正常数$L,\beta$成立以下不等式:$$|r_j(x)|\le \beta ,\Vert \nabla r_j(x)\Vert \le \beta ,$$$$|r_j(x)-r_j(\tilde{x})|\le L\Vert x-\tilde{x}\Vert ,\Vert \nabla r_j(x)-\nabla r_j(\tilde{x})\Vert \le L\Vert x-\tilde{x}\Vert ,\forall x,\tilde{x}\in \mathcal{N},$$$$j=1,2,\dots ,m.$$易知存在常数$\bar{\beta }>0$使得$\Vert J(x{)}^T\Vert =\Vert J(x)\Vert \le \bar{\beta },\forall x\in \mathcal{L}$, 以及由于$\nabla f(x)={\sum }_{j=1}^m{r}_j(x)\nabla r_j(x)$, 从而$\nabla f$是Lipschitz连续的. 因此Zoutendijk定理的条件满足.
下面我们验证搜索方向$p_k^{\mathrm{G}\mathrm{N}}$和负梯度$-\nabla f_k$的夹角${\theta }_k$一致远离$\pi /2$. 事实上, $$\cos {\theta }_k=-\frac{(\nabla f{)}^T{p}^{\mathrm{G}\mathrm{N}}}{\Vert p^{\mathrm{G}\mathrm{N}}\Vert \Vert \nabla f\Vert }=\frac{\Vert J{p}^{\mathrm{G}\mathrm{N}}{\Vert }^2}{\Vert p^{\mathrm{G}\mathrm{N}}\Vert \Vert J^{T}J{p}^{\mathrm{G}\mathrm{N}}\Vert }\ge \frac{{\gamma }^2\Vert p^{\mathrm{G}\mathrm{N}}{\Vert }^2}{{\bar{\beta }}^2\Vert p^{\mathrm{G}\mathrm{N}}{\Vert }^2}=\frac{{\gamma }^2}{{\bar{\beta }}^2}>0.$$由Zoutendijk定理知$\nabla f(x_k)\to 0$, 得证.

若$J_k$(对某个$k$)亏秩, 此时一致满秩条件不成立, 系数矩阵$J_k^T{J}_k$是奇异的. $J_k^T{J}_{k}p=-J_k^T{r}_k$仍然有解, 不过有无穷多解, 其中每个都具有形式$$p=\sum _{{\sigma }_i̸=0}-\frac{u_i^T{r}_k}{{\sigma }_i}{v}_i+\sum _{{\sigma }_i=0}{\tau }_i{v}_i,$$这里${\tau }_i$任意. 但这样一来我们便无法保证$\cos {\theta }_k$一致远离0, 从而得不到如定理1一般的全局收敛性.

当$J_k^T{J}_k$(相对第二项而言)占主时, Gauss-Newton法向解$x^{\ast }$的收敛速度可以很快. 假设$x_k$距离$x^{\ast }$充分近, 且Jacobi矩阵$J(x)$满足一致满秩条件. 类似于牛顿法的分析, 对于单位步长的Gauss-Newton步, 我们有$$\begin{array}{rl}{x}_k+p_k^{\mathrm{G}\mathrm{N}}-x^{\ast }& =x_k-x^{\ast }-[J^{T}J(x_k){]}^{-1}\nabla f(x_k)\\ & =[J^{T}J(x_k){]}^{-1}[J^{T}J(x_k)(x_k-x^{\ast })+\nabla f(x^{\ast })-\nabla f(x_k)],\end{array}$$这里$J^{T}J(x)$为$J(x{)}^{T}J(x)$的简写. 以$H(x)$表示${\nabla }^{2}f(x)$表达式中的二阶项, 由Taylor定理可知$$\begin{array}{rl}\nabla f(x_k)-\nabla f(x^{\ast })& ={\int }_0^1{J}^{T}J(x^{\ast }+t(x_k-x^{\ast }))(x_k-x^{\ast })\mathrm{d}t\\ & +{\int }_0^{1}H(x^{\ast }+t(x_k-x^{\ast }))(x_k-x^{\ast })\mathrm{d}t.\end{array}$$假定$J$在$x^{\ast }$附近有Lipschitz连续性, 则$$\begin{array}{rl}\Vert x_k+p_k^{\mathrm{G}\mathrm{N}}-x^{\ast }\Vert & \le {\int }_0^1\Vert [J^{T}J(x_k){]}^{-1}H(x^{\ast }+t(x_k-x^{\ast }))\Vert \Vert x_k-x^{\ast }\Vert \mathrm{d}t+O(\Vert x_k-x^{\ast }{\Vert }^2)\\ & \approx \Vert [J^{T}J(x^{\ast }){]}^{-1}H(x^{\ast })\Vert \Vert x_k-x^{\ast }\Vert +O(\Vert x_k-x^{\ast }{\Vert }^2).\end{array}$$因此, 若$\Vert [J^{T}J(x^{\ast }){]}^{-1}H(x^{\ast })\Vert \ll 1$, 单位步长的Gauss-Newton步的效果就很好, 从而有较好的收敛性. 特别当$H(x^{\ast })=O$(当残差为线性时), 就有二次收敛性.
而当$n,m$都很大且Jacobi矩阵$J(x)$稀疏时, 每步迭代通过分解$J_k$或$J_k^T{J}_k$精确计算步长的代价(相较于计算函数和梯度值)就会非常大. 基于此, 我们可以构造类似于第七章中非精确牛顿法的Gauss-Newton法的不精确变体. 在这些方法中, 我们直接以$J_k^T{J}_k$代替Hessian${\nabla }^{2}f(x_k)$. 与之前相同, 这一半正定近似简化了算法的许多方面.

3.2 Levenberg-Marquardt法

3.2.1 Levenberg-Marquardt法介绍

3.1中介绍的Gauss-Newton法其实就是线搜索框架下的牛顿法. 唯一的区别在于, 对Hessian我们充分挖掘了问题的内在结构, 使用了更加便利与高效的近似方式. Levenberg-Marquardt法也可用同样的Hessian近似得到, 不同在于它嵌入的是信赖域的框架. 信赖域的使用避免了Gauss-Newton法的一个缺陷, 即当Jacobi矩阵$J(x)$(接近)亏秩时往往效果不好. 由于二者使用相同的Hessian近似, 因此它们的收敛性质也是相似的.
Levenberg-Marquardt法可用第四章信赖域的框架阐明与分析. (事实上, Levenberg-Marquardt法有时也被视为是一般无约束优化信赖域算法的前身.) 我们选取球形的信赖域, 此时每步迭代的子问题为$$\underset{p}{\min }\frac{1}{2}\Vert J_{k}p+r_k{\Vert }^2,\mathrm{s}.\mathrm{t}.\Vert p\Vert \le {\Delta }_k,$$其中${\Delta }_k>0$为信赖域半径. 事实上, 我们选取的模型函数为$$m_k(p)=\frac{1}{2}\Vert r_k{\Vert }^2+p^T{J}_k^T{r}_k+\frac{1}{2}{p}^T{J}_k^T{J}_{k}p.$$下面的讨论中我们省去迭代指标$k$. 第四章的结论让我们对以上子问题的解有了如下的了解: 当Gauss-Newton法的$p^{\mathrm{G}\mathrm{N}}$严格落在信赖域中(即$\Vert p^{\mathrm{G}\mathrm{N}}\Vert <\Delta$)时, 此步$p^{\mathrm{G}\mathrm{N}}$也是子问题的解; 否则, 存在$\lambda >0$使得解$p=p^{\mathrm{L}\mathrm{M}}$满足$\Vert p\Vert =\Delta$以及$$(J^{T}J+\lambda I)p=-J^{T}r.$$注意$J^{T}J$本身半正定以及$\lambda \ge 0$保证了第四章中结论中的半正定性. 这就是下面的引理.

引理2 $p^{\mathrm{L}\mathrm{M}}$为信赖域子问题$$\underset{p}{\min }\Vert Jp+r{\Vert }^2,\mathrm{s}.\mathrm{t}.\Vert p\Vert \le \Delta$$的解当且仅当$p^{\mathrm{L}\mathrm{M}}$可行且存在标量$\lambda \ge 0$使得$$\begin{array}{rl}(J^{T}J+\lambda I)p^{\mathrm{L}\mathrm{M}}& =-J^{T}r,\\ \lambda (\Delta -\Vert p^{\mathrm{L}\mathrm{M}}\Vert )& =0.\end{array}$$

求解方程$(J^{T}J+\lambda I)p=-J^{T}r$实际上等价于求解以下线性最小二乘问题$$\underset{p}{\min }\frac{1}{2}{\Vert \left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right]p+\left[\begin{array}{c}r\\ O\end{array}\right]\Vert }^2.$$如同Gauss-Newton法中所说明, 这一等价性使我们不计算矩阵-矩阵乘积$J^{T}J$以及其Cholesky分解, 就可求解子问题.

3.2.2 Levenberg-Marquardt法的实施

1. 关于Cholesky分解.
   为求得引理2中的$\lambda$, 我们可以使用第四章中的求根算法. 这一过程是良好的: 只要当前的估计${\lambda }^{(l)}$为正, Cholesky因子$R$就一定存在. 由$B=J^{T}J$的特殊结构, 我们无需每步重新计算$B+\lambda I$的Cholesky分解.

   我们先来关注如何高效地求得系数矩阵$$\left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right]$$的QR分解:$$\left[\begin{array}{c}{R}_{\lambda }\\ O\end{array}\right]=Q_{\lambda }^T\left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right],$$其中$Q_{\lambda }$正交, $R_{\lambda }$上三角. 易知, $R_{\lambda }$即满足$R_{\lambda }^T{R}_{\lambda }=(J^{T}J+\lambda I)$.

   我们可以组合使用Householder变换和Givens变换以节省QR分解的计算时间. 假定我们使用Householder变换单独计算了$J$的QR分解$$J=Q\left[\begin{array}{c}R\\ O\end{array}\right].$$于是我们有$$\left[\begin{array}{c}R\\ O\\ \sqrt{\lambda }I\end{array}\right]=\left[\begin{array}{cc}{Q}^T& \\ & I\end{array}\right]\left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right].$$上式左端的矩阵的上半部分为上三角矩阵, 下半部分则包括$n$个非零项. 因此左端矩阵可用$n(n+1)/2$次Givens变换化为上三角阵(这里的计数包括了消除旋转过程中产生填充的过程). 具体说来, 头几步为:

   1. 旋转$R$的第$n$行与$\sqrt{\lambda }I$的第$n$行, 消去$\sqrt{\lambda }I$的$(n,n)$元;
   2. 旋转$R$的第$n-1$行与$\sqrt{\lambda }I$的第$n-1$行, 消去$\sqrt{\lambda }I$的$(n-1,n-1)$元. 这一步旋转会产生$\sqrt{\lambda }I$的$(n-1,n)$位置上的填充, 而这可通过旋转$R$的第$n$行与$\sqrt{\lambda }I$的第$n-1$行消去;
   3. 旋转$R$的第$n-2$行与$\sqrt{\lambda }I$的第$n-2$行, 消去$\sqrt{\lambda }I$的$(n-2,n-2)$元. 这一步旋转会产生$\sqrt{\lambda }I$的$(n-2,n-1),(n-2,n)$位置上的填充, 而这可通过先后旋转$R$的第$n-1$行与$\sqrt{\lambda }I$的第$n-2$行、$R$的第$n$行与$\sqrt{\lambda }I$的第$n-2$行消去.
   4. …

   依此类推. 若将所有的Givens变换汇成一个矩阵${\bar{Q}}_{\lambda }$, 我们就有$${\bar{Q}}_{\lambda }^T\left[\begin{array}{c}R\\ O\\ \sqrt{\lambda }I\end{array}\right]=\left[\begin{array}{c}{R}_{\lambda }\\ O\\ O\end{array}\right],$$因此前面的正交矩阵$Q_{\lambda }$就是$$Q_{\lambda }=\left[\begin{array}{cc}Q& \\ & I\end{array}\right]{\bar{Q}}_{\lambda }.$$这一方法的优点在于, 当我们在求根时会改变$\lambda$的值, 而这样我们就只需要再计算${\bar{Q}}_{\lambda }$而无需再管Householder变换的部分. 这在$m\gg n$时可以节省很多的计算量: 对$\lambda$计算${\bar{Q}}_{\lambda }$与$R_{\lambda }$仅需$O(n^3)$次运算, 而计算$Q$则需$O(m{n}^2)$次运算.

2. 尺度变换.
   最小二乘问题往往尺度较为恶性, 比如一些变量可能会达$10^4$量阶, 而其他的一些又会小到$10^{-6}$. 若我们忽略如此巨大的差距, 算法就会不稳定或者产生一些不好的解. 一种减缓尺度带来的问题的途径是, 选取适当的椭球型信赖域代替上述的球形信赖域. 此时信赖域子问题变为:$$\underset{p}{\min }\frac{1}{2}\Vert J_{k}p+r_k{\Vert }^2,\mathrm{s}.\mathrm{t}.\Vert D_{k}p\Vert \le {\Delta }_k,$$其中$D_k$为对角元为正的对角阵. 相应地, 解满足$$(J_k^T{J}_k+\lambda D_k^2)p_k^{\mathrm{L}\mathrm{M}}=-J_k^T{r}_k,$$这又等价于求解下面的线性最小二乘问题$$\underset{p}{\min }{\Vert \left[\begin{array}{c}{J}_k\\ \sqrt{\lambda }{D}_k\end{array}\right]p+\left[\begin{array}{c}{r}_k\\ O\end{array}\right]\Vert }^2.$$这里对角阵$D_k$可随迭代改变, 其依据为$x$的每个分量的典型范围信息. 若变动在一定范围内, 则球形情形的收敛理论就仍然适用, 其中仅需稍微做些修正. 进一步地, 以上计算$R_{\lambda }$的步骤无需改动. Seber与Wild表示可以选取$D_k^2$为$J_k^T{J}_k$的对角元, 从而使得算法在$x$的对角尺度变换下不变. 这与第四章中缩放Hessian对角元的方法类似.

3. 大规模问题.
   而对于$m,n$都较大以及$J(x)$稀疏的问题, 我们更倾向于使用第七章CG-Steihaug算法求解, 其中以$J_k^T{J}_k$代替真实的${\nabla }^2{f}_k$. $J_k^T{J}_k$的半正定性可用来简化算法, 这是因为原本算法中着重考虑的负曲率不会出现. 同时我们也不需要显式地去计算$J_k^T{J}_k$, 而是先后做两次矩阵-向量乘积.

3.2.3 Levenberg-Marquardt法的收敛性

为达全局收敛, 我们其实不必精确求解信赖域子问题. 下面的收敛性结果为第四章中定理的直接推论.

定理3 设信赖域算法中$\eta \in (0,\frac{1}{4})$, 水平集$\mathcal{L}$有界, 残差函数$r_j(\cdot ),j=1,2,\dots ,m$在$\mathcal{L}$的一个邻域$\mathcal{N}$中Lipschitz连续可微. 假设对每个$k$, 近似解$p_k$满足不等式$$m_k(0)-m_k(p_k)\ge c_1\Vert J_k^T{r}_k\Vert \min \left({\Delta }_k,\frac{\Vert J_k^T{r}_k\Vert }{\Vert J_k^T{J}_k\Vert }\right),$$其中$c_1>0,\Vert p_k\Vert \le \gamma {\Delta }_k,\gamma \ge 1$. 于是有$$\underset{k\to \infty }{\lim }\nabla f_k=\underset{k\to \infty }{\lim }{J}_k^T{r}_k=0.$$

也如第四章, 我们不需精确地计算上面不等式右端项, 而仅需要求近似解$p_k$给出的函数值下降不低于Cauchy点. 而Cauchy点可用第四章的方法方便地计算. 若使用迭代算法CG-Steihaug, 则不等式对$c_1=1/2$自动成立, 这是因为CG_Steihaug的$p_k$第一步估计就是Cauchy点, 而后面的估计只可能会给出更小的函数值.

Levenberg-Marquardt法的局部收敛性质与Gauss-Newton法类似. 在解$x^{\ast }$附近${\nabla }^{2}f(x^{\ast })$的第一项起主要作用, 此时信赖域约束不起作用, 算法将取Gauss-Newton步从而有较快的收敛速度.

3.3 大残差问题的算法

对于大残差的问题, 我们就不能再忽略${\nabla }^{2}f(x)$的第二项了. 在数据拟合问题中, 大残差的出现可能就说明模型不适合数据或者是在观测时引入了较大的误差. 尽管如此, 我们仍需要利用当前的模型和数据求解最小二乘问题, 以提出在观测的赋权、模型建立或者数据收集过程中可以做出的改进.

在大残差问题中, Gauss-Newton法与Levenberg-Marquardt法的渐进收敛速度仅为线性——这要比一些求解一般无约束问题的算法(如牛顿法、拟牛顿法)慢. 若每个Hessian阵${\nabla }^2{r}_j$容易计算, 我们不如忽略最小二乘而直接使用信赖域或线搜索框架下的牛顿法计算. 无需计算${\nabla }^2{r}_j$的拟牛顿法也是个选择. 不过话说回来, 牛顿法与拟牛顿法在迭代早期(即还未进入解的某个邻域)的表现可能并不如Gauss-Newton法与Levenberg-Marquardt法.

当然通常我们是没有问题是小的还是大残差的先验的. 因此, 使用混合算法就比较合理了. 具体说, 它们在残差较小时表现得像Gauss-Newton法或Levenberg-Marquardt法(从而也继承了相应的计算优势), 但在残差较大时转为牛顿法或拟牛顿法.

我们有很多构建混合算法的方式. 由Fletcher和Xu提出的一种方式需要保存一系列正定Hessian近似$B_k$:

• 若由$x_k$出发的Gauss-Newton步以一定的因子(如5)减小了函数值, 我们就采纳这一步并重写$B_k$为$J_k^T{J}_k$.
• 否则, 使用$B_k$计算搜索方向, 并利用线搜索得到新点$x_{k+1}$.

二者均以类似于BFGS更新公式的方式更新$B_k$得到$B_{k+1}$. 在零残差的情形, 这一策略最终总会采纳Gauss-Newton步(从而二次收敛); 在残差非零情形, 最终会约减为BFGS(从而超线性收敛).

第二种结合Gauss-Newton法与拟牛顿法的方式是仅保存Hessian二阶部分的近似. 也就是说, 我们保留矩阵列 { S k } \{S_k\} {Sk​}近似${\sum }_{j=1}^m{r}_j(x_k){\nabla }^2{r}_j(x_k)$, 然后以之估计整个Hessian$$B_k=J_k^T{J}_k+S_k,$$接着使用信赖域或线搜索计算$p_k$. 对$S_k$的更新要求近似矩阵$B_k$或其组成部分能较好地模拟刚走完的那一步的特征. 更新公式则基于割线方程. 这里有许多不同的方式可以用来定义割线方程以及设计其他的条件得到$S_k$的更新公式. 下面我们介绍Dennis, Gay和Welsch的算法.
理想状态下, $S_{k+1}$应当很接近于$x=x_{k+1}$处的二阶项:$$S_{k+1}\approx \sum _{j=1}^m{r}_j(x_{k+1})]{\nabla }^2{r}_j(x_{k+1}).$$我们不想计算右端的${\nabla }^2{r}_j$, 因此我们可代之以某个近似$(B_j{)}_{k+1}$并在$(B_j{)}_{k+1}$上提些条件. 也即$$\begin{array}{rl}(B_j{)}_{k+1}(x_{k+1}-x_k)& =\nabla r_j(x_{k+1})-\nabla r_j(x_k)\\ & =(\mathrm{r}\mathrm{o}\mathrm{w}j\mathrm{o}\mathrm{f}J(x_{k+1}){)}^T-(\mathrm{r}\mathrm{o}\mathrm{w}j\mathrm{o}\mathrm{f}J(x_k){)}^T.\end{array}$$这一条件最终推出$S_{k+1}$上的割线方程:$$\begin{array}{rl}{S}_{k+1}(x_{k+1}-x_k)& =\sum _{j=1}^m{r}_j(x_{k+1})(B_j{)}_{k+1}(x_{k+1}-x_k)\\ & =\sum _{j=1}^m{r}_j(x_{k+1})[(\mathrm{r}\mathrm{o}\mathrm{w}j\mathrm{o}\mathrm{f}J(x_{k+1}){)}^T-(\mathrm{r}\mathrm{o}\mathrm{w}j\mathrm{o}\mathrm{f}J(x_k){)}^T]\\ & =J_{k+1}^T{r}_{k+1}-J_k^T{r}_{k+1}.\end{array}$$当然割线方程还不能完全决定$S_{k+1}$. Dennis, Gay和Welsch又要求$S_{k+1}$对称并且$S_{k+1}-S_k$要在某种意义下达到极小, 从而得到了下面的更新公式:$$S_{k+1}=S_k+\frac{(y^{\#}-S_{k}s)y^T+y(y^{\#}-S_{k}s{)}^T}{y^{T}s}-\frac{(y^{\#}-S_{k}s{)}^{T}s}{(y^{T}s{)}^2}y{y}^T,$$其中$$\begin{array}{rl}s& =x_{k+1}-x_k,\\ y& =J_{k+1}^T{r}_{k+1}-J_k^T{r}_k,\\ y^{\#}& =J_{k+1}^T{r}_{k+1}-J_k^T{r}_{k+1}.\end{array}$$注意上述更新公式仅是DFP更新公式的微小改动版本. 若$y^{\#},y$相同, 则二者就完全一样了. Dennis, Gay和Welsch将他们的近似Hessian$J_k^T{J}_k+S_k$联合信赖域的框架使用, 不过其中需要加以更多的约束以提升表现. 其基本策略的一个缺陷在于, $S_k$的更新策略并不保证当迭代点趋近一个零残差解时$S_k$会消失, 因此有时难以保证超线性收敛性. 这一问题可通过在$S_k$更新之前对其缩放避免: 我们以${\tau }_k{S}_k$代替$S_k$, 其中$${\tau }_k=\min \left(1,\frac{|s^T{y}^{\#}|}{|s^T{S}_{k}s|}\right).$$最后, 若Gauss-Newton法能产生较好的下降, 我们就应省去Hessian近似中的$S_k$.

4. 正交距离回归

在例1中我们假设在抽取血液样本的时间上不存在误差, 从而模型$\varphi (x;t_j)$和观测$y_j$的差在于模型构造或$y_j$的测量误差上. 这里实际上我们假设在横坐标——时间$t_j$——上的误差远比观测上误差小从而可以忽略. 这一假设通常是合理的, 但有时若我们不考虑横坐标上的误差, 得到的结果就会有严重偏差. 将这种误差纳入考量的模型在统计上称作变量含误差模型(errors-in-variables models), 而引出的优化问题则称为总体最小二乘(线性模型下)或正交距离回归(非线性模型下).

下面我们从数学上严格地进行表述. 引入$t_j$上的扰动${\delta }_j$, $y_j$上的扰动${ϵ}_j$. 我们需要求得这$2m$个扰动的值, 使之极小化模型和观测的差异. 这里差异以加权最小二乘目标函数度量. 具体地说, 定义极小化问题为$$\underset{x,{\delta }_j,{ϵ}_j}{\min }\frac{1}{2}\sum _{j=1}^m{w}_j^2{ϵ}_j^2+d_j^2{\delta }_j^2,\mathrm{s}.\mathrm{t}.y_j=\varphi (x;t_j+{\delta }_j)+{ϵ}_j,j=1,2,\dots ,m.$$$w_i,d_i$为权重, 它们由使用者或一些自动估计误差项的相对重要性的方式决定.
当我们图示上述问题时, 我们便能知道"正价距离回归"这个称谓是怎么来的了.

若所有的$w_i,d_i$都相等, 则目标函数中求和的每一项就是点$(t_j,y_j)$与曲线$\varphi (x;t)$之间的最短距离. 而点与曲线之间的最短直线就与曲线在与直线的交点处正交.
我们可以用约束消去${ϵ}_j$, 从而得到无约束最小二乘问题$$\underset{x,\delta }{\min }F(x,\delta )=\frac{1}{2}\sum _{j=1}^m{w}_j^2[y_j-\varphi (x;t_j+{\delta }_j){]}^2+d_j^2{\delta }_j^2=\frac{1}{2}\sum _{j=1}^{2m}{r}_j^2(x,\delta ),$$其中$\delta =({\delta }_1,{\delta }_2,\dots ,{\delta }_m{)}^T$, $$r_j(x,\delta )=\{\begin{array}{ll}{w}_j[\varphi (x;t_j+{\delta }_j)-y_j],& j=1,2,\dots ,m,\\ d_{j-m}{\delta }_{j-m},& j=m+1,\dots ,2m.\end{array}$$此时的问题就是带有$2m$个残差和$m+n$个未知量的标准最小二乘问题. 我们可以使用本章中介绍的算法求解之. 不过直接使用可能会带来计算量上的困难, 这是因为此时的问题所带的参数数目与观测数可能原始的问题要大得多.

但若我们进一步深究Gauss-Newton法或Levenberg-Marquardt法中Jacobi矩阵的结构, 我们会发现: 它有许多元素都为0, 例如$$\begin{gathered} \frac{\partial r_j}{\partial {\delta }_i}=\frac{\partial [\varphi (t_j+{\delta }_j;x)-y_j]}{\partial {\delta }_i}=0,i,j=1,2,\dots ,m,i̸=j, \\ \frac{\partial r_j}{\partial x_i}=0,j=m+1,\dots ,2m,i=1,2,\dots ,n, \\ \frac{\partial r_{m+j}}{\partial {\delta }_i}=\{\begin{array}{ll}{d}_j& i=j,\\ 0& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\end{array}i,j=1,2,\dots ,m. \end{gathered}$$因此我们可以将Jacobi矩阵分块写作$$J(x,\delta )=\left[\begin{array}{cc}\hat{J}& V\\ O& D\end{array}\right],$$其中$V,D$为$m\times m$对角阵, $\hat{J}$为$m\times n$矩阵, 其中元素为$w_j\varphi (t_j+{\delta }_j;x)$对$x$的偏导. 对应地, 将$p,r$也分块为$$p=\left[\begin{array}{c}{p}_x\\ p_{\delta }\end{array}\right],r=\left[\begin{array}{c}{\hat{r}}_1\\ {\hat{r}}_2\end{array}\right],$$并将正规方程分块为$$\left[\begin{array}{cc}{\hat{J}}^T\hat{J}& {\hat{J}}^{T}V\\ V\hat{J}& V^2+D^2+\lambda I\end{array}\right]\left[\begin{array}{c}{p}_x\\ p_{\delta }\end{array}\right]=-\left[\begin{array}{c}{\hat{J}}^T{\hat{r}}_1\\ V{\hat{r}}_1+D{\hat{r}}_2\end{array}\right].$$由于右下方的子阵$V^2+D^2+\lambda I$对角, 因此我们可以方便地消去$p_{\delta }$得到用来求$p_x$的$n\times n$子系统. 这样求得一步的计算量就仅比标准最小二乘模型下的$m\times n$问题稍微大一点.

5. 再谈大规模情形

对大规模非线性最小二乘问题, Wright和Holt提出了一种非精确的Levenberg-Marquardt方法. 这一方法直接控制参数$\lambda$的变动而非借用与信赖域算法的联系. 它将(类似第七章)采纳满足以下不等式的${\bar{p}}_k$:$$\Vert (J_k^T{J}_k+{\lambda }_{k}I){\bar{p}}_k+J_k^T{r}_k\Vert \le {\eta }_k\Vert J_k^T{r}_k\Vert ,{\eta }_k\in [0,\eta ],$$其中$\eta \in (0,1)$为一常数, { η k } \{\eta_k\} {ηk​}为一强迫序列. 之后再由实际对预测的下降比值决定是否要采纳${\bar{p}}_k$. 在一定的假设条件下, 我们可以证明这一方法的全局收敛性.