约束优化问题--乘子法

乘子法

惩罚函数法再分析

主要缺点: 当惩罚因子$M$增大, 无约束问题的目标函数$F(\mathbf{x},M)$的Hesse矩阵越来越病态, 使得无约束优化方法的计算难以进行. 下面分析$M$必须无限增大的原因.

考虑等式约束优化问题
$$\begin{array}{rl}& \min f(\mathbf{x})\\ & s.t.h_i(\mathbf{x})=0,i=1,2,...,p.\end{array}$$
引入罚函数
$$F(\mathbf{x},M)=f(\mathbf{x})+\sum _{i=1}^p[h_i(\mathbf{x}){]}^2,$$
若${\mathbf{x}}^{\ast }$是问题的最优解, 则$h_i({\mathbf{x}}^{\ast })=0,i=1,2,...,p$. 于是有
$${\nabla }_{x}F({\mathbf{x}}^{\ast },M)=\nabla f({\mathbf{x}}^{\ast })+2M\sum _{i=1}^p{h}_i({\mathbf{x}}^{\ast })\nabla h_i({\mathbf{x}}^{\ast })=\nabla f({\mathbf{x}}^{\ast }).$$
另一方面, 考虑上述约束问题的$Lagrange$函数
$$L(\mathbf{x},\mu )=f(\mathbf{x})+\sum _{i=1}^p{\mu }_i{h}_i(\mathbf{x}),$$
由KKT条件有
$${\nabla }_{x}L({\mathbf{x}}^{\ast },{\mu }^{\ast })=\nabla f({\mathbf{x}}^{\ast })+\sum _{i=1}^p{\mu }_i\nabla h_i({\mathbf{x}}^{\ast })=0.$$
则可知在通常情况下, $\nabla f({\mathbf{x}}^{\ast })\ne 0$. 只有当${\mathbf{x}}^{\ast }$是$f(\mathbf{x})$无约束问题稳定点时, $\nabla f({\mathbf{x}}^{\ast })=0$才成立. 因此, 找不到一个有限的$M$, 使得${\nabla }_{x}F({\mathbf{x}}^{\ast },M)=0$成立, 只能期望${\lim }_{M\to +\infty }{\nabla }_{x}F({\mathbf{x}}^{\ast },M)=0$.

以上分析过程启发我们将$Lagrange$函数与惩罚项相结合, 构造增广$Lagrange$函数, 也就是接下来要介绍的乘子法.

Hestenes乘子法

Hestenes首先引入增广$Lagrange$函数
$$\begin{array}{rl}\phi (\mathbf{x},\mathbf{\mu })& =f(\mathbf{x})+\sum _{i=1}^p{\mu }_i{h}_i(x)+\frac{c}{2}\sum _{i=1}^p[h_i(x){]}^2\\ & =:L(\mathbf{x},\mathbf{\mu })+\frac{c}{2}\sum _{i=1}^p[h_i(x){]}^2.\end{array}$$
其中$\mathbf{\mu }=[{\mu }_1,...,{\mu }_p]$为$Lagrange$乘子, 罚因子被改写为$\frac{c}{2}$.

由最优性一阶必要条件可知, $\exists {\mathbf{\mu }}^{\ast },s.t.({\mathbf{x}}^{\ast },{\mathbf{\mu }}^{\ast })$为$L(\mathbf{x},\mathbf{\mu })$的稳定点, 即${\nabla }_{x}L({\mathbf{x}}^{\ast },{\mathbf{\mu }}^{\ast })=0$. 又因为附加项$\frac{c}{2}{\sum }_{i=1}^p[h_i(x){]}^2$在$x^{\ast }$梯度为0, 因此有
$${\nabla }_x\phi ({\mathbf{x}}^{\ast },{\mathbf{\mu }}^{\ast })=0.$$
说明$x^{\ast }$是$\phi ({\mathbf{x}}^{\ast },{\mathbf{\mu }}^{\ast })$稳定点. 则问题进一步转化为, 对某个${\mathbf{\mu }}^{\ast }$, 求$\phi ({\mathbf{x}}^{\ast },{\mathbf{\mu }}^{\ast })$的无约束极小点.

