{'fraction_non_alphanumeric': 0.04749830201494227, 'fraction_numerical': 0.01136517998641612, 'mean_word_length': 4.233649289099526, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Detection-Directed Estimation Approach for Noisy Compressive Sensing

In this paper, we investigate a Bayesian sparse reconstruction algorithm called compressive sensing via Bayesian support detection (CS-BSD). This algorithm is quite robust against measurement noise and achieves the performance of a minimum mean square error (MMSE) estimator that has support knowledge beyond a certain SNR threshold. The key idea behind CS-BSD is that reconstruction takes a detection-directed estimation structure consisting of two parts: support detection and signal value estimation. Belief propagation (BP) and a Bayesian hypothesis test perform support detection, and an MMSE estimator finds the signal values belonging to the support set. CS-BSD converges faster than other BP-based algorithms, and it can be converted to a parallel architecture to become much faster. Numerical results are provided to verify the superiority of CS-BSD compared to recent algorithms.

elements are sparsely non-zeros, called sparse signal. Then, the NCS decoder observes a measurement vector z = [z 1 , ..., z M ] ∈ R M , given as where x 0 ∈ R N is a deterministic sparse signal; Φ ∈ R M ×N is a sensing matrix whose columns represent a possibly overcomplete basis, i.e., rank(Φ) ≤ M , where M < N ;and n ∈ R M is an additive noise vector generated by a certain distribution.
Alternatively, Bayesian approaches to NCS have received attention [7]- [15]. This type of approach offers powerful mitigation of noise effects by using many existing statistical signal processing techniques and several statistical signal-noise models. In these approaches, the reconstruction problem is described as the maximum a posteriori (MAP) estimation problem as follows: where Gaussian noise is assumed and f (·) is a probability density function.
The SBL algorithm iteratively determines the posterior density of the signal on basis of a three-layer hierarchical prior model, so the prior density is a function of certain parameters. The algorithm estimates the parameters of the prior using expectation maximization (EM) and applies these parameters to finding world applications. In addition, these approaches can be made more attractive if they are applied in in conjunction with belief propagation (BP). BP replaces the reconstruction process by iterative messagepassing processes. This replacement reduces the reconstruction complexity to the O(N log N ) order. Baron et al. for the first time proposed the use of sparse matrices to the NCS setup and developed a BP-based algorithm, called CS-BP [12], [13]. CS-BP iteratively updates the signal posterior from the twostate Gaussian mixture prior via the message-passing algorithm, where the messages are the probability densities of the signal elements. In [14], Tan et al. proposed another BP-based algorithm called, BP-SBL. They applied BP to the SBL-framework in [9] to reduce the complexity of the EM algorithm. Most recently, Akcakaya et al. devised SuPrEM using an idea similar to BP-SBL, but in a different framework [15] which is based on Gaussian scale mixture [16] with a specific type of prior called the Jeffreys' prior [17]. In addition, the authors restrict the story of SuPrEM to a class of sensing matrices, called low-density frames, in which the matrices have fixed column and row weights.
In this paper, we propose a sparse reconstruction algorithm based on the Bayesian approach and the use of sparse matrices. We call our algorithm as Compressive sensing via bayesian support detection (CS-BSD). CS-BSD has the following properties: 1) Robustness against the measurement noise effects.
2) Ability to perform as the minimum mean square error (MMSE) estimator that has knowledge of the support set.
CS-BSD has a detection-directed (DD) estimation structure which consists of signal support detection and signal value estimation, as shown in Fig.1. We consider the common procedure of first using the measurements at hand to detect the signal's support set. This detected support set is then used in the model of the sparse signal, and the value estimator is built as if the detected support set is in fact the correct set. The support detection component consists of a combination of the Bayesian hypothesis test (BHT) and BP, and signal value estimation using the detected support set is achieved via an MMSE estimator. CS-BSD iterates the detection and estimation processes until the constraint in (3) is met.
The DD estimation methodology was investigated in [21] for estimation of noisy signals and have been widely applied to wireless communication systems [22], [23]. For CS, the methodology was first reported in [24], [25]; we tailor the methodology to the NCS problem by refining that work. The complexity of CS-BSD is O(N log N + KM ) whereas that of the other BP-based algorithm is O(N log N ) because CS-BSD includes the cost of MMSE in addition to that of BP. However, CS-BSD converges faster than the other BP-based algorithms; thus, its computational cost is lower in practice. In addition, CS-BSD can be much faster by converting to an parallel architecture.
The rest of the paper is organized as follows. Section II introduces the sparse sensing matrix, the prior model for our system model. The details of CS-BSD are given in Section III. A few practical issues are discussed in Section IV. We compare the numerical results of CS-BSD to the other recent CS algorithms in Section V. Section VI concludes the paper.

