Alternating Direction Method of Multipliers Based on $\ell_{2,0}$ -norm for Multiple Measurement Vector Problem

Zekun Liu, Siwei Yu
School of Mathematics, Harbin Institute of Technology Corresponding author: [email protected]

Abstract

The multiple measurement vector (MMV) problem is an extension of the single measurement vector (SMV) problem, and it has many applications. Nowadays, most studies of the MMV problem are based on the $\ell_{2,1}$ -norm relaxation, which will fail in recovery under some adverse conditions. We propose an alternating direction method of multipliers (ADMM)-based optimization algorithm to achieve a larger undersampling rate for the MMV problem. The key innovation is the introduction of an $\ell_{2,0}$ -norm sparsity constraint to describe the joint-sparsity of the MMV problem; this differs from the $\ell_{2,1}$ -norm constraint that has been widely used in previous studies. To illustrate the advantages of the $\ell_{2,0}$ -norm, we first prove the equivalence of the sparsity of the row support set of a matrix and its $\ell_{2,0}$ -norm. Then, the MMV problem based on the $\ell_{2,0}$ -norm is proposed. Next, we give our algorithm called MMV-ADMM- $\ell_{2,0}$ by applying ADMM to the reformulated problem. Moreover, based on the Kurdyka-Lojasiewicz property of objective functions, we prove that the iteration generated by the proposed algorithm globally converges to the optimal solution of the MMV problem. Finally, the performance of the proposed algorithm and comparisons with other algorithms under different conditions are studied with simulated examples. The results show that the proposed algorithm can solve a larger range of MMV problems even under adverse conditions.

1 Introduction

The multiple measurement vector (MMV) problem [1] has received considerable attention in the field of compressed sensing [2, 3]. In the single measurement vector (SMV) problem, the purpose is to recover a vector $x\in\mathbb{R}^{n}$ from $y=Ax$ , where $A\in\mathbb{R}^{m\times n}(m<n)$ and $y\in\mathbb{R}^{m}$ are given. Because $m<n$ , the problem is underdetermined. With the prior information that $x$ is sparse enough [4], we can solve

		$\displaystyle\min_{x\in\mathbb{R}^{n}}\left\\|x\right\\|_{0}$		(1.1)
		$\displaystyle s.t.\enspace y=Ax,$		(1.1)

to obtain the unique solution, where the $\ell_{0}$ -norm of a vector is the number of its nonzero elements. However, problem (1.1) is NP-hard and is difficult to solve directly. A common approach to overcome this drawback is to relax the $\ell_{0}$ -norm to the $\ell_{1}$ -norm; this is known as the basis pursuit [5]:

		$\displaystyle\min_{x\in\mathbb{R}^{n}}\left\\|x\right\\|_{1}$		(1.2)
		$\displaystyle s.t.\enspace y=Ax.$		(1.2)

Problem (1.2) is convex, implying that it can be solved much more efficiently. By solving problem (1.2) we can obtain the same solution as that for problem (1.1) under suitable conditions [3].

The MMV problem is an extension of the SMV problem. It describes the case in which multiple signals are jointly measured by the same sensing matrix from signals that are jointly sparse, that is

y^{j}=Ax^{j},j=1,2,\cdots,J,

where all of $x^{j}$ have the same locations of nonzero elements. The MMV problem occurs in many applications including hyperspectral imagery [6], computation of sparse solutions to linear inverse problems [7], neuromagnetic imaging [8], source localization [9], and equalization of sparse communication channels [10]. Given $Y=(y^{1},y^{2},\cdots,y^{J})\in\mathbb{R}^{m\times J}$ , the MMV problem is to recover $X=(x^{1},x^{2},\cdots,x^{J})$ from $Y=AX$ where only a few rows of $X$ have nonzero elements. For this purpose, the most widely used approach is $\ell_{2,1}$ -norm minimization [7, 11, 12]:

		$\displaystyle\min_{X\in\mathbb{R}^{n\times J}}\left\\|X\right\\|_{2,1}$		(1.3)
		$\displaystyle\quad s.t.\enspace Y=AX,$		(1.3)

and the definition of the generalized $\ell_{p,q}$ -norm of a matrix is

\left\|X\right\|_{p,q}=(\sum_{i=1}^{n}\left\|X_{i}\right\|_{p}^{q})^{\frac{1}{q}},

where $X_{i}$ is the $i$ th row of $X$ . Owing to the benefit of the additional joint-sparsity property, the recovery performance in the accuracy and total time spent for the MMV problem are better than those of the SMV problem [11, 12].

In recent decades, several recovery algorithms for the MMV problem have been proposed [11, 7, 13, 14]. [13] employed the greedy pursuit strategy to recover the signals because the joint-sparsity recovery problem is NP-hard. As a similar convex relaxation technique for the SMV problem (1.2), [7, 11, 12] used $\ell_{2,1}$ -norm minimization to solve the MMV problem. For hyperspectral imagery applications, [6] applied the alternating direction method of multipliers (ADMM) to solve problem (1.3). Most theoretical results on the relaxation of the SMV problem can be extended to the MMV problem [11]. However, as (1.2) is a convex relaxation of (1.1) that yields the same solution only under suitable conditions [3], (1.3) is only a convex relaxation of the $\ell_{2,0}$ -norm minimization problem for MMV. (1.3) cannot provide an accurate solution under certain unfavorable conditions, whereas the $\ell_{2,0}$ -norm minimization problem can. Therefore, some drawbacks limit the use of previous algorithms.

In this study, instead of considering the widely used convex relaxation $\ell_{2,1}$ -norm minimization problem, we directly study the original $\ell_{2,0}$ -norm minimization in the MMV problem. We first note the equivalence of the joint-sparsity property and the $\ell_{2,0}$ -norm in the MMV problem. Then, we reformulate the MMV problem via the sparsity constraint in [11]. Next, we propose our algorithm called MMV-ADMM- $\ell_{2,0}$ by applying ADMM [15] to the reformulated problem. Theoretical analysis shows that the proposed algorithm is globally convergent to the unique solution of the MMV problem. Compared with existing algorithms, the proposed algorithm can solve the MMV problem when the signals are not very sparse or the number of sensors is small, and it can achieve a larger undersampling rate.

The rest of this paper is organized as follows: In Section 2, we present some definitions used for the subsequent analysis. In Section 3, we propose the optimization problem by reformulating the MMV problem. In Section 4, we present the proposed MMV-ADMM- $\ell_{2,0}$ algorithm and discuss the subproblems in detail. In Section 5, we present the convergence results obtained with the proposed algorithm. In Section 6, we describe numerical experiments performed to verify the validity of the proposed algorithm and compare it with other algorithms. In Section 7, we present the conclusions of this study.

2 Preliminaries

We list some of the notations and definitions used for the subsequent analysis. For any vector $x=\left(x_{1},x_{2},\cdots,x_{N}\right)^{T}\in\mathbb{R}^{N}$ , its sparse support set is defined as

Supp(x)=\left\{i|x_{i}\neq 0\right\}\subseteq\left\{1,2,\cdots,N\right\}.

Note that the assumption of the MMV problem is that the set of all vectors $\left\{x^{j}\right\}_{j=1}^{J}\in\mathbb{R}^{N}$ has the same sparse support set, implying that

Supp(x^{1})=Supp(x^{2})=\cdots=Supp(x^{J}).

Definition 2.1.

For a set of vectors $\left\{x^{j}\right\}_{j=1}^{J}\in\mathbb{R}^{N}$ , the common sparse support set is defined as

Supp(\left\{x^{j}\right\}_{j=1}^{J})=\bigcup_{j=1}^{J}Supp(x^{j}).

In fact, although it is almost impossible to satisfy the above assumption with real datasets, the MMV problem works well when real data have sufficient joint-sparsity [11].

The following few definitions are used for the subsequent analysis.

Definition 2.2 ([16]).

Let $f$ be a generalized real function.

(i)

