Numerical Optimization Ch15. Fundamentals of Algorithms for Nonlinear Constrained Optimization

第十五章: 非线性约束优化算法基础

本章, 我们开启对求解一般约束优化问题 $$\underset{x\in {\mathbb{R}}^n}{\min }f(x),\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\begin{array}{ll}{c}_i(x)=0,& i\in \mathcal{E},\\ c_i(x)\ge 0,& i\in \mathcal{I}\end{array}$$算法的讨论, 这里目标函数 $f$、约束函数 $c_i$均是 ${\mathbb{R}}^n$子集上的 光滑实值函数, $\mathcal{I},\mathcal{E}$分别是不等式和等式约束的有限指标集. 在 第十二章, 我们使用这一一般形式推导了问题的 最优性条件. 这是本书剩下各类算法的 核心动机. 它们尽管各自不同, 但本质上都是 迭代方法: 它们都是产生解 $x^{\ast }$的一串估计序列, 而我们期望序列能收敛到解. 在某些情形下, 它们也产生对Lagrange乘子的估计序列. 如同在无约束优化中, 我们仅考虑求问题局部解的算法.
本章我们并不纠结于某个算法的推导, 而是企图建立基本的概念、搭建普遍的框架. 在阅读完第1,2节后, 读者可以一瞥第3-6节中的材料与示例. 在接下来几个章节的学习中, 读者可以回过头来利用这些例子进行学习.

1. 优化算法分类

我们现在罗列本书剩下部分的算法. 对于非线性优化算法, 现在还没有标准的分类方式. 但在接下来的章节中, 我们按如下方式分类归总讲述:

1. 第十六章中我们讨论求解二次规划quadratic programming问题的算法. 我们将这一类算法单独拎出来主要是因为其本身相当重要:
   • 二次规划问题本身具有独特的性质, 基于此可以设计高效的算法;
   • 在逐步二次规划算法和某些非线性规划内点法中我们需要求解二次规划子问题, 这为后续的讨论作铺垫.
     在这里, 我们探讨积极集法、内点法和梯度投影法.
2. 第十七章中我们讨论罚函数和增广Lagrange函数(penalty and augmented Lagrangian)法.
   • 结合目标函数和约束函数, 我们可以构建罚函数(penalty function)并以求解一系列无约束问题代替求解原本的约束问题. 例如, 若原始问题中仅有等式约束, 我们就可以定义二次罚函数为$$f(x)+\frac{\mu }{2}\sum _{i\in \mathcal{E}}{c}_i^2(x),$$其中$\mu >0$为惩罚因子(penalty paramater). 我们将对不断增大的$\mu$极小化这一函数直至认定我们(以一定精度)获取了约束问题的解.
   • 使用精确罚函数(exact penalty function), 我们就可以仅通过求解单个无约束问题得到原始约束问题的局部解. 例如对等式约束问题, 我们定义函数$$f(x)+\mu \sum _{i\mathbb{E}}|c_i(x)|,\mu >0充分大$$就是一个精确罚函数. 尽管通常精确罚函数不可微, 但精确罚函数可用一系列光滑的子问题逼近得到极小.
   • 在增广Lagrange函数法中, 我们定义一个结合了Lagrange函数和二次惩罚函数性质的函数, 称作增广Lagrange函数. 对于等式约束问题, 它具有以下形式:$${\mathcal{L}}_A(x,\lambda ;\mu )=f(x)-\sum _{i\in \mathcal{E}}{\lambda }_i{c}_i(x)+\frac{\mu }{2}\sum _{i\in \mathcal{E}}{c}_i^2(x).$$基于此函数的算法将先固定$\lambda$为最优Lagrange乘子的近似、$\mu$为某个正值, 求得${\mathcal{L}}_A(\cdot ,\lambda ;\mu )$的近似极小点$x$. 之后在这一新的迭代点, $\lambda ,\mu$可能会被更新, 接着重复以上过程. 此法能避免二次惩罚函数带来的一些问题.
3. 第十八章中我们介绍逐步二次规划sequential quadratic programming, SQP算法. 本质上, 它们在每个迭代点都对问题构建一个二次规划子问题并定义搜索方向为此子问题的解. 基本的SQP算法中, 我们定义$(x_k,{\lambda }_k)$处的搜索方向$p_k$为以下二次规划的解:$$\begin{array}{rl}{\min }_p& \frac{1}{2}{p}^T{\nabla }_{xx}^2\mathcal{L}(x_k,{\lambda }_k)p+\nabla f(x_k{)}^{T}p\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& \nabla c_i(x_k{)}^{T}p+c_i(x_k)=0,i\in \mathcal{E},\\ & \nabla c_i(x_k{)}^{T}p+c_i(x_k)\ge 0,i\in \mathcal{I},\end{array}$$这里$\mathcal{L}$为Lagrange函数. 子问题中的目标函数其实就是对Lagrange函数从$x_k$到$x_k+p$改变量的一个近似(少了约束的导数), 而这里的约束就是对约束问题中约束的线性化. 我们可以加上一个信赖域约束控制$p$的长度和质量; 拟牛顿近似Hessian阵可用来代替${\nabla }_{xx}^2\mathcal{L}(x_k,{\lambda }_k)$. 有一变体称为是逐步线性-二次规划(sequential linear-quadratic programming), 其$p_k$的计算要经两个阶段. 首先, 求解忽略二次项的线性规划, 此时要加上信赖域约束. 接着, 利用积极集法求解等式约束子问题, 其中积极集为在上一步线性规划解处起作用约束的指标集.
4. 第十九章, 我们讨论非线性规划的内点法(interior-point methods for nonlinear programming). 这些方法可视作第十四章线性规划原始-对偶内点法的延伸. 我们也可以将它们视作障碍函数法(barrier methods), 其中搜索方向通过求解以下以障碍函数作为目标函数的问题获得(这里以对数障碍函数为例):
   $$\begin{array}{rl}\underset{x,s}{\min }& f(x)-\mu \sum _{i=1}^m\log s_i\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& c_i(x)=0,i\in \mathcal{E},\\ & c_i(x)-s_i=0,i\in \mathcal{I},\end{array}$$其中障碍因子(barrier parameter)$\mu >0$, 变量$s_i>0$为松弛变量. 内点法是时下非线性规划的最新研究热点, 并已成为逐步二次规划算法的重要竞争对手.

