Numerical Optimization Ch11. Nonlinear Equations

第十一章: 非线性方程组

在许多应用场景下, 我们 无需显式地优化某个目标函数, 而是需要求得模型中满足一些给定关系式的变量(或参数)的值. 当这些关系式构成了 $n$个等式时——与模型中变量的 个数相同——此问题就变成了求解 非线性方程组(nonlinear equations). 数学上我们将这个问题写作 $$r(x)=0,$$其中 $r:{\mathbb{R}}^n\to {\mathbb{R}}^n$为一向量函数, 即 $$r(x)=\left[\begin{array}{c}{r}_1(x)\\ r_2(x)\\ ⋮\\ r_n(x)\end{array}\right].$$本章中, 我们假设每个函数 $r_i:{\mathbb{R}}^n\to \mathbb{R},i=1,2,\dots ,n$ 光滑. 满足 $r(x)=0$的向量 $x^{\ast }$称为非线性方程组的 解(solution). 一个简单的例子是 $$r(x)=\left[\begin{array}{c}{x}_2^2-1\\ \sin x_1-x_2\end{array}\right]=0,$$这就是个 有无穷多解的一个 $n=2$的方程组, 其中的两个解为 $x^{\ast }=(3\pi /2,-1{)}^T,x^{\ast }=(\pi /2,1{)}^T$. 一般的, 非线性方程组可以 无解、有唯一解或有多解.

求解非线性方程组的方法在动机、分析与实施上都与先前讨论的优化方法密切相关. 特别地, 在优化和非线性方程组的求解中, 牛顿法都处于许多重要算法的核心. 每步迭代的线搜索、信赖域以及线性的代数子问题的非精确求解在两个领域中都很重要. 二者也同时涉及如导数计算、全局收敛等问题.

由于许多求解非线性方程组的重要算法都是通过极小化方程离差平方和, 即求解$$\underset{x}{\min }\sum _{i=1}^n{r}_i^2(x),$$因此它们与第十章中求解非线性最小二乘问题的算法联系密切. 其中的差异在于:

1. 在求解非线性方程组时, 方程的个数与变量数是相等的;
2. 我们期望解满足所有的方程, 而不是仅仅极小化离差的平方和. 由于非线性方程组往往表征了一些生理的或经济意义上的约束, 例如保守性规则或一致性原理, 因此令解精确地满足方程组就显得尤为重要.

许多应用需要我们求解一系列紧密相关的非线性系统, 例如下面这个例子.

例1 控制学领域中的一个问题是要分析飞行器对于飞行员指令响应的稳定性. 以下是基于受力平衡方程的一个简化模型, 其中我们忽略了重力项.

某个飞行器的平衡方程组由5个方程和8个未知量组成:$$F(x)\equiv Ax+\varphi (x)=0,$$其中$F:{\mathbb{R}}^8\to {\mathbb{R}}^5$, 矩阵$A$为$$A=\left[\begin{array}{cccccccc}-3.933& 0.107& 0;126& 0& -9.99& 0& -45.83& -7.64\\ 0& -0.987& 0& -22.95& 0& -28.37& 0& 0\\ 0.002& 0& -0.235& 0& 5.67& 0& -0.921& -6.51\\ 0& 1.0& 0& -1.0& 0& -0.168& 0& 0\\ 0& 0& -1.0& 0& -0.196& 0& -0.0071& 0\end{array}\right],$$而非线性部分则定义为$$\varphi (x)=\left[\begin{array}{c}-0.727x_2{x}_3+8.39x_3{x}_4-684.4x_4{x}_5+63.5x_4{x}_2\\ 0.949x_1{x}_3+0.173x_1{x}_5\\ -0.716x_1{x}_2-1.578x_1{x}_4+1.132x_4{x}_2\\ -x_1{x}_5\\ x_1{x}_4\end{array}\right].$$前三个变量$x_1,x_2,x_3$分别表示横摇、纵摇和横摆的速度, $x_4$表示增量迎角, $x_5$表示侧滑角. 最后三个变量$x_6,x_7,x_8$分别表示控制的参数——升降舵、副翼和方向舵的偏转.

给定控制变量$x_6,x_7,x_8$, 我们就得到了具有5个方程和5个未知数的系统. 若我们希望研究飞行器的行为随控制变量的变化, 我们就需要对每一组控制参数求解相应的$x_1,x_2,\dots ,x_5$的非线性方程组.

尽管非线性方程组、无约束优化以及最小二乘问题的求解算法有许多的相似之处, 但它们之间也存在重大的差异:

1. 导数的计算. 优化中我们若要求算法具有二次收敛性, 我们往往需要目标函数的二阶导数, 而在求解非线性方程组时仅需一阶导数;
2. 价值函数的选取. 在无约束优化中, 目标函数是我们用来评估求解进程的价值函数的一个天然的选择, 但在求解非线性方程组时, 我们可以使用各种各样的价值函数, 而它们也都各自具有一些缺陷;
3. 框架的偏好. 线搜索和信赖域方法在优化中可以说是同等重要, 但在求解非线性方程组时信赖域方法在理论上显得更优.

求解非线性方程组的一些难点可用一个简单的一维例子阐明. 例如$$r(x)=\sin (5x)-x,$$它的图像如下所示.

从图中我们看到, $r(x)=0$具有三个解(也称$r$有三个根), 分别为0以及(大约)$\pm 0.519148$. 多解的情形类似于优化问题中函数具有多个局部极小点. 但二者并不等同:

• 优化中, 某个局部极小可能比其他的局部极小的函数值更小, 即局部极小之间是存在好坏之分的; 而在非线性方程组中, 这些解(从数学的角度上看)都是一样的. 当然从现实的角度或者从需求的角度, 我们会排除一些解; 甚至如果所有的解都不符合实际, 我们就需要重新建模.

本章我们先概览与牛顿法相关的算法, 并讨论它们的局部收敛性质. 这些算法除了牛顿法本身, 还包括Broyden的拟牛顿法、非精确牛顿法和张量法. 之后我们讨论全局收敛性. 最后我们讨论一种连续(或同伦)的方法. 它的主要思想是: 逐步求解较为"简单"的问题——它们的解是显然的或易求的——并逐渐将这些"简单"问题转换为$F(x)=0$. 随着问题的转换, 我们将跟踪关注解的改变, 并希望最终得到$F(x)=0$的解. 这种方法的思想在之后带约束优化中的内点法的设计上显得尤为关键.

本章中我们假设向量函数$r$在某个定义域$\mathcal{D}$上连续可微. 换句话说, 就是Jacobi矩阵$J(x)$存在且连续. 我们称满足$r(x^{\ast })=0$的$x^{\ast }$为退化解(degenerate solution), 若相应的$J(x^{\ast })$奇异; 否则称之为非退化解(nondegenerate solution).

1. 局部算法

1.1 非线性方程组的牛顿法

回忆之前在极小$f:{\mathbb{R}}^n\to \mathbb{R}$的过程中, 我们构建了$f$在当前迭代点$x_k$附近的二次模型, 之后导出了牛顿法, 牛顿步就是极小化该模型的向量. 而在非线性方程组的情形, 牛顿法也是以类似的方式导出, 不过取而代之的是线性模型. 我们先陈述一下Taylor定理的向量值函数变体.

定理1 设$r:{\mathbb{R}}^n\to {\mathbb{R}}^n$在某个开凸集$\mathcal{D}$上连续可微, $x,x+p$均在$\mathcal{D}$中. 于是我们有$$r(x+p)=r(x)+{\int }_0^{1}J(x+tp)p\mathrm{d}t.$$
因此我们可定义$r(x_k+p)$的一个线性模型$M_k(p)$:$$M_k(p)\stackrel{def}{=}r(x_k)+J(x_k)p.$$而最原始的牛顿法就是通过取搜索方向为使得$M_k(p_k)=0$的$p_k=-J(x_k{)}^{-1}r(x_k)$. 下面是正式的算法描述.

算法1 (Newton’s Method for Nonlinear Equations)
Choose $x_0$;
for $k=0,1,2,\dots$
\quad\quad Calculate a solution $p_k$ to the Newton equations$$J(x_k)p_k=-r(x_k);$$ \quad\quad $x_{k+1}\leftarrow x_k+p_k$;
end (for)

在非线性方程组的情形下, 我们不用二次模型而仅用线性模型计算牛顿步, 是因为线性模型能够正常地得到解且相应的算法具有较快的收敛速度. 事实上, 无约束优化问题的牛顿法亦可通过将算法1实施于非线性方程组$\nabla f(x)=0$得到. 在后面讨论带等式约束的逐步二次规划(Sequence Quadratic Programming, SQP)时, 我们依然能看到算法1的身影, 此时的方程对应问题的一阶最优性条件. 它与非线性最小二乘的Gauss-Newton法也有关联: 当$J(x_k)$非奇异时二者完全等价.

