1

On the Support Recovery of Jointly Sparse
Gaussian Sources using Sparse Bayesian Learning
Saurabh Khanna, Student Member, IEEE and Chandra R. Murthy, Senior Member, IEEE

arXiv:1703.04930v3 [] 26 Oct 2017

Dept. of ECE, Indian Institute of Science, Bangalore, India
{sakhanna, cmurthy}@iisc.ac.in
Abstract—In this work, we provide non-asymptotic, probabilistic guarantees for successful sparse support recovery by
the multiple sparse Bayesian learning (M-SBL) algorithm in
the multiple measurement vector (MMV) framework. For joint
sparse Gaussian sources, we show that M-SBL perfectly recovers
their common nonzero support with arbitrarily high probability
using only finitely many MMVs. In fact, the support error
probability decays exponentially fast with the number of MMVs,
with the decay rate depending on the restricted isometry property
of the self Khatri-Rao product of the measurement matrix.
Our analysis theoretically confirms that M-SBL is capable of
recovering supports of size as high as O(m2 ), where m is the
number of measurements per sparse vector. In contrast, popular
MMV algorithms in compressed sensing such as simultaneous
orthogonal matching pursuit and row-LASSO can recover only
O(m) sized supports. In the special case of noiseless measurements, we show that a single MMV suffices for perfect recovery
of the k-sparse support in M-SBL, provided any k + 1 columns
of the measurement matrix are linearly independent. Unlike
existing support recovery guarantees for M-SBL, our sufficient
conditions are non-asymptotic in nature, and do not require the
orthogonality of the nonzero rows of the joint sparse signals.
Index Terms—Sparse Signal Processing, Sparse Bayesian
Learning, Joint Sparsity, Restricted Isometry Property, KhatriRao product, Statistical Learning Theory.

I. I NTRODUCTION
In recent years, joint sparsity has emerged as one of the most
important and versatile signal structures in the field of sparse
signal processing. Two or more vectors are said to be jointly
sparse if their nonzero coefficients belong to the same index
set, i.e., they share a common nonzero support. Joint sparsity
arises naturally in multi-modal or multi-channel analysis of
signals residing in low dimensional signal subspaces. The underlying joint sparsity can be exploited to resolve ambiguities
which may arise due to erroneous estimation of the support of
the individual sparse vectors from noisy measurements. This
idea has been exploited in several practical scenarios such as
MIMO channel estimation [1]–[3], distributed source coding
[4], [5], multi-task compressive sensing [6], distributed event
localization [7], array signal processing [8], and cooperative
spectrum sensing [9]–[11].
In the sparse signal recovery literature, the estimation of
jointly sparse signals is referred to as the multiple measurement
vector (MMV) problem [12] where the signal of interest is a
matrix X ∈ Rn×L whose columns are jointly sparse vectors
in Rn . As a result, X is a row sparse matrix with only a
fraction of its rows containing nonzero elements and the rest
of the rows made up entirely of zeros. In the MMV problem,
the goal is to recover X from its noisy, linear measurements

Y ∈ Rm×L . The measurement matrix Y (each column is
called a single measurement vector (SMV)) is generated as
Y = AX + W,

(1)

where A ∈ Rm×n is a known measurement matrix and W ∈
Rm×L is the unknown noise matrix. For m < n, the above
linear system is under-determined and therefore has infinitely
many solutions for X. However, if A satisfies certain restricted
isometry properties, a unique row-sparse solution can still be
guaranteed [12]–[15].
In many applications, the performance of a joint sparse
signal recovery algorithm is judged on the basis of how
accurately it can identify the true support or the locations of
the nonzero rows of X. This gives rise to the joint sparse
support recovery (JSSR) problem where the goal is to recover
the row support of X given Y and A. Interestingly, unlike the
nonzero coefficients in a k-sparse X which can be recovered
only if m ≥ k, the nonzero support can be recovered even
from m < k measurements. In fact, for i.i.d. Gaussian entries
in both A and X, a non-iterative, correlation based algorithm
called One Step Greedy Algorithm (OSGA) [16] is capable of
recovering the true support using only m ≥ 1 measurement
per signal with probability approaching one as L → ∞.
A majority of the existing MMV algorithms [12], [14], [17]–
[21] implicitly assume that the number of nonzero rows in
X is less than the number of measurements per signal, i.e,
k < m. Currently, M-SBL [22], Co-LASSO [23] and RDCMP [24] are the only MMV algorithms which have been
theoretically or empirically shown to be capable of recovering
a k-sparse support from fewer than k measurements per signal.
However, a theoretical understanding of when these algorithms
are guaranteed to be successful is far from complete. In
this work, we derive non-asymptotic guarantees for M-SBL’s
performance in the JSSR problem. Specifically, we derive the
sufficient conditions for exact support recovery in terms of the
number of MMVs and properties of A.
M-SBL’s support recovery performance has been investigated in [25] for k < m, and in [22] for k ≥ m. In both studies, the orthogonality of the nonzero rows of X is identified
as one of the sufficient conditions for exact support recovery.
For finite L, the row-orthogonality condition is too restrictive
for a deterministic X and almost never true for a continuously
distributed random source. Thus, the existing support recovery
guarantees for M-SBL are in reality only applicable in the
asymptotic sense when L → ∞. Furthermore, the earlier analysis is restricted only to noiseless measurements. In contrast,
the new sufficient conditions derived in this work are nonasymptotic in nature and do not assume row orthogonality

2

in the signal matrix X. Our analysis also accounts for the
presence of measurement noise. While our focus here is on
sufficient conditions for exact support recovery, the necessary
conditions have been partly covered in [26] and [27].
A. Existing Theoretical Guarantees for Support Recovery
The earliest theoretical work focused on seeking guarantees
for a unique joint sparse solution to the canonical ℓ0 norm
minimization problem:
L0 :

min

X∈Rm×L

||X||0

s.t. AX = Y,

(2)

i.i.d. N (0, diag(1S )), where S denotes the unknown

 support

set. For m = Ω k log nk , it is shown that L ≥ O logloglogn n
suffices for diminishing support error probability with increasing L. One of our contributions in this paper is to extend this
result to a more general signal prior on X and show that the
support error probability vanishes even if m scales sublinearly
in the support size k.
In [32], the support recovery problem is analyzed as a
single-input-multi-output MAC communication problem. For
n
number of nonzero rows fixed to k, m = log
c(X) is shown to be
both necessary and sufficient for successful support recovery
as the problem size tends to infinity. The quantity c(X) is
a capacity like term which depends on the elements of the
nonzero rows in X and the noise power. Even fewer measurements m = Ω( Lk log n) suffices when each measurement
vector is generated using a different measurement matrix [33].
In the above discussion, in particular about algorithmspecific support recovery guarantees, the number of nonzero
rows in X is typically assumed to be less than m, the number
of measurements per MMV. In the following subsection, we
review the covariance matching framework, which is capable
of recovering supports of size k greater than m.

where ||X||0 denotes the number of nonzero rows in X.
In [12], [13], the authors showed that for spark1 (A) = m + 1
and rank(Y) ≤ m, the L0 problem admits a unique k-sparse
solution if k < ⌈(spark(A) − 1 + rank(Y)) /2⌉. This result
confirmed that the SMV bottleneck of k < m/2 for ℓ0 norm
based support recovery is not applicable when multiple measurement vectors are used. Furthermore, the sparsity bound
suggests that supports of size k < m are potentially uniquely
recoverable. Since this result, numerous JSSR algorithms have
been proposed in the quest to meet the k < m sparsity bound.
To circumvent the combinatorial hardness of the L0 problem, [17] proposed minimizing the mixed-norm ℓp,q norm
of X instead of the ℓ0 norm. Variants of the ℓp,q norm
minimization problem with different combinations of p and
q have been investigated independently in several works [12],
[13], [28], [29]. For p ≥ 1, q = 1, [13] has shown that ℓp,q
norm minimization problem has a unique k-sparse solution,
< 1, for all j ∈
/ S and for
provided A satisfies A†S aj
1
−1 T
†
T
all S ⊂ [n], |S| ≤ k, where AS = AS AS
AS . This
also serves as a sufficient condition for exact support recovery
in Simultaneous Orthogonal Matching Pursuit (SOMP) [14],
a greedy support reconstruction algorithm. In [18], support
recovery performance of various correlation based greedy and
iterative hard-thresholding type algorithms is studied in the
noiseless MMV setup. The sufficient conditions for exact
support recovery are specified in terms of the asymmetric
restricted isometry constants of the measurement matrix A.
A limitation of the aforementioned support recovery techniques and associated guarantees is that only supports of size
up to k < m/2 are uniquely recoverable. In [30], rank aware
OMP and rank aware order recursive matching pursuit are
proposed and analyzed, and it is shown that both these algorithms perfectly recover any k-sparse support from noiseless
measurements as long as k < spark(A) − 1 and rank(X) = k.
For the rank defective case, i.e., rank(X) < k, compressed
sensing MUSIC [20] and subspace-augmented MUSIC [19]
are still capable of recovering any k < spark(A) − 1 sized
support as long as k−rank(X) partial support can be estimated
by another sparse signal recovery algorithm.
In [31], the support recovery problem is formulated as a
multiple hypothesis testing problem. Necessary and sufficient
conditions for perfect support recovery with high probability
are derived under the assumption that the columns of X are

In (3), the linear constraints are the vector form of the
second order moment constraints, i.e., the covariance matching
equation: L1 YYT = Adiag(γ)AT + σ 2 Im , where σ 2 Im
denotes the noise covariance matrix. Since the constraints in
(3) comprise up to (m2 +m)/2 linearly independent equations,
sparsity levels as high as O(m2 ) are potentially recoverable.
To recover the maximum level of sparsity, k = (m2 + m)/2,
a necessary condition derived in [23], [27] is that the columnwise self Khatri-Rao product matrix A ⊙ A must have full
Kruskal rank,2 i.e., Krank(A ⊙ A) = (m2 + m)/2.
The M-SBL algorithm [25], our focus in this paper, also
imposes a common Gaussian prior N (0, diag(γ)) on the
columns of X and hence implicitly exploits the latent uncorrelatedness of the nonzero entries in X. Interestingly, similar
to Co-LASSO, the support recovery performance of M-SBL is
closely related to the properties of the self Khatri-Rao product

1 Spark of a matrix A is the smallest integer p such that there exist p
linearly dependent columns in A.

2 The Kruskal rank of an m × n matrix A is the largest integer k such that
any k columns of A are linearly independent.