在第1,3,4类中的算法均使用了消元技术(elimination techniques), 即利用约束消去了问题中的自由度. 为方便后续, 我们在第3节讨论消元. 而在之后的小节中, 我们介绍价值函数(merit functions)和滤子(filters). 它们对于非线性规划算法的全局收敛性具有重要意义.

2. 不等式约束问题的组合困难

求解非线性规划问题的主要困难之一在于处理不等式约束, 特别地, 是决定这些约束在解处是否起作用. 积极集法的本质方法就是首先对最优积极集${\mathcal{A}}^{\ast }$——在解处满足等式的约束指标——作一个估计. 我们称这一估计为工作集(working set), 记为$\mathcal{W}$. 接着利用$\mathcal{W}$求解相应的等式约束问题, 检验对于$\mathcal{W}$的解的Lagrange乘子是否满足KKT条件. 如果满足, 则接受$x^{\ast }$作为原始约束问题的局部解; 不然, 改变工作集重复以上步骤. 此法基于"求解等式约束问题要比求解一般非线性规划简单得多"这一事实.

但是工作集$\mathcal{W}$的选择可能会很大——至多到$2^{|\mathcal{I}|}$, 其中$|\mathcal{I}|$为不等式约束的数量. 因此这一上界随着不等式数量而指数增长——这种现象我们称为是非线性规划的组合困难——我们不能设计出考虑所有$\mathcal{W}$选择的实用算法.
下面的例子表明即使不等式约束较少, 决定最优积极集也不是件容易的事.

例1 考虑问题$$\begin{array}{rl}\underset{x,y}{\min }& f(x,y)\stackrel{def}{=}\frac{1}{2}(x-2{)}^2+\frac{1}{2}(y-\frac{1}{2}{)}^2\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& (x+1{)}^{-1}-y-\frac{1}{4}\ge 0,\\ & x\ge 0,\\ & y\ge 0.\end{array}$$我们依次给约束标号1-3. 下图中的虚线展示了目标函数的等高线. 可行域由曲线和两个轴封闭而成.

由图中可知, 在解$(x^{\ast },y^{\ast }{)}^T=(1.953,0.089{)}^T$处起作用的只有第一个约束.
下面我们将工作集方法应用于此, 考虑$\mathcal{W}$的所有$2^3=8$种可能的选择.

1. 解处无约束积极, 即$\mathcal{W}=\varnothing$. 由于$\nabla f=(x-2,y-1/2{)}^T$, 可见$f$的无约束极小点位于可行域外面. 因此最优积极集不可能为空.
2. 解处三个约束都积极, 即 W = { 1 , 2 , 3 } \mathcal{W}=\{1,2,3\} W={1,2,3}. 由图中可得这是不可能的, 因为这三个约束的图像并不交于一点.
   还剩下6种可能性, 我们仅详细讨论其中的3种.
3. W = { 2 } \mathcal{W}=\{2\} W={2}, 即必有$x=0$. 据此极小化$f$得到$(0,1/2{)}^T$. 易验证此点不满足KKT条件.
4. W = { 1 , 3 } \mathcal{W}=\{1,3\} W={1,3}, 则只有单个可行点$(3,0{)}^T$. 因约束2不积极, 因此${\lambda }_2=0$. 求解可得${\lambda }_1=-16,{\lambda }_3=-16.5$. 因此不满足必要性条件.
5. W = { 1 } \mathcal{W}=\{1\} W={1}. 即得$(x,y{)}^T=(1.953,0.089{)}^T$, 此时${\lambda }_1=0.411,{\lambda }_2={\lambda }_3=0$. 满足KKT条件. 进一步可验证此点处二阶充分性条件也成立.

对于这一个小例子, 我们也能发现考虑所有可能的工作集选择相当耗时耗力. 但上图表明, 某些$\mathcal{W}$的选择可用关于函数和导数的信息去除. 事实上, 第十六章的积极集法就是利用这样的信息对最优积极集作一系列估计, 避免显然不是最优积极集的$\mathcal{W}$.

一种不同的方式是第十九章中的内点(或障碍函数)法. 这些方法产生远离由不等式约束界定的可行域边界的迭代点. 随着算法逐渐趋近解, 障碍的效果将会减弱, 从而保证估计的精度. 这样, 内点法就避免了非线性规划的组合困难.

3. 消元

在处理带约束的优化问题时, 一个很自然的想法是用约束减少问题的变量个数、降低自由度, 得到更简单的问题. 不过, 消元技术必须小心使用, 因为它们可能会改变会改变问题或者让问题变得病态.

我们先以一个简单的例子开头. 对它来说, 消元是安全且便捷的. $$\min f(x)=f(x_1,x_2,x_3,x_4),\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\begin{array}{l}{x}_1+x_3^2-x_4{x}_3=0,\\ -x_2+x_4+x_3^2=0.\end{array}$$我们可以放心大胆地利用约束, 令$$x_1=x_4{x}_3-x_3^2,x_2=x_4+x_3^2,$$从而得到仅有两个变量的$$h(x_3,x_4)=f(x_4{x}_3-x_3^2,x_4+x_3^2,x_3,x_4),$$这里就可以用无约束优化的算法.

非线性消元的危害则以下例说明.

例2 考虑问题$$\min x^2+y^2,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}(x-1{)}^3=y^2.$$目标函数和约束的等高线可见下图, 易得解为$(x,y)=(1,0)$.