我们之后会证明, 当$x_k$距离某个非退化解$x^{\ast }$足够近时, 牛顿法具有超线性收敛性. 牛顿法可能存在的缺陷包括:

• 当初始点距离解较远时, 算法1的表现可能会相当病态. 特别地, 当$J(x_k)$奇异, 我们甚至无法定义牛顿步.
• 一阶导数信息(即Jacobi矩阵$J$)计算困难.
• 当$n$较大, 计算牛顿步$p_k$的计算量可能会令人难以接受.
• 解$x^{\ast }$可能是退化解, 即$J(x^{\ast })$可能是奇异的.

对于第四条, 求标量函数$r(x)=x^2$的根就是一个退化问题. 它有退化根$x^{\ast }=0$. 此时算法1从任意非零$x_0$出发, 均将产生以下迭代点:$$x_k=\frac{1}{2^k}{x}_0,$$这显然仅线性收敛到解$0$.

后面我们会介绍, 牛顿法可通过各种方式修正或增强以解决这里的大部分问题. 而这些变体就构成了可求解非线性方程组的许多软件的基础.
我们在下面的定理中将阐明算法1的局部收敛性质. 这一部分, 我们假设Jacobi矩阵Lipschitz连续, 从而存在常数${\beta }_L$使得$$\Vert J(x_0)-J(x_1)\Vert \le {\beta }_L\Vert x_0-x_1\Vert ,$$其中$x_0,x_1$为定义域中的任意点.

定理2 设$r$在一开凸集$\mathcal{D}\subset {\mathbb{R}}^n$上连续可微. $x^{\ast }\in \mathcal{D}$为$r(x)=0$的一个非退化解, { x k } \{x_k\} {xk​}为由算法1产生的迭代序列. 则当$x_k\in \mathcal{D}$距离$x^{\ast }$充分近时, 我们有$$x_{k+1}-x^{\ast }=o(\Vert x_k-x^{\ast }\Vert ),$$即有局部Q-超线性收敛. 进一步, 当$r$在$x^{\ast }$附近Lipschitz连续可微时, 我们有$$x_{k+1}-x^{\ast }=O(\Vert x_k-x^{\ast }{\Vert }^2),$$即有局部Q-二次收敛.

证明: 由于$r(x^{\ast })=0$, 因此由定理1有$$r(x_k)=r(x_k)-r(x^{\ast })=J(x_k)(x_k-x^{\ast })+w(x_k,x^{\ast }),$$其中$$w(x_k,x^{\ast })={\int }_0^1[J(x_k+t(x^{\ast }-x_k))-J(x_k)](x_k-x^{\ast })\mathrm{d}t.$$由$J$的连续性, 我们有$$\begin{array}{rl}\Vert w(x_k,x^{\ast })\Vert & =\Vert {\int }_0^1[J(x_k+t(x^{\ast }-x_k))-J(x_k)](x_k-x^{\ast })\mathrm{d}t\Vert \\ & \le {\int }_0^1\Vert J(x_k+t(x^{\ast }-x_k))-J(x_k)\Vert \Vert x_k-x^{\ast }\Vert \mathrm{d}t\\ & =o(\Vert x_k-x^{\ast }\Vert ).\end{array}$$由于$J(x^{\ast })$非奇异, 因此存在半径$\delta >0$以及正数${\beta }^{\ast }$使得对于 B ( x ∗ , δ ) = d e f { x ∣ ∥ x − x ∗ ∥ ≤ δ } \mathcal{B}(x^*,\delta)\xlongequal{def}\{x|\Vert x-x^*\Vert\le\delta\} B(x∗,δ)def {x∣∥x−x∗∥≤δ}中所有的$x$, 均有$$\Vert J(x{)}^{-1}\Vert \le {\beta }^{\ast }且x\in \mathcal{D}.$$假设$x_k\in \mathcal{B}(x^{\ast },\delta )$, 在Taylor展开式两端左乘$J(x_k{)}^{-1}$得到$$\begin{array}{rl}-p_k& =(x_k-x^{\ast })+\Vert J(x_k{)}^{-1}\Vert o(\Vert x_k-x^{\ast }\Vert ),\\ \Rightarrow x_k+p_k-x^{\ast }& =o(\Vert x_k-x^{\ast }\Vert ),\\ \Rightarrow x_{k+1}-x^{\ast }& =o(\Vert x_k-x^{\ast }\Vert ),\end{array}$$这就得到了超线性收敛.
进一步, 当$r$在$x^{\ast }$附近Lipschitz连续可微, 我们就有$\Vert w(x_k,x^{\ast })\Vert =O(\Vert x_k-x^{\ast }{\Vert }^2)$. 重复上述分析, 便得到二次收敛. 证毕.

1.2 非精确牛顿法

下面我们不要求精确求解牛顿方程, 而是仅求非精确搜索方向$p_k$, 要求它满足$$\Vert r_k+J_k{p}_k\Vert \le {\eta }_k\Vert r_k\Vert ,{\eta }_k\in [0,\eta ],$$其中$\eta \in [0,1)$为常数. 同第七章一样, 我们称 { η k } \{\eta_k\} {ηk​}为强迫序列. 不同的方法会选取不同的强迫序列, 且它们会使用不同的算法求近似解$p_k$. 此类方法的一般框架陈述如下:

框架2 (Inexact Newton for Nonlinear Equations)
Given $\eta \in [0,1)$;
Choose $x_0$;
for $k=0,1,2,\dots$
\quad\quad Choose forcing parameter ${\eta }_k\in [0,\eta ]$;
\quad\quad Find a vector $p_k$: $\Vert r_k+J_k{p}_k\Vert \le {\eta }_k\Vert r_k\Vert$;
\quad\quad $x_{k+1}\leftarrow x_k+p_k$;
end (for)

这些方法的收敛性质仅依赖于$p_k$满足的条件, 而与如何计算$p_k$无关. 其中最重要的方法使用如广义极小残量法(GMRES)或其他的Krylov子空间方法求解线性系统$Jp=-r$. 同第五章的共轭梯度法一样(CG不可以直接应用到这里, 因为$J$不保证对称正定), 这些方法在每步迭代仅需计算矩阵-向量乘积$Jd$, 并储存一定数量的维数为$n$的的工作向量. GMRES每一步都会多储存一个向量, 因此使用它时必须要周期性地重启(一般每10或20次迭代重启一次)以控制内存的使用在合理的水平.

矩阵-向量乘积$Jd$也可不需Jacobi矩阵$J$的显式表示就能计算. 第八章我们介绍了如何使用有限差分近似计算$Jd$. 精确(或至少在有限精度计算的极限内)计算$Jd$可用自动微分的前向模式, 其计算代价至多为计算一次$r(\cdot )$的计算量的小倍数.

这里我们不介绍稀疏线性系统的迭代方法. 下面我们证明非精确牛顿法的局部收敛定理.

定理3 设$r$在一开凸集$\mathcal{D}\subset {\mathbb{R}}^n$中连续可微.$x^{\ast }\in \mathcal{D}$为$r(x)=0$的非退化解, { x k } \{x_k\} {xk​}为由框架2产生的迭代序列. 则当$x_k\in \mathcal{D}$距离$x^{\ast }$充分近时, 有以下结论成立:

1. 若强迫序列的上界$\eta$充分小, 则 { x k } \{x_k\} {xk​}Q-线性收敛于$x^{\ast }$.
2. 若${\eta }_k\to 0$, 则收敛是Q-超线性的.
3. 若进一步, $J(\cdot )$在$x^{\ast }$的某个邻域内Lipschitz连续, 且${\eta }_k=O(\Vert r_k\Vert )$, 则收敛是Q-二次的.

