线性等式约束广义Lasso问题的算法研究

Algorithm Research on Generalized Lasso Problems with Linear Equality Constraints

摘要: 随着大数据时代的到来，众多研究领域都涉及到优化问题的求解，其中Lasso问题的求解尤其受到学者们的广泛研究。针对Lasso问题的求解，学者们研发出众多算法。随着应用的场景不同以及对数据的要求不同，带有约束的广义Lasso问题逐渐受到人们关注。本文将已有的快速邻近点算法结合半光滑牛顿算法，应用到对一类含线性等式约束的广义Lasso问题进行求解，并在一定的假设条件下证明了该算法的收敛性。最后，通过数值实验证实了该算法的高效性。

Abstract: With the advent of the big data era, many research fields involve solving optimization problems, among which the solution of the Lasso problem has been widely studied by scholars. Scholars have developed numerous algorithms for solving the Lasso problem. As the application scenarios vary and the data requirements differ, the generalized Lasso problem with constraints has gradually attracted attention. This paper combines existing fast proximal point algorithms with a semi-smooth Newton algorithm to solve a class of generalized Lasso problems with linear equality constraints. The convergence of the algorithm is proved under certain assumptions. Finally, the efficiency of the algorithm is verified through numerical experiments.

1. 引言

零和约束Lasso问题在生物领域有着广泛的应用背景。这种约束模型适用于许多生物学问题，尤其是在基因组学和生物信息学领域。例如：文献[1]中考虑具有零和约束的Lasso问题，如下：

$\begin{array}{l}\underset{x}{\min }\frac{1}{2}{\Vert Ax-y\Vert }^2+\lambda {\Vert x\Vert }_1\\ \text{s}\text{.t}.{\displaystyle {\sum }_{i=1}^n{x}_i}=0,\end{array}$ (1)

其中$\Vert \cdot \Vert$ 表示欧式范数，${\Vert \cdot \Vert }_1$ 表示$l_1$ 范数，$x\in R^n$ 为变量，$A\in R^{m\times n}$ ，$y\in R^m$ ，$\lambda \in R$ 均为给定量。在一些组成数据的回归模型中，即表示整体百分比或比例的数据中，需要这种特殊结构的约束。

由于约束的特殊性，该类优化问题通常出现在许多不同的领域，例如地质学、生物学、生态学和经济学。如：在微生物组分析中，数据集通常是标准化的，并产生成分数据[2] [3]。因此，设计求解该类问题的高效算法具有更加广泛的应用价值。为了使求解该类问题的算法具有广泛的应用性，现考虑求解如下更一般的等式约束广义Lasso问题：

$\begin{array}{l}\underset{x}{\min }\frac{1}{2}{\Vert Ax-y\Vert }^2+h\left(x\right)\\ \text{s}\text{.t}\text{.}\mathcal{A}\left(x\right)=b,\end{array}$ (2)

其中非光滑项函数$h\left(x\right):=\lambda {\Vert x\Vert }_1$ ，$\Vert \cdot \Vert$ 表示欧式范数，${\Vert \cdot \Vert }_1$ 表示$l_1$ 范数，$x\in R^n$ 为变量$A\in R^{m\times n}$ ，$y\in R^m$ ，$\lambda \in R$ 均为给定量；约束中的$\mathcal{A}:R^n\to R^s$ 表示线性映射；向量$b\in R^s$ 。

近年来，许多学者利用该问题的特殊结构，建立了良好的算法用来求解模型(1)。例如，Lin等人于2014年在[4]中提出了一种通过循环坐标下降的拉格朗日方法解决子问题；而Altenbuchinge等人于2017年在[5]中提出了基于变量随机选择的坐标下降策略。对于复合凸优化问题，Li [6]等人提出了一种快速邻近点算法求解复合凸优化问题模型，并基于对偶原理的半光滑牛顿算法高效稳定地求解了邻近点算法所涉及的重要子问题。

此外，考虑到约束套索的更一般形式，2016年Gaines等人在[7]中分析了一种基于二次规划的方法和ADMM方法，2020年Deng等人在[8]中提出了一种半光滑牛顿增广拉格朗日方法，文献[7] [9] [10]中设计了路径算法。

基于以上算法的启发，本文的目标是将邻近点半光滑牛顿算法推广到用于求解带线性等式约束Lasso问题(2)，并通过数值实验证实该算法求解该类带等式约束Lasso问题具有良好性能。

2. 理论知识

本节将讨论凸复合优化问题的一些稳定性性质。在之后的理论分析中可以很容易看出，这些稳定性性质是我们建立邻近点算法快速收敛的关键。

首先，给出邻近映射的定义，这将对分析邻近点算法起着至关重要的作用。

定义1 ([11] Definition 6.1) (proximal mapping)令$\mathscr{H}$ 为有限维实值希尔伯特空间，$h:R\to R\cup \left\{+\infty \right\}$ 为凸且封闭的函数。给定参数$t>0$ ，邻近映射定义如下：

${\text{prox}}_{th}\left(x\right):=\arg \underset{u\in \mathscr{H}}{\min }\left\{h\left(u\right)+\frac{1}{2t}{\Vert u-x\Vert }^2\right\},\forall x\in \mathscr{H}\text{.}$

接下来将给出Moreau包络的具体定义，如下：

定义2 ([11] Definition 6.52) (Moreau envelope)给定适当的闭凸函数$h:R\to R\cup \left\{+\infty \right\}$ 且参数$t>0$ ，则函数h的Moreau包络为：