B. Covariance Matching Based Support Recovery
A key insight was propounded in [23], that there often exists
a latent structure in the MMV problem: the nonzero entries
of X are uncorrelated. This signal structure can be enforced
by assuming each column of X to be i.i.d. N (0, diag(γ)),
where γ ∈ Rn+ is a non-negative vector of variance parameters.
Under this source model, identifying the nonzero rows of X
is tantamount to estimating the support of γ. In [23], the
Co-LASSO algorithm was proposed for the recovery of γ.
Instead of working directly with the linear observations, CoLASSO uses their covariance form, L1 YYT , as input, and γ
is recovered by solving the following non-negative ℓ1 norm
minimization problem:


1
YYT − σ 2 Im . (3)
minn ||γ||1 s.t. (A ⊙ A)γ = vec
γ∈R+
L

3

A ⊙ A. Making this relation explicit is one of the main
contributions of this paper.
C. Contributions
Following are the main contributions of this work.
1) A new interpretation of the M-SBL algorithm as a
Bregman matrix divergence minimization problem, which
opens up new avenues to exploit the vast literature on
Bregman divergence minimization for devising faster,
more robust algorithms towards minimizing the M-SBL
cost function.
2) New sufficient conditions under which M-SBL exactly
recovers the true support of the joint sparse vectors
with high probability in the noisy MMV problem. The
conditions are specified in terms of the number of MMVs
and the properties of the measurement matrix.
3) The support error probability in M-SBL is theoretically
shown to decay exponentially with the number of MMVs.
The error exponent is related to the restricted isometry
property of A ⊙ A, the self Khatri-Rao product of A.
Explicit upper bounds on the number of MMVs sufficient
for vanishing support error probability in M-SBL for both
noisy as well as noiseless measurements are derived.
A key feature of our analysis is that our sufficient conditions
are expressed in terms of the restricted isometry property (RIP)
constants of the measurement matrix A and that of its selfKhatri-Rao product A ⊙ A. This makes our results applicable
to both random as well as deterministic constructions of A.
The remainder of the paper is organized as follows. In
section II, we formulate the JSSR problem and introduce our
signal model for X. We also review the M-SBL algorithm
[25] and interpret the M-SBL cost function as a Bregman
matrix divergence. In section III, we cover some preliminary
concepts which are used while analyzing the support recovery
performance of M-SBL. In section IV, we derive an abstract
upper bound for the support error probability, which is used
in section V to derive our main result, namely, the sufficient
conditions for vanishing support error probability in M-SBL.
In section VI, we discuss the implications of the new results
in the context of several interesting special cases. Our final
conclusions are presented in section VII.
D. Notation
Throughout this paper, scalar variables are denoted by
lowercase alphabets and vectors are denoted by boldface lowercase alphabets. Matrices are denoted by uppercase boldface
alphabets and uppercase calligraphic alphabets denote sets.
Given a vector x, x(i) represents its ith entry. supp(x)
denotes the support of x, the set of indices corresponding
to nonzero entries in x. Likewise, R(X) denotes the set of
indices of all nonzero rows in X and is called the rowsupport of X. For any n ∈ N, [n] , {1, 2, . . . , N }. For any
n dimensional vector x and index set S ⊆ [n], the vector
xS is an |S| × 1 sized vector retaining only those entries
of x which are indexed by elements of S. Likewise, AS
is a submatrix comprising the columns of A indexed by S.

Null(A) and Col(A) denote the null space and column space
of the matrix A, respectively. The spectral, Frobenius and
maximum absolute row sum matrix norms of A are denoted
by |||A|||2 , ||A||F , and |||A|||∞ , respectively. P(E) denotes
the probability of event E. N (µ, Σ) denotes the Gaussian
probability density with mean µ and covariance matrix Σ.
For any square matrix C, tr(C) and |C| denote its trace and
n
determinant, respectively. Lastly, S++
denotes the set of all
n × n positive definite matrices.
II. S YSTEM M ODEL AND M-SBL A LGORITHM
A. Joint Sparse Support Recovery (JSSR)
Suppose x1 , x2 , . . . , xL are L distinct joint-sparse vectors
in Rn with a common nonzero support denoted by the index
set S ∗ ⊆ [n]. Let K be the maximum size of the common
support, i.e., |S ∗ | ≤ K. In JSSR, we are interested in recovering S ∗ from noisy underdetermined linear measurements
y1 , y2 , . . . , yL generated as
yj = Axj + wj , 1 ≤ j ≤ L.

(4)

The measurement matrix A ∈ Rm×n is assumed to be a
non-degenerate matrix, with m ≤ n. By non-degeneracy
of A, it is implied that any m columns of A are linearly
independent, or spark(A) = m + 1. The noise vector w ∈ Rm
is assumed to be zero mean Gaussian distributed with diagonal
covariance matrix σ 2 Im . The measurement matrix A and
the noise variance σ 2 are assumed to be known. The linear
measurement model in (4) can be rewritten in a compact
MMV form as Y = AX + W, where Y = [y1 , y2 , . . . , yL ],
X = [x1 , x2 , . . . , xL ] and W = [w1 , w2 , . . . , wL ] are the
observation, signal and noise matrices, respectively. Since the
columns of X are jointly sparse, with common support S ∗ , X
is a row sparse matrix with row support R(X) = S ∗ .
B. Source Signal Model
We assume that if the ith row of the unknown signal matrix
X nonzero, then it is a Gaussian ensemble of L i.i.d. zero mean
random variables with a common variance γ ∗ (i) belonging
to the interval [γ min , γ max ]. In the sequel, we refer to this as
Assumption (A1). An immediate consequence of (A1) is that
there exists a bounded, nonnegative, and at most K sparse
vector, γ ∗ ∈ Rn+ , such that the columns xj are i.i.d. N (0, Γ∗ )
with Γ∗ , diag(γ ∗ ). Furthermore, R(X) and supp(γ ∗ ) are
the same and equal to S ∗ .
C. M-SBL Algorithm
In this section, we review the M-SBL algorithm [25], a
type-II maximum likelihood (ML) procedure for estimation
of joint sparse signals from compressive linear measurements. In M-SBL, the columns of X are assumed to be
i.i.d. N (0, Γ) distributed, where Γ = diag(γ), and γ =
T
[γ(1), γ(2), . . . , γ(n)] is an n × 1 vector of non-negative
variance parameters. The elements of γ are collectively called
hyperparameters as they represent the parameters of the signal
prior. Since the hyperparameter γ(i) models the common
variance of the ith row of X, if γ(i) = 0, it drives the posterior

4

estimate of xj (i) to zero for 1 ≤ j ≤ L. Consequently, if γ
is a sparse vector, it induces joint sparsity in X.
In M-SBL, the hyperparameter vector γ is chosen to maximize the Bayesian evidence p(Y; γ), which is tantamount to
finding the maximum likelihood (ML) estimate of γ. Let γ̂ ML
denote the ML estimate of γ, i.e.,
γ̂ ML = arg max log p(Y; γ).

(5)

γ∈Rn
+

The Gaussian prior on xj combined with the linear measurement model induces Gaussian observations, i.e., p(yj ; γ) =
N (0, σ 2 Im + AΓAT ). For a fixed γ, the MMVs yj are
mutually independent. Hence, it follows that
log p(Y; γ) =

L
X

log p(yj ; γ)

j=1

2


T
,
∝ −L log |Σγ | − Tr Σ−1
γ YY

(6)

T

where Σγ = σ Im + AΓA . The log likelihood log p(Y; γ)
in (6) is a non convex function of γ and its global maximizer
γ̂ ML is not available in closed form. However, its local
maximizer can still be found via fixed point iterations or the
Expectation-Maximization (EM) procedure.
In [25], it is empirically shown that the EM procedure
faithfully recovers the true support S ∗ , provided m and L
are sufficiently large. In this paper, we derive new sufficient
conditions for supp (γ̂ ML ) = S ∗ , i.e., for perfect support
recovery via maximization of M-SBL’s log-likelihood cost.

Σγ , with point wise distances measured using the Bregman
LogDet divergence. Thus, the M-SBL algorithm, at its core,
is essentially a second order moment matching or covariance
matching procedure which selects γ such that the associated
covariance matrix Σγ is closest to the sample covariance
matrix, in the Bregman LogDet divergence sense. In a later
section, we theoretically show that if the second moment
matching equations are too ill-conditioned, then M-SBL fails
to recover the true support of the joint sparse columns of X.
This new interpretation of the M-SBL cost as a Bregman
matrix divergence beckons two interesting questions:
i Are there other matrix divergences besides LogDet Bregman matrix divergence which are better suited for covariance matching?
ii How to exploit the structural similarities between the MSBL cost and the Bregman (LogDet) matrix divergence
to devise faster and more robust techniques for the type-2
maximum likelihood procedure?
We believe that evaluating the performance of other matrix
divergences for covariance matching is a worthwhile exercise
to pursue and can lead to development of new, improved
algorithms for the JSSR problem. Our preliminary results in
this direction can be found in [24].
III. P RELIMINARY C ONCEPTS
In this section, we introduce a few key definitions and
results which are used in the later sections.

D. The M-SBL cost is a Bregman Matrix Divergence
We now present an interesting interpretation of M-SBL’s
log-marginalized likelihood cost in (6) which facilitates a
deeper understanding of what is accomplished by its maximization. We begin by introducing the Bregman matrix divergence Dϕ (X, Y) between any two n × n positive definite
matrices X and Y as
Dϕ (X, Y) , ϕ(X) − ϕ(Y) − h∇ϕ(Y), X − Yi,

(7)

n
where ϕ : S++
→ R is a convex function with ∇ϕ(Y)
as its first order derivative evaluated at point Y. In (7), the
matrix inner product hX, Yi is evaluated as tr XYT . For the
specific case of ϕ(·) = − log | · |, a strongly convex function,
we obtain the Bregman LogDet matrix divergence given by

(8)
Dlogdet (X, Y) = tr XY−1 − log XY−1 − n.

By termwise comparison of (6) and (8), we observe that the
negative log likelihood − log p(Y; γ) and Dlogdet (RY , Σγ ) are
the same up to a constant. In fact, it is shown in [34], [35] that
every regular exponential family of probability distributions is
associated with a unique Bregman divergence.
In the divergence term Dlogdet (RY , Σγ ), the first argument
RY , L1 YYT is the sample covariance matrix of the
observations Y and the second argument Σγ = σ 2 I + AΓAT
is the parameterized covariance matrix of Y. This connection between M-SBL’s log likelihood cost and the LogDet
divergence reveals that by maximizing the M-SBL cost, we
seek a γ which minimizes the distance between RY and

A. Restricted Isometry Property
A matrix A ∈ Rm×n is said to satisfy the restricted
isometry property (RIP) of order k if there exists a constant
δkA ∈ (0, 1) such that
2

2

2

(1 − δkA ) ||x||2 ≤ ||Ax||2 ≤ (1 + δkA ) ||x||2

(9)

