第十章: 最小二乘问题

文章目录

第十章: 最小二乘问题
- 1. 背景介绍
- 2. 线性最小二乘问题及求解算法
- - 2.1 基于Cholesky分解的直接法
  - 2.2 基于QR分解的直接法
  - 2.3 基于奇异值分解(SVD)的直接法
  - 2.4 三种方法的讨论
- 3. 求解非线性最小二乘问题的算法
- - 3.1 Gauss-Newton法
  - - 3.1.1 Gauss-Newton法介绍
    - 3.1.2 Gauss-Newton法的收敛性
  - 3.2 Levenberg-Marquardt法
  - - 3.2.1 Levenberg-Marquardt法介绍
    - 3.2.2 Levenberg-Marquardt法的实施
    - 3.2.3 Levenberg-Marquardt法的收敛性
  - 3.3 大残差问题的算法
- 4. 正交距离回归
- 5. 再谈大规模情形

在 最小二乘问题(Least-square problems)中, 目标函数 fff具有如下特殊形式: f(x)=12∑j=1mrj2(x),f(x)=\frac{1}{2}\sum_{j=1}^mr_j^2(x),f(x)=21j=1∑mrj2(x),其中每个 rj:Rn→Rr_j:\mathbb{R}^n\to\mathbb{R}rj:Rn→R为光滑函数. 我们称每个 rjr_jrj为 残差(residual). 本章中假设 m≥nm\ge nm≥n.

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

上述优化问题比之一般的无约束极小化问题, 具有最优值非负的特点. 其次, 上述fff的特殊结构能使得最小二乘问题比一般的无约束极小化问题更易求解. 首先将每个残差组分rjr_jrj组装成一个残差向量r:Rn→Rmr:\mathbb{R}^n\to\mathbb{R}^mr:Rn→Rm:r(x)=(r1(x),r2(x),…,rm(x))T.r(x)=(r_1(x),r_2(x),\ldots,r_m(x))^T.r(x)=(r1(x),r2(x),…,rm(x))T.利用这一表示, 我们可以将fff写作f(x)=12∥r(x)∥22f(x)=\frac{1}{2}\Vert r(x)\Vert_2^2f(x)=21∥r(x)∥22. fff的导数就可以m×nm\times nm×nJacobi矩阵J(x)J(x)J(x)的形式表示:J(x)=[∂ri∂xj]i=1,2,…,m,j=1,2,…,n=[∇r1(x)T∇r2(x)T⋮∇rm(x)T],J(x)=\left[\frac{\partial r_i}{\partial x_j}\right]_{i=1,2,\ldots,m,j=1,2,\ldots,n}=\begin{bmatrix}\nabla r_1(x)^T\\\nabla r_2(x)^T\\\vdots\\\nabla r_m(x)^T\end{bmatrix},J(x)=[∂xj∂ri]i=1,2,…,m,j=1,2,…,n=⎣⎢⎢⎢⎡∇r1(x)T∇r2(x)T⋮∇rm(x)T⎦⎥⎥⎥⎤,其中每个∇rj(x),j=1,2,…,m\nabla r_j(x),j=1,2,\ldots,m∇rj(x),j=1,2,…,m为rjr_jrj的梯度. 于是fff的梯度与Hessian矩阵为:∇f(x)=∑j=1mrj(x)∇rj(x)=J(x)Tr(x),∇2f(x)=∑j=1m∇rj(x)∇rj(x)T+∑j=1mrj(x)∇2rj(x)=J(x)TJ(x)+∑j=1mrj(x)∇2rj(x).\begin{aligned}\nabla f(x)&=\sum_{j=1}^mr_j(x)\nabla r_j(x)=J(x)^Tr(x),\\\nabla^2f(x)&=\sum_{j=1}^m\nabla r_j(x)\nabla r_j(x)^T+\sum_{j=1}^mr_j(x)\nabla^2r_j(x)\\&=J(x)^TJ(x)+\sum_{j=1}^mr_j(x)\nabla^2r_j(x).\end{aligned}∇f(x)∇2f(x)=j=1∑mrj(x)∇rj(x)=J(x)Tr(x),=j=1∑m∇rj(x)∇rj(x)T+j=1∑mrj(x)∇2rj(x)=J(x)TJ(x)+j=1∑mrj(x)∇2rj(x).一般, 残差的一阶偏导(从而Jacobi矩阵J(x)J(x)J(x))相对容易(或相对便宜)计算. 从而∇f(x)\nabla f(x)∇f(x)的表达式是可使用的. 而有了J(x)J(x)J(x), 我们也可以计算∇2f(x)\nabla^2f(x)∇2f(x)中的J(x)TJ(x)J(x)^TJ(x)J(x)TJ(x). 这一过程完全不需要计算rjr_jrj的任何二阶导数. ∇2f(x)\nabla^2f(x)∇2f(x)的这一部分"免费"可用性是最小二乘问题的关键特征. 进一步地, J(x)TJ(x)J(x)^TJ(x)J(x)TJ(x)往往要比第二项来得更重要. 这可能是由于在解附近

残差rjr_jrj接近于线性(也就是说, ∇2rj(x)\nabla^2r_j(x)∇2rj(x)相对较小); 或者
残差较小(即rj(x)r_j(x)rj(x)相对较小).

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

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

1. 背景介绍

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

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

基于我们过往的经验, 我们发现函数ϕ(x;t)\phi(x;t)ϕ(x;t)对药物在ttt时刻的凝聚度有很好的预测能力, 其中xxx为五维参数向量x=(x1,x2,x3,x4,x5)x=(x_1,x_2,x_3,x_4,x_5)x=(x1,x2,x3,x4,x5):ϕ(x;t)=x1+tx2+t2x3+x4e−x5t.\phi(x;t)=x_1+tx_2+t^2x_3+x_4e^{-x_5t}.ϕ(x;t)=x1+tx2+t2x3+x4e−x5t.我们待定xxx, 最终需要我们的模型以某种方式最佳匹配拟合我们的观测数据. 一种较好的表示预测模型值域观测值差距的方式就是如下最小二乘函数:12∑j=1m[ϕ(x;tj)−yj]2.\frac{1}{2}\sum_{j=1}^m[\phi(x;t_j)-y_j]^2.21j=1∑m[ϕ(x;tj)−yj]2.我们定义rj(x)=ϕ(x;tj)−yjr_j(x)=\phi(x;t_j)-y_jrj(x)=ϕ(x;tj)−yj. 几何上, 每个∣rj∣|r_j|∣rj∣均表示点(tj,yj)(t_j,y_j)(tj,yj)与曲线ϕ(x;t)\phi(x;t)ϕ(x;t)(视作ttt的函数, xxx为固定的参数向量)的垂直距离. 可见下图.

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

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