$M_h^t\left(x\right):=\underset{u\in \mathscr{H}}{\min }\left\{h\left(u\right)+\frac{1}{2t}{\Vert u-x\Vert }^2\right\}.$

其中参数t称为光滑参数。

接下来，将给出有关非光滑性的相关分析的重要结论。

设$\mathcal{X},\mathcal{Y}$ 为两个有限维实向量空间。设$\mathcal{O}$ 为$\mathcal{X}$ 的一个开集，且$\Phi :\mathcal{O}\to \mathcal{Y}$ 为局部Lipschitz连续函数。由Rademacher定理([12] Theorem 9.60)可知，函数$\Phi (\cdot )$ 在$\mathcal{O}$ 中几乎处处F可微。用${\mathcal{D}}_{\Phi }$ 表示集合$\mathcal{O}$ 中函数$\Phi$ 均F可微的点集。用${\partial }_B\Phi \left(x\right)$ 表示函数$\Phi$ 在点$x\in \mathcal{O}$ 处的B次微分，定义如下：

${\partial }_B\Phi \left(x\right):=\left\{\underset{k\to \infty }{\lim }\nabla \Phi \left(x^k\right):x^k\in {\mathcal{D}}_{\Phi },x^k\to x\right\}.$

进一步，函数$\Phi$ 在点x处的Clarke广义Jacobian矩阵是由${\partial }_B\Phi \left(x\right)$ 的凸包定义，即$\partial \Phi \left(x\right):=co\left({\partial }_B\Phi \left(x\right)\right)$ ，并且用$\nabla \Phi \left(x^k;\Delta x\right)$ 表示函数$\Phi$ 在点x处沿着非零方向$\Delta x\in \mathcal{X}$ 的方向导数。

下面回顾一些关于半光滑性的概念。这有助于建立算法局部收敛的超线性速度。类似于[13]和[14]，给出半光滑性的定义。

定义3 (半光滑性质) [13] [14]令$F:U\subseteq R^n\to R^m$ 是局部Lipschitz连续的。令$ℳ\left(x\right)$ 为函数族，其中每个函数$M\left(x;\cdot \right)$ 为$R^n$ 到$R^m$ 的映射，并且满足$M\left(x;0\right)=0$ 。函数F被称为在点$\overline{x}\in U$ 处关于$ℳ$ 是半光滑的，若F在点$\overline{x}\in U$ 处方向可微且对于任意$M\left(\overline{x}+\text{Δ}x;\cdot \right)\in ℳ\left(\overline{x}+\text{Δ}x\right)$有

$F\left(\overline{x}+\Delta x\right)-F\left(\overline{x}\right)-M\left(\overline{x}+\Delta x;\Delta x\right)=o\left({\Vert \Delta x\Vert }_2\right),$

其中$o\left({\Vert \Delta x\Vert }_2\right)$ 表示当$\Delta x\to 0$ 时${\Vert \Delta x\Vert }_2$ 的高阶无穷小量。

同样，函数F被称为在点$\overline{x}\in U$ 处关于$ℳ$ 是强半光滑的，若F在点$\overline{x}\in U$ 处方向可微且对于任意$M\left(\overline{x}+\text{Δ}x;\cdot \right)\in ℳ\left(\overline{x}+\text{Δ}x\right)$有

$F\left(\overline{x}+\Delta x\right)-F\left(\overline{x}\right)-M\left(\overline{x}+\Delta x;\Delta x\right)=O\left({\Vert \Delta x\Vert }_2^2\right),$

其中$O\left({\Vert \Delta x\Vert }_2^2\right)$ 表示当$\Delta x\to 0$ 时${\Vert \Delta x\Vert }_2^2$ 的同阶无穷小量。

在上述理论知识的基础上，将进行优化问题的具体算法设计。

3. 求解等式约束广义Lasso问题(2)

由于问题(2)中$l_1$ 范数的存在，使得问题不易求解，因此考虑采用邻近点算法(PPA)求解问题(2)，邻近点算法求解产生了问题(2)的PPA子问题，此时PPA子问题的目标函数中两个变量可以完全分离，又保证了其强凸性。

3.1. PPA子问题

将问题(2)的目标函数Moreau-Yosida正则化后得PPA子问题：

$\begin{array}{l}\underset{x}{\min }\frac{1}{2}{\Vert Ax-y\Vert }_2^2+h\left(x\right)+\frac{1}{2t}{\Vert x-x^{k-1}\Vert }_2^2\\ \text{s}\text{.t}\text{.}\mathcal{A}\left(x\right)=b,\end{array}$ (3)

其中$t>0$ 是正则化因子，$x^{k-1}\in R^n$ 。

由于问题(2)的目标函数是凸函数，所以加入正则项后可得目标函数

$\frac{1}{2}{\Vert Ax-y\Vert }_2^2+\lambda {\Vert x\Vert }_1+\frac{1}{2t}{\Vert x-x^{k-1}\Vert }_2^2$ 是强凸的，因此在约束准则Slater条件成立的情况下，强对偶性是成立

的，即原始问题与对偶问题的对偶间隙为零，因此可以通过求解(3)的对偶问题的最优可行解，从而得到对应的原始问题(2)的最优可行解$x^{*}$ 。

3.2. 对偶问题

首先给出PPA子问题(3)的对偶问题。通过引入松弛变量z，问题(3)等价于

