Average Case Analysis of Multichannel Sparse Recovery Using Convex Relaxation

I Introduction

Recovery of sparse signals from a small number of measurements is a fundamental problem in many different signal processing tasks such as image denoising , analog-to-digital conversion , radar, compression, inpainting, and many more. The recent framework of compressed sensing (CS), founded in the works of Donoho , Candès, Romberg and Tao , studies acquisition methods as well as efficient computational algorithms that allow reconstruction of a sparse vector $x$ from linear measurements $y=Ax$, where $A\in{\hbox{\msbm{R}}}^{n\times N}$ is referred to as the measurement matrix. The key observation is that $y$ can be relatively short, so that $n<N$, and still contain enough information to recover $x$.

In order to capture more closely the true underlying behavior of existing algorithms and observe a performance gain when using several channels, we consider an average-case analysis. In this setting, the inputs are considered to be random variables. The idea is to develop conditions on the measurement matrix $A$ such that the inputs can be recovered with high probability given a certain input distribution.

The theoretical average-case results we develop for multichannel BP are superior to the average bounds developed on thresholding and SOMP. For an equally mild or even milder condition on the sparsity and on the matrix $A$, we obtain faster exponential decay of the failure probability with respect to the number of channels. Thus, in this sense, the extension of BP to the multichannel case is superior to existing greedy algorithms, just as in the single channel setting. Moreover, our recovery results are applicable also in the single channel case whereas previous results require a large number of channels to yield meaningful (i.e., positive) probability bounds (although our new bound for thresholding generalizing the one in does not suffer from this drawback). Note, however, that in simulations SOMP often exhibits the best performance. This may be explained by the fact that the bounds are not tight (at least for SOMP).

We consider multichannel signal recovery where our goal is to recover a jointly-sparse matrix $X\in{\hbox{\msbm{C}}}^{N\times L}$ from $n$ linear measurements per channel. Here $N$ denotes the signal length and $L$ the number of channels, i.e., the number of signals. We assume that $X$ is jointly $k$-sparse, meaning that there are at most $k$ rows in the matrix $X$ that are not identically zero. More formally, we define the support of the matrix $X$ as

Our assumption is that $\|X\|_{0}:=|\operatorname{supp}X|\leq k$. The measurements are given by

which promotes joint sparsity, as argued for instance in . In the single channel case $L=1$ this is the usual BP principle. Therefore, our results can also be used to deduce the average-case behavior of the BP method. This is in contrast to , in which the recovery results derived are not applicable to the single channel case. As we discuss in Section VI, our theoretical results are superior to the previous average-case analysis of in the sense that we use an equally mild or even milder condition on the sparsity and on the matrix $A$, but at the same time get a faster exponential decay of the failure probability with respect to the number of channels $L$.

II-B Recovery Results

Recovery results for the program (5) were considered in . In particular, the lemma below is derived in and follows also from where the more general case of block sparsity is considered.

Let $S\subset\{1,\ldots,N\}$ and suppose that

with $A_{S}^{\dagger}=(A_{S}^{*}A_{S})^{-1}A_{S}^{*}$ denoting the pseudo-inverse of $A_{S}$. Then (5) recovers all $X\in{\hbox{\msbm{C}}}^{N\times L}$ with $\operatorname{supp}X=S$ from $Y=AX$.

Note, that the condition above does not depend on the number of channels. In the next section we will derive a condition similar to (6) that involves the $2$-norm instead of the $1$-norm, and is therefore weaker (namely, easier to satisfy).

The following result follows from by noting that the block coherence in this setting is equal to $\mu/d$.

Then (5) recovers all $X$ with $\|X\|_{0}\leq k$ from $Y=AX$.

Under the same conditions as in Propositions II.1 and II.2, it is shown in that BP will recover a single $k$-sparse vector. Therefore, if (6) holds, then instead of solving (5) we can use BP on each of the columns of $Y$.

The lower bound behaves like $1/\sqrt{n}$ for large $N$, which limits the Proposition II.2 to maximal sparsities $k={\cal{O}}(\sqrt{n})$. To improve on this we can generalize existing recovery results based on RIP to the multichannel setup. The restricted isometry constant $\delta_{k}$ of a matrix $A$ is defined to be the smallest constant $\delta_{k}$ such that