利用约束消去$y$, 得到$$h(x)=x^2+(x-1{)}^3.$$显然, $h(x)\to -\infty (x\to -\infty )$. 如此看来, 盲目地变换会使我们得出问题无界的错误结论. 事实上, 约束$(x-1{)}^3=y^2$隐含$x\ge 1$的条件, 而这在解处是积极的. 因此若要消去$y$, 我们应当再显式地添加界约束$x\ge 1$进去.

此例表明使用非线性方程组消元可能导致不可逆转的后果. 基于此, 大多数优化算法并不使用非线性消元, 而是先对约束做线性化再对简化的问题消元. 下面我们系统地介绍下如何用线性约束消元.

3.1 使用线性约束的简单消元

现考虑带线性等式约束的非线性目标极小化问题:$$\min f(x),\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Ax=b,$$其中$A$为$m\times n:m\le n$矩阵. 不妨设$A$行满秩. 这样就可找到$A$的$m$列线性无关. 若将这些列组合起来成$m\times m$矩阵$B$, 并定义$n\times n$的$P$为将这些列调换至$A$的前$m$列的置换阵, 则有$$AP=[B\mid N],$$这里$N$表示$A$的剩下$n-m$列. 对应地我们定义子向量$x_B\in {\mathbb{R}}^m,x_N\in {\mathbb{R}}^{n-m}$为$$\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right]=P^{T}x,$$称$x_B$为基变量, $B$为基矩阵. 注意$P{P}^T=I$, 于是我们可以将约束重新写作$$b=Ax=AP(P^{T}x)=B{x}_B+N{x}_N.$$整理可得基变量的表示:$$x_B=B^{-1}b-B^{-1}N{x}_N.$$我们将上式称为简单消元(simple elimination of variables). 这样对任意$x_N$的取值, 都有相对于约束的可行点. 于是问题等价于无约束问题$$\underset{x_N}{\min }h(x_N)\stackrel{def}{=}f\left(P\left[\begin{array}{c}{B}^{-1}b-B^{-1}N{x}_N\\ x_N\end{array}\right]\right).$$以上说明, 带线性等式约束的非线性优化问题(在数学角度上)等价于一个无约束问题.

例3 考虑问题$$\begin{array}{rl}\min & \sin (x_1+x_2)+x_3^2+\frac{1}{3}(x_4+x_5^4+x_6/2)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}& 8x_1-6x_2+x_3+9x_4+4x_5=6,\\ & 3x_1+2x_2-x_4+6x_5+4x_6=-4.\end{array}$$定义$P$为置换阵, 使得$$Px=\left[\begin{array}{c}{x}_3\\ x_6\\ x_1\\ x_2\\ x_4\\ x_5\end{array}\right],AP=\left[\begin{array}{ccccccc}1& 0& \mid & 8& -6& 9& 4\\ 0& 4& \mid & 3& 2& -1& 6\end{array}\right].$$此时基矩阵$B$为对角阵, 易于求逆. 从简单消元公式可得$$\left[\begin{array}{c}{x}_3\\ x_6\end{array}\right]=-\left[\begin{array}{cccc}8& -6& 9& 4\\ \frac{3}{4}& \frac{1}{2}& -\frac{1}{4}& \frac{3}{2}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_4\\ x_5\end{array}\right]+\left[\begin{array}{c}6\\ -1\end{array}\right].$$消去原问题的$x_3,x_6$, 问题变成$$\begin{array}{rl}\underset{x_1,x_2,x_4,x_5}{\min }& \sin (x_1+x_2)+(8x_1-6x_2+9x_4+4x_5-6{)}^2\\ & +\frac{1}{3}(x_4+x_5^4-[(1/2)+(3/8)x_1+(1/4)x_2-(1/8)x_4+(3/4)x_5]).\end{array}$$我们也可以选$A$的其他两列做基, 但这样$B^{-1}N$就不会那么好算了.

一般, 我们可通过Gauss消去法选取$m$个线性无关列. 用线性代数的语言说, 我们可以计算矩阵的行阶梯型, 选取主元所在列构作基矩阵$B$. 我们自然希望$B$容易分解且良态. 一种满足这些需求的技术是稀疏Gauss消去法. 它会在控制舍入误差的同时保留稀疏性. 这一算法的实施可用HSL库中的MA48. 但下面我们会说明, Gauss消去过程不能保证能够选取最佳的基矩阵.

简单消元法有一个有趣的解释. 为简化记号, 假设基矩阵就是$A$的前$m$列, 即$P=I$.
从以上, 我们知道任意可行点均有形式$$\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right]=x=Yb+Z{x}_N,$$其中$$Y=\left[\begin{array}{c}{B}^{-1}\\ O\end{array}\right],Z=\left[\begin{array}{c}-B^{-1}N\\ I\end{array}\right].$$注意$Z$的列向量线性无关且$AZ=O$, 因此$Z$的列向量张成$A$的核空间. 另外, $Y,Z$的列构成了一个线性无关组; $Yb$为线性约束$Ax=b$的特解.

换句话说, 简单消元法将可行点表示成$Ax=b$的特解加上沿着$A$核空间的位移. 特解$Yb$有时被称作坐标松弛步(coordinate relaxation step). 下图表示, 我们选取基矩阵$B$为$A$的第一列, 从而坐标松弛步位于横轴; 若选取$B$为$A$的第二列, 则坐标松弛步位于纵轴.

简单消元成本较低但可能会带来数值上的不稳定. 若上图中的可行集为一条几近与$x_1$轴平行的直线, 则在横轴上的坐标松弛步就会很长. 计算$x$时就会要计算两个相当大的向量的差. 此时, 选取位于$x_2$轴的坐标松弛步就更好. 因此选取最佳基并不是很直观. 为避免出现过长的坐标松弛步, 我们可以定义特解$Yb$为约束下的最小范数步. 这种方法是我们将要介绍的更广泛的消元法的特殊情形.

3.2 利用线性约束的广义消元策略