For a nonempty set $\mathcal{X}$ , we call $f$ proper to $\mathcal{X}$ if there exist $x\in\mathcal{X}$ such that $f(x)<+\infty$ and for any $x\in\mathcal{X},f(x)>-\infty$ .
(ii)

$f$ is a lower semicontinuous function if for any $x\in\mathbb{R}^{m\times n},\liminf_{y\to x}f(y)\geq f(x)$ .
(iii)

$f$ is a closed function if its epigraph

$epif=\left\{(x,t)\in\mathbb{R}^{m\times n}\times\mathbb{R}|f(x)\leq t\right\}$

is a closed set.

For a generalized real function $f$ , the lower semicontinuity of $f$ is equivalent to it being a closed function.

Definition 2.3 ([17]).

For a proper lower semicontinuous function $f$ , its critical point is one that satisfies $0\in\partial f(x)$ .

In keeping with the first-order optimal condition, if $x$ is a minimizer of $f$ , it must be a critical point of $f$ first.

Definition 2.4 ([18]).

$S\subseteq\mathbb{R}^{m\times n}$ is a semi-algebraic set if there exist polynomial functions $f_{ij},g_{ij}:\mathbb{R}^{m\times n}\to\mathbb{R}$ with $1\leq i\leq p,1\leq j\leq q$ such that

S=\bigcup_{i=1}^{p}\bigcap_{j=1}^{q}\left\{x\in\mathbb{R}^{m\times n}:f_{ij}(x)=0,g_{ij}(x)>0\right\}.

The proper function $f$ is semi-algebraic if its graph

graphf=\left\{(x,t)\in\mathbb{R}^{m\times n}\times\mathbb{R}:f(x)=t\right\}

is a semi-algebraic set of $\mathbb{R}^{m\times n}\times\mathbb{R}$ .

Semi-algebra has a few useful properties:

(i)

Semi-algebraic sets are closed under finite unions and Cartesian products.
(ii)

Many common functions are actually semi-algebraic, including real polynomials, compositions of semi-algebraic functions, indicator functions of semi-algebraic sets, and $\pi(A)$ where $A\subseteq\mathbb{R}^{m\times n}\times\mathbb{R}$ is a semi-algebraic set and $\pi:\mathbb{R}^{m\times n}\times\mathbb{R}\to\mathbb{R}^{m\times n}$ is the projection operator.

The definitions of Gâteaux differentiable, subdifferential, Karush-Kuhn-Tucker (KKT) conditions, and Kurdyka-Lojasiewicz (KL) property can be referred elsewhere [19, 17, 16, 18, 20].

The following proposition states the connection between semi-algebraic functions and KL functions.

Proposition 2.1 ([20]).

If the proper lower semicontinuous function $f$ is semi-algebraic, it must be a KL function.

3 Problem formulation

First, we state a proposition that is useful to describe the MMV problem in mathematical language.

Proposition 3.1.

For a set of vectors $\left\{s^{j}\right\}_{j=1}^{J}\in\mathbb{R}^{N}$ , denote $S=(s^{1},s^{2},\cdots,s^{J})$ . Then, $\left\{s^{j}\right\}_{j=1}^{J}$ has the common sparse support set with sparsity $k\Longleftrightarrow\left\|S\right\|_{2,0}=k$ .

Proof.

Denote the common sparse support set of $\left\{s^{j}\right\}_{j=1}^{J}$ as $Supp$ . $S$ can be expressed by its row vector as $S=(\alpha_{1},\alpha_{2},\cdots,\alpha_{N})^{T}$ . Then, we can conclude that

		$\displaystyle\left\\|S\right\\|_{2,0}=\left\\|(\left\\|\alpha_{1}\right\\|_{2},\left\\|\alpha_{2}\right\\|_{2},\cdots,\left\\|\alpha_{N}\right\\|_{2})^{T}\right\\|_{0}=k$
	$\displaystyle\Longleftrightarrow$	$\displaystyle\begin{cases}\left\\|\alpha_{j}\right\\|_{2}\neq 0\quad\text{for}\enspace j\in Supp\\ \left\\|\alpha_{j}\right\\|_{2}=0\quad\text{for}\enspace j\notin Supp\end{cases}\quad and\enspace\left\|Supp\right\|=k\quad\text{(definition of $\ell_{0}$-norm)}$
	$\displaystyle\Longleftrightarrow$	$\displaystyle\begin{cases}\alpha_{j}\neq 0\quad\text{for}\enspace j\in Supp\\ \alpha_{j}=0\quad\text{for}\enspace j\notin Supp\end{cases}\quad\text{(positivity of $\ell_{2}$-norm)}$
	$\displaystyle\Longleftrightarrow$	$\displaystyle Supp\enspace\text{is actually the common sparse support set with}\enspace\left\|Supp\right\|=k.$

This completes the whole proof. $\square$

Next, we convert the MMV problem to an optimization problem by using Proposition 3.1.

Assuming there are $J$ sensors, the sparse vectors after sparse representation are $s_{1},s_{2},\cdots,s_{J}$ , where $s_{j}\in\mathbb{R}^{N}$ for all $j\in\left\{1,2,\cdots,J\right\}$ ; we denote $S=(s_{1},s_{2},\cdots,s_{J})\in\mathbb{R}^{N\times J}$ . Assuming the sensing matrix $\Phi\in\mathbb{R}^{M\times N}(M<N)$ , the measurement vectors are $y_{1},y_{2},\cdots,y_{J}$ , where $y_{j}=\Phi s_{j}\in\mathbb{R}^{M}$ for all $j\in\left\{1,2,\cdots,J\right\}$ ; we denote $Y=(y_{1},y_{2},\cdots,y_{J})\in\mathbb{R}^{M\times J}$ . By Proposition 3.1, minimizing the joint-sparsity of $\left\{s_{j}\right\}_{j=1}^{J}$ is equivalent to minimizing $\left\|S\right\|_{2,0}$ . Therefore, the MMV problem can be described as

		$\displaystyle\min_{S\in\mathbb{R}^{N\times J}}\left\\|S\right\\|_{2,0}$		(3.1)
		$\displaystyle\quad s.t.\quad Y=\Phi S,$		(3.1)

The following theorem offers the prerequisite to make a further conversion of problem (3.1).

Theorem 3.1 ([11]).

For $\Phi\in\mathbb{R}^{M\times N}$ and $S\in\mathbb{R}^{N\times J}$ , if $Y=\Phi S$ and

\left\|S\right\|_{2,0}<\frac{Spark(\Phi)+Rank(Y)-1}{2},

(3.2)

then problem (3.1) will only possess the unique solution $S$ .

Based on Theorem 3.1, we do not require all vectors to share the same sparse support set in the MMV problem. Instead, we just require that $S$ satisfies (3.2).

To consider the odevity of $Spark(\Phi)+Rank(Y)$ and the measurement error between $Y$ and the real $\Phi S$ , set

s=\left\lfloor\frac{Spark(\Phi)+Rank(Y)-2}{2}\right\rfloor,

(3.3)

where $\left\lfloor\cdot\right\rfloor$ is the round down operator. Problem (3.1) can be converted to

		$\displaystyle\min_{S\in\mathbb{R}^{N\times J}}\left\\|Y-\Phi S\right\\|_{F}^{2}$		(3.4)
		$\displaystyle\quad s.t.\quad\left\\|S\right\\|_{2,0}\leq s.$		(3.4)

It is a nonconvex constrained optimization problem. The indicator function

\mathcal{I}_{\mathcal{M}}(X)=\begin{cases}0,&\text{ {if} }X\in\mathcal{M}\\ +\infty,&\text{ {if} }X\notin\mathcal{M}\end{cases}

is introduced to the set

\mathcal{M}=\left\{X\in\mathbb{R}^{N\times J}:\left\|X\right\|_{2,0}\leq s\right\}

to move the nonconvex constraint to the objective function, and matrix $B\in\mathbb{R}^{N\times J}$ is introduced as an auxiliary variable of $S$ to reformulate problem (3.4) as

		$\displaystyle\min_{B,S\in\mathbb{R}^{N\times J}}\left\\|Y-\Phi S\right\\|_{F}^{2}+\mathcal{I}_{\mathcal{M}}(B)$		(3.5)
		$\displaystyle\quad s.t.\quad\quad B-S=0.$		(3.5)

