Joint-sparse recovery from multiple measurements

Introduction

A problem of central importance in compressed sensing is the following: given an $m\times n$ matrix $A$, and a measurement vector $b=Ax_{0}$, recover $x_{0}$. When $m<n$, this problem is ill-posed, and it is not generally possible to uniquely recover $x_{0}$ without some prior information. In many important cases, $x_{0}$ is known to be sparse, and it may be appropriate to solve

A natural extension of the single-measurement-vector (SMV) problem just described is the multiple-measurement-vector (MMV) problem. Instead of a single measurement $b$, we are given a set of $r$ measurements

in which the vectors $x_{0}^{(k)}$ are jointly sparse—i.e., have nonzero entries at the same locations. Such problems arise in source localization , neuromagnetic imaging , and equalization of sparse-communication channels . Succinctly, the aim of the MMV problem is to recover $X_{0}$ from observations $B=AX_{0}$, where $B=[b^{(1)},\ b^{(2)},\ldots,\ b^{(r)}]$ is an $m\times r$ matrix, and the $n\times r$ matrix $X_{0}$ is row sparse—i.e., it has nonzero entries in only a small number of rows. The most widely studied approach to the MMV problem is based on solving the convex optimization problem

and X^{{j}{\scalebox{0.6}{\rightarrow}}} is the (column) vector whose entries form the $j$th row of $X$. In particular, Cotter et al. consider $p=2$, $q\leq 1$; Tropp analyzes $p=1$, $q=\infty$; Malioutov et al. and Eldar and Mishali use $p=1$, $q=2$; and Chen and Huo study $p=1$, $q\geq 1$. A different approach is given by Mishali and Eldar , who propose the ReMBo algorithm, which reduces MMV to a series of SMV problems.

The conditions under which (1.2) gives the sparsest possible solution have been studied by applying a number of different techniques. By far the most popular analytical approach is based on the restricted isometry property, introduced by Candès and Tao , which gives sufficient conditions for equivalence. Donoho obtains necessary and sufficient (NS) conditions by analyzing the underlying geometry of (1.2). Several authors characterize the NS-conditions in terms of properties of the kernel of $A$:

Fuchs and Tropp express sufficient conditions in terms of the solution of the dual of (1.2):

In this paper we are mainly concerned with the geometric and kernel conditions. We use the geometrical interpretation of the problems to get a better understanding, and resort to the null-space properties of $A$ to analyze recovery. To make the discussion more self-contained, we briefly recall some of the relevant results in the next three sections.

When we need to find the recoverability of vectors restricted to a support $\mathcal{I}$, this probability becomes

where $\mathcal{F}_{\mathcal{I}}(\mathcal{C})=2^{|\mathcal{I}|}$ denotes the number of faces in $\mathcal{C}$ formed by the convex hull of $\{\pm e_{j}\}_{i\in\mathcal{I}}$, and $\mathcal{F}_{\mathcal{I}}(A\mathcal{C})$ is the number of faces on $A\mathcal{C}$ generated by \{\pm A^{\scalebox{0.6}{\downarrow}{j}}\}_{j\in\mathcal{I}}.

Equivalence results in terms of null-space properties generally characterize equivalence for the set of all vectors $x$ with a fixed support, which is defined as

Sufficient conditions for recovery can be derived from the first-order optimality conditions necessary for $x^{*}$ and $y^{*}$ to be solutions of (1.2) and (2.1) respectively. The Karush-Kuhn-Tucker (KKT) conditions are also sufficient in this case because the problems are convex. The Lagrangian function for (1.2) is given by

where $\partial_{x}\mathcal{L}$ denotes the subdifferential of $\mathcal{L}$ with respect to $x$. The second condition reduces to

Recovery using sums-of-row norms

Our analysis of sparse recovery for the MMV problem of recovering $X_{0}$ from $B=AX_{0}$ begins with an extension of Theorem 2.1 to recovery using the convex relaxation

note that the norm within the summation is arbitrary. Define the row support of a matrix as

With these definitions we have the following result. (A related result is given by Stojnic et al. .)

For the “only if” part, suppose that there is a $Z$ with columns Z^{\scalebox{0.6}{\downarrow}{k}}\in\textrm{Ker}(A)\setminus\{0\} such that (3.2) does not hold. Now, choose X^{{j}{\scalebox{0.6}{\rightarrow}}}=Z^{{j}{\scalebox{0.6}{\rightarrow}}} for all $j\in\mathcal{I}$ and with all remaining rows zero. Set $B=AX$. Next, define $V=X-Z$, and note that $AV=AX-AZ=AX=B$. The construction of $V$ implies that \sum_{j}\|X^{{j}{\scalebox{0.6}{\rightarrow}}}\|\geq\sum_{j}\|V^{{j}{\scalebox{0.6}{\rightarrow}}}\|, and consequently $X$ cannot be the unique solution of (3.1).