for all $k$-sparse vectors $x$. The next proposition follows from .

Assume $A\in{\hbox{\msbm{C}}}^{n\times N}$ with $\delta_{2k}<\sqrt{2}-1$ Let $X\in{\hbox{\msbm{C}}}^{N\times L}$, $Y=AX$, and let $\overline{X}$ be the minimizer of (5). Then

It is well known that Gaussian and Bernoulli random matrices $A\in{\hbox{\msbm{R}}}^{n\times N}$ satisfy $\delta_{2k}\leq\sqrt{2}-1$ with high probability as long as

For random partial Fourier matrices the respective condition is $n\geq ck\log^{4}(N)$ . Therefore, Proposition II.3 allows for a smaller number of measurements. However, there is still no dependency on the number of channels. Indeed, under the same RIP condition BP will recover a single $k$-sparse vector and therefore, as before, BP may as well be applied to each of the columns of $Y$ individually.

III A Recovery Condition

Before turning to analyze the average-case behavior of (5), we first develop a new condition on $A$ that allows for perfect recovery. This formulation will be useful in deriving the average-case results.

In the following theorem we give a sufficient condition on the minimizers of (5). This theorem generalizes a result of for the $L=1$ case. To this end we denote by $\operatorname{sgn}(X)\in{\hbox{\msbm{C}}}^{N\times L}$ the matrix with entries

In this definition, each element of $X$ is normalized by the norm of the corresponding row. When $L=1$, $\operatorname{sgn}(X)$ reduces to the sign of the elements of the vector $x$.

Let $X\in{\hbox{\msbm{C}}}^{N\times L}$ with $\operatorname{supp}X=S$ and assume $A_{S}$ to be non-singular. If there exists a matrix $H\in{\hbox{\msbm{C}}}^{n\times L}$ such that

Before proving the theorem we note that the two conditions on $H$ easily imply that

Let $Y=AX$, and assume there exists a matrix $H$ such that $X,H$ satisfy (13) and (14). Let $X^{\prime}$ be an alternative matrix satisfying $Y=AX^{\prime}$. Our goal is to show that $\|X\|_{2,1}<\|X^{\prime}\|_{2,1}$. To this end, we note that

where $\operatorname{Tr}$ denotes the trace. Substituting $A_{S}^{*}H=\operatorname{sgn}(X^{S})$ into (16), and using the cyclicity of the trace we have

where we used the fact that $A_{S}X^{S}=Y=AX^{\prime}$ and $S^{\prime}$ denotes the support of $X^{\prime}$. We next rely on the following lemma.

where the second inequality is a result of applying Cauchy-Schwartz. Under the condition of the lemma, we have strict inequality in the last inequality. ∎

Thus, we have shown that $\|X^{\prime}\|_{2,1}>\|X\|_{2,1}$ for any $X^{\prime}$ such that $Y=AX^{\prime}$, and therefore (5) recovers the true sparse matrix $X$. ∎

Choosing $H=(A_{S}^{\dagger})^{*}\operatorname{sgn}(X_{S})$ in Theorem III.1 results in the following corollary.

Let $X\in{\hbox{\msbm{C}}}^{N\times L}$ with $\operatorname{supp}X=S$ and assume $A_{S}$ to be non-singular. If

This corollary will be instrumental in proving the average-case performance of (5). It can easily be seen that Corollary III.3 implies Proposition II.1. This follows from the triangle inequality,

where we used the fact that $\|\operatorname{sgn}(X^{j})\|_{2}=1$.

IV Average Case Analysis

Realizing that (5) is not more powerful than usual BP in the worst case, we seek an average-case analysis. This means that we impose a probability model on the $k$-sparse $X$. In particular, as in , we will assume that on the support $S$ of size $k$ the coefficients of $X$ are chosen at random. We then show that under a suitable probability model on the non-zero elements of $X$, the condition given by Corollary III.3 is satisfied with high probability, which depends on $L$.

We follow the probability model used in : let $S$ be the joint support of cardinality $k$. On $S$ the coefficients are given by