证明: 首先我们将$p_k$满足的条件写作$$J(x_k)p_k+r(x_k)=v_k,\Vert v_k\Vert \le {\eta }_k\Vert r(x_k)\Vert .$$由于$x^{\ast }$为一非退化解, 从而存在半径$\delta >0$以及正数${\beta }^{\ast }$使得$\Vert J(x{)}^{-1}\Vert \le {\beta }^{\ast },\forall x\in \mathcal{B}(x^{\ast },\delta )$. 上式两端左乘$J(x_k{)}^{-1}$可得$$\Vert p_k+J(x_k{)}^{-1}r(x_k)\Vert =\Vert J(x_k{)}^{-1}{v}_k\Vert \le {\beta }^{\ast }{\eta }_k\Vert r(x_k)\Vert .$$与定理2的证明同, 我们有$$r(x)=J(x)(x-x^{\ast })+w(x,x^{\ast }),$$其中$\rho (x)\stackrel{def}{=}\Vert w(x,x^{\ast })\Vert /\Vert x-x^{\ast }\Vert \to 0,x\to x^{\ast }$. 通过适当减小$\delta$, 我们就有$$\Vert r(x)\Vert \le 2\Vert J(x^{\ast })\Vert \Vert x-x^{\ast }\Vert +o(\Vert x-x^{\ast }\Vert )\le 4\Vert J(x^{\ast })\Vert \Vert x-x^{\ast }\Vert ,x\in \mathcal{B}(x^{\ast },\delta ).$$联合以上可得$$\begin{array}{rl}\Vert x_k+p_k-x^{\ast }\Vert & =\Vert p_k+J(x_k{)}^{-1}(r(x_k)-w(x_k,x^{\ast }))\Vert \\ & \le {\beta }^{\ast }{\eta }_k\Vert r(x_k)\Vert +\Vert J(x_k{)}^{-1}\Vert \Vert w(x_k,x^{\ast })\Vert \\ & \le \left[4\Vert J(x^{\ast })\Vert {\beta }^{\ast }{\eta }_k+{\beta }^{\ast }\rho (x_k)\right]\Vert x_k-x^{\ast }\Vert .\end{array}$$选取$x_k$充分靠近$x^{\ast }$使得$\rho (x_k)\le 1/(4{\beta }^{\ast })$, 选取$\eta =1/(8\Vert J(x^{\ast })\Vert {\beta }^{\ast })$, 我们就有上式方括号至多为$3/4$. 因此有Q-线性收敛. 这就证明了1.
2的证明是显然的, 因为方括号随着$x_k\to x^{\ast },{\eta }_k\to 0$也将趋近于0. 对于3, 证明手段与定理2的第二部分相同.

1.3 Broyden方法

相比于前面介绍的, 拟牛顿法并不需要计算Jacobi矩阵$J(x)$: 它们将自行构建对Jacobi的近似矩阵, 并在每步迭代更新使得近似矩阵能较好的模拟真实的Jacobi$J$在方才走过一步上的行为. 我们以$B_k$表示迭代第$k$步中的近似Jacobi矩阵, 于是有类似的线性模型:$$M_k(p)=r(x_k)+B_{k}p.$$$p_k$便是这个模型为0的解. 当$B_k$非奇异, 我们有以下显式公式:$$p_k=-B_k^{-1}r(x_k).$$以下考虑更新公式. 我们可以将近似矩阵对真实矩阵的模拟明确如下: 令$$s_k=x_{k+1}-x_k,y_k=r(x_{k+1})-r(x_k).$$由定理1, $$y_k={\int }_0^{1}J(x_k+t{s}_k)s_k\mathrm{d}t\approx J(x_{k+1})s_k+o(\Vert s_k\Vert ).$$于是我们要求更新的Jacobi近似$B_{k+1}$满足以下割线方程:$$y_k=B_{k+1}{s}_k,$$这就保证了$B_{k+1}$和$J(x_{k+1})$沿着方向$s_k$的行为是相似的. 注意到割线方程并没有提供$B_{k+1}$沿着与$s_k$正交的方向上的行为信息. 事实上, 我们可以将割线方程看成是具有$n^2$个未知量的$n$个线性方程构成的系统. 至此$B_{k+1}$并不被唯一确定. ($n=1$的标量情形就能推出标量的割线法(secant methods).)

确定$B_{k+1}$的方法有很多, 而这其中在应用上最成功的是Broyden法, 其更新公式为$$B_{k+1}=B_k+\frac{(y_k-B_k{s}_k)s_k^T}{s_k^T{s}_k}.$$Broyden更新的大致思想是, 极小化Jacobi矩阵的该变量(以欧式范数$\Vert B_k-B_{k+1}{\Vert }_2$度量), 同时要求满足割线方程.

引理4 在所有的满足割线方程$B{s}_k=y_k$的矩阵$B$中, 由Broyden公式确定的$B_{k+1}$极小矩阵改变量$\Vert B-B_k\Vert$.

证明: 设$B$为任意满足$B{s}_k=y_k$的矩阵. 由欧式范数与向量范数的相容性以及$\Vert s{s}^T/s^{T}s\Vert =1,\forall s$可得$$\begin{array}{rl}\Vert B_{k+1}-B_k\Vert & =\Vert \frac{(y_k-B_k{s}_k)s_k^T}{s_k^T{s}_k}\Vert \\ & =\Vert \frac{(B-B_k)s_k{s}_k^T}{s_k^T{s}_k}\Vert \le \Vert B-B_k\Vert \Vert \frac{s_k{s}_k^T}{s_k^T{s}_k}\Vert =\Vert B-B_k\Vert .\end{array}$$因此我们有$$B_{k+1}\in \arg \underset{B:y_k=B{s}_k}{\min }\Vert B-B_k\Vert ,$$证毕.

以下我们给出具体算法, 其中附加了沿着搜索方向$p_k$的线搜索过程, 即有$s_k=\alpha p_k,\alpha >0$.

算法3 (Broyden)
Choose $x_0$ and a nonsingular initial Jacobian approximation $B_0$;
for $k=0,1,2,\dots$
\quad\quad Calculate a solution $p_k$ to the linear equations$$B_k{p}_k=-r(x_k);$$ \quad\quad Choose ${\alpha }_k$ by performing a line search along $p_k$;
\quad\quad $x_{k+1}\leftarrow x_k+{\alpha }_k{p}_k$;
\quad\quad $s_k\leftarrow x_{k+1}-x_k$;
\quad\quad $y_k\leftarrow r(x_{k+1})-r(x_k)$;
\quad\quad Obtain $B_{k+1}$ from the update formula;
end (for)

在一定的假设条件下, Broyden法超线性收敛, 即$$\Vert x_{k+1}-x^{\ast }\Vert =o(\Vert x_k-x^{\ast }\Vert ).$$这一局部收敛速度虽不及牛顿法, 但对于大多数实际的应用已经够快了.

我们以一个例子来阐释牛顿法与Broyden法收敛速度的差别. 定义向量值函数$r:{\mathbb{R}}^2\to {\mathbb{R}}^2$为$$r(x)=\left[\begin{array}{c}(x_1+3)(x_2^3-7)+18\\ \sin (x_2{e}^{x_1}-1)\end{array}\right],$$它有一非退化根$x^{\ast }=(0,1{)}^T$. 两种方法下我们均以$x_0=(-0.5,1.4{)}^T$为初始点. 在Broyden法中选取精确的Jacobi矩阵$J(x_0)$作为初始矩阵. 结果如下表.

	xk-x*的2范数	xk-x*的2范数
Iteration $k$	Broyden	Newton
0	$0.64\times 10^0$	$0.64\times 10^0$
1	$0.62\times 10^{-1}$	$0.62\times 10^{-1}$
2	$0.52\times 10^{-3}$	$0.21\times 10^{-3}$
3	$0.25\times 10^{-3}$	$0.18\times 10^{-7}$
4	$0.43\times 10^{-4}$	$0.12\times 10^{-15}$
5	$0.14\times 10^{-6}$
6	$0.57\times 10^{-9}$
7	$0.18\times 10^{-11}$
8	$0.87\times 10^{-15}$

	r(x_k)的2范数	r(x_k)的2范数
Iteration $k$	Broyden	Newton
0	$0.74\times 10^1$	$0.74\times 10^1$
1	$0.59\times 10^0$	$0.59\times 10^0$
2	$0.20\times 10^{-2}$	$0.23\times 10^{-2}$
3	$0.21\times 10^{-2}$	$0.16\times 10^{-6}$
4	$0.37\times 10^{-3}$	$0.22\times 10^{-15}$
5	$0.12\times 10^{-5}$
6	$0.49\times 10^{-8}$
7	$0.15\times 10^{-10}$
8	$0.11\times 10^{-18}$

从每步迭代残差范数的幂次可知Newton法有Q-二次收敛性. 而Broyden法则用了牛顿法的2倍迭代步数, 且最后残差减小的速度逐渐加快. $\Vert r(x_k)\Vert$的表现与$\Vert x_k-x^{\ast }\Vert$是差不多的, 因为$$r(x_k)=r(x_k)-r(x^{\ast })\approx J(x^{\ast })(x_k-x^{\ast }),$$其中$J(x^{\ast })$非奇异.

Broyden法的收敛性分析要比牛顿法复杂得多. 我们不加证明的叙述之.
定理5 设定理2的假设条件成立. 则存在$ϵ>0,\delta >0$使得对任意满足$$\Vert x_0-x^{\ast }\Vert \le \delta ,\Vert B_0-J(x^{\ast })\Vert \le ϵ$$的初始点$x_0$和初始矩阵$B_0$, 由Broyden发产生的迭代序列 { x k } \{x_k\} {xk​}良定且Q-超线性收敛于$x^{\ast }$.