最大绝对值. max⁡j=1,2,…,m∣ϕ(x;tj)−yj∣.\max_{j=1,2,\ldots,m}|\phi(x;t_j)-y_j|.j=1,2,…,mmax∣ϕ(x;tj)−yj∣.
绝对值和. ∑j=1m∣ϕ(x;tj)−yj∣.\sum_{j=1}^m|\phi(x;t_j)-y_j|.j=1∑m∣ϕ(x;tj)−yj∣.

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

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

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

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

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

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

2.1 基于Cholesky分解的直接法

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

计算系数矩阵JTJJ^TJJTJ以及右端项JTyJ^TyJTy;
计算对称矩阵JTJJ^TJJTJ的Cholesky分解;
经前代、回代求得解x∗x^*x∗.

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

2.2 基于QR分解的直接法

第二种方法基于对矩阵JJJ的QR分解. 由于正交变换保长, 所以∥Jx−y∥=∥QT(Jx−y)∥,\Vert Jx-y\Vert=\Vert Q^T(Jx-y)\Vert,∥Jx−y∥=∥QT(Jx−y)∥,其中QQQ为任意m×mm\times mm×m的正交阵. 假设我们对JJJ有带列选主元的QR分解:JΠ=Q[RO]=[Q1Q2][RO]=Q1R,J\Pi=Q\begin{bmatrix}R\\O\end{bmatrix}=\begin{bmatrix}Q_1 & Q_2\end{bmatrix}\begin{bmatrix}R\\ O\end{bmatrix}=Q_1R,JΠ=Q[RO]=[Q1Q2][RO]=Q1R,其中

Π\PiΠ为n×nn\times nn×n的排列矩阵(因此正交);
QQQ为m×mm\times mm×m的正交阵;
Q1Q_1Q1为QQQ的前nnn列, Q2Q_2Q2为QQQ的后m−nm-nm−n列;
RRR为对角元为正的n×nn\times nn×n上三角阵;

因此我们有∥Jx−y∥22=∥[Q1TQ2T](JΠΠTx−y)∥22=∥[RO](ΠTx)−[Q1TyQ2Ty]∥2=∥R(ΠTx)−Q1Ty∥22+∥Q2Ty∥2.\begin{aligned}\Vert Jx-y\Vert_2^2&=\left\Vert\begin{bmatrix}Q_1^T\\Q_2^T\end{bmatrix}(J\Pi\Pi^Tx-y)\right\Vert_2^2\\&=\left\Vert\begin{bmatrix}R\\O\end{bmatrix}(\Pi^Tx)-\begin{bmatrix}Q_1^Ty\\Q_2^Ty\end{bmatrix}\right\Vert^2\\&=\Vert R(\Pi^Tx)-Q_1^Ty\Vert_2^2+\Vert Q_2^Ty\Vert^2.\end{aligned}∥Jx−y∥22=∥∥∥∥[Q1TQ2T](JΠΠTx−y)∥∥∥∥22=∥∥∥∥[RO](ΠTx)−[Q1TyQ2Ty]∥∥∥∥2=∥R(ΠTx)−Q1Ty∥22+∥Q2Ty∥2.对于上面的第二项我们无能为力, 不过我们可以把第一项消成0, 即令x∗=ΠR−1Q1Ty.x^*=\Pi R^{-1}Q_1^Ty.x∗=ΠR−1Q1Ty.(实际操作时, 我们先回代求解Rz=Q1TyRz=Q_1^TyRz=Q1Ty, 再排列zzz的元素得到x∗=Πzx^*=\Pi zx∗=Πz.)

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

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

在一些特殊情形下, 我们要求解对数据(J,y)(J,y)(J,y)上的扰动更加鲁棒. 此时基于JJJ的SVD的方法就可派上用场. 矩阵JJJ的SVD为:J=U[SO]VT=[U1U2][SO]VT=U1SVT,J=U\begin{bmatrix}S\\O\end{bmatrix}V^T=\begin{bmatrix}U_1&U_2\end{bmatrix}\begin{bmatrix}S\\O\end{bmatrix}V^T=U_1SV^T,J=U[SO]VT=[U1U2][SO]VT=U1SVT,其中

UUU为m×mm\times mm×m的正交阵;
U1U_1U1为UUU的前nnn列, U2U_2U2为UUU的后m−nm-nm−n列;
VVV为n×nn\times nn×n的正交阵;
SSS为n×nn\times nn×n对角阵, 其中对角元σ1≥σ2≥⋯≥σn>0\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_n>0σ1≥σ2≥⋯≥σn>0.

(注意JTJ=VS2VTJ^TJ=VS^2V^TJTJ=VS2VT, 因此VVV的列就是JTJJ^TJJTJ的对应于奇异值σj2,j=1,2,…,n\sigma_j^2,j=1,2,\ldots,nσj2,j=1,2,…,n的特征向量.) 类似地我们有∥Jx−y∥2=∥[SO](VTx)−[U1TU2T]y∥2=∥S(VTx)−U1Ty∥2+∥U2Ty∥2.\begin{aligned}\Vert Jx-y\Vert^2&=\left\Vert\begin{bmatrix}S\\O\end{bmatrix}(V^Tx)-\begin{bmatrix}U_1^T\\U_2^T\end{bmatrix}y\right\Vert^2\\&=\Vert S(V^Tx)-U_1^Ty\Vert^2+\Vert U_2^Ty\Vert^2.\end{aligned}∥Jx−y∥2=∥∥∥∥[SO](VTx)−[U1TU2T]y∥∥∥∥2=∥S(VTx)−U1Ty∥2+∥U2Ty∥2.我们令第一项为0, 即有x∗=VS−1U1Ty.x^*=VS^{-1}U_1^Ty.x∗=VS−1U1Ty.进一步地, 以ui∈Rm,vi∈Rnu_i\in\mathbb{R}^m,v_i\in\mathbb{R}^nui∈Rm,vi∈Rn分别表示U,VU,VU,V的第iii列, 我们有x∗=∑i=1nuiTyσivi.x^*=\sum_{i=1}^n\frac{u_i^Ty}{\sigma_i}v_i.x∗=i=1∑nσiuiTyvi.这一公式给予了我们关于x∗x^*x∗敏度的丰富的信息: 当σi\sigma_iσi很小, x∗x^*x∗对于yyy中(以及JJJ中)影响uiTyu_i^TyuiTy的部分就非常敏感. 这些信息在JJJ几近亏秩时尤其有用, 也就是说当σn/σ1≪1\sigma_n/\sigma_1\ll1σn/σ1≪1时. 我们可以适当地采取一些防护措施避免数值上的不稳定. 有时以SVD的高昂计算量换取这些有效信息是值当的.

