第十七章: 惩罚函数法与增广Lagrange函数法

文章目录

第十七章: 惩罚函数法与增广Lagrange函数法
- 1. 二次罚函数法
- - 1.1 动机
  - 1.2 算法框架
  - 1.3 二次罚函数法的收敛性
  - 1.4 病态与重构
- 2. 非光滑精确罚函数法
- - 2.1 实用l1l_1l1罚函数法
  - 2.2 一般的非光滑罚函数法
- 3. 增广Lagrange函数法: 等式约束情形
- - 3.1 动机和算法框架
  - 3.2 增广Lagrange函数的性质
- 4. 实用增广Lagrange函数法
- - 4.1 界约束重构
  - 4.2 约束线性化重构
  - 4.3 无约束化重构
- 5. 不同算法的比较与软件介绍

一些重要的求解约束优化的方法 将原问题替换为一系列子问题, 在子问题中原本的约束被替换成加在目标函数上的项. 本章我们介绍属于此类的三种方法.

二次罚函数法the quadratic penalty method是在目标函数上增加每个约束违反度平方的某个倍数. 这种方法比较简单、直观. 尽管它有许多的缺陷, 但实际中还是被经常使用.
非光滑精确罚函数法the nonsmooth exact penalty methods是用一个(而不是一系列)无约束问题替代原本的约束问题. 使用这种罚函数, 我们可以仅利用一次无约束极小化找到解, 但随之而来, 非光滑性可能会制造一些麻烦. 此类的常见罚函数比如l1l_1l1罚函数.
乘子法the methods of multipliers或增广Lagrange函数法the augmented Lagrangian method是另一种精确罚函数法, 其中显式估计了Lagrange乘子, 避免了二次罚函数带来的病态问题.
对数障碍函数法log-barrier method使用对数项来防止迭代点过于接近可行域边界. 此法是非线性规划内点法的部分基础. 我们放在第十九章讲述.

1. 二次罚函数法

1.1 动机

考虑使用单一个函数代替约束优化问题, 其中那个函数由以下组成:

约束优化问题的目标函数, 加上
对每个约束的一个附加项, 其在当前点违反约束时为正, 否则取0.

许多方法都会通过定义一系列这样的罚函数求解问题, 其中罚项乘上了一个正系数(惩罚因子). 此系数越大, 我们对违反约束的现象越不能容忍, 从而越强制罚函数的极小点靠近约束问题的可行域.

此类最简单的罚函数为二次罚函数(又称Courant罚函数), 其中的罚项记为约束违反度的平方. 我们先在等式约束问题下讨论这一方法.min⁡xf(x),subjecttoci(x)=0,i∈E.\min\limits_xf(x),\quad\mathrm{subject\,to\,}c_i(x)=0,\quad i\in\mathcal{E}.xminf(x),subjecttoci(x)=0,i∈E.对此, 二次罚函数为Q(x;μ)=deff(x)+μ2∑i∈Eci2(x),Q(x;\mu)\xlongequal{def}f(x)+\frac{\mu}{2}\sum_{i\in\mathcal{E}}c_i^2(x),Q(x;μ)deff(x)+2μi∈E∑ci2(x),其中μ>0\mu>0μ>0为惩罚因子(penalty parameter). 当μ\muμ趋向于∞\infty∞, 我们就对约束违反惩罚得愈加猛烈. 直观上, 我们很容易想到要构作一列{μk}:μ↑∞,k→∞\{\mu_k\}:\mu\uparrow\infty,k\to\infty{μk}:μ↑∞,k→∞, 并对每个kkk求得Q(x;μk)Q(x;\mu_k)Q(x;μk)的近似极小点xkx_kxk. 由于上述罚函数是光滑的, 因此我们可以使用无约束优化的技术求xkx_kxk. 在搜索xkx_kxk时, 我们可以使用先前的极小点xk−1,xk−2x_{k-1},x_{k-2}xk−1,xk−2作为算法初始的迭代点. 对于合适的{μk}\{\mu_k\}{μk}和初始点, 我们在每次无约束极小过程中都无需消耗过多的计算量.

例1 考虑以下等式约束问题min⁡x1+x2,subjecttox12+x22−2=0,\min x_1+x_2,\quad\mathrm{subject\,to\,}x_1^2+x_2^2-2=0,minx1+x2,subjecttox12+x22−2=0,其解为(−1,−1)T(-1,-1)^T(−1,−1)T且二次罚函数为Q(x;μ)=x1+x2+μ2(x12+x22−2)2.Q(x;\mu)=x_1+x_2+\frac{\mu}{2}(x_1^2+x_2^2-2)^2.Q(x;μ)=x1+x2+2μ(x12+x22−2)2.下图为μ=1\mu=1μ=1时, QQQ的等高线图. 从图中可见QQQ有一差不多是(−1.1,−1.1)T(-1.1,-1.1)^T(−1.1,−1.1)T的极小点, 也有一接近(0.3,0.3)T(0.3,0.3)^T(0.3,0.3)T的极大点.

而下图则对应μ=10\mu=10μ=10, 因此不再可行域x12+x22=2x_1^2+x_2^2=2x12+x22=2的点会遭受比上图更大的惩罚. 此图中QQQ的极小点愈加靠近解. 同时还有一局部极大点在(0,0)T(0,0)^T(0,0)T附近, 且QQQ在出可行域后就迅速地趋向于∞\infty∞.

但实际情形一般没有例1中那么良态. 对于一给定的μ\muμ, 即使原约束问题有唯一解, 惩罚函数也可能是下无界的. 例如,min⁡−5x12+x22,subjecttox1=1,\min -5x_1^2+x_2^2,\quad\mathrm{subject\,to\,}x_1=1,min−5x12+x22,subjecttox1=1,其解为(1,0)T(1,0)^T(1,0)T. 而对此任何小于10的惩罚因子, 其罚函数都下无界. 不过, 这点对于本章讨论的所有惩罚函数都是存在的.

对于一般的约束优化问题min⁡xf(x),subjecttoci(x)=0,i∈E,ci(x)≥0,i∈I.\min\limits_xf(x),\quad\mathrm{subject\,to\,}c_i(x)=0,i\in\mathcal{E},\quad c_i(x)\ge0,i\in\mathcal{I}.xminf(x),subjecttoci(x)=0,i∈E,ci(x)≥0,i∈I.定义二次罚函数为Q(x;μ)=deff(x)+μ2∑i∈Eci2(x)+μ2∑i∈I([ci(x)]−)2,Q(x;\mu)\xlongequal{def}f(x)+\frac{\mu}{2}\sum_{i\in\mathcal{E}}c_i^2(x)+\frac{\mu}{2}\sum_{i\in\mathcal{I}}\left([c_i(x)]^-\right)^2,Q(x;μ)deff(x)+2μi∈E∑ci2(x)+2μi∈I∑([ci(x)]−)2,这里[y]−=max⁡(−y,0)[y]^-=\max(-y,0)[y]−=max(−y,0). 此时QQQ没有光只有等式约束的惩罚函数光滑, 也不如目标函数和约束函数光滑. 例如若有一个不等式约束为x1≥0x_1\ge0x1≥0, 则min⁡(0,x1)2\min(0,x_1)^2min(0,x1)2就没有连续的二阶导数, QQQ就不再二阶连续可微.

1.2 算法框架

基于二次罚函数的一般算法框架如下.

框架1 (Quadratic Penalty Method)
Given μ0>0\mu_0>0μ0>0, a nonnegative sequence {τk}\{\tau_k\}{τk} with τk→0\tau_k\to0τk→0, and a starting point x0sx_0^sx0s;
for k=0,1,2,…k=0,1,2,\ldotsk=0,1,2,…
\quad\quadFind an approximate minimizer xkx_kxk of Q(⋅;μk)Q(\cdot;\mu_k)Q(⋅;μk), starting at xksx_k^sxks, and terminating when ∥∇xQ(x;μk)∥≤τk\Vert\nabla_xQ(x;\mu_k)\Vert\le\tau_k∥∇xQ(x;μk)∥≤τk;
\quad\quadif final convergence test satisfied
\quad\quad\quad\quadstop with approximate solution xkx_kxk;
\quad\quadend (if)
\quad\quadChoose new penalty parameter μk+1>μk\mu_{k+1}>\mu_kμk+1>μk;
\quad\quadChoose new starting point xk+1sx_{k+1}^sxk+1s;
end (for)

我们可以基于在每次迭代极小化罚函数的困难度自适应地选取惩罚因子序列{μk}\{\mu_k\}{μk}:

当极小化Q(x;μk)Q(x;\mu_k)Q(x;μk)对某一kkk来说比较昂贵, 我们就选取比μk\mu_kμk稍微大一点的μk+1\mu_{k+1}μk+1, 比如μk+1=1.5μk\mu_{k+1}=1.5\mu_kμk+1=1.5μk.
若我们比较容易就能找到Q(x;μk)Q(x;\mu_k)Q(x;μk)的近似极小点, 我们就尝试更大幅度的增长, 比如μk+1=10μk\mu_{k+1}=10\mu_kμk+1=10μk.

框架1的收敛理论给予了非负序列{τk}\{\tau_k\}{τk}选取充分的自由度: 仅需τk→0\tau_k\to0τk→0, 以保证算法随着迭代, 其精度在不断提升.

正如之前讨论过的, 我们无法保证∥∇xQ(x;μk)∥≤τk\Vert\nabla_xQ(x;\mu_k)\Vert\le\tau_k∥∇xQ(x;μk)∥≤τk一定会成立. 因此实际的实施必须包括当约束违反度下降不够快, 或者当迭代点趋向于发散时, 增大惩罚因子(或存储初始点)的防护措施.

当只有等式约束时, Q(x;μk)Q(x;\mu_k)Q(x;μk)光滑, 因此可使用无约束极小化的算法求得近似解xkx_kxk. 不过, 当μk\mu_kμk变大, Q(x;μk)Q(x;\mu_k)Q(x;μk)的极小化会变得愈发困难(求解系统的条件数变大), 除非我们使用特殊的手段计算搜索方向.

一方面, Hessian阵∇xx2Q(x;μk)\nabla_{xx}^2Q(x;\mu_k)∇xx2Q(x;μk)会在极小点附近变得愈发病态. 这一点就足以使得许多无约束算法(如拟牛顿法、共轭梯度法)效果变差. 而诸如牛顿法这般对Hessian条件数不敏感的算法, 则可能会遭遇较大μk\mu_kμk带来的另外两个问题:

在我们求解线性系统计算牛顿步时, 病态的∇xx2Q(x;μk)\nabla_{xx}^2Q(x;\mu_k)∇xx2Q(x;μk)可能会引发数值不稳定. 后面我们会指出, 这种影响并不严重, 我们可以重构牛顿方程.
即使xxx很接近Q(⋅;μk)Q(\cdot;\mu_k)Q(⋅;μk)的极小点, 但对于Q(x;μk)Q(x;\mu_k)Q(x;μk)的二阶Taylor展开仅在xxx的一个小邻域内才是合适的. 这一点可从第二张图看出来, 其中QQQ在极小点附近的等高线呈现出一种"香蕉"状而不是二次函数典型的椭圆型. 由于牛顿法是基于二次模型的, 因此其产生的迭代步可能无法快速收敛到Q(x;μk)Q(x;\mu_k)Q(x;μk)的极小点. 对此, 我们需要谨慎地选择初始点xk+1sx_{k+1}^sxk+1s或直接设置xk+1s=xkx_{k+1}^s=x_kxk+1s=xk, 并选取μk+1\mu_{k+1}μk+1只比μk\mu_kμk大一点.

1.3 二次罚函数法的收敛性

我们在下面两个定理里介绍二次罚函数方法的一些收敛性质. 我们将讨论限制在等式约束问题.

对于第一个结果, 我们假设对每个μk\mu_kμk, 罚函数Q(x;μk)Q(x;\mu_k)Q(x;μk)都有(有限个)极小点.

定理1 设xkx_kxk为框架1中Q(x;μk)Q(x;\mu_k)Q(x;μk)的精确全局极小点, 且μk↑∞\mu_k\uparrow\inftyμk↑∞. 则序列{xk}\{x_k\}{xk}的每个聚点都是原问题的全局解.

证明: 令xˉ\bar{x}xˉ为原问题的全局解. 即有f(xˉ)≤f(x),∀x:ci(x)=0,i∈E.f(\bar{x})\le f(x),\quad\forall x:c_i(x)=0,i\in\mathcal{E}.f(xˉ)≤f(x),∀x:ci(x)=0,i∈E.由于对每个kkk, xkx_kxk极小化Q(⋅;μk)Q(\cdot;\mu_k)Q(⋅;μk), 于是Q(xk;μk)≤Q(xˉ;μk)Q(x_k;\mu_k)\le Q(\bar{x};\mu_k)Q(xk;μk)≤Q(xˉ;μk), 这就得到不等式f(xk)+μk2∑i∈Eci2(xk)≤f(xˉ)+μk2∑i∈Eci2(xˉ)=f(xˉ).f(x_k)+\frac{\mu_k}{2}\sum_{i\in\mathcal{E}}c_i^2(x_k)\le f(\bar{x})+\frac{\mu_k}{2}\sum_{i\in\mathcal{E}}c_i^2(\bar{x})=f(\bar{x}).f(xk)+2μki∈E∑ci2(xk)≤f(xˉ)+2μki∈E∑ci2(xˉ)=f(xˉ).重新整理可得∑i∈Eci2(xk)≤2μk[f(xˉ)−f(xk)].\sum_{i\in\mathcal{E}}c_i^2(x_k)\le\frac{2}{\mu_k}[f(\bar{x})-f(x_k)].i∈E∑ci2(xk)≤μk2[f(xˉ)−f(xk)].假设x∗x^*x∗为{xk}\{x_k\}{xk}的一个聚点, 因此存在{1,2,…,n,…}\{1,2,\ldots,n,\ldots\}{1,2,…,n,…}的无穷子列K\mathcal{K}K使得lim⁡k→Kxk=x∗.\lim_{k\to\mathcal{K}}x_k=x^*.k→Klimxk=x∗.对上面不等式两边取极限k→∞,k∈Kk\to\infty,k\in\mathcal{K}k→∞,k∈K, 我们有∑i∈Eci2(x∗)=lim⁡k∈K∑i∈Eci2(xk)≤lim⁡k∈K2μk[f(xˉ)−f(xk)]=0.\sum_{i\in\mathcal{E}}c_i^2(x^*)=\lim_{k\in\mathcal{K}}\sum_{i\in\mathcal{E}}c_i^2(x_k)\le\lim_{k\in\mathcal{K}}\frac{2}{\mu_k}[f(\bar{x})-f(x_k)]=0.i∈E∑ci2(x∗)=k∈Klimi∈E∑ci2(xk)≤k∈Klimμk2[f(xˉ)−f(xk)]=0.因此ci(x∗)=0,i∈Ec_i(x^*)=0,i\in\mathcal{E}ci(x∗)=0,i∈E, 从而x∗x^*x∗可行. 进一步, f(x∗)≤f(x∗)+lim⁡k∈Kμk2∑i∈Eci2(xk)≤f(xˉ).f(x^*)\le f(x^*)+\lim_{k\in\mathcal{K}}\frac{\mu_k}{2}\sum_{i\in\mathcal{E}}c_i^2(x_k)\le f(\bar{x}).f(x∗)≤f(x∗)+k∈Klim2μki∈E∑ci2(xk)≤f(xˉ).由于x∗x^*x∗可行, 且其函数值不大于全局解xˉ\bar{x}xˉ处的函数值, 因此x∗x^*x∗也是全局解. 证毕.