where $\Sigma=\operatorname{diag}(\sigma_{j},j\in S)\in{\hbox{\msbm{R}}}^{k\times k}$ is an arbitrary diagonal matrix with positive diagonal elements $\sigma_{j}$. The matrix $\Phi$ will be chosen at random according to one of the following models.

Real Gaussian: each entry of $\Phi\in{\hbox{\msbm{R}}}^{k\times L}$ is chosen independently from a standard normal distribution.

Real spherical: the rows of $\Phi\in{\hbox{\msbm{R}}}^{k\times L}$ are chosen independently and uniformly at random from the real sphere $S^{L-1}$.

Complex Gaussian: the real and imaginary parts of each entry of $\Phi\in{\hbox{\msbm{C}}}^{k\times L}$ are chosen independently according to a standard normal distribution.

Complex spherical: the rows of $\Phi\in{\hbox{\msbm{C}}}^{k\times L}$ are chosen independently and uniformly at random from the complex sphere $S_{\hbox{\msbm{C}}}^{L-1}$.

Before stating the first theorem, we derive the following result on the norm of sums of independent random vectors, uniformly distributed on a sphere.

Let $a\in{\hbox{\msbm{C}}}^{k}$ and let $Z_{j}$, $j=1,\ldots,k$, be a sequence of independent random vectors which are uniformly distributed on the real sphere $S^{L-1}$. Then for any $u>1$

Theorem IV.2 generalizes the Bernstein inequality for Steinhaus sequences in [46, Theorem 13] to higher dimensions. We may extend the estimate easily to random vectors uniformly distributed on complex unit spheres.

Let $a\in{\hbox{\msbm{C}}}^{k}$ and let $Z_{j}$, $j=1,\ldots,k$, be a sequence of independent random vectors which are uniformly distributed on the complex sphere $S^{L-1}_{\hbox{\msbm{C}}}$. Then for any $u>1$

First observe that $a_{j}Z_{j}$ has the same distribution as $|a_{j}|Z_{j}$. We may therefore assume without loss of generality that $a_{j}\in{\hbox{\msbm{R}}}$. Next, a random vector $Z\in S_{\hbox{\msbm{C}}}^{L-1}$ is uniformly distributed on $S_{\hbox{\msbm{C}}}^{L-1}$ if and only if $(\operatorname{Re}(Z)^{T},\operatorname{Im}(Z)^{T})^{T}$ is uniformly distributed on the real sphere $S^{2L-1}$. Applying Theorem IV.2 with $L$ replaced by $2L$ yields the statement. ∎

With this tool at hand we can now easily prove the following average-case recovery theorem.

Let $S\subset\{1,\ldots,N\}$ be a set of cardinality $k$ and suppose

Let $X\in{\hbox{\msbm{R}}}^{N\times L}$ with $\operatorname{supp}X\subset\{1,\ldots,N\}$ such that the coefficients on $S$ are given by (21) with some diagonal matrix $\Sigma\in{\hbox{\msbm{R}}}^{k\times k}$ and $\Phi\in{\hbox{\msbm{R}}}^{k\times L}$ chosen from the real Gaussian or spherical probability. Then with probability at least

If the real probability model is replaced by one of the two complex models then $L/2$ can be replaced by $L$ in (23).

For $\alpha<1$ we are guaranteed that the exponent in (23) has a negative argument, and therefore the error decays exponentially in $L$.

The complex case follows analogously using Corollary IV.3. ∎

For $L=1$, Theorem IV.4 is contained implicitly in [46, Theorem 13]. The appearance of the $2$-norm in (24) instead of the $1$-norm as in (6) makes the condition of the theorem weaker than worst-case estimates (recall that $\|x\|_{2}\leq\|x\|_{1}\leq\sqrt{k}\|x\|_{2}$ for any length-$k$ vector $x$). In Section V this will be made more evident when we consider conditions on the coherence $\mu$ and the RIP constant to allow for recovery with high probability. The requirement we obtain on $\mu$ is weaker than that of Proposition II.2 and allows for recovery with $k$ on the order of $n$, while the worst-case results limit recovery to order $\sqrt{n}$. Furthermore, in contrast to the worst-case results which depend on $\delta_{2k}$, we will show that high-probability recovery is possible as long as $\delta_{k+1}$ is small enough.

This provides a useful average-case analysis even for $L=1$.

