Efficient Sparse Bayesian Learning via Gibbs Sampling - IEEE Xplore
EFFICIENT SPARSE BAYESIAN LEARNING VIA GIBBS SAMPLING Xing Tan and Jian Li ∗
Peter Stoica
Dept. of Electrical and Computer Eng., Univ. of Florida, Gainesville, FL, USA, email: tanxing@uÀ.edu, [email protected]À.edu
Dept. of Information Technology, Uppsala University, Uppsala, Sweden, Email: [email protected]
ABSTRACT Sparse Bayesian learning (SBL) has been used as a signal recovery algorithm for compressed sensing. It has been shown that SBL is easy to use and can recover sparse signals more accurately than the well-known Basis Pursuit (BP) algorithm. However, the computational complexity of SBL is quite high, which limits its use in large-scale problems. We propose herein an ef¿cient Gibbs sampling approach, referred to as GS-SBL, for compressed sensing. Numerical examples show that GS-SBL can be faster and perform better than the existing SBL approaches.
collecting N data samples and computing N transform coef¿cients in which only the K largest ones are retained and encoded. Fortunately, with the recent advent of compressed sensing techniques (see, e.g., [1]), we can recover the original signal s from much fewer samples. In compressed sensing, the signal is measured through a linear process that can be modeled as projecting the signal s onto a set of vectors {m 1 , m 2 , . . . , m M } to obtain the following samples: ym = mm s +em , m = 1, 2, . . . , M , where em is the measurement noise. If we de¿ne the data vector y = [y1 , y2 , . . . , yM ] , the measurement matrix M = [m 1 , m 2 , . . . , m M ] and the noise vector e = [e1 , e2 , . . . , eM ] , then y can be expressed as y = M s + e. Using the fact that s = Bx, we obtain the following data model:
Index Terms— Compressed Sensing, Sparse Bayesian Learning, Gibbs Sampling
(1)
y = Ax + e, 1. INTRODUCTION We denote vectors by boldface lowercase letters and matrices by boldface uppercase letters. Other notations that we use throughout this paper are given in Table 1.
where A M B is an M × N matrix. To obtain a sparse solution of (1), an l1 -norm based optimization approach, referred to as Basis Pursuit (BP), was proposed in [2]: ˆ = arg min ||x||1 x x
(·) (·)∗ · q μn Σk,l diag{v}
Table 1. Notations used in the text transpose of a vector or a matrix conjugate transpose of a vector or a matrix q-norm. 2-norm when the subscript is omitted. the nth element of the vector μ the (k, l)th element of the matrix Σ A diagonal matrix whose diagonal entries are the elements of the vector v
(2)
where δ is a user parameter. An alternative l1 -norm based optimization is ˆ = arg min ||y − Ax||2 + λ||x||1 . x x
In the past few decades, transform coding has been widely used for data compression. Denote the signal to be compressed as an N -dimensional real-valued vector s. Let B be an orthogonal basis matrix which satis¿es B B = BB = I and de¿ne the transform coef¿cient vector x as x = B s. If x can be approximated ˆ (i.e., x ˆ contains no more than K non-zero by a K-sparse vector x elements and usually K N ), then we only need to encode the ˆ instead of all N elK non-zero elements and their positions in x ements of s. Transform coding is inef¿cient, because it requires ∗ This material is based on research sponsored in part by the of¿ce of Naval Research (ONR) under grants No. N00014-09-1-0211 and No. N00014-07-1-0193, the Army Research Of¿ce (ARO) under Grant No. W911NF-07-1-0450, the National Science Foundation (NSF) under Grant No. ECCS-0729727, the Komen Breast Cancer Foundation under grant No. BCTR0707587, the Swedish Research Council (VR) and the European Research Council (ERC). Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the United States Government.
978-1-4244-4296-6/10/$25.00 ©2010 IEEE
subject to ||y − Ax||2 < δ,
3634
(3)
From a Bayesian point of view, (3) is equivalent to the maximum a posteriori approach for estimating x in (1), assuming that e is a normal random vector and x has a Laplacian a priori distribution. Many ef¿cient methods for solving (3) have been proposed recently, including the iterative shrinkage/thresholding (IST, see, e.g., [3]) and gradient projection for sparse reconstruction (GPSR, e.g., [4]). All these methods depend on the fast computation of the product of A and an arbitrary vector. If the complexity of the matrix vector product is O(L ), then IST and GPSR also have O(L ) complexity. Hence, when A is a fast Fourier Transform (FFT) matrix, a Discrete Cosine Transform (DCT) matrix, a Wavelet Transform matrix or a sparse matrix, then the complexities of IST and GPSR are nearly linear to N . Therefore, they are highly ef¿cient in solving largescale compressed sensing problems. Another ef¿cient sparse signal recovery method is the so-called Compressive Sampling Matching Pursuit (CoSaMP, [5]), which also has O(L ) complexity. One can also obtain a sparse estimate of x via sparse Bayesian learning (SBL), which was originally proposed in the machine learning literature [6]. The model used in SBL is a three-stage hierarchical Bayesian model. The parameters in the Bayesian model were obtained by using the maximum a posteriori (MAP) criterion. However, the MAP problem cannot be solved by setting the derivative
ICASSP 2010
of the cost function to zero, since the derivative is highly nonlinear. Hence, the author proposed two iterative approaches: a direct differentiation approach and an Expectation Maximization (EM) approach (see Appendix A of [6]). The vector x can be estimated by computing its a posteriori mean during the iterative process. A greedy SBL method, referred to as the Fast Marginal Likelihood Maximization algorithm was proposed in [7] and later on applied to compressed sensing in [8] under the name Bayesian Compressive Sensing (BCS). In [8], the authors show that in some cases, the greedy SBL gives a more accurate estimate of the signal than BP algorithm in terms of ˆ ||2 x dB. the reconstruction error de¿ned by err = 10 log10 ||s−B ||s||2 A major advantage of SBL over the l1 -norm based optimization is that the former provides not only a point estimate of x, but also the a posteriori covariance matrix of x, which measures the accuracy of the point estimate. However, both the EM-based SBL and the direct differentiation based SBL have high computational complexities of O(M 2 N ) (see [6]). The Fast Marginal Likelihood Maximization approach has a complexity of O(M N +K(L +K 2 )) (see Appendix A in [7]). This method is indeed fast when K is small. However, when K is large, it can be rather slow. Herein, we propose a novel Gibbs sampling based SBL approach which computes the a posteriori pdf of x in O(L ) Àops. We refer to this approach as GS-SBL. It is known that, in general, the Markov Chain Monte Carlo (MCMC) methods are slower than other Bayesian methods. However, with an ef¿cient sampling step, GSSBL only has a complexity of O(L ). Although we still need to generate a few hundred samples to obtain a good estimate for x, in large-scale compressed sensing problems the total time needed by GS-SBL is still much less than that required by the original SBL proposed in [6]. Furthermore, numerical examples show that GSSBL can outperform the existing SBL approaches in both small- and large-scale compressed sensing problems.
Consider the following Bayesian model:
f (p, η|y, x, a, b) ∝ f (y|x, p, η)f (x|p)f (p|a, b) N −(a+3/2) − 1 (b+ x2n ) 1 (y −Ax2 ) −(M/2) − 2η p 2 e pn e n . (5) ∝η n=1
and f (x|y, p, η, a, b) ∝ f (y|x, p, η)f (x|p)f (p|a, b) Σ−1 (x−μ)
(8)
From (6), x has a conditional normal distribution: x|y, p, η, a, b ∼ N (μ, Σ).
(9)
Using (7), (8) and (9), we can derive the Gibbs sampling algorithm (see, e.g., [9]). Assume that an initial estimate x(0) is given (see the comments on the initial estimator at the end of this section). In the tth iteration, we perform the following two steps: (t)
1. We ¿rst draw a sample of pn from the conditional pdf f (pn |y, x(t−1) , a, b) for n = 1, 2, . . . , N , and draw a sample of η (t) from the conditional pdf f (η|y, x(t−1) , a, b). By using Equation (8), we obtain ⎛
2 ⎞ (t−1) x n 1 ⎟ ⎜ , b+ p(t) ⎠ , n = 1, 2, . . . , N, n ∼ IG ⎝ a + 2 2 (10) and by using Equation (7), we obtain y − Ax(t−1) 2 M η (t) ∼ IG − 1, . 2 2
(11)
(12)
where (4)
where p = [p1 , p2 , . . . , pN ] , P = diag{p}, N (μ, Σ) represents the multivariate normal distribution with mean μ and covariance matrix Σ and IG(a, b) represents the inverse gamma distribution. We assume the inverse gamma prior for pn because we know that p is sparse and the inverse gamma prior promotes sparsity in the estimate of p, when b is small. But for η we do not have any a priori information, so it is proper to assume that f (η) ∝ 1. Based on the Bayesian model (4), we can readily obtain:
1 x2 pn |y, x, a, b ∼ IG a + , b + n . 2 2
x(t) ∼ N (μ(t) , Σ(t) ),
y|x, η ∼ N (Ax, ηI) x|p ∼ N (0, P ), f (η) ∝ 1 pn |a, b ∼ IG(a, b),
1
and
2. Then we draw a sample of x(t) from the conditional pdf of f (x|y, p(t) , η (t) , a, b). De¿ne the matrix P (t) diag{p(t) }. By Equation (9), we obtain
2. AN EFFICIENT GIBBS SAMPLING APPROACH
∝ e− 2 (x−μ)
It follows from (5) that p and η are conditionally independent and that they have the following distributions, respectively, M y − Ax2 η|y, x, a, b ∼ IG − 1, , (7) 2 2
,
(6)
where Σ = (η −1 A A + P −1 )−1 , and μ = η −1 ΣA y.
3635
Σ(t) =
η (t)
−1
−1 −1 A A + P (t)
(13)
−1 = P (t) − P (t) A η (t) I + AP (t) A AP (t) , (14) and
−1 μ(t) = P (t) A η (t) I + AP (t) A y.
(15)
Note that the matrix inversion lemma was used to obtain (14) from (13). When t is suf¿ciently large, x(t) can be viewed as a sample generated from the a posteriori pdf f (x|y, a, b). In the tth iteration (t) (t) of the GS-SBL approach, {pn }N can be generated by n=1 and η using the accept/reject algorithm (see, e.g., [10]) whose computational complexity is O(N ). The normal random vector x(t) can be generated by means of an af¿ne transformation of a standard normal random vector z: x(t) = μ(t) + Φ(t) z, where z ∼ N (0, I) and Φ(t) (Φ(t) ) = Σ(t) . Conventionally, Φ(t) is obtained by the Cholesky decomposition. However, the Cholesky decomposition has
a complexity of O(N 3 ), and therefore it is not suitable for largescale compressed sensing applications. Instead, we present an ef¿cient way to generate x(t) as follows. Let z 1 and z 2 be two independent normal random vectors of dimensions M and N , respectively. Then,
−1/2
−1/2 A z 1 + P (t) z2 (16) z 3 = η (t) is
−1 random vector with mean 0 and covariance matrix a normal (see (13)). The complexity of (16) is determined by the Σ(t) complexity of the matrix vector product A z 1 . In other words, Equation (16) has O(L ) complexity. Next, we can verify that x(t) = μ(t) + Σ(t) z 3 ,
(17)
is a normal random vector with mean μ(t) and covariance matrix Σ(t) . Using (14) and (15), we can simplify (17) as follows:
−1 x(t) = P (t) z 3 − P (t) A η (t) I + AP (t) A v, (18) where v = y − AP (t) z 3 . We can use a conjugate gradient (CG) method to compute
−1 η (t) I + AP (t) A v. CG only involves the product of the matrix η (t) I + AP (t) A and a vector. Let the vector be d. Then, (η (t) I + AP (t) A )d = η (t) d + A p(t) (A d) , (19) where is the Hadamard matrix product. Hence, if the matrix vector product has O(L ) complexity, then the computational complexity of CG for a ¿xed number of iterations is O(L ). As a result, each iteration of GS-SBL requires a complexity of O(L ) Àops. Let Tb be the time to begin sampling (also called time of burnin) and Ts be the time of stop sampling. We can estimate xn by the sample median: (Ts ) b) b +1) x ˆn = the median of {x(T , x(T , . . . , xn }. n n
(20)
When there are suf¿ciently many samples (i.e., Ts → ∞), the sample median x ˆn approachesthe a posteriori median (i.e., the u de¿ned u ∞ by −∞ f (xn |y)dxn = u f (xn |y)dxn = 0.5), which is the solution to ˆn |. min Exn |y |xn − x x ˆn
(21)
We conclude this section by making the following comments: 1. We recommend using CoSaMP [5] to obtain the initial estimate of x, because CoSaMP has O(L ) complexity and it converges faster than GS-SBL. CoSaMP requires the knowledge of the sparsity level, i.e., the number of non-zeros elements in x (i.e., K), which is unknown a priori in most applications. Since K is unknown, we simply set it to M/2. 2. Since each iteration of GS-SBL has a O(L ) complexity and the number of iterations is ¿xed to Ts (we set Ts to 500 in the numerical examples), the total computational complexity of GS-SBL is O(L ). 3. GS-SBL needs N (Ts − Tb ) memory to compute (20). To reduce the memory requirement, we can simply store the recently-generated T samples (we set T to 200 in the numerical examples) and we compute the sample mean of these T samples every Tr iterations (we set Tr to 40 in the numerical examples). Then, after Ts iterations, we select the sample mean which minimizes ||y − Ax||2 as the estimate for x.
3636
3. NUMERICAL EXAMPLES In this section, we ¿rst compare the performance of GS-SBL with those of the direct differentiation based SBL (which we simply refer to as SBL), the EM-based SBL (which we refer to as EM-SBL), the Fast Marginal Likelihood Maximization (which we refer to as BCS according to [8]) by using a small- and a large-scale compressed sensing problem. 3.1. Small-Scale Problem First, let us consider a spectral estimation example in time series analysis, which has many applications in diverse ¿elds. The measurement vector y and the spectral component x satisfy the following noisy DCT equations: ym y(tm ) = w(tm )
N n=1
xn cos
π(2n − 1)(tm − 1) + e(tm ), 2N
√ 1/ N , tm = 1; , {tm }M where w(tm ) = m=1 is a random 2/N , tm > 1, subset of the set {1, 2, . . . , N } and the noise term e(tm ) has a normal distribution with zero mean and variance η. In this case, the measurement matrix M is a random sampling matrix and the basis matrix B is a DCT matrix. Similar to the example in [8], we set K = 20, M = 128 and N = 512. The 20 non-zero elements in x are either 1 or −1 and their locations are chosen randomly. Since GS-SBL is initialized by CoSaMP, we also initialize all the algorithms in the numerical examples by CoSaMP. We run SBL and EM-SBL until they satisfy the convergence criterion ˆ (t) ||2 /||ˆ ˆ (t) is the estimate ||ˆ x(t+1) − x x (t) ||2 < 10−3 , where x in the tth iteration. The code for BCS is downloaded from the website “http://people.ee.duke.edu/ lihan/cs/”. BCS also assumes the true noise power η, since otherwise its performance can be rather poor. Note that SBL, EM-SBL and GS-SBL do not assume such a knowledge. The default stopping criteria is used for BCS. We set a = 0.5, b = 10−4 , Tb = 101 and Ts = 500 for GS-SBL. From Figure 1 we can see that GS-SBL gives the most accurate estimate of x. Figure 2(a) shows the average reconstruction errors (from 50 Monte-Carlo trials) as functions of the signal-to-noise ratio (SNR) for K = 20, M = 128 and N = 512. Here, since the signal is either 1 or −1, the SNR is simply de¿ned as the inverse noise power. We can see from Figure 2(a) that GS-SBL has better performance the existing SBL approaches. For SNR≥ 25 dB, GSSBL is 1 − 5 dB better than BCS, 1 − 8 dB better than SBL and 6 − 26 dB better than EM-SBL. Figure 2(b) shows the average reconstruction errors (based on 50 Monte-Carlo trials) as functions of M for SNR= 30 dB, K = 20, and N = 512. Again, GS-SBL has the best performance. For M ≥ 128, GS-SBL is 4−5 dB better than BCS, 5 − 6 dB better than SBL and 4 − 14 dB better than EM-SBL.
3.2. Large-Scale Problem We consider the previous spectral analysis example once again. This time, we set K = 500, M = 3000 and N = 8192. The 500 nonzero elements in x are either 1 or −1 and their locations are chosen randomly. The noise power η is set to be −30 dB. We compare the performance of GS-SBL with that of BCS in this scenario. We make use of the DCT transform to compute the matrix-vector product in the BCS code and in the GS-SBL code. We use the default stopping criteria for BCS, and we set the maximum number of iterations to
105 . As before, we set a = 0.5, b = 10−4 , Tb = 101 and Ts = 500 for GS-SBL. We can see from Figure 3 that the result generated by GS-SBL is more accurate than that of BCS, and its reconstruction error is smaller than that of BCS. The average reconstruction error of GSSBL (from 50 Monte-Carlo Trials) is −24.1 dB and the average reconstruction error of BCS is −19.7 dB. On average, GS-SBL needs 476 seconds and BCS needs 1105 seconds on a server with Intel Xeon 2.33 GHz Quad-Core Processor. Hence, GS-SBL is both faster and more accurate than BCS in this case. For problems of larger sizes, GS-SBL offers further computational savings over BCS.
0.5
0
í0.5
í1
í1.5
100
200
300
[5] D. Needell and J. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples,” Applied and Computational Harmonic Analysis, To appear, Download available at http://arxiv.org/PS cache/arxiv/pdf/0803/0803.2392v2.pdf.
500
í1
100
200
0.5
0
í0.5
í1
200
400
500
(b) 1.5
True value of x BCS
100
300 Index
1
í1.5
300
400
0.5
0
í0.5
í1
í1.5
500
True value of x GSíSBL
1
100
200
Index
300
400
500
Index
(c)
(d)
Fig. 1. The transform coef¿cients estimated by (a) SBL, (b) EMSBL, (c) BCS and (d) GS-SBL for K = 20, M = 128 and N = 512.
0
SBL EMíSBL BCS GSíSBL
Reconstruction error (dB)
í5 í10 í15 í20 í25
0
Reconstruction error (dB)
[2] S. Chen, D. Donoho, and M. Saunders, “Atomic decomposition by basis pursuit,” SIAM Journal on Scienti¿c Computing, vol. 20, pp. 33 – 61, August 1998.
[4] M. Figueiredo, R. Nowak, and S. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems,” IEEE Journal of Selected Topics in Signal Process, vol. 1, pp. 586 – 597, 2007.
0
í0.5
í1.5
Transform coefficient vector x
Transform coefficient vector x
We have proposed a Gibbs sampling based SBL approach for compressed sensing. We refer to this algorithm as GS-SBL. With an ef¿cient sampling step, GS-SBL can have a lower computational complexity than those of the direct differentiation based SBL, the EM-based SBL and the Fast Marginal Likelihood Maximization algorithms. We have also shown by using numerical examples that GS-SBL outperforms these SBL approaches.
[3] I. Daubechies, M. Defrise, and C. Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics, vol. 57, pp. 1413 – 1457, 2004.
0.5
(a)
4. CONCLUSIONS
[1] R. Baraniuk, “Compressive sensing [lecture notes],” IEEE Signal Processing Magazine, vol. 24, pp. 118–121, July 2007.
400
True value of x EMíSBL
1
Index
1.5
5. REFERENCES
1.5
True value of x SBL
1
Transform coefficient vector x
Transform coefficient vector x
1.5
SBL EMíSBL BCS GSíSBL
í5
í10
í15
í20
í30 í35 20
25
30 35 SNR (dB)
40
í25 100
45
200
300 M
400
500
Fig. 2. (a) The average reconstruction errors as functions of SNR for K = 20, M = 128 and N = 512; (b)The average reconstruction errors as functions of M for SNR = 30 dB, K = 20 and N = 512.
[6] M. E. Tipping, “Sparse Bayesian learning and the relevance vector machine,” Journal of Machine Learning Research, vol. 1, pp. 211–244, 2001.
0.5
0
í0.5
í1
í1.5
[9] G. Casella and E. George, “Explaining the gibbs sampler,” The American Statistician, vol. 46, pp. 167–174, Aug. 1992. [10] G. Casella and R. L. Berger, Statistical Inference. Grove, CA: Thomson Learning, Inc., 2001.
Paci¿c
3637
1.5
True value of x BCS
1
Transform coefficient vector x
[8] S. Ji, Y. Xue, and L. Carin, “Bayesian compressive sensing,” IEEE Transactions on Signal Processing, vol. 56, pp. 2346 – 2356, June 2008.
1.5
Transform coefficient vector x
[7] M. Tipping and A. Faul, “Fast marginal likelihood maximisation for sparse Bayesian models,” Proceedings of the Ninth International Workshop on Arti¿cial Intelligence and Statistics, pp. 3–6, Jan. 2003.
100
200
300 Index
(a)
400
500
True value of x GSíSBL
1
0.5
0
í0.5
í1
í1.5
100
200
300
400
500
Index
(b)
Fig. 3. The ¿rst 500 transform coef¿cients estimated by (a) BCS and (b) GS-SBL for SNR= 30 dB, K = 500, M = 3000 and N = 8192.
Peter Stoica
Dept. of Electrical and Computer Eng., Univ. of Florida, Gainesville, FL, USA, email: tanxing@uÀ.edu, [email protected]À.edu
Dept. of Information Technology, Uppsala University, Uppsala, Sweden, Email: [email protected]
ABSTRACT Sparse Bayesian learning (SBL) has been used as a signal recovery algorithm for compressed sensing. It has been shown that SBL is easy to use and can recover sparse signals more accurately than the well-known Basis Pursuit (BP) algorithm. However, the computational complexity of SBL is quite high, which limits its use in large-scale problems. We propose herein an ef¿cient Gibbs sampling approach, referred to as GS-SBL, for compressed sensing. Numerical examples show that GS-SBL can be faster and perform better than the existing SBL approaches.
collecting N data samples and computing N transform coef¿cients in which only the K largest ones are retained and encoded. Fortunately, with the recent advent of compressed sensing techniques (see, e.g., [1]), we can recover the original signal s from much fewer samples. In compressed sensing, the signal is measured through a linear process that can be modeled as projecting the signal s onto a set of vectors {m 1 , m 2 , . . . , m M } to obtain the following samples: ym = mm s +em , m = 1, 2, . . . , M , where em is the measurement noise. If we de¿ne the data vector y = [y1 , y2 , . . . , yM ] , the measurement matrix M = [m 1 , m 2 , . . . , m M ] and the noise vector e = [e1 , e2 , . . . , eM ] , then y can be expressed as y = M s + e. Using the fact that s = Bx, we obtain the following data model:
Index Terms— Compressed Sensing, Sparse Bayesian Learning, Gibbs Sampling
(1)
y = Ax + e, 1. INTRODUCTION We denote vectors by boldface lowercase letters and matrices by boldface uppercase letters. Other notations that we use throughout this paper are given in Table 1.
where A M B is an M × N matrix. To obtain a sparse solution of (1), an l1 -norm based optimization approach, referred to as Basis Pursuit (BP), was proposed in [2]: ˆ = arg min ||x||1 x x
(·) (·)∗ · q μn Σk,l diag{v}
Table 1. Notations used in the text transpose of a vector or a matrix conjugate transpose of a vector or a matrix q-norm. 2-norm when the subscript is omitted. the nth element of the vector μ the (k, l)th element of the matrix Σ A diagonal matrix whose diagonal entries are the elements of the vector v
(2)
where δ is a user parameter. An alternative l1 -norm based optimization is ˆ = arg min ||y − Ax||2 + λ||x||1 . x x
In the past few decades, transform coding has been widely used for data compression. Denote the signal to be compressed as an N -dimensional real-valued vector s. Let B be an orthogonal basis matrix which satis¿es B B = BB = I and de¿ne the transform coef¿cient vector x as x = B s. If x can be approximated ˆ (i.e., x ˆ contains no more than K non-zero by a K-sparse vector x elements and usually K N ), then we only need to encode the ˆ instead of all N elK non-zero elements and their positions in x ements of s. Transform coding is inef¿cient, because it requires ∗ This material is based on research sponsored in part by the of¿ce of Naval Research (ONR) under grants No. N00014-09-1-0211 and No. N00014-07-1-0193, the Army Research Of¿ce (ARO) under Grant No. W911NF-07-1-0450, the National Science Foundation (NSF) under Grant No. ECCS-0729727, the Komen Breast Cancer Foundation under grant No. BCTR0707587, the Swedish Research Council (VR) and the European Research Council (ERC). Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the United States Government.
978-1-4244-4296-6/10/$25.00 ©2010 IEEE
subject to ||y − Ax||2 < δ,
3634
(3)
From a Bayesian point of view, (3) is equivalent to the maximum a posteriori approach for estimating x in (1), assuming that e is a normal random vector and x has a Laplacian a priori distribution. Many ef¿cient methods for solving (3) have been proposed recently, including the iterative shrinkage/thresholding (IST, see, e.g., [3]) and gradient projection for sparse reconstruction (GPSR, e.g., [4]). All these methods depend on the fast computation of the product of A and an arbitrary vector. If the complexity of the matrix vector product is O(L ), then IST and GPSR also have O(L ) complexity. Hence, when A is a fast Fourier Transform (FFT) matrix, a Discrete Cosine Transform (DCT) matrix, a Wavelet Transform matrix or a sparse matrix, then the complexities of IST and GPSR are nearly linear to N . Therefore, they are highly ef¿cient in solving largescale compressed sensing problems. Another ef¿cient sparse signal recovery method is the so-called Compressive Sampling Matching Pursuit (CoSaMP, [5]), which also has O(L ) complexity. One can also obtain a sparse estimate of x via sparse Bayesian learning (SBL), which was originally proposed in the machine learning literature [6]. The model used in SBL is a three-stage hierarchical Bayesian model. The parameters in the Bayesian model were obtained by using the maximum a posteriori (MAP) criterion. However, the MAP problem cannot be solved by setting the derivative
ICASSP 2010
of the cost function to zero, since the derivative is highly nonlinear. Hence, the author proposed two iterative approaches: a direct differentiation approach and an Expectation Maximization (EM) approach (see Appendix A of [6]). The vector x can be estimated by computing its a posteriori mean during the iterative process. A greedy SBL method, referred to as the Fast Marginal Likelihood Maximization algorithm was proposed in [7] and later on applied to compressed sensing in [8] under the name Bayesian Compressive Sensing (BCS). In [8], the authors show that in some cases, the greedy SBL gives a more accurate estimate of the signal than BP algorithm in terms of ˆ ||2 x dB. the reconstruction error de¿ned by err = 10 log10 ||s−B ||s||2 A major advantage of SBL over the l1 -norm based optimization is that the former provides not only a point estimate of x, but also the a posteriori covariance matrix of x, which measures the accuracy of the point estimate. However, both the EM-based SBL and the direct differentiation based SBL have high computational complexities of O(M 2 N ) (see [6]). The Fast Marginal Likelihood Maximization approach has a complexity of O(M N +K(L +K 2 )) (see Appendix A in [7]). This method is indeed fast when K is small. However, when K is large, it can be rather slow. Herein, we propose a novel Gibbs sampling based SBL approach which computes the a posteriori pdf of x in O(L ) Àops. We refer to this approach as GS-SBL. It is known that, in general, the Markov Chain Monte Carlo (MCMC) methods are slower than other Bayesian methods. However, with an ef¿cient sampling step, GSSBL only has a complexity of O(L ). Although we still need to generate a few hundred samples to obtain a good estimate for x, in large-scale compressed sensing problems the total time needed by GS-SBL is still much less than that required by the original SBL proposed in [6]. Furthermore, numerical examples show that GSSBL can outperform the existing SBL approaches in both small- and large-scale compressed sensing problems.
Consider the following Bayesian model:
f (p, η|y, x, a, b) ∝ f (y|x, p, η)f (x|p)f (p|a, b) N −(a+3/2) − 1 (b+ x2n ) 1 (y −Ax2 ) −(M/2) − 2η p 2 e pn e n . (5) ∝η n=1
and f (x|y, p, η, a, b) ∝ f (y|x, p, η)f (x|p)f (p|a, b) Σ−1 (x−μ)
(8)
From (6), x has a conditional normal distribution: x|y, p, η, a, b ∼ N (μ, Σ).
(9)
Using (7), (8) and (9), we can derive the Gibbs sampling algorithm (see, e.g., [9]). Assume that an initial estimate x(0) is given (see the comments on the initial estimator at the end of this section). In the tth iteration, we perform the following two steps: (t)
1. We ¿rst draw a sample of pn from the conditional pdf f (pn |y, x(t−1) , a, b) for n = 1, 2, . . . , N , and draw a sample of η (t) from the conditional pdf f (η|y, x(t−1) , a, b). By using Equation (8), we obtain ⎛
2 ⎞ (t−1) x n 1 ⎟ ⎜ , b+ p(t) ⎠ , n = 1, 2, . . . , N, n ∼ IG ⎝ a + 2 2 (10) and by using Equation (7), we obtain y − Ax(t−1) 2 M η (t) ∼ IG − 1, . 2 2
(11)
(12)
where (4)
where p = [p1 , p2 , . . . , pN ] , P = diag{p}, N (μ, Σ) represents the multivariate normal distribution with mean μ and covariance matrix Σ and IG(a, b) represents the inverse gamma distribution. We assume the inverse gamma prior for pn because we know that p is sparse and the inverse gamma prior promotes sparsity in the estimate of p, when b is small. But for η we do not have any a priori information, so it is proper to assume that f (η) ∝ 1. Based on the Bayesian model (4), we can readily obtain:
1 x2 pn |y, x, a, b ∼ IG a + , b + n . 2 2
x(t) ∼ N (μ(t) , Σ(t) ),
y|x, η ∼ N (Ax, ηI) x|p ∼ N (0, P ), f (η) ∝ 1 pn |a, b ∼ IG(a, b),
1
and
2. Then we draw a sample of x(t) from the conditional pdf of f (x|y, p(t) , η (t) , a, b). De¿ne the matrix P (t) diag{p(t) }. By Equation (9), we obtain
2. AN EFFICIENT GIBBS SAMPLING APPROACH
∝ e− 2 (x−μ)
It follows from (5) that p and η are conditionally independent and that they have the following distributions, respectively, M y − Ax2 η|y, x, a, b ∼ IG − 1, , (7) 2 2
,
(6)
where Σ = (η −1 A A + P −1 )−1 , and μ = η −1 ΣA y.
3635
Σ(t) =
η (t)
−1
−1 −1 A A + P (t)
(13)
−1 = P (t) − P (t) A η (t) I + AP (t) A AP (t) , (14) and
−1 μ(t) = P (t) A η (t) I + AP (t) A y.
(15)
Note that the matrix inversion lemma was used to obtain (14) from (13). When t is suf¿ciently large, x(t) can be viewed as a sample generated from the a posteriori pdf f (x|y, a, b). In the tth iteration (t) (t) of the GS-SBL approach, {pn }N can be generated by n=1 and η using the accept/reject algorithm (see, e.g., [10]) whose computational complexity is O(N ). The normal random vector x(t) can be generated by means of an af¿ne transformation of a standard normal random vector z: x(t) = μ(t) + Φ(t) z, where z ∼ N (0, I) and Φ(t) (Φ(t) ) = Σ(t) . Conventionally, Φ(t) is obtained by the Cholesky decomposition. However, the Cholesky decomposition has
a complexity of O(N 3 ), and therefore it is not suitable for largescale compressed sensing applications. Instead, we present an ef¿cient way to generate x(t) as follows. Let z 1 and z 2 be two independent normal random vectors of dimensions M and N , respectively. Then,
−1/2
−1/2 A z 1 + P (t) z2 (16) z 3 = η (t) is
−1 random vector with mean 0 and covariance matrix a normal (see (13)). The complexity of (16) is determined by the Σ(t) complexity of the matrix vector product A z 1 . In other words, Equation (16) has O(L ) complexity. Next, we can verify that x(t) = μ(t) + Σ(t) z 3 ,
(17)
is a normal random vector with mean μ(t) and covariance matrix Σ(t) . Using (14) and (15), we can simplify (17) as follows:
−1 x(t) = P (t) z 3 − P (t) A η (t) I + AP (t) A v, (18) where v = y − AP (t) z 3 . We can use a conjugate gradient (CG) method to compute
−1 η (t) I + AP (t) A v. CG only involves the product of the matrix η (t) I + AP (t) A and a vector. Let the vector be d. Then, (η (t) I + AP (t) A )d = η (t) d + A p(t) (A d) , (19) where is the Hadamard matrix product. Hence, if the matrix vector product has O(L ) complexity, then the computational complexity of CG for a ¿xed number of iterations is O(L ). As a result, each iteration of GS-SBL requires a complexity of O(L ) Àops. Let Tb be the time to begin sampling (also called time of burnin) and Ts be the time of stop sampling. We can estimate xn by the sample median: (Ts ) b) b +1) x ˆn = the median of {x(T , x(T , . . . , xn }. n n
(20)
When there are suf¿ciently many samples (i.e., Ts → ∞), the sample median x ˆn approachesthe a posteriori median (i.e., the u de¿ned u ∞ by −∞ f (xn |y)dxn = u f (xn |y)dxn = 0.5), which is the solution to ˆn |. min Exn |y |xn − x x ˆn
(21)
We conclude this section by making the following comments: 1. We recommend using CoSaMP [5] to obtain the initial estimate of x, because CoSaMP has O(L ) complexity and it converges faster than GS-SBL. CoSaMP requires the knowledge of the sparsity level, i.e., the number of non-zeros elements in x (i.e., K), which is unknown a priori in most applications. Since K is unknown, we simply set it to M/2. 2. Since each iteration of GS-SBL has a O(L ) complexity and the number of iterations is ¿xed to Ts (we set Ts to 500 in the numerical examples), the total computational complexity of GS-SBL is O(L ). 3. GS-SBL needs N (Ts − Tb ) memory to compute (20). To reduce the memory requirement, we can simply store the recently-generated T samples (we set T to 200 in the numerical examples) and we compute the sample mean of these T samples every Tr iterations (we set Tr to 40 in the numerical examples). Then, after Ts iterations, we select the sample mean which minimizes ||y − Ax||2 as the estimate for x.
3636
3. NUMERICAL EXAMPLES In this section, we ¿rst compare the performance of GS-SBL with those of the direct differentiation based SBL (which we simply refer to as SBL), the EM-based SBL (which we refer to as EM-SBL), the Fast Marginal Likelihood Maximization (which we refer to as BCS according to [8]) by using a small- and a large-scale compressed sensing problem. 3.1. Small-Scale Problem First, let us consider a spectral estimation example in time series analysis, which has many applications in diverse ¿elds. The measurement vector y and the spectral component x satisfy the following noisy DCT equations: ym y(tm ) = w(tm )
N n=1
xn cos
π(2n − 1)(tm − 1) + e(tm ), 2N
√ 1/ N , tm = 1; , {tm }M where w(tm ) = m=1 is a random 2/N , tm > 1, subset of the set {1, 2, . . . , N } and the noise term e(tm ) has a normal distribution with zero mean and variance η. In this case, the measurement matrix M is a random sampling matrix and the basis matrix B is a DCT matrix. Similar to the example in [8], we set K = 20, M = 128 and N = 512. The 20 non-zero elements in x are either 1 or −1 and their locations are chosen randomly. Since GS-SBL is initialized by CoSaMP, we also initialize all the algorithms in the numerical examples by CoSaMP. We run SBL and EM-SBL until they satisfy the convergence criterion ˆ (t) ||2 /||ˆ ˆ (t) is the estimate ||ˆ x(t+1) − x x (t) ||2 < 10−3 , where x in the tth iteration. The code for BCS is downloaded from the website “http://people.ee.duke.edu/ lihan/cs/”. BCS also assumes the true noise power η, since otherwise its performance can be rather poor. Note that SBL, EM-SBL and GS-SBL do not assume such a knowledge. The default stopping criteria is used for BCS. We set a = 0.5, b = 10−4 , Tb = 101 and Ts = 500 for GS-SBL. From Figure 1 we can see that GS-SBL gives the most accurate estimate of x. Figure 2(a) shows the average reconstruction errors (from 50 Monte-Carlo trials) as functions of the signal-to-noise ratio (SNR) for K = 20, M = 128 and N = 512. Here, since the signal is either 1 or −1, the SNR is simply de¿ned as the inverse noise power. We can see from Figure 2(a) that GS-SBL has better performance the existing SBL approaches. For SNR≥ 25 dB, GSSBL is 1 − 5 dB better than BCS, 1 − 8 dB better than SBL and 6 − 26 dB better than EM-SBL. Figure 2(b) shows the average reconstruction errors (based on 50 Monte-Carlo trials) as functions of M for SNR= 30 dB, K = 20, and N = 512. Again, GS-SBL has the best performance. For M ≥ 128, GS-SBL is 4−5 dB better than BCS, 5 − 6 dB better than SBL and 4 − 14 dB better than EM-SBL.
3.2. Large-Scale Problem We consider the previous spectral analysis example once again. This time, we set K = 500, M = 3000 and N = 8192. The 500 nonzero elements in x are either 1 or −1 and their locations are chosen randomly. The noise power η is set to be −30 dB. We compare the performance of GS-SBL with that of BCS in this scenario. We make use of the DCT transform to compute the matrix-vector product in the BCS code and in the GS-SBL code. We use the default stopping criteria for BCS, and we set the maximum number of iterations to
105 . As before, we set a = 0.5, b = 10−4 , Tb = 101 and Ts = 500 for GS-SBL. We can see from Figure 3 that the result generated by GS-SBL is more accurate than that of BCS, and its reconstruction error is smaller than that of BCS. The average reconstruction error of GSSBL (from 50 Monte-Carlo Trials) is −24.1 dB and the average reconstruction error of BCS is −19.7 dB. On average, GS-SBL needs 476 seconds and BCS needs 1105 seconds on a server with Intel Xeon 2.33 GHz Quad-Core Processor. Hence, GS-SBL is both faster and more accurate than BCS in this case. For problems of larger sizes, GS-SBL offers further computational savings over BCS.
0.5
0
í0.5
í1
í1.5
100
200
300
[5] D. Needell and J. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples,” Applied and Computational Harmonic Analysis, To appear, Download available at http://arxiv.org/PS cache/arxiv/pdf/0803/0803.2392v2.pdf.
500
í1
100
200
0.5
0
í0.5
í1
200
400
500
(b) 1.5
True value of x BCS
100
300 Index
1
í1.5
300
400
0.5
0
í0.5
í1
í1.5
500
True value of x GSíSBL
1
100
200
Index
300
400
500
Index
(c)
(d)
Fig. 1. The transform coef¿cients estimated by (a) SBL, (b) EMSBL, (c) BCS and (d) GS-SBL for K = 20, M = 128 and N = 512.
0
SBL EMíSBL BCS GSíSBL
Reconstruction error (dB)
í5 í10 í15 í20 í25
0
Reconstruction error (dB)
[2] S. Chen, D. Donoho, and M. Saunders, “Atomic decomposition by basis pursuit,” SIAM Journal on Scienti¿c Computing, vol. 20, pp. 33 – 61, August 1998.
[4] M. Figueiredo, R. Nowak, and S. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems,” IEEE Journal of Selected Topics in Signal Process, vol. 1, pp. 586 – 597, 2007.
0
í0.5
í1.5
Transform coefficient vector x
Transform coefficient vector x
We have proposed a Gibbs sampling based SBL approach for compressed sensing. We refer to this algorithm as GS-SBL. With an ef¿cient sampling step, GS-SBL can have a lower computational complexity than those of the direct differentiation based SBL, the EM-based SBL and the Fast Marginal Likelihood Maximization algorithms. We have also shown by using numerical examples that GS-SBL outperforms these SBL approaches.
[3] I. Daubechies, M. Defrise, and C. Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics, vol. 57, pp. 1413 – 1457, 2004.
0.5
(a)
4. CONCLUSIONS
[1] R. Baraniuk, “Compressive sensing [lecture notes],” IEEE Signal Processing Magazine, vol. 24, pp. 118–121, July 2007.
400
True value of x EMíSBL
1
Index
1.5
5. REFERENCES
1.5
True value of x SBL
1
Transform coefficient vector x
Transform coefficient vector x
1.5
SBL EMíSBL BCS GSíSBL
í5
í10
í15
í20
í30 í35 20
25
30 35 SNR (dB)
40
í25 100
45
200
300 M
400
500
Fig. 2. (a) The average reconstruction errors as functions of SNR for K = 20, M = 128 and N = 512; (b)The average reconstruction errors as functions of M for SNR = 30 dB, K = 20 and N = 512.
[6] M. E. Tipping, “Sparse Bayesian learning and the relevance vector machine,” Journal of Machine Learning Research, vol. 1, pp. 211–244, 2001.
0.5
0
í0.5
í1
í1.5
[9] G. Casella and E. George, “Explaining the gibbs sampler,” The American Statistician, vol. 46, pp. 167–174, Aug. 1992. [10] G. Casella and R. L. Berger, Statistical Inference. Grove, CA: Thomson Learning, Inc., 2001.
Paci¿c
3637
1.5
True value of x BCS
1
Transform coefficient vector x
[8] S. Ji, Y. Xue, and L. Carin, “Bayesian compressive sensing,” IEEE Transactions on Signal Processing, vol. 56, pp. 2346 – 2356, June 2008.
1.5
Transform coefficient vector x
[7] M. Tipping and A. Faul, “Fast marginal likelihood maximisation for sparse Bayesian models,” Proceedings of the Ninth International Workshop on Arti¿cial Intelligence and Statistics, pp. 3–6, Jan. 2003.
100
200
300 Index
(a)
400
500
True value of x GSíSBL
1
0.5
0
í0.5
í1
í1.5
100
200
300
400
500
Index
(b)
Fig. 3. The ¿rst 500 transform coef¿cients estimated by (a) BCS and (b) GS-SBL for SNR= 30 dB, K = 500, M = 3000 and N = 8192.