此结论对带不等式的问题也成立. 但由于需要我们求每个子问题的全局极小点, 因此此性质一般并不能成立. 下一结论则允许我们不精确(但精确度要提升)极小化Q(⋅;μk)Q(\cdot;\mu_k)Q(⋅;μk). 相较于定理1, 它表示{xk}\{x_k\}{xk}可能会收敛到不可行点或是KKT点. 它也说明了在一定情形下, −μkci(xk)-\mu_kc_i(x_k)−μkci(xk)可以当做Lagrange乘子λi∗\lambda_i^*λi∗的估计. 这一点对于我们在第3节讨论增广Lagrange函数法时比较重要.

为建立定理, 我们需乐观地假设对所有kkk, ∥∇xQ(x;μk)∥≤τk\Vert\nabla_xQ(x;\mu_k)\Vert\le\tau_k∥∇xQ(x;μk)∥≤τk都会成立.

定理2 设框架1中的容忍限序列和惩罚因子序列满足τk→0,μk↑∞\tau_k\to0,\mu_k\uparrow\inftyτk→0,μk↑∞, 则若序列{xk}\{x_k\}{xk}的聚点x∗x^*x∗不可行, 它就是函数∥c(x)∥2\Vert c(x)\Vert^2∥c(x)∥2的稳定点. 若聚点x∗x^*x∗可行且约束梯度∇ci(x∗)\nabla c_i(x^*)∇ci(x∗)线性无关, 则x∗x^*x∗为原问题的KKT点. 于此, 对任一收敛到x∗x^*x∗的子列, 即∀K:lim⁡k∈Kxk=x∗\forall\mathcal{K}:\lim_{k\in\mathcal{K}}x_k=x^*∀K:limk∈Kxk=x∗, 我们有lim⁡k∈K−μkci(xk)=λi∗,∀i∈E,\lim_{k\in\mathcal{K}}-\mu_kc_i(x_k)=\lambda_i^*,\quad\forall i\in\mathcal{E},k∈Klim−μkci(xk)=λi∗,∀i∈E,这里λ∗\lambda^*λ∗为满足等式约束原问题KKT条件的乘子.

证明: 对Q(x;μk)Q(x;\mu_k)Q(x;μk)求导, 得到∇xQ(xk;μk)=∇f(xk)+∑i∈Eμkci(xk)∇ci(xk),\nabla_xQ(x_k;\mu_k)=\nabla f(x_k)+\sum_{i\in\mathcal{E}}\mu_kc_i(x_k)\nabla c_i(x_k),∇xQ(xk;μk)=∇f(xk)+i∈E∑μkci(xk)∇ci(xk),于是从终止条件, 我们有∥∇f(xk)+∑i∈Eμkci(xk)∇ci(xk)∥≤τk.\left\Vert\nabla f(x_k)+\sum_{i\in\mathcal{E}}\mu_kc_i(x_k)\nabla c_i(x_k)\right\Vert\le\tau_k.∥∥∥∥∥∇f(xk)+i∈E∑μkci(xk)∇ci(xk)∥∥∥∥∥≤τk.利用三角不等式, 有∥∑i∈Eci(xk)∇ci(xk)∥≤1μk[τk+∥∇f(xk)∥].\left\Vert\sum_{i\in\mathcal{E}}c_i(x_k)\nabla c_i(x_k)\right\Vert\le\frac{1}{\mu_k}[\tau_k+\Vert\nabla f(x_k)\Vert].∥∥∥∥∥i∈E∑ci(xk)∇ci(xk)∥∥∥∥∥≤μk1[τk+∥∇f(xk)∥].令x∗x^*x∗为迭代点序列的聚点. 于是存在无穷指标子列K\mathcal{K}K使得lim⁡k∈Kxk=x∗\lim_{k\in\mathcal{K}}x_k=x^*limk∈Kxk=x∗. 对上式取k→∞,k∈Kk\to\infty,k\in\mathcal{K}k→∞,k∈K的极限, 可得∑i∈Eci(x∗)∇ci(x∗)=0.\sum_{i\in\mathcal{E}}c_i(x^*)\nabla c_i(x^*)=0.i∈E∑ci(x∗)∇ci(x∗)=0.我们可能有ci(x∗)≠0c_i(x^*)\ne0ci(x∗)=0(若约束梯度∇ci(x∗)\nabla c_i(x^*)∇ci(x∗)线性相关), 但这起码说明x∗x^*x∗为∥c(x)∥2\Vert c(x)\Vert^2∥c(x)∥2的稳定点.

若∇ci(x∗)\nabla c_i(x^*)∇ci(x∗)线性无关, 则ci(x∗)=0,∀i∈Ec_i(x^*)=0,\forall i\in\mathcal{E}ci(x∗)=0,∀i∈E, 因此x∗x^*x∗可行. 于是KKT条件之可行性条件成立. 记约束的Jacobi阵为AAA,向量−μkc(xk)-\mu_kc(x_k)−μkc(xk)为λk\lambda_kλk, 则有A(xk)Tλk=∇f(xk)−∇xQ(xk;μk),∥∇xQ(xk;μk)∥≤τk.A(x_k)^T\lambda_k=\nabla f(x_k)-\nabla_xQ(x_k;\mu_k),\quad\Vert\nabla_xQ(x_k;\mu_k)\Vert\le\tau_k.A(xk)Tλk=∇f(xk)−∇xQ(xk;μk),∥∇xQ(xk;μk)∥≤τk.对k∈Kk\in\mathcal{K}k∈K充分大, A(xk)A(x_k)A(xk)行满秩, 于是A(xk)A(xk)TA(x_k)A(x_k)^TA(xk)A(xk)T非奇异. 上式左乘A(xk)A(x_k)A(xk)并整理可得λk=[A(xk)A(xk)T]−1A(xk)[∇f(xk)−∇xQ(xk;μk)].\lambda_k=\left[A(x_k)A(x_k)^T\right]^{-1}A(x_k)[\nabla f(x_k)-\nabla_xQ(x_k;\mu_k)].λk=[A(xk)A(xk)T]−1A(xk)[∇f(xk)−∇xQ(xk;μk)].取极限k→∞,k∈Kk\to\infty,k\in\mathcal{K}k→∞,k∈K, 可得lim⁡k∈Kλk=λ∗=[A(x∗)A(x∗)T]−1A(x∗)∇f(x∗).\lim_{k\in\mathcal{K}}\lambda_k=\lambda^*=\left[A(x^*)A(x^*)^T\right]^{-1}A(x^*)\nabla f(x^*).k∈Klimλk=λ∗=[A(x∗)A(x∗)T]−1A(x∗)∇f(x∗).又由前面取极限可知∇f(x∗)−A(x∗)Tλ∗=0,\nabla f(x^*)-A(x^*)^T\lambda^*=0,∇f(x∗)−A(x∗)Tλ∗=0,于是λ∗\lambda^*λ∗满足KKT条件之稳定性条件. 所以, x∗x^*x∗为原问题的KKT点, 且由唯一的Lagrange乘子λ∗\lambda^*λ∗. 证毕.

需重复强调的是, 若聚点x∗x^*x∗不可行, 则它至少是∥c(x)∥2\Vert c(x)\Vert^2∥c(x)∥2的稳定点. 牛顿类的算法总是会陷在这种类型的点. 这有点像我们在第十一章中讨论非线性方程组的情形. 在原问题不可行时, 我们往往可以观察到二次罚函数算法收敛于∥c(x)∥2\Vert c(x)\Vert^2∥c(x)∥2的稳定点或极小点.