为推广上一小节的内容, 我们选取$Y\in {\mathbb{R}}^{n\times m},Z\in {\mathbb{R}}^{n\times (n-m)}$:$$[Y\mid Z]\in {\mathbb{R}}^{n\times n}非奇异,AZ=O.$$这些性质表明, $Z$的列向量是$A$核空间的一组基. 由于$A$行满秩, 于是$A[Y\mid Z]=[AY\mid O]$也行满秩, 这就得出$AY$非奇异. 我们将$Ax=b$的任一解用以下式子表示:$$x=Y{x}_Y+Z{x}_Z,$$其中$x_Y\in {\mathbb{R}}^m,x_Z\in {\mathbb{R}}^{n-m}$. 代换进约束可得$$Ax=(AY)x_Y=b;$$因此由$AY$的非奇异性, $x_Y$可被显式写作$$x_Y=(AY{)}^{-1}b.$$代入$x$的表达式, 可知对任一$x_Z\in {\mathbb{R}}^{n-m}$, $$x=Y(AY{)}^{-1}b+Z{x}_Z$$均满足约束$Ax=b$. 因此, 原约束问题等价与无约束问题$$\underset{x_Z}{\min }f(Y(AY{)}^{-1}b+Z{x}_Z).$$我们自然希望$AY$尽可能良态. 我们可以用$A^T$的QR分解计算$Y,Z$:$$A^T\Pi =\left[\begin{array}{cc}{Q}_1& Q_2\end{array}\right]\left[\begin{array}{c}R\\ O\end{array}\right],$$这里$\left[\begin{array}{cc}{Q}_1& Q_2\end{array}\right]$正交, $Q_1,Q_2$分别是$n\times m,n\times (n-m)$的子阵, $R$为$m\times m$的上三角阵, $Pi$为$m\times m$的置换阵. 定义$$Y=Q_1,Z=Q_2,$$所以$Y,Z$的列就构成了${\mathbb{R}}^n$的一组正交基. 进一步整理可得$$AY=\Pi R^T,AZ=O.$$此时$Y,Z$满足条件, 且$AY$的条件数与$R$相同, 从而与$A$相同. 利用这个表示, 我们有$Ax=b$的任一解均可被表示成$$x=Q_1{R}^{-T}{\Pi }^{T}b+Q_2{x}_Z.$$这里$R^{-T}{\Pi }^{T}b$的计算成本也并不高.

注意特解$Q_1{R}^{-T}{\Pi }^{T}b$也可以写作$A^T(A{A}^T{)}^{-1}b$. 这正是下面问题的解:$$\min \Vert x{\Vert }^2,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Ax=b;$$即是$Ax=b$的最小范数解. 如下图所示.

使用正交基的消元法从数值稳定性上讲是理想的. 而这一消元策略的主要计算量在于QR分解. 不幸的是, 对于$A$较大且稀疏的问题, 稀疏QR分解要比稀疏Gauss分解成本高得多. 因此也有人提出其他的折中技术.

仍假设$P=I$. 定义$$Y=\left[\begin{array}{c}I\\ (B^{-1}N{)}^T\end{array}\right],Z=\left[\begin{array}{c}-B^{-1}N\\ I\end{array}\right].$$此时$Y,Z$的列范数不再为1, 且$AZ=O,Y^{T}Z=O$. 这说明这样的$Y,Z$是合适的. 进一步, 我们可以证明这里的$Y(AY{)}^{-1}b$也是$Ax=b$的最小范数解. 事实上$$\begin{array}{rl}Y(AY{)}^{-1}& =\left[\begin{array}{c}I\\ (B^{-1}N{)}^T\end{array}\right](B+N{N}^T{B}^{-T}{)}^{-1}\\ & =\left[\begin{array}{c}I\\ (B^{-1}N{)}^T\end{array}\right]{\left(A{A}^T{B}^{-T}\right)}^{-1}\\ & =\left[\begin{array}{c}{B}^T\\ N^T\end{array}\right](A{A}^T{)}^{-1}=A(A{A}^T{)}^{-1}.\end{array}$$这表明$Y(AY{)}^{-1}$与$B$的选取时无关的, 且其条件数仅由$A$决定. 不过仍要注意, 矩阵$Z$仍然显式依赖于$B$, 因此仔细选取$B$依然有助于计算.

3.3 不等式约束

不等式约束存在时, 消元并不总能带来积极的效果. 例如, 若例3中的问题再加上约束$x\ge 0$, 则在消去$x_3,x_6$后, 得到的等价问题则还有三个约束:$$\begin{array}{rl}(x_1,x_2,x_4,x_5)& \ge 0,\\ 8x_1-6x_2+9x_4+4x_5& \le 6,\\ (3/4)x_1+(1/2)x_2-(1/4)x_4+(3/2)x_5& \le -1.\end{array}$$尽管这样去掉了等式约束, 但却让不等式约束变得更加复杂了. 对大多数算法, 这种变换并不会带来多少益处.

但若给例3加上的是$3x_1+2x_3\ge 1$, 则消元就使得等价问题的约束只剩下$$-13x_1+12x_2-18x_4-8x_5\ge -11.$$这时不等式约束就并没有变得更加复杂. 这样消元就是值得的.

4. 价值函数和滤子

假设求解非线性规划问题的算法产生了一步减少了目标函数值但增加了约束的违反度. 我们应当接收这一步吗?

想要回答这个问题并不容易. 我们必须寻求一种权衡二者的方法. 价值函数merit functions法和滤子filters法为达成平衡的两种方法. 在约束优化算法中, 一迭代步$p$被接收仅当它带来了价值函数$\varphi$的充分下降或是它被滤子所接受. 下面我们解释这些概念.

4.1 价值函数

在无约束优化里, 目标函数$f$本身就是价值函数的天然选择. 本书中所介绍的所有无约束优化算法均需要$f$在每步(或隔几步)有所下降.

• 在无约束优化的可行方法(即初始点和后续迭代点均满足问题约束)里, 目标函数仍然是价值函数的一种适宜选择;
• 在不可行方法(及允许迭代点违反约束)中, 我们需要评估迭代步和迭代点的质量, 此时的价值函数就要结合目标函数和约束的违反度.