定理5中的第二个条件——初始Jacobi近似$B_0$必须距离$J(x^{\ast })$足够近——在实际中比较难保证. 与无约束极小问题相反, 在这里$B_0$的选取至关重要. 较为常见的选取方式为, 令$B_0=J(x_0)$或$J(x_0)$的某种差分近似.

一般地, 即使真实的Jacobi矩阵$J$很稀疏, Broyden矩阵$B_k$也可能会很稠密. 因此当$n$很大时, 完整储存与使用$B_k$的Broyden方法可能并不高效. 因此, 我们可以使用有限内存的方式, 以一定量的向量隐式地储存与代表$B_k$, 而求逆的操作可用Sherman-Morrison-Woodbury公式. 这与第七章的有限内存BFGS较为类似.

1.4 张量方法

在张量方法中, 我们需要给牛顿法中的使用的线性模型$M_k(p)$加上额外的项以捕捉$r$的一些非线性、更高阶的行为. 如此, 算法便能更快地、更可靠地收敛到退化解, 尤其是那些Jacobi矩阵$J(x^{\ast })$的秩为$n-1$或$n-2$的$x^{\ast }$. 这里我们仅简单介绍该方法.

我们记模型函数为${\hat{M}}_k(p)$:$${\hat{M}}_k(p)=r(x_k)+J(x_k)p+\frac{1}{2}{T}_{k}pp,$$其中$T_k$为张量(tensor), 具有$n^3$个元素$(T_k{)}_{ijl}$. 它在${\mathbb{R}}^n$中的任意一对向量上的作用定义为:$$(T_{k}uv{)}_i=\sum _{j=1}^n\sum _{l=1}^n(T_k{)}_{ijl}{u}_j{v}_l.$$若我们延续牛顿法的思路, 我们就会从$r$在点$x_k$处的二阶导数构造$T_k$:$$(T_k{)}_{ijl}=[{\nabla }^2{r}_i(x_k){]}_{jl}.$$例如就拿1.3节中的例子来说, 我们有$$\begin{array}{rl}(T(x)uv{)}_1& =u^T{\nabla }^2{r}_1(x)v=u^T\left[\begin{array}{cc}0& 3x_2^2\\ 3x_2^2& 6x_2(x_1+3)\end{array}\right]v\\ & =3x_2^2(u_1{v}_2+u_2{v}_1)+6x_2(x_1+3)u_2{v}_2.\end{array}$$然而, 在大多数情形下二阶导的计算并不实际. 若我们要显式地储存二阶导, 就需要储存大约$n^3/2$的内容. 这大概是牛顿法的$n$倍. 更有甚之, 可能根本不存在$p$使得${\hat{M}}_k(p)=0$, 从而这一步是无法定义的.

下面我们介绍一种仅需额外增加少量的储存就可定义$T_k$的方法. 它同时还能给予${\hat{M}}_k$一些潜在的性质. 具体说来, $T_k$的选取需要使得${\hat{M}}_k(p)$成功插值算法于前几次迭代产生的迭代点. 即$${\hat{M}}_k(x_{k-j}-x_k)=r(x_{k-j}),j=1,2,\dots ,q,q\in \mathbb{N}.$$将$T_k$的定义代入, 得到$T_k$必须满足条件$$\frac{1}{2}{T}_k{s}_{jk}{s}_{jk}=r(x_{k-j})-r(x_k)-J(x_k)s_{jk},$$其中$$s_{jk}\stackrel{def}{=}{x}_{k-j}-x_k,j=1,2,\dots ,q.$$我们可以证明当我们选取$T_k$满足$$T_{k}uv=\sum _{j=1}^q{a}_j(s_{jk}^{T}u)(s_{jk}^{T}v),$$其中$a_j,j=1,2,\dots ,q$为长度为$n$的向量时, 以上条件是满足的. 插值点数$q$通常较小, 一般小于$\sqrt{n}$. 此时$T_k$可用$2nq$个位置存储(即存储$a_j$和$s_{jk},j=1,2,\dots ,q$). 注意这一思想与Broyden方法的联系, 后者也是要求模型满足一定的插值条件(虽说仅是模型的一阶部分).

这一方法可用许多不同的方式加以完善. 例如, 选取插值点使得方向$s_{jk}$更加线性无关. 同时, 满足${\hat{M}}_k(p)$的$p$可能仍不存在, 但我们可以以极小化$\Vert {\hat{M}}_k(p){\Vert }_2^2$的$p$替代. 而这可用一种特殊的最小二乘技术得到. 由于我们无法保证这样求得的方向是价值函数$\frac{1}{2}\Vert r(x){\Vert }^2$的下降方向, 因此到时我们可以标准牛顿方向$-J_k^{-1}{r}_k$代之.

2. 实用方法

下面我们考虑以上介绍的牛顿类算法的实用变体, 其中将引入线搜索与信赖域以得到更好的全局收敛性质.

2.1 价值函数

如前所述, 不管是牛顿法还是带单位步长的Broyden方法, 我们都只能得到算法的局部收敛. 若从较远的位置出发迭代, 未知量、函数向量或是Jacobi矩阵的元素就有可能爆炸. 另外, 有种更加奇怪的算法行为是循环, 其中迭代点在不同的参数空间中移动却从不趋近某个解. 一个例子就是标量函数$$r(x)=-x^5+x^3+4x.$$它有五个非退化的根. 当我们从$x_0=1$出发, 牛顿法就会产生一列在1和-1之间振荡的迭代点.

加入线搜索与信赖域的框架, 牛顿法和Broyden方法就可变得更强健. 在介绍之前, 我们需要先定义一个标量价值函数(merit function), 它能(从迫近$r$的根的角度)指示新的迭代点是否比当前的迭代点更好, 从而指导算法接收新点. 无约束优化中, 目标函数$f$自身就是个不错的价值函数; 而在非线性方程组的求解中, 价值函数则以某种方式组合$r$的$n$个分量得到.

这里最常用的价值函数就是分量平方和:$$f(x)=\frac{1}{2}\Vert r(x){\Vert }^2=\frac{1}{2}\sum _{i=1}^n{r}_i^2(x).$$其中$\frac{1}{2}$是为了计算上的方便. 任何$r$的根$x^{\ast }$显然有$f(x^{\ast })=0$, 并且由于$f(x)\ge 0,\forall x$, 因此任意根都是$f$的最小点. 这样看, $f$的局部极小点有时就不是$r$的根. 即使这样, 价值函数仍然成功应用于许多实际问题.

$r(x)=\sin (5x)-x$的价值函数如下图所示.

图中可见, 价值函数有3个局部极小点对应3个根, 但还有许多局部极小并不对应根(例如在$\pm 1.53053$附近的极小点). 这些不是根的极小点满足一个有趣的性质: 由于$$\nabla f(x^{\ast })=J(x^{\ast }{)}^{T}r(x^{\ast })=0,$$所以若$r(x^{\ast })̸=0$就应有$J(x^{\ast })$为奇异.

由于本章中的算法可能会陷入平方和价值函数的局部极小, 因此我们这里讨论的全局收敛性要比无约束优化弱.
我们也可以使用其他的价值函数. 例如$l_1$范数价值函数:$$f_1(x)=\Vert r(x){\Vert }_1=\sum _{i=1}^m|r_i(x)|.$$我们在后续考虑约束优化时会着重讨论这一函数.

2.2 线搜索算法

我们可以将线搜索算法应用到平方和价值函数$f(x)=\frac{1}{2}\Vert r(x){\Vert }^2$上, 以得到算法的全局收敛性. 当迭代良定时, 牛顿步$p_k$:$$J(x_k)p_k=-r(x_k)$$在$r_k̸=0$时就是$f(\cdot )$的一个下降方向, 这是因为$$p_k^T\nabla f(x_k)=-p_k^T{J}_k^T{r}_k=-\Vert r_k{\Vert }^2<0.$$步长${\alpha }_k$按第三章中的某个条件选取. 最终迭代公式就是$$x_{k+1}=x_k+{\alpha }_k{p}_k,k=0,1,2,\dots .$$对于满足Wolfe条件的${\alpha }_k$, 我们有如下收敛结果. 这是Zoutendijk定理的直接推论.

定理6 设$J(x)$在水平集 L = { x : f ( x ) ≤ f ( x 0 ) } \mathcal{L}=\left\{x:f(x)\le f(x_0)\right\} L={x:f(x)≤f(x0​)}的某个邻域$\mathcal{D}$内Lipschitz连续, 且$\Vert J(x)\Vert ,\Vert r(x)\Vert$在$\mathcal{D}$上有界. 将线搜索算法应用于$f$, 其中搜索方向$p_k$满足$p_k^T\nabla f_k<0$, 步长${\alpha }_k$满足Wolfe条件. 于是Zoutendijk条件满足, 我们有$$\sum _{k\ge 0}{\cos }^2{\theta }_k\Vert J_k^T{r}_k{\Vert }^2<\infty ,$$其中$$\cos {\theta }_k=\frac{-p_k^T\nabla f(x_k)}{\Vert p_k\Vert \Vert \nabla f(x_k)\Vert }.$$