若考虑不等式约束, 则x∗x^*x∗为函数∥[c(x)]−∥2\Vert [c(x)]^-\Vert^2∥[c(x)]−∥2的稳定点, 其中[c(x)]−[c(x)]^-[c(x)]−定义为[c(x)]i−={ci(x),i∈E,[ci(x)]−,i∈I.[c(x)]_i^-=\left\{\begin{array}{ll}c_i(x), & i\in\mathcal{E},\\ [c_i(x)]^-, & i\in\mathcal{I}.\end{array}\right.[c(x)]i−={ci(x),[ci(x)]−,i∈E,i∈I.但所谓的KKT点不一定成立, 因为此处无法估计乘子的正负.

1.4 病态与重构

我们现在来验证Hessian阵∇xx2Q(x;μk)\nabla^2_{xx}Q(x;\mu_k)∇xx2Q(x;μk)的病态内在. 对这一矩阵和在其他罚函数和障碍函数法中出现的Hessian阵性质的理解, 将有助于我们选取和设计合适、高效的算法.

这里的Hessian阵为∇xx2Q(x;μk)=∇2f(x)+∑i∈Eμkci(x)∇2ci(x)+μkA(x)TA(x).\nabla^2_{xx}Q(x;\mu_k)=\nabla^2f(x)+\sum_{i\in\mathcal{E}}\mu_kc_i(x)\nabla^2c_i(x)+\mu_kA(x)^TA(x).∇xx2Q(x;μk)=∇2f(x)+i∈E∑μkci(x)∇2ci(x)+μkA(x)TA(x).当xxx接近Q(⋅;μk)Q(\cdot;\mu_k)Q(⋅;μk)的极小点且定理2的条件满足时, 我们从定理2的结论可知, 上式右端前两项就差不多是Lagrange函数的Hessian. 具体地, ∇xx2Q(x;μk)≈∇xx2L(x,λ∗)+μkA(x)TA(x).\nabla^2_{xx}Q(x;\mu_k)\approx\nabla^2_{xx}\mathcal{L}(x,\lambda^*)+\mu_kA(x)^TA(x).∇xx2Q(x;μk)≈∇xx2L(x,λ∗)+μkA(x)TA(x).依此, ∇xx2Q(x;μk)\nabla^2_{xx}Q(x;\mu_k)∇xx2Q(x;μk)就接近于以下两个部分的和:

(Lagrange项)独立于μk\mu_kμk的矩阵; 和
秩为∣E∣|\mathcal{E}|∣E∣, 且非零特征值与μk\mu_kμk同阶的矩阵.

由于约束的数量∣E∣|\mathcal{E}|∣E∣通常小于nnn, 因此第二项奇异. 因此整个矩阵的某些特征值会趋近常数, 而其他的则与μk\mu_kμk同阶. 因μk→∞\mu_k\to\inftyμk→∞, 通过条件数可知, ∇xx2Q(x;μk)\nabla^2_{xx}Q(x;\mu_k)∇xx2Q(x;μk)会变得愈发病态.

Hessian病态的一个后果是, 以下牛顿步的计算可能会不够精确:∇xx2Q(x;μk)p=−∇xQ(x;μk).\nabla^2_{xx}Q(x;\mu_k)p=-\nabla_xQ(x;\mu_k).∇xx2Q(x;μk)p=−∇xQ(x;μk).一般来说, 不论使用何种直接计算手段, 此病态系统都会导致在ppp上的巨大计算误差. 基于同样的原因, 除非有预处理的策略, 否则迭代计算的方法也不会有较好的表现.

不过, 我们可以重构以上牛顿方程以避免病态. 引入一个新的变量ζ:=μA(x)p\zeta:=\mu A(x)pζ:=μA(x)p, 我们会发现满足牛顿方程的ppp同样满足以下系统:[∇2f(x)+∑i∈Eμkci(x)∇2ci(x)A(x)TA(x)−(1/μk)I][pζ]=[−∇xQ(x;μk)0].\begin{bmatrix}\nabla^2f(x)+\sum\limits_{i\in\mathcal{E}}\mu_kc_i(x)\nabla^2c_i(x) & A(x)^T\\A(x) & -(1/\mu_k)I\end{bmatrix}\begin{bmatrix}p\\\zeta\end{bmatrix}=\begin{bmatrix}-\nabla_xQ(x;\mu_k)\\0\end{bmatrix}.[∇2f(x)+i∈E∑μkci(x)∇2ci(x)A(x)A(x)T−(1/μk)I][pζ]=[−∇xQ(x;μk)0].当xxx离x∗x^*x∗较近时, 此系统的系数矩阵并不会有与μk\mu_kμk同阶的大特征值. 事实上, 当μk→∞\mu_k\to\inftyμk→∞, 可将此系统的系数矩阵近似看做[∇xx2L(x,λ∗)A(x)TA(x)O],\begin{bmatrix}\nabla^2_{xx}\mathcal{L}(x,\lambda^*) & A(x)^T\\A(x) & O\end{bmatrix},[∇xx2L(x,λ∗)A(x)A(x)TO],进一步若假设∇xx2L(x,λ∗)\nabla^2_{xx}\mathcal{L}(x,\lambda^*)∇xx2L(x,λ∗)非奇异, 利用分块矩阵的初等变换, 可得[∇xx2L(x,λ∗)A(x)TA(x)O]∼[∇xx2L(x,λ∗)A(x)TO−A(x)[∇xx2L(x,λ∗)]−1A(x)T],\begin{bmatrix}\nabla^2_{xx}\mathcal{L}(x,\lambda^*) & A(x)^T\\A(x) & O\end{bmatrix}\sim\begin{bmatrix}\nabla^2_{xx}\mathcal{L}(x,\lambda^*) & A(x)^T\\O & -A(x)\left[\nabla^2_{xx}\mathcal{L}(x,\lambda^*)\right]^{-1}A(x)^T\end{bmatrix},[∇xx2L(x,λ∗)A(x)A(x)TO]∼[∇xx2L(x,λ∗)OA(x)T−A(x)[∇xx2L(x,λ∗)]−1A(x)T],此阵的行列式不为0. 因此此系统可视作牛顿方程的良态重构(well conditioned reformulation). 不过, 由于∇2f(x)+∑i∈Eμkci(x)∇2ci(x)\nabla^2f(x)+\sum_{i\in\mathcal{E}}\mu_kc_i(x)\nabla^2c_i(x)∇2f(x)+∑i∈Eμkci(x)∇2ci(x)对∇xx2L(x,λ∗)\nabla_{xx}^2\mathcal{L}(x,\lambda^*)∇xx2L(x,λ∗)(也即μkci(x)\mu_kc_i(x)μkci(x)对−λi∗-\lambda_i^*−λi∗)的近似可能不会很好, 所以这两个系统可能都不会得出一个比较好的搜索方向ppp. 这一事实可能就说明二次模型并不是Q(⋅;μk)Q(\cdot;\mu_k)Q(⋅;μk)的一个合适的近似模型, 而牛顿步可能本身就不适合作为搜索方向. 后续我们会讨论针对以上的补救措施.

如果要通过求解重构的问题计算迭代步, 我们需要求解一个n+∣E∣n+|\mathcal{E}|n+∣E∣维而不是nnn维的系统. 后面在下一章中讨论逐步二次规划(SQP)时我们也得求解类似的系统. 事实上,

当μk\mu_kμk很大, 上述重构可视作SQP迭代步的正则化, 其中的−(1/μk)I-(1/\mu_k)I−(1/μk)I使得迭代矩阵即使A(x)A(x)A(x)行亏秩时依然是非奇异的.
当μk\mu_kμk较小, 则重构系统的第二行说明计算出的迭代步并不能较好地满足约束的线性化. 我们是不希望这种情况发生的, 因为这样一来迭代步就不能朝着可行性改进多少, 从而全局表现不够高效.
若{μk}\{\mu_k\}{μk}趋于无穷过快, 我们就又可能丧失线性化满足时获得超线性收敛速度的机会. 具体讨论可见下一章.

总之, 重构系统可以视作无约束极小在二次惩罚函数Q(⋅;μk)Q(\cdot;\mu_k)Q(⋅;μk)上的应用, 也可以视作SQP方法的一个变体.

2. 非光滑精确罚函数法

有些罚函数是精确的(exact), 即对于特定的惩罚因子, 仅做一次极小化即可得到非线性规划问题的精确解. 这一性质是很吸引人的, 因为它使得罚函数方法不依赖于惩罚因子的更新策略. 第1节介绍的二次罚函数并不是精确的, 这是因为对任何μ\muμ的正值, 其极小点一般都不是原非线性规划的解. 本节我们讨论非光滑精确罚函数法.
对一般非线性规划的一种广受欢迎的非光滑罚函数是l1l_1l1罚函数, 定义为ϕ1(x;μ)=f(x)+μ∑i∈E∣ci(x)∣+μ∑i∈I[ci(x)]−.\phi_1(x;\mu)=f(x)+\mu\sum_{i\in\mathcal{E}}|c_i(x)|+\mu\sum_{i\in\mathcal{I}}[c_i(x)]^-.ϕ1(x;μ)=f(x)+μi∈E∑∣ci(x)∣+μi∈I∑[ci(x)]−.其称谓源于它使用的惩罚项是μ\muμ乘以约束违反度的l1l_1l1范数. 注意ϕ1(x;μ)\phi_1(x;\mu)ϕ1(x;μ)在某些xxx上并不可微.

下面的定理揭示了l1l_1l1罚函数的精确性.

定理3 设x∗x^*x∗为非线性规划问题的严格局部解, 且其满足一阶必要性条件, 有Lagrange乘子λi∗,i∈E∪I\lambda_i^*,i\in\mathcal{E}\cup\mathcal{I}λi∗,i∈E∪I. 于是x∗x^*x∗为ϕ1(x;μ),∀μ>μ∗\phi_1(x;\mu),\forall\mu>\mu^*ϕ1(x;μ),∀μ>μ∗的局部极小点, 其中μ∗=∥λ∗∥∞=max⁡i∈E∪I∣λi∗∣.\mu^*=\Vert\lambda^*\Vert_{\infty}=\max_{i\in\mathcal{E}\cup\mathcal{I}}|\lambda_i^*|.μ∗=∥λ∗∥∞=i∈E∪Imax∣λi∗∣.另外, 若在x∗x^*x∗还满足二阶充分性条件, 且μ>μ∗\mu>\mu^*μ>μ∗, 则x∗x^*x∗为ϕ1(x;μ)\phi_1(x;\mu)ϕ1(x;μ)的严格局部极小点.

证明可见 Han S P和Mangasarian O L¹的定理4.4.

粗略地说, 在非线性规划的解x∗x^*x∗处, 任何向不可行域的移动都会被严重地惩罚, 这样使得函数值会大于ϕ1(x∗;μ)=f(x∗)\phi_1(x^*;\mu)=f(x^*)ϕ1(x∗;μ)=f(x∗)(即方向导数都大于0), 强迫ϕ(⋅;μ)\phi(\cdot;\mu)ϕ(⋅;μ)的极小点落在x∗x^*x∗.

例2 考虑如下单变量问题:min⁡x,subjecttox≥1,\min x,\quad\mathrm{subject\,to\,}x\ge1,minx,subjecttox≥1,其解为x∗=1x^*=1x∗=1. 我们有ϕ1(x;μ)=x+μ[x−1]−={(1−μ)x+μ,x≤1,xx>1.\phi_1(x;\mu)=x+\mu[x-1]^-=\left\{\begin{array}{ll}(1-\mu)x+\mu, & x\le1,\\x & x>1.\end{array}\right.ϕ1(x;μ)=x+μ[x−1]−={(1−μ)x+μ,xx≤1,x>1.于是从下图可见, 当μ>1\mu>1μ>1, 罚函数有极小点x∗=1x^*=1x∗=1, 但在μ<1\mu<1μ<1时却是个单调递增函数.

由于罚函数法通过直接极小化罚函数起作用, 所以我们需要刻画ϕ1\phi_1ϕ1的稳定点. 即使ϕ1\phi_1ϕ1不可微, 但它沿任何方向具有方向导数D(ϕ1(x;μ);p)D(\phi_1(x;\mu);p)D(ϕ1(x;μ);p).

定义1 一点x^∈Rn\hat{x}\in\mathbb{R}^nx^∈Rn是罚函数ϕ1(x;μ)\phi_1(x;\mu)ϕ1(x;μ)的稳定点, 若D(ϕ1(x^;μ);p)≥0,∀p∈Rn.D(\phi_1(\hat{x};\mu);p)\ge0,\quad\forall p\in\mathbb{R}^n.D(ϕ1(x^;μ);p)≥0,∀p∈Rn.类似地, x^\hat{x}x^为不可行度量(the measure of infeasibility)h(x)=∑i∈E∣ci(x)∣+∑i∈I[ci(x)]−h(x)=\sum_{i\in\mathcal{E}}|c_i(x)|+\sum_{i\in\mathcal{I}}[c_i(x)]^-h(x)=i∈E∑∣ci(x)∣+i∈I∑[ci(x)]−的稳定点, 若D(h(x^);p)≥0,∀p∈RnD(h(\hat{x});p)\ge0,\forall p\in\mathbb{R}^nD(h(x^);p)≥0,∀p∈Rn. 若一点不可行但对不可行度量hhh是稳定的, 则我们称其为不可行稳定点(infeasible stationary point).

对于例2中的函数, 我们在x∗=1x^*=1x∗=1有D(ϕ1(x∗;μ);p)={p,p≥0,(1−μ)p,p<0.D(\phi_1(x^*;\mu);p)=\left\{\begin{array}{ll}p, & p\ge0,\\(1-\mu)p, & p<0.\end{array}\right.D(ϕ1(x∗;μ);p)={p,(1−μ)p,p≥0,p<0.于是当μ>1\mu>1μ>1, D(ϕ1(x∗;μ);p)≥0,∀p∈RD(\phi_1(x^*;\mu);p)\ge0,\forall p\in\mathbb{R}D(ϕ1(x∗;μ);p)≥0,∀p∈R.

下面的结论则在一定程度上弥补了定理3的不足, 即证明了反方向: 在一定条件下, ϕ1(x;μ)\phi_1(x;\mu)ϕ1(x;μ)的稳定点对应了原非线性规划的KKT点.

定理4 设x^\hat{x}x^为罚函数ϕ1(x;μ),∀μ>μ^>0\phi_1(x;\mu),\forall\mu>\hat{\mu}>0ϕ1(x;μ),∀μ>μ^>0的稳定点. 若x^\hat{x}x^为原非线性规划的可行点, 则其为KKT点. 若不可行, 则是不可行稳定点.

证明: 若x^\hat{x}x^可行, 于是D(ϕ1(x^;μ);p)=∇f(x^)Tp+μ∑i∈E∣∇ci(x^)Tp∣+μ∑i∈I∩A(x^)[∇ci(x^)Tp]−,???D(\phi_1(\hat{x};\mu);p)=\nabla f(\hat{x})^Tp+\mu\sum_{i\in\mathcal{E}}\left|\nabla c_i(\hat{x})^Tp\right|+\mu\sum_{i\in\mathcal{I}\cap\mathcal{A}(\hat{x})}\left[\nabla c_i(\hat{x})^Tp\right]^-,???D(ϕ1(x^;μ);p)=∇f(x^)Tp+μi∈E∑∣∣∇ci(x^)Tp∣∣+μi∈I∩A(x^)∑[∇ci(x^)Tp]−,???其中A(x^)\mathcal{A}(\hat{x})A(x^)为在x^\hat{x}x^处的积极集. 考虑任何线性化可行方向集F(x^)\mathcal{F}(\hat{x})F(x^)中的ppp, 由F(x^)\mathcal{F}(\hat{x})F(x^)的定义, 我们有∣∇ci(x^)Tp∣+∑i∈I∩A(x^)[∇ci(x^)Tp]−=0,\left|\nabla c_i(\hat{x})^Tp\right|+\sum_{i\in\mathcal{I}\cap\mathcal{A}(\hat{x})}\left[\nabla c_i(\hat{x})^Tp\right]^-=0,∣∣∇ci(x^)Tp∣∣+i∈I∩A(x^)∑[∇ci(x^)Tp]−=0,因此由稳定性假设, 我们有0≤D(ϕ1(x^;μ);p)=∇f(x^)Tp,∀p∈F(x^).0\le D(\phi_1(\hat{x};\mu);p)=\nabla f(\hat{x})^Tp,\quad\forall p\in\mathcal{F}(\hat{x}).0≤D(ϕ1(x^;μ);p)=∇f(x^)Tp,∀p∈F(x^).利用Farkas引理, 可得存在λ^i:λ^i≥0,i∈I∩A(x^)\hat{\lambda}_i:\hat{\lambda}_i\ge0,i\in\mathcal{I}\cap\mathcal{A}(\hat{x})λ^i:λ^i≥0,i∈I∩A(x^), 使得∇f(x^)=∑i∈A(x^)λ^i∇ci(x^).\nabla f(\hat{x})=\sum_{i\in\mathcal{A}(\hat{x})}\hat{\lambda}_i\nabla c_i(\hat{x}).∇f(x^)=i∈A(x^)∑λ^i∇ci(x^).这说明x^\hat{x}x^为KKT点.

若x^\hat{x}x^不可行, 反证, 假设存在ppp使得D(h(x^);p)<0D(h(\hat{x});p)<0D(h(x^);p)<0. 由0≤D(ϕ1(x^;μ);p)=∇f(x^)Tp+μD(h(x^);p),∀μ充分大,0\le D(\phi_1(\hat{x};\mu);p)=\nabla f(\hat{x})^Tp+\mu D(h(\hat{x});p),\forall\mu充分大,0≤D(ϕ1(x^;μ);p)=∇f(x^)Tp+μD(h(x^);p),∀μ充分大,可推得∇f(x^)Tp→∞\nabla f(\hat{x})^Tp\to\infty∇f(x^)Tp→∞, 矛盾! 证毕.

例3 再考虑例1中的问题, 这里使用l1l_1l1罚函数,ϕ1(x;μ)=x1+x2+μ∣x12+x22−2∣.\phi_1(x;\mu)=x_1+x_2+\mu|x_1^2+x_2^2-2|.ϕ1(x;μ)=x1+x2+μ∣x12+x22−2∣.下图描绘了ϕ(x;2)\phi(x;2)ϕ(x;2)的等高线, 其极小点可见即为原问题的解x∗=(−1,−1)Tx^*=(-1,-1)^Tx∗=(−1,−1)T.

事实上由定理3, 我们知道对于所有的μ>∣λ∗∣=0.5\mu>|\lambda^*|=0.5μ>∣λ∗∣=0.5, ϕ1(x;μ)\phi_1(x;\mu)ϕ1(x;μ)的极小点都是x∗x^*x∗. 而图中等高线的尖角则体现了l1l_1l1罚函数的不光滑性.

这些结论为基于l1l_1l1罚函数的算法框架提供了动机.

框架2 (Classical l1l_1l1 Penalty Method)
Given μ0>0\mu_0>0μ0>0, tolerance τ>0\tau>0τ>0, starting point x0sx_0^sx0s;
for k=0,1,2,…k=0,1,2,\ldotsk=0,1,2,…
\quad\quadFind an approximate minimizer xkx_kxk of ϕ1(x;μk)\phi_1(x;\mu_k)ϕ1(x;μk), starting at xksx_k^sxks;
\quad\quadif h(xk)≤τh(x_k)\le\tauh(xk)≤τ
\quad\quad\quad\quadstop with approximate solution xkx_kxk;
\quad\quadend (if)
\quad\quadChoose new penalty parameter μk+1>μk\mu_{k+1}>\mu_kμk+1>μk;
\quad\quadChoose new starting point xk+1sx_{k+1}^sxk+1s;
end (for)

这其中对于ϕ1(x;μk)\phi_1(x;\mu_k)ϕ1(x;μk)的极小因非光滑性而变得难以实现. 但若使用ϕ1(x;μk)\phi_1(x;\mu_k)ϕ1(x;μk)的光滑模型近似, 再对光滑模型做极小, 事情就变得好理解多了. 我们将在下一小节介绍, 这有点类似于逐步二次规划算法.

更新惩罚因子μk\mu_kμk的最简单格式是乘以某个常数(比如5或10). 此格式有时很好用, 但也可能会影响效率; 若初始的惩罚因子μ0\mu_0μ0过小, 则框架2可能需要更多的迭代更新才能取到一个比较合适的因子; 在迭代早期, 迭代点可能会逐渐远离解x∗x^*x∗, 这时就需要尽早终止ϕ1(x;μk)\phi_1(x;\mu_k)ϕ1(x;μk)的极小化过程, xksx_k^sxks也应该重新设置为一个先前的迭代点; 若μk\mu_kμk相当大, 则极小化罚函数的困难度也会增加, 需要的迭代书也会变多. 之后我们会讨论如何选取惩罚因子.

2.1 实用l1l_1l1罚函数法

如前所述, ϕ1(x;μ)\phi_1(x;\mu)ϕ1(x;μ)在任何x:ci(x)=0,∃i∈E∪Ix:c_i(x)=0,\exists i\in\mathcal{E}\cup\mathcal{I}x:ci(x)=0,∃i∈E∪I都不可微. 考虑到此函数不可微的特殊性, 我们不考虑非光滑优化的技术(比如捆集法(bundle methods)), 而是充分考虑特殊性, 重新设计算法. 在无约束优化中, 我们构建ϕ1(x;μ)\phi_1(x;\mu)ϕ1(x;μ)的简单模型并通过极小化模型得到迭代步. 而这一简单模型可以通过线性化约束cic_ici并以二次函数近似非线性fff得到:q(p;μ)=f(x)+∇f(x)Tp+12pTWp+μ∑i∈E∣ci(x)+∇ci(x)Tp∣+μ∑i∈I[ci(x)+∇ci(x)Tp]−,q(p;\mu)=f(x)+\nabla f(x)^Tp+\frac{1}{2}p^TWp+\mu\sum_{i\in\mathcal{E}}|c_i(x)+\nabla c_i(x)^Tp|+\mu\sum_{i\in\mathcal{I}}[c_i(x)+\nabla c_i(x)^Tp]^-,q(p;μ)=f(x)+∇f(x)Tp+21pTWp+μi∈E∑∣ci(x)+∇ci(x)Tp∣+μi∈I∑[ci(x)+∇ci(x)Tp]−,其中WWW为对称阵, 其通常包含了fff和ci,i∈E∪Ic_i,i\in\mathcal{E}\cup\mathcal{I}ci,i∈E∪I的二阶导数信息. 模型q(p;μ)q(p;\mu)q(p;μ)仍不光滑, 但我们可以添加人工变量, 构造光滑的二次规划问题:min⁡p,r,s,tf(x)+12pTWp+∇f(x)Tp+μ∑i∈E(ri+si)+μ∑i∈Itisubjectto∇ci(x)Tp+ci(x)=ri−si,i∈E,∇ci(x)Tp+ci(x)≥−ti,i∈I,r,s,t≥0.\begin{array}{rl}\min\limits_{p,r,s,t} & f(x)+\frac{1}{2}p^TWp+\nabla f(x)^Tp+\mu\sum\limits_{i\in\mathcal{E}}(r_i+s_i)+\mu\sum\limits_{i\in\mathcal{I}}t_i\\\mathrm{subject\,to} & \nabla c_i(x)^Tp+c_i(x)=r_i-s_i,\quad i\in\mathcal{E},\\& \nabla c_i(x)^Tp+c_i(x)\ge-t_i,\quad i\in\mathcal{I},\\ & r,s,t\ge0.\end{array}p,r,s,tminsubjecttof(x)+21pTWp+∇f(x)Tp+μi∈E∑(ri+si)+μi∈I∑ti∇ci(x)Tp+ci(x)=ri−si,i∈E,∇ci(x)Tp+ci(x)≥−ti,i∈I,r,s,t≥0.此子问题可用标准的二次规划求解器求解. 即使再添加一个盒型(box-shaped)信赖域约束∥p∥∞≤Δ\Vert p\Vert_{\infty}\le\Delta∥p∥∞≤Δ, 它仍是一个二次规划. 这样极小化ϕ1\phi_1ϕ1的方法与逐步二次规划联系紧密, 我们将在下一章讨论之.

选取和更新惩罚因子μk\mu_kμk的策略对迭代的成功与否很关键. 我们曾提到过一种简单但不一定高效的方法, 是先选取初值, 再重复地增大它直到可行性条件满足. 在此方法的一些变体中, 每步迭代惩罚因子的选取满足μk>∥λk∥∞\mu_k>\Vert\lambda_k\Vert_{\infty}μk>∥λk∥∞, 其中λk\lambda_kλk为xkx_kxk处Lagrange乘子的估计. 此策略是基于定理2的, 但由于乘子估计可能不够精确, 且在离解不够近时往往无法得到较好的μk\mu_kμk, 因此这种策略也并不总是能成功.

因选取合适μk\mu_kμk较为困难, 非光滑罚函数法在上世纪九十年代逐渐淡出, 从而刺激了对滤子法的发展. 相比较而言, 滤子法不需要选取惩罚因子. 不过近些年来, 又出现了对罚函数法新一轮的研究兴趣, 部分由于它们能够处理退化的问题. 而新的更新惩罚因子的方法在某些特定的情形下, 能够克服先前的困难. 具体可见下一章的算法5.

在选取初始点xk+1sx_{k+1}^sxk+1s时需要尤其注意. 若μk\mu_kμk对当前合适, 则我们可以设xk+1sx_{k+1}^sxk+1s为ϕ1(x;μk)\phi_1(x;\mu_k)ϕ1(x;μk)的极小点xkx_kxk. 否则应考虑选取更早些迭代的极小点.

2.2 一般的非光滑罚函数法

精确非光滑罚函数可以定义为除了l1l_1l1范数以外的其他范数形式. 我们可以写作ϕ(x;μ)=f(x)+μ∥cE(x)∥+μ∥[cI(x)]−∥,\phi(x;\mu)=f(x)+\mu\Vert c_{\mathcal{E}}(x)\Vert+\mu\Vert[c_{\mathcal{I}}(x)]^-\Vert,ϕ(x;μ)=f(x)+μ∥cE(x)∥+μ∥[cI(x)]−∥,其中∥⋅∥\Vert\cdot\Vert∥⋅∥为任意向量范数, 且所有的等式、不等式约束都分别存放于向量值函数cE,cIc_{\mathcal{E}},c_{\mathcal{I}}cE,cI中. 框架2可用于任何此类罚函数. 我们仅需重新定义不可行性度量为h(x)=∥cE(x)∥+∥[cI(x)]−∥h(x)=\Vert c_{\mathcal{E}}(x)\Vert+\Vert[c_{\mathcal{I}}(x)]^-\Verth(x)=∥cE(x)∥+∥[cI(x)]−∥. 事实上可以证明, 这些范数所对应的不同罚函数的极小点是一一对应的, 证明可见Han S P和Mangasarian O L²的定理4.2. 最常用的范数为l1,l∞,l2l_1,l_{\infty},l_2l1,l∞,l2(无平方)范数. 对l∞l_{\infty}l∞范数容易得到类似于上一小节的线性化重构.

对l1l_1l1罚函数成立的理论性质对一般形式也成立. 在定理3中, 我们把不等式换成μ≥μ∗=∥λ∗∥D,\mu\ge\mu^*=\Vert\lambda^*\Vert_D,μ≥μ∗=∥λ∗∥D,其中∥⋅∥D\Vert\cdot\Vert_D∥⋅∥D是∥⋅∥\Vert\cdot\Vert∥⋅∥的对偶范数. 而定理4则无需修改直接可应用.

下说明此类的罚函数若是精确的必是非光滑的. 为说明简单, 我们仅考虑只有一个等式约束c1(x)=0c_1(x)=0c1(x)=0的情形, 考虑罚函数ϕ(x;μ)=f(x)+μh(ci(x)),\phi(x;\mu)=f(x)+\mu h(c_i(x)),ϕ(x;μ)=f(x)+μh(ci(x)),其中h:R→Rh:\mathbb{R}\to\mathbb{R}h:R→R为满足h(y)≥0,∀y∈R;h(0)=0h(y)\ge0,\forall y\in\mathbb{R};h(0)=0h(y)≥0,∀y∈R;h(0)=0的函数. 若hhh连续可微, 则hhh有一极小点在000. 于是∇h(0)=0\nabla h(0)=0∇h(0)=0. 若x∗x^*x∗为问题的局部解, 则ci(x∗)=0⇒∇h(c1(x∗))=0c_i(x^*)=0\Rightarrow\nabla h(c_1(x^*))=0ci(x∗)=0⇒∇h(c1(x∗))=0. 若x∗x^*x∗为ϕ(x;μ)\phi(x;\mu)ϕ(x;μ)的局部极小点(只需μ\muμ足够大), 则有0=∇ϕ(x∗;μ)=∇f(x∗)+μ∇ci(x∗)∇h(c1(x∗))=∇f(x∗).0=\nabla\phi(x^*;\mu)=\nabla f(x^*)+\mu\nabla c_i(x^*)\nabla h(c_1(x^*))=\nabla f(x^*).0=∇ϕ(x∗;μ)=∇f(x∗)+μ∇ci(x∗)∇h(c1(x∗))=∇f(x∗).然而, fff在约束优化问题中的解一般不是稳定点, 因此hhh连续可微的假设错误, ϕ(⋅;μ)\phi(\cdot;\mu)ϕ(⋅;μ)不可能光滑.

非光滑罚函数也可用作其他机制下计算迭代步时的价值函数. 具体可见第十五章第四节的一般性讨论和第十八和第十九章中的具体实施.

3. 增广Lagrange函数法: 等式约束情形

下面我们讨论一种被称为乘子法the method of multipliers或增广Lagrange函数法augmented Lagragian methods的方法. 此法与第1节中的二次罚函数法相关, 但它通过显式引入Lagrange乘子估计至目标函数, 减小了问题的病态程度. 引入显式Lagrange乘子估计后的函数就被称为增广Lagrange函数(augmented Lagrangian function). 与第2节中讨论的罚函数相反, 增广Lagrange函数极大地保留了光滑性, 因此具体实施可调用标准的无约束或带界约束的优化软件求解.
本节我们在Lagrange乘子估计上, 用上标表示迭代指标, 下标表示分量. 而对于其他的变量, 使用下标表示迭代指标.

3.1 动机和算法框架

我们首先考虑等式约束问题. 二次罚函数Q(x;μ)Q(x;\mu)Q(x;μ)以不可行度量的平方惩罚违反行为. 但我们从定理2中也看到了, Q(x;μk)Q(x;\mu_k)Q(x;μk)的极小点xkx_kxk并不满足可行性条件ci(x)=0,i∈Ec_i(x)=0,i\in\mathcal{E}ci(x)=0,i∈E, 而是近似地有ci(xk)≈−λi∗/μk,∀i∈E.c_i(x_k)\approx-\lambda_i^*/\mu_k,\quad\forall i\in\mathcal{E}.ci(xk)≈−λi∗/μk,∀i∈E.当然, ci(xk)→0,μk→∞c_i(x_k)\to0,\mu_k\to\inftyci(xk)→0,μk→∞, 但考虑到条件数的问题, 我们可能会考虑能否修改函数Q(x;μk)Q(x;\mu_k)Q(x;μk)的形式来避免这种系统上的偏差, 即使得近似极小点无需μk\mu_kμk充分大就差不多满足等式约束ci(x)=0c_i(x)=0ci(x)=0.

增广Lagrange函数LA(x,λ;μ)\mathcal{L}_A(x,\lambda;\mu)LA(x,λ;μ)正满足我们的需求. 基于定理2, 它向目标加入了Lagrange乘子λ\lambdaλ的显式估计, 即有LA(x,λ;μ)=deff(x)−∑i∈Eλici(x)+μ2∑i∈Eci2(x).\mathcal{L}_A(x,\lambda;\mu)\xlongequal{def}f(x)-\sum_{i\in\mathcal{E}}\lambda_ic_i(x)+\frac{\mu}{2}\sum_{i\in\mathcal{E}}c_i^2(x).LA(x,λ;μ)deff(x)−i∈E∑λici(x)+2μi∈E∑ci2(x).以上形式, 与标准的Lagrange函数相比多了平方项, 而与二次罚函数相比则多了关于λ\lambdaλ的求和项. 从这个角度上看, 它是Lagrange函数和二次罚函数的结合.

下面我们设计一种对xxx极小化的算法, 其中固定惩罚因子μ\muμ为μk>0\mu_k>0μk>0, 固定λ\lambdaλ为λk\lambda^kλk. 用xkx_kxk表示LA(x,λk;μk)\mathcal{L}_A(x,\lambda^k;\mu_k)LA(x,λk;μk)的近似极小点, 于是由无约束极小的最优性条件可知0=∇xLA(xk,λk;μk)=∇f(xk)−∑i∈E[λik−μkci(xk)]∇ci(xk).0=\nabla_x\mathcal{L}_A(x_k,\lambda^k;\mu_k)=\nabla f(x_k)-\sum_{i\in\mathcal{E}}[\lambda_i^k-\mu_kc_i(x_k)]\nabla c_i(x_k).0=∇xLA(xk,λk;μk)=∇f(xk)−i∈E∑[λik−μkci(xk)]∇ci(xk).比较约束优化的一阶最优性条件, 若∇ci(xk)\nabla c_i(x_k)∇ci(xk)线性无关, 则有λi∗≈λik−μkci(xk),∀i∈E.\lambda_i^*\approx\lambda_i^k-\mu_kc_i(x_k),\quad\forall i\in\mathcal{E}.λi∗≈λik−μkci(xk),∀i∈E.整理可得ci(xk)≈−1μk(λi∗−λik),∀i∈E,c_i(x_k)\approx-\frac{1}{\mu_k}(\lambda_i^*-\lambda_i^k),\quad\forall i\in\mathcal{E},ci(xk)≈−μk1(λi∗−λik),∀i∈E,从这可看出, 若λk\lambda^kλk接近于最优乘子λ∗\lambda^*λ∗, 则xkx_kxk处的不可行性度量将会小于(1/μk)(1/\mu_k)(1/μk), 而不是成比例(见定理2). 而上面的式子也为我们提供了更新当前估计λk\lambda^kλk的公式:λik+1=λik−μkci(xk),∀i∈E.\lambda_i^{k+1}=\lambda_i^k-\mu_kc_i(x_k),\quad\forall i\in\mathcal{E}.λik+1=λik−μkci(xk),∀i∈E.

于是我们得到下述算法框架.

框架3 (Augmented Lagrangian Method-Equality Constraints)
Given μ0>0\mu_0>0μ0>0, tolerance τ0>0\tau_0>0τ0>0, starting points x0sx_0^sx0s and λ0\lambda^0λ0;
for k=0,1,2,…k=0,1,2,\ldotsk=0,1,2,…
\quad\quadFind an approximate minimizer xkx_kxk of LA(⋅,λk;μk)\mathcal{L}_A(\cdot,\lambda^k;\mu_k)LA(⋅,λk;μk), starting at xksx_k^sxks, and terminating when ∥∇xLA(xk,λk;μk)∥≤τk\Vert\nabla_x\mathcal{L}_A(x_k,\lambda^k;\mu_k)\Vert\le\tau_k∥∇xLA(xk,λk;μk)∥≤τk;
\quad\quadif a convergence test for the original problem is satisfied
\quad\quad\quad\quadstop with approximate solution xkx_kxk;
\quad\quadend (if)
\quad\quadUpdate Lagrange multipliers to obtain λk+1\lambda^{k+1}λk+1;
\quad\quadChoose new penalty perameter μk+1≥μk\mu_{k+1}\ge\mu_kμk+1≥μk;
\quad\quadSet starting point for the next iteration to xk+1s=xkx_{k+1}^s=x_kxk+1s=xk;
\quad\quadSelect tolerance τk+1\tau_{k+1}τk+1;
end (for)

我们通过一个例子来阐释此法可以无需将μ\muμ增加到∞\infty∞即可保证收敛. 于是病态的情况要比框架1有所缓解; 初始点xk+1sx_{k+1}^sxk+1s的选取也就不那么关键; 容忍限τk\tau_kτk可依据不可行性∑i∈E∣c(xk)∣\sum_{i\in\mathcal{E}}|c(x_k)|∑i∈E∣c(xk)∣选取; 若当前迭代下不可行性度量没有充分下降, 则惩罚因子μ\muμ可能就需要增大一些.

例4 再次考虑例1中的问题, 其增广Lagrange函数为LA(x,λ;μ)=x1+x2−λ(x12+x22−2)+μ2(x12+x22−2)2.\mathcal{L}_A(x,\lambda;\mu)=x_1+x_2-\lambda(x_1^2+x_2^2-2)+\frac{\mu}{2}(x_1^2+x_2^2-2)^2.LA(x,λ;μ)=x1+x2−λ(x12+x22−2)+2μ(x12+x22−2)2.如前所述, 原问题解为x∗=(−1,−1)Tx^*=(-1,-1)^Tx∗=(−1,−1)T, 最优Lagrange乘子为λ∗=−0.5\lambda^*=-0.5λ∗=−0.5.
假设在第kkk步迭代, μk=1\mu_k=1μk=1, λk=−0.4\lambda^k=-0.4λk=−0.4. 增广Lagrange函数的等高线可见下图. 与之对比, μk=1\mu_k=1μk=1的二次罚函数等高线图可见图1.

注意到等高线之间的跨度说明此问题的条件数与二次罚函数Q(x;1)Q(x;1)Q(x;1)差不多. 但此时的极小点xk≈(−1.02,−1.02)Tx_k\approx(-1.02,-1.02)^Txk≈(−1.02,−1.02)T则更靠近解x∗=(−1,−1)Tx^*=(-1,-1)^Tx∗=(−1,−1)T. 这个例子表明, 将Lagrange乘子估计项包含进LA(x,λ;μ)\mathcal{L}_A(x,\lambda;\mu)LA(x,λ;μ)中可以得到相较于二次罚函数法的精度提升.

3.2 增广Lagrange函数的性质

下面我们证明关于增广Lagrange函数使用和等式约束问题乘子法的两条结论.

第一条证实, 若我们有精确Lagrange乘子λ∗\lambda^*λ∗的信息, 则原问题的解x∗x^*x∗是LA(x,λ∗;μ),μ\mathcal{L}_A(x,\lambda^*;\mu),\muLA(x,λ∗;μ),μ充分大的严格极小点. 尽管我们一般不能预先知道准确的λ∗\lambda^*λ∗, 但此结论和其证明表明只要λ\lambdaλ是λ∗\lambda^*λ∗的一个较好的估计, 我们可以通过在μ\muμ不是特别大时, 极小化LA(x,λ;μ)\mathcal{L}_A(x,\lambda;\mu)LA(x,λ;μ)得到关于x∗x^*x∗的一个好估计.

定理5 令x∗x^*x∗为原问题的局部解, 且LICQ成立, 对λ=λ∗\lambda=\lambda^*λ=λ∗有二阶充分性条件成立. 则存在μˉ\bar{\mu}μˉ使得对∀μ≥μˉ\forall\mu\ge\bar{\mu}∀μ≥μˉ, x∗x^*x∗为LA(x,λ∗;μ)\mathcal{L}_A(x,\lambda^*;\mu)LA(x,λ∗;μ)的严格局部极小点.

证明: 我们证明x∗x^*x∗对μ\muμ充分大, 满足关于LA(x,λ∗;μ)\mathcal{L}_A(x,\lambda^*;\mu)LA(x,λ∗;μ)的二阶充分性条件, 即∇xLA(x∗,λ∗;μ)=0,∇xx2LA(x∗,λ∗;μ)正定.\nabla_x\mathcal{L}_A(x^*,\lambda^*;\mu)=0,\quad\nabla^2_{xx}\mathcal{L}_A(x^*,\lambda^*;\mu)正定.∇xLA(x∗,λ∗;μ)=0,∇xx2LA(x∗,λ∗;μ)正定.由于x∗x^*x∗为原问题局部解, 且LICQ成立, 于是由KKT条件知∇xL(x∗,λ∗)=0,ci(x∗)=0,∀i∈E\nabla_x\mathcal{L}(x^*,\lambda^*)=0,c_i(x^*)=0,\forall i\in\mathcal{E}∇xL(x∗,λ∗)=0,ci(x∗)=0,∀i∈E, 从而∇xLA(x∗,λ∗;μ)=∇f(x∗)−∑i∈E[λi∗−μci(x∗)]∇ci(x∗)=∇f(x∗)−∑i∈Eλi∗∇ci(x∗)=∇xL(x∗,λ∗)=0,\begin{aligned}\nabla_x\mathcal{L}_A(x^*,\lambda^*;\mu)&=\nabla f(x^*)-\sum_{i\in\mathcal{E}}[\lambda_i^*-\mu c_i(x^*)]\nabla c_i(x^*)\\&=\nabla f(x^*)-\sum_{i\in\mathcal{E}}\lambda_i^*\nabla c_i(x^*)=\nabla_x\mathcal{L}(x^*,\lambda^*)=0,\end{aligned}∇xLA(x∗,λ∗;μ)=∇f(x∗)−i∈E∑[λi∗−μci(x∗)]∇ci(x∗)=∇f(x∗)−i∈E∑λi∗∇ci(x∗)=∇xL(x∗,λ∗)=0,这一点与μ\muμ无关. 对于二阶充分性条件, 定义AAA为约束梯度矩阵, 于是∇xx2LA(x∗,λ∗;μ)=∇xx2L(x∗,λ∗)+μATA.\nabla^2_{xx}\mathcal{L}_A(x^*,\lambda^*;\mu)=\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*)+\mu A^TA.∇xx2LA(x∗,λ∗;μ)=∇xx2L(x∗,λ∗)+μATA.若对于充分大μ\muμ二阶充分性条件不成立, 则对每个k≥1k\ge1k≥1, 我们可选取向量wk:∥wk∥=1w_k:\Vert w_k\Vert=1wk:∥wk∥=1使得0≥wkT∇xx2LA(x∗,λ∗;k)wk=wkT∇xx2L(x∗,λ∗)wk+k∥Awk∥22,0\ge w_k^T\nabla^2_{xx}\mathcal{L}_A(x^*,\lambda^*;k)w_k=w_k^T\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*)w_k+k\Vert Aw_k\Vert^2_2,0≥wkT∇xx2LA(x∗,λ∗;k)wk=wkT∇xx2L(x∗,λ∗)wk+k∥Awk∥22,因此∥Awk∥22≤−(1/k)wkT∇xx2L(x∗,λ∗)wk→0,k→∞.\Vert Aw_k\Vert_2^2\le-(1/k)w_k^T\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*)w_k\to0,\quad k\to\infty.∥Awk∥22≤−(1/k)wkT∇xx2L(x∗,λ∗)wk→0,k→∞.由于{wk}\{w_k\}{wk}位于一紧集中, 于是存在聚点www. 而上面的式子说明Aw=0Aw=0Aw=0. 这说明w∈F(x∗)w\in\mathcal{F}(x^*)w∈F(x∗). 进一步, wkT∇xx2L(x∗,λ∗)wk≤−k∥Awk∥22≤0,w_k^T\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*)w_k\le-k\Vert Aw_k\Vert_2^2\le0,wkT∇xx2L(x∗,λ∗)wk≤−k∥Awk∥22≤0,取极限有wT∇xx2L(x∗,λ∗)w≤0w^T\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*)w\le0wT∇xx2L(x∗,λ∗)w≤0. 这与原问题的二阶充分性矛盾. 证毕.

第二条结论则立意于更实际的情形λ≠λ∗\lambda\ne\lambda^*λ=λ∗. 它给出了LA(x,λ;μ)\mathcal{L}_A(x,\lambda;\mu)LA(x,λ;μ)极小点靠近x∗x^*x∗的条件, 且给出了xkx_kxk和更新的乘子估计λk+1\lambda^{k+1}λk+1分别与x∗,λ∗x^*,\lambda^*x∗,λ∗的误差界.

定理6 设定理5的条件在x∗,λ∗x^*,\lambda^*x∗,λ∗上成立, 且令μˉ\bar{\mu}μˉ为定理5中的阈值. 则存在正标量δ,ϵ,M\delta,\epsilon,Mδ,ϵ,M使得以下断言成立:

对所有满足∥λk−λ∗∥≤μkδ,μk≥μˉ\Vert\lambda^k-\lambda^*\Vert\le\mu_k\delta,\quad\mu_k\ge\bar{\mu}∥λk−λ∗∥≤μkδ,μk≥μˉ的λk,μk\lambda^k,\mu_kλk,μk, 问题min⁡xLA(x,λk;μk),subjectto∥x−x∗∥≤ϵ\min_x\mathcal{L}_A(x,\lambda^k;\mu_k),\quad\mathrm{subject\,to\,}\Vert x-x^*\Vert\le\epsilonxminLA(x,λk;μk),subjectto∥x−x∗∥≤ϵ有唯一解xkx_kxk. 进一步地, ∥xk−x∗∥≤M∥λk−λ∗∥/μk.\Vert x_k-x^*\Vert\le M\Vert\lambda^k-\lambda^*\Vert/\mu_k.∥xk−x∗∥≤M∥λk−λ∗∥/μk.
对所有满足∥λk−λ∗∥≤μkδ,μk≥μˉ\Vert\lambda^k-\lambda^*\Vert\le\mu_k\delta,\quad\mu_k\ge\bar{\mu}∥λk−λ∗∥≤μkδ,μk≥μˉ的λk,μk\lambda^k,\mu_kλk,μk, 我们有∥λk+1−λ∗∥≤M∥λk−λ∗∥/μk,\Vert\lambda^{k+1}-\lambda^*\Vert\le M\Vert\lambda^k-\lambda^*\Vert/\mu_k,∥λk+1−λ∗∥≤M∥λk−λ∗∥/μk,其中λk+1\lambda^{k+1}λk+1由公式λik+1=λik−μkci(xk),∀i∈E\lambda_i^{k+1}=\lambda_i^k-\mu_kc_i(x_k),\quad\forall i\in\mathcal{E}λik+1=λik−μkci(xk),∀i∈E给出.
对所有满足∥λk−λ∗∥≤μkδ,μk≥μˉ\Vert\lambda^k-\lambda^*\Vert\le\mu_k\delta,\quad\mu_k\ge\bar{\mu}∥λk−λ∗∥≤μkδ,μk≥μˉ的λk,μk\lambda^k,\mu_kλk,μk, 矩阵∇xx2LA(xk,λk;μk)\nabla^2_{xx}\mathcal{L}_A(x_k,\lambda^k;\mu_k)∇xx2LA(xk,λk;μk)正定, 约束梯度∇ci(xk),i∈E\nabla c_i(x_k),i\in\mathcal{E}∇ci(xk),i∈E线性无关.

证明: 以下以μ,λ,λ~\mu,\lambda,\tilde\lambdaμ,λ,λ~分别表示μk,λk,λk+1\mu_k,\lambda^k,\lambda^{k+1}μk,λk,λk+1. 对μ>0\mu>0μ>0, 考虑以下以(x,λ~,λ,μ)(x,\tilde\lambda,\lambda,\mu)(x,λ~,λ,μ)为变量的系统:∇f(x)−∇h(x)λ~=0,h(x)+(λ−λ~)/μ=0,\nabla f(x)-\nabla h(x)\tilde\lambda=0,\quad h(x)+(\lambda-\tilde\lambda)/\mu=0,∇f(x)−∇h(x)λ~=0,h(x)+(λ−λ~)/μ=0,其中h(x)=[c1(x)c2(x)⋮cm(x)],∇h(x)=[∇c1(x)∇c2(x)⋯∇cm(x)].h(x)=\begin{bmatrix}c_1(x)\\c_2(x)\\\vdots\\c_m(x)\end{bmatrix},\quad \nabla h(x)=\begin{bmatrix}\nabla c_1(x) & \nabla c_2(x) & \cdots & \nabla c_m(x)\end{bmatrix}.h(x)=⎣⎢⎢⎢⎡c1(x)c2(x)⋮cm(x)⎦⎥⎥⎥⎤,∇h(x)=[∇c1(x)∇c2(x)⋯∇cm(x)].引入变量t∈Rm,γ∈Rt\in\mathbb{R}^m,\gamma\in\mathbb{R}t∈Rm,γ∈R:t=(λ−λ∗)/μ,γ=1/μ,t=(\lambda-\lambda^*)/\mu,\quad\gamma=1/\mu,t=(λ−λ∗)/μ,γ=1/μ,上述系统可写作∇f(x)−∇h(x)λ~=0,h(x)+t+γλ∗−γλ~=0.\nabla f(x)-\nabla h(x)\tilde\lambda=0,\quad h(x)+t+\gamma\lambda^*-\gamma\tilde\lambda=0.∇f(x)−∇h(x)λ~=0,h(x)+t+γλ∗−γλ~=0.对t=0,γ∈[0,1/μˉ]t=0,\gamma\in[0,1/\bar{\mu}]t=0,γ∈[0,1/μˉ], 上述系统有解x=x∗,λ~=λ∗x=x^*,\tilde\lambda=\lambda^*x=x∗,λ~=λ∗. 在此对变量(x,λ~)(x,\tilde\lambda)(x,λ~)的Jacobi矩阵为[∇xx2L(x∗,λ∗)∇h(x∗)−∇h(x∗)T−γI].\begin{bmatrix}\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*) & \nabla h(x^*)\\-\nabla h(x^*)^T & -\gamma I\end{bmatrix}.[∇xx2L(x∗,λ∗)−∇h(x∗)T∇h(x∗)−γI].我们证明, 这一矩阵对所有γ∈[0,1/μˉ]\gamma\in[0,1/\bar\mu]γ∈[0,1/μˉ]可逆. 显然γ=0\gamma=0γ=0时这是成立的. 为证明这对于∀γ∈(0,1/μˉ]\forall\gamma\in(0,1/\bar\mu]∀γ∈(0,1/μˉ]亦成立, 设对某个z∈Rn,w∈Rmz\in\mathbb{R}^n,w\in\mathbb{R}^mz∈Rn,w∈Rm, 我们有[∇xx2L(x∗,λ∗)∇h(x∗)−∇h(x∗)T−γI][zw]=0\begin{bmatrix}\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*) & \nabla h(x^*)\\-\nabla h(x^*)^T & -\gamma I\end{bmatrix}\begin{bmatrix}z\\w\end{bmatrix}=0[∇xx2L(x∗,λ∗)−∇h(x∗)T∇h(x∗)−γI][zw]=0或等价地, ∇xx2L(x∗,λ∗)z+∇h(x∗)w=0,−∇h(x∗)Tz−γw=0.\begin{aligned}\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*)z+\nabla h(x^*)w&=0,\\-\nabla h(x^*)^Tz-\gamma w&=0.\end{aligned}∇xx2L(x∗,λ∗)z+∇h(x∗)w−∇h(x∗)Tz−γw=0,=0.消去www, 得到[∇xx2L(x∗,λ∗)−(1/γ)∇h(x∗)∇h(x∗)T]z=0.[\nabla^2_{xx}\mathcal{L}(x^*,\lambda^*)-(1/\gamma)\nabla h(x^*)\nabla h(x^*)^T]z=0.[∇xx2L(x∗,λ∗)−(1/γ)∇h(x∗)∇h(x∗)T]z=0.对于γ=1/μ,μ≥μˉ\gamma=1/\mu,\mu\ge\bar\muγ=1/μ,μ≥μˉ, 这就得出∇xx2LA(x∗,λ∗;μ)z=0\nabla^2_{xx}\mathcal{L}_A(x^*,\lambda^*;\mu)z=0∇xx2LA(x∗,λ∗;μ)z=0, 而由定理5知∇xx2LA(x∗,λ∗;μ)>0,μ≥μˉ\nabla^2_{xx}\mathcal{L}_A(x^*,\lambda^*;\mu)>0,\mu\ge\bar\mu∇xx2LA(x∗,λ∗;μ)>0,μ≥μˉ, 所以z=0⇒w=0z=0\Rightarrow w=0z=0⇒w=0. 这就说明矩阵对所有γ∈[0,1/μˉ]\gamma\in[0,1/\bar\mu]γ∈[0,1/μˉ]可逆.
由隐函数定理, 存在ϵ>0,δ>0\epsilon>0,\delta>0ϵ>0,δ>0和定义在S(K;δ)S(K;\delta)S(K;δ)(其中K={(0,γ)∣γ∈[0,1/μˉ]}K=\{(0,\gamma)\mid\gamma\in[0,1/\bar\mu]\}K={(0,γ)∣γ∈[0,1/μˉ]})上的连续可微函数x^(t,γ),λ^(t,γ)\hat x(t,\gamma),\hat\lambda(t,\gamma)x^(t,γ),λ^(t,γ)使得(∣x^(t,γ)−x∗∣2+∣λ^(t,γ)−λ∗∣2)1/2<ϵ,∀(t,γ)∈S(K;δ),\left(|\hat{x}(t,\gamma)-x^*|^2+|\hat{\lambda}(t,\gamma)-\lambda^*|^2\right)^{1/2}<\epsilon,\forall(t,\gamma)\in S(K;\delta),(∣x^(t,γ)−x∗∣2+∣λ^(t,γ)−λ∗∣2)1/2<ϵ,∀(t,γ)∈S(K;δ),且满足∇f[x^(t,γ)]−∇h[x^(t,γ)]λ^(t,γ)=0,h[x^(t,γ)]+t+γλ∗−γλ^(t,γ)=0.\begin{aligned}\nabla f[\hat{x}(t,\gamma)]-\nabla h[\hat{x}(t,\gamma)]\hat{\lambda}(t,\gamma)&=0,\\h[\hat{x}(t,\gamma)]+t+\gamma\lambda^*-\gamma\hat\lambda(t,\gamma)&=0.\end{aligned}∇f[x^(t,γ)]−∇h[x^(t,γ)]λ^(t,γ)h[x^(t,γ)]+t+γλ∗−γλ^(t,γ)=0,=0.显然δ,ϵ\delta,\epsilonδ,ϵ的选取可使∇h[x^(t,γ)]\nabla h[\hat{x}(t,\gamma)]∇h[x^(t,γ)]秩为mmm且∇xx2L[x^(t,γ),λ^(t,γ)]−μ∇h[x^(t,γ)]∇h[x^(t,γ)]T>0,∀(t,γ)∈S(K;δ),μ≥μˉ.\nabla^2_{xx}\mathcal{L}[\hat{x}(t,\gamma),\hat{\lambda}(t,\gamma)]-\mu\nabla h[\hat x(t,\gamma)]\nabla h[\hat x(t,\gamma)]^T>0,\forall (t,\gamma)\in S(K;\delta),\mu\ge\bar\mu.∇xx2L[x^(t,γ),λ^(t,γ)]−μ∇h[x^(t,γ)]∇h[x^(t,γ)]T>0,∀(t,γ)∈S(K;δ),μ≥μˉ.对μ≥μˉ,∥λ−λ∗∥≤δμ\mu\ge\bar\mu,\Vert\lambda-\lambda^*\Vert\le\delta\muμ≥μˉ,∥λ−λ∗∥≤δμ, 定义x(λ,μ)=x^(λ−λ∗μ,1μ),λ~(λ,μ)=λ^(λ−λ∗μ,1μ).x(\lambda,\mu)=\hat{x}\left(\frac{\lambda-\lambda^*}{\mu},\frac{1}{\mu}\right),\quad\tilde\lambda(\lambda,\mu)=\hat\lambda\left(\frac{\lambda-\lambda^*}{\mu},\frac{1}{\mu}\right).x(λ,μ)=x^(μλ−λ∗,μ1),λ~(λ,μ)=λ^(μλ−λ∗,μ1).于是对所有满足∥λ−λ∗∥≤μδ,μ≥μˉ\Vert\lambda-\lambda^*\Vert\le\mu\delta,\quad\mu\ge\bar{\mu}∥λ−λ∗∥≤μδ,μ≥μˉ的λk,μk\lambda^k,\mu_kλk,μk, 有∇f[x(λ,μ)]−∇h[x(λ,c)]λ~(λ,μ)=0,λ~(λ,μ)=λ−μh[x(λ,μ)],∇xx2L[x(λ,μ),λ~(λ,μ)]−μ∇h[x(λ,μ)]∇h[x(λ,μ)]T=∇xx2LA[x(λ,μ),λ]>0.\begin{aligned}\nabla f[x(\lambda,\mu)]-\nabla h[x(\lambda,c)]\tilde\lambda(\lambda,\mu)&=0,\\\tilde\lambda(\lambda,\mu)&=\lambda-\mu h[x(\lambda,\mu)],\\\nabla^2_{xx}\mathcal{L}[x(\lambda,\mu),\tilde\lambda(\lambda,\mu)]-\mu\nabla h[x(\lambda,\mu)]\nabla h[x(\lambda,\mu)]^T=\nabla^2_{xx}\mathcal{L}_A[x(\lambda,\mu),\lambda]&>0.\end{aligned}∇f[x(λ,μ)]−∇h[x(λ,c)]λ~(λ,μ)λ~(λ,μ)∇xx2L[x(λ,μ),λ~(λ,μ)]−μ∇h[x(λ,μ)]∇h[x(λ,μ)]T=∇xx2LA[x(λ,μ),λ]=0,=λ−μh[x(λ,μ)],>0.这就证明了第三条. 而为了证明第一条和第二条, 对∇f[x^(t,γ)]−∇h[x^(t,γ)]λ^(t,γ)=0,h[x^(t,γ)]+t+γλ∗−γλ^(t,γ)=0.\begin{aligned}\nabla f[\hat{x}(t,\gamma)]-\nabla h[\hat{x}(t,\gamma)]\hat{\lambda}(t,\gamma)&=0,\\h[\hat{x}(t,\gamma)]+t+\gamma\lambda^*-\gamma\hat\lambda(t,\gamma)&=0.\end{aligned}∇f[x^(t,γ)]−∇h[x^(t,γ)]λ^(t,γ)h[x^(t,γ)]+t+γλ∗−γλ^(t,γ)=0,=0.中的t,γt,\gammat,γ求导, 得到[∇tx^(t,γ)T∇γx^(t,γ)T∇tλ^(t,γ)T∇γλ^(t,γ)T]=A(t,γ)[OO−Iλ^(t,γ)−λ∗],\begin{bmatrix}\nabla_t\hat x(t,\gamma)^T & \nabla_{\gamma}\hat x(t,\gamma)^T\\\nabla_t\hat\lambda(t,\gamma)^T & \nabla_{\gamma}\hat\lambda(t,\gamma)^T\end{bmatrix}=A(t,\gamma)\begin{bmatrix}O & O\\-I & \hat\lambda(t,\gamma)-\lambda^*\end{bmatrix},[∇tx^(t,γ)T∇tλ^(t,γ)T∇γx^(t,γ)T∇γλ^(t,γ)T]=A(t,γ)[O−IOλ^(t,γ)−λ∗],其中A(t,γ)=[∇xx2L[x^(t,γ),λ^(t,γ)]−∇h[x^(t,γ)]∇h[x^(t,γ)]T−γI]−1.A(t,\gamma)=\begin{bmatrix}\nabla^2_{xx}\mathcal{L}[\hat x(t,\gamma),\hat\lambda(t,\gamma)] &- \nabla h[\hat x(t,\gamma)]\\\nabla h[\hat x(t,\gamma)]^T & -\gamma I\end{bmatrix}^{-1}.A(t,γ)=[∇xx2L[x^(t,γ),λ^(t,γ)]∇h[x^(t,γ)]T−∇h[x^(t,γ)]−γI]−1.对∀(t,γ):∣t∣<δ,γ∈[0,1/μˉ]\forall(t,\gamma):|t|<\delta,\gamma\in[0,1/\bar\mu]∀(t,γ):∣t∣<δ,γ∈[0,1/μˉ], 我们有[x^(t,γ)−x∗λ^(t,γ)−λ∗]=[x^(t,γ)−x^(0,0)λ^(t,γ)−λ^(0,0)]=∫01A(ζt,ζγ)[OO−Iλ^(ζt,ζγ)−λ∗][tγ]dζ.\begin{aligned}\begin{bmatrix}\hat x(t,\gamma)-x^*\\\hat\lambda(t,\gamma)-\lambda^*\end{bmatrix}&=\begin{bmatrix}\hat x(t,\gamma)-\hat x(0,0)\\\hat\lambda(t,\gamma)-\hat\lambda(0,0)\end{bmatrix}\\&=\int_0^1A(\zeta t,\zeta\gamma)\begin{bmatrix}O & O\\-I & \hat\lambda(\zeta t,\zeta\gamma)-\lambda^*\end{bmatrix}\begin{bmatrix}t\\\gamma\end{bmatrix}\,\mathrm{d}\zeta.\end{aligned}[x^(t,γ)−x∗λ^(t,γ)−λ∗]=[x^(t,γ)−x^(0,0)λ^(t,γ)−λ^(0,0)]=∫01A(ζt,ζγ)[O−IOλ^(ζt,ζγ)−λ∗][tγ]dζ.由前面证得的矩阵对所有γ∈[0,1/μˉ]\gamma\in[0,1/\bar\mu]γ∈[0,1/μˉ]可逆, 于是对充分小的δ\deltaδ, A(t,γ)A(t,\gamma)A(t,γ)在{(t,γ)∣∣t∣<δ,γ∈[0,1/μˉ]}\{(t,\gamma)\mid|t|<\delta,\gamma\in[0,1/\bar\mu]\}{(t,γ)∣∣t∣<δ,γ∈[0,1/μˉ]}上已知有界. 令Mˉ\bar MMˉ为界, 即∥A(t,γ)∥≤Mˉ,∀∣t∣<δ,γ∈[0,1/μˉ]\Vert A(t,\gamma)\Vert\le\bar M,\forall|t|<\delta,\gamma\in[0,1/\bar\mu]∥A(t,γ)∥≤Mˉ,∀∣t∣<δ,γ∈[0,1/μˉ], 且进一步可取δ\deltaδ充分小以保证Mˉδ<1\bar M\delta<1Mˉδ<1. 于是从上式得到(∣x^(t,γ)−x∗∣2+∣λ^(t,γ)−λ∗∣2)1/2≤Mˉ(∣t∣+max⁡0≤ζ≤1∣λ^(ζt,ζγ)−λ∗∣γ).\begin{aligned}&\left(|\hat x(t,\gamma)-x^*|^2+|\hat\lambda(t,\gamma)-\lambda^*|^2\right)^{1/2}&\le\bar M\left(|t|+\max_{0\le\zeta\le1}|\hat\lambda(\zeta t,\zeta\gamma)-\lambda^*|\gamma\right).\end{aligned}(∣x^(t,γ)−x∗∣2+∣λ^(t,γ)−λ∗∣2)1/2≤Mˉ(∣t∣+0≤ζ≤1max∣λ^(ζt,ζγ)−λ∗∣γ).由此, 对∀(t,γ):∣t∣<δ,γ∈[0,1/μˉ],γ<δ\forall(t,\gamma):|t|<\delta,\gamma\in[0,1/\bar\mu],\gamma<\delta∀(t,γ):∣t∣<δ,γ∈[0,1/μˉ],γ<δ,∣λ^(t,γ)−λ∗∣≤Mˉ∣t∣+Mˉγmax⁡0≤ζ≤1∣λ^(ζt,ζγ)−λ∗∣.|\hat\lambda(t,\gamma)-\lambda^*|\le\bar M|t|+\bar M\gamma\max_{0\le\zeta\le1}|\hat\lambda(\zeta t,\zeta\gamma)-\lambda^*|.∣λ^(t,γ)−λ∗∣≤Mˉ∣t∣+Mˉγ0≤ζ≤1max∣λ^(ζt,ζγ)−λ∗∣.将上式中的t,γt,\gammat,γ换成ζt,ζγ,ζ∈[0,1]\zeta t,\zeta\gamma,\zeta\in[0,1]ζt,ζγ,ζ∈[0,1], 得到max⁡0≤ζ≤1∣λ^(ζt,ζγ)−λ∗)∣≤Mˉ1−Mˉγ∣t∣.\max_{0\le\zeta\le1}|\hat\lambda(\zeta t,\zeta\gamma)-\lambda^*)|\le\frac{\bar M}{1-\bar M\gamma}|t|.0≤ζ≤1max∣λ^(ζt,ζγ)−λ∗)∣≤1−MˉγMˉ∣t∣.组合起来就有, 对∀(t,γ):∣t∣<δ,γ∈[0,1/μˉ],γ<δ\forall(t,\gamma):|t|<\delta,\gamma\in[0,1/\bar\mu],\gamma<\delta∀(t,γ):∣t∣<δ,γ∈[0,1/μˉ],γ<δ,(∣x^(t,γ)−x∗∣2+∣λ^(t,,γ)−λ∗∣2)1/2≤(μ+μ2γ1−μγ)∣t∣≤Mˉ1−Mˉδ∣t∣.\left(|\hat x(t,\gamma)-x^*|^2+|\hat\lambda(t,,\gamma)-\lambda^*|^2\right)^{1/2}\le\left(\mu+\frac{\mu^2\gamma}{1-\mu\gamma}\right)|t|\le\frac{\bar M}{1-\bar M\delta}|t|.(∣x^(t,γ)−x∗∣2+∣λ^(t,,γ)−λ∗∣2)1/2≤(μ+1−μγμ2γ)∣t∣≤1−MˉδMˉ∣t∣.令δ\deltaδ充分小, 我们有(∣x^(t,γ)−x∗∣2+∣λ^(t,,γ)−λ∗∣2)1/2≤2Mˉ∣t∣.\left(|\hat x(t,\gamma)-x^*|^2+|\hat\lambda(t,,\gamma)-\lambda^*|^2\right)^{1/2}\le2\bar M|t|.(∣x^(t,γ)−x∗∣2+∣λ^(t,,γ)−λ∗∣2)1/2≤2Mˉ∣t∣.由ttt的定义, 并令x(λ,μ)=x^(t,γ),λ~(λ,μ)=λ^(t,γ)x(\lambda,\mu)=\hat x(t,\gamma),\tilde\lambda(\lambda,\mu)=\hat\lambda(t,\gamma)x(λ,μ)=x^(t,γ),λ~(λ,μ)=λ^(t,γ), 我们有对∀(λ,μ):∣λ−λ∗∣<δμ,μ>max⁡{μˉ,1/δ}\forall(\lambda,\mu):|\lambda-\lambda^*|<\delta\mu,\mu>\max\{\bar\mu,1/\delta\}∀(λ,μ):∣λ−λ∗∣<δμ,μ>max{μˉ,1/δ}, 有∣x(λ,μ)−x∗∣≤2Mˉ∣λ−λ∗∣/μ,∣λ~(λ,μ)−λ∗∣≤2Mˉ∣λ−λ∗∣/μ.|x(\lambda,\mu)-x^*|\le2\bar M|\lambda-\lambda^*|/\mu,\quad|\tilde\lambda(\lambda,\mu)-\lambda^*|\le2\bar M|\lambda-\lambda^*|/\mu.∣x(λ,μ)−x∗∣≤2Mˉ∣λ−λ∗∣/μ,∣λ~(λ,μ)−λ∗∣≤2Mˉ∣λ−λ∗∣/μ.令M=2MˉM=2\bar MM=2Mˉ即可得第一条和第二条对∀(λ,μ):∣λ−λ∗∣<δμ,μ>max⁡{μˉ,1/δ}\forall(\lambda,\mu):|\lambda-\lambda^*|<\delta\mu,\mu>\max\{\bar\mu,1/\delta\}∀(λ,μ):∣λ−λ∗∣<δμ,μ>max{μˉ,1/δ}是成立的. 而由于x(⋅,⋅)x(\cdot,\cdot)x(⋅,⋅)连续可微, 因此我们也可找到MMM使得第一条和第二条对∀(λ,μ):∣λ−λ∗∣<δμ,≤ˉμ≤max⁡{μˉ,1/δ}\forall(\lambda,\mu):|\lambda-\lambda^*|<\delta\mu,\bar\le\mu\le\max\{\bar\mu,1/\delta\}∀(λ,μ):∣λ−λ∗∣<δμ,≤ˉμ≤max{μˉ,1/δ}成立. 证毕.