$\begin{array}{l}\underset{x,z}{\min }f\left(x,z\right)=\frac{1}{2}{\Vert z\Vert }_2^2+h\left(x\right)+\frac{1}{2t}{\Vert x-x^{k-1}\Vert }_2^2\\ \text{s}\text{.t}.Ax-y-z=0,\\ \mathcal{A}\left(x\right)=b.\end{array}$ (4)

其拉格朗日函数为：

$ℒ\left(x,z;{\xi }_1,{\xi }_2\right):=\frac{1}{2}{\Vert z\Vert }_2^2+h\left(x\right)+\frac{1}{2t}{\Vert x-x^{k-1}\Vert }_2^2+\langle Ax-y-z,{\xi }_1\rangle +\langle \mathcal{A}\left(x\right)-b,{\xi }_2\rangle ,$ (5)

其中${\xi }_1$ ，${\xi }_2$ 分别是$Ax-y-z=0$ ，$\mathcal{A}\left(x\right)=b$ 对应的拉格朗日乘子。关于$x,z$ 极小化$ℒ$ 得：

$\begin{array}{l}{D}_t\left({\xi }_1,{\xi }_2\right):=\underset{x,z}{\inf }ℒ\left(x,z;{\xi }_1,{\xi }_2\right)\\ =\underset{x,z}{\inf }\frac{1}{2}{\Vert z\Vert }_2^2-\langle z,{\xi }_1\rangle +h\left(x\right)+\frac{1}{2t}{\Vert x-x^{k-1}\Vert }_2^2+\langle x,A^{\text{T}}{\xi }_1\rangle -\langle y,{\xi }_1\rangle +\langle \mathcal{A}\left(x\right)-b,{\xi }_2\rangle \\ =\frac{1}{2}{\Vert \overline{z}\Vert }_2^2-\langle \overline{z},{\xi }_1\rangle +h\left(\overline{x}\right)+\frac{1}{2t}{\Vert \overline{x}-x^{k-1}\Vert }_2^2+\langle \overline{x},A^{\text{T}}{\xi }_1\rangle -y,{\xi }_1+\langle \mathcal{A}\left(\overline{x}\right)-b,{\xi }_2\rangle \\ =\frac{1}{2}{\Vert \overline{z}\Vert }_2^2-\langle \overline{z},{\xi }_1\rangle +h\left(\overline{x}\right)+\frac{1}{2t}{\Vert \overline{x}-x^{k-1}\Vert }_2^2+\langle \overline{x},A^{\text{T}}{\xi }_1\rangle -\langle y,{\xi }_1\rangle +\langle \mathcal{A}\left(\overline{x}\right),{\xi }_2\rangle -\langle b,{\xi }_2\rangle \\ =\frac{1}{2}{\Vert \overline{z}\Vert }_2^2-\langle \overline{z},{\xi }_1\rangle +h\left(\overline{x}\right)+\frac{1}{2t}{\Vert \overline{x}-x^{k-1}\Vert }_2^2+\langle \overline{x},A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\rangle -\langle y,{\xi }_1\rangle -\langle b,{\xi }_2\rangle ,\end{array}$(6)

其中$\overline{x},\overline{z}$ 是极小化拉格朗日函数(5)中的变量$x,z$ 得到的最优解，$\overline{x},\overline{z}$ 的具体形式将通过如下的引理给出。

于是得到问题(4)的对偶问题：

$\underset{{\xi }_1,{\xi }_2}{\max }{D}_t\left({\xi }_1,{\xi }_2\right)$ (7)

引理1 设$\left({\xi }_1,{\xi }_2\right)$ 是问题(4)的对偶问题(7)的一组最优解，则存在$\left(\overline{x},\overline{z}\right)$ 是PPA子问题的等价问题(4)的一组最优解，且满足如下等式：

$\{\begin{array}{l}\overline{x}={\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right),\\ \overline{z}={\xi }_1.\end{array}$(8)

证明：下面将给出具体的求解过程，分别极小化(5)中的广义拉格朗日函数中的变量x和z，可以得到

$\begin{array}{c}\overline{x}=\arg \underset{x}{\min }ℒ\left(x,z;{\xi }_1,{\xi }_2\right)\\ =\arg \underset{x}{\min }h\left(x\right)+\frac{1}{2t}{\Vert x-x^{k-1}\Vert }_2^2+\langle x,A^{\text{T}}{\xi }_1\rangle +\langle x,{\mathcal{A}}^{*}\left({\xi }_2\right)\rangle \\ =\arg \underset{x}{\min }\left\{h\left(x\right)+\frac{1}{2t}{\Vert x+t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)-x^{k-1}\Vert }_2^2\right\}\\ ={\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right),\end{array}$(9)

$\begin{array}{c}\overline{z}=\arg \underset{z}{\min }ℒ\left(x,z;{\xi }_1,{\xi }_2\right)\\ =\arg \underset{z}{\min }\left\{\frac{1}{2}{\Vert z\Vert }_2^2-\langle z,{\xi }_1\rangle \right\}\\ ={\xi }_1,\end{array}$ (10)

其中(9)式中倒数第二个等式的解是$h\left(x\right)$ 在$x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)$ 处的邻近算子。

□

由引理1可得原始变量$x,z$ 都可以由最优对偶变量${\xi }_1,{\xi }_2$ 表示，所以问题(4)的对偶问题(7)的目标函数$D_t\left({\xi }_1,{\xi }_2\right)$ 中关于变量x的部分

$\begin{array}{l}h\left(\overline{x}\right)+\frac{1}{2t}{\Vert \overline{x}+t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)-x^{k-1}\Vert }_2^2\\ =h\left({\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)\right)\\ +\frac{1}{2t}{\Vert {\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)+t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)-x^{k-1}\Vert }_2^2\\ =M_h^t\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)\end{array}$(11)