holds for any k-sparse vector x ∈ Rn . The smallest such δkA
is called the k th order restricted isometry constant (RIC) of A.
B. ǫ-Cover, ǫ-Net and Covering Number
Suppose T is a set equipped with a pseudo-metric d. For
any set A ⊆ T , its ǫ-cover is defined as the coverage of A with
open balls of radius ǫ and centers in T . The set Aǫ comprising
the centers of these covering balls is called an ǫ-net of A. The
minimum number of ǫ-balls which can cover A is called the
ǫ-covering number of A, and is given by
ǫ
Ncov
(A, d) = min {|Aǫ | : Aǫ is an ǫ-net of A}

(10)

In computational theory of learning, ǫ-net constructs are often
useful in converting a union over the elements of a continuous
set to a finite sized union.
Proposition 1 ( [36]). Let B(0, 1) be a unit ball in Rn centered
at 0. Then, its ǫ-covering number with respect to the standard
n
ǫ
Euclidean metric is bounded as Ncov
(B(0, 1), ||·||2 ) ≤ (3/ǫ) .

5

C. Strong Convexity of − log det
A differentiable function f is called strongly convex with
strong convexity constant mf > 0 if the following inequality
holds for all points x, y in the domain.
mf
2
f (y) ≥ f (x) + h∇f (x), y − xi +
||y − x||2 ,
(11)
2
where ||·||2 is the Frobenius norm.

m
Proposition 2. Let MK ⊆ S++
be the collection of ma2
trices of form Σγ = σ I + AΓAT , where Γ = diag(γ),
γ ∈ Rn+ , ||γ||0 ≤ K, and A ∈ Rm×n . The matrix function ψ(·) = − log det(·) is strongly convex in the domain
MK with strong convexity constant mψ given by mψ =
−1
supγ∈Rn+ ,||γ||0 ≤2K |||Σγ ◦ Σγ |||2
, where ◦ denotes the
Hadamard product.

Proof. See Appendix A.
D. Rényi Divergence
Let (X , F ) be a measurable space and P and Q be two
probability measures on F with densities p and q, respectively,
with respect to the dominating Lebesgue measure µ on F .
For α ∈ R+ \1, the Rényi divergence of order α between P
and Q, denoted Dα (p||q), is defined as
Z
1
p(x)α q(x)1−α µ(dx).
(12)
log
Dα (p||q) =
α−1
X

For p = N (0, Σ1 ) and q = N (0, Σ2 ), the α-Rényi
divergence Dα (p||q) is available in closed form [37],
Dα (p||q) =

1
|(1 − α)Σ1 + αΣ2 |
.
log
1−α
α
2(1 − α)
|Σ1 |
|Σ2 |

(13)

E. A Lower Bound on the α-Rényi Divergence Between Multivariate Gaussian Distributions
The strong convexity of − log det(·) can be used to derive a
lower bound on the Rényi divergence between two multivariate
Gaussian distributions, stated below as Theorem 1.
Theorem 1. Let p1 and p2 be two multivariate Gaussian
distributions with zero mean and positive definite covariance
matrices Σ1 and Σ2 , respectively. Then, for α ∈ [0, 1), the
α-Rényi divergence between p1 and p2 satisfies
1
2
(14)
Dα (p1 , p2 ) ≥ αm∗ ||Σ2 − Σ1 ||F ,
4
where m∗ is the strong convexity constant of the matrix
function ψ(Σ) = − log |Σ| over the domain Θ(Σ1 , Σ2 ), the
collection of all convex combinations of Σ1 and Σ2 .
Proof. See Appendix B.

We now proceed with the analysis of support recovery performance of M-SBL algorithm. We derive an M-SBL specific
Chernoff bound for the support error probability in the JSSR
problem under assumption (A1). We begin by introducing
some of the frequently used notation in the following table.
K
Maximum number of nonzero rows in X
Sk
Collection of all support sets of k or lesser
size, i.e., Sk = {S ⊆ [n], |S| ≤ k}
Θ(S) A bounded hyperparameter set associated with
the
set S, formally defined as Θ(S) ,
 support
γ ∈ Rn+ : supp(γ) = S, γ min  γ S  γ max .
ΘK
Collection of all nonnegative K or less sparse
vectors in Rn+ whose nonzeroSelements belong
to [γ min , γ max ]. Also, ΘK = S∈SK Θ(S).
S∗
True row support of X.
γ∗
Principal diagonal of the common covariance
matrix Γ∗ of the i.i.d. columns in X.
Also, supp(γ ∗ ) = S ∗ .

By assumption (A1) on X, we have γ ∗ ∈ ΘK . Thus, in
order to recover γ ∗ from Y, we consider solving a constrained
version of the M-SBL optimization problem in (5):

C
ǫ2

log 2δ .

γ̂ = arg max L(Y; γ).

cM-SBL(K):

(16)

γ∈ΘK

In (16), the objective L(Y; γ) , log p(Y; γ) is the same as
the M-SBL log-likelihood cost in (6). The row support of X
is then estimated as supp(γ̂), where γ̂ is the solution of (16).
We define the set of bad MMVs,

(17)
ES ∗ , Y ∈ Rm×L : supp (γ̂) 6= S ∗ ,

which result in erroneous estimation of the true support. We
are interested in identifying the conditions under which P(ES ∗ )
can be made arbitrarily small.
Note that ES ∗ is the collection of undesired MMVs for
which the log-likelihood objective is globally maximized by
γ ∈ Θ(S), for S 6= S ∗ . Based on this interpretation, ES ∗ can
be rewritten as