Let $S\subset\{1,\ldots,N\}$ be a set of cardinality $k$, and let $X\in{\hbox{\msbm{R}}}^{N\times L}$ be random sparse coefficients with $\operatorname{supp}X=S$ given by the real Gaussian probability model. If

and $\Gamma$ denotes the Gamma function, then with probability at least

It follows from Stirling’s formula $\Gamma(z)\sim\sqrt{2\pi z}z^{z-1/2}e^{-z}$, that

Moreover, for all $L\geq 1$ it holds that $\sqrt{L}\geq A_{L}\geq\sqrt{\frac{2}{\pi}}\sqrt{L}\approx 0.797\sqrt{L}$.

Note that $\gamma=\frac{A_{L}}{3\sqrt{L}+2\sqrt{k}}$ is monotonically increasing in $L$. In addition, the probability $P$ is also increasing (towards $1$) in $L$. Therefore, more channels increase the probability of success and in addition relax the requirements on the matrix $A$.

To prove the theorem we show that if (24) is satisfied, then condition (20) of Corollary III.3 holds with probability $P$.

By the assumption of the theorem $\|b_{j}\|_{2}<\gamma$ where $\gamma$ is defined by (24). It therefore remains to bound $\|\Phi\|_{2}$ and $\|D\|_{2}$. ¿From [10, equation (4.35)], see also , the operator norm of $\Phi$ satisfies

with probability at least $1-\exp(-r^{2}/2)$.

Next we consider $\|D\|_{2}$. Observe that the $s_{j}^{2}$ are $\chi^{2}(L)$ distributed. Therefore, denoting a $\chi^{2}(L)$-variable by $Y$,

As a function of $\Phi^{j}$ the $s_{j}$ are Lipschitz continuous, i.e., $s_{j}(\Phi^{j}-\Psi^{j})\leq\|\Phi^{j}-\Psi^{j}\|_{2}$. Using these two observations we rely on the following standard concentration of measure result, see e.g. [28, eq. (2.35)] or [29, eq. (1.6)].

Let $f$ be a Lipschitz function on ${\hbox{\msbm{R}}}^{L}$, i.e., $|f(x)-f(y)|\leq B\|x-y\|_{2}$ for all $x,y\in{\hbox{\msbm{R}}}^{L}$. Further assume that $Z=(Z_{1},Z_{2},\ldots,Z_{L})$ is a vector of independent standard Gaussian random variables. Then

Our goal is to show that $\|D\|_{2}$ is bounded from above, which is equivalent to bounding the smallest value of $s_{j}$ from below. Applying Theorem IV.6 to $s_{j}$,

where we used the fact that $B=1$ and ${\hbox{\msbm{E}}}[s_{j}]=A_{L}$. Using a union bound over all $j$, we obtain

Assuming that $\min_{j\in S}s_{j}\geq A_{L}(1-t)$ holds, $\|D\|_{2}\leq 1/(A_{L}(1-t))$. Combining this bound with (26) for $r=\sqrt{L}s$ we have

¿From (27) and Corollary III.3, $X$ is recoverable using (5).

The probability that (27) does not hold can be computed by applying a union bound to the probabilities that the spectral norms of each of the matrices $\Phi$ and $D$ are not bounded. This shows that (27) does not hold with probability at most $\exp(-L/8)+k\exp(-A_{L}^{2}/8)$ completing the proof of the theorem. ∎

V Bounded Norm Condition

Let $A\in{\hbox{\msbm{C}}}^{n\times N}$ have unit-norm columns and coherence $\mu$, and let $S\subset\{1,\ldots,N\}$ be a set of cardinality $k$. Assume that

Gershgorin’s disk theorem implies that the smallest eigenvalue $\lambda_{\min}$ of $A_{S}^{*}A_{S}$ is bounded from below by $1-(k-1)\mu$. In particular, $A_{S}^{*}A_{S}$ is invertible provided $(k-1)\mu<1$. Further,

where the last inequality follows from the fact that (28) implies $\delta>\sqrt{k}/(1-(k-1)\mu)^{-1}$. ∎

