有效集方法求解欠定线性方程组的稀疏非负解

Sparse Nonnegative Solution of Underdetermined Linear Equations by Active Set Method

摘要: 针对欠定线性方程组稀疏非负解的求解问题，本文首先将原问题松弛为l₀正则优化模型。随之提出有效集方法识别严格L-稳定点邻域内的零分量，基于这种快速识别技术，设计了有效集Barzilar-Borwein算法求解l₀正则极小化模型。最后的数据实验证明该算法可以快速有效地求解欠定线性方程组的稀疏非负解。

Abstract: In this paper, for acquiring the sparse nonnegative solution of underdetermined linear equations, the original problem is relaxed into a l₀ regularized optimization model. An active set identification technique is developed to accurately identify the zero components in a neighbourhood of the strict L-stationary point. Based on the active set identification technique, we propose an active set Barzilar-Borwein method to solve a l₀ regularized minimization model. Numerical results show that the algorithm can effectively solve the sparse nonnegative solutions of underdetermined linear equations.

1. 引言

基于欠定线性方程组 $Ax=b$，这里 $b\in R^m$， $x\in R^n$，A是 $m\times n\left(m<n\right)$ 的矩阵，b是已知的，x是未知的，本文只考虑非负解 $x\ge 0$ 的部分。寻找欠定线性方程组的稀疏非负解，就是求解如下的优化问题：

$\begin{array}{l}\min {\Vert x\Vert }_0\\ \text{s}\text{.t}\text{.}Ax=b,x\ge 0.\end{array}$ (1)

其中 ${\Vert x\Vert }_0$ 指向量x中非零分量的个数。这类问题在信号处理，图像恢复和多光谱数据处理等领域引起了广泛的研究 [1] [2] [3]。

在现有的研究成果中，主要用贪婪法和松弛法来解决问题(1)。贪婪法主要包括正交匹配追踪算法 [4]，正则化正交匹配追踪算法 [5] 以及压缩采样匹配追踪 [6] 和子空间追踪 [7]。这种类型的算法由于其相对较快的计算速度，而在实践中得到了广泛的应用。然而，当信号不是很稀疏或观测噪声水平较高时，贪婪法的性能无法得到保证。

松弛法是通过用一定的非负连续函数代替 $l_0$ 范数，将 $l_0$ 极小化模型转化为更易于处理的模型。常用于逼近 $l_0$ 极小化模型的是 $l_{\text{1}}$ 极小化模型。 $l_{\text{1}}$ 极小化是一个凸优化问题，在矩阵A满足适当的条件下， $l_{\text{1}}$ 极小化可以精确地恢复稀疏信号。因此， $l_{\text{1}}$ 极小化被认为是解决稀疏性问题的非常有用的工具并得到了广泛的运用。此外，一些非凸函数 $l_p\left(0<p<1\right)$ 范数被提出作为 $l_0$ 范数的松弛。与 $l_{\text{1}}$ 极小化方法相比，非凸松弛模型通常具有更好的稀疏性，但通常较难求解。

求解约束优化问题(1)主要有两种算法。第一种算法是迭代重加权算法，其中两个最重要的迭代重加权算法是重加权 $l_{\text{1}}$ 极小化 [8] 和迭代重加权最小二乘算法 [9]。另一种方法通常称为正则化方法，它通过引入正则化参数将问题(1)转化为下面的优化问题：

$\begin{array}{l}\min \varphi \left(x\right):=f\left(x\right)+\mu {\Vert x\Vert }_0\\ \text{s}\text{.t}\text{.}x\ge 0.\end{array}$ (2)

其中 $f\left(x\right)=\frac{1}{2}{\Vert Ax-b\Vert }_2^2$，且 $\mu >0$ 是给定的正则化参数。

本文的其他内容如下：首先给出一些符号定义和关于问题(2)的一阶优化条件。基于Cheng和Chen [10] 的研究，我们在加入非负约束 $x\ge 0$ 后，重新定义了四种关于问题(2)的稳定点：基本稳定点，强基本稳定点，L-稳定点，严格L-稳定点，并说明了它们之间的关系。随之提出了有效集方法求解问题(2)，这种方法能够精准地找出L-稳定点邻域内的值为零的分量。接着提出了有效集Barzilar-Borwein算法求解问题(2)。最后进行了数据实验，证明了该算法可以有效的求解欠定线性方程组的稀疏非负解。

