Bayesian Compressive Sensing Using Normal Product Priors
Zhou Zhou Kaihui Liu Jun Fang














Bayesian Compressive Sensing Using Normal Product Priors
    Zhou Zhou, Kaihui Liu, and Jun Fang,
Member, IEEE
Zhou Zhou, Kaihui Liu, and Jun Fang are with the National Key Laboratory on Communications,
University of Electronic Science and Technology of China, Chengdu
611731, China, Email: ZhouZhou@std.uestc.edu.cn; kaihuil@std.uestc.edu.cn; JunFang@uestc.edu.cn


This work was supported in part by the National
Science Foundation of China under Grant 61172114.
    Received: date / Accepted: date
=======================================================================================================================================================================================================================================================================================================================================================================================================================






In this paper, we introduce a new sparsity-promoting prior,
namely, the “normal product” prior, and develop an efficient
algorithm for sparse signal recovery under the Bayesian framework.
The normal product distribution is the distribution of a product
of two normally distributed variables with zero means and possibly
different variances. Like other sparsity-encouraging distributions
such as the Student's t-distribution, the normal product
distribution has a sharp peak at origin, which makes it a suitable
prior to encourage sparse solutions. A two-stage normal
product-based hierarchical model is proposed. We resort to the
variational Bayesian (VB) method to perform the inference.
Simulations are conducted to illustrate the effectiveness of our
proposed algorithm as compared with other state-of-the-art
compressed sensing algorithms.


Compressed Sensing, sparse Bayesian learning, normal product
prior.





















§ INTRODUCTION


Compressed sensing
<cit.> is a new data
acquisition technique that has attracted much attention over the
past decade. Existing methods for compressed sensing can generally
be classified into the following categories, i.e. the greedy
pursuit approach <cit.>, the convex
relaxation-type approach <cit.> and the nonconvex
optimization method <cit.>. Another class of
compressed sensing techniques that have received increasing
attention are Bayesian methods. In the Bayesian framework, the
signal is usually assigned a sparsity-encouraging prior, such as
the Laplace prior and the Gaussian-inverse Gamma prior
<cit.>, to encourage sparse solutions. It has
been shown in a series of experiments that Bayesian compressed
sensing techniques <cit.> demonstrate superiority
for sparse signal recovery as compared with the greedy methods and
the basis pursuit method. One of the most popular prior model for
Bayesian compressed sensing is the Gaussian-inverse Gamma prior
proposed in <cit.>. The Gaussian-inverse Gamma
prior is a two-layer hierarchical model in which the first layer
specifies a Gaussian prior to the sparse signal and an inverse
Gamma priori is assigned to the parameters characterizing the
Gaussian prior. As discussed in <cit.>, this
two-stage hierarchical model is equivalent to imposing a Student's
t-distribution on the sparse signal. Besides the
Gaussian-inverse Gamma prior, authors in
<cit.> employ Laplace priors for the sparse
signal to promote sparse solutions.

In this paper, we introduce a new sparsity-encouraging prior,
namely, the normal product (NP) prior, for sparse signal recovery.
A two-stage normal product-based hierarchical model is
established, and we resort to the variational Bayesian (VB) method
to perform the inference for the hierarchical model. Our
experiments show that the proposed algorithm achieves similar
performance as the sparse Bayesian learning method, while with
less computational complexity.







§ THE BAYESIAN NETWORK



In the context of compressed sensing, we are given a noise corrupted linear measurements of a vector x_0

    y = Ax_0 + n

Here A∈R^M × N is an underdetermined
measurement matrix with low coherence of columns, n∈R^M represents the acquisition noise vector and x_0 a sparse signal. For the inverse problem, both x_0
and n are unknown, and our goal is to recover x_0
from y. We formulate the observation noise as zero mean
independent Gaussian variable with variance ξ ^2 i.e.
y∼ N(A x_0,ξ ^2I) and
seek a sparse signal xΔ =
(x_1,x_2, …x_N) for x_0.In this section,
we utilize a two-stage hierarchical bayesian prior for the signal
model. In the first layer of signal model, we use the Normal
Product (NP) distribution as sparseness prior:

    x∼ NP(0,Σ)