A. Sparse Sensing Matrix Φ
For signal sensing, we employ sparse-Bernoulli matrices Φ ∈ {0, 1, −1} M ×N , which have been successfully used in CS recently [12]- [14]. In the matrix, sparsely nonzero elements are equiprobably equal to 1 or −1. We set the sparsity of Φ using the fixed column weight L. Because the column weight rather than the row weight is fixed, all elements of x 0 have an even chance of being sensed. In addition, the fixed column weight unifies the energy of the basis of the measurement space spanned by the column vectors of Φ.
With the sparse-Bernoulli matrix, the linear system z = Φx 0 + n can be represented over a bipartite for all j ∈ C, respectively. Note that |N V (i)| = L for all i ∈ V under our assumption regarding Φ. Fig.2 depicts a simple example of the graphical representation corresponding to N = 6, M = 4, L = 2.

B. Prior Model
We limit our discussion to the random vector x whose elements are i.i.d. random variables. This assumption is commonly used in many papers [7]- [15]. We characterize the signal sparsity in a probabilistic manner, called sparsity rate. The sparsity rate q is defined as q := Pr{x i = 0} for all i ∈ V. Namely, each signal element independently belongs to the signal support set with the rate q. The supportiveness of each signal element is represented by a state variable s i , defined as Hence, we model the prior density of x i using a spike-and-slab model originating in a two-state mixture density as follows: where δ(x) indicates a Dirac distribution having nonzero value between x ∈ [0−, 0+] and δ(x)dx = 1.
In the prior density, we use Gaussian density N (x; 0, σ 2 x ) for f x (x|s = 1) although it includes the probability mass at x i = 0. The reason is the probability mass at x i = 0 is very small and Gaussian densities are mathematically tractable. In addition, we drop the index i from the prior density under the assumption of i.i.d. elements. The spike-and-slab prior has been widely employed in Bayesian inference problems [26]- [28] and was recently applied to CS [11] as well.

III. PROPOSED ALGORITHM
In this section, we discuss the details of the proposed algorithm based on the DD estimation structure.
The proposed algorithm, CS-BSD, is an iterative algorithm that repeats the support detection and signal value estimation processes until E Φx − z 2 ≤ is met.