显然成立，其中(11)式中的第二个等式是根据定义2得到的。因此，函数(6)式可以简化为：

$\begin{array}{l}{D}_t\left({\xi }_1,{\xi }_2\right):=M_h^t\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)+\langle A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right),x^{k-1}\rangle \\ -\frac{1}{2}{\Vert {\xi }_1\Vert }_2^2-\frac{t}{2}{\Vert A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\Vert }_2^2-\langle y,{\xi }_1\rangle -\langle b,{\xi }_2\rangle .\end{array}$

由引理1可知，PPA子问题的等价问题(4)的解$x,z$ 可以由对偶变量${\xi }_1,{\xi }_2$ 表示，因此接下来将给出求解问题(4)的邻近点算法的框架以及算法的收敛性分析。

3.3. 邻近点算法设计

下面算法将给出求解问题(4)的邻近点算法。

算法1. 邻近点算法(PPA)求解问题(4)

1	初始化：$\left(x^0,z^0\right)\in R^n\times R^n,ϵ>0,t^0>0$ ；
2	for $j=0,1,2,\cdots$ do
3	计算：$\left({\xi }_1^{j+1},{\xi }_2^{j+1}\right)\approx \arg \max D_{t^j}\left({\xi }_1,{\xi }_2\right)$ ；
4	计算：$\left(x^{j+1},z^{j+1}\right)$ 的更新由(8)式得到；
5	if $\max \left\{{\Vert x^{j+1}-x^j\Vert }_2^2,{\Vert z^{j+1}-z^j\Vert }_2^2\right\}\le ϵ$ then
6	输出：$x^{j+1},z^{j+1}$ ；
7	else
8	更新：$t^{j+1}$ ；
9	end
10	end

对于算法1中$x^{j+1},z^{j+1}$ 的更新的具体形式已在引理1中详细给出。

3.4. 邻近点算法收敛性分析

本小节将给出算法1的收敛结果。在此之前，提出一些假设和命题。将与问题(4)相关的Moreau包络定义为

$f_t\left(x^{k-1}\right):={\min }_{x,z}\left\{\frac{1}{2}{\Vert z\Vert }_2^2+\lambda {\Vert x\Vert }_1+\frac{1}{2t}{\Vert x-x^{k-1}\Vert }_2^2:Ax-y-z=0,\mathcal{A}\left(x\right)=b\right\}.$ (12)

命题1 ([15] Proposition 2.2)当对偶问题(7)的极大值在点$\left({\xi }_1^{*},{\xi }_2^{*}\right)$ 处取得，则以下条件等价：

1) $\left(x^{*},z^{*}\right)$ 是优化问题(4)的极小值点；

2) $f\left(x^{*},z^{*}\right)=f_t\left(x^{k-1^{*}}\right)$ ，其中$f\left(x^{*},z^{*}\right)$ 和$f_t\left(x^{k-1^{*}}\right)$ 分别是问题(4)和(12)式的极小值。

命题1的成立，说明问题(4)取到最优点$\left(x^{*},z^{*}\right)$ 和(12)式取到极值点$x^{k-1^{*}}$ 是同时成立的，且此时原问题和对偶问题的目标函数值也相等。

在分析算法1的收敛性和收敛速度之前，首先给出如下假设。

假设1 问题(4)的最优解集$V^{*}$ 非空，${\displaystyle {\sum }_{j=0}^{\infty }{t}^j}=+\infty$ 且存在参数$\eta >0$ ，$ϵ>0$ 和$\omega \ge 1$ 使得

$f\left(x^{*},z^{*}\right)+\eta {\left(d\left(x,z\right)\right)}^{\omega }\le f\left(x,z\right),\forall \left(x,z\right)\in R^n\times R^n,$

且

$d\left(x,z\right)\le ϵ,$

其中

$d\left(x,z\right)=\underset{\left(x^{*},z^{*}\right)\in V^{*}}{\min }{\left({\Vert x-x^{*}\Vert }_2^2+{\Vert z-z^{*}\Vert }_2^2\right)}^{\frac{1}{2}}.$

定理1 [16]设$\left\{\left(x^j,z^j\right)\right\}$ 是由邻近点算法1产生的序列。若$t^j$ 满足${\displaystyle {\sum }_{j=0}^{\infty }{t}^j}=+\infty$ ，则有

$f\left(x^j,z^j\right)\downarrow f\left(x^{*},z^{*}\right),$

且若问题(4)的最优解集$V^{*}$ 非空，则序列$\left\{\left(x^j,z^j\right)\right\}$ 收敛到问题(4)的最优解$\left(x^{*},z^{*}\right)$ 。

假设1成立时，由定理1可知算法1产生的序列$\left\{\left(x^j,z^j\right)\right\}$ 收敛到$\left(x^{*},z^{*}\right)$ ，如下的定理可说明算法1是线性收敛的。

定理2 [16]令假设1中的$\omega =2$ 且对于任意的j，都有$\left(x^j,z^j\right)\ne \left(x^{*},z^{*}\right)$ 成立，于是有：

1) 若${\lim }_{j\to \infty }{t}^j=c$ ，$c\in \left(0,+\infty \right)$ ，则

$\lim {\sup }_{j\to \infty }\frac{d\left(x^{j+1},z^{j+1}\right)}{d\left(x^j,z^j\right)}\le \frac{1}{1+\eta c};$

