Compressive Sensing via Variational Bayesian Inference under Two Widely Used Priors: Modeling, Comparison and Discussion

Abstract

Compressive sensing is a sub-Nyquist sampling technique for efficient signal acquisition and reconstruction of sparse or compressible signals. In order to account for the sparsity of the underlying signal of interest, it is common to use sparsifying priors such as Bernoulli–Gaussian-inverse Gamma (BGiG) and Gaussian-inverse Gamma (GiG) priors on the components of the signal. With the introduction of variational Bayesian inference, the sparse Bayesian learning (SBL) methods for solving the inverse problem of compressive sensing have received significant interest as the SBL methods become more efficient in terms of execution time. In this paper, we consider the sparse signal recovery problem using compressive sensing and the variational Bayesian (VB) inference framework. More specifically, we consider two widely used Bayesian models of BGiG and GiG for modeling the underlying sparse signal for this problem. Although these two models have been widely used for sparse recovery problems under various signal structures, the question of which model can outperform the other for sparse signal recovery under no specific structure has yet to be fully addressed under the VB inference setting. Here, we study these two models specifically under VB inference in detail, provide some motivating examples regarding the issues in signal reconstruction that may occur under each model, perform comparisons and provide suggestions on how to improve the performance of each model.

1. Introduction

Compressive sensing (CS) involves efficient signal acquisition and reconstruction techniques in a sub-Nyquist sampling sense. The CS framework can capture the vital information of the underlying signal via a small number of measurements while retaining the ability to reconstruct the signal. CS operates under the assumption that the signal is compressible or sparse, and the number and location of dominating components are unknown in most cases [1,2,3]. Compressibility or sparsity means that the signal has few dominating elements under some proper basis. CS has been used in a variety of applications such as the single-pixel camera, missing pixels and inpainting removal of images, biomedical such as heart rate estimation, internet of things (IoT), geostatistical data analysis, seismic tomography, communications such as blind multi-narrowband signals sampling and recovery, the direction of arrival (DoA) estimation, spectrum sharing of radar and communication signals, wireless networks and many more [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. In the linear CS framework, the problem is posed as

$$\mathbf{y}=A{\mathbf{x}}_s+\mathbf{e},$$

(1)

where $\mathbf{y}\in {\mathbb{R}}^M$ contains the measurements, ${\mathbf{x}}_s\in {\mathbb{R}}^N$ is the sparse signal of interest, $\mathbf{e}$ is the noise representing either the measurement noise or the insignificant coefficients of ${\mathbf{x}}_s$ and, generally, $M\ll N$ [1,2]. The measurement matrix can be defined as $A=\Phi \Psi$, where $\Phi$ is the sensing design matrix and $\Psi$ is a proper sparsifying basis. There exist various approaches to solve for ${\mathbf{x}}_s$ in (1) including greedy-based, convex-based, thresholding-based and sparse Bayesian learning (SBL) algorithms [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64]. Typically, the performance of CS reconstruction is determined in terms of the mean-squared reconstruction error. In this paper, we are also interested in the more demanding requirements of the probability of detection and the false alarm of the nonzero components. This is of more interest to CS applications such as blind multinarrowband signals, spectrum sharing RADAR, etc. [11,12,13,14,15].

The focus of this paper is on sparse Bayesian learning (SBL) for the CS problem. Bayesian learning models are flexible in incorporating prior knowledge of the characteristics of the underlying signal into the model. Bayesian learning also provides a distribution of the hidden variables, which is more informative than the point estimate approaches. A prior favoring the sparsity or compressibility in ${\mathbf{x}}_s$ can be represented in the SBL framework via Gaussian-inverse Gamma (GiG), Laplace-inverse Gamma (LiG), Bernoulli–Gaussian-inverse Gamma (BGiG), often referred to as spike-and-slab prior, etc. [27,46,47,48,49,50,51,52,53,54,55,56,57,58,59]. The inference on parameters and hidden variables in these models is usually made using Markov chain Monte Carlo (MCMC) and variational Bayes (VB) [27,45,46,47,48,49,50,51,52]. In this paper, we focus on the two most commonly used SBL prior models for solving the inverse problem of compressive sensing: Bernoulli–Gaussian-inverse Gamma (BGiG) prior and Gaussian-inverse Gamma (GiG). These models have been widely used, along with some additional priors, for sparse recovery of signals or images with block-sparsity/clustering patterns, sparse signals with correlated coefficients or other structured patterns [26,27,48,49,50,51,62,63].

We use VB inference to estimate the variables and parameters of the SBL model. VB is preferred over MCMC because MCMC is computationally expensive, though it can numerically approximate exact posterior distributions with a sufficient amount of computation. The convergence diagnostic of MCMC requires additional work, such as measuring the potential scale reduction factor (PSRF) for all the hidden variables and parameters of the model or monitoring their trace plots [45,50,51,52,65]. In contrast, VB inference can lead to a reasonable approximation of the exact posteriors, using less computation than MCMC and less effort to monitor the convergence [45,51,66,67]. In this paper, we present the derivation of the update rules of the parameters and variables using VB inference for both the BGiG and GiG models. (Portions of this derivation have been previously presented in [68,69]). Although these prior models have been widely used in various applications of compressive sensing, the study of the overall performance of these models under VB inference has yet to be thoroughly investigated. The preference for one model over the other becomes crucial when dealing with moderate or low sampling ratios, which we discuss in this paper. Here, we study the issues associated with each model via some motivational examples. Pre-/postprocessing approaches will then be discussed to tackle the issues. Finally, the overall performance of BGiG and GiG is compared.

The remainder of this work is organized as follows. In Section 2, we present a brief background on VB inference. We study Bernoulli–Gaussian-inverse Gamma modeling for CS using VB in Section 3. Some motivational examples are provided to show the issues with this approach. Section 4 represents Gaussian-inverse Gamma modeling, the associated update rules using VB inference and a motivational example of the issue that may occur using this approach. In Section 5, we study the improvement of the performances of the models after some pre-/postprocessing along with simulation results and comparisons. Section 6 concludes this work.

2. Variational Bayesian Inference

Variational Bayes (VB) is an effective approach to approximate intractable integrals that may arise in Bayesian inference. The main idea behind variational methods is to use a family of distributions over the latent variables with their own variational parameters. VB is a fast alternative to sampling methods such as Markov chain Monte Carlo (MCMC) and Sequential Monte Carlo (SMC) for performing approximate Bayesian inference [70,71]. For a probabilistic model with unknown parameters $\theta$ and hidden variables $\mathbf{x}$, the posterior distribution of the unknowns, given a set of observations $\mathbf{y}$, can be written as $p(\mathbf{x},\theta |\mathbf{y})=p(\mathbf{x},\theta ,\mathbf{y})/p(\mathbf{y})$. Finding the exact posterior in closed form to perform the inference would be a challenge, as the marginal distribution $p\left(\mathbf{y}\right)=\int p(\mathbf{y},\mathbf{x},\theta )d\mathbf{x}d\theta$ is often intractable. As an efficient approximation method for such inference problems, VB provides an analytical approximation to the posterior $p(\mathbf{x},\theta |\mathbf{y})$. VB approximates the joint density $p(\mathbf{x},\theta |\mathbf{y})$ via a variational distribution $Q_{x,\theta }(\mathbf{x},\theta )$, i.e., $p(\mathbf{x},\theta |\mathbf{y})\approx Q_{x,\theta }(\mathbf{x},\theta )$. VB assumes that the distribution Q can be fully factorized with respect to the unknown parameters and hidden variables, i.e.,

$$\begin{array}{cc}\hfill Q_{x,\theta }(\mathbf{x},\theta )& =q_x\left(\mathbf{x}\right)q_{\theta }\left(\theta \right)\hfill \\ & =\prod _{i=1}^I{q}_x\left(x_i\right)\prod _{j=1}^J{q}_{\theta }\left({\theta }_j\right),\hfill \end{array}$$

where I and J are the number of unknown parameters and hidden variables, respectively. This independence assumption in VB further simplifies the search for a closed-form solution to the approximation of the actual posterior. We desire to select the variational distribution $Q_{x,\theta }^{\ast }(\mathbf{x},\theta )$ as close as possible to $p(\mathbf{x},\theta |\mathbf{y})$, where the closeness metric for distribution $Q_{x,\theta }(\mathbf{x},\theta )$ is formulated as minimizing the Kullback–Leibler (KL) divergence of the approximation $Q_{x,\theta }(\mathbf{x},\theta )$ and the true posterior $p(\mathbf{x},\theta |\mathbf{y})$ as

$$\begin{array}{cc}\hfill Q_{x,\theta }^{🟉}(\mathbf{x},\theta )& =\underset{Q_{x,\theta }(\mathbf{x},\theta )}{\arg \min }\mathrm{KL}\left(Q_{x,\theta }(\mathbf{x},\theta )\left|\right|p(\mathbf{x},\theta |\mathbf{y})\right)\hfill \\ & =\underset{Q_{x,\theta }(\mathbf{x},\theta )}{\arg \min }\int Q_{x,\theta }(\mathbf{x},\theta )log\frac{Q_{x,\theta }(\mathbf{x},\theta )}{p(\mathbf{x},\theta |\mathbf{y})}d\mathbf{x}d\theta .\hfill \end{array}$$

The quantity $logp\left(\mathbf{y}\right)$ can be written as $logp\left(\mathbf{y}\right)=log\{\int p(\mathbf{x},\theta ,\mathbf{y}\left)d\mathbf{x}d\theta \right\}$. Then, defining

$$F\left(Q_{x,\theta }(\mathbf{x},\theta )\right)=\int Q_{x,\theta }(\mathbf{x},\theta )log\frac{p(\mathbf{x},\theta ,\mathbf{y})}{Q_{x,\theta }(\mathbf{x},\theta )}d\mathbf{x}d\theta .$$

It is straightforward to show that

$$logp\left(\mathbf{y}\right)=F\left(Q_{x,\theta }(\mathbf{x},\theta )\right)+\mathrm{KL}\left(Q_{x,\theta }(\mathbf{x},\theta ),p(\mathbf{x},\theta |\mathbf{y})\right).$$

Since (by Jensen’s inequality) $KL(Q_{x,\theta }(\mathbf{x},\theta )\ge 0$, $log(p\left(\mathbf{y}\right)\ge F(Q_{x,\theta }(\mathbf{x},\theta )$. Since $log\left(p\right(\mathbf{y})$ is constant with respect to $Q_{x,\theta }$, minimizing the KL-divergence between the actual posterior distribution and the variational distribution is equivalent to maximizing the lower bound $F(·)$ [66,67]. Since the term $p\left(\mathbf{y}\right)$ in $p(\mathbf{x},\theta |\mathbf{y})=p(\mathbf{x},\theta ,\mathbf{y}\left)p\right(\mathbf{y})$ does not involve the variational distribution $Q_{x,\theta }(\mathbf{x},\theta )$, this term can be ignored when maximizing $F(·)$. The lower bound $F(·)$ on the model log-marginal likelihood can be iteratively optimized until the convergence by the following update rules [66,72].

VB-E step: 

$$q_x^{[t+1]}\left(\mathbf{x}\right)\propto exp\left\{{\mathrm{E}}_{q_{\theta }^{\left[t\right]}}[logp(\mathbf{x},\mathbf{y}|\theta )]\right\}$$

(2)

VB-M step: 

$$q_{\theta }^{[t+1]}\left(\theta \right)\propto p\left(\theta \right)exp\left\{{\mathrm{E}}_{q_x^{[t+1]}}[logp(\mathbf{x},\mathbf{y}|\theta )]\right\}$$

(3)

This results in an iterative algorithm analogous to the expectation-maximization (EM) approach.

3. Bernoulli–Gaussian-Inverse Gamma Modeling and SBL(BGiG) Algorithm

In the inverse problem of CS defined in (1), the goal is to recover the sparse vector ${\mathbf{x}}_s$. In the Bernoulli–Gaussian-inverse Gamma model, the sparse solution is defined as

$${\mathbf{x}}_s=(\mathbf{s}\circ \mathbf{x}),$$

(4)

where $\mathbf{s}$ is a binary support vector indicating the non-zero locations in the solution, $\mathbf{x}$ represents values of the solution and ∘ is Hadamard (element-by-element) product [47]. We refer to the algorithm associated with this Bayesian modeling based on VB inference as SBL(BGiG). SBL using VB inference for the clustered pattern of sparse signals has already been investigated in the recent literature [45,50,51,58]. In this paper, however, we intend to focus on the ordinary SBL using VB inference modeling without promoting any structure on the supports other than sparsity itself. We show that when the sampling ratio is moderate or low (with respect to the sparsity level), the reconstruction performance becomes sensitive to selecting the support-related hyperparameters.

We define a set of priors as follows [47,68,69]. We model the elements of vector $\mathbf{s}$ as

$$s_n\sim \mathrm{Bernoulli}\left({\gamma }_n\right),{\gamma }_n\sim \mathrm{Beta}({\alpha }_0,{\beta }_0),\forall n,$$

(5)

where ${\alpha }_0$ and ${\beta }_0$ are the support-related hyperparameters. Setting ${\alpha }_0$ and ${\beta }_0$ to small values and with ${\alpha }_0\ll {\beta }_0$ encourages $\mathbf{s}$ to be sparse on average. The prior on the solution value vector is defined as

$$\mathbf{x}\sim \mathcal{N}(0,{\tau }^{-1}{I}_N),\tau \sim \mathrm{Gamma}(a_0,b_0).$$

(6)

Here, $\tau$ is the precision value. Finally, the prior on the noise is

$$\mathbf{e}\sim \mathcal{N}(0,{ϵ}^{-1}{I}_M),ϵ\sim \mathrm{Gamma}({\theta }_0,{\theta }_1),$$

(7)

where ${\theta }_0$ and ${\theta }_1$ are set to small positive values.

3.1. Update Rules of SBL(BGiG) Using VB Inference

According to the VB algorithm defined in (2) and (3), the update rule of the variables and parameters of the BGiG model can be simplified as follows [68]. The details of these derivations appear in Appendix A.1.

• Update rule for the support vector $\mathbf{s}$

$$q\left(s_n\right|-)\sim \mathrm{Bernoulli}\left(\frac{1}{1+c_n{\kappa }_n}\right),\forall n=1,\dots ,N,$$

where conditioning on − denotes conditioning on all relevant variables and observations. Therefore,

$${\tilde{s}}_n=\frac{1}{1+c_n{\kappa }_n},\forall n=1,\dots ,N,$$

(8)

where

$$\begin{array}{cc}\hfill c_n& :=e^{\psi \left({\beta }_{1,n}\right)-\psi \left({\alpha }_{1,n}\right)},\hfill \\ \hfill {\kappa }_n& :=e^{\frac{1}{2}\widetilde{ϵ}\left(\parallel {\mathbf{a}}_n{\parallel }_2^2({\widetilde{x_n}}^2+{\sigma }_{\widetilde{x_n}}^2)-2\widetilde{x_n}{\mathbf{a}}_n^T{\tilde{\mathbf{y}}}^{-n}\right)},\hfill \\ \hfill {\tilde{y}}_m^{-n}& :=y_m-\sum _{l\ne n}^N{a}_{ml}{\tilde{s}}_l{\tilde{x}}_l.\hfill \end{array}$$

(9)

Here, $\tilde{\mathbf{x}}:=<\mathbf{x}{>}_{q_x}$, $\psi$ is the digamma function (the logarithmic derivative of the gamma function), and ${\tilde{\mathbf{y}}}^{-n}={[{\tilde{y}}_1^{-n},\dots ,{\tilde{y}}_M^{-n}]}^T$.

• Update rule for the solution value matrix $\mathbf{x}$

$$q\left(\mathbf{x}\right|-)\sim \mathcal{N}(\tilde{\mathbf{x}},{{\Sigma }}_{\tilde{x}}),$$

where

$${{\Sigma }}_{\tilde{x}}={\left(\tilde{\tau }{I}_N+\widetilde{ϵ}\tilde{\Phi }\right)}^{-1}\mathrm{and}\tilde{\mathbf{x}}=\widetilde{ϵ}{{\Sigma }}_{\tilde{x}}\mathrm{diag}\left(\tilde{\mathbf{s}}\right)A^T\mathbf{y},$$

(10)

and where $\mathrm{diag}\left(\mathbf{s}\right)$ denotes a diagonal matrix with the components of $\mathbf{s}$ on its main diagonal, and

$$\tilde{\Phi }:=\left[\left(A^{T}A\right)\circ \left(\tilde{\mathbf{s}}{\tilde{\mathbf{s}}}^T+\mathrm{diag}(\tilde{\mathbf{s}}\circ (1-\tilde{\mathbf{s}}))\right)\right].$$

(11)

• Update rule for ${\gamma }_n$

$$q\left({\gamma }_n\right|-)\sim \mathrm{Beta}({\alpha }_{1,n},{\beta }_{1,n}),\forall n=1,\dots ,N.$$

Therefore,

$$\widetilde{{\gamma }_n}=\frac{{\alpha }_{1,n}}{{\alpha }_{1,n}+{\beta }_{1,n}},\forall n=1,\dots ,N,$$

(12)

where

$${\alpha }_{1,n}:={\alpha }_0+{\tilde{s}}_n\mathrm{and}{\beta }_{1,n}:={\beta }_0+1-{\tilde{s}}_n.$$

(13)

• Update rule for the solution precision $\tau$

$$q\left(\tau \right|-)\sim \mathrm{Gamma}\left(a_0+\frac{N}{2},b_0+\frac{1}{2}(\parallel \tilde{\mathbf{x}}{\parallel }_2^2+Tr\left({{\Sigma }}_{\tilde{x}}\right))\right),$$

where ${{\Sigma }}_{\tilde{x}}=diag({\sigma }_{{\tilde{x}}_1}^2,\dots ,{\sigma }_{{\tilde{x}}_N}^2)$ and $Tr\left(A\right)$ is the trace of matrix A. Thus

$$\tilde{\tau }=\frac{a_0+\frac{N}{2}}{b_0+\frac{1}{2}\left(\parallel \tilde{\mathbf{x}}{\parallel }_2^2+{\sum }_{n=1}^N{\sigma }_{{\tilde{x}}_n}^2\right)}.$$

(14)

• Update rule for the noise precision $ϵ$

$$q\left(ϵ\right|-)\sim \mathrm{Gamma}({\theta }_0+\frac{M}{2},{\theta }_1+\frac{1}{2}\tilde{\Psi }),$$

where

$$\tilde{\Psi }:=\left({\mathbf{y}}^T\mathbf{y}-2{(\tilde{\mathbf{x}}\circ \tilde{\mathbf{s}})}^T{A}^T\mathbf{y}+Tr\left((\tilde{\mathbf{x}}{\tilde{\mathbf{x}}}^T+{{\Sigma }}_{\tilde{x}})\tilde{\Phi }\right)\right).$$

(15)

This yields to the following update rule for the precision of the noise component

$$\widetilde{ϵ}=\frac{{\theta }_0+\frac{M}{2}}{{\theta }_1+\frac{1}{2}\tilde{\Psi }}.$$

(16)

The stopping criterion of the algorithm is made based on the log-marginalized likelihood. We define the stopping condition in terms of $L:=log\left\{p\right(\mathbf{y}|\mathbf{s},ϵ,\tau )\}$. The marginalized likelihood can be written as

$$p\left(\mathbf{y}\right|\mathbf{s},ϵ,\tau )=\int p\left(\mathbf{y}\right|\mathbf{x},\mathbf{s},ϵ)p\left(\mathbf{x}\right|\tau I_N)d\mathbf{x}.$$

After some simplification, the negative log-likelihood is proportional to

$$-L\propto log|{{\Sigma }}_0^{-1}|+{\mathbf{y}}^T{{\Sigma }}_0\mathbf{y},$$

where

$${{\Sigma }}_0={({\widetilde{ϵ}}^{-1}{I}_M+{\tilde{\tau }}^{-1}A{\tilde{S}}^2{A}^T)}^{-1}$$

(17)

and $\tilde{S}:=\mathrm{diag}\left\{\tilde{\mathbf{s}}\right\}$. Therefore, the stopping condition can be made as

$${\Delta }{L}_n^{\left[t\right]}:=|{\Delta }{L}^{\left[t\right]}|/|L^{[t-1]}|\le T_0,$$

(18)

for some small value of threshold $T_0$ [50], where

$$L^{\left[t\right]}:=log{{\Sigma }}_0^{\left[t\right]}-{\mathbf{y}}^T{{\Sigma }}_0^{\left[t\right]}\mathbf{y}.$$

(19)

and

$$\begin{array}{cc}\hfill {\Delta }{L}^{\left[t\right]}:& =L^{\left[t\right]}-L^{[t-1]}\hfill \\ \hfill & =log|\frac{{{\Sigma }}_0^{\left[t\right]}}{{{\Sigma }}_0^{[t-1]}}|+{\mathbf{y}}^T({{\Sigma }}_0^{[t-1]}-{{\Sigma }}_0^{\left[t\right]})\mathbf{y}.\hfill \end{array}$$

(20)

Figure 1 illustrates the graphical Bayesian representation of the BGiG model, which is an undirected graph. The shaded node $\mathbf{y}$ shows the observations (measurements), and the small solid nodes represent the hyperparameters. Each unshaded node denotes a random variable (or a group of random variables).

The flowchart representation of the algorithm is shown in Figure 2 motivated by the graphical approach in [47,73]. According to the pseudocode in Algorithm 1 and the flowchart in Figure 2, first, the hyperparameters of the model are set. The support-related hyperparameters ${\alpha }_0$ and ${\beta }_0$ are suggested to be set to small values with ${\alpha }_0\ll {\beta }_0$ to encourage $\mathbf{s}$ to be sparse on the average. The hyperparameters $a_0$ and $b_0$ on the precision of the solution-value vector are also initialized and suggested to be small not to bias the estimation when the measurements are incorporated. The hyperparameters ${\theta }_0$ and ${\theta }_1$ on the precision of the noise are recommended to be of order ${10}^{-6}$ for high SNRs. For moderate and low SNRs, higher values are recommended. In the next step, all the main variables of the model are drawn i.i.d. from their corresponding prior distributions defined in (5)–(7). Then, the stopping condition is computed based on the log-marginalized likelihood in (19). In the main loop, all of the main variables of the model are updated via the expected values obtained from the VB inference. Specifically, we first update the support vector and the solution value components; then, the precisions of the solution vector and the noise are updated. Finally, the stopping criterion is computed through the measure of the log-marginalized likelihood of the observations. The pseudocode of the algorithm is provided below.

Algorithm 1: SBL(BGiG) Algorithm

${\widehat{\mathbf{x}}}_s=\tilde{\mathbf{x}}\circ \tilde{\mathbf{s}}$ $\left[\tilde{\mathbf{x}},\tilde{\mathbf{s}}\right]=\mathbf{SBL}-\mathbf{BGiG}(Y,A)$ Set the hyperparameters, i.e., (${\alpha }_0,{\beta }_0$), ($a_0,b_0$), and (${\theta }_0,{\theta }_1$) % Variables Initialization Draw $\tilde{\mathbf{s}}$ and $\tilde{\gamma }$ from (5) Draw $\tilde{\mathbf{x}}$ and $\tilde{\tau }$ from (6) Draw $\widetilde{ϵ}$ from (7) t = 1     % Iterator Compute $L^{\left[\mathrm{t}\right]}$ from (19) and set $L^{\left[0\right]}=0$ % Main Loop for Estimations While$\frac{|L^{\left[\mathrm{t}\right]}-L^{[\mathrm{t}-1]}|}{|L^{[\mathrm{t}-1]}|}\ge T_0$. For example $T_0={10}^{-6}$.        Compute ${\tilde{s}}_n$ from (8), $\forall n=1,\dots ,N$                         % (Support vector component )        Compute ${{\Sigma }}_{\tilde{x}}$ and $\tilde{\mathbf{x}}$ from (10)                                   % (Solution-value matrix component)        Compute ${\alpha }_{1,n}$ and ${\beta }_{1,n}$ from (13) $\forall n=1,\dots ,N$        % (Parameters of the hyperprior $\mathbf{\gamma }$)        t Compute $\tilde{\tau }$ from (14)                                         % (Precision on the solution)        Compute $\widetilde{ϵ}$ from (16)                                            % (Precision on the noise)        Compute $L^{\left[\mathrm{t}\right]}$ from (19) and then t = t + 1 End While

3.2. Issues with SBL(BGiG)

In this section, we show that the estimated solution using SBL(BGiG) algorithm is sensitive to support-related hyperparameters, i.e., ${\alpha }_0$ and ${\beta }_0$ in (5). We provide an example under three cases to demonstrate this issue. We generated a random scenario, where the true solution ${\mathbf{x}}_s\in {\mathbf{R}}^{100}$ has the sparsity level of $k=25$, that is, the true $\mathbf{x}$ (or $\mathbf{s}$) has k active elements. The active elements of $\mathbf{s}$ were drawn randomly. The nonzeros of ${\mathbf{x}}_s$, corresponding to the active locations of $\mathbf{s}$, were drawn from $\mathcal{N}(0,{\sigma }_x^2)$, with ${\sigma }_x^2=1$. Each entry of the sensing matrix A was drawn i.i.d. from the Gaussian distribution $\mathcal{N}(0,1)$, then normalized, so each column has the Euclidian norm of 1. The elements of measurement noise were drawn from $\mathcal{N}(0,{\sigma }^2)$ with $SNR=25$ dB, where $SNR:=20{log}_{10}({\sigma }_x/\sigma )$. The hyperparameters of $\tau$ and $ϵ$ were set to $a_0=b_0={10}^{-3}$ and ${\theta }_0={\theta }_1={10}^{-6}$, respectively. In Cases 1–3, we set the pair $({\alpha }_0,{\beta }_0)$ with low emphasis on the prior $(0.01,0.99)$, moderate emphasis $(0.1,0.9)$ and fairly high emphasis $(1.4,2)$, respectively.

From the top to the bottom row of Figure 3, Figure 4 and Figure 5, we illustrate the estimated results with the number of measurements set to 80, 60 and 40 (that is, the sample ratio $\lambda$ is 0.80, 0.60, and 0.40), respectively. In each row of Figure 3, Figure 4 and Figure 5 from left to right, we show the comparison between the measurements $\mathbf{y}$ and the computed measurements based on $\widehat{\mathbf{y}}=A(\tilde{\mathbf{s}}\circ \tilde{\mathbf{x}})$, the true signal ${\mathbf{x}}_s=\mathbf{s}\circ \mathbf{x}$ and the reconstructed signal ${\widehat{\mathbf{x}}}_s=\tilde{\mathbf{s}}\circ \tilde{\mathbf{x}}$, the true support vector $\mathbf{s}$ and the estimated support vector $\tilde{\mathbf{s}}$ and the evolution of the estimated supports with respect to the iterations in the SBL(BGiG) algorithm.

According to Figure 3, the setting for $({\alpha }_0,{\beta }_0)$ in Case 1 fails to provide perfect results even for high sampling ratios. Similarly, Figure 4 shows that the settings for $({\alpha }_0,{\beta }_0)$ in Case 2 do not provide encouraging results even for high sampling ratios. Specifically, it turns out that Case 1 and Case 2 provide sparse solutions for the sampling ratios within the range $[0,1]$, where $\lambda =1$ means $M=N$.

According to Figure 5, setting $({\alpha }_0,{\beta }_0)$ to $(1.4,2)$ seems to be a reasonable choice for high sampling ratios (over $70\%$), while it is not a good choice for the lower sampling ratios. This issue can be seen in the supports plot in the 2nd and 3rd row of Figure 5. One may argue that the estimated support vector $\widehat{\mathbf{s}}$ can be filtered via some threshold value (such as $0.3$) for $\lambda =0.6$. However, thresholding will adversely affect the detection rate, and setting the threshold depends on our understanding of the signal characteristics. Furthermore, we should account for the effect of the filtered supports since their corresponding estimated components in ${\widehat{\mathbf{x}}}_s$ contribute to fitting the model to the measurements.

In 

Table 1. Performance results of SBL(BGiG) for Cases 1–3.

Case 1: $({\mathit{\alpha }}_0=0.01,{\mathit{\beta }}_0=0.99)$				Case 2: $({\mathit{\alpha }}_0=0.1,{\mathit{\beta }}_0=0.9)$				Case 3: $({\mathit{\alpha }}_0=1.4,{\mathit{\beta }}_0=2.0)$
$\mathbf{\lambda }$	${\mathit{P}}_{\mathit{D}}$	${\mathit{P}}_{\mathit{FA}}$	NMSE (dB)	$\mathbf{\lambda }$	${\mathit{P}}_{\mathit{D}}$	${\mathit{P}}_{\mathit{FA}}$	NMSE (dB)	$\mathbf{\lambda }$	${\mathit{P}}_{\mathit{D}}$	${\mathit{P}}_{\mathit{FA}}$	NMSE (dB)
0.8	0.20	0	−2.367	0.8	0.24	0	−3.109	0.8	0.72	0	−16.264
0.6	0.08	0	−1.326	0.6	0.16	0	−2.197	0.6	1	0	−5.226
0.4	0.08	0	−1.181	0.4	0.08	0	−1.181	0.4	1	1	−0.088

, we summarize the performance of the generated example for Cases 1–3, where $P_D$, $P_{FA}$ and $NMSE$ denote the detection rate and false alarm rate in support recovery and the normalized mean-squared error between the true and the estimated sparse signal. This also shows that the algorithm fails to provide reasonable results for the sampling ratio of $\lambda =0.4$.

These experiments suggest that there is no fixed setting for $({\alpha }_0,{\beta }_0)$ capable of performing reasonably well for all sampling ratios and thus, selecting the hyperparameters $({\alpha }_0,{\beta }_0)$ should be made with care.

Continuing this examination, in Figure 6, Figure 7 and Figure 8, we illustrate the negative log-marginalized likelihood, the noise precision estimation and the estimated precision on the generated true solution in Cases 1–3, respectively. The horizontal axis shows the iterations until the stopping rule is met.

As expected, as the sampling ratio increases, the algorithm requires fewer iterations to meet its stopping condition. This can be seen on the negative log-marginalized likelihood plots in Figure 6, Figure 7 and Figure 8. In these experiments, the actual precision of the solution components was set to $\tau =1$, and the actual noise precision was set to $ϵ=316.2$.

For Cases 1 and 2, according to Figure 6, Figure 7 and Figure 8, the estimated precisions on both the noise and solution components were far off from the actual ones even for $\lambda =0.8$. Thus, it resulted in poor performance in signal recovery for Cases 1 and 2 (see Figure 3 and Figure 4).

For Case 3, the estimated precisions on the noise and the solution components were acceptable for $\lambda =0.8$ but far off from the actual ones for lower sampling ratios (see Figure 8). The main issue of the failures can be found in the update rule of the support learning vector $\tilde{\mathbf{s}}$ defined in (8). It is important to balance between the terms $c_n$ and ${\kappa }_n$, where $c_n$ imposes the effect of hyper-prior on $\mathbf{s}$ accompanied by the current estimate of $s_n$. In contrast, ${\kappa }_n$ imposes the contribution of the current estimates of noise precision, solution and other supports in fitting the model to the measurements. Therefore, if we impose a substantial weight on the sparsity via $c_n$, the solution tends to neglect the effect of ${\kappa }_n$ and vice versa. This is why we had sparse (with poor performance) in Cases 1 and 2 for all the represented sampling ratios and nonsparse (with poor performance) for moderate and lower sampling ratios in Case 3. These results suggest that the algorithm and its update rules are sensitive to the selection of hyperparameters on the Gamma prior on the support vector $\mathbf{s}$. The main issue can be seen in (9), where the selection of the hyperparameters ${\alpha }_0$ and ${\beta }_0$ resulted in a large or small value in $c_n$ due to the digamma function.

4. Gaussian-Inverse Gamma Modeling and SBL(GiG) Algorithm

In this section, we consider the Gaussian-inverse Gamma (GiG) model. In this model, each component $x_n$ of the solution is modeled by zero-mean Gaussian with the precision ${\tau }_n$. The main difference between this model and the model defined in Section 3 is that the GiG model does not have the support vector $\mathbf{s}$; instead, different precisions are considered on the components of the solution vector ${\mathbf{x}}_s$ in (1). A simpler version of GiG can also be used by defining the same precision $\tau$ for all the components of ${\mathbf{x}}_s$.

Here, we rather use different precisions to make the GiG model have almost the same complexity as the BGiG model in terms of the parameters to be learned. The set of priors in this model is defined as follows.

$$x_n\sim \mathcal{N}(0,{\tau }_n^{-1}),{\tau }_n\sim \mathrm{Gamma}(a_0,b_0),\forall n,$$

(21)

where $a_0$ and $b_0$ denote the shape and rate of the Gamma distribution, respectively. The entries of the noise component $\mathbf{e}$ are defined the same as (7), i.e., 

$$\mathbf{e}\sim \mathcal{N}(0,{ϵ}^{-1}{I}_M),ϵ\sim \mathrm{Gamma}({\theta }_0,{\theta }_1),$$

where ${\theta }_0$ and ${\theta }_1$ are set to small positive values. The estimation of the parameters in this model is carried out using VB inference, as discussed below.

4.1. Update Rules of SBL(GiG) Using VB Inference

According to the VB algorithm described in (2) and (3), the update rule of the variables and parameters of the GiG model can be simplified as follows. The details of these derivations appear in Appendix A.2.

• Update rule for the precision ${\tau }_n$ on $x_n$ using VB

$$q\left({\tau }_n\right)\sim \mathrm{Gamma}\left(a_0+\frac{1}{2},b_0+\frac{1}{2}({\tilde{x}}_n^2+{\sigma }_{{\tilde{x}}_n}^2)\right),\forall n=1,\dots ,N.$$

Thus,

$$\widetilde{{\tau }_n}=\frac{a_0+\frac{1}{2}}{b_0+\frac{1}{2}({\tilde{x}}_n^2+{\sigma }_{{\tilde{x}}_n}^2)},\forall n=1,2,\dots ,N.$$

(22)

• Update rule for the noise precision $ϵ$ using VB

$$q\left(ϵ\right)\sim \mathrm{Gamma}({\theta }_0+\frac{M}{2},b_0+\frac{1}{2}\tilde{\Psi })$$

which yields

$$\widetilde{ϵ}=\frac{{\theta }_0+\frac{M}{2}}{{\theta }_1+\frac{1}{2}\tilde{\Psi }},$$

(23)

where

$$\tilde{\Psi }:={\mathbf{y}}^T\mathbf{y}-2{\tilde{\mathbf{x}}}^T{A}^T\mathbf{y}+Tr\left((\tilde{\mathbf{x}}{\tilde{\mathbf{x}}}^T+{{\Sigma }}_{\tilde{x}})A^{T}A\right).$$

(24)

• Update rule for the solution vector $\mathbf{x}$ using VB

$$q_x\left(\mathbf{x}\right)\sim \mathcal{N}(\tilde{\mathbf{x}},{{\Sigma }}_{\tilde{x}}),$$

(25)

where

$${{\Sigma }}_{\tilde{x}}:={(\tilde{T}+\widetilde{ϵ}{A}^{T}A)}^{-1}\mathrm{and}\tilde{\mathbf{x}}:=\widetilde{ϵ}{{\Sigma }}_{\tilde{x}}{A}^T\mathbf{y},$$

(26)

and

$$\tilde{T}:=\mathrm{diag}\left\{[{\tilde{\tau }}_1,\dots ,{\tilde{\tau }}_N]\right\}.$$

We set the stopping rule of the algorithm using the marginalized likelihood (evidence) defined as

$$p\left(\mathbf{y}\right|ϵ,\tau )=\int p(\mathbf{y}|\mathbf{x},ϵ,\tau )p\left(\mathbf{x}\right|\tau )d\mathbf{x}.$$

After simplification and for the comparison purposes of $L^{\left[t\right]}$ with $L^{[t-1]}$ in the updating process, we have

$$L^{\left[t\right]}\propto log|{{\Sigma }}_0^{\left[t\right]}|-{\mathbf{y}}^T{{\Sigma }}_0^{\left[t\right]}\mathbf{y},$$

where ${{\Sigma }}_0$ is defined as

$${{\Sigma }}_0:={({\widetilde{ϵ}}^{-1}{I}_M+{\tilde{T}}^{-1}A{A}^T)}^{-1}.$$

(27)

Therefore, similar to SBL(BGiG), the stopping condition can be made as

$${\Delta }{L}_n^{\left[t\right]}:=|{\Delta }{L}^{\left[t\right]}|/|L^{[t-1]}|\le T_0,$$

(28)

for some small value of threshold $T_0$.

Figure 9 illustrates the graphical Bayesian representation of the GiG model, which is an undirected graph. Similar to Figure 1, the shaded node $\mathbf{y}$ shows the observations, the small solid nodes represent the hyperparameters and the unshaded nodes denote the random variables.

The flowchart representation of the algorithm is shown in Figure 10. According to the pseudocode in Algorithm 2 and the flowchart in Figure 10, first, the hyperparameters of the model are set. The hyperparameters $a_0$ and $b_0$ on the precision of the solution-value vector are initialized and suggested to be small. Similar to SBL(BGiG), the hyperparameters ${\theta }_0$ and ${\theta }_1$ on the precision of the noise are recommended to be of order ${10}^{-6}$ for high SNRs. All the main variables of the model are drawn i.i.d. from their corresponding prior distributions defined in (22)–(26). Then, the stopping condition is computed based on (28). In the main loop, all the main variables of the model are updated via the expected values obtained from the VB inference through (22)–(26). The pseudocode of the algorithm is provided below.      

Algorithm 2: SBL(GiG) Algorithm

${\tilde{\mathbf{x}}}_s=\mathbf{SBL}-\mathbf{GiG}(Y,A)$                 Set the hyperparameters, i.e., ($a_0,b_0$) and (${\theta }_0,{\theta }_1$) % Variables’ Initialization Draw ${\tilde{\mathbf{x}}}_s$ and $\widetilde{\mathsf{ø}}$ from (21) Draw $\widetilde{ϵ}$ from (7) t = 1     % Iterator Compute ${\tilde{L}}^{\left[\mathrm{t}\right]}$ from (28) and (27), and set ${\tilde{L}}^{\left[0\right]}=0$ % Main Loop for Estimations t = 1 While$\frac{|L^{\left[\mathrm{t}\right]}-L^{[\mathrm{t}-1]}|}{|L^{[\mathrm{t}-1]}|}\ge T_0$. For example $T_0={10}^{-6}$.        Compute ${{\Sigma }}_{\tilde{x}}$ and ${\tilde{\mathbf{x}}}_s$ from (26)          % (Solution-value matrix component)        Compute $\tilde{T}$ from (22)                    % (Precisions on the solution)        Compute $\widetilde{ϵ}$ from (23)                     % (Precision on the noise)        Compute $L^{\left[\mathrm{t}\right]}$ from (28) and (27), and then t = t + 1 End While

4.2. Issues with SBL(GiG)

An issue with the SBL(GiG) algorithm is that the solution becomes nonsparse since it does not incorporate a binary vector $\mathbf{s}$ (hard-thresholding or soft-thresholding if the expected value is used) as we had in SBL(BGiG). This may have no major effect on the signal reconstruction for high sampling ratios. However, the nonsparseness effect appears in low sampling ratios by misleading the algorithm to wrongly activate many components in the estimated signal yet providing a good fit of the model to the measurements. Here, we use the same example as we made for the SBL(BGiG) model with the same sensing matrix A, measurement vector $\mathbf{y}$ and noise $\mathbf{e}$. Notice that in the SBL(BGiG) model, we considered the same precision $\tau$ on all the components of the solution value vector $\mathbf{x}$ support vector $\mathbf{s}$. In contrast, the SBL(GiG) model does not have the support learning vector; instead, we assume that each component of the solution vector has different precision ${\tau }_n$. It turns out that SBL(GiG) is not very sensitive to the selection of the hyperparameters as the SBL(BGiG). Thus, here, we show the results for one case scenario for the hyperparameters. We use the same setting for the parameters of $ϵ$ in the hyper prior as before, i.e., ${\theta }_0={\theta }_1={10}^{-6}$, and the same parameters for all the precisions ${\tau }_n$ of the solution component, i.e., $a_0=b_0={10}^{-3}$. In Figure 11 and Figure 12, we illustrate the results after applying the SBL(GiG) algorithm. In Figure 11, from left to right, we show the results for sampling ratios of $\lambda =$0.8, 0.6, and 0.40, respectively. The first row shows the comparison of $\mathbf{y}$ with $\widehat{\mathbf{y}}=A{\tilde{\mathbf{x}}}_s$, the second row shows the true solution ${\mathbf{x}}_s$ and the estimated solution ${\tilde{\mathbf{x}}}_s$, and the third row demonstrates the estimated precisions on the solution components. In Figure 12, we demonstrate the negative log-marginalized likelihood comparison and the estimated noise precision against the true noise precision for the sampling ratios of $\lambda =$ 0.8, 0.6 and 0.4.

From the results shown in Figure 11 and Figure 12, we observe that the recovered signal tends to become nonsparse. This effect is illustrated in the second row of Figure 11. This can also be observed in the precision estimations of the solution components. More specifically, the true nonzero components in our simulations were drawn from a zero-mean Gaussian with the precision of ${\tau }_n=1$. Thus, the ideal precision estimation would be within the two classes of values of 1 and infinity or very large values. However, the estimated results in our simulation do not show such a classification. As the sampling ratio decreases, the solution estimate has poor performance, due not only to the reduction in the number of measurements but also the nonsparseness behavior.

5. Preprocessing versus Postprocessing and Simulations

In this section, we show that in order to improve the performance of Bernoulli–Gaussian-inverse Gamma modeling using the SBL(BGiG) algorithm, we need to perform a preprocessing step. The results in Section 4 suggest one can perform some postprocessing for the SBL(GiG) algorithm to improve the reconstruction performance. Below, we provide more details for each of these algorithms.

5.1. Pre-Processing for the SBL(BGiG) Algorithm

Based on the observations made on the performance of SBL(BGiG) in Section 3.2, we showed that the pair of hyperparameters $({\alpha }_0$,${\beta }_0)$ should be selected with care. In other words, obtaining good performance with this algorithm needs some preprocessing to assess an appropriate setting for the parameters. For a more rigorous study, here, we perform a grid search on the hyperparameters $({\alpha }_0,{\beta }_0)$ to see whether we can find some common pattern in selecting these parameters for all sampling ratios. The grid search runs the algorithm for different values of ${\alpha }_0$ and ${\beta }_0$ with the search range of $[0.1,2]$ with the resolution of $0.1$. For each $({\alpha }_0,{\beta }_0)$ within this range, we ran 200 random trials and then averaged the results. The settings of these trials are represented in

Table 2. Settings for preprocessing analysis and simulations on SBL(BGiG).

${\alpha }_0$	${\beta }_0$	$a_0$	$b_0$	${\theta }_0$	${\theta }_1$	Sparsity	$\gamma$	N
$[0.1,2]$	$[0.1,2]$	${10}^{-3}$	${10}^{-3}$	${10}^{-6}$	${10}^{-6}$	25	(5)	100
$\mathbf{s}$	$\tau$	$\mathbf{x}$	${\mathbf{x}}_s$	M	$ϵ$	$\mathbf{e}$	A	$\mathbf{y}$
(5)	(6)	(6)	${\mathbf{x}}_s=\mathbf{x}\circ \mathbf{s}$	$5:N$	316	(7)	${\left[A\right]}_{mn}\sim \mathcal{N}(0,1)$	$A{\mathbf{x}}_s+\mathbf{e}$

.

We generated a random scenario, where the true solution ${\mathbf{x}}_s\in {\mathbf{R}}^{100}$ has the sparsity level of $k=25$. The active elements of $\mathbf{s}$ were drawn randomly. The nonzeros of ${\mathbf{x}}_s$ were drawn from $\mathcal{N}(0,{\sigma }_x^2)$, with ${\sigma }_x^2=1$. Each entry of the sensing matrix A was drawn i.i.d. from the Gaussian distribution $\mathcal{N}(0,1)$, then normalized. The elements of measurement noise were drawn from $\mathcal{N}(0,{\sigma }^2)$ with $SNR=25$ dB. The results were examined to see what values of $({\alpha }_0,{\beta }_0)$ provided the highest performance in the detection rate vs. and false alarm rate. The simulation was executed for a range of sampling ratios in the range $[0.05,1]$ with the step size of 0.05. The results are demonstrated in Figure 13. In this figure, we also provide the results of performing a random Sobol search for $({\alpha }_0,{\beta }_0)$. A Sobol sequence is a low discrepancy quasirandom sequence. The two right plots in Figure 13 show the results for the best setting of $({\alpha }_0,{\beta }_0)$.

It should be clear from Figure 13 that there is no fixed setting for these parameters in order to get the best performance for all sampling ratios. The two plots on the right of Figure 13 illustrate the performance based on the best values of these hyperparameters, which provided the best performance, i.e., tuned hyperparameters. We also examined the grid search results for the top 10 highest performances for each sampling ratio, where performance is in terms of $P_D-P_{FA}$ and the normalized mean-squared error (NMSE). In Figure 14a, we demonstrate the top 10 highest performances based on NMSE and $P_D-P_{FA}$ for different sampling ratios. In Figure 14b,c, we illustrate the values of $({\alpha }_0,{\beta }_0)$, which led to the performances shown in Figure 14a for different sampling ratios. Figure 15 details the top 10 values of $({\alpha }_0,{\beta }_0)$ vs. sampling ratio.

According to Figure 14b,c, there is no specific pattern for these hyperparameters. Figure 15 also shows that hyperparameters need to be carefully selected.

5.2. Post-Processing for the SBL(GiG) Algorithm

Since the SBL(GiG) algorithm does not include the binary support vector $\mathbf{s}$, as SBL(BGiG) possesses, the resulting solution tends to become nonsparse. This leads to a high detection rate for the location of active supports and a high false alarm rate. Thus, as the sampling ratio decreases, there is a high chance that this algorithm overwhelms the locations of the true solution. Therefore, SBL(GiG) requires some postprocessing to discard the components with low amplitudes. This problem becomes of great importance for applications where detecting the correct nonzero components is more crucial than the magnitudes of the nonzeros in the signal. This effect can be seen in Figure 16b. The curves with solid lines in this plot show the detection and false alarm rate in support recovery and the difference between the rates. This issue can be resolved by some postprocessing such as data-driven threshold tuning. That way, the amplitudes in the reconstructed signal with lower values than the threshold can be discarded. For this purpose, we set up 200 random trials, the same way as the one explained for SBL(BGiG), and then evaluate the performance in terms of NMSE by varying the threshold. Figure 16b shows the averaged results of 200 trials. The settings of these trials are represented in

Table 3. Settings for preprocessing analysis and simulations on SBL(GiG).

$a_0$	$b_0$	${\theta }_0$	${\theta }_1$	Sparsity	N
${10}^{-3}$	${10}^{-3}$	${10}^{-6}$	${10}^{-6}$	25	100
${\tau }_n$	${\mathbf{x}}_s$	M	$\mathbf{e}$	A	$\mathbf{y}$
(22)	(25)	$5:N$	(23)	${\left[A\right]}_{mn}\sim \mathcal{N}(0,1)$	$A{\mathbf{x}}_s+\mathbf{e}$

.

In Figure 16a, we observe that the postprocessing does not benefit us so much in terms of the reconstruction error for low and moderate sampling ratios. However, there is a threshold of around 0.25, for which the postprocessing step reduced the reconstruction error by approximately 3 dB. We set the threshold to 0.25 and ran 200 random trials by applying SBL(GiG) and evaluating the performance based on the detection and false alarm rate in support recovery. According to Figure 16b, the additional post-processing step provides reasonable performance.

Finally, in Figure 17, we compare the performance of the SBL(BGiG) algorithm (with performing the preprocessing step) with the SBL(GiG) algorithm (after performing postprocessing). We see that Bernoulli–Gaussian-inverse Gamma implemented via SBL(BGiG) provides better performance for low and high sampling ratios. In contrast, Gaussian-inverse Gamma modeling implemented via SBL(GiG) performs much better for the moderate sampling ratios.

6. Conclusions

We investigated solving the inverse problem of compressive sensing using VB inference for two sparse Bayesian models of Bernoulli–Gaussian-inverse Gamma (BGiG) and Gaussian-inverse Gamma (GiG). The issues of each approach were discussed and the performance between the two models was compared. Specifically, we showed the behavior of these models and algorithms when the sampling ratio is low and moderate as well as the importance of selecting the hyperparameters of BGiG model with care. We further provided some intuition for performing additional pre/post-processing steps, depending on the selected model for better performance.

Based on our study on the synthetic data and considering the overall performance of both algorithms and the complexity in additional pre-/postprocessing, we observed that for moderate sampling ratios, SBL(GiG) is performing better than SBL(BGiG) modeling when using VB for sparse signals with no specific pattern in the supports. In contrast, SBL(BGiG) provided better perfomance for low and high sampling ratios. Finally, a rigorous comparison is required to study in the future under real-world scenarios and various applications. The MATLAB codes for GiG and BGiG modeling are available at https://github.com/MoShekaramiz/Compressive-Sensing-GiG-versus-BGiG-Modeling.git, accessed on 15 December 2022.