Distribution-aware Block-sparse Recovery via Convex Optimization

We study the problem of reconstructing a block-sparse signal from compressively sampled measurements. In certain applications, in addition to the inherent block-sparse structure of the signal, some prior information about the block support, i.e. blocks containing non-zero elements, might be available. Although many block-sparse recovery algorithms have been investigated in Bayesian framework, it is still unclear how to incorporate the information about the probability of occurrence into regularization-based block-sparse recovery in an optimal sense. In this work, we bridge between these fields by the aid of a new concept in conic integral geometry. Specifically, we solve a weighted optimization problem when the prior distribution about the block support is available. Moreover, we obtain the unique weights that minimize the expected required number of measurements. Our simulations on both synthetic and real data confirm that these weights considerably decrease the required sample complexity.


I. INTRODUCTION
C OMPRESSED Sensing (CS) has emerged in the past decade as a modern technique for recovering a sparse vector x P R n from compressed measurements (see [1], [2] for more explanations about this field). In this work, we consider signals x P R n that have a block-sparse structure, namely their non-zero entries appear in blocks. This property has been referred in the literature as block sparsity. It is common to use the following optimization problem for recovering the signal x from compressive measurements. P η 1,2 : min zPR n }z} 1,2 :" where A P R mˆn represents a fat measurement matrix with m ! n, tV b u q b"1 are the default disjoint blocks of size tk b u q b"1 that partition the set t1, ..., qu, y :" Ax`e P R m is the observation vector 1 , e is the noise term which is considered to be i.i.d. Gaussian with variance σ 2 , and η is an upper-bound for }e} 2 . Most of the earlier literature in block-sparse recovery is focused on the case of single constraint }Az´y} 2 ď η. However, in many applications such as DNA micro-arrays [3], [4], computational neuroscience [5] multi-band signal reconstruction, multiple measurement vector (MMV) problem S. Daei and F. Haddadi are with the School of Electrical Engineering, Iran University of Science & Technology. A. Amini is with EE department, Sharif University of Technology. 1 Our analysis holds for both real and complex-valued signals and measurements. [6], and the reconstruction of signals in union of subspaces [7], [8] [9], [10], there exist additional information (or alternatively additional constraints in P η 1,2 ) about the signal of interest. In this work, we explore the benefits of having access to extra information about the distribution of the block support (blocks containing non-zero elements) on the required number of measurements. To this end, we propose the optimization problem P η 1,2,w : min zPR n }z} 1,2,w :" where the quantities w b , b " 1, ..., q are some positive scalars, w :" rw 1 , ..., w q s T and M is some predefined model that restricts the feasible set of the solution. Specifically, we consider two new models for prior information that are of practical interest: ‚ Model 1 (Prior distribution) : We assume that the prior distribution of the block support is available. Under this setting, there are known probabilities associated with each block index b P t1, ..., qu. Namely, where bsuppp¨q returns the block support of a vector. ‚ Model 2 (Multiple Bayesian estimates): We consider L independent Bayesian estimates of bsupppxq i.e. tP i u L i"1 with accuracies We assume, without loss of generality, α i 's sum to unity. To each subset P i Ă t1, ..., qu, we assign a fixed weight ω i . In fact, it holds that w " where 1 C is the indicator function of the set C. In this work, we obtain the weights w˚P R q and ω˚P R L that minimize a threshold m 0 describing the expected number of required measurements for Models 1 and 2, respectively. Our approach is to find a suitable upper-bound for m 0 . The bound is not necessarily tight but leads to closed-form expressions for w P R q and ω P R L in Models 1 and 2, respectively. Availability of the parameters (p b 's and α i 's) associated with Models 1 and 2 is feasible in certain practical scenarios. For instance, they might be obtained from previous arXiv:1809.03005v1 [cs.IT] 9 Sep 2018 signal estimation in a dynamic scenario or from the statistic of previous training data. We further show that how Model 1 and 2 apply to brain functional magnetic resonance imaging ( fMRI) reconstruction and direction of arrival (DOA) estimation, respectively (see sections I-B,III). It often occurs in practice that the parameters p b 's and α i 's are not accurately accessible, but only an approximation is available. Our last contribution in this work is about this challenge. In fact, we show that our method of obtaining the optimal weights wå nd ω˚is robust to small changes of p b and α i as long as p b Ç 1 10 and α i Ç 1 10 . It is worth mentioning that in this work, two structures are used to construct the measurement matrix A mˆn :randomly subsampled discrete Fourier transform (DFT) 2 and i.i.d. standard normal distribution. The former lies in the structure of DOA estimation (see subsection I-B) while we use the latter for fMRI reconstruction. Although, at first look, might be different, these two distributions have a common property making them usable in our conic geometry approach [11,Theorem 1]; in fact, for both, null space of A mˆn is drawn uniformly at random from Grassmann manifold of subspaces of dimension n´m.

A. Related works and Key Differences
CS in presence of prior information has been studied in different signal models. While a large part of research (see for example [12]- [17]) deals with deterministic signal models, only a few works (see [18]- [20]) have investigated random signal models with Bayesian information. In the deterministic model, the ground truth signal has intersected with a few sets which called support estimates. The contributing level of each set to the support is available to the experimenter [12], [13]. This exact situation is investigated in [17]. They propose a non-uniform model for capturing deterministic prior information. The work [19] considers a probabilistic model where there is a continuous shape function describing the probability of contributing each element to the support. The authors obtain an upper-bound for failure probability of weighted 1 minimization. Their approach is based on calculating the internal and external angles of a weighted cross polytope. With a different approach, [20] has investigated a discrete measure for Bayesian information (a special case of Model 1 with k " 1). They relate the weights of weighted Fig. 1. Schematic diagram of DOA estimation of far-field sources. The angular half-space is divided into q " 30 angular clusters tV b u q b"1 with equal length. The associated parameters corresponding to Model 2 are L " 5,

B. Broadband DOA Estimation
Suppose that s far-field broadband signals tx i pt, f qu s i"1 in the frequency range f P rf L , f H s incident on an q-element uniform linear array (ULA). The received signal in m sensors at time t and l-th frequency bin can be expressed as: where a i pf l q :" r1, ..., e´j 2πf l pq´1qd sinpθiq c s T is the steering vector, c is the propagation velocity, d is the inter-sensor spacing and ept, f l q is a Gaussian noise term with variance σ 2 . In practice, one has to take several snapshots typt, f l qu k t"1 . This temporal redundancy is crucial in practice since the array size is limited due to physical constraints [21], [22]. Consequently, one may write Y pf l q :" ryp1, f l q, ..., ypk, f l qs " Apf l qXpf l q`Epf l q P C mˆk , where Xpf l q " rxp1, f l q, ..., xpq, f l qs P C qˆk , Apf l q " ra 1 pf l q, ..., a q pf l qs P C mˆq and Epf l q P C mˆk is defined similarly. If the sources is time-invariant over the period of snapshotting, then for all t " 1, ..., k the non-zero dominant peaks in xpt, f l q occur at the same locations corresponding to the ground truth DOAs. Hence, DOA estimation can be cast as recovering a joint sparse signal Xpf l q from Y pf l q. In addition, it is realistic for a radar engineer to know the probability of appearing the ground truth DOAs in some angular bands [23] (see Figure 1 for a schematic model of this scenario). Notation. Throughout, scalars are denoted by lowercase letters, vectors by lowercase boldface letters, and matrices by uppercase boldface letters. The ith element of a vector x is shown either by xpiq or x i . C˝denotes the polar of a cone C. We show sets (e.g. B) by calligraphic uppercase letters. B is used to represent the complement t1, ..., nuzB of a set B Ă t1, ..., nu. g P R n is a vector with standard normal distribution. The performance of an image x X P R n1ˆn2 is measured by peak signal to noise ratio (PSNR) given by: where X is the true image and }¨} 8 returns the maximum intensity.