2.4 三种方法的讨论

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

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

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

3.1 Gauss-Newton法

3.1.1 Gauss-Newton法介绍

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

∇2fk≈JkTJk\nabla^2f_k\approx J_k^TJ_k∇2fk≈JkTJk的近似省去了我们计算每个残差的Hessian∇2rj,j=1,2,…,m\nabla^2r_j,j=1,2,\ldots,m∇2rj,j=1,2,…,m的功夫. 事实上, 若我们在计算梯度∇fk=JkTrk\nabla f_k=J_k^Tr_k∇fk=JkTrk的过程中就已经计算了Jacobi矩阵JkJ_kJk的话, 这一近似是根本不会牵涉到任何的导数计算的. 这在某些应用上可以节省大量的计算时间.
实际上我们经常看到第一项JTJJ^TJJTJ(相对于第二项)占主的场景(当然这得离解x∗x^*x∗足够近), 从而近似是得当的, Gauss-Newton法的收敛速度也并不会逊于Newton法太多. 具体说, 比如第二项中每一小项的范数(即∣rj(x)∣∥∇2rj(x)∥|r_j(x)|\Vert\nabla^2r_j(x)\Vert∣rj(x)∣∥∇2rj(x)∥)比JTJJ^TJJTJ的特征值要小得多. 在之前我们提到一种情形: 当残差rjr_jrj较小或者它们近似于线性时. 实际中, 许多最小二乘问题在解附近均有较小的残差, 从而保证了Gauss-Newton法的收敛速度.
只要JkJ_kJk列满秩以及梯度∇fk\nabla f_k∇fk非零, 得到的方向pkGNp_k^{\mathrm{GN}}pkGN就是下降方向, 从而可用于线搜索中. 事实上, (pkGN)T∇fk=(pkGN)TJkTrk=−(pkGN)TJkTJkpkGN=−∥JkpkGN∥2≤0.(p_k^{\mathrm{GN}})^T\nabla f_k=(p_k^{\mathrm{GN}})^TJ_k^Tr_k=-(p_k^{\mathrm{GN}})^TJ_k^TJ_kp_k^{\mathrm{GN}}=-\Vert J_kp_k^{\mathrm{GN}}\Vert^2\le0.(pkGN)T∇fk=(pkGN)TJkTrk=−(pkGN)TJkTJkpkGN=−∥JkpkGN∥2≤0.其中最后一个不等式只有在JkpkGN=0J_kp_k^{\mathrm{GN}}=0JkpkGN=0的时候取等, 此时也应当有JkTrk=∇fk=0J_k^Tr_k=\nabla f_k=0JkTrk=∇fk=0. 这就是说xkx_kxk已经是个稳定点了.
Gauss-Newton法所解的方程与线性情形的正规方程方程有一定的相似度. 具体说来, pkGNp_k^{\mathrm{GN}}pkGN可以看做是(也实际上就是)以下线性最小二乘问题的解:min⁡p12∥Jkp+rk∥2.\min_p\frac{1}{2}\Vert J_kp+r_k\Vert^2.pmin21∥Jkp+rk∥2.因此, 我们可以用求解线性最小二乘问题的算法求解以上子问题. 进一步, 若我们使用基于QR分解或SVD的算法, 我们甚至不需要显式地计算出Hessian的近似JkTJkJ_k^TJ_kJkTJk. 若使用共轭梯度法也是一样: 我们只需计算矩阵JkTJkJ_k^TJ_kJkTJk-向量乘积, 而这一步可以通过先后对Jk,JkTJ_k,J_k^TJk,JkT操作得到.
大规模情形. 若残差的数量mmm很大而变量数nnn相对较小, 此时显式地存储JJJ似乎就显得不太符合情理. 不过我们可以通过连续计算rj,∇rj,j=1,2,…,mr_j,\nabla r_j,j=1,2,\ldots,mrj,∇rj,j=1,2,…,m再求和得到:JTJ=∑j=1m(∇rj)(∇rj)T,JTr=∑j=1mrj(∇rj).J^TJ=\sum_{j=1}^m(\nabla r_j)(\nabla r_j)^T,\quad J^Tr=\sum_{j=1}^mr_j(\nabla r_j).JTJ=j=1∑m(∇rj)(∇rj)T,JTr=j=1∑mrj(∇rj).
以上子问题还启发我们给出得到Gauss-Newton法的另外一种途径. 我们利用Taylor展开得到向量值函数的近似r(xk+p)≈rk+Jkpr(x_k+p)\approx r_k+J_kpr(xk+p)≈rk+Jkp. 因此f(xk+p)=12∥r(xk+p)∥2≈12∥Jkp+rk∥2,f(x_k+p)=\frac{1}{2}\Vert r(x_k+p)\Vert^2\approx\frac{1}{2}\Vert J_kp+r_k\Vert^2,f(xk+p)=21∥r(xk+p)∥2≈21∥Jkp+rk∥2,再选取pkGNp_k^{\mathrm{GN}}pkGN作为这一近似模型的极小点.

Gauss-Newton法的实施通常需要沿着pkGNp_k^{\mathrm{GN}}pkGN做线搜索, 这其中需要步长参数αk\alpha_kαk满足第三章中提到的条件, 例如Armijo条件、Wolfe条件.