where 0∈R^N stands for its mean and Σ given as a diagonal matrix diag(σ _1^2,σ _2^2, … ,σ _N^2) ∈R^N × N represents its variance. For each element x_i of x , its probability density distribution is K_0( |x_i|/σ _i)/
 .
 -πσ _i, which exhibits a sharp peak at the origin and heavy tails, where K_0( · ) is the zero order modified Bessel function of the second kind<cit.> thus the probability density function (PDF) of x is
 
    1/π ^N∏_i = 1^N 1/σ _iK_0(| x_i|/σ _i)
    .

  A NP distributed scaler variate x_i ∈R can be decomposed as the product of two independent normally distributed variables i.e. if x_i ∼ NP(0,σ_i^2 ), then there are a_i ∼ N(0,κ_i^2 ) and b_i ∼ N(0,γ_i^2 )
 satisfying x_i = a_ib_i with moment relationship σ_i  = γ_i κ_i<cit.>. We call this property as the generating rule of Normal Product distribution. Similarly for the vector x,
 we can decompose it into the Hadamard product of two independent virtual normally distributed vector variables a and b whose
variance are diagonal matrix κ^2Δ =diag(κ _1^2,κ _2^2, … ,κ _N^2) and γ^2Δ = diag(γ _1^2,γ _2^2, … ,γ _N^2) respectively i.e. x  = a ∘b and Σ = κ^2∘γ^2, where ∘ denote the Hadamard product.
Finally, we set σ _i^ - 2 as a realization of Gamma hyperprior Γ (α ,β ) and choose α→ 0,β→ 0 to construct a sharper distribution of x during the variational procedure <cit.> for the second layer of signal.According to the generating rule of Normal Product mentioned above, we can add two parallel nodes a and b to construct one latent
 layer. Then the posterior distribution can be expressed as:

    P_x(x,Σ|y;α ,β ) = ∫P(x,a,b,Σ|y;α ,β )dadb

where

    [                               P(x,a,b,Σ|y;α ,β );       ∝ P(y|x)P(x|a,b)P(a|κ^2)P(b|γ^2)P(Σ|α ,β ); = P(y - Ax)δ (x - a∘b)P(a|κ^2)P(b|γ^2)P(Σ|α ,β ) ]


and δ(·) denotes the Dirac Delta function. In this expression we know that the value of x must
be equal to a∘b while keeping the value of P(x,a,b,Σ|y;α ,β ) nonzero.
Consider the Bayesian Risk function:

    [ E(L(x̂,x)) = ∫L(x̂,x)P_x(x,Σ|y;α ,β )dxdΣ;     ∝∫L(x̂,a∘b)P_a,b(a,b,Σ|y;α ,β )dadbdΣ ]

where

    [                      P_a,b(a,b,Σ|y;α ,β ); = P(y - A(a∘b))P(a|κ^2)P(b|γ^2)P(Σ|α ,β ) ]

and L(x̂,x) represents the loss function.
Thus, we can replace P_x(x,Σ|y;α ,β ) with P_a,b(a,b,Σ|y;α ,β )
 in the inference procedure while maintaining the value of the Bayesian risk function at a consistent level. So the modified Bayesian Network is:


    y∼ N(A(a∘b),ξ ^2I)


    a∼ N(0,κ^2)


    b∼ N(0,γ^2)


    σ _i^ - 2 = γ _i^ - 2κ _i^ - 2∼Γ (α ,β )

as depicted in Fig. <ref>.




§ THE VARATIONAL BAYESIAN INFERENCE

In order to infer the bayesian network, we use the mean-field
variational Bayesian method to analytically approximate the entire
posterior distribution of latent variable. In this theorem, it was
assumed that the approximate variational distribution
q(z|y) where z stands for the latent
variable in the model could be factorized into the product
∏_i = 1^l q_i(z_i|y) for the
partition z_1,…, z_l of z. It
could be shown that the solution of the variational method for
each factor q_i(z_i|y) could be written as:

    q_i^*(z_i|y) = e^ < ln p(z|y)) > _z_i/∫e^ < ln p(z|y)) > _z_idz_i