A. Detection of Support Set
The decoder detects the signal support in each element unit. Namely, the supportive state of each signal element is detected independently and converted to the support set information for the signal. First, the following simple hypothesis test can be considered for the state detection of x i : where H 0 and H 1 are two possible hypotheses. If we marginalize over s i , the left hand side of (6) becomes where The right hand side of (6) is where From (7) and (9), the hypothesis test in (6) is refined as Here where the posterior density, f xi (x|z, s i = 1), is Gaussian because the signal and noise are assumed to be Gaussian (see p.326 in [29]). The term of Pr{x i = 0|z, s i = 1} in right hand side of (12) is caused by the use of Gaussian density N (x; 0, σ 2 x ) for the prior of nonzero x i . Because the variance of f xi (x|z, s i = 1) is a function of the noise variance, the probability Pr{x i = 0|z, s i = 1} is very small, and it approaches zero as the SNR increases. Therefore, we suggest setting the threshold of the hypothesis test in (12) to 1. This implies that the hypothesis test can detect the supportive state of the signal elements with a high success probability if SNR is sufficiently high.
We now describe how to obtain the probability ratio, Pr{si=0|z} Pr{si=1|z} . By factorizing over x i , the ratio becomes where f xi (x|z) denotes the posterior density of x i given z. The signal elements are not i.i.d. anymore given z. In (14), Pr{s i |z, x i } = Pr{s i |x i } holds true since the measurements z does not provide any additional information on the state given x i . Using the Bayesian rule and the prior information, we finally obtain the hypothesis test as the following form: Since we know the prior of the state P r{s} from the sparsity rate, i.e., Pr{s = 1} = q, we can move the prior term to the right side, and then treat it as a threshold for the hypothesis test γ := Pr{s=1} Pr{s=0} = q (1−q) . Therefore, the state of each elements can be sensed from the corresponding posterior and prior densities.
Definition 1 (BHT for state detection): Let s i denote the detected state of x i ; f xi (x|z) indicates the posterior density of x i ;, and f x (x|s) denotes the conditional prior density of a signal element given the state. Then, state detection for all i ∈ V is performed by choosing the hypothesis that result from where

B. Belief Propagation for Posterior Update
The posterior density used for the BHT is obtained and updated at every iteration via BP. Our BP process is similar to that in [12], [13] and was independently devised from [14], [15]. Distinctively, our BP process uses the information on the noise distributions f nj (n) = N (n; 0, σ 2 n ) under the i.i.d. zero-mean Gaussian noise assumption.
Using Bayesian rule, we can represent the posterior density of x i in the form of Posterior = Prior × Likelihood Evidence , given as If the sensing matrix Φ is sufficiently sparse such that the corresponding bipartite graph is tree-like, we postulate that the elements of z associated with x i are independent of each other given x i [19]. Under the tree-like assumption, we can decompose the likelihood density f z (z|x i ) to the product of densities: We call each decomposition of the likelihood, f zj (z|x i ) the measurement density. Theorem 1 below demonstrates that the measurement density can be composed of the densities of the associated signal elements.
Theorem 1 (Measurement density in BP): The measurement density f zj (z|x i ) is expressed as the linear convolution of all the associated distributions of the signal elements and the corresponding noise distribution f nj (n) as follows: for all (i, j) ∈ E, where ⊗ and are the operator for linear convolution and the linear convolution of a sequence of functions, respectively Proof : See Appendix A.
Therefore, the essence of the BP-process is to update the signal and measurement densities by exchanging probability density messages, associated with the neighboring relation in the bipartite graph.
Let a i→j denote the message from the i-th signal element to the j-th measurement element, called the signal message; b j→i is the message from the j-th measurement element to the i-th signal element, called the measurement message. The signal message is an approximation of the density of the signal element, i.e., a i→j ≈ f xi (x|z) and it is obtained from (19) simply by replacing the measurement density with the measurement message of the previous iteration. Note that in BP-process the message coming from the j-th measurement is excluded in the calculation of a i→j . Thus, the signal message at the l-th iteration is expressed as for all (i, j) ∈ E, where η[·] is the normalization function to make a i→j dx = 1. Similarly, the measurement message approximates the measurement density, i.e., b j→i ≈ f zj (z|x i ), and it is obtained from the expression of (20) by replacing the associated densities of signal elements f xk (x) with the signal messages for the purpose of iteration , that is, The convolution operations in (22) can be efficiently computed by using the Fast fourier transform (FFT). Therefore, we express for the measurement message calculation as where F ∈ C Nd×Nd denotes a Fourier matrix of size N d . In fact, the use of the FFT brings a small calculation gap between this result and that of (20) since the FFT-based calculation performs a circular convolution that produces output having a heavy tail, as shown in Fig.3. The heaviness increases as the corresponding row weights in Φ increase. However, the difference is can be ignored, especially when the densities are bell-shaped distributions.
Finally, the update of the posterior density of x i at the l-th iteration is provided as given in Definition 2.
Definition 2 (Posterior update in BP): Let b l j→i denote a measurement message at the l-th iteration for all (i, j) ∈ E. Then, the posterior density of x i at the l-th iteration is calculated by where η[·] is the normalization function that makes f x l i (x|z)dx = 1.