2) 若${\lim }_{j\to \infty }{t}^j=+\infty$ ，则

$\lim {\sup }_{j\to \infty }\frac{d\left(x^{j+1},z^{j+1}\right)}{d\left(x^j,z^j\right)}=0.$

4. PPA子问题(3)的求解

已知问题(2)的求解实际上等价于求解一系列PPA子问题(3)，以及PPA算法具有线性收敛性。由于PPA子问题(3)的等价问题(4)和其对偶问题(7)是零对偶间隙的，且等价问题(4)的解$x,z$ 可以由对偶变量${\xi }_1,{\xi }_2$ 表示。因此，设计高效的算法求解对偶问题(7)非常关键，本节将进一步分析并设计求解对偶问题(7)的高效算法。

为便于分析，首先将对偶问题(7)写成如下等价形式：

$\begin{array}{c}\underset{{\xi }_1,{\xi }_2}{\min }\psi \left({\xi }_1,{\xi }_2\right):=-D_t\left({\xi }_1,{\xi }_2\right)\\ =-M_{\lambda {\Vert \cdot \Vert }_1}^t\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)+\frac{1}{2}{\Vert {\xi }_1\Vert }_2^2-\langle A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right),x^{k-1}\rangle \\ +\frac{t}{2}{\Vert A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\Vert }_2^2+\langle y,{\xi }_1\rangle +\langle b,{\xi }_2\rangle .\end{array}$ (13)

注意到，此时问题(13)的目标函数$\psi \left({\xi }_1,{\xi }_2\right)$ 是凸的且连续可微。

4.1. 梯度

为方便表示，记$\left({\xi }_1,{\xi }_2\right)=u$ ，$\left(\Delta {\xi }_1,\Delta {\xi }_2\right)=\Delta u$ 。

问题(13)的目标函数$\psi \left(u\right)$ 关于变量${\xi }_1,{\xi }_2$ 的梯度${\nabla }_{{\xi }_1}\psi \left(u\right),{\nabla }_{{\xi }_2}\psi \left(u\right)$ 具体形式为：

$\begin{array}{c}{\nabla }_{{\xi }_1}\psi \left(u\right)=A\left[x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)-{\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)\right]\\ +{\xi }_1-A{x}^{k-1}+tA{A}^{\text{T}}{\xi }_1+tA\left({\mathcal{A}}^{*}\left({\xi }_2\right)\right)+y\\ =-A{\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)+{\xi }_1+y,\end{array}$ (14)

$\begin{array}{c}{\nabla }_{{\xi }_2}\psi \left(u\right)=\mathcal{A}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)-{\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)\right)\\ -\mathcal{A}\left(x^{k-1}\right)+t\mathcal{A}\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)+b\\ =-\mathcal{A}{\text{prox}}_{th}\left(x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right)\right)+b,\end{array}$ (15)

因此问题(13)的目标函数$\psi \left(u\right)$ 的梯度$\nabla \psi \left(u\right)$ 为：

$\nabla \psi \left(u\right):=\left[{\nabla }_{{\xi }_1}\psi \left(u\right);{\nabla }_{{\xi }_2}\psi \left(u\right)\right].$ (16)

4.2. 广义Hessian矩阵

根据已知梯度，现在进一步分析问题(13)的目标函数$\psi \left(u\right)$ 的广义Hessian矩阵。

考虑到$\nabla \psi \left(u\right)$ 是Lipschitz连续的，则在Clarke意义下$\nabla \psi$ 在点$\left({\xi }_1,{\xi }_2\right)$ 处的广义Jacobian矩阵是B-次微分${\partial }_B\left(\nabla \psi \left(u\right)\right)$ 的凸包，即

$\partial \left(\nabla \psi \left(u\right)\right)=co\left({\partial }_B\left(\nabla \psi \left(u\right)\right)\right).$

记

$\partial \left(\nabla \psi \left(u\right)\right)\left[\left(\Delta u\right)\right]:=\left(\partial \left({\nabla }_{{\xi }_1}\psi \left(u\right)\right)\left[\left(\Delta u\right)\right],\partial \left({\nabla }_{{\xi }_2}\psi \left(u\right)\right)\left[\left(\Delta u\right)\right]\right).$ (17)

对(17)式的分析如下：

为了方便表示，记

$H_{{\xi }_1,{\xi }_2}:=A^{\text{T}}\Delta {\xi }_1+{\mathcal{A}}^{*}\left(\Delta {\xi }_2\right),$ (18)

$\widehat{x}\left({\xi }_1,{\xi }_2\right):=x^{k-1}-t\left(A^{\text{T}}{\xi }_1+{\mathcal{A}}^{*}\left({\xi }_2\right)\right),$(19)

由于${\Vert x\Vert }_1$ 在$x=0$ 处不可导，因此邻近算子的导数需要使用次梯度来表示。当x不等于0时，邻近算子的次梯度为其导数。当$x=0$ 时，邻近算子的次梯度可以是任何介于$\left[-1,1\right]$ 的值。因此，此处有