where < ln p(z|y)) > _z_i means taking the expectation over the set of z using the variational distribution except z_i<cit.> .
    Applying the variational Bayesian method upon the bayesian model mentioned in secetion <ref>, the posterior approximation of Σ,a,b are respectively:

    ln (q^*(a)) ∝  < ln P(a,b,Σ|y;α,β) > _q^*(b),q^*(Σ)


    ln (q^*(b)) ∝  < ln P(a,b,Σ|y;α,β) > _q^*(a),q^*(Σ)


    ln (q^*(Σ)) ∝  < ln P(a,b,Σ|y;α,β) > _q^*(b),q^*(a)




 §.§ posterior approximation of a and b

Substituting (<ref>), (<ref>), (<ref>) and (<ref>) into (<ref>), we can arrive at:

    q^*(b)∼ N(Γ'^ - 1c',Γ'^ - 1)

where

    Γ' = 1/ξ ^2(A^TA) ∘  < aa^T > _q^*(a) +  < γ^ - 2 > _q^*(γ)
 and

    c' = 1/ξ ^2diag(  < a > _q^*(a))A^Ty
.
From the principle of variational inference, we know that q( b) is an approximation of P(b|y) i.e. P(b|y) ≈ q( b) and in order to get a more concise iteration formula, we choose the relaxation of the secondary moment as:


    < aa^T > _q^*(a)≈  < a > _q^**(a) < a^T > _q^**(a)


    < γ^ - 2 > _q^*(γ)≈  < γ > _q^**(γ)^ - 2


This relaxation can be interpreted as ignoring the posteriori variance during the updating, which is inspired by the fact that after the learning procedure
 the posteriori variance always approaches zero to ensure the posteriori mean's concentration on the estimated value.

Then the corrected posterior approximation is q^**(b)∼ N(Γ^ - 1c,Γ^ - 1) where

    [ Γ^ - 1 = {1/ξ ^2( Adiag(  < a > _q^**(a)))^TAdiag(  < a > _q^**(a));                                       +  < γ > _q^**(γ)^ - 2} ^ - 1 ]

and

    c = 1/ξ ^2diag(  < a > _q^**(a))A^Ty

In the noise free case, using Woodbury identity we have:

    [                                               lim_ξ→ 0Γ^ - 1c; =  < γ > _q^**(γ)(Adiag( < a > _q^**(a)) < γ > _q^**(γ))^ + y ]

and

    [                                                      lim _ξ→ 0Γ^ - 1; = (I -  < γ > _q^**(γ)( Adiag(  < b > _q^**(b)) < γ > _q^**(γ))^ + );                            ×Adiag(  < a > _q^**(a)) < γ > _q^**(γ)^2 ]

Similarly, for a, substituting (<ref>), (<ref>), (<ref>) and (<ref>) into (<ref>), we have:

    q^*(a)∼ N( H'^ - 1f', H'^ - 1)

where

    H' = 1/ξ ^2(A^TA) ∘  < bb^T > _q^*(b) +  < κ^ - 2 > _q^*(κ)

and

    f' = 1/ξ ^2diag(  < b > _q^*(b))A^Ty

So the corrected approximation is:

    q^**(a)∼ N( H^ - 1f, H^ - 1)

Then again, in the noiseless case, we have

    [                                             lim _ξ→ 0 H^ - 1f; =  < κ > _q^**(κ)(Adiag( < b > _q^**(b)) < κ > _q^**(κ))^ + y ]

and

    [                                                      lim _ξ→ 0H^ - 1; = (I -  < κ > _q^**(κ)( Adiag(  < b > _q^**(b)) < κ > _q^**(κ))^ + );                            ×Adiag(  < b > _q^**(b)) < κ > _q^**(κ)^2 ]






 §.§ posterior approximation of Σ 


As discussed in section <ref>, we have:

    [ ln (q^*(Σ)) = ln (q^*(κ^2∘γ^2)) = ln (q^*(κ^2)) + ln (q^*(γ^2));                 ∝  < ln P(a,b,κ^2∘γ^2|y;α,β) > _q^**(b),q^**(a) ]

Substituting (<ref>), (<ref>), (<ref>) into (<ref>) and using the separability of ln (q^*(κ^2)), it can be shown that:

    [ lnq( κ _i^2) ∝  - 1/2(  < a_i^2 > _q^**(a) + 2β  < γ _i^ - 2 > _q^*(γ _i^2))κ _i^ - 2;                                                        + (α  + 1/2 - 1)ln( κ _i^ - 2) ]

which means

    κ _i^ - 2∼Γ (1/2 + α ,1/2( < a_i^2 > _q^**(a) + 2β  < γ _i^ - 2 > _q^*(γ^2)))

and

    < κ _i^ - 2 > _q^*(κ ^2) = 1/2 + α/1/2( < a_i^2 > _q^**(a) + 2β  < γ _i^ - 2 > _q^*(γ^2))


Similarly, we have:

    γ _i^ - 2∼Γ (1/2 + α ,1/2( < b_i^2 > _q^**(b) + 2β  < κ _i^ - 2 > _q^*(κ^2)))

and

    < γ _i^ - 2 > _q^*(γ^2 ) = 1/2 + α/1/2( < b_i^2 > _q^**(b) + 2β  < κ _i^ - 2 > _q^*(κ^2))





 §.§ The Proposed Algorithm 

As a result, we summarize the procedure above as two algorithms
named “NP-0” and “NP-1” which represent the inference results
for the one layer signal model and two layer signal model
respectively. The difference between Np-0 and NP-1 is whether the
learning process updates the precision of Normal Product.


The algorithm NP-0 is similar to the FOCUSS algorithm<cit.>. The difference between them is that NP-0 uses < a>^(t) and < b>^(t), the decomposed component of x^(t) in an interleaved way to update x^(t+1) while FOCUSS uses the whole x^(t) to directly update x^(t+1). Furthermore, we set the initial value of the algorithm as a constant to avoid the local minimum results being returned.





The algorithm NP-1 looks like coupling two sparse Bayesian learning (SBL) <cit.> procedures together and we will experientially
prove that the MSE via NP-1 descends faster than SBL in section <ref>.



§ NUMERICAL RESULTS

In this section, we compare our proposed algorithms among SBL, iterative reweighted least squares (IRLS)<cit.> and basis pursuit (BP)<cit.>.
In the following results, we set the dimension of original signal to be 100 and in each experiment every entries of the sensing matrix A is uniformly distributed as standard normal random variable. The stop criterion parameters ε and t_max in both algorithm NP-0 and NP-1 are set as 10^-3 and 300 respectively.The success rate in the Fig. <ref> is defined by the average probability of successful restoration, which is counted as the ratio between the number of successful trials and the total number of the independent experiments. It collects the result of 300 Monte Carlo trials and
we assume successful trial if x̂ - x_2/
 .
 -x_2 < 10^ - 3, where x̂ is the estimate of the original signal x.
The mean squared error (MSE) in Fig. <ref> which measures the average of the squares of the difference between the estimator and estimated signal is calculated as follow:

    MSE(dB) = 20log (x̂ - x_2/
     .
     -x_2)
Fig. <ref> demonstrats the superior performance of the NP-1 algorithm in a few measurements. It shows the success rate of respective algorithms versus the number of measurement M and sparse level K of the signal. We can see that the reconstruction performance of the one layer normal product algorithm(NP-0) is almost the same as BP, which is the MAP estimation of unknowns in the one layer Laplace prior framework. The comparisons also show that the two layers normal product algorithm(NP-1) requires as few number of CS measurements as SBL while inheriting the similar reconstruction precision.

In Fig. <ref>, the comparison among the computational cost of the three algorithms on the reconstruction performance are performed on a personal PC with dual-core 3.1 GHz CPU and 4GB RAM. It is interesting to note that the MSE of NP-1 descends faster than SBL in Fig. 4(a) while it has the same sparsity-undersampling tradeoff performance as SBL. It should also be noticed that we haven't use any basis pruning technology as in <cit.> to reduce the computational complexity. To understand the reason of low computational time from Fig. 4(b), we observe that the NP-1 requires only a few numbers of iterations to arrive at stop condition. In summary, these experimental results confirm that the NP-1 is a fast Sparse Bayesian algorithm.















§ CONCULSION


In this letter, we formulated a Normal Product prior based Bayesian framework to solve the compressed sensing problem in noise free case. Using this framework, we derived two algorithms named NP-0 and NP-1 and compared them with different algorithms. We have shown that our algorithm NP-1 has the similar reconstruction performance with SBL while the interleaved updating procedure provides improved performance in computational times.

























IEEEtran