3.1.2 Gauss-Newton法的收敛性

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

定理1 设每个残差函数rjr_jrj在水平集的一个邻域N\mathcal{N}N内Lipschitz连续可微, Jacobi矩阵J(x)J(x)J(x)在N\mathcal{N}N上满足一致满秩条件. 则由Gauss-Newton法产生的迭代序列{xk}\{x_k\}{xk}(其中步长参数αk\alpha_kαk满足Wolfe条件)成立lim⁡k→∞JkTrk=0.\lim_{k\to\infty}J_k^Tr_k=0.k→∞limJkTrk=0.
证明: 首先注意有界水平集L\mathcal{L}L的邻域可以选取得足够小使得对正常数L,βL,\betaL,β成立以下不等式:∣rj(x)∣≤β,∥∇rj(x)∥≤β,|r_j(x)|\le\beta,\quad\Vert\nabla r_j(x)\Vert\le\beta,∣rj(x)∣≤β,∥∇rj(x)∥≤β,∣rj(x)−rj(x~)∣≤L∥x−x~∥,∥∇rj(x)−∇rj(x~)∥≤L∥x−x~∥,∀x,x~∈N,|r_j(x)-r_j(\tilde{x})|\le L\Vert x-\tilde{x}\Vert,\quad \Vert\nabla r_j(x)-\nabla r_j(\tilde{x})\Vert\le L\Vert x-\tilde{x}\Vert, \quad \forall x,\tilde{x}\in\mathcal{N},∣rj(x)−rj(x~)∣≤L∥x−x~∥,∥∇rj(x)−∇rj(x~)∥≤L∥x−x~∥,∀x,x~∈N,j=1,2,…,m.j=1,2,\ldots,m.j=1,2,…,m.易知存在常数βˉ>0\bar{\beta}>0βˉ>0使得∥J(x)T∥=∥J(x)∥≤βˉ,∀x∈L\Vert J(x)^T\Vert=\Vert J(x)\Vert\le\bar{\beta},\forall x\in\mathcal{L}∥J(x)T∥=∥J(x)∥≤βˉ,∀x∈L, 以及由于∇f(x)=∑j=1mrj(x)∇rj(x)\nabla f(x)=\sum_{j=1}^mr_j(x)\nabla r_j(x)∇f(x)=∑j=1mrj(x)∇rj(x), 从而∇f\nabla f∇f是Lipschitz连续的. 因此Zoutendijk定理的条件满足.
下面我们验证搜索方向pkGNp_k^{\mathrm{GN}}pkGN和负梯度−∇fk-\nabla f_k−∇fk的夹角θk\theta_kθk一致远离π/2\pi/2π/2. 事实上, cos⁡θk=−(∇f)TpGN∥pGN∥∥∇f∥=∥JpGN∥2∥pGN∥∥JTJpGN∥≥γ2∥pGN∥2βˉ2∥pGN∥2=γ2βˉ2>0.\cos\theta_k=-\frac{(\nabla f)^Tp^{\mathrm{GN}}}{\Vert p^{\mathrm{GN}}\Vert\Vert\nabla f\Vert}=\frac{\Vert Jp^{\mathrm{GN}}\Vert^2}{\Vert p^{\mathrm{GN}}\Vert\Vert J^TJp^{\mathrm{GN}}\Vert}\ge\frac{\gamma^2\Vert p^{\mathrm{GN}}\Vert^2}{\bar{\beta}^2\Vert p^{\mathrm{GN}}\Vert^2}=\frac{\gamma^2}{\bar{\beta}^2}>0.cosθk=−∥pGN∥∥∇f∥(∇f)TpGN=∥pGN∥∥JTJpGN∥∥JpGN∥2≥βˉ2∥pGN∥2γ2∥pGN∥2=βˉ2γ2>0.由Zoutendijk定理知∇f(xk)→0\nabla f(x_k)\to 0∇f(xk)→0, 得证.

若JkJ_kJk(对某个kkk)亏秩, 此时一致满秩条件不成立, 系数矩阵JkTJkJ_k^TJ_kJkTJk是奇异的. JkTJkp=−JkTrkJ_k^TJ_kp=-J_k^Tr_kJkTJkp=−JkTrk仍然有解, 不过有无穷多解, 其中每个都具有形式p=∑σi≠0−uiTrkσivi+∑σi=0τivi,p=\sum_{\sigma_i\ne0}-\frac{u_i^Tr_k}{\sigma_i}v_i+\sum_{\sigma_i=0}\tau_iv_i,p=σi̸=0∑−σiuiTrkvi+σi=0∑τivi,这里τi\tau_iτi任意. 但这样一来我们便无法保证cos⁡θk\cos\theta_kcosθk一致远离0, 从而得不到如定理1一般的全局收敛性.

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

3.2 Levenberg-Marquardt法

3.2.1 Levenberg-Marquardt法介绍