1. 对于非线性规划问题的一类广受欢迎的价值函数是$l_1$罚函数:$${\varphi }_1(x;\mu )=f(x)+\mu \sum _{i\in \mathcal{E}}|c_i(x)|+\mu \sum _{i\in \mathcal{I}}[c_i(x){]}^{-},$$其中 [ z ] − = max ⁡ { 0 , − z } [z]^{-}=\max\{0,-z\} [z]−=max{0,−z}. 正标量$\mu$为惩罚因子(penalty parameter), 它表现了我们给约束违反度所赋的权重. 注意$l_1$价值函数${\varphi }_1$因绝对值和$[\cdot {]}^{-1}$的存在并不可微, 但它是个精确(exact)的价值函数.
   定义1 (Exact Merit Function) 我们称价值函数$\varphi (x;\mu )$是精确的, 若存在正标量${\mu }^{\ast }$使得对$\forall \mu >{\mu }^{\ast }$, 非线性规划问题的任何局部解都是$\varphi (x;\mu )$的局部极小点.
   我们将在后面证明, 在一定假设条件下, $l_1$价值函数${\varphi }_1(x;\mu )$是精确的, 且阈值${\mu }^{\ast }$由下式给定: μ ∗ = max ⁡ { ∣ λ i ∗ ∣ , i ∈ E ∪ I } , \mu^*=\max\{|\lambda_i^*|,i\in\mathcal{E}\cup\mathcal{I}\}, μ∗=max{∣λi∗​∣,i∈E∪I},其中${\lambda }_i^{\ast }$表示对应于最优解$x^{\ast }$的Lagrange乘子. 但由于最优Lagrange乘子并不预先已知, 基于$l_1$价值函数的算法都会包含调整惩罚因子的策略. 这些策略依赖于所选取的优化算法, 我们将在后面的章节里详细讨论.

2. 另一种常用的价值函数是精确$l_2$函数, 它对不等式约束问题有形式$${\varphi }_2(x;\mu )=f(x)+\mu \Vert c(x){\Vert }_2.$$由于$2-$范数没有平方, 因此它并不可微.

3. 有一些函数既光滑又精确. 为保证函数同时拥有两种性质, 我们必须在价值函数中增加额外的项. 对于等式约束问题, Fletcher增广Lagrange函数Fletcher’s augmented Lagrangian为$${\varphi }_F(x;\mu )=f(x)-\lambda (x{)}^{T}c(x)+\frac{1}{2}\mu \sum _{i\in \mathcal{E}}{c}_i(x{)}^2,$$其中$\mu >0$为惩罚因子, $$\lambda (x)=[A(x){]}^{†}\nabla f(x)=[A(x)A(x{)}^T{]}^{-1}A(x)\nabla f(x).$$这里$A(x)$表示$c(x)$的Jacobi矩阵. 尽管这一函数有一些有趣的理论性质, 但其在实际使用上有所限制. 例如计算$\lambda (x)$和$\lambda (x)$的导数.

4. 一种稍微不同的价值函数为标准增广Lagrange函数. 对不等式约束问题有形式$${\mathcal{L}}_A(x,\lambda ;\mu )=f(x)-{\lambda }^{T}c(x)+\frac{1}{2}\mu \Vert c(x){\Vert }_2^2.$$对此, 我们将通过比较${\mathcal{L}}_A(x^{+},{\lambda }^{+};\mu ),{\mathcal{L}}_A(x,\lambda ;\mu )$的值评估试探点$(x^{+},{\lambda }^{+})$是否可以接受. 严格来讲, ${\mathcal{L}}_A$不是个精确价值函数, 这是因为非线性规划问题的解往往并不是${\mathcal{L}}_A$的极小点, 而只是个稳定点. 尽管一些逐步二次规划算法通过自适应地调整$\mu ,\lambda$将${\mathcal{L}}_A$成功地应用成了价值函数, 但今后我们不再考虑它作为价值函数的使用. 我们仅考虑非光滑的精确罚函数${\varphi }_1,{\varphi }_2$.

由线搜索得到的试探点$x^{+}=x+\alpha p$若能够产生价值函数$\varphi (x;\mu )$上的充分下降, 它就会被算法所接受. 一种定义充分下降的方式类似于无约束优化, 即真实下降量相对于预测下降量不算太小. 尽管$l_1,l_2$函数并不可微, 但它们都有方向导数. 我们将$\varphi (x;\mu )$沿方向$p$的方向导数写作$$D(\varphi (x;\mu );p).$$

1. 线搜索算法中, 充分下降条件要求步长$\alpha >0$充分小以满足不等式$$\varphi (x+\alpha p;\mu )\le \varphi (x;\mu )+\eta \alpha D(\varphi (x;\mu );p),\eta \in (0,1).$$这里$\varphi (x;\mu )-\varphi (x+\alpha p;\mu )$是真实下降, $-\alpha D(\varphi (x;\mu );p)$就是预测下降.
2. 信赖域算法则一般使用一个二次模型$q(p)$估计价值函数$\varphi$的值. 此时充分下降条件可以模型的方式陈述:$$\varphi (x+p;\mu )\le \varphi (x;\mu )-\eta (q(0)-q(p)),\eta \in (0,1).$$这里预测下降变成了$q(0)-q(p)$.

4.2 滤子

滤子法是一种基于多目标优化的迭代步接收机制. 我们知道非线性规划有两个目标: 极小化目标函数、满足约束. 若定义约束违反度为$$h(x)\sum _{i\in \mathcal{E}}|c_i(x)|+\sum _{i\in \mathcal{I}}[c_i(x){]}^{-1},$$我们就可将这两个目标写作$$\underset{x}{\min }f(x),\underset{x}{\min }h(x).$$价值函数可以看做是将两个问题以一定的权重合成为一个极小化问题, 而滤子法则保持这两个目标是分离的. 若$(f(x^{+}),h(x^{+}))$不被先前算法所生成的$(f_l,h_l)=(f(x_l),h(x_l))$占优, 则滤子法就会接收试探点$x^{+}$. “占优”、"滤子"以及"接收"的概念定义如下.