2. 符号说明

在本节中，我们将定义一些符号说明。假设 $x^{*}$ 是问题(2)的最优解，我们定义有效集 $\gamma \left(x^{*}\right)$ 是 $x^{*}$ 中零分量的下标集合，无效集 $\tau \left(x^{*}\right)$ 是 $x^{*}$ 中非零分量的下标集合：

$\gamma \left(x^{*}\right)=\left\{i:x_i^{*}=0\right\}和\tau \left(x^{*}\right)=\left\{i:x_i^{*}>0\right\}.$

定义2.1 [11] 令 $\nu \in R^n$， $\alpha >0$，若向量值函数 $h_0\left(\nu ,\alpha \right)$ 满足

$h_0{\left(\nu ,\alpha \right)}_i=\{\begin{array}{l}{\nu }_i,\left|{\nu }_i\right|>\sqrt{2\alpha },\\ 0,\left|{\nu }_i\right|\le \sqrt{2\alpha }.\end{array}i=1,2,\cdots ,n.$

则称 $h_0\left(\nu ,\alpha \right)$ 为硬阈值算子。

注2.1 硬阈值算子是 $l_0$ 极小化问题的最优解，即

$h_0\left(\nu ,\alpha \right)=\underset{u\in R^n}{\arg \min }\left\{{\Vert u\Vert }_1:u=\underset{x\in R^n}{\arg \min }\left\{\alpha {\Vert x\Vert }_0+\frac{1}{2}{\Vert x-\nu \Vert }^2\right\}\right\}.$

在最近的研究中，Beck和Hallak [12] 提出了下面问题的一些基本定义

(3)

其中 $\mu >0$， $B\subseteq R^n$。

定义2.2 (基本稳定点)如果 $B=R_{+}^n$，则有：

$g_i\left(x^{*}\right)=0,\forall i\in \tau \left(x^{*}\right).$ (4)

其中 $g_i\left(x^{*}\right)$ 表示 $f\left(x^{*}\right)$ 在第i个分量的梯度。那我们称满足条件(4)的 $x^{*}$ 为问题(2)的基本稳定点。

定义2.3 (L-稳定点)向量 $x^{*}$ 称作问题(3)的L-稳定点需满足：

$\begin{array}{l}{x}^{*}\in pro{x}_{\frac{\mu }{L}{g}_B}\left(x^{*}-\frac{1}{L}g\left(x^{*}\right)\right)\\ =\underset{x\in B}{\arg \min }\left\{\frac{\mu }{L}{\Vert x\Vert }_0+\frac{1}{2}{\Vert x-\left(x^{*}-\frac{1}{L}g\left(x^{*}\right)\right)\Vert }^2\right\}.\end{array}$

其中 $L>0$， $g_B\left(x\right)={\delta }_B\left(x\right)+{\Vert x\Vert }_0$ 且 ${\delta }_B\left(x\right)$ 是集合B的指标函数，而指标函数也就是

${\delta }_B\left(x\right)=\{\begin{array}{l}0,x\in B;\\ 1,x\notin B.\end{array}$

当 $B=R_{+}^n$ 时，问题 等价于问题 ，则上述条件等价于下面的条件：

$\{\begin{array}{l}\text{0}\le g_i\left(x^{*}\right)\le \sqrt{2\mu L},i\in \gamma \left(x^{*}\right);\\ x_i^{*}\ge \sqrt{\frac{2\mu }{L}}且g_i\left(x^{*}\right)=0,i\in \tau \left(x^{*}\right).\end{array}$

向量 $x^{*}$ 满足如下条件就称为问题(2)的严格L-稳定点：

$\{\begin{array}{l}\text{0}\le g_i\left(x^{*}\right)\le \sqrt{2\mu L},i\in \gamma \left(x^{*}\right);\\ x_i^{*}>\sqrt{\frac{2\mu }{L}}且g_i\left(x^{*}\right)=0,i\in \tau \left(x^{*}\right).\end{array}$ (5)

向量 $x^{*}$ 满足如下条件就称为问题(2)的强基本稳定点：

$\{\begin{array}{l}\text{0}\le g_i\left(x^{*}\right)\le \sqrt{2\mu L},i\in \gamma \left(x^{*}\right);\\ g_i\left(x^{*}\right)=0,i\in \tau \left(x^{*}\right).\end{array}$

由上述条件不难得出由问题(2)的严格L-稳定点可以推出问题(2)的L-稳定点，强基本稳定点和基本稳定点。

对于函数 $\varphi \left(x\right)$，称 $d\in R^n$ 是函数 $\varphi \left(x\right)$ 在x处的下降方向，如果常数 $\delta >0$，使得

$\varphi \left(x+\alpha d\right)<\varphi \left(x\right),\forall \alpha \in \left(0,\delta \right).$

我们用 $D\left(x\right)$ 表示函数 $\varphi \left(x\right)$ 的所有下降方向的集合：

$D\left(x\right)=\left\{d\in R^n|\varphi \left(x+\alpha d\right)<\varphi \left(x\right),\forall \alpha \in \left(0,\delta \right)\right\}.$

类似于文献 [10] 的定理2.1，下面的定理将给出集合 $D\left(x\right)$ 更精确的表达形式，这将有利于我们设计线性搜索方法解决问题(2)。

定理2.1 假设 $x>0$，由于 $f\left(x\right)$ 是连续可微的，集合 $D\left(x\right)$ 可表示为：

$D\left(x\right)=\left\{d\in R^n|g_{\tau \left(x\right)}^{\text{T}}{d}_{\tau \left(x\right)}<0,d_{\gamma \left(x\right)}=0\right\}.$

3. 有效集估计

通常来说，有效集的效率很大程度上取决于如何找到最优点处的有效约束。接下来，将介绍有效集的估计并用它来逼近原问题的严格L-稳定点。我们用集合 $A\left(x\right)$ 逼近 $\gamma \left(x^{*}\right)$，集合 $F\left(x\right)$ 逼近 $\tau \left(x^{*}\right)$ ：

$A\left(x\right)=\left\{i:0\le x_i\le \frac{1}{L}{g}_i\left(x\right)+\sqrt{\frac{2\mu }{L}}\right\},$

和

$F\left(x\right)=\left\{i:x_i>\frac{1}{L}{g}_i\left(x\right)+\sqrt{\frac{2\mu }{L}}\right\}.$

下面的定理表明，当点 $x^k$ 足够接近严格L-稳定点时，则上述有效集的逼近集合 $A\left(x^k\right)$ 是集合 $\gamma \left(x^{*}\right)$ 的子集。

定理3.1 假设 $x^{*}$ 是问题(2)的严格L-稳定点，序列 $\left\{x^k\right\}$ 收敛于 $x^{*}$，则对于足够大的k有：

$A\left(x^k\right)\subseteq \gamma \left(x^{*}\right)且\tau \left(x^{*}\right)\subseteq F\left(x^k\right).$

证明：假设 $x^{*}$ 是问题(2)的严格L-稳定点，如果 $\tau \left(x^{*}\right)=\varnothing$，由于

$A\left(x^k\right)\cup F\left(x^k\right)=\gamma \left(x^{*}\right)\cup \tau \left(x^{*}\right)=\left\{1,2,\cdots ,n\right\}.$

所以结论成立。假设 $\tau \left(x^{*}\right)\ne \varnothing$，对于任意的 $i\in \tau \left(x^{*}\right)$，由(5)可知：

$x_i^{*}-\frac{1}{L}{g}_i\left(x^{*}\right)=x_i^{*}>\sqrt{\frac{2\mu }{L}}.$

因为 $g\left(x\right)$ 是连续的，故有：

$x_i^k-\frac{1}{L}{g}_i\left(x^k\right)>\sqrt{\frac{2\mu }{L}}.$

因此， $i\in F\left(x^k\right)$ 也就是 $\tau \left(x^{*}\right)\subseteq F\left(x^k\right)$。

另一方面，对于任意的 $i\in A\left(x^k\right)$，有

$0\le x_i^k\le \frac{1}{L}{g}_i\left(x^k\right)+\sqrt{\frac{2\mu }{L}}.$

由于序列 $\left\{x^k\right\}$ 收敛于 $x^{*}$，可得：

$0\le x_i^{*}\le \frac{1}{L}{g}_i\left(x^{*}\right)+\sqrt{\frac{2\mu }{L}}.$

对于足够大的k，可得 $x_i^{*}=0$。再由最优性条件，我们有 $0\le g_i\left(x^{*}\right)\le \sqrt{2\mu L}$。因此， $i\in \gamma \left(x^{*}\right)$ 也就是 $A\left(x^k\right)\subseteq \gamma \left(x^{*}\right)$，证毕。

下面的定理表明，当 $x^k$ 是问题(2)的严格L-稳定点时，我们可以得到集合 $A\left(x^k\right)=\gamma \left(x^k\right)$。

定理3.2 假设是问题的严格L-稳定点，则有：

$A\left(x^k\right)=\gamma \left(x^k\right),F\left(x^k\right)=\tau \left(x^k\right).$

证明：假设 $x^k$ 是问题的严格L-稳定点，对于任意的 $i\in \gamma \left(x^k\right)$，由条件(5)可得 $i\in A\left(x^k\right)$。而对于任意的 $i\in A\left(x^k\right)$，我们有 $x_i^k=0$。但事实上，如果存在 $x_i^k\ne 0$，由条件(5)可得 $g_i\left(x^k\right)=0$ 且有：

$x_i^k\le \sqrt{\frac{2\mu }{L}}\Rightarrow x_i^k=0.$

由此可得 $i\in \gamma \left(x^k\right)$，故 $A\left(x^k\right)=\gamma \left(x^k\right)$。同理可得 $F\left(x^k\right)=\tau \left(x^k\right)$，证毕。

4. 有效集算法

在本节中，首先给出搜索方向的定义和一些相关性质，然后介绍一种有效集Barzilar-Borwein算法来求解问题(2)，并分析了该算法的全局收敛性。

4.1. 搜索方向的一些性质

令 $x^k\in \Omega$ 是第k个迭代点，为了简化表述，令

$A^k=A\left(x^k\right),F^k=\left\{i\in F\left(x^k\right):x_i>0\right\}.$

为了定义搜索方向，进一步将 $A^k$ 分成两个部分：

$A_1^k=\left\{i\in A^k:x_i^k{g}_i\left(x^k\right)\ge 0\right\},A_2^k=\left\{i\in A^k:x_i^k{g}_i\left(x^k\right)<0\right\}.$

搜索方向 $d^k=\left(d_{A_1^k}^k,d_{A_2^k}^k,d_{F^k}^k\right)$，定义 $d_{F^k}^k$ 为：

$d_i^k=-{\lambda }^k{g}_i^k,i\in F^k.$ (6)

其中 $0<{\lambda }_{\min }\le {\lambda }^k<{\lambda }_{\max }$ 是Barzilar-Borwein算法的步长，其精确的定义将在4.2节中给出。集合 $A_1^k$ 表示零分量的下标集合，且包含满足条件(5)的变量的下标集合。集合 $A_2^k$ 包含不满足条件(5)的变量的下标集合。因此，给出定义：

$d_i^k=-x_i^k,i\in A_1^k.$ (7)

$d_i^k=-g_i^k,i\in A_2^k.$ (8)

下面的定理表明 $d^k=0$ 当且仅当 $x^k$ 是问题(2)的一个严格L-稳定点。

定理4.1 假设 $d^k$ 由上述定义式(6)~(8)所表示。 $d^k=0$ 当且仅当 $x^k$ 是问题(2)的一个严格L-稳定点。

证明：令 $d^k=0$，由 $A_2^k$ 的定义和搜索方向 $d_i^k=-g_i^k=0$，可得 $A_2^k=\varnothing$。如果 $i\in A_i^k$，由 $A_1^k$ 的定义和 $d_i^k=-x_i^k=0$ 可得：

$x_i^k=0且g_i\left(x^k\right)\le \sqrt{2\mu L}.$

如果 $i\in F^k$，由 $F^k$ 的定义可知：

$x_i^k>\sqrt{\frac{2\mu }{L}}且g_i\left(x^k\right)=0.$

反过来，假设 $x^k$ 是问题(2)的严格L-稳定点。如果 $i\in A^k$，我们称 $x_i^k=0$， $d_i^k=-x_i^k=0$ 且 $A_2^k=\varnothing$。事实上，如果 $x_i^k\ne 0$，由式(5)可得：

$g_i\left(x^k\right)=0且x_i^k\le \sqrt{\frac{2\mu }{L}}\Rightarrow x_i^k=0.$

这与假设条件矛盾。因此，对 $\forall i\in F^k$ 有：

$g_i\left(x^k\right)=0且d_i^k=0.$

证毕。

定理4.2 假设 $x^k\to x^{*}$ 且 $d^k\to 0$，其中 $d^k$ 由式(6)~(8)所表示，则 $x^{*}$ 是问题(2)的一个L-稳定点。

证明：对于 $\forall i\in A_1^k$，当 $x^k\to x^{*}$ 时，由于 $d_i^k=-x_i^k\to 0$，可得 $x_i^{*}=0$。再由 $A_1^k$ 的定义，我们可得：

$x_i^{*}=0且g_i\left(x^{*}\right)\le \sqrt{2\mu L},$

对于 $\forall i\in A_2^k$，由 $d_i^k=-g_i^k\to 0$，可得：

$g_i\left(x^{*}\right)=0.$

对于 $\forall i\in F^k$，由于 $g\left(x\right)$ 的连续性和 $d^k\to 0$，有

$x_i^{*}\ge \sqrt{\frac{2\mu }{L}}且g_i\left(x^{*}\right)=0.$

证毕。

下面的定理表明，如果 $x^k$ 是问题(2)的强基本稳定点，由(6)~(8)定义的 $d^k$ 是f在 $x^k$ 处的下降方向。

定理4.3 假设 $d^k$ 由(6)~(8)式定义，可得

${\left(g^k\right)}^{\text{T}}{d}^k\le 0.$ (9)

此外，等号成立当且仅当 $x^k$ 是问题(2)的强基本稳定点。

证明：由 $A_1^k$ 的定义可知：

$g_i^k{d}_i^k=-g_i^k{x}_i^k\le 0,\forall i\in A_i^k.$

再由 $A_2^k$ 的定义可知：

$g_i^k{d}_i^k=-{\left(g_i^k\right)}^2<0.$

对于任意的 $i\in F^k$，我们有：

${\left(g_i^k\right)}^{\text{T}}{d}_i^k\le -{\lambda }_{\min }{\Vert {\left(g_i^k\right)}^{\text{T}}\Vert }^2\le 0.$

综上所述，我们可以得到 ${\left(g^k\right)}^{\text{T}}{d}^k\le 0$。

假设 $x^k$ 是问题(2)的强基本稳定点。对于任意的 $i=1,2,\cdots ,n$，有 $x_i^k{g}_i\left(x^k\right)=0$。由此，对于任意的 $i\in A_1^k$，有 $d_i^k=-x_i^k$，则有 ${\left(g^k\right)}^{\text{T}}{d}^k=0$。由 $A_2^k$ 的定义可知 $A_2^k=\varnothing$。对于任意的 $i\in F^k$，我们有 $g_i^k=0$。相反地，假设不等式(9)的左边等于0，我们有：

$g_i\left(x^k\right)=0且x_i^k>\sqrt{\frac{2\mu }{L}},\forall i\in F^k.$

由 $A_2^k$ 的定义可知 $A_2^k=\varnothing$。由的定义可得：

$x_i^k=0,0\le g_i^k\le \sqrt{2\mu L}或0<x_i^k\le \sqrt{\frac{2\mu }{L}},g_i^k=0.$

由此可知， $x^k$ 是问题 的强基本稳定点，证毕。

4.2. 算法

基于上述讨论和非单调线性搜索 [13]，类似于文献 [10]，我们提出了一种有效集Barzilar-Borwein算法来求解问题(2)。但不同于文献 [10] 的是我们在加入非负约束后，重新定义了有效集和有效集部分的搜索方向，产生了不同的迭代点，加快了算法的收敛速度。

算法4.1 (有效集Barzilar-Borwein算法，简记ABB算法)

步0 初始化 $x^0\in R_{+}^n$，取 $\eta ,\delta \in \left(0,1\right)$， ${\lambda }^0=1$ 且 ${\lambda }_{\max }>{\lambda }_{\min }>0$，正整数M，取 $k=0$。

步1 进行收敛性检验，如果满足停止条件，则用近似解 $x^k$ 终止。

步2 计算 ${\lambda }^k:{\lambda }_{\min }\le {\lambda }^k\le {\lambda }_{\max }$ 且 $d^k$ 由式(6)~(8)确定。

步3 如果，满足不等式

$\varphi \left(x^k+{\alpha }^k{d}^k\right)\le \varphi \left(x^k\right)-\delta {\left({\alpha }^k\Vert d^k\Vert \right)}^2.$ (10)

执行步5，否则执行步4。

步4 确定 ${\alpha }^k:=\max \left\{{\eta }^j,j=0,1,\cdots \right\}$ 满足不等式

$f\left(x^k+{\alpha }^k{d}^k\right)\le \underset{0\le j\le m\left(k\right)}{\max }f\left(x^{k-j}\right)-\delta {\left({\alpha }^k\Vert d^k\Vert \right)}^2.$ (11)

其中 $m\left(0\right)=0$， $0\le m\left(k\right)\le \min \left[m\left(k-1\right)+1,M\right]$，执行步5。

步5 计算 $x^{k+1}:=x^k+{\alpha }^k{d}^k$，取 $k:=k+1$，返回步1。

接下来，将给出 ${\lambda }^k$ 的显性形式，运用了Barzilar-Borwein (BB)算法 [14] 的基本思想。类似于Barzilar-Borwein步的思想，通过求解下面的优化问题：

$\underset{\lambda \in R}{\min }{\Vert D\left(\lambda \right)s_{F^k}^{k-1}-y_{F^k}^{k-1}\Vert }^2.$

其中 $s_{F^k}^{k-1}=x_{F^k}^k-x_{F^k}^{k-1}$ 且 $y_{F^k}^{k-1}=g_{F^k}^k-g_{F^k}^{k-1}$，我们得到Barzilar-Borwein步长如下：

${\lambda }_{BB}^k=\frac{{\Vert s_{F^k}^{k-1}\Vert }^2}{{\left(s_{F^k}^{k-1}\right)}^{\text{T}}{y}_{F^k}^{k-1}}.$

为了控制 ${\lambda }_{BB}^k$ 的大小值，我们用区间 $\left[{\lambda }_{\min },{\lambda }_{\max }\right]$ 来限制，其中 ${\lambda }_{\min }<{\lambda }_{\max }$ 是给定的正常数。因此，我们设

${\lambda }^k=\max \left\{{\lambda }_{\min },\min \left\{{\lambda }_{\max },{\lambda }_{BB}^k\right\}\right\}$ (12)

为了简单起见，我们将步长(12)代入式(6)，简化了算法4.1。

4.3. 全局收敛性分析

下面将分析算法4.1的全局收敛性。类似于文献 [13] 中的分析，可以得到了如下的引理。

引理4.1 假设序列是由算法4.1产生。如果存在一个正整数 $k_1$，对于所有 $k>k_1$，满足线性搜索(11)，我们有：

$\underset{k\to \infty }{\lim }{\alpha }^k{d}^k=0.$ (13)

下面的定理和定理4.2表明序列 $\left\{x^k\right\}$ 的任意聚点都是问题(2)的L-稳定点和强基本稳定点。

定理4.4 假设序列 $\left\{x^k\right\}$ 是由算法4.1产生。如果对于所有的k都有 $d^k\ne 0$，则序列 $\left\{x^k\right\}$ 的任意聚点 $x^{*}$ 都是问题(2)的L-稳定点或强基本稳定点。

证明：假设对于所有的k都有 $d^k\ne 0$，且 $x^{*}$ 为序列 $\left\{x^k\right\}$ 的任意聚点。则存在一个无限指标集合K，满足 $k\in K$ 使得：

$\underset{k\to \infty }{\lim }{x}^k=x^{*}.$

考虑到不同集合 $A_1^k,A_2^k,F^k$ 的数量是有限的，对于任意的 $k\in K$，我们考虑两种情况。情形(1)：如果存在一个无限指标集合 $K_1\subseteq K$，对于所有的 $k\in K_1$，线搜索(11)都满足。通过对不等式(10)求和，可得

则由定理4.2可知结论成立。情形(2)：如果存在一个正整数 $k_1$，对于所有的 $k>k_1$，线搜索(11)都满

足。如果序列 ${\left\{{\alpha }^k\right\}}_{k\in K}$ 有一个非零的极限值，我们由(13)式可知

$\underset{k\to \infty }{\lim }{d}^k=0.$

则由定理4.2可得知结论成立。否则，如果我们有 $\underset{k\to \infty }{\lim }{\alpha }^k=0$，对于所有的 $k\ge 1$，由线搜索(11)规则可得

$f\left(x^k+\frac{{\alpha }^k}{\eta }{d}^k\right)-f\left(x^k\right)\ge f\left(x^k+\frac{{\alpha }^k}{\eta }{d}^k\right)-\underset{0\le j\le m\left(k\right)}{\max }f\left(x^{k-j}\right)>-\delta {\left(\frac{{\alpha }^k}{\eta }\Vert d^k\Vert \right)}^2.$ (14)

另一方面，对于充分大的 $k\in K$，由中值定理可得

$f\left(x^k+\frac{{\alpha }^k}{\eta }{d}^k\right)-f\left(x^k\right)=\frac{{\alpha }^k}{\eta }{\left(g\left(x^k+{\theta }^k\frac{{\alpha }^k}{\eta }{d}^k\right)-g\left(x^k\right)\right)}^{\text{T}}{d}^k+\frac{{\alpha }^k}{\eta }g{\left(x^k\right)}^{\text{T}}{d}^k.$

其中 ${\theta }^k\in \left(0,1\right)$。将上述等式带入(14)，可得

${\left(g\left(x^k+{\theta }^k\frac{{\alpha }^k}{\eta }{d}^k\right)-g\left(x^k\right)\right)}^{\text{T}}{d}^k+g{\left(x^k\right)}^{\text{T}}{d}^k>-\delta \frac{{\alpha }^k}{\eta }{\Vert d^k\Vert }^2.$

由于序列是有界的，不失一般性地，存在 $d^{*}$ 使得 $\underset{k\to \infty }{\lim }{d}^k=d^{*}$，结合引理4.1可知， $\underset{k\to \infty }{\lim }{\alpha }^k=0$。再由最后一个不等式取极限得 $g{\left(x^{*}\right)}^{\text{T}}{d}^{*}>0$。由定理4.3得出矛盾，故结论成立，证毕。

5. 数值实验

下面将进行一些数值实验来验证算法4.1对于求解欠定线性方程组的稀疏非负解是有效的。所有的实验结果都是在联想笔记本电脑(2.50 GHz, 8.00 GB of RAM)上实现的，数学软件采用MATLAB R2016a。

基于文献 [15] 的实验证实了 $l_0$ 优化问题可连续的有效性，我们将算法中嵌入连续化过程。具体方法就是在考虑求解问题(2)前，先求解如下一系列问题

$x_{{\mu }_k}^{*}:=\underset{x\in R_{+}^n}{\arg \min }\left\{{\varphi }_{{\mu }_k}\left(x\right)=f\left(x\right)+{\mu }_k{\Vert x\Vert }_0\right\}.$

其中 ${\mu }_{\text{0}}>{\mu }_1>\cdots >{\mu }_N$ 且N是一个正整数。我们用问题的解 $x_{{\mu }_{k-1}}^{*}$ 作为下一个问题的初始迭代点。取 ${\mu }_0=\frac{1}{2}{\Vert A^{\text{T}}b\Vert }_{\infty }^2$， ${\mu }_N={10}^{-15}{\mu }_0$，区间 $\left[{\mu }_N,{\mu }_0\right]$ 经过取对数化处理后，再将其分成N个均匀分布的子区间。我们选择序列 $\left\{{\mu }_k\right\}$ 作为子区间的端点，并对 ${\mu }_k$ 做了一个处理，即令 ${\mu }_k=\max \left\{{\mu }_k,0.005\right\}$。算法中涉及的参数设计如下：

${\lambda }^0=1,M=10,{\lambda }_{\min }={10}^{-3},{\lambda }_{\max }=10,\delta ={10}^{-2},\eta =0.5$。首先用标准高斯分布的独立样本填充 $m\times n$ 的度量矩阵A，并且对这些行进行正交化。而真实解 $x^{*}$ 是一个具有有效集 $A^{*}$ 的T-稀疏信号，观测向量b定义为：

$b=A{x}^{*}+e$

其中e是根据均值为零且方差为10^-2的高斯分布绘制的。

在实验中，我们研究了参数N和L对算法的影响，设 $n=5\times {10}^3$ 和 $n={10}^4$， $m=round\left(0.2\ast n\right)$，并考虑了稀疏度的范围： $x^{*}$ 的非零分量的数量T在1到60都是可行的。算法的初始点是零向量，当满足 ${\Vert g{\left(x^k\right)}_{\tau \left(x^k\right)}\Vert }_{\infty }\le {10}^{-5}$ 或 ${\Vert g{\left(x^k\right)}_{\gamma \left(x^k\right)}\Vert }_{\infty }\le \text{0}\text{.05}$ 条件时，停止运算。下面将从平均迭代步数(Iter)、平均误差( ${\Vert x^k-x^{*}\Vert }_{\infty }$ )和运算时间(Time)三个方面来对ABB算法的性能进行评估。

Table 1. n = 5 × 10 ³ , L = 1 / 5 , the results of ABB algorithm to solve the problem

表1. n = 5 × 10 ³ , L = 1 / 5 时ABB算法求解问题的计算结果

Table 2. n = 5 × 10 ³ , L = 1 / 4 , the results of ABB algorithm to solve the problem

表2. $n=5\times {10}^3,L=1/4$ 时ABB算法求解问题的计算结果

Table 3. n = 5 × 10 ³ , L = 1 / 2 , the results of ABB algorithm to solve the problem

表3. $n=5\times {10}^3,L=1/2$ 时ABB算法求解问题的计算结果

Table 4. n = 5 × 10 ³ , L = 1 , the results of ABB algorithm to solve the problem

表4. $n=5\times {10}^3,L=1$ 时ABB算法求解问题的计算结果

由上述表1~4可以看出，当问题的规模取 $n=5\times {10}^3$ 时， $T=31:60,N=15$，ABB算法产生了较大的平均误差。在 $T=1:30,N=10$ 时，ABB算法产生的平均误差较小。而L的取值变化对于ABB算法产生的平均误差影响很小。

对于上述表格中计算的数据，我们将绘图分析在 $T=1:30$ 时，L和N的取值变化对于ABB算法求解问题的平均迭代次数和运算时间产生的影响。

Figure 1. The influence of the change of N value on the average number of iterations at different L levels

图1. 不同L水平下，N取值变化对平均迭代次数的影响

Figure 2. The influence of the change of N value on operation time at different L levels

图2. 不同L水平下，N取值变化对运算时间的影响

Figure 3. The influence of the change of L value on the average number of iterations at different N levels

图3. 不同N水平下，L取值变化对平均迭代次数的影响

Figure 4. The influence of the change of L value on operation time at different N levels

图4. 不同N水平下，L取值变化对运算时间的影响

由图1和图2可以明显看出，当L的值一定时，随着N的取值变大，迭代步数和运算时间都会增加。同样的，由图3和图4可以明显看出，当N的取值一定时，可以看出 $L=1/4$ 时，平均迭代步数和运算时间都是最小的。由此可知，当 $L=1/4,N=10$ 时，可以得到最优的实验结果。

Table 5. n = 10 4 , L = 1 / 5 , the results of ABB algorithm to solve the problem

表5. $n={10}^4,L=1/5$ 时ABB算法求解问题的计算结果

Table 6. n = 10 4 , L = 1 / 4 , the results of ABB algorithm to solve the problem

表6. $n={10}^4,L=1/4$ 时ABB算法求解问题的计算结果

Table 7. n = 10 4 , L = 1 / 2 , the results of ABB algorithm to solve the problem

表7. $n={10}^4,L=1/2$ 时ABB算法求解问题的计算结果

Table 8. n = 10 4 , L = 1 , the results of ABB algorithm to solve the problem

表8. $n={10}^4,L=1$ 时ABB算法求解问题的计算结果

由表5~8中的结果可知，当问题规模达到 $n={10}^4$ 时，我们依旧可以得到相同的结论。当L的值一定时，随着N的取值变大，迭代步数和运算时间都会增加。当N的取值一定时，可以看出 $L=1/4$ 时，平均迭代步数和运算时间都是最小的。由此可知，当 $L=1/4,N=10$ 时，可以得到最优的实验结果。综上可知，ABB算法可以快速有效地求解欠定线性方程组的稀疏非负解。

6. 小结

本文主要讨论了求解欠定线性方程组的稀疏非负解问题。首先将线性约束条件惩罚到目标函数上，给出了 $l_0$ 正则优化模型的严格L-稳定点定义。随之提出有效集方法识别严格L-稳定点邻域内的零分量，基于这种快速精准的识别技术，设计了有效集Barzilar-Borwein算法，并且证明了由该算法产生的序列的任意聚点都是强稳定点。最后进行了数据实验，证明了该算法可以快速有效地求解欠定线性方程组的稀疏非负解。

NOTES

^*第一作者。

^#通讯作者。