Now we have obtained the final problem formulation (3.5) to describe the MMV problem.

4 Algorithm

Problem (3.5) is a two block nonconvex optimization problem with a linear equality constraint; it can be solved using ADMM [15].

The augmented Lagrangian function associated with problem (3.5) is defined as

\mathcal{L}_{\rho}(B,S,L)=\left\|Y-\Phi S\right\|_{F}^{2}+\mathcal{I}_{\mathcal{M}}(B)+\left\langle L,B-S\right\rangle+\frac{\rho}{2}\left\|B-S\right\|_{F}^{2},

(4.1)

where $L\in\mathbb{R}^{N\times J}$ is the Lagrangian multiplier of the equation constraint, $\left\langle\cdot,\cdot\right\rangle$ is the inner product of matrices, and $\rho>0$ is the penalty parameter.

By using the ADMM scheme, we can obtain an approximate solution of problem (3.5) by minimizing the variables one by one with others fixed, as follows:

\displaystyle\left\{\begin{matrix}B^{k+1}=&\mathop{\arg\min}\limits_{B\in\mathbb{R}^{N\times J}}\mathcal{L}_{\rho}(B,S^{k},L^{k}),\\ S^{k+1}=&\mathop{\arg\min}\limits_{S\in\mathbb{R}^{N\times J}}\mathcal{L}_{\rho}(B^{k+1},S,L^{k}),\\ L^{k+1}=&L^{k}+\rho(B^{k+1}-S^{k+1}).\end{matrix}\right.

(4.2)

Next, the subproblems in Equation (4.2) are solved one by one.

4.1 Update B

Fix $S$ and $L$ ,

\begin{split}B^{k+1}&=\mathop{\arg\min}\limits_{B\in\mathbb{R}^{N\times J}}\mathcal{L}_{\rho}(B,S^{k},L^{k})\\ &=\mathop{\arg\min}\limits_{B\in\mathbb{R}^{N\times J}}\left\{\mathcal{I}_{\mathcal{M}}(B)+\left\langle L^{k},B-S^{k}\right\rangle+\frac{\rho}{2}\left\|B-S^{k}\right\|_{F}^{2}\right\}\\ &=\mathop{\arg\min}\limits_{B\in\mathcal{M}}\left\|B-S^{k}+\frac{L^{k}}{\rho}\right\|_{F}^{2}\\ &=\mathcal{P}_{\mathcal{M}}(S^{k}-\frac{L^{k}}{\rho}),\end{split}

(4.3)

where $\mathcal{P}_{\mathcal{M}}(\cdot)$ is the projection operator on set $\mathcal{M}=\left\{X\in\mathbb{R}^{N\times J}:\left\|X\right\|_{2,0}\leq s\right\}$ .

Actually, $\mathcal{P}_{\mathcal{M}}(\cdot)$ is the hard thresholding operator [21]. If $\left\|S^{k}-\frac{L^{k}}{\rho}\right\|_{2,0}\leq s$ , then $B^{k+1}=S^{k}-\frac{L^{k}}{\rho}$ ; otherwise, truncate $S^{k}-\frac{L^{k}}{\rho}$ with the rows whose $\ell_{2}$ -norm is the top $s$ large preserved.

4.2 Update S

Fix $B$ and $L$ ,

\begin{split}S^{k+1}&=\mathop{\arg\min}\limits_{S\in\mathbb{R}^{N\times J}}\mathcal{L}_{\rho}(B^{k+1},S,L^{k})\\ &=\mathop{\arg\min}\limits_{S\in\mathbb{R}^{N\times J}}\left\{\left\|Y-\Phi S\right\|_{F}^{2}+\left\langle L^{k},B^{k+1}-S\right\rangle+\frac{\rho}{2}\left\|B^{k+1}-S\right\|_{F}^{2}\right\}\\ &=\mathop{\arg\min}\limits_{S\in\mathbb{R}^{N\times J}}\left\{\left\|Y-\Phi S\right\|_{F}^{2}+\frac{\rho}{2}\left\|B^{k+1}-S+\frac{L^{k}}{\rho}\right\|_{F}^{2}\right\},\end{split}

(4.4)

Denote $f(X)=\left\|Y-\Phi X\right\|_{F}^{2}+\frac{\rho}{2}\left\|B^{k+1}-X+\frac{L^{k}}{\rho}\right\|_{F}^{2}$ . The first-order optimal condition for unconstrained differentiable optimization problem gives

S^{k+1}=\mathop{\arg\min}\limits_{S\in\mathbb{R}^{N\times J}}f(S)\Longleftrightarrow\nabla f(S^{k+1})=0.

It is easy to calculate

\nabla f(S^{k+1})=(2\Phi^{T}\Phi+\rho I)S^{k+1}-2\Phi^{T}Y-\rho B^{k+1}-L^{k}.

Therefore, by solving the linear equation $\nabla f(S^{k+1})=0$ , we obtain

S^{k+1}=(2\Phi^{T}\Phi+\rho I)^{-1}(2\Phi^{T}Y+\rho B^{k+1}+L^{k}).

(4.5)

4.3 Update L

Fix $B$ and $S$ . The Lagrangian multiplier $L$ can be updated as

L^{k+1}=L^{k}+\rho(B^{k+1}-S^{k+1}).

(4.6)

4.4 Convergence criterion

The convergence criterion of the proposed algorithm is given below, and the reasons will be discussed in the following section.

Assuming $\left\{(B^{k},S^{k},L^{k})\right\}$ is the sequence generated using the ADMM procedure (4.2), the convergence criterion is given as

\begin{cases}\left\|B^{k}-S^{k}\right\|_{F}\to 0,\\ \left\|S^{k+1}-S^{k}\right\|_{F}\to 0,\\ \left\|L^{k}\right\|_{F}\to 0.\end{cases}

(4.7)

We call our algorithm for solving problem (3.5) the MMV-ADMM- $\ell_{2,0}$ ; it is summarized in Algorithm 1. Section 6 describes the parameter initialization.

Input: The measurement data matrix

Y

, the sensing matrix

\Phi

;

Output: The reconstructed sparse matrix

\hat{S}

;

1 initialization:

S^{0},L^{0},s,\rho

, and let

k=0

;

2do

3 Update

B

by:

B^{k+1}=\mathcal{P}_{\mathcal{M}}(S^{k}-\frac{L^{k}}{\rho})

;

4 Update

S

by:

S^{k+1}=(2\Phi^{T}\Phi+\rho I)^{-1}(2\Phi^{T}Y+\rho B^{k+1}+L^{k})

;

5 Update the Lagrangian multiplier

L

L^{k+1}=L^{k}+\rho(B^{k+1}-S^{k+1})

;

6 Update

k:k=k+1

;

7while not satisfy the convergence criterion;

return

\hat{S}=B^{k}

Algorithm 1 MMV-ADMM-

\ell_{2,0}

Notably, $2\Phi^{T}\Phi+\rho I$ is unchanged in the iteration; therefore, its inverse needs to be calculated only once. Algorithm 1 can be accelerated in an obvious way. When updating $S$ by (4.5), we need to calculate the inverse of an $N\times N$ matrix; the computational time required for this purpose will be high when $N$ is large. Because $M<N$ in compressed sensing, we can use the Sherman-Morrison-Woodbury (SMW) formula [22, 23] to simplify the calculation of the inverse. Specifically, based on the SMW-formula, we have

(2\Phi^{T}\Phi+\rho I)^{-1}=\frac{I}{\rho}-\frac{2\Phi^{T}(I+\frac{2\Phi\Phi^{T}}{\rho})^{-1}\Phi}{\rho^{2}}.

(4.8)

It converts the calculation of an $N\times N$ matrix inverse to the inverse of an $M\times M$ matrix. Because $M$ is much smaller than $N$ in compressed sensing, this will reduce the computational time. We call our algorithm MMV-ADMM- $\ell_{2,0}$ -SMW when using (4.8) to update $S$ .

5 Convergence analysis

Algorithm 1 is a two-block ADMM for a nonconvex problem; its convergence analysis remains an open problem. However, owing to the KL property of objective functions in problem (3.5), Algorithm 1 has global convergence.

First, we prove some properties of objective functions in problem (3.5).

Proposition 5.1.

The $\ell_{2,0}$ -norm of a matrix is proper and lower semicontinuous.

Proof.

First, we prove that the $\ell_{0}$ -norm of a vector is lower semicontinuous. Denote

	$\displaystyle f:$	$\displaystyle\mathbb{R}^{N}\to\mathbb{R}$
		$\displaystyle x\mapsto\left\\|x\right\\|_{0}.$

For any $x=(x_{1},x_{2},\cdots,x_{N})^{T}\in\mathbb{R}^{N}$ , assume

f(x)=\left\|x\right\|_{0}=s,

Supp(x)=\left\{l_{1},l_{2},\cdots,l_{s}\right\}\subseteq\left\{1,2,\cdots,N\right\}.

When $y^{k}=(y_{1}^{k},y_{2}^{k},\cdots,y_{N}^{k})^{T}\overset{k\to\infty}{\longrightarrow}x$ , for

\varepsilon=\frac{\min_{j=1,2,\cdots,s}\left|x_{l_{j}}\right|}{2}>0,

there exist $K\in\mathbb{N}$ such that

\left|y_{i}^{k}-x_{i}\right|<\varepsilon,i=1,2,\cdots,N\text{ for all }k>K.

Therefore, when $k>K$ ,

y_{l_{j}}^{k}\neq 0,j=1,2,\cdots,s,

implying $f(y^{k})=\left\|y^{k}\right\|_{0}\geq s.$

According to the sign-preserving property of limit, we have $\liminf_{k\to\infty}f(y^{k})\geq s=f(x)$ ; therefore, $f$ is lower semicontinuous.

For any $S\in\mathbb{R}^{N\times J}$ , $S$ can be expressed by its row vector as $S=(\alpha_{1},\alpha_{2},\cdots,\alpha_{N})^{T}$ . Denote

	$\displaystyle g:$	$\displaystyle\mathbb{R}^{N\times J}\to\mathbb{R}^{N}$
		$\displaystyle S\mapsto(\left\\|\alpha_{1}\right\\|_{2},\left\\|\alpha_{2}\right\\|_{2},\cdots,\left\\|\alpha_{N}\right\\|_{2})^{T}.$

Then, the $\ell_{2,0}$ -norm of a matrix

	$\displaystyle h:$	$\displaystyle\mathbb{R}^{N\times J}\to\mathbb{R}$
		$\displaystyle S\mapsto\left\\|S\right\\|_{2,0}$

is the composition of $f$ and $g$ , that is, $h=f\circ g$ .

It is apparent that $g$ is continuous. Upon combining $f$ , a lower semicontinuous function, we can conclude that for any $S\in\mathbb{R}^{N\times J}$ , $S^{k}\to S$ ,

\liminf_{k\to\infty}h(S^{k})=\liminf_{k\to\infty}f(g(S^{k}))\geq f(\lim_{k\to\infty}g(S^{k}))=f(g(S))=h(S).

Therefore, $h$ is lower semicontinuous.

Apparently, the $\ell_{2,0}$ -norm is proper. Therefore, it is proper and lower semicontinuous. $\square$

Based on the equivalence of the semicontinuous function and closed function, Corollary 5.1 holds immediately.

Corollary 5.1.

The $\ell_{2,0}$ -norm of a matrix is a closed function.

Corollary 5.2.

The indicator function

\mathcal{I}_{\mathcal{M}}(X)=\begin{cases}0,&\text{ {if} }X\in\mathcal{M}\\ +\infty,&\text{ {if} }X\notin\mathcal{M}\end{cases}

to the set $\mathcal{M}=\left\{X\in\mathbb{R}^{N\times J}:\left\|X\right\|_{2,0}\leq s\right\}$ is also proper and lower semicontinuous.

Proof.

We conclude that $\mathcal{I}_{\mathcal{M}}(X)$ is a closed function. In fact, for any sequence $(A_{k},t_{k})\in epi\mathcal{I}_{\mathcal{M}}$ , we have $t_{k}\geq\mathcal{I}_{\mathcal{M}}(A_{k})$ . According to Definition 2.2, we only need to prove if $(A_{k},t_{k})\to(A,t)$ ; then, $t\geq\mathcal{I}_{\mathcal{M}}(A)$ .

Because $t_{k}\to t$ , $k$ is large enough in only two cases:

If $\left\|A_{k}\right\|_{2,0}>s$ , then $t_{k}=\mathcal{I}_{\mathcal{M}}(A_{k})=+\infty$ . From $t_{k}\to t$ , we have $t=+\infty$ ; therefore, $t\geq\mathcal{I}_{\mathcal{M}}(A)$ .

If $\left\|A_{k}\right\|_{2,0}\leq s$ , then $t_{k}\geq\mathcal{I}_{\mathcal{M}}(A_{k})=0$ . From $t_{k}\to t$ , we have $t\geq 0$ . Through Corollary 5.1, we know that the $\ell_{2,0}$ -norm of a matrix is a closed function. Therefore,

$\left.\begin{aligned} &(A_{k},s)\in epi\left\|\cdot\right\|_{2,0}\\ &(A_{k},s)\to(A,s)\end{aligned}\right\}\Longrightarrow(A,s)\in epi\left\|\cdot\right\|_{2,0},$