[
Y : max L(Y; γ) ≥ ′ max ∗ L(Y; γ ′ ) .
ES ∗ =
S∈SK \{S ∗ }

Since

γ∈Θ(S)

γ ∈Θ(S )

(18)

max ∗ L(Y; γ ′ ) ≥ L(Y; γ ∗ ), we have
′

γ ∈Θ(S )

ES ∗ ⊆

[

S∈SK \{S ∗ }

=

F. Concentration of Sample Covariance Matrix
Proposition 3 (Vershynin [38]). Let y1 , y2 , . . . , yL ∈ Rm
beP
L independent samples from N (0, Σ), and let ΣL =
L
1
T
j=1 yj yj denote the sample covariance matrix. Then, for
L
any ǫ > 0,
|||ΣL − Σ|||2 ≤ ǫ |||Σ|||2 ,
(15)
holds with probability exceeding 1 − δ provided L ≥
Here, C is an absolute constant.

IV. S UPPORT E RROR P ROBABILITY A NALYSIS

[

S∈SK

\{S ∗ }



Y : max L(Y; γ) ≥ L(Y; γ ∗ )
γ∈Θ(S)

[

{Y : L(Y; γ) − L(Y; γ ∗ ) ≥ 0} .(19)

γ∈Θ(S)

The continuous union over infinitely many elements of Θ(S)
in (19) can be relaxed to a finite sized union by using the
following ǫ-net argument.
Consider Θǫ (S), a finite sized ǫ-net of the hyperparameter
set Θ(S), such that for any γ ∈ Θ(S), there exists an element
γ ′ ∈ Θǫ (S) so that |L(Y; γ) − L(Y; γ ′ )| ≤ ǫ. Proposition 4
gives an upper bound on the size of such an ǫ-net.

6

Proposition 4. Given a support set S ⊆ [n], there exists a
finite set Θǫ (S) ⊂ Θ(S) such that it simultaneously satisfies

(i) For any γ ∈ Θ(S), there exists a γ ′ ∈ Θǫ (S) such that
|L(Y; γ) − L(Y;γ ′ )| ≤ ǫ.

|S| 
p
ǫ
(ii) |Θ (S)| ≤ max 1, 3CL,S (γ max − γ min ) |S|/ǫ
,

where CL,S is the Lipschitz constant of L(Y; γ) with
respect to γ in the bounded domain Θ(S).
The set Θǫ (S) is an ǫ-net of Θ(S).
Proof. See Appendix C.

From Proposition 4-(ii), we observe that the construction
of Θǫ (S) depends on the Lipschitz continuity of the loglikelihood L(Y; γ) with respect to γ over the domain Θ(S).
By virtue of the data-dependent nature of L(Y; γ), its Lipschitz constant CL,S depends on the instantaneous value of
Y. To make the rest of the analysis independent of Y, we
introduce a new MMV set G such that CL,S can be uniformly
bounded solely in terms of the statistics of Y. A possible
choice of G could be




1
T
T
m×L
≤ 2E y1 y1
. (20)
YY
G, Y⊂R
:
L
2

For Y ∈ G, L(Y; γ) is uniformly Lipschitz continuous
irrespective of Y, and hence the ǫ-net Θǫ (S) can now be
constructed entirely independent of Y.
Since for arbitrary sets A and B, A ⊆ (A ∩ B) ∪ B c , the
RHS in (19) relaxes as



 [
[
{L(Y; γ) − L(Y; γ ∗ ) ≥ 0} ∩ G ∪ G c .
ES ∗ ⊆


∗
S∈SK \S γ∈Θ(S)

(21)
Let Θǫ (S)|G denote an ǫ-net of Θ(S) constructed under the
assumption that Y ∈ G. Then, the continuous union over Θ(S)
relaxes to a finite sized union over Θǫ (S)|G as shown below.


 [

[
ES ∗ ⊆
{L(Y; γ)−L(Y; γ ∗ ) ≥ −ǫ} ∩ G ∪ G c


S∈SK \S ∗ γ∈Θǫ (S)|G



 [
[
{L(Y; γ) − L(Y; γ ∗ ) ≥ −ǫ} ∪ G c .
⊆


ǫ
∗
S∈SK \S

γ∈Θ (S)|G

By applying the union bound, we obtain
X
X
P (L(Y; γ)−L(Y; γ ∗ ) ≥ −ǫ)
P (ES ∗ ) ≤ P (G c )+
S∈SK \S ∗ γ∈Θǫ (S)|G

(22)
From (22), P(ES ∗ ) will be small if each of the constituent
probabilities P (L(Y; γ) − L(Y; γ ∗ ) ≥ −ǫ), γ ∈ Θǫ (S)|G are
sufficiently small so that their collective contribution to the
double summation in (22) remains small. In Theorem 2, we
show that each event within the summation is a large deviation
event which occurs with an exponentially decaying probability.
Theorem 2. Let pγ (y) denote the marginal probability density
of the columns of Y induced by the joint sparse X with
columns drawn independently from N (0, diag(γ)). Then, the

P
log-likelihood L(Y; γ) = L
j=1 log pγ (yj ) satisfies the following large deviation property.

P (L(Y; γ) − L(Y; γ ∗ ) ≥ −ǫ) ≤ exp (−Lψ ∗ (−ǫ/L)),
(23)
where ψ ∗ (·) is the Legendre transform3 of ψ(t) , (t −
1)Dt (pγ , pγ ∗ ), and Dt is the t-Rényi divergence (of order
t > 0) between the probability densities pγ and pγ ∗ .
Proof. See Appendix D.

Note that, when the measurement noise is Gaussian, the
marginal density pγ (yj ) of the individual observations is
also Gaussian with zero mean and covariance matrix Σγ =
σ 2 Im +AΓAT . If σ 2 > 0, both marginals pγ and pγ ∗ are nondegenerate and hence the Rényi divergence term Dt (pγ , pγ ∗ )
in Theorem 2 is well defined. We now restate Theorem 2 as
Corollary 1, which is the final form of the large deviation
result for L(Y; γ) that will be used in bounding P(ES ∗ ).

Corollary 1. For an arbitrary γ ∈ Rn+ , and the true variance
parameters γ ∗ , let the associated marginal densities pγ and
pγ ∗ be as defined in Theorem 2, and suppose σ 2 > 0. Then, the
log-likelihood L(Y; γ) satisfies the large deviation property

P L(Y; γ) − L(Y; γ ∗ ) ≥ −LD1/2 (pγ , pγ ∗ )/2

≤ exp −LD1/2 (pγ , pγ ∗ )/4 . (24)

Proof. The
 large deviation result is obtained by replacing
ψ ∗ − Lǫ in Theorem 2 by its lower bound − tǫ
L − ψ (t),
followed by setting t = 1/2 and ǫ = LD1/2 (pγ , pγ ∗ )/2.

Note that, in the above, we have used the sub-optimal
choice t = 1/2 for the Chernoff parameter t, since its
optimal value is not available in closed form. However, this
suboptimal selection of t is inconsequential as it figures only
as a multiplicative constant in the final sample complexity. By
using Corollary 1 in (22), we can upper bound P (ES ∗ ) as


X
∗
κcov exp −LD1/2
(S)/4 + P (G c ) , (25)
P(ES ∗ ) ≤
S∈SK \S ∗

∗
D1/2
(S)

where
is the minimum value of D1/2 (pγ , pγ ∗ ) across
all ǫ-net points γ in Θǫ (S)|G . A new constant κcov is introduced in (25) which denotes the maximum size of the ǫ-nets
Θǫ (S)|G , S ∈ SK , i.e.,
κcov , max Θǫ (S)|G .
S∈SK

(26)

∗
Proposition 5 presents a lower bound for D1/2
(S) for S ∈ SK .

Proposition 5. Let pγ denote the parameterized multivariate
Gaussian density with zero mean and covariance matrix Σγ =
σ 2 I + AΓAT , Γ = diag(γ). For any pair γ 1 , γ 2 ∈ ΘK such
that supp(γ 1 ) differs from supp(γ 2 ) in exactly kd locations,
the α-Rényi divergence between pγ 1 and pγ 2 satisfies

2
⊙
 α
kd (1 − δ2K
)
γ min
,
Dα pγ 1 , pγ 2 ≥
T
2
2
4 σ + γ max Ka∗ sup
AS AS 2
S:|S|=2K

3 For any convex function f : X → R on a convex set X ⊆ Rn , its
Legendre transform is the function f ∗ defined by

f ∗ (z) = sup (hz, xi − f (x)) .
x∈X

7

⊙
where δ2K
is the 2K th order restricted isometry constant of
the self Khatri-Rao product A ⊙ A, and a∗ = ||vec(A)||∞ .

Proof. See Appendix E.
Using Proposition 5 with α = 1/2 in (25), we obtain



X
∗
η
log κcov
exp −LkdS,S
P(ES ∗ ) ≤
+ P (G c ) ,
−
4
L
∗
S∈SK \S

with η =

1
8



(27)

γ min
σ 2 + γ max

2

Ka2∗

⊙
1 − δ2K
ATS AS
sup

S:|S|=2K

.

(28)

2

Theorem 3. Suppose S ∗ is the true row support of the
unknown X satisfying assumption A1 and with |S ∗ | ≤ K.
Then, for any δ ∈ (0, 1), P(ES ∗ ) ≤ 2δ, if





8
2
1+δ
8
, log κcov , C log
,
log 3enK
L ≥ max
η
δ
η
δ
(29)
where κcov and η are as defined in (26) and (28), respectively, and C > 0 is a universal numerical constant.
c

Proof. Since L ≥ C log (2/δ), by Proposition 3, P(G ) ≤ δ.
Combined with L ≥ 8 log κcov /η, (27) can be rewritten as


X
∗
exp −ηLkdS,S /8 + δ.
P(ES ∗ ) ≤
(30)
S∈SK \S ∗

The total number of support sets belonging to SK \S ∗
∗
which differ
the true support
S in exactly kd locaPkd from
j+|S ∗ |
n−|S ∗ |
tions is j=0
, which is further upper
j
min(kd ,j+|S ∗ |)
k

bounded by (3enK) d . Thus, we can rewrite (30) as
P(ES ∗ )

≤ δ +

K+|S ∗ |

X

kd =1
∗

≤ δ+

K+|S |

X

X



ηLkd
exp −
8

S∈SK \S ∗ ,
|(S\S ∗ )∪(S ∗ \S)|=kd
kd

(3enK)

kd =1

 ηL kd
e− 8
.

P(ES ∗ )

≤ δ+

X

kd =1

δ
1+δ

kd

by noting that for all S ⊆ [n], |S| = 2K,
q

ATS AS 2 ≤ tr AS ATS AS ATS = vec(AS ATS ) 2
q
√
⊙
= (AS ⊙ AS )1|S| 2 ≤ 1 + δ2K
2K. (33)
B. Upper Bound for κcov
∗
By setting ǫ = LD1/2
/2 in Proposition 4, we have
 
K 
√
K
κcov ≤ max 1, 6 K(γ max − γ min )CL /ηL
,

(31)

≤ δ + δ = 2δ.

In Theorem 3, we finally have an abstract bound on the
sufficient number of MMVs, L, which guarantee vanishing
support error probability in cM-SBL(K), for a fixed true
support S ∗ . The MMV bound is meaningful only when η is
strictly positive. We now proceed to show that by imposing
certain isometry conditions on the self Khatri-Rao product
A ⊙ A, the following can be ensured.
i. η > 0 (a positive support error exponent),
ii. η and κcov scale favorably with the system dimensions.

(34)

where CLK = maxS∈SK CL,S denotes the Lipschitz constant
of L(Y; γ) with respect to γ over the bounded domain ΘK .
Proposition 6 characterizes the Lipschitz property of L(Y, γ).
Proposition 6. For S ∈ SK , the log-likelihood L(Y; γ) :
Θ(S) → R is Lipschitz continuous in γ as shown below.


|||RY |||2
LK
|L(Y,γ 2 )−L(Y,γ 1 )| ≤
1+
||γ 2 − γ 1 ||2 ,
γ min
σ2
for any γ 1 , γ 2 ∈ Θ(S). Here, RY ,

1
T
L YY .

Proof. See Appendix F.
Under 2K-RIP compliance of A ⊙ A√and for Y ∈ G
 (as
⊙ 1/2
defined in (20)), |||RY |||2 ≤ 2 σ 2 +γ max K(1+δK
)
. By
using Proposition 6 in (34), combined with the lower bound
for η in (32) yields the following upper bound for κcov .
 
K 
7/2 2
2
, (35)
κcov ≤ max 1, K a∗ ζ(γ max − γ min )/σ
√  1+δ⊙   σ2 +γ 2  3σ2 +2γ 
max
max
.
where ζ = 48 2 1−δ2K
⊙
γ
γ
2K


, P(ES ∗ ) can be upper
Since L ≥ η8 log 3enK 1+δ
δ
bounded by a geometric series as
K+|S ∗ | 

A. Lower Bound for η
From (28), it is evident that if the self Khatri-Rao product
⊙
A ⊙ A satisfies the RIP condition: δ2K
< 1, then η is strictly
positive. Furthermore, η can be lower bounded as

2
⊙
1 − δ2K
γ min
1
√
q
,
(32)
η≥
2
⊙
8 2K 3/2 a2∗ σ + γ max
1 + δ2K

min

min

V. S UFFICIENT C ONDITIONS FOR S UPPORT R ECOVERY
Equipped with explicit bounds for η and κcov , we now
state the sufficient conditions for vanishing support error
probability in cM-SBL(K) as the following theorem.

Theorem 4. Under assumption (A1), suppose X has row
support S ∗ , |S ∗ | ≤ K. Let γ̂ be the solution of cMSBL(K). Then, supp (γ̂) = S ∗ with probability exceeding
1 − 2δ, for any δ ∈ (0, 1), provided
C1. A is non-degenerate, and A ⊙ A satisfies the RIP
⊙
of order 2K, i.e., δ2K
< 1.
C2. The number of MMVs, L, satisfies
L ≥ K 3/2 a2∗ ξ max {log (3enK(1 + δ)/δ) , log κcov }
 2
2
⊙ 1/2
√
σ + γ max
(1 + δ2K
)
where ξ = 64 2
,
⊙
γ min
1 − δ2K

8

a∗ = ||vec(A)||∞ and κcov is bounded as in (35).
Proof. C1 ensures that η in (28) is strictly positive. Further,
C2 ensures that the abstract MMV bound in (29) is satisfied.
Therefore, by Theorem 3, P(ES ∗ ) ≤ 2δ.
In the following corollary, we state an extra condition
besides C1 and C2 which extends the above cM-SBL result
to the M-SBL cost function.

B. Continuous versus Binary Hyperparameters
In [31], Tang and Nehorai formulated support recovery
using MMVs as a multiple hypothesis testing problem by
assuming that each column of X is i.i.d. N (0, diag(1S ∗ )),
where S ∗ is the k-sparse support of X. For this choice of
signal prior, finding the true support via type-2 likelihood
maximization as in (5) is no longer a continuous variable
optimization but rather a combinatorial search over all k-sparse
n
vertices of the hypercube {0, 1} , as described below.
γ̂ =

Corollary 2. For the same setting as Theorem 4, let γ̂ be
the output of M-SBL optimization in (5). If γ̂ ∈ ΘK , and
the conditions C1 and C2 hold for a given δ ∈ (0, 1),
then supp(γ̂) = S ∗ with probability exceeding 1 − 2δ.
Proof. Since γ̂ belongs to ΘK and also maximizes the MSBL objective log p(Y; γ), it follows that γ̂ is a solution to
the constrained cM-SBL(K) optimization in (16). Hence, the
statement of Corollary 2 follows directly from Theorem 4.
According to Corollary 2, M-SBL has vanishing support error
probability under conditions C1 and C2, however only in a
retrospective sense, i.e., when the M-SBL output γ̂ belongs to
the bounded hyperparameter set ΘK .
So far, we have analyzed the support error probability for the
∗
case where the true row support of X
 is nfixed
 to S n. In reality,
n
R(X) can assume any one of the 1 + 2 +· · ·+ K possible
supports in the collection SK . The support error probability
averaged over all possible supports can be evaluated as
X
err
=
Pavg
P(R(X) = S ∗ )P (ES ∗ ) ,
(36)
S ∗ ∈SK

where ES ∗ is as defined in (17). For example, if all supports
err
in SK are equiprobable, then Pavg
≤ 2δ under C1 and C2.
VI. D ISCUSSION
In this section, we interpret the sufficiency conditions C1
and C2 in Theorem 4 in the context of various interesting
cases.
A. Column Normalized Measurement Matrices
According to C2 in Theorem 4, the number of sufficient
MMVs grows quadratically with a∗ , the largest absolute value
of entries in A. If A has unit norm columns, then
 a∗ scales
inversely with the number of rows in A as O √1m . For
example, if A is a random matrix with i.i.d. N (0, √1m ) entries,

q
log n
with high probability. Likewise, if
then |a∗ | ≤ O
m

A is a partial DFT sensing matrix, then |a∗ | = √1m .
Often, column normalization is essential for accurate modelling of the norm-preserving nature of the application specific
measurement modality. In light of this, if A has normalized
columns, then one must account for O( √1m ) scaling of a∗
while interpreting the sufficient number of MMVs in C2.

arg max

log p (Y; γ).

(37)

γ∈{0,1}n ,||γ||0 =k

Binary valued hyperparameters can be accommodated as a
special case by setting γ min = γ max = 1 in our source signal
model. For γ min = γ max , according to Proposition 4, the ǫnet Θǫ (S) collapses to a single point for all S ∈ SK , which
ultimately amounts to κcov = 1. Thus, for binary γ and column
normalized A, by setting κcov = 1 and a∗ = √1m in C2, we
can conclude that the support error probability decays

 3/2 exponentially fast with increasing L provided L ≥ O K mlog n .
For continuous valued γ, i.e. γ(i) ∈ [γ min , γ max ] ∪ {0},
the number
MMVs
must scale as L ≥
 3/2of sufficient

5/2
O max K mlog n , K mlog K .
C. On the Restricted Isometry of A⊙A and Support Recovery
Guarantees for m ≤ K

After L surpasses the MMV threshold in C2, the support
error probability is guaranteed to decay exponentially with L,
provided η is strictly positive. The positivity of η (as defined
in (28)) is contingent upon the self Khatri-Rao product A ⊙ A
⊙
satisfying the 2K-RIP condition: δ2K
< 1.
Suppose A is an m × n sized random matrix comprising
i.i.d. subgaussian entries with zero mean and variance equal to
√1 . Then, according to [39, Theorem 3], the self Khatri-Rao
m
⊙
product A ⊙ A satisfies the RIP condition: δ2K
high

≤√ t with
K log3/2 n
.
probability for any t ∈ (0, 1) provided m ≥ O
t
Thus,
for
a
subgaussian
A,
it
is
sufficient
for
m
to
scale
as
√
O( K) in order to guarantee vanishing error probability for
M-SBL based reconstruction of all K or less sparse supports.
Furthermore, in [39, Theorem 3], it is also shown that the
self Khatri-Rao
A ⊙ A satisfies K-RIP for K as
 2 product

m
high as O log
.
This
implies that for a fixed n, O(m2 )
3n
sized supports can be perfectly recovered by M-SBL with high
probability using only finitely many MMVs. The exact number
of sufficient MMVs is given by C2. An order-wise sufficient
MMV condition is reported in Table I.
D. Support Recovery from Noiseless Measurements
For K < spark(A) − 1, it can be shown that as the noise
variance σ 2 → 0, the support error exponent D1/2 (pγ , pγ ∗ ) in
(24) grows in an unbounded fashion for all γ ∈ ΘK \Θ(S ∗ ).
This is formally proved below.
The 1/2-Rényi divergence between two multivariate Gaussian densities pγ i (y) = N (0, Σγ i ), i = 1, 2 is given by
D1/2 pγ 1 , pγ 2



=

log

Σγ 1 + Σγ 2
1
− log Σγ 1 Σγ 2
2
2

9

=
1
2


 1
1
log 0.5 H 2 + H− 2 .

(38)

1
2

where H , Σγ 1 Σ−1
γ 2 Σγ 1 is referred to as the discrimination
matrix. Since H is a normal matrix, it is unitarily diagonalizable. Let H = UΛUT , where Λ = diag(λ1 , . . . , λm ) with
λi ’s being the strictly positive eigenvalues of H, and U being
a unitary matrix with the eigenvectors of H as its columns.
The 1/2-Rényi divergence can be expressed in terms of λi as
m
 

X

−1/2
1/2
/2
log λi + λi
=
D1/2 pγ 1 , pγ 2
i=1
 

1
(λmax (H))1/2 + (λmax (H))−1/2 . (39)
≥ log
2
The above inequality is obtained by dropping all positive terms
in the summation except the one term which corresponds to
λmax (H), the maximum eigenvalue of H. Proposition 7 below
relates λmax (H) to the noise variance σ 2 .

Proposition 7. If K < spark(A) − 1, then, for any pair
γ 1 , γ 2 ∈ ΘK such that supp(γ 1 )\supp(γ 2 ) 6= φ, the maxi1
1
c1
2
mum eigenvalue of H , Σγ2 1 Σ−1
γ 2 Σγ 1 satisfies λmax (H) ≥ σ2
2
for some constant c1 > 0 independent of σ .
Proof. See Appendix G.
According to Proposition 7, in the limit σ 2  → 0,
λmax (H) → ∞, and consequently, D1/2 pγ 1 , pγ 2 grows
unbounded (due to (39)) whenever supp(γ 1 ) 6= supp(γ 2 ) and
K < spark(A) − 1. Based on this observation, we now state
Theorem 5 which lays forward the sufficient conditions for
exact support recovery in the noiseless case.

Theorem 5. Consider the noiseless MMV problem, with
observations Y = AX corresponding to an unknown
X satisfying assumption (A1). Suppose S ∗ is the true
nonzero row support of X with |S ∗ | ≤ K. Then, if a
solution γ̂ to the M-SBL optimization problem in (5)
satisfies γ̂ ∈ ΘK , then supp(γ̂) = S ∗ almost surely,
provided that K < spark(A) − 1. This result holds even
in the SMV case, i.e., when L = 1.
Proof. Under assumption A1, there exists a γ ∗ ∈ ΘK
such that every column in X is i.i.d. N (0, diag(γ ∗ )), and
supp(γ ∗ ) = S ∗ . Since γ̂ globally maximizes the M-SBL
objective L(Y; γ), it follows that L(Y; γ̂) ≥ L(Y; γ ∗ ) if
γ̂ 6= γ ∗ . Moreover, the following chain of implications holds.
{supp(γ̂) 6= S ∗ } =
⊆

{supp(γ̂) 6= supp(γ ∗ )} ⊆ {γ̂ 6= γ ∗ }
{L (Y; γ̂) ≥ L (Y; γ ∗ )} .

By applying Corollary 1, this further implies that
P (supp(γ̂) 6= S ∗ ) ≤ P (L (Y; γ̂) ≥ L (Y; γ ∗ ))


≤ exp −LD1/2 (pγ̂ , pγ ∗ ) /4 .

By using the lower bound (39) for D1/2 (pγ̂ , pγ ∗ ), we have
"
!#− L4
1 p
1
∗
P (supp(γ̂) 6= S ) ≤
λmax (H) + p
,
2
λmax (H)
(40)

1/2

1/2

∗
where H = Σγ̂ Σ−1
γ ∗ Σγ̂ . Since both γ̂ and γ belong to
ΘK , as long as K < spark(A) − 1, by Proposition 7, σ 2 → 0
results in λmax (H) → ∞ which in turn drives the RHS in (40)
to zero for L ≥ 1.

From Theorem 5, we conclude that, in the noiseless scenario
and for X satisfying assumption A1, M-SBL requires only a
single measurement vector (L = 1) to perfectly recover any
K < spark(A)−1 sized support. This result is in line with the
sufficient conditions identified for successful support recovery
by M-SBL in [25, Theorem 1]. However, unlike in [25], the
nonzero rows of X need not be orthogonal. Also, our result
improves over the k ≤ m/2 condition shown in [31].
E. Impact of Measurement Noise on Sufficient MMVs
∗
For σ 2 > 0, the error exponent term D1/2
in (25) is always
bounded. This implies that unlike in the noiseless case, a single
MMV is no longer sufficient, and multiple MMVs are needed
to drive the error probability close to zero.
A close inspection of the abstract MMV bound in (29)
reveals that the noise variance influences the error probability
in a twofold manner: (i) through η, and (ii) through the ǫ-net
cardinality bound κcov . As σ 2 increases, η decreases polynomially (see (28)) while log κcov increases at most logarithmically
(see (35)). The overall effect is captured by condition C2
in Theorem 4, which suggests that if the noise variance is
very high relative to γmax , it is sufficient to have a roughly
quadratically larger number of MMVs to guarantee the desired
probability of error. As the noise variance approaches zero, the
MMV bound in C2 loosens and is not informative.

F. Support Recovery using RIP Compliant A
If the measurement matrix A is itself 2K-RIP compliant,
then the support error probability vanishes for fewer MMVs
than specified by C2. In [31], the authors showed that for 2KRIP
 compliant
 A, and binary hyperparameters γ in M-SBL,
O logloglogn n MMVs suffice for error probability to vanish.
Here, we consider the general case of continuous valued γ.
⊙
A 2
If A has unit norm columns, then δ2K
≤ (δ2K
) by [39,
Theorem 1], which results in tighter bounds for η and κcov :


2 
A 2
1
γ min
1 − (δ2K
)
, (41)
η ≥
A
8Ka2∗ σ 2 + γ max
1 + δ2K
( 
K )
K 5/2 a2∗ (γ max − γ min )ζ ′
κcov ≤ max 1,
, (42)
A )
γ min (1 − δ2K

 2
A
1+δK
(σ +γ max )3
where ζ ′ = γ144
2
A
σ2 +γ min . By using (41) and (42) in
1−δ
min
K
Theorem 3, we obtain the following condition on the number
of MMVs to guarantee P(ES ∗ ) ≤ 2δ:


2

 
3enK(1+δ)
64Ka2∗ σ 2 +γ max
, log κcov ,
max log
L≥
A )
γ min
δ
(1−δ2K
(43)
where κcov is bounded as in (42).
From (43), we conclude that for A with unit norm columns,
P(ES ∗ ) can be driven arbitrarily close to zero if L exceeds

10

TABLE I
S UFFICIENT CONDITIONS FOR EXACT SUPPORT RECOVERY OF K OR LESS
SPARSE SUPPORT IN M-SBL WHEN ITS OUTPUT γ̂ ∈ ΘK .
General A with
unit norm columns
√
|a∗ | = O(1/ m)
A < 1
δ2K

Sufficient number of MMVs
Binary γ

σ2

=0 L≥1

σ2 > 0 L ≥ O

⊙
δ2K
<1

σ2 ≥ 0 L ≥ O

A with i.i.d.
N (0, √1m ) entries
m = Ω(K log n)

Noise

Noise



Continuous γ

K log n
m



3





K 2 log n
m

L≥1

L≥O
L≥O





K 2 log K
m
5

K 2 log K
m



Binary γ

Continuous γ
L≥1

σ2 > 0 L ≥ O (log n)
L ≥ O (K log K)

 2
√

√
3
log K
m = O( K log 2 n) σ2 ≥ 0 L ≥ O K log n L ≥ O K 3/2
log

n




n K 2 log K
,
. For A with i.i.d. N (0, √1m )
O max K log
m
m
entries, the sufficient MMV condition is reported in Table I.
VII. C ONCLUSIONS
In this work, we presented new sufficient conditions for
exact support recovery by M-SBL in the joint sparse support
recovery problem. The sufficient conditions are expressed in
terms of the number of MMVs and the restricted isometry
constants of the measurement matrix and its self Khatri-Rao
product. The new conditions cater to a wider, more useful
class of Gaussian signals, and dispenses with the restrictive
row-orthogonality condition on the signal matrix required in
the previous results [22], [25].
We have shown that M-SBL is capable of recovering ksparse support from fewer than k measurements per MMV
for both noiseless and noisy measurements. In the noiseless
case, a single MMV suffices for perfect recovery of any ksparse support, provided that spark(A) > k + 1. In case of
noisy measurements, M-SBL can perfectly recover any k-sized
support with vanishing error probability using finitely many
⊙
MMVs, provided that A⊙A satisfies the RIP condition: δ2k
<
⊙
1. For random subgaussian A, the RIP condition
δ
<
1
 is
 2k
m2
[39],
satisfied with high probability for k as high as O log
3n
which suggests that M-SBL can successfully recover O(m2 )
sized supports for fixed n.
There still remain the following open questions about MSBL’s support recovery performance:
(i) What are necessary conditions for exact support recovery
in terms of the number of required MMVs?
(ii) Under what conditions do all local maxima of the M-SBL
objective yield the correct support estimate?
(iii) How is the support recovery performance impacted by
inter and intra vector correlations in the signals?
A PPENDIX
A. Proof of Proposition 2
Let Σγ 1 , Σγ 2 ∈ MK . From Taylor’s expansion of
ψ(Σγ 2 ) = − log det Σγ 2 around Σγ 1 , we have
ψ(Σγ 2 )

= ψ(Σγ 1 ) + h∇ψ(Σγ 1 ), Σγ 2 − Σγ 1 i

∇2 ψ(Σγ t )[Σγ 2 − Σγ 1 ]

(44)

with Σγ t , tΣγ 1 + (1 − t)Σγ 2 for some t ∈ [0, 1].
2
The
 Hessian operator ∇ ψ(Σγ t )[Σγ 2− Σγ 1 ] is evaluated as
−1
tr Σ−1
γ t (Σγ 2−Σγ 1 )Σγ t (Σγ 2−Σγ 1 ) . To prove the strong
convexity of ψ over Mk , one needs to show the existence
of a positive constant mψ such that



Sufficient number of MMVs

σ2 = 0 L ≥ 1

+

∇2 ψ(Σγ t )[Σγ 2 − Σγ 1 ] ≥ mψ Σγ 2 − Σγ 1

2
F

(45)

for all Σγ 1 , Σγ 2 ∈ MK . Note that Σγ t = σ 2 Im + AΓt AT ,
where Γt = diag(γ t ) and γ t = tγ 1 + (1 − t)γ 2 is at most
2K-sparse. Further, let S denote the nonzero support of ∆γ =
γ 2 − γ 1 , and let ∆Γ = diag(∆γ). Then, we have Σγ 2 −
Σγ 1 = AS ∆ΓS ATS . Using these new definitions, the Hessian
term in (44) can be lower bounded as


−1
)
−
Σ
)Σ
(Σ
−
Σ
tr Σ−1
(Σ
γ
γ
γ
γ
γt
γt
1
2
1
2


−1
T −1
= tr Σγ t AS ∆ΓS AS Σγ t AS ∆ΓS ATS


(a)
T
T −1
= vec (∆ΓS ) ATS Σ−1
γ t AS ⊗ AS Σγ t AS vec (∆ΓS )

T

(b)
−1
−1
−1
−1
= ∆γ TS Σγ t2 AS ⊙ Σγ t2 AS
Σγ t2 AS ⊙ Σγ t2 AS ∆γ S


(c)
T −1
= ∆γ TS ATS Σ−1
A
◦
A
Σ
A
S
S ∆γ S
γt
S γt
−1
(d)
T
(AS ◦ AS ) ∆γ S
= ∆γ TS (AS ◦ AS ) Σγ t ◦ Σγ t


T
T
≥ ∆γ S (AS ◦AS ) (AS ◦AS ) ∆γ S / Σγ t◦Σγ t 2
2

2
AS ∆ΓS ATS F
||(AS ⊙ AS ) ∆γ S ||2
=
Σγ t ◦ Σγ t 2
Σγ t ◦ Σγ t 2
−1

≥
sup
|||Σγ ◦Σγ |||2
Σγ 2−Σγ 1
(e)

≥

γ∈Rn
+ ,||γ||0 ≤2K

2
F

. (46)

In above, steps (a) and (d) follow from facts 7.4.9 and 8.21.49
in [40], respectively. Step (b) is due to Columnwise KhatroRao product (⊙) being a submatrix of the Kronecker product
(⊗). Step (c) is a consequence of the property: (A⊙A)T (A⊙
A) = AT A ◦ AT A satisfied by any matrix A. Step (e) is due
to Hadamard product (◦) being a submatrix of the columnwise
Khatri-Rao product matrix (⊙).
−1

From (46), mψ = supγ∈Rn+ ,||γ||0 ≤2K |||Σγ ◦Σγ |||2
>0
is a valid choice for the strong convexity constant of ψ over
domain MK , thereby establishing the Proposition.
B. Proof of Theorem 1
From Proposition 2, − log det(·) is a strongly convex function in Θ(Σ1 , Σ2 ) with strong convexity constant m∗ . Any
strongly convex function f : X → R with strong convexity
constant mf satisfies tf (x)+(1−t)f (y) ≥ f (x + (1 − t)y)+
0.5mf t(1 − t)kx − yk22 for t ∈ [0, 1] and x, y ∈ X . Thus, for
Σ1 , Σ2 ∈ Θ(Σ1 , Σ2 ),
− log |(1 − α)Σ1 + αΣ1 | ≤ −(1 − α) log |Σ1 | − α log |Σ2 |
1
2
− m∗ α(1 − α) ||Σ2 − Σ1 ||F .
2

11

By rearranging terms, we obtain
1
1
|(1 − α)Σ1 + αΣ1 |
≥ m∗ α ||Σ2 − Σ1 ||2F .
log
1−α
α
2(1 − α)
4
|Σ1 |
|Σ2 |
The left hand side in the above equation is precisely
Dα (p1 , p2 ) Thus, we have the desired result.
C. Proof of Proposition 4
The following stepwise procedure shows how to construct a
δ-net of Θ(S) (with respect to the Euclidean distance metric)
which is entirely contained in Θ(S).
1) Consider an δ-blow up of Θ(S), denoted by Θ↑δ (S) ,
{x : ∃x′ ∈ Θ(S) such that kx − x′ k2 ≤ δ}.
2) Let Θδ↑δ (S) be a δ-net of Θ↑δ (S). Some points in Θδ↑δ (S)
may lie outside Θ(S).
3) Let P denote the set containing the projections of all
points in Θδ↑δ (S) ∩ Θ(S)c onto the set Θ(S). By conδ
c
struction, P ⊂ Θ(S),
and |P| ≤ |Θ

 ↑δ (S) ∩ Θ(S) |.

4) Then, Θδ (S) , Θδ↑δ (S) ∩ Θ(S) ∪ P is a valid δ-net
of Θ(S) which is entirely contained in Θ(S).
To prove the validity of the above δ-net construction, we
need to show that for any γ ∈ Θ(S), there exists an element a
in Θδ (S) such that kγ−ak2 ≤ δ. Let γ be an arbitrary element
in Θ(S). Then, γ also belongs to the larger set Θ↑δ (S), and
consequently, there exists γ ′ ∈ Θδ↑δ (S) such that kγ − γ ′ k2 ≤
′
δ. Now, there are two
and (ii) γ ′ ∈
/ Θ(S).

 cases. (i) γ ∈ Θ(S),
′
δ
In case (i), γ ∈ Θ↑δ (S) ∩ Θ(S) , and hence also belongs
to Θδ (S). Further, kγ − γ ′ k2 ≤ δ. Hence a = γ ′ will work.
In case (ii), γ ′ ∈ Θδ↑δ (S) ∩ Θ(S)c . Let γ ′′ be the projection
of γ ′ onto Θ(S), then γ ′′ must belong to P, and hence must
also belong to Θδ (S). Note that since γ ′′ is the projection
of γ ′ onto the convex set Θ(S), for any γ ∈ Θ(S), we have
hγ − γ ′′ , γ ′ − γ ′′ i ≤ 0. Further, we have
δ

≥

=

≥

2

2

||γ − γ ′ ||2 = ||(γ − γ ′′ ) + (γ ′′ − γ ′ )||2
2

2

||γ − γ ′′ ||2 + ||γ ′′ − γ ′ ||2 + 2hγ − γ ′′ , γ ′′ − γ ′ i
2

||γ − γ ′′ ||2 .

(47)

The last inequality is obtained by dropping the last two
nonnegative terms in the RHS. From (47), a = γ ′′ will work.
Since case (i) and (ii) together are exhaustive, Θδ (S) in
step-4 is a valid δ-net of Θ(S) which is entirely inside Θ(S).
δ
pCardinality of Θ (S): The diameter of Θ(S) is
|S|(γ max − γ min ). Based on the construction in step4, the cardinality of Θδ (S) can be upper bounded as:
|Θδ (S)|

= |Θδ↑δ (S) ∩ Θ(S)| + |Θδ↑δ (S) ∩ Θ(S)c |

= |Θδ↑δ (S)| ≤ |δ-net of ℓ2 ball of radius
p
|S|(γ max − γ min ) in R|S|
 
|S| 
p
. (48)
≤ max 1, 3 |S|(γ max − γ min )/δ

The last step is an extension of the volumetric arguments in
[36] to show that the δ-covering number of a unit ball B1 (0) in
Rk with respect to the standard Euclidean norm ||·||2 satisfies

δ
Ncov
(B1 (0), ||·||2 ) ≤ (3/δ)k . The max operation with unity
covers the case when δ is larger than the diameter of Θ(S).
Now consider the modified net Θǫ/CL,S (S) obtained by
ǫ
setting δ = CL,S
in steps 1-4, where CL,S is the Lipschitz
constant of L(Y, γ) with respect to γ ∈ Θ(S). We claim
ǫ/C
that ΘS L,S is the desired set which simultaneously satisfies
conditions (i) and (ii) stated in Proposition 4.
ǫ/C
To show condition (i), we observe that since ΘS L,S is
an (ǫ/CL,S )-net of Θ(S) with respect to ||·||2 , for any γ ∈
Θ(S), there exists a γ ′ ∈ Θǫ/CL,S (S) such that ||γ − γ ′ ||2 ≤
ǫ/CL,S . Since L(Y, γ) is CL,S -Lipschitz in Θ(S), it follows
that |L(Y, γ) − L(Y, γ ′ )| ≤ CL,S ||γ − γ ′ ||2 ≤ ǫ.
Condition (ii) follows from (48) by setting δ = ǫ/CL,S .

D. Proof of Theorem 2
For continuous probability densities pγ and pγ ∗ defined on
the observation space Rm , the tail probability of the random
variable log (pγ (Y)/pγ ∗ (Y)) has a Chernoff upper bound
with parameter t > 0 as shown below.




L
X
pγ (Y)
p
(y
)
γ
j
P log
≥ −ǫ
= P
log
≥ −ǫ
∗ (yj )
pγ ∗ (Y)
p
γ
j=1



L
X
pγ (yj ) 
exp (tǫ)
log
≤ Epγ ∗ exp t
p
γ ∗ (yj )
j=1
L



pγ (y)
exp (tǫ)
= Epγ ∗ exp t log
pγ ∗ (y)
"
t #!L
pγ (y)
exp (tǫ)
= Epγ ∗
pγ ∗ (y)
Z
L
=
pγ (y)t pγ ∗ (y)1−t dy
exp (tǫ)
y∈Y
i

h  ǫ
− (t − 1)Dt (pγ , pγ ∗ ) .(49)
= exp −L t −
L
In the above, the first and third steps follow from the independence of yj . The second step is the application of Chernoff
bound. The last step is obtained by using the definition of the
Rényi divergence and rearranging the terms in the exponent.
By introducing the function ψ(t) = (t − 1)Dt (pγ , pγ ∗ ), the
Chernoff bound (49) can be restated as


i

h  ǫ
pγ (Y)
− ψ(t) .
≥ −ǫ ≤ exp −L t −
P log
pγ ∗ (Y)
L
(50)


For t = arg supt>0 t − Lǫ − ψ(t) , the upper bound in (50)
attains its tightest value exp (−Lψ ∗ (−ǫ/L)), where ψ ∗ is the
Legendre transform of ψ.
E. Proof of Proposition 5

∗
2
Σγ 2 − Σγ 1 F .
From Theorem 1, Dα pγ 1 , pγ 2 ≥ αm
4
In the following arguments, we derive separate lower bounds
2
for m∗ and Σγ 2 − Σγ 1 F .
2
The term Σγ 2 − Σγ 1 F can be lower bounded as:
Σγ 2 − Σγ 1

2
F

= A (Γ2 − Γ1 ) AT

2
F

12

⊙
= ||(A ⊙ A) (γ 2 − γ 1 )||22 ≥ (1 − δ2K
)kd γ 2min .(51)


−1
≤ L |||RY |||2 tr ATS Σ−1
γ Σγ AS


−1/2
−1/2
−1
= L |||RY |||2 tr ΓS ÃTS Σ−1
γ Σγ ÃS ΓS


−1
≤ (L |||RY |||2 /γ min ) tr ÃTS Σ−1
γ Σγ ÃS

 

T −1
tr
Ã
Σ
Ã
≤ L |||RY |||2 Σ−1
/γ
S
S
γ
γ
min
2

(a) 
≤ L |||RY |||2 Σ−1
min
(m,
|S|)
/γ min
γ
2

The last inequality is a consequence of A⊙A satisfying RIP of
2
⊙
order 2K with RIC δ2K
and the fact that ||γ 2 −γ 1 ||2 ≥ kd γ 2min .
Next, we derive a lower bound for m∗ . We note that
(σ 2 Im + AΓAT ) ◦ (σ 2 Im + AΓAT )

|||Σγ ◦ Σγ |||2 =

≤ σ 4 + 2σ 2 γ max
+ γ 2max

AS ATS ◦ Im

AS ATS ◦ AS ATS

2

2
2

,

(52)

where S = supp(γ) and |S| ≤ 2K. The last step follows from
the triangle inequality. Let a∗ denote the maximum absolute
value of elements in A. Then, it can be shown that
AS ATS ◦ Im

AS ATS

◦

AS ATS

2
2

≤

≤

a2∗ K,

a2∗ K

(53)
ATS AS

2

. (54)

Using (53) and (54) in (52), and by Proposition 2, we have
−1

2
ATS AS 2
.
(55)
m∗ ≥ Ka2∗ σ 2 + γ max
Finally, we get the desired result by using Theorem 1, and
2
combining the lower bounds for m∗ and Σγ 2 − Σγ 1 F in
(55) and (51), respectively.
F. Proof of Proposition 6
The log-likelihood L(Y; γ) can be expressed as the sum
f (γ) + g(γ) with f (γ) = −L log |Σγ | and g(γ) =
2
T
−Ltr Σ−1
γ RY . Here, Σγ = σ Im + AΓA .
First, we derive an upper bound for the Lipschitz constant of
f (γ) = −L log |Σγ | for γ ∈ Θ(S). By mean value theorem,
any upper bound for k∇γ f (γ)k2 also serves as an upper bound
for the Lipschitz constant of f . So, we derive an upper bound
(γ)
for k∇γ f (γ)k2 . Note that ∂f
= L(ai Σ−1
γ ai ) for i ∈ S,
γ(i)
th
and 0 otherwise. Here, ai denotes the i column of A. Then,
k∇γ f (γ)k2 can be upper bounded as shown below.
X
k∇γ f (γ)k2 ≤ k∇γ f (γ)k1 = L
aTi Σ−1
γ ai
u∈S


= L tr ATS (σ 2 Im + AS ΓS ATS )−1 AS

 
−1/2
−1/2
= L tr ΓS ÃTS (σ 2 Im + ÃS ÃTS )−1 ÃS ΓS


(a)
tr ÃTS (σ 2 Im + ÃS ÃTS )−1 ÃS
= L Γ−1
S
2
(b)

≤ (L min (m, |S|)) /γ min .
1/2
AS ΓS .

(56)

In the above, step (a) follows from
where ÃS =
the trace inequality tr(AB) ≤ k|Ak|2 tr(B) for any positive
definite matrices A and B. Step (b) follows from the observation that input argument of trace operator has min (m, |S|)
nonzero eigenvalues, all of them less than unity.
We now shift focus to the second term g(γ) of the loglikeT −1
−1
lihood. Note that ∂g(γ)
∂γ(i) = L(ai Σγ RY Σγ ai ) for i ∈ S,
and 0 otherwise. Then, k∇γ g(γ)k2 can be upper bounded as
X
−1
k∇g(γ)k2 ≤ k∇γ g(γ)k1 = L
aTi Σ−1
γ RY Σγ ai
= L tr

i∈S

−1
ATS Σ−1
R
Σ
Y
γ
γ AS

(b)

≤ (L |||RY |||2 min (m, |S|)) /γ min σ 2 .

(57)

1/2

where ÃS , AS ΓS . The inequality in (57a) follows from

ÃTS Σγ−1 ÃS having min (m, |S|) nonzero eigenvalues, all
of them less than unity. The last inequality in (57b) is due to
Σ−1
≤ 1/σ 2 . Finally, the Lipschitz constant CL,S can
γ
2
be bounded as CL,S ≤ k∇γ f (γ)k2 + k∇γ g(γ)k2 . Thus, by
combining (56) and (57), and noting that min(m, |S|) ≤ K,
we obtain the desired result.
G. Proof of Proposition 7
1

1

2
Let µ∗ be the largest eigenvalue of Σγ2 1 Σ−1
γ 2 Σγ 1 . Then,


i

tr Σγ−1
Σ
γ
1
1 h 2  −1 
2
T
µ∗ ≥
σ Tr Σγ 2 + Tr Σ−1
=
γ 2 AΓ1 A
m
m

1  −1
(58)
tr Σγ 2 AΓ1 AT .
≥
m

Here, the second step is setting Σγ 1 = σ 2 Im + AΓ1 AT . The
last inequality

 is obtained by dropping the strictly positive
σ2
−1
tr
Σ
term.
γ2
m
Let S1 and S2 be the nonzero supports of γ 1 and γ 2 , respectively. Further, let the eigendecomposition of Σγ 2 be UΛUT ,
where Λ = diag(λ1 , . . . , λm ), λi ’s are the eigenvalues of Σγ 2
and U is a unitary matrix with columns as the eigenvectors
of Σγ 2 . Then, U can be partitioned as [U2 U2⊥ ], where the
columns of U2 and U2⊥ span the orthogonal complementary
subspaces Col(AS2 ) and Col(AS2 )⊥ , respectively. Further, let
Λ2 and Λ2⊥ be |S2 | × |S2 | and ((m − |S2 |) × (m − |S2 |))
sized diagonal matrices containing the eigenvalues in Λ corresponding to the eigenvectors in U2 and U2⊥ , respectively.
We observe that Λ2⊥ = σ 2 Im−|S2 | .
−1 T
−1 T
By setting Σ−1
γ 2 = U2 Λ2 U2 + U2⊥ Λ2⊥ U2⊥ in (58), we
get

1
T
T
tr U2 Λ−1
+ tr U2⊥ Λ−1
UT
µ∗ ≥
2 U2 AΓ1 A
2⊥ 2⊥
m


1
T
T
AΓ1 AT ≥ tr Λ−1
⊥ U2⊥ AΓ1 A U2⊥ ,
2
m
where the last inequality is due to nonnegativity of the dropped
first term. Since UT2⊥ AS2 = 0 by construction of U2⊥ ,

1 
T
T
c
c
µ∗ ≥ tr Λ−1
⊥ U2⊥ AS2 Γ1,S2 AS c U2⊥
2
2
m
m−|S2 |
1 X
(u ⊥ )T AS2c Γ1,S2c ATS2c u2⊥ ,i
=
mσ 2 i=1 2 ,i
=

m−|S2 |
1 X
(u ⊥ )T AS1 \S2 Γ1,S1 \S2 ATS1 \S2 u2⊥ ,i . (59)
mσ 2 i=1 2 ,i

13

In the above, u2⊥ ,i denotes the ith column of U2⊥ . The last
equality is obtained by observing that the nonzero elements of
γ 1,S2c are located in the index set S1 \S2 .
We now prove that if K < spark(A) − 1, then there exists
at least one strictly positive term in the above summation. Let
us assume the contrary, i.e., let each term in the summation
in (59) be equal to zero. This implies that the columns of
1/2
U2⊥ belong to Null(Γ1,S1 \S2 ATS1 \S2 ), which means that they
also belong to Null(AS1 \S2 Γ1,S1 \S2 ATS1 \S2 ). Since, for a
symmetric matrix, the row and column spaces are equal and
orthogonal to the null space of the matrix, it follows that
1/2
Col(AS1 \S2 Γ1,S1 \S2 ATS1 \S2 ) (same as Col(AS1 \S2 Γ1,S1 \S2 ))
is spanned by the columns of U2 , or equivalently by columns
of AS2 . Thus, every column in AS1 \S2 can be expressed
as a linear combination of columns in AS2 . This contradicts
our initial assumption that K + 1 < spark(A), implying that
any K + 1 or fewer columns of A are linearly independent.
Therefore, we conclude that there is at least one strictly
positive term in the summation in (59), and consequently there
exists a constant c1 > 0 such that µ∗ ≥ c1 /σ 2 .
ACKNOWLEDGEMENT
The authors gratefully acknowledge Heng Qiao from UC
San Diego for helpful comments and discussions.
R EFERENCES
[1] M. Masood, L. Afify, and T. Al-Naffouri, “Efficient coordinated recovery
of sparse channels in massive MIMO,” IEEE Trans. Signal Process.,
vol. 63, no. 1, pp. 104–118, 2015.
[2] Y. Barbotin, A. Hormati, S. Rangan, and M. Vetterli, “Estimation of
sparse MIMO channels with common support,” IEEE Trans. Commun.,
vol. 60, no. 12, pp. 3705–3716, 2012.
[3] R. Prasad, C. Murthy, and B. Rao, “Joint channel estimation and
data detection in MIMO-OFDM systems: A sparse Bayesian learning
approach,” IEEE Trans. Signal Process., vol. 63, no. 20, pp. 5369–5382,
Oct. 2015.
[4] D. Baron, M. F. Duarte, M. B. Wakin, S. Sarvotham, and R. G.
Baraniuk, “Distributed compressive sensing,” CoRR, vol. abs/0901.3403,
2009. [Online]. Available: http://arxiv.org/abs/0901.3403
[5] H. Zhu, G. Giannakis, and A. Cano, “Distributed in-network channel
decoding,” IEEE Trans. Signal Process., vol. 57, no. 10, pp. 3970–3983,
Oct. 2009.
[6] S. Ji, D. Dunson, and L. Carin, “Multitask compressive sensing,” IEEE
Trans. Signal Process., vol. 57, no. 1, pp. 92–106, Jan. 2009.
[7] T. Erseghe, “A distributed and maximum-likelihood sensor network
localization algorithm based upon a nonconvex problem formulation,”
IEEE Trans. Signal Inf. Process. Netw., vol. 1, no. 4, pp. 247–258, Dec.
2015.
[8] Z. Tan, P. Yang, and A. Nehorai, “Joint sparse recovery method for
compressed sensing with structured dictionary mismatches,” IEEE Trans.
Signal Process., vol. 62, no. 19, pp. 4997–5008, Oct. 2014.
[9] Z. Fanzi, C. Li, and Z. Tian, “Distributed compressive spectrum sensing
in cooperative multihop cognitive networks,” IEEE J. Sel. Topics Signal
Process., vol. 5, no. 1, pp. 37–48, Feb. 2011.
[10] Y. Wang, A. Pandharipande, Y. L. Polo, and G. Leus, “Distributed
compressive wide-band spectrum sensing,” in Information Theory and
Applications Workshop, Feb. 2009, pp. 178–183.
[11] J. Bazerque and G. Giannakis, “Distributed spectrum sensing for cognitive radio networks by exploiting sparsity,” IEEE Trans. Signal Process.,
vol. 58, no. 3, pp. 1847–1862, Mar. 2010.
[12] S. Cotter, B. Rao, K. Engan, and K. Kreutz-Delgado, “Sparse solutions
to linear inverse problems with multiple measurement vectors,” IEEE
Trans. Signal Process., vol. 53, no. 7, pp. 2477–2488, Jul. 2005.
[13] J. Chen and X. Huo, “Theoretical results on sparse representations of
multiple-measurement vectors,” IEEE Trans. Signal Process., vol. 54,
no. 12, pp. 4634–4643, Dec. 2006.

[14] J. Tropp, A. Gilbert, and M. Strauss, “Simultaneous sparse approximation via greedy pursuit,” in Proc. ICASSP, vol. 5, 2005, pp. 721–724.
[15] S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive
Sensing. Birkhäuser Basel, 2013.
[16] M. Duarte, S. Sarvotham, D. Baron, M. Wakin, and R. Baraniuk,
“Distributed compressed sensing of jointly sparse signals,” in Proc.
Asilomar Conf. on Signals, Syst., and Comput., Oct. 2005, pp. 1537–
1541.
[17] B. D. Rao, K. Engan, and S. Cotter, “Diversity measure minimization
based method for computing sparse solutions to linear inverse problems
with multiple measurement vectors,” in Proc. ICASSP, vol. 2, May 2004,
pp. 369–372.
[18] J. D. Blanchard, M. Cermak, D. Hanle, and Y. Jing, “Greedy algorithms
for joint sparse recovery,” IEEE Trans. Signal Process., vol. 62, no. 7,
pp. 1694–1704, Apr. 2014.
[19] K. Lee, Y. Bresler, and M. Junge, “Subspace methods for joint sparse
recovery,” IEEE Trans. Inf. Theory, vol. 58, no. 6, pp. 3613–3641, Jun.
2012.
[20] J. M. Kim, O. K. Lee, and J. C. Ye, “Compressive music: Revisiting the
link between compressive sensing and array signal processing,” IEEE
Trans. Inf. Theory, vol. 58, no. 1, pp. 278–301, Jan. 2012.
[21] M. Mishali and Y. C. Eldar, “Reduce and boost: Recovering arbitrary
sets of jointly sparse vectors,” IEEE Trans. Signal Process., vol. 56,
no. 10, pp. 4692–4702, Oct. 2008.
[22] O. Balkan, K. Kreutz-Delgado, and S. Makeig, “Localization of more
sources than sensors via jointly-sparse Bayesian learning,” IEEE Signal
Process. Lett., vol. 21, no. 2, pp. 131–134, Feb. 2014.
[23] P. Pal and P. P. Vaidyanathan, “Pushing the limits of sparse support
recovery using correlation information,” IEEE Trans. Signal Process.,
vol. 63, no. 3, pp. 711–726, Feb. 2015.
[24] S. Khanna and C. R. Murthy, “Rényi divergence based covariance
matching pursuit of joint sparse support,” in Sig. Process. Advances
in Wireless Commun. (SPAWC), Jul. 2017, pp. 1 – 5.
[25] D. Wipf and B. Rao, “An empirical Bayesian strategy for solving
the simultaneous sparse approximation problem,” IEEE Trans. Signal
Process., vol. 55, no. 7, pp. 3704–3716, Jul. 2007.
[26] S. E. Chiu and B. Rao, “Correlation learning on joint support recovery
for more sources than measurements,” in IEEE Sensor Array and
Multichannel Signal Processing Workshop (SAM), Jul. 2016, pp. 1–5.
[27] P. Pal and P. P. Vaidyanathan, “Parameter identifiability in sparse
Bayesian learning,” in Proc. ICASSP, May 2014, pp. 1851–1855.
[28] E. van den Berg and M. P. Friedlander, “Theoretical and empirical results
for recovery from multiple measurements,” IEEE Trans. Inf. Theory,
vol. 56, no. 5, pp. 2516–2527, May 2010.
[29] Y. C. Eldar and M. Mishali, “Robust recovery of signals from a
structured union of subspaces,” IEEE Trans. Inf. Theory, vol. 55, no. 11,
pp. 5302–5316, Nov. 2009.
[30] M. E. Davies and Y. C. Eldar, “Rank awareness in joint sparse recovery,”
IEEE Trans. Inf. Theory, vol. 58, no. 2, pp. 1135–1146, Feb. 2012.
[31] G. Tang and A. Nehorai, “Performance analysis for sparse support
recovery,” IEEE Trans. Inf. Theory, vol. 56, no. 3, pp. 1383–1399, Mar.
2010.
[32] Y. Jin and B. D. Rao, “Support recovery of sparse signals in the presence
of multiple measurement vectors,” IEEE Trans. Inf. Theory, vol. 59,
no. 5, pp. 3139–3157, May 2013.
[33] S. Park, N. Y. Yu, and H. N. Lee, “An information-theoretic study for
joint sparsity pattern recovery with different sensing matrices,” IEEE
Trans. Inf. Theory, vol. PP, no. 99, pp. 1–1, 2017.
[34] A. Banerjee, S. Merugu, I. S. Dhillon, and J. Ghosh, “Clustering with
Bregman divergences,” J. Mach. Learn. Res., vol. 6, pp. 1705–1749,
Dec. 2005.
[35] I. S. Dhillon and J. A. Tropp, “Matrix nearness problems with Bregman divergences,” SIAM Journal on Matrix Analysis and Applications,
vol. 29, no. 4, pp. 1120–1146, 2008.
[36] R. Vershynin, “On the role of sparsity in compressed sensing and random
matrix theory,” in Proc. CAMSAP, Dec. 2009, pp. 189–192.
[37] M. Gil, F. Alajaji, and T. Linder, “Rényi divergence measures for commonly used univariate continuous distributions,” Information Sciences,
vol. 249, pp. 124 – 131, Nov. 2013.
[38] R. Vershynin, “Introduction to the non-asymptotic analysis of random
matrices,” CoRR, vol. abs/1011.3027, 2011. [Online]. Available:
http://arxiv.org/abs/1011.3027
[39] S. Khanna and C. R. Murthy, “On the restricted isometry of the
columnwise Khatri-Rao product,” CoRR, vol. abs/1709.05789, 2017.
[Online]. Available: http://arxiv.org/abs/1709.05789
[40] D. S. Bernstein, Matrix Mathematics: Theory, Facts, and Formulas
(Second Edition). Princeton University Press, 2009.