定义2

1. 称二元组$(f_k,h_k)$占优二元组$(f_l,h_l)$若$f_k\le f_l,h_k\le h_l$同时成立.
2. 滤子是一些二元组$(f_l,h_l)$构成的集合, 其中每两个元组互不占优.
3. 称迭代点$x_k$可被当前滤子所接收, 若$(f_k,h_k)$不占优当前滤子中任一二元组.

当迭代点$x_k$被滤子所接收, 我们(通常)将$(f_k,h_k)$加到滤子中, 并移除被其占优的元组. 下图就表示了一个滤子, 滤子中的每一个二元组$(f_l,h_l)$用黑点表示.

如图所示, 滤子中的每个元组都对应了一个无限长方形区域, 而这些区域的并就构成了当前滤子的不可接收域. 具体地说, 试探点$x^{+}$可被滤子接收当且仅当$(f^{+},h^{+})$位于上图实线的左下方.

为比较滤子法和价值函数法, 我们在下图中用虚线画出使得$f+\mu h=f_k+\mu h_k$的所有二元组.

虚线左端的区域就是可降低价值函数$\varphi (x;\mu )=f(x)+\mu h(x)$的点的全体. 显然这与滤子法的可接收域完全不同. 可以看出, 一般一个试探点被滤子法所接收的可能性大于被价值函数法所接收的可能性.

在滤子法中也有线搜索和信赖域. 若线搜索得到的试探点$x^{+}=x_k+{\alpha }_k{p}_k$对应的二元组$(f^{+},h^{+})$为滤子所接收, 则置$x_{k+1}=x^{+}$; 否则启动回溯线搜索. 而在信赖域算法, 若试探点不被滤子接收, 则算法缩小信赖域, 重新计算试探点.

为使滤子法获得全局收敛和较好的实用表现, 我们需对算法做许多强化.

1. 首先, 我们需要保证算法不会接收特别接近当前$(f_k,h_k)$或滤子中二元组的$(f,h)$. 我们可以通过修改接收的判定准则, 加上一个充分下降条件达到目的: 试探点$x^{+}$对滤子可接收, 若对滤子中的所有$(f_j,h_j)$, 都有$$f(x^{+})\le f_j-\beta h_j或h(x^{+})\le h_j-\beta h_j,$$其中$\beta \in (0,1)$. 尽管这种条件很实用(例如取$\beta =10^{-5}$), 但为分析方便将第一个不等式替换为$$f(x^{+})\le f_j-\beta h^{+}.$$
2. 第二个加强方面是要解决滤子法一些内在的问题. 在一定条件下, 由线搜索产生的搜索方向可能需要滤子能接收任意小的步长${\alpha }_k$. 这种现象可能会导致算法停滞甚至失效. 为此, 若回溯线搜索产生的步长小于给定阈值${\alpha }_{\min }$, 则算法就切换到可行性恢复阶段(feasibility restoration phase), 我们等会儿再介绍. 类似地, 在信赖域算法中, 若一串试探点被滤子连续地拒绝, 则信赖域半径可能会过小使得子问题变得不可行. 此时也要调用可行性恢复阶段. 当然我们也可以使用其他的手段, 但我们下面会说到, 可行性恢复阶段对算法还有他用.

可行性恢复阶段专门用于减小约束违反度. 即要求得以下问题的近似解:$$\underset{x}{\min }h(x).$$尽管之前的定义$h(x)$不光滑, 但我们将在十七章说明如何通过求解光滑约束优化子问题解决上述. 此阶段终止于具有充分小$h$值且与滤子兼容的迭代点.

下面我们给出信赖域算法下滤子法的框架. 在之后的章节中我们会讨论约束优化的信赖域算法.

算法1 (General Filter Method)
Choose a starting point $x_0$ and an initial trust-region radius ${\Delta }_0$;
Set $k\leftarrow 0$;
repeat until a convergence test is satisfied
\quad\quad if the step-generation subproblem is infeasible
\quad\quad\quad\quad Compute $x_{k+1}$ using the feasibility restoration phase;
\quad\quad else
\quad\quad\quad\quad Compute a trial iterate $x^{+}=x_k+p_k$;
\quad\quad\quad\quad if $(f^{+},h^{+})$ is acceptable to the filter
\quad\quad\quad\quad\quad\quad Set $x_{k+1}=x^{+}$ and add $(f_{k+1},h_{k+1})$ to the filter;
\quad\quad\quad\quad\quad\quad Choose ${\Delta }_{k+1}$ such that ${\Delta }_{k+1}\ge {\Delta }_k$;
\quad\quad\quad\quad\quad\quad Remove all pairs from the filter that are dominated by $(f_{k+1},h_{k+1})$;
\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad Reject the step, set $x_{k+1}=x_k$;
\quad\quad\quad\quad\quad\quad Choose ${\Delta }_{k+1}<{\Delta }_k$;
\quad\quad\quad\quad end if
\quad\quad end if
\quad\quad $k\leftarrow k+1$;
end repeat

其他的加强措施则依赖于之后章节中的具体算法.

5. Maratos效应

基于价值函数或滤子的一些算法可能无法快速收敛, 是因为它们拒绝了好的试探点. 这种现象被称为Maratos效应(Maratos effect), 这由Maratos首次发现. Maratos效应可用下面的例子阐明, 其中本可以带来二次收敛的$p_k$因带来了目标函数和约束违反度的上升而被拒绝.

例4 考虑问题$$\min f(x_1,x_2)=2(x_1^2+x_2^2-1)-x_1,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{x}_1^2+x_2^2-1=0.$$其图像可见下图.

可以验证最优解为$x^{\ast }=(1,0{)}^T$. 对应的Lagrange乘子为${\lambda }^{\ast }=\frac{3}{2}$, ${\nabla }_{xx}^2\mathcal{L}(x^{\ast },{\lambda }^{\ast })=I$.