Implying that $\left\|A\right\|_{2,0}\leq s$ . Therefore, $t\geq 0=\mathcal{I}_{\mathcal{M}}(A)$ .

Further, $t\geq\mathcal{I}_{\mathcal{M}}(A)$ . This proves that $\mathcal{I}_{\mathcal{M}}(X)$ is closed. It is apparent that $\mathcal{I}_{\mathcal{M}}(X)$ is proper. Therefore, $\mathcal{I}_{\mathcal{M}}(X)$ is also proper and lower semicontinuous. $\square$

Proposition 5.2.

The $\ell_{2,0}$ -norm of a matrix is a KL function.

Proof.

[20] proved that both $\left\|\cdot\right\|_{0}$ and $\left\|\cdot\right\|_{p}$ with $p\in\mathbb{Q}_{+}$ are semi-algebraic. Denote functions $f$ , $g$ , and $h$ in the same way as the definitions in Proposition 5.1. For any $S\in\mathbb{R}^{N\times J}$ , $S$ can be expressed by its row vector as $S=(\alpha_{1},\alpha_{2},\cdots,\alpha_{N})^{T}$ :

g(S)=(\left\|\alpha_{1}\right\|_{2},\left\|\alpha_{2}\right\|_{2},\cdots,\left\|\alpha_{N}\right\|_{2})^{T}=\prod_{i=1}^{N}\left\{\left\|\alpha_{i}\right\|_{2}\right\}

is the Cartesian product of the $\ell_{2}$ -norm. Because the semi-algebra is stable under a Cartesian product, $g$ is also semi-algebraic. Therefore, the composition $h=f\circ g$ is semi-algebraic.

Proposition 5.1 proves that $h$ is proper lower semicontinuous. Therefore, $h$ is a KL function according to Proposition 2.1. $\square$

Corollary 5.3.

The indicator function

\mathcal{I}_{\mathcal{M}}(X)=\begin{cases}0,&\text{ {if} }X\in\mathcal{M}\\ +\infty,&\text{ {if} }X\notin\mathcal{M}\end{cases}

to the set $\mathcal{M}=\left\{X\in\mathbb{R}^{N\times J}:\left\|X\right\|_{2,0}\leq s\right\}$ is also a KL function.

Proof.

First, we prove that $\mathcal{M}$ is semi-algebraic. Because $\left\|\cdot\right\|_{2,0}\in\mathbb{N}$ , $\mathcal{M}$ can be expressed as