Condition (28) is slightly weaker than (8) as long as $\delta>1/\sqrt{k}$. This follows from the $2$-norm that replaced the $1$-norm in the upper bound. However, (28) still suffers the square-root bottleneck $k={\cal O}(\sqrt{n})$. To improve on this result, we next provide a condition based on the following refinement of the RIP of $A$. For a set $S\subset\{1,\ldots,N\}$ we let

The restricted isometry constant $\delta_{k}$ of (10) satisfies $\delta_{k}=\max_{|S|\leq k}\|A_{S}^{*}A_{S}-I\|_{2}$ so that if $S$ has cardinality $k$ then $\delta(S)\leq\delta_{k}$. We further define

Clearly, $\delta(S)\leq\delta^{*}(S)\leq\delta_{k+1}$. Finally, we make use of the following “local” $2$-coherence function,

If $A$ satisfies $\delta^{*}(S)\leq\delta<1/2$ then

If $A$ satisfies $\delta(S)\leq\delta<1$ and $\mu_{2}(S)\leq\eta$ then

Denoting by $\lambda$ an eigenvalue of $A_{S}^{*}A_{S}$, the definition of $\delta(S)\leq\delta^{*}(S)\leq\delta$ implies that $|1-\lambda|\leq\delta$. Consequently, the smallest eigenvalue of $A_{S}^{*}A_{S}$ is bounded from below by $1-\delta$ and therefore

Proposition V.2 applies if $\delta_{k+1}$ is small while in contrast Theorem II.3 works with $\delta_{2k}$, which is generally larger than $\delta_{k+1}$. By (11) the condition $\delta_{k+1}\leq\delta$ can be satisfied if $n\geq C_{\delta}k\log(N/k)$. Working with $\delta^{*}(S)$ instead of $\delta_{k+1}$ allows to improve on the bound (11) for Gaussian, Bernoulli and random spherical matrices.

Let $S\subset\{1,\ldots,N\}$ be a set of cardinality $k$ and suppose that $A=\frac{1}{\sqrt{n}}\Phi\in{\hbox{\msbm{R}}}^{n\times N}$, where $\Phi$ is drawn at random according to a standard Gaussian or Bernoulli distribution (with expectation and variance $1/n$). Then $\delta^{*}(S)\leq\delta$ with probability at least $1-\epsilon$ provided that

The same statement holds (with possibly a different constant) for a random matrix whose columns are chosen independently at random according to the uniform distribution on a sphere.

A straightforward extension of the proof, as in , also shows that a random matrix $A\in{\hbox{\msbm{R}}}^{n\times N}$ with independent columns drawn from the uniform distribution on the sphere satisfies RIP, $\delta_{k}\leq\delta$ with probability at least $1-\epsilon$ provided $n\geq C\delta^{-2}(k\log(N/k)+\log(\epsilon^{-1}))$. Although this fact seems to be known, we are not aware of reference where this is rigorously stated.

The next result relies on a theorem by Tropp [46, Theorem B] that uses random support sets $S$ and allows to work with the coherence $\mu$ alone. Note that choosing $S$ at random is perfectly in line with an average-case analysis.

Let $A\in{\hbox{\msbm{C}}}^{n\times N}$ have unit norm columns and coherence $\mu$. Let $S\subset\{1,\ldots,N\}$ be a set of cardinality $k\geq 4$ chosen uniformly at random. Let $\delta,\epsilon\in(0,1)$ and assume that

where $c=\frac{\log(2)e^{-1/2}}{4\cdot 144\log(3)}\approx 6.64\cdot 10^{-4}$. Then

The proof relies on [46, Theorem 12]. The formulation below follows from by setting $s=\log(\epsilon^{-1})/\log(k/2)$ and estimating $\log(k/2+1)/\log(k/2)\leq\log(3)/\log(2)$ for $k\geq 4$.

Assume $A\in{\hbox{\msbm{C}}}^{n\times N}$ has unit norm columns and coherence $\mu$. Let $S\subset\{1,\ldots,N\}$ be a set of cardinality $k\geq 4$ chosen uniformly at random. The condition

Using (34) and the value of $c$, the square-root in (36) becomes $\delta/(2e^{1/4})$. Combining this with (35) shows that (36) is satisfied. Therefore, $\|A_{S}^{*}A_{S}-I\|_{2}\leq\delta$ with probability at least $1-\epsilon$, which implies that