定理是证明只需我们验证$\nabla f$在$\mathcal{D}$上Lipschitz连续且$f$在$\mathcal{D}$上下有界(0)即可.

若迭代序列满足$$\cos {\theta }_k\ge \delta ,\delta \in (0,1),k充分大,$$则定理6告诉我们$J_k^T{r}_k\to 0$, 表明迭代点会趋近于价值函数$f$的稳定点. 进一步地, 若有$\Vert J(x_k{)}^{-1}\Vert$有界, 则就能得到$r_k\to 0$.

下面我们来讨论由牛顿法和非精确牛顿法产生的方向对应的$\cos {\theta }_k$.

1. 对于精确牛顿步, 我们有$$\cos {\theta }_k=-\frac{p_k^T\nabla f(x_k)}{\Vert p_k\Vert \Vert \nabla f(x_k)\Vert }=\frac{\Vert r_k{\Vert }^2}{\Vert J_k^{-1}{r}_k\Vert \Vert J_k^T{r}_k\Vert }\ge \frac{1}{\Vert J_k^T\Vert \Vert J_k^{-1}\Vert }=\frac{1}{\kappa (J_k)}.$$
2. 当$p_k$为非精确牛顿步, 即有$$\begin{array}{rl}\Vert r_k+J_k{p}_k{\Vert }^2\le {\eta }_k^2\Vert r_k{\Vert }^2& \Rightarrow 2p_k^T{J}_k^T{r}_k+\Vert r_k{\Vert }^2+\Vert J_k{p}_k{\Vert }^2\le {\eta }^2\Vert r_k{\Vert }^2\\ & \Rightarrow p_k^T\nabla f_k=p_k^T{J}_k^T{r}_k\le [({\eta }^2-1)/2]\Vert r_k{\Vert }^2.\end{array}$$与此同时, $$\Vert p_k\Vert \le \Vert J_k^{-1}\Vert [\Vert r_k+J_k{p}_k\Vert +\Vert r_k\Vert ]\le \Vert J_k^{-1}\Vert (\eta +1)\Vert r_k\Vert ,$$以及$$\Vert \nabla f_k\Vert =\Vert J_k^T{r}_k\Vert \le \Vert J_k\Vert \Vert r_k\Vert .$$结合以上估计, 我们有$$\cos {\theta }_k=-\frac{p_k^T\nabla f_k}{\Vert p_k\Vert \Vert \nabla f_k\Vert }\ge \frac{1-{\eta }^2}{2\Vert J_k\Vert \Vert J_k^{-1}\Vert (1+\eta )}\ge \frac{1-\eta }{2\kappa (J_k)}.$$

因此若条件数$\kappa (J_k)$有界, 则精确与非精确的牛顿法均是全局收敛的. 我们观察到这里牛顿法的收敛似乎变得与条件数相关了. 但注意到我们这里讨论的是全局收敛性, 而不是局部收敛性. 因此牛顿法的尺度不变性依然是成立的.

不过, 当$\kappa (J_k)$很大时, 这两个下界就会接近于0, 使用牛顿方向就会使算法产生不好的表现. 事实上, 下面的例子就表明$\cos {\theta }_k$是可以收敛到0的, 从而导致算法无法继续. 这个例子也说明了线搜索算法的一个基本的缺点.

例2 考虑求以下非线性系统的根:$$r(x)=\left[\begin{array}{c}{x}_1\\ \frac{10x_1}{(x_1+0.1)}+2x_2^2\end{array}\right].$$已知它有唯一解$x^{\ast }=0$. 我们以带精确线搜索的牛顿迭代求解该问题. 可以证明, 从$(3,1{)}^T$出发, 迭代点将收敛于$(1.8016,0{)}^T$. 但这并不是根, 甚至连$f$的稳定点都不是! 从这点出发沿着负梯度$-\nabla f$走将能使得$r$在两个分量上都有所下降. 下验证之. 注意$r$的Jacobi矩阵为$$J(x)=\left[\begin{array}{cc}1& 0\\ \frac{1}{(x_1+0.1{)}^2}& 4x_2\end{array}\right],$$它在$x_2=0$的点处奇异. 对这样的点, 我们有$$\nabla f(x)=\left[\begin{array}{c}{x}_1+\frac{10x_1}{(x_1+0.1{)}^3}\\ 0\end{array}\right],$$从而若$x_1>0$则梯度将指向$x_1$轴的正向. 因此$(1.8016,0{)}^T$不是$f$的稳定点.

对于这个例子, 计算可知从一个靠近(但不在)$x_1$轴的点出发的牛顿方向趋于平行于$x_2$轴, 差不多与梯度$\nabla f(x)$正交. 因此, $\cos {\theta }_k$将以任意精度接近0.

在这个例子中, 带精确线搜索的牛顿法陷于一个Jacobi矩阵奇异的点. 由于非线性系统常包含奇异点, 因此这样的现象是需要我们去关注的.

为防止这一现象发生, 保证$\cos {\theta }_k\ge \delta >0$, 我们需要修正牛顿方向.

一种可能的方式是, 在$J_k^T{J}_k$上加上某个纯量阵${\lambda }_{k}I$, 定义$p_k$为$$p_k=-(J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}{J}_k^T{r}_k.$$对任意${\lambda }_k>0$, 括号里的矩阵都是非奇异的, 且若${\lambda }_k$远离0值, 则就满足条件$\cos {\theta }_k\ge \delta >0,\exists \delta >0$. 事实上, 设$J_k^T{J}_k$有谱分解$$J_k^T{J}_k=Q\Lambda Q^T=\sum _{j=1}^n{\mu }_j{q}_j{q}_j^T,$$其中$Q$为正交阵, ${\mu }_j$是$J_k^T{J}_k$的特征值, 且满足${\mu }_1\ge {\mu }_2\ge \cdots \ge {\mu }_n\ge 0,{\mu }_1>0$. 利用不等式$$\begin{gathered} \Vert J_k^T{J}_k+{\lambda }_{k}I\Vert \ge \Vert J_k^T{J}_k\Vert , \\ \Vert (J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}\Vert \Vert J_k^T{J}_k+{\lambda }_{k}I\Vert =\kappa (J_k^T{J}_k+{\lambda }_{k}I)\le \delta ,\exists \delta >0, \end{gathered}$$我们有$$\begin{array}{rl}\cos {\theta }_k& =-\frac{p_k^T\nabla f_k}{\Vert p_k\Vert \Vert \nabla f_k\Vert }=\frac{r_k^T{J}_k(J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}{J}_k^T{r}_k}{\Vert J_k^T{r}_k\Vert \Vert (J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}{J}_k^T{r}_k\Vert }\\ & \ge \frac{r_k^T{J}_k(J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}{J}_k^T{r}_k}{\Vert J_k^T{r}_k\Vert \Vert (J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}\Vert \Vert J_k^T{r}_k\Vert }\\ & \ge \frac{r_k^T{J}_k(J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}{J}_k^T{r}_k\Vert J_k^T{J}_k\Vert }{\Vert J_k^T{r}_k{\Vert }^2\Vert (J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}\Vert \Vert J_k^T{J}_k+{\lambda }_{k}I\Vert }\\ & \ge \frac{1}{\delta }\frac{(r_k^T{J}_k)\left({\sum }_{j=1}^n\frac{{\mu }_1}{{\mu }_j+{\lambda }_k}{q}_j{q}_j^T\right)(J_k^T{r}_k)}{\Vert J_k^T{r}_k{\Vert }^2}\\ & \ge \frac{1}{\delta }\frac{\frac{{\mu }_1}{{\mu }_1+{\lambda }_k}(r_k^T{J}_k)Q{Q}^T(J_k^T{r}_k)}{\Vert J_k^T{r}_k{\Vert }^2}\\ & =\frac{{\mu }_1}{({\mu }_1+{\lambda }_k)\delta }.\end{array}$$因此只要对于所有$J_k^T{J}_k$, ${\mu }_1>0$且${\lambda }_k$远离0值, 就有全局收敛性. 因此, 一些使用算法自适应地选取${\lambda }_k$以保证矩阵的非奇异性(同时也不接近奇异). 这一方法类似于第十章的经典Levenberg-Marquardt算法. 为使算法免于显式计算$J_k^T{J}_k$并对矩阵$(J_k^T{J}_k+\lambda I)$做Cholesky分解, 我们可以使用先前在最小二乘情形下讨论的方法. 它主要利用的是, $(J_k^T{J}_k+\lambda I)$的Cholesky因子就等于$R^T$, 其中$R$为$$\left[\begin{array}{c}{J}_k\\ \sqrt{\lambda }I\end{array}\right]$$QR分解的上三角矩阵. 我们同样可以组合使用Householder和Givens变换, 从而在改变${\lambda }_k$的取值时减少计算量.