Applying the reverse triangle inequality, $\|a+b\|-\|b\|\geq-\|a\|$, to the summation over $j\in\mathcal{I}$ and reordering exactly gives condition (3.2). ∎

For uniform recovery on support $\mathcal{I}$ to hold it follows from Theorem 3.1 that for any matrix $Z$ with columns Z^{\scalebox{0.6}{\downarrow}{k}}\in\textrm{Ker}(A)\setminus\{0\}, property (3.2) holds. In particular it holds for $Z$ with Z^{\scalebox{0.6}{\downarrow}{k}}={\bar{z\mkern 2.8mu}\mkern-2.8mu}{} for all $k$, with ${\bar{z\mkern 2.8mu}\mkern-2.8mu}{}\in\textrm{Ker}(A)\setminus\{0\}$. Note that for these matrices there exist a norm-dependent constant $\gamma$ such that

where $\Braket{V,W}\mathrel{\mathop{:}}=\mathop{\hbox{\rm trace}}(V^{T}\!W)$ is an inner-product defined over real matrices. The dual is then given by maximizing

over $Y$. (Because the primal problem has only linear constraints, there necessarily exists a dual solution $Y^{*}$ that maximizes this expression [25, Theorem 28.2].) To simplify the supremum term, we note that for any convex, positively homogeneous function $f$ defined over an inner-product space,

To derive these conditions, note that positive homogeneity of $f$ implies that $f(0)=0$, and thus $w\in\partial f(0)$ implies that $\Braket{w,v}\leq f(v)$ for all $v$. Hence, the supremum is achieved with $v=0$. If on the other hand $w\not\in\partial f(0)$, then there exists some $v$ such that $\Braket{w,v}>f(v)$, and by the positive homogeneity of $f$, $\Braket{w,\alpha v}-f(\alpha v)\to\infty$ as $\alpha\to\infty$. Applying this expression for the supremum to (3.5), we arrive at the necessary condition

Combining this expression with (3.6), we arrive at the dual of (3.3):

The following conditions are therefore necessary and sufficient for a primal-dual pair $(X^{*},Y^{*})$ to be optimal for (3.3) and its dual (3.8):

The existence of a matrix $Y^{*}$ that satisfies (3.9) provides a certificate that the feasible matrix $X^{*}$ is an optimal solution of (3.3). However, it does not guarantee that $X^{*}$ is also the unique solution. The following theorem gives sufficient conditions, similar to those in Section 2.3, that also guarantee uniqueness of the solution.

The first three conditions clearly imply that $(X,Y)$ primal and dual feasible, and thus satisfy (3.9a) and (3.9b). Conditions (3.10b) and (3.10c) together imply that

The first and last identities above follow directly from the definitions of the matrix trace and of the norm $\|\cdot\|_{1,2}$, respectively; the middle equality follows from the standard Cauchy inequality. Thus, the zero-gap requirement (3.9c) is satisfied. The conditions (3.10a)–(3.10c) are therefore sufficient for $(X,Y)$ to be an optimal primal-dual solution of (3.3). Because $Y$ determines the support and is a Lagrange multiplier for every solution $X$, this support must be unique. It then follows from condition (3.10d) that $X$ must be unique. ∎

2 Counter examples

3 Experiments

Recovery using ReMBo

which is the number of nonzeros of the sparsest vector in the kernel of $A$; any vector $x_{0}$ with fewer than $\textrm{Spark}(A)/2$ nonzeros is the unique sparsest solution of $Ax=Ax_{0}=b$ . Unfortunately, the spark is prohibitively expensive to compute, but under the assumption that $A$ is in general position, $\textrm{Spark}(A)=m+1$. Note that choosing a higher value can help to recover signals with row sparsity exceeding $m/2$. However, in this case it can no longer be guaranteed to be the sparsest solution.

The term within brackets denotes the probability of failure and the fraction represents the success rate, which is given by the ratio of the number of faces $\mathcal{F}_{\mathcal{I}}(A\mathcal{C})$ that survived the mapping to the total number of faces to consider. The total number reduces by two at each trial because we can exclude the face $f$ we just tried, as well as $-f$. The factor of two in $C(|\mathcal{I}|,r)/2$ is also due to this symmetryHenceforth we use the convention that the uniqueness of a sign pattern is invariant under negation..