Let us now compare worst-case and average results based on the coherence $\mu$, by relying on Theorem V.4. For simplicity, we consider the case in which $A$ is a unit-norm tight frame, for which $\|A\|^{2}_{2}=\frac{N}{n}$. In this case, (35) is equivalent to $k\leq\frac{\delta}{4e^{1/4}}n$. If additionally $\mu=c/\sqrt{n}$, then conditions (34) and (35) are both satisfied for fixed $\epsilon,\delta$ provided

This beats the square-root bottleneck and even removes the $\log$-factor present in estimates for the restricted isometry constants, see (11). Moreover, we have the additional advantage that the coherence is much easier to estimate than the restricted isometry constants.

Combining Theorem V.4 with the average-case analysis of Theorems IV.4 and IV.5 shows that for a unit norm tight frame $A$ of coherence $\mu$ multichannel sparse recovery by (5) can be ensured in the average-case provided $k\leq C\mu^{-2}$, which can be as small as $k\leq Cn$. Moreover, the failure probability decays exponentially in the number of channels.

In the next sections we provide further examples when we discuss particular choices of the matrix $A$.

VI Comparison with Multichannel Greedy Algorithms

In $p$-thresholding, we select a set $S$ of $k$ indices whose $p$-correlation with $Y$ are among the $k$ largest:

After the support $S$ is determined, the non-zero coefficients of $\hat{X}$ are computed via an orthogonal projection: $\hat{X}^{S}=A_{S}^{\dagger}Y$.

Using the probability model (21) average-case recovery theorems for $p$-thresholding and $p$-SOMP have been proven in [25, 24, Theorems 4,6,7,8]. We improve slightly on these in the following. (Note, however, that also treats the noisy case.) Our first result generalizes the one in to the multichannel setup.

Let $A\in{\hbox{\msbm{C}}}^{n\times N}$ have unit norm columns and local $2$-coherence function $\mu_{2}(S)$ defined in (30). Let $X\in{\hbox{\msbm{R}}}^{N\times L}$ with $\operatorname{supp}X\subset S$ where $S\subset\{1,\ldots,N\}$, and such that the coefficients on $S$ are given by (21), $X^{S}=\Sigma\Phi$, where we choose the real spherical model for $\Phi$. Set $Y=AX$ and $R=\max_{j}\sigma_{j}/\min_{j}\sigma_{j}$. If

then the probability that $2$-thresholding applied to $Y$ fails to recover $X$ is bounded by

If we use the complex spherical model instead of the real spherical model then $L/2$ in the above probability estimate may be replaced by $L$.

We proceed similarly as in . We denote by $\Theta$ the event that $2$-thresholding fails. Clearly,

where $\rho$ will be specified later. Denote by $Z_{j}$, $j\in S$, a sequence of independent random vectors which are uniformly distributed on the unit sphere of ${\hbox{\msbm{R}}}^{L}$. Then,

Choosing $\rho=\sigma_{\min}/2$ and applying Theorem IV.2 we obtain

where we used the definition of $\theta$ and $\mu_{2}(S)$. Similarly we estimate

Combining the two estimates completes the proof for the real case. Choosing the vectors $Z_{j}$, $j\in S$, from the complex unit sphere $S_{\hbox{\msbm{C}}}^{L}$ and using Corollary IV.3 yields the statement for the complex case. ∎

We now state the corresponding result for $2$-SOMP, which slightly improves the one in for the noiseless case. (Note that we restrict to $p=2$ here, although the theorem is easily extended to general values of $p$.)

Let $A$ be a matrix with unit norm columns and constants $\delta(S),\mu_{2}(S)<1$ where $S\subset\{1,\ldots,N\}$. Assume that

for some $\epsilon\in(0,1)$. Let $X$ be a random coefficient matrix with support $S$ that is selected according to the real Gaussian probability model, see (21), and let $Y=AX$. Then $2$-SOMP applied to $Y$ recovers $X$ in $k$ steps with probability at least

where $A_{L}\sim\sqrt{L}$ is given by (25).

If we use the complex Gaussian model instead of the real Gaussian model then the same conclusion holds with $L$ replaced by $2L$ in (42).