3.1中介绍的Gauss-Newton法其实就是线搜索框架下的牛顿法. 唯一的区别在于, 对Hessian我们充分挖掘了问题的内在结构, 使用了更加便利与高效的近似方式. Levenberg-Marquardt法也可用同样的Hessian近似得到, 不同在于它嵌入的是信赖域的框架. 信赖域的使用避免了Gauss-Newton法的一个缺陷, 即当Jacobi矩阵J(x)J(x)J(x)(接近)亏秩时往往效果不好. 由于二者使用相同的Hessian近似, 因此它们的收敛性质也是相似的.
Levenberg-Marquardt法可用第四章信赖域的框架阐明与分析. (事实上, Levenberg-Marquardt法有时也被视为是一般无约束优化信赖域算法的前身.) 我们选取球形的信赖域, 此时每步迭代的子问题为min⁡p12∥Jkp+rk∥2,s.t. ∥p∥≤Δk,\min_p\frac{1}{2}\Vert J_kp+r_k\Vert^2,\quad \mathrm{s.t.\,}\Vert p\Vert\le\Delta_k,pmin21∥Jkp+rk∥2,s.t.∥p∥≤Δk,其中Δk>0\Delta_k>0Δk>0为信赖域半径. 事实上, 我们选取的模型函数为mk(p)=12∥rk∥2+pTJkTrk+12pTJkTJkp.m_k(p)=\frac{1}{2}\Vert r_k\Vert^2+p^TJ_k^Tr_k+\frac{1}{2}p^TJ_k^TJ_kp.mk(p)=21∥rk∥2+pTJkTrk+21pTJkTJkp.下面的讨论中我们省去迭代指标kkk. 第四章的结论让我们对以上子问题的解有了如下的了解: 当Gauss-Newton法的pGNp^{\mathrm{GN}}pGN严格落在信赖域中(即∥pGN∥<Δ\Vert p^{\mathrm{GN}}\Vert<\Delta∥pGN∥<Δ)时, 此步pGNp^{\mathrm{GN}}pGN也是子问题的解; 否则, 存在λ>0\lambda>0λ>0使得解p=pLMp=p^{\mathrm{LM}}p=pLM满足∥p∥=Δ\Vert p\Vert=\Delta∥p∥=Δ以及(JTJ+λI)p=−JTr.(J^TJ+\lambda I)p=-J^Tr.(JTJ+λI)p=−JTr.注意JTJJ^TJJTJ本身半正定以及λ≥0\lambda\ge0λ≥0保证了第四章中结论中的半正定性. 这就是下面的引理.

引理2 pLMp^{\mathrm{LM}}pLM为信赖域子问题min⁡p∥Jp+r∥2,s.t. ∥p∥≤Δ\min_p\Vert Jp+r\Vert^2,\quad\mathrm{s.t.\,}\Vert p\Vert\le\Deltapmin∥Jp+r∥2,s.t.∥p∥≤Δ的解当且仅当pLMp^{\mathrm{LM}}pLM可行且存在标量λ≥0\lambda\ge0λ≥0使得(JTJ+λI)pLM=−JTr,λ(Δ−∥pLM∥)=0.\begin{aligned}(J^TJ+\lambda I)p^{\mathrm{LM}}&=-J^Tr,\\\lambda(\Delta-\Vert p^{\mathrm{LM}}\Vert)&=0.\end{aligned}(JTJ+λI)pLMλ(Δ−∥pLM∥)=−JTr,=0.

求解方程(JTJ+λI)p=−JTr(J^TJ+\lambda I)p=-J^Tr(JTJ+λI)p=−JTr实际上等价于求解以下线性最小二乘问题min⁡p12∥[JλI]p+[rO]∥2.\min_p\frac{1}{2}\left\Vert\begin{bmatrix}J\\\sqrt{\lambda}I\end{bmatrix}p+\begin{bmatrix}r\\O\end{bmatrix}\right\Vert^2.pmin21∥∥∥∥[JλI]p+[rO]∥∥∥∥2.如同Gauss-Newton法中所说明, 这一等价性使我们不计算矩阵-矩阵乘积JTJJ^TJJTJ以及其Cholesky分解, 就可求解子问题.

3.2.2 Levenberg-Marquardt法的实施

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

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

我们可以组合使用Householder变换和Givens变换以节省QR分解的计算时间. 假定我们使用Householder变换单独计算了JJJ的QR分解J=Q[RO].J=Q\begin{bmatrix}R\\O\end{bmatrix}.J=Q[RO].于是我们有[ROλI]=[QTI][JλI].\begin{bmatrix}R\\O\\\sqrt{\lambda}I\end{bmatrix}=\begin{bmatrix}Q^T&\\& I\end{bmatrix}\begin{bmatrix}J\\\sqrt{\lambda}I\end{bmatrix}.⎣⎡ROλI⎦⎤=[QTI][JλI].上式左端的矩阵的上半部分为上三角矩阵, 下半部分则包括nnn个非零项. 因此左端矩阵可用n(n+1)/2n(n+1)/2n(n+1)/2次Givens变换化为上三角阵(这里的计数包括了消除旋转过程中产生填充的过程). 具体说来, 头几步为:
1. 旋转RRR的第nnn行与λI\sqrt{\lambda}IλI的第nnn行, 消去λI\sqrt{\lambda}IλI的(n,n)(n,n)(n,n)元;
2. 旋转RRR的第n−1n-1n−1行与λI\sqrt{\lambda}IλI的第n−1n-1n−1行, 消去λI\sqrt{\lambda}IλI的(n−1,n−1)(n-1,n-1)(n−1,n−1)元. 这一步旋转会产生λI\sqrt{\lambda}IλI的(n−1,n)(n-1,n)(n−1,n)位置上的填充, 而这可通过旋转RRR的第nnn行与λI\sqrt{\lambda}IλI的第n−1n-1n−1行消去;
3. 旋转RRR的第n−2n-2n−2行与λI\sqrt{\lambda}IλI的第n−2n-2n−2行, 消去λI\sqrt{\lambda}IλI的(n−2,n−2)(n-2,n-2)(n−2,n−2)元. 这一步旋转会产生λI\sqrt{\lambda}IλI的(n−2,n−1),(n−2,n)(n-2,n-1),(n-2,n)(n−2,n−1),(n−2,n)位置上的填充, 而这可通过先后旋转RRR的第n−1n-1n−1行与λI\sqrt{\lambda}IλI的第n−2n-2n−2行、RRR的第nnn行与λI\sqrt{\lambda}IλI的第n−2n-2n−2行消去.
4. …
依此类推. 若将所有的Givens变换汇成一个矩阵Qˉλ\bar{Q}_{\lambda}Qˉλ, 我们就有QˉλT[ROλI]=[RλOO],\bar{Q}_{\lambda}^T\begin{bmatrix}R\\O\\\sqrt{\lambda}I\end{bmatrix}=\begin{bmatrix}R_{\lambda}\\O\\O\end{bmatrix},QˉλT⎣⎡ROλI⎦⎤=⎣⎡RλOO⎦⎤,因此前面的正交矩阵QλQ_{\lambda}Qλ就是Qλ=[QI]Qˉλ.Q_{\lambda}=\begin{bmatrix}Q&\\&I\end{bmatrix}\bar{Q}_{\lambda}.Qλ=[QI]Qˉλ.这一方法的优点在于, 当我们在求根时会改变λ\lambdaλ的值, 而这样我们就只需要再计算Qˉλ\bar{Q}_{\lambda}Qˉλ而无需再管Householder变换的部分. 这在m≫nm\gg nm≫n时可以节省很多的计算量: 对λ\lambdaλ计算Qˉλ\bar{Q}_{\lambda}Qˉλ与RλR_{\lambda}Rλ仅需O(n3)O(n^3)O(n3)次运算, 而计算QQQ则需O(mn2)O(mn^2)O(mn2)次运算.
尺度变换.
最小二乘问题往往尺度较为恶性, 比如一些变量可能会达10410^4104量阶, 而其他的一些又会小到10−610^{-6}10−6. 若我们忽略如此巨大的差距, 算法就会不稳定或者产生一些不好的解. 一种减缓尺度带来的问题的途径是, 选取适当的椭球型信赖域代替上述的球形信赖域. 此时信赖域子问题变为:min⁡p12∥Jkp+rk∥2,s.t. ∥Dkp∥≤Δk,\min_p\frac{1}{2}\Vert J_kp+r_k\Vert^2,\quad \mathrm{s.t.\,}\Vert D_kp\Vert\le\Delta_k,pmin21∥Jkp+rk∥2,s.t.∥Dkp∥≤Δk,其中DkD_kDk为对角元为正的对角阵. 相应地, 解满足(JkTJk+λDk2)pkLM=−JkTrk,(J_k^TJ_k+\lambda D_k^2)p_k^{\mathrm{LM}}=-J_k^Tr_k,(JkTJk+λDk2)pkLM=−JkTrk,这又等价于求解下面的线性最小二乘问题min⁡p∥[JkλDk]p+[rkO]∥2.\min_p\left\Vert\begin{bmatrix}J_k\\\sqrt{\lambda}D_k\end{bmatrix}p+\begin{bmatrix}r_k\\O\end{bmatrix}\right\Vert^2.pmin∥∥∥∥[JkλDk]p+[rkO]∥∥∥∥2.这里对角阵DkD_kDk可随迭代改变, 其依据为xxx的每个分量的典型范围信息. 若变动在一定范围内, 则球形情形的收敛理论就仍然适用, 其中仅需稍微做些修正. 进一步地, 以上计算RλR_{\lambda}Rλ的步骤无需改动. Seber与Wild表示可以选取Dk2D_k^2Dk2为JkTJkJ_k^TJ_kJkTJk的对角元, 从而使得算法在xxx的对角尺度变换下不变. 这与第四章中缩放Hessian对角元的方法类似.
大规模问题.
而对于m,nm,nm,n都较大以及J(x)J(x)J(x)稀疏的问题, 我们更倾向于使用第七章CG-Steihaug算法求解, 其中以JkTJkJ_k^TJ_kJkTJk代替真实的∇2fk\nabla^2f_k∇2fk. JkTJkJ_k^TJ_kJkTJk的半正定性可用来简化算法, 这是因为原本算法中着重考虑的负曲率不会出现. 同时我们也不需要显式地去计算JkTJkJ_k^TJ_kJkTJk, 而是先后做两次矩阵-向量乘积.