II. MAIN RESULTS
Before stating our main results, we introduce two concepts from conic integral geometry that are used in our analysis.
Descent cone: Let x be a vector and f be its structure inducing convex function. Then, the set of decent directions forms a convex set defined by: Dpf, xq " ď tě0 tz P C n : f px`tzq ď f pxqu¨ (8) Statistical dimension: It is shown in [11] that statistical dimension of the above decent cone defined by δpDpf, xqq :" Edist 2 pg, Dpf, xq˝q, specifies the required number of measurements (i.e. m) that needs for perfect recovery. Define the expected number of measurements needed for P 1,2,w as s m :" E x δpDp}¨} 1,2,w , xqq.
Then, we call the weights that minimize s m optimal in sense of expected sample complexity. In the following propositions, we obtain closed-form solutions for optimal weights in case that x satisfies Model 1 and 2. The proofs are provided in Appendix A. Proposition 1. Let x P C n satisfies Model 1 with parameter p " rp 1 , ..., p q s T . Then, the optimal weights w˚" rw1 , ..., wq s T in P 1,2,w are obtained by solving the following equations simultaneously: Proposition 2. Let x P C n be decomposed into q blocks V b of equal size k. Assume that there exist L independent Bayesian estimates tP i u L i"1 of bsupppxq with parameter α " rα 1 , ..., α L s T . Then, the optimal weights ω˚" rω1 , ..., ωLs T in Model 2 are obtained by In applications, it is of important practical value to know how the inaccuracies of p b and α i in Model 1 and 2, respectively, affect the optimal weights. The following theorem is about this challenge. Theorem 1. Assume that p and p 1 be the true and approximate estimate of bsupppxq, respectively. Let w and w 1 be the Fig. 2. The middle image compares the success rate of P 0 1,2,w˚w ith that in P 0 1,2 when there exists a prior distribution p P R q (depicted in the left image) about the block indexes t1, ..., 50u. The optimal weight w˚is obtained by (12). The right image represents cpk b , p b q in Theorem 11 for different k b 's. In the right image, recovery with P η 1,2 is shown in the dashed line while recovery using P η 1,2,w˚i s shown in the solid line. The prior distribution p P R q that is used to calculate w˚via (12) is depicted in the left image. The used parameters are k b " 1 for all b and q " 128 2 " 16384. Selected true frames are also displayed in the right image, with green markers showing where these frames correspond to the PSNR of P η 1,2 .
optimal weights corresponding to p, and p 1 . Then, there exists a constant cpk b , p b q such that where γpa, zq :" ş 8 z u a´1 e´udu is incomplete gamma function and hp¨q : p b Ñ w b is the nonlinear function in (12).
From the above theorem and the right image of Figure 2, one can infer that the method of obtaining wb and ωi is robust to slight changes of p b and α i as long as p b Ç 1 10 and α i Ç 1 10 .