考虑迭代点$x_k=(\cos \theta ,\sin \theta {)}^T$, 这对$\theta$的任意值均是可行的. 假设算法计算出的迭代步为$$p_k=\left[\begin{array}{c}{\sin }^2\theta \\ -\sin \theta \cos \theta \end{array}\right],$$则有试探点$$x_k+p_k=\left[\begin{array}{c}\cos \theta +{\sin }^2\theta \\ \sin \theta (1-\cos \theta )\end{array}\right].$$利用三角恒等式, 我们有$$\Vert x_k+p_k-x^{\ast }{\Vert }_2=2{\sin }^2(\theta /2),\Vert x_k-x^{\ast }{\Vert }_2=2|\sin (\theta /2)|,$$因此有$$\frac{\Vert x_k+p_k-x^{\ast }{\Vert }_2}{\Vert x_k-x^{\ast }{\Vert }_2^2}=\frac{1}{2}.$$因此这一步具有二次收敛速度. 但此时我们有$$\begin{array}{rl}f(x_k+p_k)& ={\sin }^2\theta -\cos \theta >-\cos \theta =f(x_k),\\ c(x_k+p_k)& ={\sin }^2\theta >0=c(x_k),\end{array}$$从图像中亦可看到, 目标函数值和约束违反度在这一步都增加了. 这对任何$\theta ̸=0$都成立, 即使初始点无限靠近解.

对于以上例子, 任何基于具有形式$$\varphi (x;\mu )=f(x)+\mu h(c(x))$$(其中$h(\cdot )$为满足$h(0)=0$的非负函数)的价值函数的算法都会拒绝好的试探点. 除此之外, 滤子法也是同样. 因此, 这些方法都会受Maratos效应的影响.

如果不采取补救措施, Maratos效应就可能会影响优化算法的收敛速度. 避免Maratos效应的策略包括以下:

1. 我们可以使用不会被Maratos效应影响的价值函数. 例如Fletcher增广Lagrange函数.
2. 使用二阶校正步(second-order correction step): 在$c(x_k+p_k)$上计算得到${\hat{p}}_k$, 用于减小约束违反度, 再在$p_k$上加上${\hat{p}}_k$.
3. 允许价值函数$\varphi$在一定数目的迭代上增加; 即使用非单调策略(nonmonotone strategy).

我们在下一节讨论后面两种方法.

6. 二阶校正和非单调技术

6.1 二阶校正

通过增加校正项减小约束违反度, 许多算法就能够克服Maratos效应. 我们以等式约束为例介绍这一技术.

给定$p_k$, 二阶校正步(the second-order correction step)${\hat{p}}_k$定义为$${\hat{p}}_k=-A_k^T(A_k{A}_k^T{)}^{-1}c(x_k+p_k),$$其中$A_k=A(x_k)$为$c$在$x_k$处的Jacobi矩阵. 注意${\hat{p}}_k$满足$c$在$x_k+p_k$的一个线性化表示, 即$$A_k{\hat{p}}_k+c(x_k+p_k)=0.$$这一定程度上体现了${\hat{p}}_k$降低约束违反度的作用. 事实上, ${\hat{p}}_k$是这一方程的最小范数解. 后面在第十八章中我们会给出二阶校正的一种不同解释.

在$p_k$满足$A_k{p}_k+c(x_k)=0$时, ${\hat{p}}_k$的作用是将$\Vert c(x)\Vert$减小至$\Vert x_k-x^{\ast }{\Vert }^3$的量阶. 这一估计表明从$x_k$到$x_k+p_k+{\hat{p}}_k$将至少会在解附近减小价值函数取值. 校正的额外成本包括计算$c$在$x_k+p_k$处的取值以及计算${\hat{p}}_k$.

下面给出在线搜索中使用价值函数和二阶校正步的算法. 假设搜索方向$p_k$和惩罚因子${\mu }_k$满足$p_k$为价值函数的下降方向, 即$D(\varphi (x_k;\mu );p_k)<0$. 在第十八章和第十九章, 我们将讨论如何实现这些目标.

这一算法的关键特征是, 若满步${\alpha }_k=1$没有提供价值函数的充分下降, 则算法在沿着原始方向$p_k$回溯之前先尝试二阶校正.

算法2 (Generic Algorithm with Second-Order Correction)
Choose parameters $\eta \in (0,0.5)$ and ${\tau }_1,{\tau }_2:0<{\tau }_1<{\tau }_2<1$;
Choose initial point $x_0$; set $k\leftarrow 0$;
repeat until a convergence test is satisfied:
\quad\quad Compute a search direction $p_k$;
\quad\quad Set ${\alpha }_k\leftarrow 1$, newpoint$\leftarrow$false
\quad\quad while newpoint = false
\quad\quad\quad\quad if $\varphi (x_k+{\alpha }_k{p}_k;\mu )\le \varphi (x_k;\mu )+\eta {\alpha }_{k}D(\varphi (x_k;\mu );p_k)$
\quad\quad\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k+{\alpha }_k{p}_k$;
\quad\quad\quad\quad\quad\quad Set newpoint$\leftarrow$true;
\quad\quad\quad\quad else if ${\alpha }_k=1$
\quad\quad\quad\quad\quad\quad Compute the second-order correction step ${\hat{p}}_k$;
\quad\quad\quad\quad\quad\quad if $\varphi (x_k+p_k+{\hat{p}}_k;\mu )\le \varphi (x_k;\mu )+\eta D(\varphi (x_k;\mu );p_k)$
\quad\quad\quad\quad\quad\quad\quad\quad Set $x_{k+1}\leftarrow x_k+p_k+{\hat{p}}_k$;
\quad\quad\quad\quad\quad\quad\quad\quad Set newpoint$\leftarrow$true;
\quad\quad\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad\quad\quad Choose new ${\alpha }_k$ in $[{\tau }_1{\alpha }_k,{\tau }_2{\alpha }_k]$;
\quad\quad\quad\quad\quad\quad end
\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad Choose new ${\alpha }_k$ in $[{\tau }_1{\alpha }_k,{\tau }_2{\alpha }_k]$;
\quad\quad\quad\quad end
\quad\quad end while
end repeat