此定理解释了增广Lagrange函数法的一些显著的性质.

∥xk−x∗∥≤M∥λk−λ∗∥/μk\Vert x_k-x^*\Vert\le M\Vert\lambda^k-\lambda^*\Vert/\mu_k∥xk−x∗∥≤M∥λk−λ∗∥/μk表示若λk\lambda^kλk精确或惩罚因子μk\mu_kμk很大, 则xkx_kxk就离得x∗x^*x∗很近. 因此这就给出了两种改进xkx_kxk精确度的方法, 而二次罚函数法只给出了一种方法——增大惩罚因子μk\mu_kμk.
∥λk+1−λ∗∥≤M∥λk−λ∗∥/μk\Vert\lambda^{k+1}-\lambda^*\Vert\le M\Vert\lambda^k-\lambda^*\Vert/\mu_k∥λk+1−λ∗∥≤M∥λk−λ∗∥/μk则表明, 局部上我们可以通过选取充分大的μk\mu_kμk保证乘子上精确度的提升.
无约束极小的二阶充分条件在第kkk个子问题下依然是成立的, 因此使用标准的无约束极小技术是有望得到好的结果的.

4. 实用增广Lagrange函数法

本节我们讨论实用的增广Lagrange函数法, 特别地, 是要处理带不等式约束的情形. 我们分别讨论三种方式——界约束重构、约束线性化重构和无约束重构.