This model would be a bound for the average performance of ReMBo if the sign patterns generated would be randomly sampled from the space of all sign patterns on the given support. However, because it is generated from the orthant intersections with a hyperplane, the actual pattern is highly structured. Indeed, it is possible to imagine a situation where the $(s-1)$-faces in $\mathcal{C}$ that perish in the mapping to $A\mathcal{C}$ have sign patterns that are all contained in the set generated by a single hyperplane. Any other set of sign patterns would then necessarily include some faces that survive the mapping and by trying all patterns in that set we would recover $X_{0}$. In this case, the average recovery over all $X_{0}$ on that support could be much higher than that given by (5.1). We do not yet fully understand how the surviving faces of $\mathcal{C}$ are distributed. Due to the simplicial structure of the facets of $\mathcal{C}$, we can expect the faces that perish to be partially clustered (if a $(d-2)$-face perishes, then so will the two $(d-1)$-faces whose intersection gives this face), and partially unclustered (the faces that perish while all their sub-faces survive). Note that, regardless of these patterns, recovery is guaranteed in the limit whenever the number of unique sign patterns tried exceeds half the number of faces lost, $(|\mathcal{F}_{\mathcal{I}}(\mathcal{C})|-|\mathcal{F}_{\mathcal{I}}(\mathcal{AC})|)/2$.

where the last inequality follows from (5.3). Consequently, all inequalities hold with equality. ∎

Given $d\leq n$, then $C(n,d)=2^{n}-C(n,n-d)$, and $C(2d,d)=2^{2d-1}$.

2 Practical considerations

In practice it is generally not feasible to generate all of the $C(|\mathcal{I}|,r)/2$ unique sign patterns. This means that we would have to replace this term in (5.1) by the number of unique patterns actually tried. For a given $X_{0}$ the actual probability of recovery is determined by a number of factors. First of all, the linear combinations of the columns of the nonzero part of $\bar{X}$ prescribe a hyperplane and therefore a set of possible sign patterns. With each sign pattern is associated a face in $\mathcal{C}$ that may or may not map to a face in $A\mathcal{C}$. In addition, depending on the probability distribution from which the weight vectors $w$ are drawn, there is a certain probability for reaching each sign pattern. Summing the probability of reaching those patterns that can be recovered gives the probability $P(A,\mathcal{I},X_{0})$ of recovering with an individual random sample $w$. The probability of recovery after $t$ trials is then of the form

3 Experiments

In this section we illustrate the theoretical results from Section 5 and examine some practical considerations that affect the performance of ReMBo. For all experiments that require the matrix $A$, we use the same $20\times 80$ matrix that was used in Section 4, and likewise for the supports $\mathcal{I}_{s}$. To solve (1.2), we again use CVX in conjunction with SDPT3. We consider $x_{0}$ to be recovered from $b=Ax_{0}=AX_{0}w$ if $\|x^{*}-x_{0}\|_{\infty}\leq 10^{-5}$, where $x^{*}$ is the computed solution.

The experiments that are concerned with the number of unique sign patterns generated depend only on the $s\times r$ matrix $\bar{X}$ representing the nonzero entries of $X_{0}$. Because an initial reordering of the rows does not affect the number of patterns, those experiments depend only on $\bar{X}$, $s=|\mathcal{I}|$, and the number of observations $r$; the exact indices in the support set $\mathcal{I}$ are irrelevant for those tests.

The practical performance of ReMBo depends on its ability to generate as many different sign patterns using the columns in $X_{0}$ as possible. A natural question to ask then is how the number of such patterns grows with the number of randomly drawn samples $w$. Although this ultimately depends on the distribution used for generating the entries in $w$, we shall, for sake of simplicity, consider only samples drawn from the normal distribution. As an experiment we take a $10\times 5$ matrix $\bar{X}$ with normally-distributed entries, and over $10^{8}$ trials record how often each sign-pattern (or negation) was reached, and in which trial they were first encountered. The results of this experiment are summarized in Figure 7. From the distribution in Figure 7(b) it is clear that the occurrence levels of different orthants exhibits a strong bias. The most frequently visited orthant pairs were reached up to $7.3\times 10^{6}$ times, while others, those hard to reach using weights from the normal distribution, were observed only four times over all trials. The efficiency of ReMBo depends on the rate of encountering new sign patterns. Figure 7(c) shows how the average rate changes over the number of trials. The curves in Figure 7(d) illustrate the theoretical probability of recovery in (5.1), with $C(n,d)/2$ replaced by the number of orthant pairs at a given iteration, and with face counts determined as in Section 4, for three instances with support cardinality $s=10$, and observations $r=5$.

3.2 Role of X¯¯𝑋\bar{X}.

3.3 Limiting the number of iterations

The number of iterations used in the previous experiments greatly exceeds that what is practically feasible: we cannot afford to run ReMBo until all possible sign patterns have been tried, even if there was a way detect that the limit had been reached. Realistically, we should set the number of iterations to a fixed maximum that depends on the computational resources available, and the problem setting.

Conclusions

All of the numerical experiments in this paper are reproducible. The scripts used to run the experiments and generate the figures can be downloaded from

Acknowledgments

The authors would like to give their sincere thanks to Özgür Yılmaz and Rayan Saab for their thoughtful comments and suggestions during numerous discussions.

References