Due to the factor $2^{k}$ the probability bound (42) becomes effective only when the number of channels becomes comparable to the sparsity $k$. This drawback is very likely due to the analysis and is not observed in practice. However, it seems to be very difficult to remove this factor by a more sophisticated proof technique.

We require $\epsilon<1$, so that the probability decay of (42) is potentially slower than that given by Theorem IV.4.

With $\delta=\epsilon=1/2$ condition (41) is satisfied if $\mu_{2}(\Lambda)\leq 1/7$ while the probability estimate (42) behaves like $1-N2^{k}\exp(-L/4)$.

With the estimates $\delta(S)\leq\delta^{*}(S)$ and $\mu_{2}(S)\leq\delta^{*}(S)$, (41) with $\epsilon=3/11$ is implied by

VI-B Comparison

Time-Frequency shifts of the Alltop window.

We now compare this result with the condition of Theorem VI.1 concerning thresholding. As noted in (32), $\mu_{2}(S)\leq\delta^{*}(S)$. Therefore, by Proposition V.3 we have

with probability at least $1-\epsilon$ provided

and the failure probability of thresholding is bounded by $N\exp(-L/2(\theta^{-2}-\log(\theta^{-2})-1))+\epsilon$.

Let us finally consider Theorem VI.2 for SOMP. By Proposition V.3 the condition $\delta^{*}(S)<1/3$ in Remark VI.3 is satisfied with probability at least $1-\epsilon$ provided

and the failure probability of SOMP is bounded by

with $A_{L}^{2}\sim L$ if the real Gaussian probability model is used.

VI-B2 Union of Dirac and Fourier

If $S$ is chosen at random then a much better bound (up to constants) is obtained using Theorem V.4. In our special case, however, further improvement is possible. A reformulation of a result of , see also [46, Proposition 3] shows the following. If the support $S$ consists of $k_{1}$ arbitrary elements of $\{1,\ldots,n\}$ and $k_{2}$ random elements of $\{n+1,\ldots,2n\}$ then with probability at least $1-\epsilon$ we have $\delta(S)\leq 1/2$ provided

with $c=0.25$. In particular $k\leq n/4$ and the same reasoning as in the proof of Theorem V.4 yields

To compute the performance of thresholding, note that condition (39), $2R\mu_{2}(S)\leq 2R\mu\sqrt{k}\leq\theta<1,$ is satisfied provided

Assuming that the non-zero rows of the matrix $\Phi$ in the probability model (21) on the coefficients are independent and uniformly distributed on the complex unit sphere $S_{\hbox{\msbm{C}}}^{L-1}$, the failure probability of thresholding is bounded by $N\exp(-L(\theta^{-2}-\log(\theta^{-2})-1))$.

Assuming $\delta(S)\leq\delta=1/2$ and $\mu\sqrt{k}\leq 1/7$, i.e.,

VI-B3 Time-Frequency shifts of Alltop window

As in the Fourier-Dirac case, under condition (48) and the complex probability model of Theorem VI.1, thresholding fails with probability at most $N\exp(-L(\theta^{-2}-\log(\theta^{-2})-1))$.

with a constant $C$ (which also implies (35)) we have

For the analysis of SOMP we choose $\delta=1/2$ in Theorem V.5. Assuming that the square-root in (36) is less than $\frac{9}{10}e^{-1/4}\frac{1}{2}$ is equivalent to

with an appropriate $C$, and condition (36) is satisfied. Then with probability at least $1-\epsilon$ we have $\delta^{*}(S)\leq 1/2$. Furthermore, as suggested by Remark VI.3(b) the condition $\mu_{2}(S)\leq 1/12$ is also implied by (50) since $\mu_{2}(S)\leq\sqrt{k}\mu=\sqrt{\frac{k}{n}}$. Assuming the complex Gaussian probability model on the non-zero coefficients of $X$ the failure probability of SOMP is bounded by $N2^{k}\exp(-A_{2L}^{2}/2)+\epsilon$ due to Theorem VI.2.

VII Numerical Simulations

$\Phi$ is chosen to be a real Gaussian random matrix (i.e., all entries independent and standard normally distributed); $\Sigma$ has independent diagonal entries with standard normal distribution.