C. Detection-Directed Estimation of Signal Values
We now describe signal value estimation based on the DD estimation structure. The DD estimator is basically an estimator that determines how to act on the input data directed by the information from the detector. In CS-BSD, the detector provides the support information s l , and the value estimator then finds the signal values as if the detected support set is the correct set at each iteration. That is, where the estimator decides that x l i = 0 for all i ∈ V : s l i = 0. From the argument in (12), the DD estimate converges to the true signal x 0 since the detected support set becomes the correct set as SNR and the number of iterations l increases. This DD methodology makes no general claim regarding optimality of the solution; however, it is common and often successful. Let x l supp ∈ R || s l ||0 denote a random vector consisting of the elements with s l i = 1. Then, the problem in (25) is reduced to Since x l supp and the noise elements are assumed to be zero-mean i.i.d. Gaussian random variables with variance σ 2 x and σ 2 n respectively, the MAP estimation in (26) is recast as where Φ l supp denotes a submatrix of Φ corresponding to i ∈ V : s i = 1. In addition, the MAP and MMSE estimates are identical, assuming the signal and noise are Gaussian (see p.358 in [29]). Therefore, the estimate x l supp can be obtained by the MMSE estimator To combine the support information s l and the estimated values x l supp , we define an index set U l := {1, ..., s l 0 } corresponding to the elements x l supp = [x supp,1 , ..., x supp, s l 0 ] and a bijective mapping function h : {i ∈ V| s l i = 1} → U l . Then, the reconstruction at each iteration is readily obtained from for all i ∈ V. CS-BSD is summarized in Algorithm 1.

A. Complexity
We implement the BP process in CS-BSD based on the sampled-message approach in [13]. The density messages a i→j , b j→i are vectors of size N d where N d is chosen to be power of two for efficient use of FFT. Next, we analyze the complexity of CS-BSD by considering each part seperately.  (16), the decoder requires O(N d ) flops to calculate a likelihood ratio. The cost for the hypothesis test is much smaller than that of BP; therefore, it is ignored.  [20] such that its complexity is O (N log N + KM log N ). If we fix the number of iteration N iter empirically, we can remove the MMSE operation from the iterations. In that case, the complexity is reduced to O (N log N + KM ).

B. Parallelization of Belief Propagation
The BP process for finding the posterior finding can be implemented using a parallel architecture.
Indeed, many parallelized BP algorithms, with applications to LDPC codes, have demonstrated superior performance in [31]- [33]. The graph representation of the sparse sensing matrix shows that the dependencies of the message calculations for any signal elements (or measurement elements) depend only upon the corresponding measurement elements (or signal elements). This allows all messages in BP to be computed in a parallel manner. Therefore, implementing BP on a parallel architecture for BP yields low power consumption, high-speed decoding, and simple logic [31].

V. NUMERICAL RESULTS
We demonstrate the advantages of CS-BSD using simulation results in several different settings. To show its average performance, we take 200 Monte Carlo trials for each point in the simulation. In each trial, we generate the deterministic sparse signal x 0 with N = 1024 and σ x = 10 whose values are represented with finite precision. The finite precision is provided by 6-bit quantization such that each signal value has 64 levels. This assumption of finite precision for the signal values is reasonable in terms of digital signal processing and implementation. In addition, we restrict the magnitude level of the signal elements to |x i | ≤ 3σ x for the same reason. We define the SNR as SNR : = 10 log 10 E Φx 2 2 M σ 2 n dB (30) and M/N as the undersampling ratio for signal acquisition.