\mathcal{M}=\left\{X\in\mathbb{R}^{N\times J}:\left\|X\right\|_{2,0}\leq s\right\}=\bigcup_{t=0}^{s}\left\{X\in\mathbb{R}^{N\times J}:\left\|X\right\|_{2,0}=t\right\}.

Proposition 5.2 proves that the $\ell_{2,0}$ -norm is a semi-algebraic function. Therefore, we conclude through Definition 2.4 that

graph\left\|\cdot\right\|_{2,0}=\left\{(X,t):\left\|X\right\|_{2,0}=t\right\}

is a semi-algebraic set.

Denote the projection operator

	$\displaystyle\pi:$	$\displaystyle\mathbb{R}^{N\times J}\times\mathbb{R}\to\mathbb{R}^{N\times J}$
		$\displaystyle(X,t)\mapsto X.$

Then,

\mathcal{M}=\bigcup_{t=0}^{s}\pi(graph\left\|\cdot\right\|_{2,0}).

Because the semi-algebra is stable under finite union and the projection operator $\pi$ , $\mathcal{M}$ is semi-algebraic.

$\mathcal{I}_{\mathcal{M}}(X)$ is the indicator function to $\mathcal{M}$ ; it is also semi-algebraic. Corollary 5.2 proves that it is also proper lower semicontinuous. Therefore, by Proposition 2.1, $\mathcal{I}_{\mathcal{M}}(X)$ is a KL function. $\square$

Next, we prove that Algorithm 1 has global convergence based on [24].

In [24] the authors discuss the following optimization problem:

		$\displaystyle\min_{x,y}f(x)+g(y)$		(5.1)
		$\displaystyle s.t.\quad Ax+y=b,$		(5.1)

where $f$ is proper and lower semicontinuous, and $g$ is Gradient- $L$ -Lipschitz continuous.

The authors make the following assumptions:

(i)

The two minimization subproblems in (4.2) possess solutions, and
(ii)

The penalty parameter $\rho>2L$ ,
(iii)

$A^{T}A\succeq\mu I$ for some $\mu>0$ .

Set $A=-I$ , $x=B$ , $y=S$ , $b=0$ , $f(B)=\mathcal{I}_{\mathcal{M}}(B)$ , and $g(S)=\left\|Y-\Phi S\right\|_{F}^{2}$ in (5.1); this gives problem (3.5). It is apparent that problem (3.5) satisfies the form of (5.1) and the above assumptions.

Next, we refer to the sufficient conditions in [24] to guarantee that the sequence $\left\{(B^{k},S^{k},L^{k})\right\}$ generated by Algorithm 1 is bounded.

Lemma 5.1 ([24]).

By applying the classic ADMM to (5.1), we obtain a sequence $\left\{w^{k}\right\}$ . Suppose that

\bar{g}:=\inf_{y}\left\{g(y)-\frac{1}{2L}\left\|\nabla g(y)\right\|^{2}\right\}>-\infty.

(5.2)

$\left\{w^{k}\right\}$ is bounded as long as one of the following two statements holds:

(i)

$\liminf_{\left\|x\right\|\to+\infty}f(x)=+\infty$ , and
(ii)

$\inf_{x}f(x)>-\infty$ and $\liminf_{\left\|y\right\|\to+\infty}g(y)=+\infty$ .

The main result in [24] is shown below.

Theorem 5.1 ([24]).

Suppose that Lemma 5.1 holds and that $f$ and $g$ are KL functions. Then, $\left\{w^{k}\right\}$ converges to the KKT point of (5.1).

As an application of the above theorem, Algorithm 1 has global convergence.

Theorem 5.2.

Algorithm 1 has global convergence, and the generated sequence $\left\{(B^{k},S^{k},L^{k})\right\}$ converges to the KKT point of problem (3.5).

Proof.

First, we prove that $\left\{(B^{k},S^{k},L^{k})\right\}$ is bounded.

For $g(S)=\left\|Y-\Phi S\right\|_{F}^{2}$ , we have $\nabla g(S)=2\Phi^{T}(\Phi S-Y)$ and $L=2\left\|\Phi\right\|_{F}^{2}$ . For any $S\in\mathbb{R}^{N\times J}$ ,

\left\|Y-\Phi S\right\|_{F}^{2}-\frac{1}{2L}\left\|2\Phi^{T}(\Phi S-Y)\right\|_{F}^{2}\geq\left\|Y-\Phi S\right\|_{F}^{2}-\frac{4\left\|\Phi\right\|_{F}^{2}}{2L}\left\|Y-\Phi S\right\|_{F}^{2}=0.

Therefore, (5.2) is satisfied. Moreover, it is not hard to see that (ii) in Lemma 5.1 holds. Therefore, $\left\{(B^{k},S^{k},L^{k})\right\}$ is bounded.

Corollary 5.3 proves that $f(B)=\mathcal{I}_{\mathcal{M}}(B)$ is a KL function. Because $g(S)=\left\|Y-\Phi S\right\|_{F}^{2}$ is a real polynomial, it is also a KL function.

All the conditions in Theorem 5.1 are satisfied; therefore, $\left\{(B^{k},S^{k},L^{k})\right\}$ globally converges to the KKT point of (3.5). $\square$

Finally, we state the theorem about settings for the convergence criterion (4.7).

Theorem 5.3.

Algorithm 1 is nearly convergent to the optimal solution of problem (3.5) if and only if its iteration sequence $\left\{(B^{k},S^{k},L^{k})\right\}$ satisfies the convergence criterion (4.7).

Proof.

In [25], the authors note that for the following problem

		$\displaystyle\min_{x}f(x)$		(5.3)
		$\displaystyle s.t.\enspace x\in\mathcal{X},\left\\|x\right\\|_{0}\leq k,$		(5.3)

where $\mathcal{X}$ denotes a convex polyhedron, the optimal solution must also satisfy the KKT conditions.

It is easy to see that problem (3.4) satisfies the form of (5.3). Because problem (3.5) is equivalent to problem (3.4), the optimal point of (3.5) is also its KKT point.

Therefore, if $\left\{(B^{k},S^{k},L^{k})\right\}$ converges to the optimal solution $(B^{\ast},S^{\ast},L^{\ast})$ , then $(B^{\ast},S^{\ast},L^{\ast})$ is a KKT point of (3.5). The Lagrangian function of (3.5) is

\mathcal{L}(B,S,L)=\left\|Y-\Phi S\right\|_{F}^{2}+\mathcal{I}_{\mathcal{M}}(B)+\left\langle L,B-S\right\rangle,

(5.4)

it satisfies the KKT conditions at $(B^{\ast},S^{\ast},L^{\ast})$ . Therefore,

\begin{cases}0=\nabla_{S}\mathcal{L}(B^{\ast},S^{\ast},L^{\ast})=2\Phi^{T}(\Phi S^{\ast}-Y)-L^{\ast},\\ 0\in\partial_{B}\mathcal{L}(B^{\ast},S^{\ast},L^{\ast})=\partial_{B}\mathcal{I}_{\mathcal{M}}(B^{\ast})+L^{\ast},\\ 0=B^{\ast}-S^{\ast}.\end{cases}

(5.5)

From $(B^{k},S^{k},L^{k})\to(B^{\ast},S^{\ast},L^{\ast})$ and the third formula of (5.5), we can conclude that $\left\|B^{k}-S^{k}\right\|_{F}\to 0$ .

Denote $f$ in the same way as in the definition in Subsection 4.2. When updating $S$ , we have

\nabla f(S^{k+1})=(2\Phi^{T}\Phi+\rho I)S^{k+1}-2\Phi^{T}Y-\rho B^{k+1}-L^{k}=0,

Let $k\to+\infty$ . We obtain $2\Phi^{T}(\Phi S^{\ast}-Y)-L^{\ast}=0$ , implying that the first formula of (5.5) is satisfied in Algorithm 1 without any further conditions.