4.1 界约束重构

给定一般的非线性规划问题, 我们可以增加松弛变量sis_isi将问题转化为带等式约束的问题, 即ci(x)−si=0,si≥0,∀i∈I.c_i(x)-s_i=0,\quad s_i\ge0,\quad\forall i\in\mathcal{I}.ci(x)−si=0,si≥0,∀i∈I.而界约束l≤x≤ul\le x\le ul≤x≤u则无需变换. 这样, 我们就可以将非线性规划写作min⁡x∈Rnf(x),subjecttoci(x)=0,i=1,2,…,m,l≤x≤u.\min\limits_{x\in\mathbb{R}^n}f(x),\quad\mathrm{subject\,to\,}c_i(x)=0,i=1,2,\ldots,m,l\le x\le u.x∈Rnminf(x),subjecttoci(x)=0,i=1,2,…,m,l≤x≤u.其中下界向量lll的一些分量可能为−∞-\infty−∞; uuu同理.

界约束重构的增广Lagrange函数法the bound-constrained Lagrangian methods, BCL methods只包含等式约束, 即LA(x,λ;μ)=f(x)−∑i=1mλici(x)+μ2∑i=1mci2(x).\mathcal{L}_A(x,\lambda;\mu)=f(x)-\sum_{i=1}^m\lambda_ic_i(x)+\frac{\mu}{2}\sum_{i=1}^mc_i^2(x).LA(x,λ;μ)=f(x)−i=1∑mλici(x)+2μi=1∑mci2(x).而界约束则显式地加在子问题上, 于是有子问题min⁡xLA(x,λ;μ),subjectto≤x≤u.\min_x\mathcal{L}_A(x,\lambda;\mu),\quad\mathrm{subject\,to\,}\le x\le u.xminLA(x,λ;μ),subjectto≤x≤u.在此问题被近似求解后, 我们就更新乘子λ\lambdaλ和惩罚因子μ\muμ, 之后重复.