这一方法的缺陷在于, 对我们来说, ${\lambda }_k$的选取相当困难: 如果选的太大, 我们就会破坏牛顿法的快速收敛. 这是因为当${\lambda }_k\uparrow \infty$时趋向于$-J_k^T{r}_k$的倍数, 从而步长将变得很小, 搜索方向也倾向于变成$f$的最速下降方向. 反之, 如果选的太小, 算法就不能很好地处理Jacobi矩阵的奇异性.
一种更加直观、合理的方式, 是将其嵌入信赖域的框架. 它能间接地选取${\lambda }_k$. 我们将在下一小节讨论之.

下面我们具体地表述基于牛顿类的搜索方向和线搜索的算法, 它将在必要时给予步长计算以正则化. 其中一些细节留作读者思考或自学.

算法4 (Line Search Newton-like Method)
Given $c_1,c_2$ with $0<c_1<c_2<\frac{1}{2}$;
Choose $x_0$;
for $k=0,1,2,\dots$
\quad\quad Calculate a Newton-like step from $J(x_k)p_k=-r(x_k)$, (regularizing with $p_k=-(j_k^T{J}_k+{\lambda }_{k}I{)}^{-1}{J}_k^T{r}_k$if $J_k$ appears to be near-singular), or $\Vert r_k+J_k{p}_k\Vert \le {\eta }_k\Vert r_k\Vert$, or $p_k=-B_k^{-1}r(x_k)$;
\quad\quad if $\alpha =1$ satisfies the Wolfe conditions
\quad\quad\quad\quad Set ${\alpha }_k=1$;
\quad\quad else
\quad\quad\quad\quad Perform a line search to find ${\alpha }_k>0$ that satisfies Wolfe conditions;
\quad\quad end (if)
\quad\quad $x_{k+1}\leftarrow x_k+{\alpha }_k{p}_k$;
end (for)

2.3 信赖域算法

对于非线性方程组最常使用的信赖域算法就是直接将第四章的算法应用于价值函数$f(x)=\frac{1}{2}\Vert r(x){\Vert }_2^2$, 其中使用$B_k=J(x_k{)}^{T}J(x_k)$作为模型函数$m_k(p)$的近似Hessian矩阵, $m_k(p)$定义为$$m_k(p)=\frac{1}{2}\Vert r_k+J_k{p}_k{\Vert }_2^2=f_k+p^T{J}_k^T{r}_k+\frac{1}{2}{p}^T{J}_k^T{J}_k^T{p}_k.$$于是$p_k$就由近似求解下面的信赖域子问题获得:$$\underset{p}{\min }{m}_k(p),\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\Vert p\Vert \le {\Delta }_k,$$其中${\Delta }_k$为信赖域的半径. 实际下降与预测下降的比值${\rho }_k$在许多信赖域算法中扮演很关键的角色:$${\rho }_k=\frac{\Vert r(x_k){\Vert }^2-\Vert r(x_k+p_k){\Vert }^2}{\Vert r(x_k){\Vert }^2-\Vert r(x_k)+J(x_k)p_k{\Vert }^2},$$它将辅助我们决定是否接受新的迭代点以及如何处理信赖域的半径.

算法5 (Trust-Region Method for Nonlinear Equations)
Given $\bar{\Delta }>0,{\Delta }_0\in (0,\bar{\Delta })$, and $\eta \in [0,\frac{1}{4})$:
for $k=0,1,2,\dots$
\quad\quad Calculate $p_k$ as an (approximate) solution of the TR subproblem;
\quad\quad Evaluate ${\rho }_k$
\quad\quad if ${\rho }_k<\frac{1}{4}$
\quad\quad\quad\quad ${\Delta }_{k+1}=\frac{1}{4}\Vert p_k\Vert$;
\quad\quad else
\quad\quad\quad\quad if ${\rho }_k>\frac{3}{4}$ and $\Vert p_k\Vert ={\Delta }_k$
\quad\quad\quad\quad\quad\quad ${\Delta }_{k+1}=\min (2{\Delta }_k,\bar{\Delta })$;
\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad ${\Delta }_{k+1}={\Delta }_k$;
\quad\quad\quad\quad end (if)
\quad\quad end (if)
\quad\quad if ${\rho }_k>\eta$
\quad\quad\quad\quad $x_{k+1}=x_k+p_k$;
\quad\quad else
\quad\quad\quad\quad $x_{k+1}=x_k$;
\quad\quad end (if)
end (for).

折线算法是信赖预算法的一个特殊情形. 它基于Cauchy点$p_k^C$和$m_k(p)$的无约束极小点构建子问题的近似解. 这里Cauchy点为$$p_k^C=-{\tau }_k({\Delta }_k/\Vert J_k^T{r}_k\Vert )J_k^T{r}_k,$$其中 τ k = min ⁡ { 1 , ∥ J k T r k ∥ 3 / ( Δ k r k T J k ( J k T J k ) J k T r k ) } ; \tau_k=\min\{1,\Vert J_k^Tr_k\Vert^3/(\Delta_kr_k^TJ_k(J_k^TJ_k)J_k^Tr_k)\}; τk​=min{1,∥JkT​rk​∥3/(Δk​rkT​Jk​(JkT​Jk​)JkT​rk​)};与之前的一般情形比较后我们发现, 我们不需要考虑$m_k(p)$中会出现不定Hessian近似的情形. 因此$m_k(p)$的无约束极小点在$J_k$非奇异时是唯一的. 此解我们以$p_k^J$表示:$$p_k^J=-(J_k^T{J}_k{)}^{-1}(J_k^T{r}_k)=-J_k^{-1}{r}_k.$$折线算法中$p_k$的选取方式如下:

程序6 (Dogleg)
Calculate $p_k^C$;
if $\Vert p_k^C\Vert ={\Delta }_k$
\quad\quad $p_k\leftarrow p_k^C$;
else
\quad\quad Calculate $p_k^J$;
\quad\quad $p_k\leftarrow p_k^C+\tau (p_k^J-p_k^C)$,
\quad\quad where $\tau$ is the largest value in $[0,1]$ such that $\Vert p_k\Vert \le {\Delta }_k$;
end (if)

由第四章中的引理, 我们知道当$J_k$非奇异时, 上述程序中选出的$p_k$就是$m_k$沿着从原点途径Cauchy点再到无约束极小点$p_k^J$这一分段直线路径上的极小点. 因此, 模型函数上的下降至少不会低于Cauchy点带来的下降, 这在最小二乘情形下可写作$$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$是某个正常数.

我们还知道, 子问题的精确解具有形式$$p_k=-(J_k^T{J}_k+{\lambda }_{k}I{)}^{-1}{J}_k^T{r}_k,\exists {\lambda }_k\ge 0,$$并且当无约束解$p_k^J$满足$\Vert p_k^J\Vert \le {\Delta }_k$时可选取${\lambda }_k=0$. 另外注意上式恰好等于第十章中的$p_k^{\mathrm{L}\mathrm{M}}$. 事实上, 求解非线性方程组的Levenberg-Marquardt法就是求解非线性最小二乘问题同样算法的特殊情形. Levenberg-Marquardt法使用信赖域的迭代方法求${\lambda }_k$. 信赖域算法中的精确算法是基于Cholesky分解的, 但在第十章, 我们可代之以计算$$\left[\begin{array}{c}{J}_k\\ \sqrt{\lambda }I\end{array}\right]$$QR分解的算法. 即使求不到子问题的精确${\lambda }_k$, 由上式求得的$p_k$若对某个$c_1$满足下降条件, 且有$$\Vert p_k\Vert \le \gamma {\Delta }_k,\gamma \ge 1,$$则仍能使算法具有全局收敛性.

折线算法在每次迭代中仅需要求解一个线性系统, 而反观精确求解需要求解多个这样的线性系统. 但也如第四章说明的, 我们要权衡在每步迭代上的耗费和需要计算的函数值域导数值的总数.

我们也可以考虑基于不同价值函数的信赖域算法以及信赖域的不同定义. 例如, 一种基于$l_1$价值函数和$l_{\infty }$信赖域的算法源自于如下形式的子问题:$$\underset{p}{\min }\Vert J_{k}p+r_k{\Vert }_1\mathrm{s}.\mathrm{t}.\Vert p{\Vert }_{\infty }\le \Delta ,$$而这可用线性规划的方法重构和求解. 我们将在后面讨论非线性规划的S$l_1$QP和SLQP算法时谈及.