From the rule of updating $B$ in (4.3), we know that $\mathcal{I}_{\mathcal{M}}(B^{k})=0$ is always satisfied during iterations in Algorithm 1. According to Corollary 5.2, $\mathcal{I}_{\mathcal{M}}(X)$ is proper and lower semicontinuous. Thus,

0\leq\mathcal{I}_{\mathcal{M}}(B^{\ast})\leq\liminf_{k\to+\infty}\mathcal{I}_{\mathcal{M}}(B^{k})=0.

Therefore, when limiting the space to all iteration points generated by Algorithm 1 and their accumulations, we have $\mathcal{I}_{\mathcal{M}}(X)\equiv 0$ . Thus, $\partial_{B}\mathcal{I}_{\mathcal{M}}(B^{\ast})=0$ . Therefore, the second formula of (5.5) holds if and only if $L^{\ast}=0$ , which is equivalent to $\left\|L^{k}\right\|_{F}\to 0$ .

Thus far we have proved that the KKT conditions (5.5) are equivalent to the first and third formulas in (4.7). From $S^{k}\to S^{\ast}$ , we have $\left\|S^{k+1}-S^{k}\right\|_{F}\to 0$ , which is the second formula in (4.7). Therefore, the optimal point $(B^{\ast},S^{\ast},L^{\ast})$ must satisfy the convergence criterion (4.7). This completes the necessity.

For sufficiency, the convergence criterion for a generalized constrained optimization problem [16] is that the distance between the adjacent iteration points tends to zero and the KKT conditions at the current iteration point are nearly satisfied. We have proved that the first and third formulas in (4.7) are equivalent to the KKT conditions. Therefore, the convergence criterion (4.7) ensures that the KKT conditions are satisfied. From combining the second formula in (4.7), which implies that the distance between the adjacent iteration points tends to zero, we know that the convergence criterion (4.7) guarantees $\left\{(B^{k},S^{k},L^{k})\right\}$ to converge to a KKT point of (3.5). Denote the KKT point as $(B^{\ast},S^{\ast},L^{\ast})$ ; we say that it is also the optimal point of (3.5).

In fact, at the KKT point $(B^{\ast},S^{\ast},L^{\ast})$ , we have proved that $L^{\ast}=0$ and (5.5) are satisfied. Therefore, from the first formula of (5.5), we can conclude that

\Phi^{T}(\Phi S^{\ast}-Y)=0.

(5.6)

The sensing matrix $\Phi\in\mathbb{R}^{M\times N}$ in the MMV problem is known to be row full rank. Thus, from (5.6), we know that $\Phi S^{\ast}-Y=0$ . We have proved $\mathcal{I}_{\mathcal{M}}(B^{\ast})=0$ , implying that $\left\|B^{\ast}\right\|_{2,0}\leq s$ . Because $B^{\ast}=S^{\ast}$ , $\left\|S^{\ast}\right\|_{2,0}\leq s$ too. Therefore, the KKT point $(B^{\ast},S^{\ast},L^{\ast})$ satisfies $Y=\Phi S^{\ast},\left\|S^{\ast}\right\|_{2,0}\leq s$ . From Theorem 3.1, $S^{\ast}$ is the unique optimal point of (3.1). According to the equivalence of (3.1) and (3.5), we can conclude that $(B^{\ast},S^{\ast},L^{\ast})$ is also the unique optimal point of (3.5).

Therefore, when the convergence criterion is satisfied, Algorithm 1 converges to the optimal solution of (3.5).

This completes the proof. $\square$

In fact, it is not hard to find that all of the discussions above are also correct when we replace the $\ell_{2,0}$ -norm with the $\ell_{p,0}$ -norm with $p\in\mathbb{Q}_{+}$ , implying that the proposed algorithm is also applicable to the $\ell_{p,0}$ -norm case with $p\in\mathbb{Q}_{+}$ . However, for the convenience of calculation, we just consider the $\ell_{2,0}$ -norm case.

6 Numerical simulations

We designed numerical experiments to verify the validity of the proposed algorithm and theorems and then explored the performance and advantages of our MMV-ADMM- $\ell_{2,0}$ by comparing it with existing algorithms.

6.1 Experiment setup

All experiments were based on synthetic data and were performed on a PC with an Intel Core i7-8565U 1.8 GHz CPU and 20 GB RAM. The results in Subsections 6.4-6.7 are averaged over 100 independent experiments. Our datasets are generated as follows: sensing matrix $\Phi\in\mathbb{R}^{M\times N}$ is an $i.i.d.$ Gaussian random matrix with unit-norm columns, and the ground truth of the recovery problem $S\in\mathbb{R}^{N\times J}$ is generated in two steps. First, $K$ rows were randomly selected as nonzero rows, with the remaining rows all being zero. Then, the elements of the selected $K$ rows were independently generated from $\mathcal{N}(0,1)$ . The solution obtained from the different recovery algorithms was denoted as $\hat{S}$ . The root-mean-square error (RMSE), given as

RMSE=\frac{\left\|\hat{S}-S\right\|_{F}}{\sqrt{NJ}},

was used to compare the recovery quality. The average running time was used to measure the computational cost. Recovery is considered successful when $RMSE<10^{-5}$ .

For comparison, we test SOMP [13] and MFOCUSS [7], both of which are traditional algorithms that show good performance for the MMV problem. We also include the MMV-SPG algorithm based on the $\ell_{2,1}$ -norm, which is derived from a solver for large-scale sparse reconstruction called SPGL1 [26], to solve the MMV problem. Because the projection operator $\mathcal{P}_{\mathcal{M}}(\cdot)$ in (4.3) is a hard thresholding operator, we also check the SNIHT [14] for the MMV problem. As the most direct comparison, MMV-ADMM- $\ell_{2,1}$ [6] is the same ADMM scheme algorithm based on the $\ell_{2,1}$ -norm. We introduce it in the experiments to compare the changes in performance after the $\ell_{2,1}$ -norm is replaced by the $\ell_{2,0}$ -norm.

The parameters of all algorithms are set as shown in Table 1. For MMV-ADMM- $\ell_{2,0}$ and MMV-ADMM- $\ell_{2,1}$ , $S^{0}$ and $L^{0}$ are separately initialized as an $i.i.d.$ Gaussian random matrix and a null matrix. For MMV-ADMM- $\ell_{2,0}$ , if we have the prior information that the joint-sparsity of $S$ is about $K$ , then we set $s$ to be slightly larger than $K$ . Otherwise, we use (3.3) to initialize $s$ . The spark of the matrix $\Phi$ is not easy to calculate accurately; however, an estimation given by [4] is enough for the proposed algorithm.

Algorithm	Parameters
SOMP[13]	Sparsity $K$
MFOCUSS[7]	$\lambda=10^{-10}$
MMV-SPG[26]	Not required
SNIHT[14]	Sparsity $K$ , $MaxIter=10^{3}$
MMV-ADMM- $\ell_{2,1}$ [6]	$\lambda=10^{-6},\rho=10^{-5},MaxIter=10^{3}$
MMV-ADMM- $\ell_{2,0}$	$\rho=1,MaxIter=10^{3}$

Table 1: Parameter settings for the tested algorithms

6.2 Validity of MMV-ADMM- $\ell_{2,0}$

In the first experiment, we observed the recovery quality of our MMV-ADMM- $\ell_{2,0}$ . MMV-ADMM- $\ell_{2,0}$ was applied without a convergence criterion ( $i.e.$ , iterates to $MaxIter=10^{3}$ ) on datasets of size $N=500,M=150,K=50,J=10$ over 10 random repetitions. As shown in Table 2, the $RMSE$ of MMV-ADMM- $\ell_{2,0}$ is close to the computer precision, and the running time is less than 1 s; this indicates that the proposed algorithm can recover the sparse signals precisely and efficiently.

Data	RMSE	Time(s)
1	8.9904e-16	0.6610
2	9.1155e-16	0.6290
3	9.1562e-16	0.5604
4	8.4978e-16	0.5965
5	9.1531e-16	0.6480
6	9.0638e-16	0.6051
7	9.9356e-16	0.5715
8	9.1655e-16	0.5679
9	7.7866e-16	0.6154
10	7.6170e-16	0.6352
Mean	8.8484e-16	0.6090
Std	6.9661e-17	0.0350