A. SER Performance of Support Detector
To determine the performance of the support detector in CS-BSD, we defined the state error rate (SER) as: where s 0,i is the state variable corresponding to the true signal value x 0,i . We simulate the SER performance as a function of the SNR for a variety of undersampling ratio M/N . In this simulation, we set q = 0.05, N d = 64, and L = 4. In addition. we compare the SER performance to a theoretical limit on the support recovery given by Fletcher et al. [35]. They found a necessary condition for maximum-likelihood (ML) estimation to asymptotically recover the support set if the sensing matrix has i.i.d. Gaussian entries.
The ML estimation is described as where the signal sparsity x 0 is assumed to be known, J ⊆ V is a subset of the index set of the signal, and P J z denotes the orthogonal projection of z onto the subspace spanned by columns of Φ corresponding to J . Namely, the ML estimate is a subset of V such that the subspace spanned by the corresponding columns of Φ contain the maximum energy of z. We rewrite the necessary condition in terms of SNR such that where minimum-to-average ratio (MAR) is defined as In this comparison, we used 200 Monte Carlo trials to find the average SNR limit , i.e., SNR limit := avg[SNR limit ]. In Fig.4, the SER curves show a waterfall behavior; the curves decline rapidly to less than 10 −5 beyond a certain threshold SNR. This behavior supports the argument in (12)

B. MSE Performance Comparison
We consider the reconstruction performance in terms of normalized means square error (MSE), which is defined as We compare our algorithm to several recent CS reconstruction algorithms: 1) CS-BP [12], [13], 2) L1-DS via linear programming [5], 3) Bayesian CS (BCS) [9], 4) CoSaMP [34], and 5) SuPrEM (reweighted version) [15]. For BCS and SuPrEM, we obtained the source code from each author's webpage; for CoSaMP we used Stephen Becker's code (available at http://www.ugcs.caltech.edu/ ∼ srbecker/algorithms. shtml). L1-DS is provided by the L1-MAGIC package (available at http://users.ece.gatech.edu/ ∼ justin/ December 21, 2013 DRAFT l1magic/).We implemented CS-BP algorithm by using the sampled-message approach and upgrading the original algorithm to use the noise information. For CS-BP, we used the sparse-Bernoulli sensing matrix with L = 4; for SuPrEM, we use a sensing matrix generated from a low-density frame [15] with the same parameters (N , M , L). L1-DS, CoSaMP and BCS were used with a Gaussian sensing matrix having the same column energy as the sparse-Bernoulli matrix, for fairness, i.e., φ j,Gaussian  Table I, with respect to thier complexity, type of sensing matrix, prior type, and algorithm type.
where Tr[·] denote the matrix trace operation. Beyond SNR=31 dB, since the SER of CS-BSD is almost error-free, the MSE performance achieves MSE * at M/N ≥ 0.5. Surprisingly, this result is superior to that of the l 1 norm based approach, which is known as an optimal algorithm in the noiseless case. The gap between the two algorithms is caused by the reconstruction error over the non-supporting elements.
CS-BSD completely removes the error from the non-supporting elements whereas the l 1 norm based approach leaves a certain amount of the reconstruction error on the non-supporting elements.
In the low SNR regime, it is noteworthy that CS-BSD works well although the proposed algorithm was originally targeted at a reasonable system having high SNR. x k = x NC(j),g(k) for all k ∈ N C (j).
By marginalizing over n j to f zj (z|x i ), we obtain where f nj (n|x i ) = f nj (n) since n is independent of x. By further marginalizing over elements of x NC(j) , we rewrite the expression in (38) as where we assume x i = x NC(j),1 without loss of generality. In addition, f zj (z|x NC(j) , n j ) = δ(z − z j ) holds true since knowing x NC(j) is equivalent to knowing (Φx) row(j) ; thus, there is no uncertainty in z j = (Φx) row(j) + n j . Since the elements of x are assumed be independent, we replace (39) with the product of the probability densities.
The expression in (40) can be represented by a sequence of convolutions of probability densities, as given in (20).