当$p_k$满足下降条件以及长度限制时, 算法5的全局收敛性可见下面的定理. 类似于第四章, 第一个结论针对$\eta =0$, 这样算法将接受所有能为价值函数$f_k$带来下降的搜索方向; 而第二个则针对$\eta >0$, 从而能得到更强的结论.

定理7 设在水平集 L = { x : f ( x ) ≤ f ( x 0 ) } \mathcal{L}=\{x:f(x)\le f(x_0)\} L={x:f(x)≤f(x0​)}的一个邻域$\mathcal{D}$内, $J(x)$Lipschitz连续, $\Vert J(x)\Vert$上有界. 另设所有子问题的近似解$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),$$以及长度限制$$\Vert p_k\Vert \le \gamma {\Delta }_k,\gamma \ge 1.$$则当$\eta =0$时, 算法5有$$\underset{k\to \infty }{\mathrm{liminf}}\Vert J_k^T{r}_k\Vert =0;$$而当$\eta \in (0,\frac{1}{4})$, 有$$\underset{k\to \infty }{\lim }\Vert J_k^T{r}_k\Vert =0.$$

下面我们转向考虑信赖域子问题在精确求解时算法的局部收敛性质. 我们假设序列 { x k } \{x_k\} {xk​}收敛到非线性方程组$r(x)=0$的一个非退化解$x^{\ast }$. 下面的结果, 意义在于在良好设计的算法中, 为达全局收敛而加强的条件并不干涉在第1节中叙述的局部快速收敛.

定理8 设由算法5产生的序列 { x k } \{x_k\} {xk​}收敛到$r(x)=0$的一个退化解$x^{\ast }$. 设$J(x)$在$x^{\ast }$的一个开邻域$\mathcal{D}$内Lipschitz连续, 且对于充分大的$k$, 信赖域子问题被精确求解. 于是序列 { x k } \{x_k\} {xk​}二次收敛于$x^{\ast }$.

证明: 我们证明存在指标$K$使得信赖域半径在第$K$步迭代后不再减小, 即${\Delta }_k\ge {\Delta }_K,\forall k\ge K$. 于是算法终将采纳标准的牛顿步, 从而由定理2得知有局部二次收敛.
令$p_k$为子问题的精确解. 注意若此步在信赖域内, 则它就是无约束牛顿步$-J_k^{-1}{r}_k$. 否则我们有$\Vert J_k^{-1}{r}_k\Vert >{\Delta }_k$, 而$\Vert p_k\Vert \le \Delta$. 因此, 不论哪种情形, 我们都有$$\Vert p_k\Vert \le \Vert J_k^{-1}{r}_k\Vert .$$考虑比值$rh{o}_k$. 直接由定义我们有$$|1-{\rho }_k|\le \frac{\left|\Vert r_k+J_k{p}_k{\Vert }^2-\Vert r(x_k+p_k){\Vert }^2\right|}{\Vert r(x_k){\Vert }^2-\Vert r(x_k)+J(x_k)p_k{\Vert }^2}.$$由定理1, 分子的第二项可以写成$$\Vert r(x_k+p_k){\Vert }^2=\Vert [r(x_k)+J(x_k)p_k]+w(x_k,x_k+p_k){\Vert }^2.$$由$J$的Lipschitz连续性(Lipschitz常数为${\beta }_L$), 我们有$$\begin{array}{rl}\Vert w(x_k,x_k+p_k)\Vert & \le {\int }_0^1\Vert J(x_k+t{p}_k)-J(x_k)\Vert \Vert p_k\Vert \mathrm{d}t\\ & \le {\int }_0^1{\beta }_L\Vert p_k{\Vert }^2\mathrm{d}t=({\beta }_L/2)\Vert p_k{\Vert }^2,\end{array}$$再由于$$\Vert r_k+J_k{p}_k\Vert \le \Vert r_k\Vert =f(x_k{)}^{1/2}$$我们有$$\begin{array}{rl}& \left|\Vert r_k+J_k{p}_k{\Vert }^2-\Vert r(x_k+p_k){\Vert }^2\right|\\ & \le 2\Vert r_k+J_k{p}_k\Vert \Vert w(x_k,x_k+p_k)\Vert +\Vert w(x_k,x_k+p_k){\Vert }^2\\ & \le f(x_k{)}^{1/2}{\beta }_L\Vert p_k{\Vert }^2+({\beta }_L/2{)}^2\Vert p_k{\Vert }^4\\ & \le ϵ(x_k)\Vert p_k{\Vert }^2,\end{array}$$其中我们定义$$ϵ(x_k)=f(x_k{)}^{1/2}{\beta }_L+({\beta }_L/2{)}^2\Vert p_k{\Vert }^2.$$由于$x_k\to x^{\ast }$, 因此$f(x_k)\to 0,\Vert r_k\Vert \to 0$. 由于$x^{\ast }$为非退化解, 所以对充分大的$k$我们有$\Vert J(x_k{)}^{-1}\Vert \le {\beta }^{\ast }$. 所以$$\Vert p_k\Vert \le \Vert J_k^{-1}{r}_k\Vert \le {\beta }^{\ast }\Vert r_k\Vert \to 0.$$因此, $ϵ(x_k)\to 0$.
下面考察分母. 定义${\bar{p}}_k$为与牛顿步$-J_k^{-1}{r}_k$同向且与$p_k$同长的搜索方向, 即$${\bar{p}}_k=-\frac{\Vert p_k\Vert }{\Vert J_k^{-1}{r}_k\Vert }{J}_k^{-1}{r}_k.$$由于${\bar{p}}_k$为子问题可行点, 所以$$\begin{array}{rl}\Vert r_k{\Vert }^2-\Vert r_k+J_k{p}_k{\Vert }^2& \ge \Vert r_k{\Vert }^2-{\Vert r_k-\frac{\Vert p_k\Vert }{\Vert J_k^{-1}{r}_k\Vert }{r}_k\Vert }^2\\ & =2\frac{\Vert p_k\Vert }{\Vert J_k^{-1}{r}_k\Vert }\Vert r_k{\Vert }^2-\frac{\Vert p_k{\Vert }^2}{\Vert J_k^{-1}{r}_k{\Vert }^2}\Vert r_k{\Vert }^2\\ & \ge \frac{\Vert p_k\Vert }{\Vert J_k^{-1}{r}_k\Vert }\Vert r_k{\Vert }^2\ge \frac{1}{{\beta }^{\ast }}\Vert p_k\Vert \Vert r_k\Vert .\end{array}$$代入就有$$|1-{\rho }_k|\le \frac{{\beta }^{\ast }ϵ(x_k)\Vert p_k{\Vert }^2}{\Vert p_k\Vert \Vert r_k\Vert }\le ({\beta }^{\ast }{)}^2ϵ(x_k)\to 0.$$因此对于充分大的$k$我们有${\rho }_k>\frac{1}{4}$, 从而信赖域半径${\Delta }_k$在此之后不会再减小. 因此存在指标$K$使得$${\Delta }_k\ge {\Delta }_K,\forall k\ge K.$$由于$\Vert J_k^{-1}{r}_k\Vert \le {\beta }^{\ast }\Vert r_k\Vert \to 0$, 因此牛顿步$-J_k^{-1}{r}_k$的长度终将小于${\Delta }_K$(以及${\Delta }_k$), 从而最终成为子问题的解. 二次收敛性则由定理2推得. 证毕.

我们亦可以将$x_k\to x^{\ast }$的假设代换为非退化解$x^{\ast }$仅为序列的一个聚点. 事实上由Cauchy序列的性质, 这样也有$x_k\to x^{\ast }$.

3. 连续/同伦方法

3.1 动机

以上提及的基于牛顿法的算法均有一个缺点: 除非$J(x)$在指定区域内非奇异——这一条件往往不能保证——否则它们就可能会收敛到价值函数的局部极小点, 而不是非线性系统的解. 而本节中我们要介绍的连续方法则更倾向于在这种麻烦的情形下收敛到$r(x)=0$的解. 它们潜在的动机非常简单: 我们并不直接处理原始问题$r(x)=0$, 而是代之以某个"简单"的、易求解的方程组. 之后我们再将这一简单的系统逐渐转换成原始的系统$r(x)$, 同时跟踪解的动向: 从简单问题的解移动到原问题的解.

一种定义所谓同伦(homotopy)映射$H(x,\lambda )$的方式如下:$$H(x,\lambda )=\lambda r(x)+(1-\lambda )(x-a),$$其中$\lambda$为一标量参数, $a\in {\mathbb{R}}^n$为固定的向量. 当$\lambda =0$时, 上式就人为地定义了一个简单系统$H(x,0)=x-a$, 其解显然是$x=a$. 当$\lambda =1$, 我们就有原系统$H(x,1)=r(x)$.