${\left[\partial {\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}\left(\widehat{x}\left({\xi }_1,{\xi }_2\right)\right)\left[H_{{\xi }_1,{\xi }_2}\right]\right]}_i=\{\begin{array}{ll}{\left[H_{{\xi }_1,{\xi }_2}\right]}_i,\hfill & {\left[\widehat{x}\left({\xi }_1,{\xi }_2\right)\right]}_i>t\lambda ,\hfill \\ a{\left[H_{{\xi }_1,{\xi }_2}\right]}_i且a\in \left[0,1\right],\hfill & {\left[\widehat{x}\left({\xi }_1,{\xi }_2\right)\right]}_i=t\lambda ,\hfill \\ 0，\hfill & {\left[\widehat{x}\left({\xi }_1,{\xi }_2\right)\right]}_i<t\lambda .\hfill \end{array}$

其中${\left[H_{{\xi }_1,{\xi }_2}\right]}_i$ 和${\left[\widehat{x}\left({\xi }_1,{\xi }_2\right)\right]}_i$ 中的下标$i=1,\cdots ,n$ 表示对应向量的第i个元素。于是有：

$\partial \left(\nabla \psi \left(u\right)\right)\left[\left(\Delta u\right)\right]=\left[tA\partial {\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}\left(\widehat{x}\left({\xi }_1,{\xi }_2\right)\right)\left[H_{{\xi }_1,{\xi }_2}\right]+\Delta {\xi }_1,t\mathcal{A}\partial {\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}\left(\widehat{x}\left({\xi }_1,{\xi }_2\right)\right)\left[H_{{\xi }_1,{\xi }_2}\right]\right].$

4.3. 最优性条件

已知问题(13)是一个凸的无约束优化问题，现通过如下的定理给出该问题的一阶最优性条件。

定理3 假设最优点对$\left({\xi }_1^{*},{\xi }_2^{*}\right):=u^{*}$ 是问题(13)的全局最优解，则满足如下条件：

$\{\begin{array}{l}{\nabla }_{{\xi }_1}\psi \left(u^{*}\right)=0,\\ {\nabla }_{{\xi }_2}\psi \left(u^{*}\right)=0.\end{array}$ (20)

注意到(20)式可以写成

$\varphi \left(u^{*}\right):=\left({\nabla }_{{\xi }_1}\psi \left(u^{*}\right),{\nabla }_{{\xi }_2}\psi \left(u^{*}\right)\right)=0,$ (21)

在点$\left({\xi }_1^k,{\xi }_2^k\right):=u^k$ 处，定义问题(13)的KKT误差为

${\eta }_{kkt}\left(u^k\right):={\Vert \varphi \left(u^k\right)\Vert }_{\infty }.$ (22)

因此，问题(13)一阶最优性条件成立等价于

${\eta }_{kkt}\left(u^k\right)=0,$

所以可以用KKT误差${\eta }^k:={\eta }_{kkt}\left(u^k\right)$ 来衡量当前迭代点的最优性。

接下来给出问题(13)的二阶充分性条件。

命题2 [17]对于函数$\nabla \psi$ 在点u处的广义Jacobian矩阵$\widehat{H}\in {\partial }_B\left(\nabla \psi \left(u^{*}\right)\right)$ 。

1) 若$u^{*}$ 为问题(13)的局部极小值点，则对任意非零s有

$s^{\text{T}}\widehat{H}s\ge 0;$

2) 令$u^{*}$ 满足KKT条件(20)。若对任意非零s有二阶充分性条件

$s^{\text{T}}\widehat{H}s\ge 0$

成立，则$u^{*}$ 为问题(13)的严格局部极小值点。

4.4. 半光滑牛顿算法

由于半光滑牛顿算法的局部快速收敛性，极大提高了子问题的求解效率，进一步可以提高算法1的整体效率。在本小节中，考虑采用半光滑牛顿算法求解子问题(13)。

由于$l_1$ 范数的邻近算子是局部Lipschitz连续的，因此可以定义${\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}$ 的Clarke广义Jacobian，记为$\partial {\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}$ 。进一步定义目标函数$\psi$ 的广义Hessian，即$\nabla \psi$ 的Clarke广义Jacobian，记为${\partial }^2\psi$ 。由于${\partial }^2\psi$ 没有显式表达式，因此定义如下的替代映射：

${\widehat{\partial }}^2\psi \left({\xi }_1,{\xi }_2\right):=\left[tA\partial {\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}\left(\widehat{x}\left({\xi }_1,{\xi }_2\right)\right)+\Delta {\xi }_1,t\mathcal{A}\partial {\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}\left(\widehat{x}\left({\xi }_1,{\xi }_2\right)\right)\right].$

对于给定的$\left({\xi }_1,{\xi }_2\right)$ ，有

${\partial }^2\psi \left({\xi }_1,{\xi }_2\right)\left(d\right)={\widehat{\partial }}^2\psi \left({\xi }_1,{\xi }_2\right)\left(d\right),$

所以，${\widehat{\partial }}^2\psi$ 可以被认为是${\partial }^2\psi$ 的一个较好的替代映射。对于给定的$\left({\xi }_1,{\xi }_2\right)$ ，令

$\widehat{H}=\left(\begin{array}{c}tAU{A}^{\text{T}}+\Delta {\xi }_1\\ t\mathcal{A}U{\mathcal{A}}^{*}\end{array}\right)$

其中$U\in \partial {\text{prox}}_{t\lambda {\Vert \cdot \Vert }_1}\left(\widehat{x}\left({\xi }_1,{\xi }_2\right)\right)$，则$\widehat{H}\in {\widehat{\partial }}^2\psi$ 。

具体来说，广义牛顿方向$\widehat{d}:=\left({\widehat{d}}_{{\xi }_1},{\widehat{d}}_{{\xi }_2}\right)$ 可以通过求解无约束的二次规划

${\min }_{\widehat{d}}\langle \nabla \psi \left(u^k\right),\widehat{d}\rangle +\frac{1}{2}\langle d,{\widehat{H}}^k\widehat{d}\rangle$ (23)