3.2.3 Levenberg-Marquardt法的收敛性

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

定理3 设信赖域算法中η∈(0,14)\eta\in(0,\frac{1}{4})η∈(0,41), 水平集L\mathcal{L}L有界, 残差函数rj(⋅),j=1,2,…,mr_j(\cdot),j=1,2,\ldots,mrj(⋅),j=1,2,…,m在L\mathcal{L}L的一个邻域N\mathcal{N}N中Lipschitz连续可微. 假设对每个kkk, 近似解pkp_kpk满足不等式mk(0)−mk(pk)≥c1∥JkTrk∥min⁡(Δk,∥JkTrk∥∥JkTJk∥),m_k(0)-m_k(p_k)\ge c_1\Vert J_k^Tr_k\Vert\min\left(\Delta_k,\frac{\Vert J_k^Tr_k\Vert}{\Vert J_k^TJ_k\Vert}\right),mk(0)−mk(pk)≥c1∥JkTrk∥min(Δk,∥JkTJk∥∥JkTrk∥),其中c1>0,∥pk∥≤γΔk,γ≥1c_1>0,\Vert p_k\Vert\le\gamma\Delta_k,\gamma\ge1c1>0,∥pk∥≤γΔk,γ≥1. 于是有lim⁡k→∞∇fk=lim⁡k→∞JkTrk=0.\lim_{k\to\infty}\nabla f_k=\lim_{k\to\infty}J_k^Tr_k=0.k→∞lim∇fk=k→∞limJkTrk=0.

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

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

3.3 大残差问题的算法

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

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

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

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

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

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