Table 2: The average recovery quality over 10 random experiments (

N=500,M=150,K=50,J=10

)

6.3 Test for convergence criterion

In the second experiment, we tested the convergence criterion (4.7) given by Theorem 5.3. We set $N=500,M=150,K=50,J=10$ . The aim was to observe the average change in the convergence criterion (4.7) in iterations over 10 random repetitions. The average change of the convergence criterion in 200 iterations is shown in Figure 1. All three variables in (4.7) tend to zero and reach an order of magnitude of $-16$ at $MaxIter=10^{3}$ ; this satisfies our analysis in Theorem 5.3.

Refer to caption — Figure 1: The average change of convergence criterion in iterations over 10 experiments ( $N=500,M=150,K=50,J=10$ )

Then, we studied the influence of the convergence criterion (4.7) on our algorithm. When all three variables in the convergence criterion (4.7) were less than $10^{-6}$ , the proposed algorithm was considered to have found an approximate optimal solution, and the iteration was stopped. As shown in Table 3, with the convergence criterion, the proposed algorithm can successfully recover the sparse signals in a short time.

In the following experiments, we applied MMV-ADMM- $\ell_{2,0}$ with the proposed convergence criterion.

data	Without convergence criterion		With convergence criterion
	RMSE	Time(s)	RMSE	Time(s)
1	8.2155e-16	0.5765	9.4005e-10	0.1560
2	9.0799e-16	0.5193	7.4102e-10	0.1179
3	7.6515e-16	0.4940	6.9234e-10	0.1116
4	8.4516e-16	0.5109	8.8110e-10	0.1132
5	8.5083e-16	0.5278	7.4768e-10	0.1269
6	9.4163e-16	0.5697	8.0561e-10	0.1508
7	8.2227e-16	0.5152	6.5987e-10	0.1108
8	8.9286e-16	0.5053	7.4714e-10	0.1069
9	7.6474e-16	0.4863	8.1573e-10	0.0949
10	8.0768e-16	0.4796	7.8128e-10	0.1024
Mean	8.4199e-16	0.5185	7.8118e-10	0.1191
Std	5.8515e-17	0.0324	8.4008e-11	0.0200

Table 3: The average recovery quality over 10 random experiments with and without convergence criterion(

N=500,M=150,K=50,J=10

)

6.4 Performance with different sparsity $K$

In the third experiment, we studied how sparsity influenced the recovery quality of all these algorithms. Set $N=500,M=150,J=10$ and let $K$ range from 25 to 150 with step size 25. Observe how the percentage of successful recovery changes with $\frac{K}{N}$ when applying different algorithms. The experiment results are shown in Figure 2.

It is not hard to see that all these algorithms fail to recover the original signals when the sparsity of signals is large ( $i.e.$ , $\frac{K}{N}\geq 0.2$ ). However, when the sparsity is not so large, that is, $0.1\leq\frac{K}{N}\leq 0.15$ , the proposed algorithm can successfully recover all signals over 100 experiments, whereas the other algorithms have poor performance. In addition, in the case of successful recovery, all these algorithms are fast enough. Therefore, we can conclude that for the MMV problem, when the original signals are not sparse enough, the proposed algorithm is more suitable.

6.5 Performance with different undersampling rate $\frac{N}{M}$

The main purpose of compressed sensing is to reduce the number of measurements while ensuring good recovery quality. To illustrate that the proposed algorithm has better performance with a larger undersampling rate compared with that of other algorithms, we studied the percentage of successful recovery with different undersampling rates. We set $N=500,K=50,J=10$ and let the undersampling rate $\frac{N}{M}$ range from 1.6 to 8 with a step size of 1.6.

As shown in Figure 3, when the undersampling rate $\frac{N}{M}\geq 6.4$ , all algorithms fail to recover the original signals. However, when $\frac{N}{M}=4.8$ , the proposed algorithm can still recover all signals, whereas MFOCUSS can only recover approximately 94% of the signals and other algorithms perform poorly. With regard to the running time, we just need to consider the case of successful recovery. We see that the proposed algorithm has average performance among all algorithms. This implies that the proposed algorithm has a larger undersampling rate compared with that of the other algorithms; this is the key of compressed sensing.

6.6 Performance with different number of sensors $J$

In this experiment, we studied how the recovery quality was affected by the number of sensors $J$ . We set $N=500,M=150,K=50$ and let $J$ range from 1 (multiplied by 2 per step) to 32.

Figure 4 shows that when the number of sensors is 1 ( $i.e.$ , a SMV problem), all algorithms cannot successfully recover all signals. However, the proposed algorithm can recover approximately 82% of signals, whereas the other algorithms perform poorly. When $J\geq 2$ , the proposed algorithm can successfully recover all signals. However, the other algorithms can only recover all signals when $J\geq 4$ . The running time of the proposed algorithm is the average of that of all the other algorithms, all of which are fast enough to recover the original signals. This experiment proves that the proposed algorithm performs well even though with a small number of sensors, whereas the other algorithms cannot.

6.7 Performance with different sample dimensions $N$

We tested the recovery quality of our algorithm for varying numbers of dimensions. We set $J=10$ and tested $N=100,500,1000,1500,3000$ with $\frac{N}{M}=3$ and $\frac{K}{N}=0.1$ .

Figure 5 shows that MFOCUSS, MMV-ADMM- $\ell_{2,1}$ , and the proposed algorithm all perform well irrespective of how the dimension $N$ changes. MFOCUSS and the proposed algorithm have similar running times; however, the running time of MMV-ADMM- $\ell_{2,1}$ is higher when $N$ is large. This experiment proves that the proposed algorithm is also suitable for a large number of dimensions.

6.8 Performance with SMW-formula

Finally, we designed an experiment to illustrate the advantage of MMV-ADMM- $\ell_{2,0}$ with the SMW formula when $N$ is large. We set $N=5000,M=1500,K=500,J=10$ and observe the RMSE and running time over 10 random experiments.

As shown in Table 4, with the SMW-formula, MMV-ADMM- $\ell_{2,0}$ -SMW requires less time to obtain solutions with the same precision as that of MMV-ADMM- $\ell_{2,0}$ .

data	MMV-ADMM- $\ell_{2,0}$		MMV-ADMM- $\ell_{2,0}$ -SMW
	RMSE	Time(s)	RMSE	Time(s)
1	2.5689e-10	12.5324	2.3941e-10	11.5070
2	2.6148e-10	13.4957	2.4949e-10	12.0098
3	2.5017e-10	13.3999	2.6850e-10	12.7097
4	2.5103e-10	13.8419	2.5707e-10	12.2158
5	2.4778e-10	13.6656	2.4574e-10	12.1301
6	2.4846e-10	13.5132	2.5370e-10	12.0958
7	2.4829e-10	13.5649	2.4438e-10	12.4199
8	2.5644e-10	14.2626	2.5308e-10	12.4762
9	2.6028e-10	13.6527	2.5655e-10	12.3948
10	2.4507e-16	13.8025	2.5998e-10	12.3597
Mean	2.5259e-10	13.5732	2.5279e-10	12.2319
Std	5.7252e-12	0.4394	8.4480e-12	0.3283

Table 4: The average recovery quality over 10 random experiments with and without SMW-formula(

N=5000,M=1500,K=500,J=10

)

7 Conclusions

In this study, we propose MMV-ADMM- $\ell_{2,0}$ , an ADMM algorithm based on $\ell_{2,0}$ -norm for the MMV problem. Instead of using the relaxation $\ell_{2,1}$ -norm, we directly solve the $\ell_{2,0}$ -norm minimization problem; this is the main difference between our study and other works [7, 11, 12]. We prove that the proposed algorithm has global convergence and determine its convergence criterion (4.7). Through numerical simulations, we test the validity of the proposed algorithm and its convergence criterion. The experiment results are consistent with Theorems 5.2 and 5.3. In comparisons against other algorithms, MMV-ADMM- $\ell_{2,0}$ shows better performance with high sparsity $K$ and few sensors $J$ , and it is insensitive to the dimension $N$ . Moreover, MMV-ADMM- $\ell_{2,0}$ has a larger undersampling rate compared with those of other algorithms, especially MMV-ADMM- $\ell_{2,1}$ , which is the key of compressed sensing.