一种求解带界约束重构非线性规划的技术是非线性投影梯度法, 可见下一章. 考虑子问题的KKT条件, 我们会发现xxx为子问题解的一阶必要条件是x−P(x−∇xLA(x,λ;μ),l,u)=0,x-P(x-\nabla_x\mathcal{L}_A(x,\lambda;\mu),l,u)=0,x−P(x−∇xLA(x,λ;μ),l,u)=0,其中P(g,l,u)P(g,l,u)P(g,l,u)为将g∈Rng\in\mathbb{R}^ng∈Rn投影至长方体[l,u][l,u][l,u]上的投影算子, 其定义为P(g,l,u)i={li,gi≤li,gi,gi∈(li,ui),ui,gi≥ui,i=1,2,…,n.P(g,l,u)_i=\left\{\begin{array}{ll}l_i, & g_i\le l_i,\\g_i, & g_i\in(l_i,u_i),\\u_i, & g_i\ge u_i,\end{array}\right.i=1,2,\ldots,n.P(g,l,u)i=⎩⎨⎧li,gi,ui,gi≤li,gi∈(li,ui),gi≥ui,i=1,2,…,n.以下为算法表述.

算法4 (Bound-Constrained Lagrangian Method)
Choose an initial point x0x_0x0 and initial multipliers λ0\lambda^0λ0;
Choose convergence tolerances η∗\eta_*η∗ and w∗w_*w∗;
Set μ0=10,w0=1μ0,η0=1/μ00.1\mu_0=10,w_0=1\mu_0,\eta_0=1/\mu_0^{0.1}μ0=10,w0=1μ0,η0=1/μ00.1;
for k=0,1,2,…k=0,1,2,\ldotsk=0,1,2,…
\quad\quadFind an approximate solution xkx_kxk of the subproblem such that∥xk−P(xk−∇xLA(xk,λk;μk),l,u)∥≤wk;\Vert x_k-P(x_k-\nabla_x\mathcal{L}_A(x_k,\lambda^k;\mu_k),l,u)\Vert\le w_k;∥xk−P(xk−∇xLA(xk,λk;μk),l,u)∥≤wk;\quad\quadif ∥c(xk)∥≤ηk\Vert c(x_k)\Vert\le\eta_k∥c(xk)∥≤ηk
\quad\quad\quad\quad(test for convergence)
\quad\quad\quad\quadif ∥c(xk)∥≤η∗\Vert c(x_k)\Vert\le\eta_*∥c(xk)∥≤η∗ and ∥xk−P(xk−∇xLA(xk,λk;μk),l,u)∥≤w∗\Vert x_k-P(x_k-\nabla_x\mathcal{L}_A(x_k,\lambda^k;\mu_k),l,u)\Vert\le w_*∥xk−P(xk−∇xLA(xk,λk;μk),l,u)∥≤w∗
\quad\quad\quad\quad\quad\quadstop with approximate solution xkx_kxk;
\quad\quad\quad\quadend (if)
\quad\quad\quad\quad(update multipliers, tighten tolerances)
\quad\quad\quad\quad λk+1=λk−μkc(xk)\lambda^{k+1}=\lambda^k-\mu_kc(x_k)λk+1=λk−μkc(xk);
\quad\quad\quad\quadμk+1=μk\mu_{k+1}=\mu_kμk+1=μk;
\quad\quad\quad\quadηk+1=ηk/μk+10.9\eta_{k+1}=\eta_k/\mu_{k+1}^{0.9}ηk+1=ηk/μk+10.9;
\quad\quad\quad\quadwk+1=wk/μk+1w_{k+1}=w_k/\mu_{k+1}wk+1=wk/μk+1;
\quad\quadelse
\quad\quad\quad\quad(increase penalty parameter, tighten tolerances)
\quad\quad\quad\quad λk+1=λk\lambda^{k+1}=\lambda^kλk+1=λk;
\quad\quad\quad\quad μk+1=100μk\mu_{k+1}=100\mu_kμk+1=100μk;
\quad\quad\quad\quad ηk+1=1/μk+10.1\eta_{k+1}=1/\mu_{k+1}^{0.1}ηk+1=1/μk+10.1;
\quad\quad\quad\quad wk+1=1/μk+1w_{k+1}=1/\mu_{k+1}wk+1=1/μk+1;
\quad\quadend (if)
end (for)

算法主要分岔位于近似求解子问题后, 算法将验证约束违反是否充分下降, 即验证条件∥c(xk)∥≤ηk.\Vert c(x_k)\Vert\le\eta_k.∥c(xk)∥≤ηk.若此条件成立, 则不改变惩罚因子, 而Lagrange乘子估计则按公式更新, 在下一次迭代前减小wk,ηkw_k,\eta_kwk,ηk. 若条件不成立, 则增大惩罚因子保证下一次迭代求解子问题时会更加“关照”约束违反度. 此时就不更新Lagrange乘子, 毕竟主要任务是改进可行性.

算法4中出现的常数0.1,0.9,100在某种程度上可以任取. 软件LANCELOT使用带信赖域的投影梯度法求解界约束非线性子问题. 此时, 投影梯度法构建增广Lagrange函数LA\mathcal{L}_ALA的一个二次模型, 并通过近似求解以下信赖域问题得到迭代步ddd:min⁡d12dT[∇xx2L(xk,λk)+μkAkTAk]d+∇xLA(xk,λk;μk)Tdsubjecttol≤xk+d≤u,∥d∥∞≤Δ,\begin{array}{rl}\min\limits_d & \frac{1}{2}d^T\left[\nabla^2_{xx}\mathcal{L}(x_k,\lambda^k)+\mu_kA_k^TA_k\right]d+\nabla_x\mathcal{L}_A(x_k,\lambda^k;\mu_k)^Td\\\mathrm{subject\,to} & l\le x_k+d\le u,\quad\Vert d\Vert_{\infty}\le\Delta,\end{array}dminsubjectto21dT[∇xx2L(xk,λk)+μkAkTAk]d+∇xLA(xk,λk;μk)Tdl≤xk+d≤u,∥d∥∞≤Δ,其中Ak=A(xk)A_k=A(x_k)Ak=A(xk), Δ\DeltaΔ为信赖域半径. 我们可以将信赖域约束写作界约束−Δe≤d≤Δe-\Delta e\le d\le\Delta e−Δe≤d≤Δe. 求解此子问题的算法每步迭代由两个阶段组成:

第一阶段, 做投影梯度线搜索以决定ddd的哪些分量需设为它们的边界值.
第二阶段, 使用CG极小化子问题, 其中不再牵涉第一阶段中被固定的分量.

这有点类似于第十六章二次规划的投影梯度法. 第二阶段相当于在可行多面体的某个面上做的极小化. 这一算法的好处在于, 无需计算KKT矩阵或者约束Jacobi矩阵AkA_kAk的分解. CG迭代仅需计算矩阵-向量乘积.

需要说明的是, Lagrange函数的Hessian阵∇xx2L(xk,λk)\nabla^2_{xx}\mathcal{L}(x_k,\lambda^k)∇xx2L(xk,λk)可以用基于BFGS或SR1更新公式的拟牛顿近似矩阵替代. LANCELOT包也被设计成可充分利用目标函数与约束中的部分可分结构, 以高效计算Lagrange函数的Hessian阵或拟牛顿近似矩阵, 可见第七章.

4.2 约束线性化重构

约束线性化重构的增广Lagrange函数法linearly constrained Lagragian methods, LCL methods的主要想法是, 在线性化原约束的限制下, 极小化Lagrange函数(或增广Lagrange函数), 产生迭代步. 对于上一小节中的一般非线性规划min⁡x∈Rnf(x),subjecttoci(x)=0,i=1,2,…,m,l≤x≤u,\min_{x\in\mathbb{R}^n}f(x),\quad\mathrm{subject\,to\,}c_i(x)=0,i=1,2,\ldots,m,l\le x\le u,x∈Rnminf(x),subjecttoci(x)=0,i=1,2,…,m,l≤x≤u,LCL方法的子问题就具有以下形式:min⁡xFk(x)subjecttoc(xk)+Ak(x−xk)=0,l≤x≤u.\begin{array}{rl}\min\limits_x & F_k(x)\\\mathrm{subject\,to} & c(x_k)+A_k(x-x_k)=0,\quad l\le x\le u.\end{array}xminsubjecttoFk(x)c(xk)+Ak(x−xk)=0,l≤x≤u.对Fk(x)F_k(x)Fk(x)我们有许多不同的选择. 早期的LCL方法定义Fk(x)=f(x)−∑i=1mλikcˉik(x),F_k(x)=f(x)-\sum_{i=1}^m\lambda_i^k\bar{c}_i^k(x),Fk(x)=f(x)−i=1∑mλikcˉik(x),其中λk\lambda^kλk为当前的Lagrange乘子估计, cˉik(x)\bar{c}_i^k(x)cˉik(x)为ci(x)c_i(x)ci(x)和它在xkx_kxk处线性化的差函数, 即cˉik(x)=ci(x)−ci(xk)−∇ci(xk)T(x−xk).\bar{c}_i^k(x)=c_i(x)-c_i(x_k)-\nabla c_i(x_k)^T(x-x_k).cˉik(x)=ci(x)−ci(xk)−∇ci(xk)T(x−xk).可以证明, 随着xkx_kxk收敛于解x∗x^*x∗, 对应于子问题中等式约束的Lagrange乘子将收敛于最优乘子. 因此, 我们可设λk\lambda^kλk为前一次迭代中对应于等式约束的Lagrange乘子.

现在LCL方法定义FkF_kFk为增广Lagrange函数Fk(x)=f(x)−∑i=1mλikcˉik(x)+μ2∑i=1m[cˉik(x)]2.F_k(x)=f(x)-\sum_{i=1}^m\lambda_i^k\bar{c}_i^k(x)+\frac{\mu}{2}\sum_{i=1}^m[\bar c_i^k(x)]^2.Fk(x)=f(x)−i=1∑mλikcˉik(x)+2μi=1∑m[cˉik(x)]2.这一定义与之前相比, 在实际中具有更好的全局收敛效果.

可以看到, 上述FkF_kFk与之前定义的增广Lagrange函数具有很大的相似度, 而不同在于原本的约束ci(x)c_i(x)ci(x)被替换为了cˉik(x)\bar c_i^k(x)cˉik(x), 后者仅捕捉cic_ici二阶及以上的信息. 线性化约束重构子问题与原本的增广Lagrange函数子问题的不同则在于, 前者要求新迭代点xxx精确满足等式约束的线性化, 同时又将约束的线性部分从目标函数中抽出(这就是用cˉik\bar c_i^kcˉik换cic_ici). 类似于算法4的程序可用于更新惩罚因子μ\muμ以及调整容忍限.

因cˉik(x)\bar c_i^k(x)cˉik(x)在x=xkx=x_kx=xk处梯度为0, 于是∇Fk(xk)=∇f(xk)\nabla F_k(x_k)=\nabla f(x_k)∇Fk(xk)=∇f(xk). 我们也可验证FkF_kFk的Hessian与Lagrange函数或增广Lagrange函数的Hessian也是紧密相关的. 基于这些性质, 线性化重构子问题就类似于第十八章中的SQP子问题.

MINOS包使用增广形式的模型函数, 用既约投影梯度法求解线性化子问题, 其中对FkF_kFk的既约Hessian用拟牛顿近似. 为确保等式约束Lagrange乘子的精度, MINOS中每一次子问题的求解需要更多的目标函数值和约束函数(以及它们的梯度)值. 这要多于SQP算法或内点法. 不过, 这样一来所需求解的子问题数一般比其他方法要少一些.

4.3 无约束化重构

利用近邻点proximal point方法, 我们可以得到带不等式约束问题的无约束化重构增广Lagrange函数子问题. 为说明简便, 假设问题无等式约束, 于是我们可以将原问题等价地写作无约束优化问题min⁡x∈RnF(x),\min\limits_{x\in\mathbb{R}^n}F(x),x∈RnminF(x),其中F(x)=max⁡λ≥0{f(x)−∑i∈Iλici(x)}={f(x),xisfeasible,∞,otherwise.F(x)=\max_{\lambda\ge0}\left\{f(x)-\sum_{i\in\mathcal{I}}\lambda_ic_i(x)\right\}=\left\{\begin{array}{ll}f(x), & x\,is\,feasible,\\\infty, & otherwise.\end{array}\right.F(x)=λ≥0max{f(x)−i∈I∑λici(x)}={f(x),∞,xisfeasible,otherwise.由于FFF不光滑, 直接极小化FFF并不实际. 于是再次想到用光滑的近似模型. 我们以F^(x;λk,μk)\hat F(x;\lambda^k,\mu_k)F^(x;λk,μk)代替FFF, 前者依赖于惩罚因子μk\mu_kμk和Lagrange乘子估计λk\lambda^kλk. 其定义为F^(x;λk,μk)=max⁡λ≥0{f(x)−∑i∈Iλici(x)−12μk∑i∈I(λi−λik)2}.\hat F(x;\lambda^k,\mu_k)=\max_{\lambda\ge0}\left\{f(x)-\sum_{i\in\mathcal{I}}\lambda_ic_i(x)-\frac{1}{2\mu_k}\sum_{i\in\mathcal{I}}(\lambda_i-\lambda_i^k)^2\right\}.F^(x;λk,μk)=λ≥0max{f(x)−i∈I∑λici(x)−2μk1i∈I∑(λi−λik)2}.其中最后一项会对λ\lambdaλ任何偏离前一次估计λk\lambda^kλk的移动做出惩罚, 即, 它会使新的乘子λ\lambdaλ与前一次估计λk\lambda^kλk尽可能地相近. 因F^\hat FF^本身的定义就是一个关于λ\lambdaλ的界约束二次规划, 因此我们有显式解λi={0,−ci(x)+λik/μk≤0,λik−μkci(x),otherwise.\lambda_i=\left\{\begin{array}{ll}0, & -c_i(x)+\lambda_i^k/\mu_k\le0,\\\lambda_i^k-\mu_kc_i(x), & otherwise.\end{array}\right.λi={0,λik−μkci(x),−ci(x)+λik/μk≤0,otherwise.注意到这可以用来更新乘子. 将此代入F^\hat FF^得到F^(x;λk,μk)=f(x)+∑i∈Iψ(ci(x),λik;μk),\hat F(x;\lambda^k,\mu_k)=f(x)+\sum_{i\in\mathcal{I}}\psi(c_i(x),\lambda_i^k;\mu_k),F^(x;λk,μk)=f(x)+i∈I∑ψ(ci(x),λik;μk),其中函数ψ\psiψ具有三个标量参数:ψ(t,σ;μ)=def{−σt+μ2t2,t−σ/μ≤0,−12μσ2,otherwise.\psi(t,\sigma;\mu)\xlongequal{def}\left\{\begin{array}{ll}-\sigma t+\frac{\mu}{2}t^2, & t-\sigma/\mu\le0,\\-\frac{1}{2\mu}\sigma^2, & otherwise.\end{array}\right.ψ(t,σ;μ)def{−σt+2μt2,−2μ1σ2,t−σ/μ≤0,otherwise.因此, 对xxx极小化F^(x;λk,μk)\hat F(x;\lambda^k,\mu_k)F^(x;λk,μk)的xkx_kxk, 并依前面的显式表达更新Lagrange乘子. 比较框架3, 我们会发现FFF扮演了LA\mathcal{L}_ALA的角色. 本节中介绍的格式就是等式约束增广Lagrange函数法向不等式约束情形的推广. 不过, 由于其性质还未得到验证, 此无约束化重构并不如界约束重构和约束线性化重构一般, 具有常用软件包的内嵌.

5. 不同算法的比较与软件介绍

在约束数目较小时, 人们通常使用二次罚函数法. 事实上, 有时仅需对一个较大的μ\muμ做一次Q(x;μ)Q(x;\mu)Q(x;μ)的极小化. 若μ\muμ选得不够好, 得到的解就可能不会很精确. 由于求解无约束优化的主要软件并不内嵌二次罚函数法, 因此很少有关于更新惩罚因子、调整容忍限以及选取初始点的讨论.

尽管二次罚函数更直观、简洁, 但第3节与第4节的增广Lagrange函数法通常更受青睐. 子问题的求解通常并不太难, 且乘子估计的引入避免了病态子问题的出现. 不过二次罚函数仍然作为许多算法正则化(regularization)的一种重要手段, 譬如SQP.

一般情形的l1l_1l1罚函数法由Fletcher在上世纪八十年代提出. 由于它与SQP算法的共同点, 它也被称为Sl1l_1l1QP算法.最近, 软件包KNITRO加入了使用线性规划子问题的l1l_1l1罚函数法. 这两种方法我们将在下一章讨论.

近些年l1l_1l1罚函数受到了大量的关注. 它已被成功用于求解一些困难的问题, 例如带互补约束的数学规划(mathematical programs with complementarity constraints, MPCCs), 其中约束不满足标准的约束规范. 将这些问题约束吸收为罚项而不是线性化, 之后后再使用诸如SQP或内点法, 我们能够扩展这些其他方法的应用面. 软件包SNOPT在SQP算法中使用l1l_1l1罚函数法作为一种防护策略, 以防二次模型不可行或无界, 亦或有无界乘子.

增广Lagrange函数法由于其简洁性, 已受多年关注. MINOS和LANCELOT则是实施增广Lagrange函数法最好的软件包. 它们也适用于大型非线性规划问题. 一般层面上, MINOS的LCL和LANCELOT的BCL有共同的性质. 但它们在计算迭代步的子问题和求解子问题的技术上具有较大差异:

MINOS使用既约空间法处理线性化约束, 并且使用(稠密)拟牛顿近似代替Lagrange函数的Hessian. 因此, MINOS对具有较小自由度的问题是最合适的.
LANCELOT则更适用于约束较少的情形. 正如第4节所说, LANCELOT无需约束Jacobi矩阵AAA的分解, 这也增强了它在大型问题上的应用能力. 它同时也提供多种Hessian近似和预处理子的选择.
PENNON软件包也基于增广Lagrange函数法, 它在处理半定矩阵约束时较有优势.

界约束重构和无约束重构Lagrange函数法的缺陷是, 它们增加了约束的平方, 让问题变得更复杂. 这时, 我们只有通过极小化增广Lagrange函数才能获得可行性的改善. 相反, LCL则在约束上计算牛顿类的迭代步, 稳步改善可行性. 因此, 不难预见, 在带线性约束的问题上MINOS比LANCELOT更具优势.

从第3节的增广Lagrange函数也可构造光滑精确罚函数(smooth exact penalty functions), 但这些就过于复杂了. 例如我们提过等式约束的Fletcher光滑精确罚函数:ϕF(x;μ)=f(x)−λ(x)Tc(x)+μ2∑i∈Eci(x)2.\phi_F(x;\mu)=f(x)-\lambda(x)^Tc(x)+\frac{\mu}{2}\sum_{i\in\mathcal{E}}c_i(x)^2.ϕF(x;μ)=f(x)−λ(x)Tc(x)+2μi∈E∑ci(x)2.其中Lagrange乘子估计λ(x)\lambda(x)λ(x)为最小二乘近似:λ(x)=[A(x)A(x)T]−1A(x)∇f(x)=(AT)†∇f(x).\lambda(x)=[A(x)A(x)^T]^{-1}A(x)\nabla f(x)=\left(A^T\right)^{\dagger}\nabla f(x).λ(x)=[A(x)A(x)T]−1A(x)∇f(x)=(AT)†∇f(x).这里的ϕF\phi_FϕF光滑且精确. ϕF\phi_FϕF的缺点在于需要计算λ(x)\lambda(x)λ(x). 当A(x)A(x)A(x)不满秩时λ(x)\lambda(x)λ(x)不是唯一的, 且当A(x)A(x)A(x)接近奇异时λ\lambdaλ的估计可能会变得很差.

Han S P , Mangasarian O L . Exact penalty functions in nonlinear programming[J]. Mathematical Programming, 1979, 17(1):251-269. ↩︎
Han S P , Mangasarian O L . Exact penalty functions in nonlinear programming[J]. Mathematical Programming, 1979, 17(1):251-269. ↩︎

Numerical Optimization Ch17. Penalty and Augmented Lagragian Methods相关推荐

【Numerical Optimization】17 Penalty and Augmented Lagrangian Methods
2019.02.25 [ 罚方法:两个重要定理.的取值问题以及在ALM中的改善罚方法中的坏条件] 求解思想是将约束函数罚到原函数上形成新的无约束优化问题: 二次惩罚+非平滑惩罚+增广拉格朗日 1.二 ...
Numerical Optimization - my afterword
历时六个月, 从第一篇(2018.9.30)到第十九篇(2019.3.17), 感谢各位博友的支持. 就个人而言, 其实这本书早在2018.11就看完了. 写博客纯粹是为了加深自己的印象.锻炼自己的英 ...
数值优化（Numerical Optimization）学习系列-惩罚和增广拉格朗日方法（Augmented Lagrangian Methods）
原文地址为: 数值优化(Numerical Optimization)学习系列-惩罚和增广拉格朗日方法(Augmented Lagrangian Methods) 概述求解带约束的最优化问题,一类很 ...
Numerical Optimization和Convex optimization 两本书的选择?
Numerical Optimization和Convex optimization 两本书的选择? - 知乎https://www.zhihu.com/question/49689245 Numer ...
Note of Numerical Optimization Ch.3
目录 Numerical Optimization Ch.3 Line Search Methods Step Length Convergence of Line Search Methods Ra ...
数值优化（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 ...
Self-adaptive Differential Evolution Algorithm for Numerical Optimization
0.论文背景本文在DE算法的基础上,提出了一种自适应的差分进化算法(SaDE).其中学习策略的选择和两个控制参数F和CR不需要预先指定. Qin A K, Suganthan P N. Self-a ...

Numerical Optimization Ch17. Penalty and Augmented Lagragian Methods

第十七章: 惩罚函数法与增广Lagrange函数法

文章目录

1. 二次罚函数法

1.1 动机

1.2 算法框架

1.3 二次罚函数法的收敛性

1.4 病态与重构

2. 非光滑精确罚函数法

2.1 实用l1l_1l1罚函数法

2.2 一般的非光滑罚函数法

3. 增广Lagrange函数法: 等式约束情形

3.1 动机和算法框架

3.2 增广Lagrange函数的性质

4. 实用增广Lagrange函数法

4.1 界约束重构

4.2 约束线性化重构

4.3 无约束化重构

5. 不同算法的比较与软件介绍

Numerical Optimization Ch17. Penalty and Augmented Lagragian Methods相关推荐

最新文章

热门文章

Numerical Optimization Ch17. Penalty and Augmented Lagragian Methods

第十七章: 惩罚函数法与增广Lagrange函数法

文章目录

1. 二次罚函数法

1.1 动机

1.2 算法框架

1.3 二次罚函数法的收敛性

1.4 病态与重构

2. 非光滑精确罚函数法

2.1 实用l1l_1l1​罚函数法

2.2 一般的非光滑罚函数法

3. 增广Lagrange函数法: 等式约束情形

3.1 动机和算法框架

3.2 增广Lagrange函数的性质

4. 实用增广Lagrange函数法

4.1 界约束重构

4.2 约束线性化重构

4.3 无约束化重构

5. 不同算法的比较与软件介绍

Numerical Optimization Ch17. Penalty and Augmented Lagragian Methods相关推荐

最新文章

热门文章

2.1 实用l1l_1l1罚函数法