第二种结合Gauss-Newton法与拟牛顿法的方式是仅保存Hessian二阶部分的近似. 也就是说, 我们保留矩阵列{Sk}\{S_k\}{Sk}近似∑j=1mrj(xk)∇2rj(xk)\sum_{j=1}^mr_j(x_k)\nabla^2r_j(x_k)∑j=1mrj(xk)∇2rj(xk), 然后以之估计整个HessianBk=JkTJk+Sk,B_k=J_k^TJ_k+S_k,Bk=JkTJk+Sk,接着使用信赖域或线搜索计算pkp_kpk. 对SkS_kSk的更新要求近似矩阵BkB_kBk或其组成部分能较好地模拟刚走完的那一步的特征. 更新公式则基于割线方程. 这里有许多不同的方式可以用来定义割线方程以及设计其他的条件得到SkS_kSk的更新公式. 下面我们介绍Dennis, Gay和Welsch的算法.
理想状态下, Sk+1S_{k+1}Sk+1应当很接近于x=xk+1x=x_{k+1}x=xk+1处的二阶项:Sk+1≈∑j=1mrj(xk+1)]∇2rj(xk+1).S_{k+1}\approx\sum_{j=1}^mr_j(x_{k+1})]\nabla^2r_j(x_{k+1}).Sk+1≈j=1∑mrj(xk+1)]∇2rj(xk+1).我们不想计算右端的∇2rj\nabla^2r_j∇2rj, 因此我们可代之以某个近似(Bj)k+1(B_j)_{k+1}(Bj)k+1并在(Bj)k+1(B_j)_{k+1}(Bj)k+1上提些条件. 也即(Bj)k+1(xk+1−xk)=∇rj(xk+1)−∇rj(xk)=(row j of J(xk+1))T−(row j of J(xk))T.\begin{aligned}(B_j)_{k+1}(x_{k+1}-x_k)&=\nabla r_j(x_{k+1})-\nabla r_j(x_k)\\&=(\mathrm{row\,}j\mathrm{\,of\,}J(x_{k+1}))^T-(\mathrm{row\,}j\mathrm{\,of\,}J(x_k))^T.\end{aligned}(Bj)k+1(xk+1−xk)=∇rj(xk+1)−∇rj(xk)=(rowjofJ(xk+1))T−(rowjofJ(xk))T.这一条件最终推出Sk+1S_{k+1}Sk+1上的割线方程:Sk+1(xk+1−xk)=∑j=1mrj(xk+1)(Bj)k+1(xk+1−xk)=∑j=1mrj(xk+1)[(row j of J(xk+1))T−(row j of J(xk))T]=Jk+1Trk+1−JkTrk+1.\begin{aligned}S_{k+1}(x_{k+1}-x_k)&=\sum_{j=1}^mr_j(x_{k+1})(B_j)_{k+1}(x_{k+1}-x_k)\\&=\sum_{j=1}^mr_j(x_{k+1})[(\mathrm{row\,}j\mathrm{\,of\,}J(x_{k+1}))^T-(\mathrm{row\,}j\mathrm{\,of\,}J(x_k))^T]\\&=J_{k+1}^Tr_{k+1}-J_k^Tr_{k+1}.\end{aligned}Sk+1(xk+1−xk)=j=1∑mrj(xk+1)(Bj)k+1(xk+1−xk)=j=1∑mrj(xk+1)[(rowjofJ(xk+1))T−(rowjofJ(xk))T]=Jk+1Trk+1−JkTrk+1.当然割线方程还不能完全决定Sk+1S_{k+1}Sk+1. Dennis, Gay和Welsch又要求Sk+1S_{k+1}Sk+1对称并且Sk+1−SkS_{k+1}-S_kSk+1−Sk要在某种意义下达到极小, 从而得到了下面的更新公式:Sk+1=Sk+(y#−Sks)yT+y(y#−Sks)TyTs−(y#−Sks)Ts(yTs)2yyT,S_{k+1}=S_k+\frac{(y^\#-S_ks)y^T+y(y^\#-S_ks)^T}{y^Ts}-\frac{(y^\#-S_ks)^Ts}{(y^Ts)^2}yy^T,Sk+1=Sk+yTs(y#−Sks)yT+y(y#−Sks)T−(yTs)2(y#−Sks)TsyyT,其中s=xk+1−xk,y=Jk+1Trk+1−JkTrk,y#=Jk+1Trk+1−JkTrk+1.\begin{aligned}s&=x_{k+1}-x_k,\\y&=J_{k+1}^Tr_{k+1}-J_k^Tr_k,\\y^\#&=J_{k+1}^Tr_{k+1}-J_k^Tr_{k+1}.\end{aligned}syy#=xk+1−xk,=Jk+1Trk+1−JkTrk,=Jk+1Trk+1−JkTrk+1.注意上述更新公式仅是DFP更新公式的微小改动版本. 若y#,yy^\#,yy#,y相同, 则二者就完全一样了. Dennis, Gay和Welsch将他们的近似HessianJkTJk+SkJ_k^TJ_k+S_kJkTJk+Sk联合信赖域的框架使用, 不过其中需要加以更多的约束以提升表现. 其基本策略的一个缺陷在于, SkS_kSk的更新策略并不保证当迭代点趋近一个零残差解时SkS_kSk会消失, 因此有时难以保证超线性收敛性. 这一问题可通过在SkS_kSk更新之前对其缩放避免: 我们以τkSk\tau_kS_kτkSk代替SkS_kSk, 其中τk=min⁡(1,∣sTy#∣∣sTSks∣).\tau_k=\min\left(1,\frac{|s^Ty^\#|}{|s^TS_ks|}\right).τk=min(1,∣sTSks∣∣sTy#∣).最后, 若Gauss-Newton法能产生较好的下降, 我们就应省去Hessian近似中的SkS_kSk.

4. 正交距离回归

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

下面我们从数学上严格地进行表述. 引入tjt_jtj上的扰动δj\delta_jδj, yjy_jyj上的扰动ϵj\epsilon_jϵj. 我们需要求得这2m2m2m个扰动的值, 使之极小化模型和观测的差异. 这里差异以加权最小二乘目标函数度量. 具体地说, 定义极小化问题为min⁡x,δj,ϵj12∑j=1mwj2ϵj2+dj2δj2,s.t. yj=ϕ(x;tj+δj)+ϵj,j=1,2,…,m.\min_{x,\delta_j,\epsilon_j}\frac{1}{2}\sum_{j=1}^mw_j^2\epsilon_j^2+d_j^2\delta_j^2,\quad\mathrm{s.t.\,}y_j=\phi(x;t_j+\delta_j)+\epsilon_j,\quad j=1,2,\ldots,m.x,δj,ϵjmin21j=1∑mwj2ϵj2+dj2δj2,s.t.yj=ϕ(x;tj+δj)+ϵj,j=1,2,…,m.wi,diw_i,d_iwi,di为权重, 它们由使用者或一些自动估计误差项的相对重要性的方式决定.
当我们图示上述问题时, 我们便能知道"正价距离回归"这个称谓是怎么来的了.

若所有的wi,diw_i,d_iwi,di都相等, 则目标函数中求和的每一项就是点(tj,yj)(t_j,y_j)(tj,yj)与曲线ϕ(x;t)\phi(x;t)ϕ(x;t)之间的最短距离. 而点与曲线之间的最短直线就与曲线在与直线的交点处正交.
我们可以用约束消去ϵj\epsilon_jϵj, 从而得到无约束最小二乘问题min⁡x,δF(x,δ)=12∑j=1mwj2[yj−ϕ(x;tj+δj)]2+dj2δj2=12∑j=12mrj2(x,δ),\min_{x,\delta}F(x,\delta)=\frac{1}{2}\sum_{j=1}^mw_j^2[y_j-\phi(x;t_j+\delta_j)]^2+d_j^2\delta_j^2=\frac{1}{2}\sum_{j=1}^{2m}r_j^2(x,\delta),x,δminF(x,δ)=21j=1∑mwj2[yj−ϕ(x;tj+δj)]2+dj2δj2=21j=1∑2mrj2(x,δ),其中δ=(δ1,δ2,…,δm)T\delta=(\delta_1,\delta_2,\ldots,\delta_m)^Tδ=(δ1,δ2,…,δm)T, rj(x,δ)={wj[ϕ(x;tj+δj)−yj],j=1,2,…,m,dj−mδj−m,j=m+1,…,2m.r_j(x,\delta)=\left\{\begin{array}{ll}w_j[\phi(x;t_j+\delta_j)-y_j], & j=1,2,\ldots,m,\\d_{j-m}\delta_{j-m}, & j=m+1,\ldots,2m.\end{array}\right.rj(x,δ)={wj[ϕ(x;tj+δj)−yj],dj−mδj−m,j=1,2,…,m,j=m+1,…,2m.此时的问题就是带有2m2m2m个残差和m+nm+nm+n个未知量的标准最小二乘问题. 我们可以使用本章中介绍的算法求解之. 不过直接使用可能会带来计算量上的困难, 这是因为此时的问题所带的参数数目与观测数可能原始的问题要大得多.

但若我们进一步深究Gauss-Newton法或Levenberg-Marquardt法中Jacobi矩阵的结构, 我们会发现: 它有许多元素都为0, 例如∂rj∂δi=∂[ϕ(tj+δj;x)−yj]∂δi=0,i,j=1,2,…,m,i≠j,∂rj∂xi=0,j=m+1,…,2m,i=1,2,…,n,∂rm+j∂δi={dji=j,0otherwise,i,j=1,2,…,m.\frac{\partial r_j}{\partial \delta_i}=\frac{\partial[\phi(t_j+\delta_j;x)-y_j]}{\partial\delta_i}=0,\quad i,j=1,2,\ldots,m,i\ne j,\\\frac{\partial r_j}{\partial x_i}=0,\quad j=m+1,\ldots,2m,\quad i=1,2,\ldots,n,\\\frac{\partial r_{m+j}}{\partial\delta_i}=\left\{\begin{array}{ll}d_j & i=j,\\0 & \mathrm{otherwise},\end{array}\right.\quad i,j=1,2,\ldots,m.∂δi∂rj=∂δi∂[ϕ(tj+δj;x)−yj]=0,i,j=1,2,…,m,i̸=j,∂xi∂rj=0,j=m+1,…,2m,i=1,2,…,n,∂δi∂rm+j={dj0i=j,otherwise,i,j=1,2,…,m.因此我们可以将Jacobi矩阵分块写作J(x,δ)=[J^VOD],J(x,\delta)=\begin{bmatrix}\hat{J} & V\\O & D\end{bmatrix},J(x,δ)=[J^OVD],其中V,DV,DV,D为m×mm\times mm×m对角阵, J^\hat{J}J^为m×nm\times nm×n矩阵, 其中元素为wjϕ(tj+δj;x)w_j\phi(t_j+\delta_j;x)wjϕ(tj+δj;x)对xxx的偏导. 对应地, 将p,rp,rp,r也分块为p=[pxpδ],r=[r^1r^2],p=\begin{bmatrix}p_x\\p_{\delta}\end{bmatrix},\quad r=\begin{bmatrix}\hat{r}_1\\\hat{r}_2\end{bmatrix},p=[pxpδ],r=[r^1r^2],并将正规方程分块为[J^TJ^J^TVVJ^V2+D2+λI][pxpδ]=−[J^Tr^1Vr^1+Dr^2].\begin{bmatrix}\hat{J}^T\hat{J} & \hat{J}^TV\\V\hat{J} & V^2+D^2+\lambda I\end{bmatrix}\begin{bmatrix}p_x\\p_{\delta}\end{bmatrix}=-\begin{bmatrix}\hat{J}^T\hat{r}_1\\V\hat{r}_1+D\hat{r}_2\end{bmatrix}.[J^TJ^VJ^J^TVV2+D2+λI][pxpδ]=−[J^Tr^1Vr^1+Dr^2].由于右下方的子阵V2+D2+λIV^2+D^2+\lambda IV2+D2+λI对角, 因此我们可以方便地消去pδp_{\delta}pδ得到用来求pxp_xpx的n×nn\times nn×n子系统. 这样求得一步的计算量就仅比标准最小二乘模型下的m×nm\times nm×n问题稍微大一点.

5. 再谈大规模情形

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

Numerical Optimization Ch10. Least-Squares Problems相关推荐

METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 翻译（七）
METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS(七) 3.7. 最后的话我们已经讨论了许多解决非线性最小二乘问题的算法.它们都出现在任何好的程序库中,并且 ...
Numerical Optimization - my afterword
历时六个月, 从第一篇(2018.9.30)到第十九篇(2019.3.17), 感谢各位博友的支持. 就个人而言, 其实这本书早在2018.11就看完了. 写博客纯粹是为了加深自己的印象.锻炼自己的英 ...
Note of Numerical Optimization Ch.3
目录 Numerical Optimization Ch.3 Line Search Methods Step Length Convergence of Line Search Methods Ra ...
数学知识--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 ...
数值优化（Numerical Optimization）学习系列-文件夹
概述数值优化对于最优化问题提供了一种迭代算法思路,通过迭代逐渐接近最优解,分别对无约束最优化问题和带约束最优化问题进行求解. 该系列教程能够參考的资料有 1. <Numerical Optim ...
支持向量机：Numerical Optimization
支持向量机:Numerical Optimization by pluskid, on 2010-09-15, in Machine Learning 15 comments 本文是&q ...
数值优化（Numerical Optimization）学习系列-目录
概述数值优化对于最优化问题提供了一种迭代算法思路,通过迭代逐渐接近最优解,分别对无约束最优化问题和带约束最优化问题进行求解. 该系列教程可以参考的资料有 1. <Numerical Optim ...

Numerical Optimization Ch10. Least-Squares Problems