得到，其中$u^k=\left({\xi }_1^k,{\xi }_2^k\right)$ ，${\widehat{H}}^k$ 为问题(13)的目标函数$\psi$ 在点$u^k$ 处的广义Hessian矩阵。

根据目标函数的性质以及算法的理论，将问题转化为求解如下牛顿方程：

${\widehat{H}}^{k}d=-\nabla \psi \left(u^k\right).$ (24)

由于${\widehat{\partial }}^2\psi \left(u^k\right)$ 中的映射是对称和半正定的，但它们仍然可能是奇异的。对此，设置

$H^k={\widehat{H}}^k+{\nu }^{k}id,$ (25)

其中${\widehat{H}}^k\in {\widehat{\partial }}^2\psi \left(u^k\right)$ ，$id:R^m\times R^s\to R^m\times R^s$ 为恒等映射，${\nu }^k$ 定义为${\nu }^k:=\min \left(\overline{\delta },\delta {\eta }^k\right)$ ，$\overline{\delta }$ 和$\delta$ 为两个正参数，并且${\eta }^k:={\eta }_{kkt}\left(u^k\right)$ 为(21)式中定义的在点$u^k$ 处的KKT误差。显然，当${\eta }^k\ne 0$ 时，$H^k$ 为正定的。

考虑使用半光滑牛顿算法求解子问题(13)，生成迭代为

$u^{k+1}=u^k-{\alpha }^k{\left(H^k\right)}^{-1}\left[\left(\nabla \psi \left(u^k\right)\right)\right],$ (26)

其中${\alpha }^k>0$ 为第k步合适的步长，${\left(H^k\right)}^{-1}$ 为广义Hessian矩阵$H^k$ 的逆。

记

$\left(\Delta {\xi }_1{}_N^k,\Delta {\xi }_2{}_N^k\right):=-{\left(H^k\right)}^{-1}\left[\left({\nabla }_{{\xi }_1}\psi \left(u^k\right),{\nabla }_{{\xi }_2}\psi \left(u^k\right)\right)\right],$

则$\left(\Delta {\xi }_1{}_N^k,\Delta {\xi }_2{}_N^k\right)$ 为在$u^k$ 处得到的正则化广义牛顿方向，它可以通过求解广义牛顿方程：

$H^k\left[\left(\Delta {\xi }_1{}_N^k,\Delta {\xi }_2{}_N^k\right)\right]=-\left[\left({\nabla }_{{\xi }_1}\psi \left(u^k\right),{\nabla }_{{\xi }_2}\psi \left(u^k\right)\right)\right]$ (27)

得到。

对此，考虑到计算成本，通常采用共轭梯度(CG)算法来求解广义牛顿方程(27)的近似解。现在给出半光滑牛顿算法求解系统(24)的基本框架。

算法2. 半光滑牛顿算法求解问题(24)

1	Input：$\epsilon >0,\sigma \in \left(0,1\right),\tau \in \left(0,1\right)$ ；
2	Given：$\left({\xi }_1^0,{\xi }_2^0\right)\in R^m\times R^s$ ；
3	Compute $d_{{\xi }_1}^0$ and $d_{{\xi }_2}^0$ from (23)
4	for $k=0,1,\cdots$ do
5	if $\max \left\{{\eta }^k,ga{p}^k\right\}\le \epsilon$ then
6	return ${\xi }_1^k$ and ${\xi }_2^k$ ；
7	end；
8	Compute $d_{{\xi }_1}^k$ and $d_{{\xi }_2}^k$ by solving (23)；
9	Set $\alpha =1$ ；
10	while $\psi \left({\xi }_1^k+{\alpha }^k{d}_{{\xi }_1}^k,{\xi }_2^k+{\alpha }^k{d}_{{\xi }_2}^k\right)-\psi \left({\xi }_1^k,{\xi }_2^k\right)>\alpha \sigma {\Delta }^k$ do
11	Set $\alpha =\tau \alpha$ ；
12	end
13	Set ${\alpha }^k=\alpha$ and ${\xi }_1^{k+1}={\xi }_1^k+{\alpha }^k{d}_{{\xi }_1}^k$ and ${\xi }_2^{k+1}={\xi }_2^k+{\alpha }^k{d}_{{\xi }_2}^k$ ；
14	end

在下文中，交代了算法的终止条件。与[18]类似，考虑使用gap来表示原问题和对偶问题的间隙，以便于测量原目标与对偶目标之间的差距，即

$gap=\frac{\left|probj-dobj\right|}{1+\left|probj\right|+\left|dobj\right|},$

其中probj和dobj分别表示原问题目标函数值和对偶问题目标函数值。

当$\max \left\{{\eta }^k,ga{p}^k\right\}\le \epsilon$ 时，算法终止，其中$\epsilon >0$ 为给定误差，${\eta }^k$ 为(22)式中定义的KKT误差，gap为点$u^k$ 处的间隙。因此，在算法2中给出了求解问题(23)的算法。

4.5. 半光滑牛顿算法收敛性分

本小节将分析半光滑牛顿算法的收敛性，为方便表示，记$u:=\left({\xi }_1,{\xi }_2\right)$ 。

定义4 假设$\nabla \psi$ 是局部Lipschitz连续的，则其广义Jacobian存在。取在$u^k$ 点任意的广义Jacobian矩阵${\widehat{H}}_k\in \partial \left(\nabla \psi \left(u^k\right)\right)$ ，若${\widehat{H}}_k$ 可逆，其基本的迭代格式为