在上述算法中, 若${\hat{p}}_k$没有带来价值函数的一定量下降, 就会被舍弃. 我们也并不沿着$p_k+{\hat{p}}_k$回溯, 这是因为这个方向并不一定是价值函数的下降方向. 此算法的一种变体则仅在充分下降条件因约束违反度增加被打破时应用二阶校正步.

二阶校正策略在实际中很高效. 它带来了强健性和效率. 因此额外的计算成本是值得的.

6.2 非单调策略

由Maratos效应带来的不便可通过偶尔接受增加价值函数的试探点避免. 这样的迭代步称作松弛步(relaxed steps). 不过接受是有一定限度的. 如果接连几次(例如$\hat{t}$次)松弛步都没有带来价值函数的充分下降, 我们就回到松弛步之前的迭代点, 实施常规迭代, 利用线搜索或其他技术强制价值函数下降.

二阶校正仅改进约束的满意度, 而非单调策略则意在同时改善可行性与最优性. 我们希望价值函数的增加只是暂时的, 后续的迭代会得到极大的改良.

下面介绍一种非单调方法的特殊形式, 称作是watchdog策略. 置$\hat{t}=1$, 则我们仅允许价值函数在下降之前仅有一次增加. 如6.1, 我们给出线搜索中利用非光滑价值函数$\varphi$的算法. 假设惩罚因子$\mu$在完成循环之前都不变. 为记号简单, 我们忽略$\varphi$对$\mu$的依赖, 将价值函数和方向导数分别写作$\varphi (x),D(\varphi (x);p_k)$.

算法3 (Watchdog)
Choose a constant $\eta \in (0,0.5)$ and an initial point $x_0$;
Set k ← 0 , S ← { 0 } k\leftarrow0,\mathcal{S}\leftarrow\{0\} k←0,S←{0};
repeat until a termination test is satisfied
\quad\quad Compute a step $p_k$;
\quad\quad Set $x_{k+1}\leftarrow x_k+p_k$;
\quad\quad if $\varphi (x_{k+1})\le \varphi (x_k)+\eta D(\varphi (x_k);p_k)$
\quad\quad\quad\quad k ← k + 1 , S ← S ∪ { k } k\leftarrow k+1,\mathcal{S}\leftarrow\mathcal{S}\cup\{k\} k←k+1,S←S∪{k};
\quad\quad else
\quad\quad\quad\quad Compute a search direction $p_{k+1}$ from $x_{k+1}$;
\quad\quad\quad\quad Find ${\alpha }_{k+1}$ such that
\quad\quad\quad\quad\quad\quad $\varphi (x_{k+2})\le \varphi (x_{k+1})+\eta {\alpha }_{k+1}D(\varphi (x_{k+1});p_{k+1})$;
\quad\quad\quad\quad Set $x_{k+2}\leftarrow x_{k+1}+{\alpha }_{k+1}{p}_{k+1}$;
\quad\quad\quad\quad if $\varphi (x_{k+1})\le \varphi (x_k)$ or $\varphi (x_{k+1})\le \varphi (x_k)+\eta D(\varphi (x_k);p_k)$
\quad\quad\quad\quad k ← k + 2 , S ← S ∪ { k } k\leftarrow k+2,\mathcal{S}\leftarrow\mathcal{S}\cup\{k\} k←k+2,S←S∪{k};
\quad\quad\quad\quad else if $\varphi (x_{k+2})>\varphi (x_k)$
\quad\quad\quad\quad\quad\quad (return to $x_k$ and search along $p_k$)
\quad\quad\quad\quad\quad\quad Find ${\alpha }_k$ such that $\varphi (x_{k+3})\le \varphi (x_k)+\eta {\alpha }_{k}D(\varphi (x_k);p_k)$;
\quad\quad\quad\quad\quad\quad Compute $x_{k+3}=x_k+{\alpha }_k{p}_k$;
\quad\quad\quad\quad\quad\quad k ← k + 3 , S ← S ∪ { k } k\leftarrow k+3,\mathcal{S}\leftarrow\mathcal{S}\cup\{k\} k←k+3,S←S∪{k};
\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad Compute a direction $p_{k+2}$ from $x_{k+2}$;
\quad\quad\quad\quad\quad\quad Find ${\alpha }_{k+2}$ such that
\quad\quad\quad\quad\quad\quad\quad\quad $\varphi (x_{k+3})\le \varphi (x_{k+2})+\eta {\alpha }_{k+2}D(\varphi (x_{k+2});p_{k+2});$
\quad\quad\quad\quad\quad\quad Set $x_{k+3}\leftarrow x_{k+2}+{\alpha }_{k+2}{p}_{k+2}$;
\quad\quad\quad\quad\quad\quad k ← k + 3 , S ← S ∪ { k } k\leftarrow k+3,\mathcal{S}\leftarrow\mathcal{S}\cup\{k\} k←k+3,S←S∪{k};
\quad\quad\quad\quad end
\quad\quad end
end (repeat)

其中集合$\mathcal{S}$并没有实质作用, 仅用于储存充分下降价值函数的迭代指标. 注意依以上, 至少有三分之一的指标在$\mathcal{S}$里. 利用这一事实, 我们可以证明许多使用watchdog技术的约束优化算法是全局收敛的. 我们也可以证明, 对于充分大的$k$, 步长终将满足${\alpha }_k=1$, 收敛速度超线性.

实际中, 允许$\hat{t}>1$可能会更好. $\hat{t}$的经典取法为5或8. 基于此, 要较好地使用watchdog技术是需要一定的复杂度的, 但这是值得的. watchdog技术相对于二阶校正策略的一个潜在的优势是, 它可能需要计算更少的约束函数值. 最好情形下, 所有迭代步均是满步, 且基本不需要再返回到之前的迭代点.