III. SIMULATIONS
In the first experiment, we construct a random blocksparse x P R 250 , whose building blocks V b have equal size k " 5. The probability of activating each block is taken from the vector p depicted in the left image of Figure 2 . Then, this signal is observed through a measurement matrix A P R mˆ250 , of which the elements are drawn from i.i.d. standard normal distribution. We obtain the optimal weights w˚corresponding to p by solving (12). In the middle image of Figure 2, we plot the success rate as a function of m. For each m, we average over 100 realizations of A and x. As turns out, P 1,2,w˚n eeds less measurements than P 1,2 . Fig. 4. Broadband DOA estimation for 10 sources using an m " 15 sensor ULA. The associated parameters are k " 10, q " 100, d " 5, c " 3ˆ10 8 , and σ " 1. The accuracies corresponding to three Bayesian estimates are α 1 " 4 5 , α 2 " 2 3 , and α 3 " 0. The top and bottom images correspond to the recovery using P η 1,2 and P η 1,2,Dω˚, respectively, where D :" r1 P 1 , 1 P 2 , 1 P 3 s and ω˚is calculated using (13).
In the second experiment, we test brain fMRI reconstruction of n 1ˆn2 images when there exists a two dimensional prior distribution (the left image in Figure 3). We utilize 27 frames of a MRI video scan, each of which contains 128ˆ128 pixels (k b " 1 @b and q " 128 2 " 16384), and consider only the gray scale components. The measurement matrix A P R 9830ˆ16384 for each frame is chosen from an i.i.d. standard normal distribution. Also, the additive noise e is chosen to be i.i.d. Gaussian with variance 1 100 . We implement two recovery programs P η 1,2 and P η 1,2,w˚( w˚is obtained by (12)) using TFOCS package [24]. From the right image of Figure 3, one can see that, for a fixed m " 9830 and with our optimal weighting strategy, P η 1,2,w˚c learly acts better than P η 1,2 (leading to PSNR values about 30dB higher than that in P η 1,2 ). In the last experiment, we test DOA estimation using broadband signals (see Subsection I-B and Figure 1). The angular half-space´r90˝, 90˝q is divided into q " 100 angular clusters (grids). For each frequency bin f l P r0, 5s GHz, we take k " 10 snapshots. We assume that there exist three (L " 3) Bayesian estimates P i with accuracies α 1 " 4 5 , α 2 " 2 3 , and α 3 " 0. Also, we use m " 15 sensors for recovering s " 10 ground truth sources and implement the optimization problems P η 1,2 and P η 1,2,w˚w here w˚" ř L i"1 ωi 1 Pi . ωi 's are obtained using (13). As turns out from Figure 4, while P η 1,2 detects many nonexisting sources with spurious DOAs, P η 1,2,w˚l ocates the ground truth sources correctly. This in turn suggests that our optimal weighting strategy considerably decreases the required number of sensors.
APPENDIX A. Proof of Propositions 1 and 2 Proof. Since δpDp}¨} 1,2,w , xqq only depends on bsupppxq and w, we show it by δpDpB, wqq. It holds that: We proceed by using an upper-bound for δpDpB, wqq a special case of which is obtained in [17,Lemma 1]. Lemma 1. The statistical dimension of descent cone of any vector x nˆ1 with bsupppxq " B satisfies: with φ B pz, kq :" Moreover, the minimum achieved at a unique t ě 0. The inequality is in fact equality in the asymptotic case (q Ñ 8) [17, Proposition 3].
Thus, by the aid of Lemma 1, it holds that By using p b :" Ppb P bsupppxqq " ř BQb q B and 1´p b " ř BSb q B , we have:

16)
Moreover, by the inequality I and (4), one may write where we used α i " PpP i Ă bsupppxqq " ř BĄPi q B and 1´α i " PpP i Ć bsupppxqq " p ř BČPi q B q. Since scaling is effectless on P 1,2,w˚, by minimizing the expressions in parentheses in (16) and (17) with respect to w b and ω i , one can get to the expressions (12) and (13).