上海理工大学学报  2022, Vol. 44 Issue (4): 373-380   PDF    
多面体约束非光滑复合函数的序列有效集方法
时闪闪, 宇振盛     
上海理工大学 理学院,上海 200093
摘要: 对带多面体约束的非光滑复合函数问题的求解进行了研究。针对非光滑复合函数问题,首先,构造光滑函数来逼近非光滑目标函数,通过求解光滑近似问题来达到求解原问题的目的。在此基础上,考虑多面体约束的特殊结构,运用序列二次规划算法的思想,利用有效集策略,通过逐次求解一系列仅含等式约束的二次规划问题来逼近搜索方向的最优解,再通过线搜索求得步长,进而得到下一步的迭代点。最后,从理论上证明了算法的全局收敛性,并进行了初步的数值实验。将该算法与光滑序列投影收缩算法作对比,结果表明,该算法在迭代次数和计算时间上都有一定的优势。
关键词: 光滑近似     多面体约束     复合函数     序列二次规划     有效集    
A sequential active set algorithm for nonsmooth composite functions with polyhedral constraints
SHI Shanshan, YU Zhensheng     
College of Science, University of Shanghai for Science and Technology, Shanghai 200093, China
Abstract: The solution of nonsmooth composite function problem with polyhedral constraints was studied. For the problem of nonsmooth composite function, a smooth function was firstly constructed to approximate the nonsmooth objective function, and the purpose of solving the original problem was then achieved by handling the smooth approximation problem. On this basis, considering the special structure of polyhedral constraints, the idea of sequential quadratic programming algorithm was adopted. A quadratic programming subproblem was solved to obtain the search direction. The active set strategy was adopted by the solution of quadratic programming to approach the optimal solution of the search direction by solving a series of quadratic programming problems with equality constraints step by step. After that, the step size was obtained by line search. The next iteration point was then obtained. Finally, the global convergence of the algorithm was proved theoretically. A preliminary numerical experiment was carried out to compare the algorithm with the smoothed sequential projection shrinkage algorithm. The results show that the algorithm has certain advantages not only in iteration times but also in computing time.
Key words: smooth approximation     polyhedral constraint     composite function     sequential quadratic programing     active set    
1 问题的提出

本文考虑如下复合函数优化问题