The proposed algorithm has some drawbacks. Although its running time is roughly the average of those of the other algorithms, its efficiency remains lower than those of SOMP and SNIHT. We propose MMV-ADMM- $\ell_{2,0}$ -SMW to reduce the computational cost, and it runs faster. However, this improvement is not significant until the dimension $N$ becomes large.

The research in this paper can still be studied through a deep dive. First, Algorithm 1 can be extended to the generalised JSM-2 problem, where multiple vectors are measured by different sensing matrices. And the derivation of the algorithm for the generalized JSM-2 problem is nearly the same as discussions in Section 3 and 4 of this paper. Second, Algorithm 1 can accelerate even further. The most time cost step in the iterations of Algorithm 1 is updating $S$ by (4.5). The key is to solve a high-dimensional system of positive-definite linear equations. Drawing on recent advances in the field of numerical algebra, $S$ can be updated more efficiently, thus Algorithm 1 can converge faster. We will focus on the above two points in the future research; give a more efficient algorithm for the MMV problem, and extend it to the generalized JSM-2 problem.

References

[1] Ewout van den Berg and Michael P. Friedlander. Theoretical and empirical results for recovery from multiple measurements. IEEE Transactions on Information Theory, 56(5):2516–2527, 2010.
[2] D.L. Donoho. Compressed sensing. IEEE Transactions on Information Theory, 52(4):1289–1306, 2006.
[3] E.J. Candes, J. Romberg, and T. Tao. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory, 52(2):489–509, 2006.
[4] David L. Donoho and Michael Elad. Optimally sparse representation in general (nonorthogonal) dictionaries via $\ell_{1}$ minimization. Proceedings of the National Academy of Sciences of the United States of America, 100:2197 – 2202, 2003.
[5] Scott Saobing Chen, David L. Donoho, and Michael A. Saunders. Atomic decomposition by basis pursuit. SIAM J. Sci. Comput., 20:33–61, 1998.
[6] Qing Qu, Nasser M. Nasrabadi, and Trac D. Tran. Abundance estimation for bilinear mixture models via joint sparse and low-rank representation. IEEE Transactions on Geoscience and Remote Sensing, 52(7):4404–4423, 2014.
[7] Shane F. Cotter, Bhaskar D. Rao, Kjersti Engan, and Kenneth Kreutz-Delgado. Sparse solutions to linear inverse problems with multiple measurement vectors. IEEE Transactions on Signal Processing, 53:2477–2488, 2005.
[8] Irina F. Gorodnitsky, John S. George, and B.D. Rao. Neuromagnetic source imaging with focuss: a recursive weighted minimum norm algorithm. Electroencephalography and clinical neurophysiology, 95 4:231–51, 1995.
[9] Dmitry M. Malioutov, Müjdat Çetin, and Alan S. Willsky. Source localization by enforcing sparsity through a laplacian prior: an svd-based approach. IEEE Workshop on Statistical Signal Processing, 2003, pages 573–576, 2004.
[10] S.F. Cotter and B.D. Rao. Sparse channel estimation via matching pursuit with application to equalization. IEEE Transactions on Communications, 50(3):374–377, 2002.
[11] Jie Chen and Xiaoming Huo. Theoretical results on sparse representations of multiple-measurement vectors. IEEE Transactions on Signal Processing, 54(12):4634–4643, 2006.
[12] Yonina C. Eldar and Moshe Mishali. Robust recovery of signals from a structured union of subspaces. IEEE Transactions on Information Theory, 55(11):5302–5316, 2009.
[13] Joel A. Tropp, Anna C. Gilbert, and Martin Strauss. Algorithms for simultaneous sparse approximation. part i: Greedy pursuit. Signal Process., 86:572–588, 2006.
[14] Jeffrey D. Blanchard, Michael Cermak, David Hanle, and Yirong Jing. Greedy algorithms for joint sparse recovery. IEEE Transactions on Signal Processing, 62(7):1694–1704, 2014.
[15] Stephen P. Boyd, Neal Parikh, Eric King wah Chu, Borja Peleato, and Jonathan Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn., 3:1–122, 2011.
[16] Wen Zaiwen, Hu Jiang, Li Yongfeng, and Liu Haoyang. Optimization: Modeling, Algorithm and Theory(in Chinese). Higher Education Press, 2020.
[17] Stephen P. Boyd and Lieven Vandenberghe. Convex optimization. IEEE Transactions on Automatic Control, 51:1859–1859, 2004.
[18] Attouch H, Bolte J, and Svaiter B. F. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods. Mathematical Programming, 2013.
[19] Boris S. Mordukhovich. Variational analysis and generalized differentiation I: Basic theory. Grundlehren der mathematischen Wissenschaften. Springer Berlin, Heidelberg, 2006.
[20] Bolte J, Sabach S, and Teboulle M. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Math. Program. 146, 459–494, 2014.
[21] Simon Foucart. Hard thresholding pursuit: An algorithm for compressive sensing. SIAM J. Numer. Anal., 49:2543–2563, 2011.
[22] Jack Sherman and Winifred J. Morrison. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. Annals of Mathematical Statistics, 21:124–127, 1950.
[23] M Woodbury. Inverting Modified Matrices. Statistical Research Group, Princeton University, Princeton, 1950.
[24] Ke Guo, Deren Han, and Tingting Wu. Convergence of alternating direction method for minimizing sum of two nonconvex functions with linear constraints. International Journal of Computer Mathematics, 94:1653 – 1669, 2017.
[25] Oleg P. Burdakov, Christian Kanzow, and Alexandra Schwartz. Mathematical programs with cardinality constraints: Reformulation by complementarity-type conditions and a regularization method. SIAM J. on Optimization, 26(1):397–425, jan 2016.
[26] Ewout van den Berg and Michael P. Friedlander. Sparse optimization with least-squares constraints. SIAM J. Optim., 21:1201–1229, 2011.

Alternating Direction Method of Multipliers Based on ℓ2,0\ell_{2,0}-norm for Multiple Measurement Vector Problem

Abstract

1 Introduction

2 Preliminaries

Definition 2.1.

Definition 2.2 ([16]).

Definition 2.3 ([17]).

Definition 2.4 ([18]).

Proposition 2.1 ([20]).

3 Problem formulation

Proposition 3.1.

Proof.

Theorem 3.1 ([11]).

4 Algorithm

4.1 Update B

4.2 Update S

4.3 Update L

4.4 Convergence criterion

5 Convergence analysis

Proposition 5.1.

Proof.

Corollary 5.1.

Corollary 5.2.

Proof.

Proposition 5.2.

Proof.

Corollary 5.3.

Proof.

Lemma 5.1 ([24]).

Theorem 5.1 ([24]).

Theorem 5.2.

Proof.

Theorem 5.3.

Proof.

6 Numerical simulations

6.1 Experiment setup

6.2 Validity of MMV-ADMM-ℓ2,0\ell_{2,0}

6.3 Test for convergence criterion

6.4 Performance with different sparsity KK

6.5 Performance with different undersampling rate NM\frac{N}{M}

6.6 Performance with different number of sensors JJ

6.7 Performance with different sample dimensions NN

6.8 Performance with SMW-formula

7 Conclusions

References

Alternating Direction Method of Multipliers Based on $\ell_{2,0}$ -norm for Multiple Measurement Vector Problem

6.2 Validity of MMV-ADMM- $\ell_{2,0}$

6.4 Performance with different sparsity $K$

6.5 Performance with different undersampling rate $\frac{N}{M}$

6.6 Performance with different number of sensors $J$

6.7 Performance with different sample dimensions $N$