求解思路: 首先给定一个足够大的正数$c$, 在迭代过程中逐步调整${\mathbf{\mu }}^{(k)}$.

假设已经求得${\mathbf{\mu }}^{(k)}$, 求解$\min \phi (\mathbf{x},{\mathbf{\mu }}^{(k)})$得最优解${\mathbf{x}}^{(k)}$, 则有
$${\nabla }_x\phi ({\mathbf{x}}^{(k)},{\mathbf{\mu }}^{(k)})=\nabla f({\mathbf{x}}^{(k)})+\sum _{i=1}^p{\mu }_i^{(k)}\nabla h_i({\mathbf{x}}^{(k)})+c\sum _{i=1}^p{h}_i({\mathbf{x}}^{(k)})\nabla h_i({\mathbf{x}}^{(k)})=0.$$
即
$$\nabla f({\mathbf{x}}^{(k)})+\sum _{i=1}^p({\mu }_i^{(k)}+c{h}_i({\mathbf{x}}^{(k)})\nabla h_i({\mathbf{x}}^{(k)})=0.$$
将其与$\nabla L({\mathbf{x}}^{(k)},{\mathbf{\mu }}^{(k)})$比较可得
$${\mu }_i^{(k+1)}={\mu }_i^{(k)}+c{h}_i({\mathbf{x}}^{(k)}),i=1,...,p.$$
迭代终止条件？

**Theorem ** 设${\mathbf{x}}^{(k)}$是$\min \phi (\mathbf{x},{\mathbf{\mu }}^{(k)})$的最优解, 则${\mathbf{x}}^{(k)}$是原约束优化问题的最优解$\iff h_i({\mathbf{x}}^{(k)})=0,i=1,...,p.$

因此, 迭代终止条件可取为: 给定计算精度$\epsilon >0,$ 若${\mathbf{x}}^{(k)}$满足$\Vert h({\mathbf{x}}^{(k)})\Vert <\epsilon$, 则停止迭代, 输出最优解${\mathbf{x}}^{(k)}$. 此外, 若迭代过程中, ${\mathbf{\mu }}^{(k)}$收敛太慢, 可适当增大因子$c$.

算法流程

步骤1 给定初始点$x^{(0)}$, 初始乘子向量${\mathbf{\mu }}^{(1)}$(可取${\mathbf{\mu }}^{(1)}=0$), $\epsilon >0,c>0,0<r<1,\alpha >1$, 令$k=1$;

步骤2 以${\mathbf{x}}^{(k-1)}$为初始点, 求解$\min \phi (\mathbf{x},{\mathbf{\mu }}^{(k)})$, 得${\mathbf{x}}^{(k)}$;

步骤3 若$\Vert h({\mathbf{x}}^{(k)})\Vert <\epsilon$, 迭代终止, ${\mathbf{x}}^{(k)}$为最优解; 否则, 计算$\beta =\frac{\Vert h({\mathbf{x}}^{(k)})\Vert }{\Vert {\mathbf{x}}^{(k-1)}\Vert }$, 若$\beta >r$, 则$c=\alpha c$;

步骤4 ${\mu }_i^{(k+1)}={\mu }_i^{(k-1)}+c{h}_i({\mathbf{x}}^{(k)}),i=1,...,p+1$, 令$k=k+1$, 返回步骤1.

Rockafellar乘子法

Rockafellar将乘子法推广到不等式约束的优化问题
$$\begin{array}{rl}& \min f(x)\\ & s.t.g_i(x)\ge 0,i=1,...,p.\end{array}$$

通过引入松弛变量$z_i(i=1,...,p)$, 将不等式约束转化为等式约束
$$g_i(x)-z_i^2=0,i=1,...,p.$$
考虑其对应的增广$Lagrange$函数
$$\Phi (\mathbf{x},\mathbf{z},\mathbf{\mu })=f(\mathbf{x})+\sum _{i=1}^p{\mu }_i[g_i(\mathbf{x})-z_i^2]+\frac{c}{2}\sum _{i=1}^p[g_i(\mathbf{x})-z_i^2{]}^2.$$
将$\Phi (\mathbf{x},\mathbf{z},\mathbf{\mu })$关于$\mathbf{z}$求极小,
$$\phi (\mathbf{x},\mathbf{\mu })=\underset{\mathbf{z}}{\min }\Phi (\mathbf{x},\mathbf{z},\mathbf{\mu }),{\nabla }_{\mathbf{z}}\Phi (\mathbf{x},\mathbf{z},\mathbf{\mu })=0,$$
得
$$z_i[c{z}_i^2-({\mu }_i+c{g}_i(\mathbf{x}))]=0,i=1,...,p.$$
若${\mu }_i+c{g}_i(\mathbf{x})\le 0$, 则$z_i^2=0$; 若${\mu }_i+c{g}_i(\mathbf{x})>0$, 则$z_i^2=\frac{1}{c}({\mu }_i+c{g}_i(x))$, 即
$$g_i(\mathbf{x})-z_i^2=\{\begin{array}{l}{g}_i(\mathbf{x}),{\mu }_i+c{g}_i(\mathbf{x})\le 0,\\ -\frac{{\mu }_i}{c},{\mu }_i+c{g}_i(\mathbf{x})>0.\end{array}$$
当${\mu }_i+c{g}_i(\mathbf{x})\le 0$时,
$${\mu }_i[g_i(\mathbf{x})-z_i^2]+\frac{c}{2}[g_i(\mathbf{x})-z_i^2{]}^2={\mu }_i{g}_i(\mathbf{x})+\frac{c}{2}[g_i(\mathbf{x}){]}^2=\frac{1}{2c}[({\mu }_i+c{g}_i(\mathbf{x}){)}^2-{\mu }_i^2],$$
当${\mu }_i+c{g}_i(\mathbf{x})>0$时,
$${\mu }_i[g_i(\mathbf{x})-z_i^2]+\frac{c}{2}[g_i(\mathbf{x})-z_i^2{]}^2=-\frac{{\mu }_i^2}{c}+\frac{c}{2}(-\frac{{\mu }_i}{c}{)}^2=-\frac{1}{2c}{\mu }_i^2.$$
故
$$\phi (\mathbf{x},\mathbf{\mu })=f(\mathbf{x})+\frac{1}{2c}\sum _{i=1}^p\{[\min (0,{\mu }_i+c{g}_i(\mathbf{x})){]}^2-{\mu }_i^2\}.$$
可得乘子迭代公式
μ i ( k + 1 ) = min ⁡ { 0 , μ i ( k ) + c g i ( x ( k ) ) } , \mu_i^{(k+1)}=\min\{0,\mu_i^{(k)}+cg_i(\boldsymbol x^{(k)})\}, μi(k+1)​=min{0,μi(k)​+cgi​(x(k))},
迭代终止准则
$$(\sum _{i=1}^p[\min (g_i(\mathbf{x}),-\frac{{\mu }_i}{c}){]}^2{)}^{\frac{1}{2}}<\epsilon .$$

应用举例

例1 等式约束问题

求解
$$\begin{array}{rl}& \min f(x_1,x_2)=x_1^2+x_2^2\\ & s.t.x_1+x_2-2=0.\end{array}$$
易知其最优解为$(1,1{)}^T$.

解: 增广$Lagrange$函数为
$$\phi (x_1,x_2,\mu ,c)=x_1^2+x_2^2+\mu (x_1+x_2-2)+\frac{c}{2}(x_1+x_2-2{)}^2.$$
令$\frac{\partial \phi }{\partial x_1}=\frac{\partial \phi }{\partial x_2}=0$, 得
$$x_1=x_2=\frac{2c-\mu }{2c+2}.$$
则有
$${\mu }^{(k+1)}=\frac{1}{c+1}{\mu }^{(k)}-\frac{2c}{c+1}.$$
取$c=10,$ 则${\mu }^{(k+1)}=\frac{1}{11}{\mu }^{(k)}-\frac{20}{11}$. 设${\mu }^{(k)}\to {\mu }^{\ast }$, 可得
$${\mu }^{\ast }=\frac{1}{11}{\mu }^{\ast }-\frac{20}{11}\Rightarrow u^{\ast }=-2.$$
则原问题最优解${\mathbf{x}}^{\ast }=(1,1{)}^T.$

例2 不等式约束优化问题

求解
$$\begin{array}{rl}& \min f(\mathbf{x})=x_1^2+x_2^2\\ & s.t.x_1-1\ge 0.\end{array}$$
易知最优解为$(1,0{)}^T.$

解: 增广$Lagrange$函数为
$$\begin{gathered} \phi (x_1,x_2,\mu ,c)=x_1^2+x_2^2+\frac{1}{2c}\{[\min (0,\mu +c(x_1-1)){]}^2-{\mu }^2\} \\ =\{\begin{array}{ll}& x_1^2+x_2^2+\mu (x_1-1)+\frac{c}{2}(x_1-1{)}^2,x_1\le 1-\frac{\mu }{c},\\ & x_1^2+x_2^2-\frac{{\mu }^2}{2c},x_1>1-\frac{\mu }{c}.\end{array} \end{gathered}$$

$$\frac{\partial \phi }{\partial x_1}=\{\begin{array}{l}2x_1+\mu +c(x_1-1),x_1\le 1-\frac{\mu }{c}\\ 2x_1,x_1>1-\frac{\mu }{c}\end{array},\frac{\partial \phi }{\partial x_2}=2x_2.$$

令$\nabla \phi (\mathbf{x},\mathbf{\mu })=0$, 得极小点
$$x_1=\frac{c-\mu }{c+2},x_2=0.$$
取$c=4,$有
x ( k ) = ( 4 − μ ( k ) 6 , 0 ) T , μ ( k + 1 ) = min ⁡ { 0 , μ ( k ) + 4 ( x 1 ( k ) − 1 ) } = μ ( k ) − 4 3 . \begin{align*} & \boldsymbol x^{(k)}=(\frac{4-\mu^{(k)}}{6},0)^T,\\ & \mu^{(k+1)}=\min\{0,\mu^{(k)}+4(x_1^{(k)}-1)\}=\frac{\mu^{(k)}-4}{3}. \end{align*} ​x(k)=(64−μ(k)​,0)T,μ(k+1)=min{0,μ(k)+4(x1(k)​−1)}=3μ(k)−4​.​
令$k\to +\infty ,{\mu }^{(k)}\to -2,{\mathbf{x}}^{(k)}\to (1,0{)}^T=x^{\ast }$.

Code

以例1中的等式约束问题为练习, 给出简易MATLAB代码, 可以看出数值结果与理论相符.

% 等式约束的乘子法
function multph
global c miu
x0=[2,1]';
x1=x0;
miu=0;
eps=10^(-4);
c=10;
r=0.4;
a=2;
k=1;
while norm(h(x0))>=eps
    x1=fminsearch('phi',x0);
    beta=norm(h(x1))/norm(h(x0));
    x0=x1;
    if beta>r
        c=a*c;
    end
    miu=miu+c*h(x1);
    k=k+1;
end
format long
disp('最优解')
disp(x1)
disp('乘子')
disp(miu)
disp('迭代次数')
disp(k)
disp('终止条件')
disp(h(x0))
end

function r=h(x)
% 等式约束对应函数
r=x(1)+x(2)-2;
end

function r=phi(x)
global miu c
r=x(1)^2+x(2)^2+miu*(x(1)+x(2)-2)+c/2*(x(1)+x(2)-2)^2;
end

最优解
   1.000053801333254
   0.999968661617259

乘子
  -1.999646326786757

迭代次数
     6

终止条件
     2.246295051211433e-05