$u^{k+1}=u^k-{\widehat{H}}_k^{-1}\nabla \psi \left(u^k\right).$ (28)

假设2 映射$\nabla \psi$ 在最优点$u^{*}$ 处是半光滑和所有Jacobian矩阵是非奇异的。

引理2 如果假设2成立，则存在常数$c>0,\kappa >0$ 和一个小邻域$N\left(u^{*},{ϵ}_0\right)$ 使得对于任意的$u\in N\left(u^{*},{ϵ}_0\right)$ 和${\widehat{H}}_k\in \partial \left(\nabla \psi \left(u^k\right)\right)$ ，下面结论成立：

1) $u^{*}$ 是一个孤立解；

2) ${\widehat{H}}_k$ 是非奇异的并且$\Vert {\widehat{H}}_k^{-1}\Vert \le c$ ；

3) 局部误差届条件对于$\nabla \psi \left(u\right)$ 在邻域$N\left(u^{*},{ϵ}_0\right)$ 上成立，也就是说$\Vert u-u^{*}\Vert \le \kappa \Vert \nabla \psi \left(u\right)\Vert$ 。

如下定理给出了半光滑牛顿算法的局部二次收敛性。

定理4 设$\nabla \psi \left({\xi }_1,{\xi }_2\right)$ 具有半光滑性并且$u^{*}$ 是优化问题(13)的最优解。那么迭代(28)是良定义的，且存在一个小邻域$N\left(u^{*},ϵ\right)$ ，使得对于任意的k有$u^k\in N\left(u^{*},ϵ\right)$ ，迭代(28)是超线性收敛的。如果$\nabla \psi \left(u\right)$ 是强半光滑的，迭代(28)是二次收敛的。

证明：根据引理2，迭代(28)式是良定义的。可以推出：

$\begin{array}{c}\Vert u^{k+1}-u^k\Vert =\Vert u^k-{\widehat{H}}_k^{-1}\nabla \psi \left(u^k\right)-u^{*}\Vert \\ =\Vert {\widehat{H}}_k^{-1}\Vert \Vert \nabla \psi \left(u^k\right)-\nabla \psi \left(u^{*}\right)-{\widehat{H}}_k\left(u^k-u^{*}\right)\Vert \\ =O\left({\Vert u^k-u^{*}\Vert }^2\right),\end{array}$

其中最后一个等式来源于强半光滑性。

□

5. 数值实验

本节将通过仿真数值实验来测试本文使用的邻近点对偶半光滑牛顿算法(记为PPA_DSN)算法求解含等式约束广义Lasso问题的实际效果。算法采用MatlabR2017a编程实现，且所有实验均在AMD Ryzen 7 6800HS Creator Edition 3.20 GHz，16.0 GB内存的个人笔记本电脑上运行。同时，考虑选择邻近梯度算法(记为PPA_DG)和交替方向乘子法算法(记为ADMM_GLP)作为对比算法。算法PPA_DSN所需参数分别为：设计矩阵A的规模：$m\times n$ ；观测向量y的规模：m维；最大迭代时间600秒；最大误差$ϵ=1\times {10}^{-6}$ 。

对于所有算法，若下面任一条件满足，则终止程序：

1) 算法残差${\eta }_{kkt}\left(u^k\right)\le ϵ$ 。

2) 达到最大迭代数maxit = 200。

为进行实验观测算法PPA_DSN求解含等式约束广义Lasso问题的效率，首先设置含等式约束广义Lasso问题的样本数量为m，特征数量为n，设计矩阵A的每个分量$A_{ij}$ 均随机生成；设向量$y=A{x}^{*}+\epsilon$ ($\epsilon$ 表示随机扰动向量)，正则化参数$\lambda =0.01$ 。

Figure 1. In correspondence with the experimental data of three algorithms

图1. 三种算法对比的实验数据

为了验证PPA_DSN算法的有效性和稳定性，我们进行了不同参数设置下的实验。下文提供PPA_DSN算法、邻近梯度算法(记为PPA_DG)和交替方向乘子算法(记为ADMM_GLP)的数值结果。我们针对零和约束Lasso问题(1)，取变量维数$n=\left[1000,1500,2000\right]$ ，$m=\left[100,200,500\right]$ ，正则化参数$\lambda =\left[0.001,0.01,0.1\right]$ ，最大迭代时间600秒，最大误差。在这些设置下，共进行了27组实验。考虑到实验数据均为随机生成的，每组实验重复了3次，并取数值结果的平均值。具体结果详见图1，其中Res、Iter、Iterf、Iterd和Time分别表示残差${\eta }_{kkt}\left(u^k\right)$ 、迭代步数、目标函数计算次数、搜索方向计算次数和CPU运行时间(以秒为单位)。从图1可以看出，算法PPA_DSN成功求解了所有测试问题，并且在大多数情况下，CPU耗时不超过1秒。

6. 总结

本文主要对一类含线性等式约束的广义Lasso问题进行算法研究。由于大多实际问题的求解过程中会涉及对样本数据的限制，如零和约束Lasso问题，因此考虑求解一类含线性等式约束的广义Lasso问题。受众多研究者先前对Lasso模型的算法研究的启发，本文考虑将邻近点半光滑牛顿算法推广到用于求解含线性等式约束Lasso问题，并通过数值实验与邻近梯度算法以及交替方向乘子法进行对比，证实了该算法求解该类含线性等式约束Lasso问题具有良好性能。