为求解$r(x)=0$, 考虑如下算法: 首先, 置$\lambda =0$并置$x=a$. 之后将$\lambda$从0逐渐以小量增加到1, 并对每个$\lambda$的取值, 我们都计算相应系统$H(x,\lambda )=0$的解. 最终对应于$\lambda =1$的$x$就是原问题的$r(x)=0$的解.
这一简单的方法听上去似乎可行, 下图也阐释了一种算法可以顺利进行的情形.

图中, 对于每个$\lambda$的取值, $H(x,\lambda )=0$都只有唯一的解. 满足$H(x,\lambda )=0$的$(x,\lambda )$构成的轨迹称为零路径(zero paths).

然而, 这样的方法往往会失败, 如下图所示.

此时, 算法在下分支从$\lambda =0$移动到$\lambda ={\lambda }_T$, 之后除非它足够幸运, 能够跳到路径的上分支, 否则算法就会失效. ${\lambda }_T$称作是拐点(turning point), 这是因为在这点处我们只有减小$\lambda$才能继续跟随路径. 事实上, 实用的连续方法将考虑上图中的情形, 即他们会显式地跟进零路径, 即使这意味着需要有时减小$\lambda$.

3.2 实用连续方法

3.2.1 方法一

在一种实用的方法中, 我们以弧长参数(arc length parameter)$s$对零路径参数化. 即$(x(s),\lambda (s))$为我们从初始点$(x(0),\lambda (0))=(a,0)$出发走过长度$s$后所达到的点. 由于$$H(x(s),\lambda (s))=0,\forall s\ge 0,$$因此我们可对$s$求导:$$\frac{\partial }{\partial x}H(x,\lambda )\dot{x}+\frac{\partial }{\partial \lambda }H(x,\lambda )\dot{\lambda }=0,(\dot{x},\dot{\lambda })=\left(\frac{\mathrm{d}x}{\mathrm{d}s},\frac{\mathrm{d}\lambda }{\mathrm{d}s}\right).$$向量$(\dot{x}(s),\dot{\lambda }(s))$为零路径的切向量, 如3.1节的第二张图所示. 从上式可得, 切向量位于$n\times (n+1)$矩阵$$\left[\begin{array}{cc}\frac{\partial }{\partial x}H(x,\lambda )& \frac{\partial }{\partial \lambda }H(x,\lambda )\end{array}\right]$$的零空间. 当此矩阵(行)满秩, 它的零空间就是1维的. 为了在这种情形下完善$(\dot{x},\dot{\lambda })$的定义, 我们需要赋给它长度和方向. 长度可用归一化条件决定:$$\Vert \dot{x}(s){\Vert }^2+|\dot{\lambda }(s){|}^2=1,\forall s,$$这保证了$s$是从$(0,a)$到$(x(s),\lambda (s))$的弧长参数. 我们还需要确定符号以保证我们是沿着零路径的正向移动的. 一种比较好的启发式方法是选取符号, 使得当前$s$的切向量$(\dot{x},\dot{\lambda })$与先前$s$的切向量夹角小于$\pi /2$.

我们将计算$(\dot{x},\dot{\lambda })$的程序概述如下:

程序7 (Tangent Vector Calculation)
Compute a vector in the null space by performing a QR factorization with column pivoting,$$Q^T\left[\begin{array}{cc}\frac{\partial }{\partial x}H(x,\lambda )& \frac{\partial }{\partial \lambda }H(x,\lambda )\end{array}\right]\Pi =\left[\begin{array}{cc}R& w\end{array}\right],$$where $Q$ is $n\times n$ orthogonal, $R$ is $n\times n$ upper triangular, $\Pi$ is an $(n+1)\times (n+1)$ permutation matrix, and $w\in {\mathbb{R}}^n$.
Set $$v=\Pi \left[\begin{array}{c}{R}^{-1}w\\ -1\end{array}\right];$$Set $(\dot{x},\dot{\lambda })=\pm v/\Vert v{\Vert }_2$, where the sign is chosen to satisfy the angle criterion mentioned above.

由于我们有在任意给定点$(x,\lambda )$处的切向量以及初始点$(x(0),\lambda (0))=(a,0)$, 因此我们可调用标准一阶常微分方程初值问题求解器得到零路径, 求解器将在它找到使得$\lambda (s)=1$的$s$后终止.

3.2.2 方法二

第二种跟进零路径的方式与刚才介绍的很相似, 只不过它是从代数角度出发而不是从微分方程的角度. 给定当前点$(x,\lambda )$, 我们依照上面的方式计算切向量$(\dot{x},\dot{\lambda })$, 并沿着这一方向走一小步(长度为较小的$ϵ$)得到一个"预估"(predictor)点$(x^p,{\lambda }^p)$, 即$$(x^p,{\lambda }^p)=(x,\lambda )+ϵ(\dot{x},\dot{\lambda }).$$通常这一预估点并不精确地落在零路径上, 从而我们再需要一些"校正"(corrector)步将它带回路径上, 从而最终得到满足$H(x^{+},{\lambda }^{+})=0$的新点$(x^{+},{\lambda }^{+})$. 这一过程如下图所示.

在校正步中, 我们将预估点$(x^p,{\lambda }^p)$相对原点$(x,\lambda )$变化最大的那个分量固定住. 若这一分量的指标为$i$, 且若我们在校正步中使用标准的牛顿法(这通常是适定的, 因为$(x^p,{\lambda }^p)$一般与目标点$(x^{+},{\lambda }^{+})$离的很近), 则校正步就有形式:$$\left[\begin{array}{c}\frac{\partial H}{\partial x}\mid \frac{\partial H}{\partial \lambda }\\ e_i\end{array}\right]\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]=\left[\begin{array}{c}-H\\ 0\end{array}\right],$$其中$\partial H/\partial x,\partial H/\partial \lambda ,H$均在校正步之前最近的迭代点处取值. 这一系统的最后一行是要保证$(\delta x,\delta \lambda )$的第$i$个分量为0. 例如在上图中, 我们固定了$\lambda$的值. 但在下一步迭代, 我们可能就需要固定$x$. (这是因为快到到达$\lambda$的拐点了.)

3.2.3 总结

以上两种路径跟随的变体都能够跟随具有拐点的曲线. 不过, 它们需要$n\times (n+1)$矩阵对于所有的$(x,\lambda )$行满秩, 从而切向量才是良定的. 下面的结论表明, 行满秩在一定条件下是成立的.

定理9 设$r$二次连续可微. 于是对几乎所有的$a\in {\mathbb{R}}^n$, 存在从$(a,0)$出发的零路径, 且路径上$n\times (n+1)$系数矩阵有行满秩. 若这一路径对$\lambda \in [0,1)$有界, 则它就有聚点$(\bar{x},1)$使得$r(\bar{x})=0$. 进一步地, 若Jacobi矩阵$J(\bar{x})$非奇异, 则$(a,0)$与$(\bar{x},1)$之间的零路径可求长.

这一定理表明, 除非我们选到了不好的$a$, 以上的算法才可能会得到发散的路径. 我们以一个简单的例子说明零路径的发散是可能发生的.

例3 考虑$r(x)=x^2-1$, 其有两个非退化解+1和-1. 选取$a=-2$并尝试应用连续方法于函数$$H(x,\lambda )=\lambda (x^2-1)+(1-\lambda )(x+2)=\lambda x^2+(1-\lambda )x+(2-3\lambda ).$$这一函数的零路径如下图所示.

从图中可以看到在$(-2,0)$和$(1,1)或者(-1,1)$之间不存在任何的零路径, 从而此例连续方法会失效. 我们可求得对$\lambda$的取值, $x$取什么范围内的值时连续方法会失效(使用二次求根公式):$$x=\frac{-(1-\lambda )\pm \sqrt{(1-\lambda {)}^2-4\lambda (2-3\lambda )}}{2\lambda }.$$当根号里为负, 意味着没有$x$的实根. 此时$$\lambda \in \left(\frac{5-2\sqrt{3}}{13},\frac{5+2\sqrt{3}}{13}\right)\approx (0.118,0.651).$$这就解释了为什么从$(-2,0)$出发的零路径会无界.

这个例子说明, 即使是对一个相当简单的非线性系统, 连续方法依然可能会失效. 然而一般来说, 它们要比本章先前介绍的基于价值函数的方法要可靠得多. 不过这鲁棒性是基于一定的代价的: 连续方法往往比基于价值函数的方法需要过多的计算量. 连续方法也将成为我们求解约束优化问题的内点法的动机.