$\Phi$ is chosen to be a complex Gaussian random matrix (i.e., the real and imaginary parts of each entry are chosen independently according to a standard normal distribution); $\Sigma$ is equal to the identity.

In the following figures the results of various simulation runs are plotted (we always used $100$ simulations for each choice of parameters).

In Fig. 1 we plot the results when choosing $A$ from a random spherical ensemble of size $n=32$ columns and $N=256$ rows for $L=1,2,4$. The matrix $X$ was generated according to model (1). The improvement with increasing $L$ is clearly evident.

Finally, in Fig. 3 we plot the results when using time-frequency shifts of the Alltop window with $n=29$ and $N=29^{2}=841$. Here the results of thresholding are extremely poor and therefore not plotted.

VIII Conclusion

Appendix A Proof of Theorem IV.2

The proof uses the following extension of Khintchine’s inequality to higher dimensions stated in ,

for all $p\geq 2$ and all vectors $a\in{\hbox{\msbm{R}}}^{k}$. By splitting in real and imaginary parts it easily follows that this inequality also holds for all $a\in{\hbox{\msbm{C}}}^{k}$. We may assume without loss of generality that $\|a\|_{2}=1$. Then an application of Markov’s inequality yields

where $(a)_{i}=a(a+1)(a+2)\cdots(a+i-1)$ denotes the Pochhammer symbol. The last equation is due to the fact that $\sum_{i=0}^{\infty}\frac{(a)_{i}}{i!}\lambda^{i}$ is the Taylor series of $(1-\lambda)^{-a}$, which converges for $\lambda<1$. Minimizing (51) with respect to $\lambda$ gives $\lambda=1-u^{-2}$. Inserting this value yields the statement of the theorem.

Appendix B Proof of Proposition V.3

Now consider a random matrix $\Psi\in{\hbox{\msbm{R}}}^{n\times N}$ with independent columns that are uniformly distributed on the sphere $S^{n-1}$. Then $\Psi$ has the same distribution as $DA$, where $A$ is Gaussian matrix as above, $D=\operatorname{diag}(s_{1}^{-1},\ldots,s_{N}^{-1})$ and $s_{j}=\sqrt{n}\|\Phi_{j}\|_{2}$ where $\Phi_{j}\in{\hbox{\msbm{R}}}^{n}$ is a vector of independent standard normally-distributed random variables. We now use the following measure concentration inequality [3, Corollary (2.3)] or [4, eq. (2.6)] for a standard Gaussian vector $Z\in{\hbox{\msbm{R}}}^{n}$,

Appendix C Proof of Theorem VI.2

We assume that until a certain step SOMP has selected only correct indices, collected in $J\subset S$. Let us first estimate the probability that it selects a correct element of $S\setminus J$ also in the next step.

We denote by $P_{J}=A_{J}A_{J}^{\dagger}$ the orthogonal projection onto the span of the columns of $A$ in $J$, and $Q_{J}=I-P_{J}$. The residual at the current iteration is given by $Y_{M}=Q_{J}Y=Q_{J}A_{S}X=Q_{J}A_{S}\Sigma\Phi$. SOMP selects a correct index in $S\setminus J$ in the next step if

By Theorem 11 in (which is proven using Theorem IV.6; note that there is a slight error in in the computation of the constant $A_{L}$) we have the following concentration of measure inequalities

where $A_{L}$ is the constant in (25) and $C_{2}(L)={\hbox{\msbm{E}}}\|Z\|_{2}$ with $Z=(Z_{1},\ldots,Z_{L})$ being a vector of independent standard normal variables. Now we assume that

Then by the above and a union bound the probability that SOMP fails can be bounded by

where we used the fact that $A_{S\setminus J}^{*}A_{J}$ is a submatrix of $A_{S}^{*}A_{S}-I$.

Next we consider the maximum on the left hand side of (54). We can estimate

Combining the above estimates, condition (54) is satisfied if

In order to complete the proof, we note that OMP successfully recovers the correct signal if (54) holds for all $J\subset S$. By a union bound of (55) over all those $2^{k}$ subsets this is true with probability at least $1-N2^{k}\exp(-\epsilon^{2}A_{L}^{2})$ provided condition (41) holds.

The extension to the complex valued case is straightforward.

References