$ \left\{ \begin{array}{*{20}{l}} \mathop {\min }\limits_x \quad \phi \left( {\boldsymbol{x}} \right): = f\left( {\boldsymbol{x}} \right) + g\left( {\boldsymbol{x}} \right) \\ {\text{s}}{\text{.t}}{\text{.}}\quad \;\;{\boldsymbol{x}} \in \varOmega \\ \end{array} \right.$ (1)

式中: $ f:{\mathbb{R}^n} \to \mathbb{R} $ 是连续可微函数; $ g:{\mathbb{R}^n} \to \mathbb{R} $ 是凸的非光滑函数;可行集 $\varOmega = \{ {\boldsymbol{x}} \in {\mathbb{R}^n}:{{\boldsymbol{a}}_i}^{\text{T}}{\boldsymbol{x}} = {{\boldsymbol{b}}_i},i \in$ $ {\mathcal{M}_E}, $ ${{\boldsymbol{a}}_i}^{\rm{T}}x \geqslant {{\boldsymbol{b}}_i},{\text{ }}i \in {\mathcal{M}_I}\}$ ${\mathcal{M}_E} = \left\{ {1,\cdots,m} \right\}$ ${\mathcal{M}_I} = \left\{ {m + 1,\cdots,p} \right\}$ ${{\boldsymbol{a}}_i} \in {\mathbb{R}^n}$ ${{\boldsymbol{b}}_i} \in \mathbb{R}$ $ i = 1,\cdots,p $

从现实生活中的问题抽象出来的数学模型,可以分为光滑和非光滑两大类。对于光滑问题,由于其梯度及Hessian矩阵易于求出,可以使用经典的优化算法求解。而对非光滑问题,函数是不可微的,因此无法直接使用经典的优化算法求解。此外,非光滑复合函数的极小化问题在信息理论[1]、无线通信[2-4]、图像恢复[5-7]、信号处理[8-9]、变量选择[4, 10]等领域都有广泛的应用。

当可行集 $ \varOmega = {\mathbb{R}^n} $ 时,问题(1)就是无约束优化问题。目前,已经有很多学者提出了各种非常有效的算法[11-18]解决无约束复合函数问题。比较成熟的算法有迭代收缩阈值算法[12-13]、临近牛顿法[14]、分裂算法[15-16]、邻近点算法[17-18]等。

实际生活中的很多问题都要求满足一些特定的条件,所以研究带约束的复合函数问题十分有必要。例如,压缩感知理论中的稀疏信号重构问题就是经典的复合函数优化问题。稀疏重构主要解决这两类问题:在对非负信号处理时,有非负约束的限制;进行图像恢复时,灰度强度取值在0~255之间。在信号恢复、图像恢复模型里加入约束,可以提高恢复的准确率。同样地,对于实际问题,考虑约束限制,能更加精确地求出符合实际条件的解,从而节约成本。

当可行集 $\varOmega = \left\{ {{\boldsymbol{l}} \leqslant {\boldsymbol{x}} \leqslant {\boldsymbol{u}}} \right\}$ 时,问题(1)就是界约束的复合函数问题。Bian等[6]提出光滑二次正则化方法,用SQR表示解决图像恢复中出现的界约束非Lipschitz连续优化问题。有效集方法是求解约束优化问题的常用方法,采用有效集策略,可以加快算法的收敛速度。张涛[19]提出了一种求解界约束优化问题的新积极集信赖域算法。

生产生活中的大多数问题,往往具有更加复杂的约束,因此研究复杂约束的复合函数优化问题更具有现实意义。当可行集 $ \varOmega $ 为多面体约束时,问题(1)就是多面体约束复合函数优化问题。简单约束复合函数优化问题的理论已经发展得十分成熟,但是关于多面体约束的复合函数优化问题,研究工作仍处于基础阶段,还有很多问题尚未解决。Liu等[4]提出了光滑序列二次规划框架(SSQP)求解一类多面体上复合函数 $ {L_q} $ 极小化问题。实际上,有效集方法已成功地应用于大规模的线性约束光滑优化问题的求解。William等[20]提出了求解多面体约束优化问题的新的有效集算法 (PASA)。Zhang等[1]提出了求解线性约束非凸非Lipschitz优化的光滑有效集算法 (SASA),且用于求解简单约束非光滑复合函数优化问题效果良好,即

$ \left\{ \begin{array}{*{20}{l}} \mathop {\min }\limits_x \quad \phi \left( {\boldsymbol{x}} \right): = f\left( {\boldsymbol{x}} \right) + \left\| {\boldsymbol{x}} \right\|_q^q \\ {\text{s}}{\text{.t}}{\text{.}}\quad \;\;{\boldsymbol{x}} \in \varOmega \\ \end{array} \right.$ (2)

也就是本文实际求解的模型:

$ \left\{ \begin{array}{*{20}{l}} \mathop {\min }\limits_x \quad \phi \left( {\boldsymbol{x}} \right): = f\left( {\boldsymbol{x}} \right) + \left\| {{\boldsymbol{b}} - {\boldsymbol{Ax}}} \right\|_q^q \\ {\bf{s}}{\bf{.t.}}\quad \;\;{\boldsymbol{x}} \in \varOmega \\ \end{array} \right.$ (3)

它是Liu等[4]在求解一类多面体上复合函数 $ {L_q} $ 极小化问题中求解的模型。

这两个问题之间的一个明显区别是复合项 $ \left\| {\, \cdot \,} \right\|_q^q $ 不同。具体来说:问题(2)试图找到一个解 ${{\boldsymbol{x}}^*}$ ,使 ${\boldsymbol{x}}$ 的非零项的数量尽可能得少;而问题(3)尝试找到一个解 ${{\boldsymbol{x}}^*}$ ,使向量 ${\boldsymbol{b}} - {\boldsymbol{Ax}}$ 正分量的数量尽可能少。问题(2)本质上是一般约束稀疏优化,问题(3)本质上是不等式约束稀疏优化。

以上算法在求解大多数问题时表现良好,但它们不能很好地解决一般多面体约束复合函数问题,例如,SQR算法[6]不能求解问题中 $f({\boldsymbol{x}})$ 的复合项 $ {L_q} $ 和一般多面体约束问题。SSQP 算法[4]的目标函数不具一般性,且由于投影收缩算法是不精确求解子问题,这会使得算法收敛速度慢。PASA算法[20]需要求解投影梯度信息,SASA算法[1]对线性约束优化器有特殊条件限制,这两种算法都是两阶段算法,程序实现难度较大。

序列二次规划(sequential quadratic programming, SQP)是求解非线性约束优化问题的主要方法之一,它在求解中等规模和小规模的非线性问题[21- 22]中表现不俗。有效集算法 (active set algorithm, ASA) 在求解不等式约束优化问题时十分有效,通过有效集策略,每次只需求解仅含等式约束的优化问题,从而可以降低维度,达到提高算法收敛速度的目的。

本文提出了将ASA和SQP结合的新算法,该算法主要借鉴文献[4]中提到的求解多面体上最小化复合函数的光滑近似技巧及SSQP框架,将非光滑目标函数转化成光滑函数,并利用有效集策略求解二次规划子问题。

2 光滑近似

由于问题(1)中涉及非光滑项 $\left\| {\boldsymbol{x}} \right\|_q^q$ ,导致无法直接使用经典的优化算法来求解。因此先给出近似光滑函数这一重要概念。

定义1  [23]  函数 $ {\tilde g_\mu }:{\mathbb{R}^n} \times {\mathbb{R}_ + } \to \mathbb{R} $ 称为非光滑函数 $ g:{\mathbb{R}^n} \to \mathbb{R} $ 的近似光滑函数,若 $ \forall \mu \gt 0 $ $ {\tilde g_\mu } $ 是连续可微的,且对 $\forall \mu \in {\mathbb{R}^n}$ ,有 $\mathop {\lim }\limits_{{\boldsymbol{z}} \to {\boldsymbol{x}},\;\mu \downarrow 0} \tilde g\left( {{\boldsymbol{z}},\mu } \right) = g\left( {\boldsymbol{x}} \right)$

非光滑优化问题的光滑近似已经得到了广泛的应用[23- 24]。本文对绝对值函数

$ \theta \left( t \right) = \left| t \right| $

构造光滑函数[25]

$ {\theta _\mu }\left( {t, \mu } \right) = \left\{ \begin{array}{*{20}{l}} t,& t > \dfrac{\mu }{2} \\ \dfrac{{{t^2}}}{\mu } + \dfrac{\mu }{4}, &- \dfrac{\mu }{2} \leqslant t \leqslant \dfrac{\mu }{2} \\ - t, & t < - \dfrac{\mu }{2} \end{array} \right. $ (4)

基于式(4),可以构造非光滑函数 $ \phi $ 的近似光滑函数 $ \tilde \phi $ 。定义

$ \left\{ \begin{array}{*{20}{l}} \mathop {\min }\limits_x \;{\text{ }}\tilde \phi ({\boldsymbol{x}},\mu ): = f\left( {\boldsymbol{x}} \right) + \tilde g\left( {{\boldsymbol{x}},\mu } \right) \\ {\text{s}}{\text{.t}}{\text{.}}\;\;\;\,{\boldsymbol{x}} \in \varOmega \\ \end{array} \right. $ (5)

则式(5)是问题(1)的光滑近似函数。

不难得到 $ {\theta _\mu }\left( {t, \mu } \right) $ 具有的特殊性质。显然地,对任意固定的 $ \mu > 0 $ ,有

$ {\theta _\mu }(t,\mu ) = \theta (t),\;\forall t \geqslant \mu $

$ \theta (t,\mu ) \geqslant \dfrac{\mu }{4},\;\,\forall t $

此外, $ {\theta ^q}(t,\mu ) $ 是连续可微的。除在 $ t = \dfrac{\mu }{2} $ $ t = - \dfrac{\mu }{2} $ 外, $ {\theta ^q}(t,\mu ) $ 在其他处都是二次连续可微的。 $ {\theta ^q}(t,\mu ) $ 关于 $ t $ 的一阶导数和二阶导数分别为

$ {\left[ {{\theta ^q}(t,\mu )} \right]^\prime } = \left\{ {\begin{array}{*{20}{l}} {q{t^{q - 1}},{\text{ }}\quad \quad \quad \quad \;t > \dfrac{\mu }{2}} \\ {2q{\theta ^{q - 1}}\left( {t,\mu } \right)\dfrac{t}{\mu },{\text{ }} - \dfrac{\mu }{2} \leqslant t \leqslant \dfrac{\mu }{2}} \\ { - q{t^{q - 1}},\quad \quad \quad \quad \;t < - \dfrac{\mu }{2}} \end{array}} \right.\quad $
$ {\left[ {{\theta ^q}(t,\mu )} \right]^{\prime \prime }} = \left\{ {\begin{array}{*{20}{l}} {q\left( {q - 1} \right){t^{q - 2}},\qquad \qquad t > \dfrac{\mu }{2}}\\ \begin{array}{l} q{\theta ^{q - 1}}\left( {t,\mu } \right)\dfrac{2}{\mu } + q\left( {q - 1} \right)\cdot \\\quad {\theta ^{q - 2}}(t,\mu )\dfrac{{2{t^2}}}{{{\mu ^2}}},\qquad - \dfrac{\mu }{2} \leqslant t \leqslant \dfrac{\mu }{2} \end{array}\\ { - q\left( {q - 1} \right){t^{q - 2}}, \qquad \qquad t < - \dfrac{\mu }{2}} \end{array}} \right. $

接下来,给出一个非常重要的引理。

引理1 对 $ \forall q \in \left( {0,{\text{ }}1} \right) $ $ \mu \in (0,\,\, + \infty ) $

a. $0 \leqslant {\theta ^q}\left( {t,\mu } \right) - {\theta ^q}\left( t \right) \leqslant {\left( {\dfrac{\mu }{4}} \right)^q},\forall t \leqslant \dfrac{\mu }{2}$ ,且有 ${\left[ {{\theta ^q}(t,\mu )} \right]^{\prime \prime }} \leqslant 16q{\mu ^{q - 2}}$ $\forall t \notin \left\{ - \dfrac{\mu }{2},\dfrac{\mu }{2}\right\}$

b. 定义

$ \kappa (t,\mu )\left\{\begin{array}{l}16q{\mu }^{q-2},\text{ }0\leqslant t\leqslant \mu \\ 0,\text{ }其他 \end{array}\right. $

$ \forall t $ $ \hat t $ ,如果 $ \hat t > \mu $ 使得 $t - \hat t \geqslant - \dfrac{{\hat t}}{2}$ ,或者 $0 \leqslant t \leqslant \mu$ 时, $ t \in ( - \infty ,\infty ) $ ,或者 $ \hat t < 0 $ 时, $ t - \hat t \leqslant 0 $ ,有

$ {\theta ^q}\left( {t,\mu } \right) \leqslant {\theta ^q}\left( {\hat t,\mu } \right) + {\left[ {{\theta ^q}(\hat t,\mu )} \right]^\prime }(t - \hat t) + \frac{{\kappa (\hat t,\mu )}}{2}{(t - \hat t)^2} $

成立。

引理1的详细证明参见文献[4]中的引理 4.1。

3 算 法

本文受文献[4]的启发,采用该文献中的SSQP框架,提出了序列有效集算法 (SQP-ASA),在每一迭代步中利用有效集策略求解一个凸二次规划(QP)子问题。此QP子问题的目标函数是问题(1)光滑近似后式(5)的目标函数 $ \tilde \phi (x,\mu ) $ 的局部上界。光滑参数的更新依赖于问题(1)的残差。

对固定的光滑参数 $ \mu \gt 0 $ ,定义非光滑函数 $ g(\; \cdot \;) $ 的光滑近似函数 $ \tilde g( \cdot \,,\mu ) $ ${{\boldsymbol{x}}_k}$ 处的二次近似为

$ \begin{aligned} {Q_1}\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu } \right) =& \tilde g\left( {{{\boldsymbol{x}}_k},\mu } \right) + \nabla \tilde g{\left( {{{\boldsymbol{x}}_k},\mu } \right)^{\text{T}}}\left( {{\boldsymbol{x}} - {{\boldsymbol{x}}_k}} \right) + \\&\frac{1}{2}{({\boldsymbol{x}} - {{\boldsymbol{x}}_k})^{\text{T}}}\tilde B({{\boldsymbol{x}}_k},\mu )({\boldsymbol{x}} - {{\boldsymbol{x}}_k}) \end{aligned} $

式中, $\tilde g({{\boldsymbol{x}}_k},\mu ) = \displaystyle \sum\limits_{i = 1}^n {{\theta ^q}} ({{\boldsymbol{x}}_k},\mu )$ $\tilde B({\boldsymbol{x}},\mu ) = {\text{Diag}}(\kappa ({{\boldsymbol{x}}_1},\mu ),\cdots, \kappa ({{\boldsymbol{x}}_n},\mu ))$

类似地,定义光滑函数 $ f\left( {\boldsymbol{x}} \right) $ $ {{\boldsymbol{x}}_k} $ 处的二次近似为

$ {Q_2}({\boldsymbol{x}},{{\boldsymbol{x}}_k}) = f({{\boldsymbol{x}}_k}) + \nabla f{({{\boldsymbol{x}}_k})^{\rm{T}}}({\boldsymbol{x}} - {{\boldsymbol{x}}_k}) + \frac{1}{2}{L_k}{\left\| {{\boldsymbol{x}} - {{\boldsymbol{x}}_k}} \right\|^2} $

假设f(x)的梯度是Lipschitz连续的,即

$ {\left\| {\nabla f({\boldsymbol{x}}) - \nabla f({\boldsymbol{y}})} \right\|_2} \leqslant L{\left\| {{\boldsymbol{x}} - {\boldsymbol{y}}} \right\|_2},\;\forall {\boldsymbol{x}},{\boldsymbol{y}} \in \varOmega $

定义

$ Q({\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu ) = {Q_1}({\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu ) + {Q_2}({\boldsymbol{x}},{{\boldsymbol{x}}_k}) $ (6)

由于 $ {Q_1}({\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu ) $ $ {Q_2}({\boldsymbol{x}},{{\boldsymbol{x}}_k}) $ 分别是 $ \tilde g( \cdot \,,\mu ) $ $ f\left( {\boldsymbol{x}} \right) $ $ {{\boldsymbol{x}}_k} $ 处的二次近似,因此,可以把 $ {Q_1}({\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu ) $ 看作是 $ \tilde \phi ({\boldsymbol{x}},\mu ) $ $ {{\boldsymbol{x}}_k} $ 处的二次近似。

下面的引理给出了凸二次规划(6)的目标函数 $ Q({\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu ) $ 是式(5)的目标函数 $ \tilde \phi (x,\mu ) $ 之间的大小关系,即 $ Q({\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu ) $ $ \tilde \phi ({\boldsymbol{x}},\mu ) $ $ {{\boldsymbol{x}}_k} $ 处的一个上界,这保证了问题(6)的最小值是问题(5)最小值的上界。

引理2 对任意的 ${{\boldsymbol{x}}_k}$ ${\boldsymbol{x}}$ ,使得

$ \begin{gathered} {\left( {{{\boldsymbol{x}}_k} - {\boldsymbol{x}}} \right)_e} \leqslant 0,\,\;e \in \mathcal{I}_{{{\boldsymbol{x}}_k}}^\mu \\ {\left( {{{\boldsymbol{x}}_k} - {\boldsymbol{x}}} \right)_e} \geqslant - \dfrac{{{{({{\boldsymbol{x}}_k})}_e}}}{2},\;e \in \mathcal{J}_{{{\boldsymbol{x}}_k}}^\mu \\ \end{gathered} $

其中

$ \begin{gathered} \mathcal{I}_{{{\boldsymbol{x}}_k}}^\mu = \left\{ {e|{{({{\boldsymbol{x}}_k})}_e} < 0} \right\} \\ \mathcal{J}_{{{\boldsymbol{x}}_k}}^\mu = \left\{ {e|{{({{\boldsymbol{x}}_k})}_e} > \mu } \right\} \\ \end{gathered} $

如果 $ {L_k} \gt 0 $ $ f({\boldsymbol{x}}) \leqslant {Q_2}({\boldsymbol{x}},{{\boldsymbol{x}}_k}) $ ,那么 $\tilde \phi ({\boldsymbol{x}},\mu ) \leqslant Q({\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu )$

该引理的证明可参考文献[4]中的引理 4.2。

基于引理2,提出序列有效集二次规划算法。

本文研究的问题(2)光滑后可以表示为

$ \left\{\begin{gathered} \mathop {\min }\limits_x \quad {\mkern 1mu} \tilde \phi ({\boldsymbol{x}},\mu ): = f({\boldsymbol{x}}) + \tilde g({\boldsymbol{x}},\mu ) \\ {\text{s}}{\text{.t}}{\text{.}}\quad \;\; {\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{x}} = {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_E} \\ {\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{x}} \geqslant {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_I} \end{gathered}\right. $ (7)

式中:函数 $ \tilde \phi $ 是连续可微的;指标集 ${\mathcal{M}_E} = \left\{ {1,\cdots,m} \right\}, {\mathcal{M}_I} = \left\{ {m + 1,\cdots,p} \right\}$ ${{\boldsymbol{a}}_i} \in {\mathbb{R}^n}$ $ {{\boldsymbol{b}}_i} \in \mathbb{R} $ $i = 1,\cdots,p$

那么,问题(7)在 $ {{\boldsymbol{x}}_k} $ 处的QP子问题为

$ \left\{\begin{array}{l} \min \quad {\mkern 1mu} Q\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu } \right) \\ {\text{s}}{\text{.t}}{\text{.}}\quad \;\; {{{\left( {{{\boldsymbol{x}}_k} - {\boldsymbol{x}}} \right)}_e} \leqslant 0,\;e \in \mathcal{I}_{{{\boldsymbol{x}}_k}}^\mu } \\ \quad\quad\;\;\;{{{\left( {{{\boldsymbol{x}}_k} - {\boldsymbol{x}}} \right)}_e} \geqslant - \dfrac{{{{({{\boldsymbol{x}}_k})}_e}}}{2},\;e \in \mathcal{J}_{{{\boldsymbol{x}}_k}}^\mu } \\ \quad\quad\;\;\; {{\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{x}} = {b_i},\;i \in {\mathcal{M}_E}} \\ \quad\quad\;\;\; {{\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{x}} \geqslant {b_i},\;i \in {\mathcal{M}_I}} \end{array}\right. $ (8)

$ {{\boldsymbol{d}}_k} = {\boldsymbol{x}} - {{\boldsymbol{x}}_k} $ $ {{\boldsymbol{H}}_k} = {L_k}{{\boldsymbol{I}}_n} + \tilde {\boldsymbol{B}}({{\boldsymbol{x}}_k},\mu ) $ ,则 $ Q\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_k},\mu } \right) = \nabla \tilde \phi {({{\boldsymbol{x}}_k},\mu )^{\text{T}}}{{\boldsymbol{d}}_k} + \dfrac{1}{2}{\boldsymbol{d}}_k^{\text{T}}{{\boldsymbol{H}}_k}{{\boldsymbol{d}}_k} $

求解问题(7)的算法思想为:在每一迭代步中,首先使用有效集方法求出其二次规划子问题(8)的最优解 $ {{\boldsymbol{d}}^*} $ ,并将它作为问题(7)下一步的搜索方向 $ {{\boldsymbol{d}}_k} $ ,再利用价值函数通过Armijo线搜索求出步长 $ {\alpha _k} $ ,进而得到问题(7)的下一步迭代点 $ {{\boldsymbol{x}}_{k + 1}} $ ,重复这一过程,直到求得问题(7)的最优解 $ {{\boldsymbol{x}}^*} $ 。具体求解过程如算法1所示。

算法1  序列有效集算法第 $ k $ 步迭代

步骤1 选择初始可行点 $ {{\boldsymbol{x}}_k} $ 和参数 $\epsilon \in (0,1]$ $ 0 \lt \sigma \lt 1 $ $ \eta \gt 1 $ $ {L^{\max }} \geqslant {L_0} \geqslant {L^{\min }} \gt 0 $ 。令 $ z = 0 $ (外部迭代索引), $ k = 0 $ (内部迭代索引), ${\mu }_{0}=\dfrac{\epsilon }{{\sigma }^{\lfloor {\mathrm{log}}_{\sigma }\epsilon\rfloor }}\in (\sigma ,1]$

步骤2 令 $ \mu = {\mu _z} $ ,如果 $ \left\| {{{\boldsymbol{d}}_k}} \right\| \leqslant \mu $ ,执行步骤5;否则执行步骤3。

步骤3 参照文献[26]中的有效集方法求解QP子问题,得到 $ {{\boldsymbol{d}}_k} $ 。通过价值函数 $ \psi ({\boldsymbol{x}},\mu ,\vartheta ) $ ,根据Armijo线搜索求出步长 $ {\beta _k} $ ,进而得到 $ {{\boldsymbol{x}}_{k + 1}} = {{\boldsymbol{x}}_k} + {\beta _k}{{\boldsymbol{d}}_k} $ $ {r_k} $ ,其中

${r_k} = \dfrac{{f({{\boldsymbol{x}}_{k + 1}}) - f({{\boldsymbol{x}}_k}) - \nabla f{{({{\boldsymbol{x}}_k})}^{\text{T}}}({{\boldsymbol{x}}_{k + 1}} - {{\boldsymbol{x}}_k})}}{{\dfrac{1}{2}{L_k}{{\left\| {{{\boldsymbol{x}}_{k + 1}} - {{\boldsymbol{x}}_k}} \right\|}^2}}}$

步骤4 令 ${{\boldsymbol{s}}_{k + 1}} = {{\boldsymbol{x}}_{k + 1}} - {{\boldsymbol{x}}_k}$ ${{\boldsymbol{y}}_{k + 1}} = \nabla f({{\boldsymbol{x}}_{k + 1}}) - \nabla f({{\boldsymbol{x}}_k})$ ,如果 $ {r_k} \leqslant 1 $ $ {L_{k + 1}} = $ $\max \left\{ {\rm{min}}\left\{ {L^{\max }},\dfrac{{{\boldsymbol{s}}_{k + 1}^{\text{T}}{{\boldsymbol{y}}_{k + 1}}}}{{{{\left\| {{{\boldsymbol{s}}_{k + 1}}} \right\|}^2}}}, \; {L^{\min }}\right\} \right\}$ $ k = k + 1 $ ,执行步骤2,否则 $ {L_k} = \eta {L_k} $ ,执行步骤3。

步骤5 如果 $\mu \leqslant \epsilon$ ,则终止该算法;否则,转到步骤6。

步骤6 令 $ {\mu _{z + 1}} = \sigma {\mu _z} $ $ z = z + 1 $ ${{\boldsymbol{x}}_0} = {{\boldsymbol{x}}_k}$ $k = 0$ ,然后转到步骤2。

${\boldsymbol{x}} - {{\boldsymbol{x}}_k} = {\boldsymbol{d}}$ ,则可将QP子问题(8)转化与 $ {\boldsymbol{d}} $ 有关的子问题求解。通过积极集策略,会有越来越多的 $ {{\boldsymbol{x}}_k} $ 逐渐满足 ${\boldsymbol{a}}_i^{\rm{T}}{{\boldsymbol{x}}_k} = {{\boldsymbol{b}}_i},\;i \in {{{\mathcal{S}}}_k}$ 。这样问题就转化为了等式约束优化问题,再利用拉格朗日乘子算法对其进行求解。因此,本文的算法可以在大大减少问题维数的同时提高计算的效率。

接下来,利用文献[26]中的有效集方法求解QP子问题,如算法2所示。

算法2  有效集方法第 $ k $ 步迭代

步骤1 选择初始可行点 $\forall {{\boldsymbol{x}}_0},\;{{\boldsymbol{d}}_0} \in \varOmega$ ,对称正定矩阵 $ {{\boldsymbol{H}}_0} $ ,令 $\forall x \in \varOmega$

步骤2 确定相应的有效集 ${{\rm{{\cal S}}}_k}$ ,求解子问题

$ \begin{gathered} \min \quad {\mkern 1mu} {q_k}\left( {\boldsymbol{t}} \right) = \frac{1}{2}{{\boldsymbol{t}}^{\text{T}}}{{\boldsymbol{H}}_k}{\boldsymbol{t}} + {\boldsymbol{C}}_k^{\text{T}}{\boldsymbol{t}} \\ {\text{s}}{\text{.t}}{\text{.}}\quad \;{\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{t}} = 0,\;i \in {\mathcal{S}_k} \\ \end{gathered} $

得极小点 $ {{\boldsymbol{d}}_k} $ 和拉格朗日乘子 $\lambda_k$ ,其中 ${{\boldsymbol{C}}_k} = {{\boldsymbol{H}}_0}{{\boldsymbol{d}}_k} + {{\boldsymbol{g}}_0}$ ,如果 $ {{\boldsymbol{t}}_k} = 0 $ ,执行步骤3;否则执行步骤4。

步骤3 令 $ {({\lambda _k})_l} = \mathop {\min }\limits_{i \in {{\cal{S}}_k}} \{ {({\lambda _k})_l}\} $ ,如果 ${({\lambda _k})_l} \geqslant 0$ ,则 $d_k$ 是全局极小点,停算。否则,若 ${({\lambda _k})_l} <0$ ,则令 $ {\mathcal{S}_k}: = {\mathcal{S}_k}\backslash \left\{ l \right\} $ ,转步骤2。

步骤4 令 ${\alpha _k} = \min \{ 1,{{\bar \alpha }_k}\} $ ,其中

${\bar \alpha _k} = \mathop {\min }\limits_{i \in {\mathcal{S}_k}} \left\{ {\dfrac{{{c_i} - {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{d}}_k}}}{{{\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{t}}_k}}}\Bigg| {{\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{t}}_k} < 0} } \right\}$

${{\boldsymbol{d}}_{k + 1}}: = {{\boldsymbol{d}}_k} + {\alpha _k}{{\boldsymbol{t}}_k}$ ,其中 ${c_i} = {{\boldsymbol{b}}_i} - {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_0}$

步骤5 若 ${\alpha _k} = 1$ ,则令 $ {\mathcal{S}_{k + 1}}: = {\mathcal{S}_k} $ ;否则,若 ${\alpha _k} < 1$ ,则令 $ {\mathcal{S}_{k + 1}}: = {\mathcal{S}_k} \cup \left\{ {{j_k}} \right\} $ ,其中 $ \left\{ {{j_k}} \right\} $ 满足 ${\bar \alpha _k} = \dfrac{{{c_{{j_k}}} - {\boldsymbol{a}}_{{j_k}}^{\text{T}}{{\boldsymbol{d}}_k}}}{{{\boldsymbol{a}}_{{j_k}}^{\text{T}}{{\boldsymbol{t}}_k}}}$

步骤6  $k=k+1$ ,转步骤 2。

此处的价值函数 $\psi \left( {{\boldsymbol{x}},\mu ,\vartheta } \right) = \nabla \tilde \phi \left( {{\boldsymbol{x}},\mu } \right)$ 。这是因为由有效集方法得到的方向 ${{\boldsymbol{d}}_k}$ 一定是可行方向。

4 收敛性分析

为证明算法1的全局收敛,现给出如下基本假设:

假设1 对 $\forall {\boldsymbol{x}} \in \varOmega$ ,线性独立约束限定 (LICQ) 成立,即梯度 ${{\boldsymbol{a}}_i}(\,\;i \in {\mathcal{S}_k} )$ 线性无关。

$ {{\boldsymbol{H}}_k} $ 的构造可知, $ {{\boldsymbol{H}}_k} $ 一定是对称正定的。本文提出的SQP-ASA算法通过求解如下形式的二次规划问题

$\left\{ \begin{gathered} \min \quad \nabla \tilde \phi {({{\boldsymbol{x}}_k},\mu )^{\text{T}}}{\boldsymbol{d}} + \frac{1}{2}{{\boldsymbol{d}}^{\rm{T}}}{{\boldsymbol{H}}_k}{\boldsymbol{d}} \\ {\text{s}}{\text{.t}}{\text{.}}\quad \; {\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{d}} + {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} = {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_E} \\ \quad\quad {\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{d}} + {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} \geqslant {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_I} \end{gathered}\right. $ (9)

得到解 $ {{\boldsymbol{d}}^*} $ ,作为原问题变量 $ {\boldsymbol{x}} $ 在第 $ k $ 次迭代过程中的搜索方向 $ {{\boldsymbol{d}}_k} $ 。在求解过程中利用有效集策略,通过逐次求解仅含带等式约束的相关二次规划问题(9)来逼近其最优解。

下面证明SQP-ASA算法的全局收敛性,即证该算法产生的点列 $ \left\{ {{{\boldsymbol{x}}_k}} \right\} $ 的任何聚点 $ {{\boldsymbol{x}}^*} $ 都是问题(7)的KT 点。

引理3  由算法2中的有效集方法求得的QP子问题方向 $ {{\boldsymbol{d}}_k} $ 一定是原问题(7)的可行下降方向。

证明 先证明 $ {{\boldsymbol{d}}_k} $ 是可行方向。依据有效集的求解过程可知: $ {{\boldsymbol{d}}_k} $ 是QP子问题(9)的可行方向,从而有 $ {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{d}}_k} \geqslant {c_i} $ $ {c_i} = {{\boldsymbol{b}}_i} - {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} $ ${\boldsymbol{ a}}_i^{\text{T}}{{\boldsymbol{d}}_k} = 0 $ $ i \in {\mathcal{S}_k} $ 。又因初始点 $ {{\boldsymbol{x}}_k} $ 是可行点,则有 $ {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} \geqslant {{\boldsymbol{b}}_i} $ ${\boldsymbol{ a}}_i^{\text{T}}{{\boldsymbol{x}}_k} = {{\boldsymbol{b}}_i}$ $ i \in {\mathcal{S}_k} $ 。设由序列二次规划求出的下一步迭代点是 $ {{\boldsymbol{x}}_{k + 1}} $ ,则

$ \begin{aligned} {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_{k + 1}} =& {\boldsymbol{a}}_i^{\text{T}}({{\boldsymbol{x}}_k} + {\beta _k}{{\boldsymbol{d}}_k}) = {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} + {\beta _k}{\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{d}}_k} =\\& {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} + 0 = {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_E} \end{aligned} $
$ \begin{gathered} {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_{k + 1}} = {\boldsymbol{a}}_i^{\text{T}}({{\boldsymbol{x}}_k} + {\beta _k}{{\boldsymbol{d}}_k}) = {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} + {\beta _k}{{\boldsymbol{a}}_i}{{\boldsymbol{d}}_k} \geqslant \\ {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} + {\beta _k}({{\boldsymbol{b}}_i} - {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k}) \geqslant (1 - {\beta _k}){\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} + \\ {\beta _k}{{\boldsymbol{b}}_i} \geqslant (1 - {\beta _k}){{\boldsymbol{b}}_i} + {\beta _k}{{\boldsymbol{b}}_i} = {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_I} \\ \end{gathered} $

因此, $ {{\boldsymbol{d}}_k} $ 一定是原问题(7)的可行方向。

其次,证明方向 $ {{\boldsymbol{d}}_k} $ 是下降的。对 $ \forall \mu \in (0,1] $ ,利用QP子问题(9)的KT条件,有

$ \left\{ \begin{gathered} \nabla \tilde \phi ({{\boldsymbol{x}}_k},\mu ) + {{\boldsymbol{H}}_k}{{\boldsymbol{d}}_k} = \mathop \sum \limits_{i \in {\mathcal{M}_I}} {\lambda _i}{{\boldsymbol{a}}_i} + \mathop \sum \limits_{i \in {\mathcal{M}_E}} {\lambda _i}{{\boldsymbol{a}}_i} \\ {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{d}}_k} = 0,\;i \in {\mathcal{M}_E} \\ {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{d}}_k} \geqslant {c_i},\;{c_i} = {{\boldsymbol{b}}_i} - {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k},\;i \in {\mathcal{M}_I} \\ {\lambda _i} \geqslant 0,\;i \in {\mathcal{M}_I};{\lambda _i} = 0,\;i \in \mathcal{M}\backslash {\mathcal{S}_k} \\ \end{gathered} \right. $ (10)

由式(10)的第一个式子可得

$ \begin{aligned} \nabla \tilde \phi ({{\boldsymbol{x}}_k},\mu ){{\boldsymbol{d}}_k} =& - {\boldsymbol{d}}_k^{\text{T}}{{\boldsymbol{H}}_k}{{\boldsymbol{d}}_k} - \mathop \sum \limits_{i \in {\mathcal{M}_I}} {\lambda _i}({\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} - {{\boldsymbol{b}}_i}) - \\ &\mathop \sum \limits_{i \in {\mathcal{M}_E}} {\lambda _i}({\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} - {{\boldsymbol{b}}_i}) = - {\boldsymbol{d}}_k^{\rm{T}}{{\boldsymbol{H}}_k}{{\boldsymbol{d}}_k} -\\ & \mathop \sum \limits_{i \in {\mathcal{S}_k}} {\lambda _i}({\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} - {{\boldsymbol{b}}_i}) - \mathop \sum \limits_{i \in \mathcal{M}\backslash {\mathcal{S}_k}} {\lambda _i}({\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} - {{\boldsymbol{b}}_i}) \\ \end{aligned} $

$ i \in {\mathcal{S}_k} $ 时, ${\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} - {{\boldsymbol{b}}_i} = 0$ ;当 $ i \in \mathcal{M}\backslash {\mathcal{S}_k} $ 时, $ {\lambda _i} = 0 $ ,从而 $\displaystyle \mathop \sum \limits_{i \in \mathcal{M}} {\lambda _i}({\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}_k} - {{\boldsymbol{b}}_i}) = 0$ ,即 $\nabla \tilde \phi ({{\boldsymbol{x}}_k},\mu ){{\boldsymbol{d}}_k} = - {\boldsymbol{d}}_k^{\text{T}}{{\boldsymbol{H}}_k}{{\boldsymbol{d}}_k}$ 。又因为 ${{\boldsymbol{H}}_k}$ 是对称正定矩阵,所以对任何非负向量 ${{\boldsymbol{d}}_k}$ 都有 ${\boldsymbol{d}}_k^{\text{T}}{{\boldsymbol{H}}_k}{{\boldsymbol{d}}_k} > 0$ ,故 $\nabla \tilde \phi ({{\boldsymbol{x}}_k},\mu ){{\boldsymbol{d}}_k} < 0$

引理4 对于上述二次规划子问题(9),设其KT点对为 $({\boldsymbol{d}}({\boldsymbol{x}}),\tilde \lambda ({\boldsymbol{x}}))$ 。如果 ${\boldsymbol{d}}({\boldsymbol{x}}) = 0$ ,则 ${\boldsymbol{x}}$ 为问题(7)的 KT点。

证明 不妨设

$ \lambda ({\boldsymbol{x}}) = \left( {{\lambda _i}({\boldsymbol{x}}),i \in \mathcal{M}} \right) \text{,} {\lambda _i}({\boldsymbol{x}}) = \left\{ \begin{gathered} {{\tilde \lambda }_i}({\boldsymbol{x}}),\;i \in \mathcal{S}({\boldsymbol{x}}) \\ 0,\;\;\quad \;i \in \mathcal{M}\backslash {\mathcal{S}}({\boldsymbol{x}}) \\ \end{gathered} \right. $

如果 $ {\boldsymbol{d}}({\boldsymbol{x}}) = 0 $ ,依据有效集方法求解方向 $ {\boldsymbol{d}}({\boldsymbol{x}}) $ 的过程,结合问题(7)的KT条件可知, $ {\tilde \lambda _i}({\boldsymbol{x}}) \geqslant 0 $ ,否则, $ {\boldsymbol{d}}({\boldsymbol{x}}) $ 不是二次规划问题(9)的最优解,与KT点对这个条件相矛盾。因此

$ \left\{ \begin{gathered} \nabla \tilde \phi ({\boldsymbol{x}},\mu ) = \mathop \sum \limits_{i \in {\mathcal{M}_I}} {\lambda _i}{{\boldsymbol{a}}_i} + \mathop \sum \limits_{i \in {\mathcal{M}_E}} {\lambda _i}{{\boldsymbol{a}}_i} \\ {\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{x}} = {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_E} \\ {\boldsymbol{a}}_i^{\text{T}}{\boldsymbol{x}} \geqslant {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_I} \\ {\lambda _i} \geqslant 0,\;i \in {\mathcal{M}_I};{\lambda _i} = 0,\;i \in \mathcal{M}\backslash \mathcal{S}(x) \\ \end{gathered} \right. $

这说明 ${\boldsymbol{x}}$ 是问题(7)的一个KT点。

定理1 由序列有效集算法生成的序列 $\left\{ {{{\boldsymbol{x}}_k}} \right\}$ 的任何聚点 ${{\boldsymbol{x}}^*}$ 都是原问题的KT点。

证明 对 $ \forall \mu \in (0,1] $ ,当该算法是有限步终止时,由引理4可知,算法产生的点列 $ \left\{ {{{\boldsymbol{x}}_k}} \right\} $ 的任何聚点 $ {{\boldsymbol{x}}^*} $ 都是问题(7)的KT点。

下设该算法产生无穷点列 $ \left\{ {{{\boldsymbol{x}}_k}} \right\} $ $ {{\boldsymbol{x}}^*} $ 为其任一给定聚点。由于 $ {\mathcal{S}_k} \in \mathcal{M} $ 为有限集,不妨设存在无穷子集 $ K $ ,使 $ {{\boldsymbol{x}}_k} \to {{\boldsymbol{x}}^ * } $ $ {\mathcal{S}_k} \equiv \mathcal{S} $ $ k \in K $ 。为不失一般性,不妨设由SQP-ASA算法产生的迭代点 ${{\boldsymbol{x}}_{k + 1}} = {{\boldsymbol{x}}_k} + {\beta _k}{{\boldsymbol{d}}_k}$ $ (\forall $ $ k \in K) $ ,而由引理3可知, $ \tilde \phi ({{\boldsymbol{x}}_k},\mu ) $ 单调下降,从而由 $ {{\boldsymbol{x}}_k} \to {{\boldsymbol{x}}^ * } $ $k \in K$ 可以得到 $\tilde \phi ({{\boldsymbol{x}}_k},\mu ) \to \tilde \phi ({{\boldsymbol{x}}^ * },\mu )$ $k \to \infty$

一方面,由于步长 $ \beta $ 是由Armijo线搜索得到的,所以 $ \beta \in \left( {0,\left. 1 \right]} \right. $ 。另一方面,通过引理3知 ${{\boldsymbol{d}}_k}$ 是下降方向,不难得到

$ \begin{aligned} 0 = &\mathop {\lim }\limits_{k \in K} \;(\tilde \phi ({{\boldsymbol{x}}_{k + 1}},\mu ) - \tilde \phi ({{\boldsymbol{x}}_k},\mu )) \leqslant \mathop {\lim }\limits_{k \in K} \,\,\beta \,\nabla \tilde \phi {({{\boldsymbol{x}}_k},\mu )^{\text{T}}}{{\boldsymbol{d}}_k} \leqslant \\& \mathop {\lim }\limits_{k \in K} \,\,\beta \,( - {\boldsymbol{d}}_k^{\text{T}}{{\boldsymbol{H}}_k}{{\boldsymbol{d}}_k}) \leqslant 0 \end{aligned} $

故可知 ${{\boldsymbol{d}}_k} \to 0$ $ k \in K $ 。从而 ${\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}^ * } - {{\boldsymbol{b}}_i} = 0$ $ i \in \mathcal{S} $ $\mathcal{S} \subseteq \mathcal{M}({{\boldsymbol{x}}^ * })$ 。那么

$ \left\{ \begin{gathered} \nabla \tilde \phi ({{\boldsymbol{x}}^ * },\mu ) = \mathop \sum \limits_{i \in {\mathcal{M}_I}} \lambda _i^ * {{\boldsymbol{a}}_i} + \mathop \sum \limits_{i \in {\mathcal{M}_E}} \lambda _i^ * {{\boldsymbol{a}}_i} \\ {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}^ * } = {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_E} \\ {\boldsymbol{a}}_i^{\text{T}}{{\boldsymbol{x}}^ * } \geqslant {{\boldsymbol{b}}_i},\;i \in {\mathcal{M}_I} \\ \lambda _i^ * \geqslant 0,\;i \in {\mathcal{M}_I};\lambda _i^ * = 0,\;i \in \mathcal{M}\backslash \mathcal{S} \\ \end{gathered} \right. $

这说明 ${{\boldsymbol{x}}^*}$ 是问题(7)的一个KT点。

5 数值计算

利用Matlab R2016a实现算法SQP-ASA。其中试验环境为:Intel(R) Core(TM) i5-5200U CPU@ 2.20GHz 4.00GB RAM。通过求解下列问题的稀疏解来检验算法的有效性。

例 1  非负箱约束稀疏信号恢复的模型如下

$ \begin{aligned} \mathop {\min }\limits_{\boldsymbol{x}}& \;\;\phi (x) = \left\| {{\boldsymbol{A}}x - {\boldsymbol{b}}} \right\|_2^2 \\& {\text{s}}{\text{.t}}{\text{.}}\quad 0 \leqslant {\boldsymbol{x}} \leqslant 255 \\ \end{aligned} $

寻求该模型的稀疏解,可以转化为求解模型:

$\begin{aligned} \mathop {\min }\limits_x & \;\;\phi ({\boldsymbol{x}}) = \left\| {{\boldsymbol{A}}{\boldsymbol{x}} - {\boldsymbol{b}}} \right\|_2^2 + \left\| {\boldsymbol{x}} \right\|_q^q \\ & {\text{s}}{\text{.t}}{\text{.}}\quad 0 \leqslant {\boldsymbol{x}} \leqslant 255 \\ \end{aligned} $

分别运用本文提出的序列有效集算法(SQP-ASA)和文献[4]中的序列投影收缩算法(SQP- PG),求出该模型的稀疏解。相关参数取法如下: $ {\mu _0} = 0.1 $ ${\sigma _2} = 0.000\;1$ $ {\boldsymbol{A}} = randi ([ - 10, $ $ ,10],50,50) $ ${\boldsymbol{b}} = randi([ - 10,10],50,1)$ ${L_0} = \left\| 2 \times \right. \left. {{({\boldsymbol{A}})}^{\text{T}}}{\boldsymbol{A}} \right\|$ ${L^{\max }} = 3\;000$ ${L^{\min }} = 1\;137.5$ ,精度参数 $\epsilon = 10^{{-5}}$ ,初始点 ${{\boldsymbol{d}}_0} = 254 \times randi([0,1],50,1)$ ,最大迭代次数为1 000。目标函数与迭代次数的关系如图1所示。


图 1 迭代次数比较 Fig. 1 Comparison of iteration times

这两种方法的迭代次数和迭代时间的比较见表1。从表1可以看出,对于同样的初始点,SQP-ASA算法比SQP-PG算法的迭代次数少,运行效率也更高。


表 1 例1的数值结果 Table 1 Numerical results of example 1

例 2 寻求如下函数

$ \begin{aligned} \mathop {\min }\limits_x &\;\;\phi ({\boldsymbol{x}}) = f({\boldsymbol{x}}) \\& {\text{s}}{\text{.t}}{\text{.}}\;\;\;{\boldsymbol{x}} \in \varOmega \\ \end{aligned}$

的稀疏解。可以转化为求解模型

$ \begin{aligned} \mathop {\min }\limits_x & \;\;\phi ({\boldsymbol{x}}) = f({\boldsymbol{x}}) + \left\| {\boldsymbol{x}} \right\|_q^q \\& {\text{s}}{\text{.t}}{\text{.}}\;\;\;{\boldsymbol{x}} \in \varOmega \\ \end{aligned}$

本文通过求解以下两个模型在不同约束条件下的稀疏解,来验证算法的有效性。两个模型均来自文献[26]。

模型1  $f(x) = \dfrac{1}{2}x_1^2 - {x_1}{x_2} + x_2^2 - 6{x_1} - 2{x_2}$

$ \varOmega = \left\{ \begin{gathered} - 2{x_1} - {x_2} + 3 \geqslant 0 \\ {x_1} - {x_2} + 1 \geqslant 0 \\ - {x_1} - 2{x_2} + 2 \geqslant 0 \\ {x_1},{x_2} \geqslant 0 \\ \end{gathered} \right. $

分别使用SQP-ASA算法和SQP- PG算法,求出该模型的稀疏解。这里的相关参数取法如下: ${\;\mu _0} = 0.01$ $ \xi = 0.5 $ $ \sigma = 0.1 $ $ \eta = 1.2 $ $ {L^{\max }} = 20 $ $ {L^{\min }} = 3 $ $ {L_0} = 3 $ ,精度参数 $\epsilon{\text{=10}}^{{-5}}$ ,初始 ${{\boldsymbol{d}}_0} =$ ${[0.5\;\;\;-1]^{\text{T}}}$ ,最大迭代次数为1 000。迭代过程如图2所示。


图 2 迭代次数比较 Fig. 2 Comparison of iteration times

这两种方法的迭代次数和迭代时间的比较见表2。从表2可以看出,对于同样的初始点,SQP-ASA算法比SQP-PG算法的迭代次数少,运行效率也更高。


表 2 模型1的数值结果 Table 2 Numerical results of model 1

模型2  $f(x) = x_1^2 - {x_1}{x_2} + 2x_2^2 - {x_1} - 10{x_2}$

$ \varOmega = \left\{ \begin{gathered} - 3{x_1} - 2{x_2} + 6 \geqslant 0 \\ {x_1},{x_2} \geqslant 0 \\ \end{gathered} \right. $

分别使用SQP-ASA算法和SQP- PG算法,求出该模型的稀疏解。此时的相关参数取法如下: ${\;\mu _0} = 0.01$ $ \xi = 0.5 $ $ \sigma = 0.1 $ $ \eta = 1.2 $ $ {L^{\max }} = 20 $ $ {L^{\min }} = 2.181 $ $ {L_0} = 2.618 $ ,精度参数 $\epsilon{\text{=10}}^{{-5}}$ ,初始 ${{\boldsymbol{d}}_0 }= {[0 \;\;\; -3]^{\text{T}}}$ ,最大迭代次数为1 000。目标函数与迭代次数的关系如图3所示。


图 3 迭代次数比较 Fig. 3 Comparison of iteration times

这两种方法的迭代次数和迭代时间的比较见表3。从表3可以看出,对于同样的初始点,SQP-ASA算法比SQP-PG算法的迭代次数少,运行效率也更高。


表 3 模型2的数值结果 Table 3 Numerical results of model 2
6 结 论

提出了一种新的求解多面体约束非光滑复合函数优化问题的序列有效集算法。首先,将非光滑目标函数转化为光滑目标函数;然后,在每一迭代步中用有效集方法求解一个二次规划子问题得到方向,利用价值函数通过线搜索取得步长,重复这一过程直到求得原问题的解;最后,证明了该算法的全局收敛性。此外,通过实例进行了数值验证,结果表明,该方法较序列投影收缩方法具有一定优势。

参考文献
[1]
ZHANG C, CHEN X J. A smoothing active set method for linearly constrained non-lipschitz nonconvex optimization[J]. SIAM Journal on Optimization, 2020, 13(1): 1-30.
[2]
COMBETTES P L, PESQUET J C. Proximal splitting methods in signal processing[M]//BAUSCHKE H H, BURACHIK R S, COMBETTES P L, et al. Fixed-point Algorithms for Inverse Problems in Science and Engineering. New York: Springer, 2011.
[3]
LIU Y F, DAI Y H, LUO Z Q. Joint power and admission control via linear programming deflation[J]. IEEE Transactions on Signal Processing, 2013, 61(6): 1327-1338. DOI:10.1109/TSP.2012.2236319
[4]
LIU Y F, MA S Q, DAI Y H, et al. A smoothing SQP framework for a class of composite Lq minimization over polyhedron [J]. Mathematical Programming, 2016, 158(1/2): 467-500.
[5]
NESTEROV Y. Gradient methods for minimizing composite functions[J]. Mathematical Programming, 2013, 140(1): 125-161. DOI:10.1007/s10107-012-0629-5
[6]
BIAN W, CHEN X J. Worst-case complexity of smoothing quadratic regularization methods for non-lipschitzian optimization[J]. SIAM Journal on Optimization, 2013, 23(3): 1718-1741. DOI:10.1137/120864908
[7]
WANG W N, WU C L, TAI X C. A globally convergent algorithm for a constrained non-lipschitz image restoration model[J]. Journal of Scientific Computing, 2020, 83(1): 14. DOI:10.1007/s10915-020-01190-4
[8]
HUANG Y K, LIU H W. A Barzilai-borwein type method for minimizing composite functions[J]. Numerical Algorithms, 2015, 69(4): 819-838. DOI:10.1007/s11075-014-9927-8
[9]
MANJU V N, FRED A L. Sparse decomposition technique for segmentation and compression of compound images[J]. Journal of Intelligent Systems, 2019, 29(1): 515-528.
[10]
HU J H, JIAN H, FENG Q. A group adaptive elastic-net approach for variable selection in high-dimensional linear regression[J]. Science China Mathematics, 2018, 61(1): 173-188. DOI:10.1007/s11425-016-0071-x
[11]
SARABI M E. Primal superlinear convergence of SQP methods in piecewise linear-quadratic composite optimization[J]. Set-Valued and Variational Analysis, 2022, 30(1): 1-37. DOI:10.1007/s11228-021-00580-6
[12]
BECH A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM Journal of Imaging Sciences, 2009, 2(1): 183-202. DOI:10.1137/080716542
[13]
李启朋. 非光滑凸优化问题的快速迭代收缩阈值算法研究[D]. 西安: 西安电子科技大学, 2018.
[14]
KANZOW C, LECHNER T. Globalized inexact proximal newton-type methods for nonconvex composite functions[J]. Computational Optimization and Applications, 2021, 78(2): 377-410. DOI:10.1007/s10589-020-00243-6
[15]
CRUZ J Y B. On proximal subgradient splitting method for minimizing the sum of two nonsmooth convex functions[J]. Set-Valued and Variational Analysis, 2017, 25(2): 245-263. DOI:10.1007/s11228-016-0376-5
[16]
PATRASCU A, IROFTI P. Stochastic proximal splitting algorithm for composite minimization[J]. Optimization Letters, 2021, 15(6): 2255-2273. DOI:10.1007/s11590-021-01702-7
[17]
ZHANG X Y, BARRIO R, MARTÍNEZ M A, et al. Bregman proximal gradient algorithm with extrapolation for a class of nonconvex nonsmooth minimization problems[J]. IEEE Access, 2019, 7: 126515-126529. DOI:10.1109/ACCESS.2019.2937005
[18]
NECOARA I. General convergence analysis of stochastic first-order methods for composite optimization[J]. Journal of Optimization Theory and applications, 2021, 189(1): 66-95. DOI:10.1007/s10957-021-01821-2
[19]
张涛. 一种针对盒子约束优化问题带有新积极集策略的信赖域算法[D]. 北京: 北京交通大学, 2016.
[20]
HAGER W W W, ZHANG H C. An active set algorithm for nonlinear optimization with polyhedral constraints[J]. Science China Mathematics, 2016, 59(8): 89-106.
[21]
FLETCHER R . Practical methods of optimization[J]. SIAM Review, 1984, 26(1):143–144.
[22]
EVTUŠENKO J G. Numerical methods for nonlinear programming[J]. Doklady Akademii Nauk SSSR, 1975, 221(5): 1016-1019.
[23]
WEI B, CHEN X J, YE Y Y. Complexity analysis of interior point algorithms for non-Lipschitz and nonconvex minimization[J]. Mathematical Programming, 2014, 149(1/2): 301-327.
[24]
NESTEROV Y. Smooth minimization of non-smooth functions[J]. Mathematical Programming, 2005, 103(1): 127-152. DOI:10.1007/s10107-004-0552-5
[25]
王传芳. 解非光滑优化问题的光滑技术及理论[D]. 南京: 南京航空航天大学, 2003.
[26]
马昌凤. 最优化方法及其Matlab程序设计[M]. 北京: 科学出版社, 2010.