# Bayesian Sparse Estimation Using Double Lomax Priors

Authors: Xiaojing Gu; Henry Leung; Xingsheng Gu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Sparsity-promoting prior along with Bayesian inference is an effective approach in solving sparse linear models (SLMs). In this paper, we first introduce a new sparsity-promoting prior coined as Double Lomax prior, which corresponds to a three-level hierarchical model, and then we derive a full variational Bayesian (VB) inference procedure. When noninformative hyperprior is assumed, we further show that the proposed method has one more latent variable than the canonical automatic relevance determination (ARD). This variable has a smoothing effect on the solution trajectories, thus providing improved convergence performance. The effectiveness of the proposed method is demonstrated by numerical simulations including autoregressive (AR) model identification and compressive sensing (CS) problems.


---

## Body

## 1. Introduction

Consider the following linear model:(1)y=Φx+e,
where Φ∈ℜN×M is the designed matrix and y=[y1,…,yN]T is the measurement vector. e is the measurement noise, and x=[x1,…,xM]T is the unknown coefficient vector. This linear model is called underdetermined if M≤N and overdetermined otherwise. In the underdetermined case, there are in general numerous solutions, and additional assumptions need to be assumed on x for getting the desired one. In the applications we consider, the sparsity of x is an important assumption: most elements of x should be tiny or equal to zero when they are not required to tell the data, but few should be allowed to be large if necessary. A linear model with such an assumption on x is referred to the sparse linear model (SLM). SLMs can represent various types of signals in different transform domains [1–5]. Obtaining such a representation is equivalent to solving for a linear inverse problem with sparsity-promoting regularization. It is a fundamental issue of model (basis) selection [5, 6], compressive sensing [7, 8], source localization [9, 10], and data fusion [11, 12].Under the sparsity assumption, one seeks for a representation ofy with a minimal number of basis vectors in Φ. This is equivalent to seeking for the sparsest solution of x, which has a maximal number of elements equal to zero. The sparsest solution can be induced by the improper ℓ0-norm of x [13]: ∥x∥0≜∑m1{xm≠0}. The associated maximum a posteriori (MAP) estimation (assuming that e is white Gaussian noise) is a penalized least squares problem:
(2)minx∥y-Φx∥22+ρ∥x∥0,
where ρ is the model parameter that controls the relative importance applied to error term and sparseness term. The challenge in solving for (2) is that it requires a combinatorial search which is NP-hard in general. A natural approach is to replace ∥·∥0 by a tractable approximation. Among all the tractable approximations, ℓ1-norm is the tightest convex relaxation of ℓ0-norm; it leads to the successful Lasso algorithm [1, 3, 4, 14–16]. However, a nonconvex relaxation is able to recover sparsity in a more efficient way. It is proved that a nonconvex relaxation that closely resembles the ℓ0-norm is able to find the correct solution with fewer measurements [17–20].One difficulty with nonconvex relaxation is the existence of many local minima. Point-estimation methods such as MAP may easily get stuck in poor local solutions [21]. A posterior-estimation Bayesian approach can nicely alleviate this problem. Given the sparsity-promoting prior p(x), which is the probabilistic representation of a particular ℓ0-norm relaxation, the Bayesian inference can be expressed by
(3)p(x∣y)=p(y∣x)p(x)∫p(y∣x)p(x)dx,
where p(x∣y) is the posterior. The Bayesian approach computes the posterior mean E[x∣y]=∫xp(x∣y)dx instead of its mode, and it integrates instead of maximizing over p(x∣y). In case of many local minima, the mass of the probability density tends to concentrate at a strong minimum, which is more likely to be the global optimum. Although E[x∣y] is not necessarily the exact global optimum in general, it is proved that [22], by taking a “zero temperature limit,” E[x∣y] is able to converge to the global optimum. In other words, as -lnp(x) approaches the ℓ0-norm, E[x∣y] converges to the ℓ0 solution.In Bayesian formulation, the choice of priorp(x) is important for both solution sparsity and computational tractability. The Student-t prior is a relatively standard prior that stands out by being non-log-concave (results in nonconvex penalty) and having a Gaussian-integral representation for efficient processing. With 𝒮(·) denoting the Student-t prior, the distribution is defined as
(4)𝒮(x∣0,η,f)=Γ((f+1)/2)Γ(f/2)(ηπf)1/2[1+ηfx2]-(f+1)/2,η>0,f>0,
where f controls the degree of sparsity enforced and is known as shape parameter and η is the scale parameter. Additionally, denoting the Gaussian prior with variance τ by 𝒩(x∣0,τ) and the Gamma prior by 𝒢, defined as 𝒢(z∣a,b)=(ba/Γ(a))za-1exp(-bz), the Student-t prior corresponds to a Gaussian integral as follows:
(5)𝒮(x∣0,η,f)=∫0∞𝒩(x∣0,(ηξ)-1)𝒢(ξ∣f2,f2)dξ,η>0,f>0,
where ξ is the latent variable modifying the precision of the Gaussian mixture. The representation in (5) can be treated as an integral over scale parameter of the Gaussian density, thus is referred to as Gaussian scale mixtures (GSM) [23], which naturally invokes a two-level hierarchical Bayes model.Specifically, settingf→+0 leads to the noninformative-hyperpriors case of Student-t prior and gives the well known automatic relevance determination (ARD) prior (the scale parameter η is also assumed to be unity) [24] as follows:
(6)ARD(x)=𝒮(x∣0,1,f),f⟶0+.
As the limited case, the ARD prior provides the strongest sparsity-promoting potential among all Student-t priors with different values of f (see the remark of Proposition 1), and it has drawn tremendous attention in many problems [5, 8, 25]. Several methods have been proposed to solve for SLMs incorporating ARD priors, including variational Bayes (VB) inference [25], type II maximum likelihood (the “evidence framework”) [26], and a series of reweighted ℓ1 problems [27].In this paper, we formulate the SLM problem from a Bayesian perspective. Despite of the two-level Student-t prior, we propose the use of a new non-log-concave sparsity prior referred to as Double Lomax prior, which corresponds to a three-level hierarchical Bayes model. Additionally, a full VB inference procedure is derived to solve for the SLM using Double Lomax priors. As shown here, when noninformative hyperpriors are assumed (f→+0), the updating procedure of proposed method includes the canonical ARD updating procedure as a special case, but it has a better global convergence performance that can result in improved sparse estimations.The organization of this paper is as follows. In Section2, we present the proposed Double Lomax prior and demonstrate its properties. In Section 3, we formulate the Bayesian model of SLM using Double Lomax priors based on a generalized noise model. In Section 4, a full VB inference procedure is developed. In Section 5, we discuss the relationship between the proposed method and the ARD approach. Experimental results are presented in Section 6. Two SLM applications with different Φ matrix properties are considered in this section. One application indicates order determination and coefficients estimation for autoregressive (AR) models, where the rows of Φ matrix are strongly correlated. The other application is sparse signal recovery from the standard compressive sensing (CS) setup, where the rows of Φ matrix are uncorrelated. Conclusions are given in Section 7.


## 2. Double Lomax Prior and Its Properties

The proposed Double Lomax prior can be derived from two back-to-back spliced Lomax distributions (also known as Pareto II distribution [28]). That is,
(7)𝒟ℒ(x∣0,η,f)=η2(1+η|x|f)-(f+1),η>0,f>0.
It can be shown that 𝒟ℒ(x∣0,η,f)>0, ∫-∞+∞𝒟ℒ(x∣0,η,f)=1, and 𝒟ℒ(-x∣0,η,f)=𝒟ℒ(x∣0,η,f). Thus (7) is a probability distribution function of x defined on (-∞,+∞), symmetrical about zero. η is the scale parameter, f is the shape parameter, and 0 refers to zero mean.Moreover, note that(8)limf→+∞𝒟ℒ(x∣0,η,f)=η2exp(-η|x|),limf→+0𝒟ℒ(x∣0,η,f)={0,x≠0,η2,x=0.
Thus when the shape parameter goes to infinite, Double Lomax prior asymptotically becomes the Laplace prior with inverse scale η; when the shape parameter goes to zero, Double Lomax prior tends to strongly concentrate its mass on the zero value, distributing with the same summit but much longer and thinner tails. This behavior can be observed in Figure 1, which shows the Double Lomax priors with different values of f.Figure 1
Double Lomax prior at zero position, withη=1 for various values of f.Proposition 1.
Double Lomax prior in (7) is log-convex.Proof.
The logarithm of density function equation (7) is
(9)ln𝒟ℒ(x∣0,η,f)=ln(η2)-(f+1)ln(1+η|x|f).
It is continuous on (-∞+∞), and the derivatives are
(10)ddxln𝒟ℒ(x∣0,η,f)=-(f+1)ηf+η|x|,x≠0,d2dx2ln𝒟ℒ(x∣0,η,f)=(f+1)η2(f+η|x|)2{>0,f↛+∞,x≠0,=0,f⟶+∞,x≠0.
Thus, Double Lomax prior is log-convex for any f>0 and is strictly log-convex for f↛+∞, even though it does not have first and second derivative at the point x=0.Remark 2.
Proposition1 shows that the implicit penalty induced by Double Lomax prior, -ln𝒟ℒ(x∣0,η,f), is (strictly) concave, spanning from the convex ℓ1-norm penalty (when f→+∞) to the strongly concave ℓ0-norm-like penalty (when f→+0). Figure 2 is plotted -ln𝒟ℒ(x∣0,η,f) for different values of f.
For Student-t prior in (4), (d2/dx2)ln𝒮(x∣0,η,f)=η(f+1)(ηx2-f)/(f+ηx2)2, thus, the implicit penalty induced by Student-t prior, -ln𝒮(x∣0,η,f), has inconsistent concavity, that is, convex on the central interval [-f/ηf/η] and concave on each of the outer intervals, (-∞-f/η] and [f/η+∞). When f→+∞, the outer concave intervals are vanished, and the Student-t prior becomes the Gaussian prior with precision η and leads to the ℓ2-norm penalty that does not induce sparsity at all; when f→+0, the central convex interval is vanished, and the Student-t prior is now referred to as ARD prior, which becomes proportional to the noninformative Jeffreys prior p(x)∝1/|x| and leads to the log-sum penalty ∑log(|x|), providing similar sparsity-promoting effect as the Double Lomax prior does.Figure 2
Examples of the negative logarithm of Double Lomax prior for different values off, normalized by ln(1+η/f). The ℓ1-norm and ℓ0-norm are included for comparison.Proposition 3.
Double Lomax prior in (7) can be represented as a GSM. Given Gaussian prior 𝒩(·), exponential prior ℰ(·), and Gamma prior 𝒢(·), one has
(11)𝒟ℒ(x∣0,η,f)=∬0∞𝒩(x∣0,α)ℰ(α∣η2υ22)𝒢(υ∣f,f)dαdυ.Proof.
Equation (7) can be expressed as
(12)𝒟ℒ(x∣0,η,f)=η2ff+1(f+η|x|)-(f+1)=Γ(f+1)Γ(f)η2ff(f+η|x|)-(f+1)=ηff2Γ(f)∫0∞υfexp{-υ(f+η|x|)}dυ=∫0∞ηυ2exp{-ηυ|x|}ffΓ(f)υf-1exp{-fυ}dυ=∫0∞ℒ(x∣0,ηυ)𝒢(υ∣f,f)dυ.
Meanwhile, note that the Laplace prior has a hierarchical representation as follows:
(13)ℒ(x∣0,λ)=∫0∞𝒩(x∣0,α)ℰ(α∣λ22)dα=λ2exp(-λ|x|),
where ℰ(z∣τ)=τexp(-τz) is an exponential prior. Substituting (13) into (12) gives (11). This proves Proposition 3.Remark 4.
Equation (11) corresponds to a three-level hierarchical Bayes model with two latent variables α and υ. On the first level, random variables are generated from zero-mean Gaussian priors with independent variance α. On the second level, a set of α is generated from exponential priors with independent inverse scales, which are formed by (1/2)η2 multiplying a set of independent variables υ2. On the third level, the positive squares of independent parameters υ are generated from a Gamma prior with two parameters of the same value.
Moreover, besides this GSM model, Double Lomax prior can be considered as a Laplacian scaled mixture as well, and the inverse scaleυm can be thought as the width-controller of the corresponding Laplace prior: if υm is large, the corresponding distribution is narrow and the regularizing effect towards zero is strong.


## 3. Sparse Linear Model Using Double Lomax Priors

We assume independent Double Lomax priors over the coefficients and use the model in (1). Without loss of generality, the scale parameter η is assumed to be unity. Using (11), the coefficient model of (1) can be written as
(14)p(x∣α)=∏m=1M𝒩(xm∣0,αm)=𝒩(x∣0,Λ-1),(15)p(αm∣υm)=ℰ(αm∣υm22),(16)p(υm∣f)=𝒢(υm∣f,f),
where αm is the latent variable modifying the variance of the GSM model, υm is the latent variable modifying the inverse scale of exponential prior, and Λ=diag(αm-1) denotes the precision matrix.Traditionally the noise of SLM is presumed to be Gaussian distributed, which means that the measurements are also Gaussian distributed. However, in many real applications, measurements are distributed with tails that decay more slowly than Gaussian, which can not be described well by Gaussian noise model. Here we allow for very slow decay of tails by modeling the noise model using Student-t distribution. As a generalized noise model, Student-t can model a wide range of noises which from Gaussian to very heavy tailed. It also makes a connection of the robust statistic and the sparsity priors, under the view of the capability of modeling extreme events: those variables who have small probabilities but a large value range.Using a hierarchical Bayes model which is similar to (5), denoting the shape parameter by d (also called degrees of freedom), the scale parameter by β (also called precision parameter), and the latent variable by δn, and imposing two additional Gamma priors on β and d, the hierarchical form of Student-t distributed data can be written as
(17)p(y∣x,β,δ)=∏n=1N𝒩(yn∣ϕnx,(βδn)-1)=𝒩(y∣Φx,β-1L-1),(18)p(δn∣d)=𝒢(δn∣d2,d2),(19)p(β)=𝒢(β∣aβ,bβ),(20)p(d)=𝒢(d∣ad,bd),
where ϕn is the nth row of Φ and L=diag(δn) is the precision matrix.Figure3 shows the complete Bayesian model and the interdependencies among the model parameters. All hyperparameters in this graph are set to small values for getting flat hyperpriors (e.g., aβ=bβ=ad=bd=f=10-4). It can be observed from the graph that the joint distribution of the observations and parameters Θ={β,δ,d,x,α,υ} can be expressed as
(21)p(y,Θ)=p(y,β,δ,d,x,α,υ)=p(y∣x,β,δ)p(x∣α)p(α∣υ)×p(υ∣f)p(β)p(δ∣d)p(d),
where the p(y∣x,β,δ), p(x∣α), p(α∣υ), p(υ∣f), p(β), p(δ∣d) and p(d) are given in (17), (14), (15), (16), (19), (18), and (20), respectively.Figure 3
Directed graph representing the Bayesian modeling.


## 4. Variational Bayesian Sparse Regression

It is usually difficult to find the analytical results of the integrals in Bayesian estimation, hence approximation methods are required. There are two major classes of the approximations. One is Markov chain Monte Carlo (MCMC), which represents the posterior by samples; unbiased results can be obtained in infinite time. However, MCMC is relatively time consuming [29]. The other is variational relaxations [22], which relaxes the Bayesian estimation into a series of feasible optimization problems, providing a good tradeoff between computation accuracy and efficiency. VB inference [30–32] is one approach of this class.VB approximates the posterior densityp(Θ∣y) by introducing a factorable distribution q(Θ). Using Jensen’s inequality, we get the lower bound of the log probability of p(y):
(22)lnp(y)≥∫q(Θ)lnp(y,Θ)q(Θ)dΘ.
In this work, we factorize the approximation distribution q(Θ) over groups of parameters Θi as follows:
(23)q(Θ)=∏iq(Θi)=q(x)q(α)q(υ)q(β)q(δ)q(d).
Maximizing the right-hand side of (22) with respect to q(Θi), we can get a general expression of the optimal solution q^(Θi),
(24)q^(Θi)=exp〈lnp(y,Θ)〉k≠i∫exp〈lnp(y,Θ)〉k≠idΘi,
where 〈·〉k≠i denotes the expectation operator with respect to the distributions q(Θk) for all k≠i. ∫exp〈lnp(y,Θ)〉k≠idΘi is the integration over Θi which is a normalization constant. If the conjugate prior exists and is chosen for each group, the approximate posterior has the same functional form as the prior [30]. Hence the normalization constant does not need to be considered.Now we derive a full VB inference procedure for the above problem. We use the general optimal solution in (24) to obtain the factorized approximate posterior densities q(Θi). However, different from most VB updates in the literature, the challenge here is, for several random variables in the proposed model, that there are no conjugate priors for their likelihoods, and hence, in order to obtain an iterative algorithm with closed form, extra integration analysis has to be incorporated in our derivation.


### 4.1. Coefficients,x

To computeq(x), we only need to consider those terms that are functionally dependent on x. Because p(x∣α) is the conjugate prior for x, q^ would still be a Gaussian. It is convenient to ignore the normalization constant. The optimal solution is given by
(25)q^(x)∝exp〈lnp(y∣x,β,δ)+lnp(x∣α)〉β,δ,α∝exp〈ln𝒩(y∣Φx,β-1L-1)+ln𝒩(x∣0,Λ-1)〉β,δ,α∝𝒩(y∣Φx,〈β〉-1〈L〉-1)𝒩(x∣0,〈Λ〉-1)=𝒩(x∣μx,Σx),
where 〈L〉=diag(〈δn〉) and 〈Λ〉=diag(〈αm-1〉). Comparing the square term with standard Gaussian function, we can identify the mean and variance of this Gaussian, and get
(26)μx=〈β〉ΣxΦT〈L〉y,Σx=(〈β〉ΦT〈L〉Φ+〈Λ〉)-1.


### 4.2. Coefficient’s Latent Variables,α

From (14) and (15) we observe that the terms involving α are composed of M-order multiplicative terms over αm, xm, and υm, a further factorization of this approximate posterior density becomes
(27)q(α)=∏m=1Mq(αm),q^(αm)∝exp〈lnp(xm∣αm)+lnp(αm∣υm)〉xm,υm∝exp〈ln𝒩(xm∣0,αm)+lnℰ(αm∣υm22)〉xm,υm∝𝒩(〈xm2〉∣0,αm)ℰ(αm∣〈υm2〉2),
where
(28)〈xm2〉=μx(m)2+Σx(mm).·(m) and ·(mm) denote the mth element of a vector and the mth diagonal element of a matrix, respectively.However, since exponential distribution is not the conjugate prior for the variance of Gaussian distribution, the normalization term must be reinstated. Integrating overαm, the normalization constant of this posteriori becomes
(29)𝒞q^(αm)=∫0∞𝒩(〈xm2〉∣0,αm)ℰ(αm∣〈υm2〉2)dαm=ℒ(〈xm2〉∣0,〈υm2〉).
Therefore, the complete form of this posteriori is
(30)q^(αm)=𝒩(〈xm2〉∣0,αm)ℰ(αm∣〈υm2〉/2)ℒ(〈xm2〉∣0,〈υm2〉)=〈υm2〉2πexp(-〈υm2〉〈xm2〉)αm-1/2×exp(-12(〈xm2〉αm+〈υm2〉αm)).
It is noted that the posteriori density in (30) corresponds to a probability function known as generalized inverse Gaussian distribution (𝒢ℐ𝒢), which has the probability function [33]
(31)𝒢ℐ𝒢(z∣ω,χ,ψ)=(ψ/χ)ω/22Kω(χψ)z(ω-1)exp(-12(χz+ψz)),
where Kω(·) is the modified Bessel function of the third kind with index ω, one of its integral representations is given by
(32)Kω(τ)=12∫0∞τω-1exp(-12ω(1τ+τ))dω.
From (32), we have
(33)K1/2(τ)=π2τ-1/2exp(-τ).
Substituting (33) into (30) and comparing with (31), it can be found that
(34)q^(αm)=𝒢ℐ𝒢(αm∣12,〈xm2〉,〈λx〉).
The rth order moment of a 𝒢ℐ𝒢 distribution is given by [33]
(35)〈zr〉=Kω+r(χψ)Kω(χψ)(χψ)r/2.
In addition, from (32), it can be derived that
(36)Kω(τ)=K-ω(τ),Kω+1(τ)=2ωτKω(τ)+Kω-1(τ).
Therefore, for the specific 𝒢ℐ𝒢 in (34), using (35) with (36), we have:
(37)〈αm-1〉=〈υm2〉〈xm2〉,〈αm〉=(1〈xm2〉〈υm2〉+1)〈xm2〉〈υm2〉=〈αm-1〉-1+〈υm2〉-1.


### 4.3. Coefficient’s Latent Variables,υ

The approximate posterior density ofυ can be factorized as
(38)q(υ)=∏m=1Mq(υm),(39)q^(υm)∝exp〈lnp(αm∣υm)+lnp(υm∣f)〉αm∝exp〈lnℰ(αm∣υm22)+ln𝒢(υm∣f,f)〉αm∝exp((f+1)lnυm-〈αm〉2υm2-fυm).
Equation (39) contains logarithmic, square, and linear terms, thus it is not in correspondence with any standard distribution. Using a definite integration result in [34], the normalization term of q^(υm) can be expressed as
(40)𝒞q^(υm)=〈αm〉-(f+2)/2Γ(f+2)×exp(f24〈αm〉)D-(f+2)(f〈αm〉),
where Dω(·) is the parabolic cylinder function, whose computation can be implemented [35]. Thereby,
(41)〈υm2〉=1Cq^(υm)∫0∞υm2·υmf+1exp(-〈αm〉2υm2-fυm)dυm=1〈αm〉Γ(f+4)Γ(f+2)D-(f+4)(f/〈αm〉)D-(f+2)(f/〈αm〉).However, when noninformative hyperprior is assumed,f is set to zero and q^(υm) becomes the Rayleigh distribution who has the probability function
(42)ℛ(z∣λ)=zλ2exp(-z22λ2).
Here q^(υm) is ℛ(υm|1/〈xm2〉)=〈αm〉υmexp(-(〈αm〉/2)υm2). Therefore,
(43)〈υm2〉=∫0∞υm2·〈αm〉υmexp(-〈αm〉2υm2)dυm=〈αm〉∫0∞υm3exp(-〈αm〉2υm2)dυm=2〈αm〉.
Equation (43) indicates that the value of 〈υm2〉 is simply inversely proportional to 〈αm〉 when f→+0.


### 4.4. Noise Precision Parameter,β

Since Gamma distribution is conjugate to the precision of Gaussian function, we still obtain a Gamma distribution forβ as follows:
(44)q^(β)=𝒢(β∣Aβ,Bβ),
where
(45)Aβ=N2+aβ,Bβ=12trace(〈L〉〈(y-Φx)(y-Φx)T〉)+bβ,〈(y-Φx)(y-Φx)T〉=yyT-2ΦμxyT+Φ(μxμxT+Σx)ΦT.


### 4.5. Noise’s Latent Variables,δ

It is observed that the approximate posterior density ofδ can be further factorized as
(46)q(δ)=∏n=1Nq(δn).
And q^(δn) is still a Gamma distribution as follows:
(47)q^(δn)=𝒢(δn∣Aδ,Bδn),
where
(48)Aδ=12+〈d〉2,Bδn=〈β〉2〈(yn-ϕnx)2〉+〈d〉2,〈(yn-ϕnx)2〉=yn2-2ynϕnμx+ϕn(Σx+μxμxT)ϕnT.


### 4.6. Noise’s Degrees of Freedom,d

It is noted that(49)q^(d)∝exp(d2(∑n=1N‍〈lnδn〉-∑n=1N‍〈δn〉-2bd)-NlnΓ(d2)+N2dlnd2+(ad-1)lnd+d2(∑n=1N‍〈lnδn〉-∑n=1N‍〈δn〉-2bd)),
which is not in correspondence with any standard distribution. This problem has been solved by expanding lnΓ(d/2) using Sterling’s approximation [6], lnΓ(d/2)≈(d/2-1/2)ln(d/2)-d/2. Thus,
(50)q^(d)∝exp(∑n=1N〈lnδn〉-(N2+ad-1)lnd+12(N+∑n=1N〈lnδn〉-∑n=1N〈δn〉-2bd)d)=𝒢(d∣Ad,Bd),
where
(51)Ad=N2+ad,Bd=-12(N+∑n=1N‍〈lnδn〉-∑n=1N‍〈δn〉)+bd,
where
(52)〈lnδn〉=ψ(Aδ)-ln(Bδn),ψ(·) is the digamma function, defined by ψ(z)=(d/dz)lnΓ(z).We summarize the VB inference of Double Lomax formulation (noninformative hyperprior) as follows.For coefficient model, consider(53)q^(x)=𝒩(x∣μx,Σx),(54)μx=〈β〉ΣxΦT〈L〉y,(55)Σx=(〈β〉ΦT〈L〉Φ+〈Λ〉)-1,(56)q^(αm)=𝒢ℐ𝒢(αm∣12,〈xm2〉,〈λx〉),(57)〈αm-1〉=〈υm2〉〈xm2〉,(58)q^(υm)=ℛ(υm∣1〈xm2〉),(59)〈υm2〉=2〈αm〉=2〈αm-1〉-1+〈υm2〉-1.
For noise (observation) model, consider
(60)q^(β)=𝒢(β∣Aβ,Bβ),(61)〈β〉=AβBβ,(62)q^(δn)=𝒢(δn∣Aδ,Bδn),(63)〈δn〉=AδBδn,(64)q^(d)=𝒢(d∣Ad,Bd),(65)〈d〉=AdBd.
When initial expectations are assumed, the iterative updating can be carried out in closed form.


## 4.1. Coefficients,x

To computeq(x), we only need to consider those terms that are functionally dependent on x. Because p(x∣α) is the conjugate prior for x, q^ would still be a Gaussian. It is convenient to ignore the normalization constant. The optimal solution is given by
(25)q^(x)∝exp〈lnp(y∣x,β,δ)+lnp(x∣α)〉β,δ,α∝exp〈ln𝒩(y∣Φx,β-1L-1)+ln𝒩(x∣0,Λ-1)〉β,δ,α∝𝒩(y∣Φx,〈β〉-1〈L〉-1)𝒩(x∣0,〈Λ〉-1)=𝒩(x∣μx,Σx),
where 〈L〉=diag(〈δn〉) and 〈Λ〉=diag(〈αm-1〉). Comparing the square term with standard Gaussian function, we can identify the mean and variance of this Gaussian, and get
(26)μx=〈β〉ΣxΦT〈L〉y,Σx=(〈β〉ΦT〈L〉Φ+〈Λ〉)-1.


## 4.2. Coefficient’s Latent Variables,α

From (14) and (15) we observe that the terms involving α are composed of M-order multiplicative terms over αm, xm, and υm, a further factorization of this approximate posterior density becomes
(27)q(α)=∏m=1Mq(αm),q^(αm)∝exp〈lnp(xm∣αm)+lnp(αm∣υm)〉xm,υm∝exp〈ln𝒩(xm∣0,αm)+lnℰ(αm∣υm22)〉xm,υm∝𝒩(〈xm2〉∣0,αm)ℰ(αm∣〈υm2〉2),
where
(28)〈xm2〉=μx(m)2+Σx(mm).·(m) and ·(mm) denote the mth element of a vector and the mth diagonal element of a matrix, respectively.However, since exponential distribution is not the conjugate prior for the variance of Gaussian distribution, the normalization term must be reinstated. Integrating overαm, the normalization constant of this posteriori becomes
(29)𝒞q^(αm)=∫0∞𝒩(〈xm2〉∣0,αm)ℰ(αm∣〈υm2〉2)dαm=ℒ(〈xm2〉∣0,〈υm2〉).
Therefore, the complete form of this posteriori is
(30)q^(αm)=𝒩(〈xm2〉∣0,αm)ℰ(αm∣〈υm2〉/2)ℒ(〈xm2〉∣0,〈υm2〉)=〈υm2〉2πexp(-〈υm2〉〈xm2〉)αm-1/2×exp(-12(〈xm2〉αm+〈υm2〉αm)).
It is noted that the posteriori density in (30) corresponds to a probability function known as generalized inverse Gaussian distribution (𝒢ℐ𝒢), which has the probability function [33]
(31)𝒢ℐ𝒢(z∣ω,χ,ψ)=(ψ/χ)ω/22Kω(χψ)z(ω-1)exp(-12(χz+ψz)),
where Kω(·) is the modified Bessel function of the third kind with index ω, one of its integral representations is given by
(32)Kω(τ)=12∫0∞τω-1exp(-12ω(1τ+τ))dω.
From (32), we have
(33)K1/2(τ)=π2τ-1/2exp(-τ).
Substituting (33) into (30) and comparing with (31), it can be found that
(34)q^(αm)=𝒢ℐ𝒢(αm∣12,〈xm2〉,〈λx〉).
The rth order moment of a 𝒢ℐ𝒢 distribution is given by [33]
(35)〈zr〉=Kω+r(χψ)Kω(χψ)(χψ)r/2.
In addition, from (32), it can be derived that
(36)Kω(τ)=K-ω(τ),Kω+1(τ)=2ωτKω(τ)+Kω-1(τ).
Therefore, for the specific 𝒢ℐ𝒢 in (34), using (35) with (36), we have:
(37)〈αm-1〉=〈υm2〉〈xm2〉,〈αm〉=(1〈xm2〉〈υm2〉+1)〈xm2〉〈υm2〉=〈αm-1〉-1+〈υm2〉-1.


## 4.3. Coefficient’s Latent Variables,υ

The approximate posterior density ofυ can be factorized as
(38)q(υ)=∏m=1Mq(υm),(39)q^(υm)∝exp〈lnp(αm∣υm)+lnp(υm∣f)〉αm∝exp〈lnℰ(αm∣υm22)+ln𝒢(υm∣f,f)〉αm∝exp((f+1)lnυm-〈αm〉2υm2-fυm).
Equation (39) contains logarithmic, square, and linear terms, thus it is not in correspondence with any standard distribution. Using a definite integration result in [34], the normalization term of q^(υm) can be expressed as
(40)𝒞q^(υm)=〈αm〉-(f+2)/2Γ(f+2)×exp(f24〈αm〉)D-(f+2)(f〈αm〉),
where Dω(·) is the parabolic cylinder function, whose computation can be implemented [35]. Thereby,
(41)〈υm2〉=1Cq^(υm)∫0∞υm2·υmf+1exp(-〈αm〉2υm2-fυm)dυm=1〈αm〉Γ(f+4)Γ(f+2)D-(f+4)(f/〈αm〉)D-(f+2)(f/〈αm〉).However, when noninformative hyperprior is assumed,f is set to zero and q^(υm) becomes the Rayleigh distribution who has the probability function
(42)ℛ(z∣λ)=zλ2exp(-z22λ2).
Here q^(υm) is ℛ(υm|1/〈xm2〉)=〈αm〉υmexp(-(〈αm〉/2)υm2). Therefore,
(43)〈υm2〉=∫0∞υm2·〈αm〉υmexp(-〈αm〉2υm2)dυm=〈αm〉∫0∞υm3exp(-〈αm〉2υm2)dυm=2〈αm〉.
Equation (43) indicates that the value of 〈υm2〉 is simply inversely proportional to 〈αm〉 when f→+0.


## 4.4. Noise Precision Parameter,β

Since Gamma distribution is conjugate to the precision of Gaussian function, we still obtain a Gamma distribution forβ as follows:
(44)q^(β)=𝒢(β∣Aβ,Bβ),
where
(45)Aβ=N2+aβ,Bβ=12trace(〈L〉〈(y-Φx)(y-Φx)T〉)+bβ,〈(y-Φx)(y-Φx)T〉=yyT-2ΦμxyT+Φ(μxμxT+Σx)ΦT.


## 4.5. Noise’s Latent Variables,δ

It is observed that the approximate posterior density ofδ can be further factorized as
(46)q(δ)=∏n=1Nq(δn).
And q^(δn) is still a Gamma distribution as follows:
(47)q^(δn)=𝒢(δn∣Aδ,Bδn),
where
(48)Aδ=12+〈d〉2,Bδn=〈β〉2〈(yn-ϕnx)2〉+〈d〉2,〈(yn-ϕnx)2〉=yn2-2ynϕnμx+ϕn(Σx+μxμxT)ϕnT.


## 4.6. Noise’s Degrees of Freedom,d

It is noted that(49)q^(d)∝exp(d2(∑n=1N‍〈lnδn〉-∑n=1N‍〈δn〉-2bd)-NlnΓ(d2)+N2dlnd2+(ad-1)lnd+d2(∑n=1N‍〈lnδn〉-∑n=1N‍〈δn〉-2bd)),
which is not in correspondence with any standard distribution. This problem has been solved by expanding lnΓ(d/2) using Sterling’s approximation [6], lnΓ(d/2)≈(d/2-1/2)ln(d/2)-d/2. Thus,
(50)q^(d)∝exp(∑n=1N〈lnδn〉-(N2+ad-1)lnd+12(N+∑n=1N〈lnδn〉-∑n=1N〈δn〉-2bd)d)=𝒢(d∣Ad,Bd),
where
(51)Ad=N2+ad,Bd=-12(N+∑n=1N‍〈lnδn〉-∑n=1N‍〈δn〉)+bd,
where
(52)〈lnδn〉=ψ(Aδ)-ln(Bδn),ψ(·) is the digamma function, defined by ψ(z)=(d/dz)lnΓ(z).We summarize the VB inference of Double Lomax formulation (noninformative hyperprior) as follows.For coefficient model, consider(53)q^(x)=𝒩(x∣μx,Σx),(54)μx=〈β〉ΣxΦT〈L〉y,(55)Σx=(〈β〉ΦT〈L〉Φ+〈Λ〉)-1,(56)q^(αm)=𝒢ℐ𝒢(αm∣12,〈xm2〉,〈λx〉),(57)〈αm-1〉=〈υm2〉〈xm2〉,(58)q^(υm)=ℛ(υm∣1〈xm2〉),(59)〈υm2〉=2〈αm〉=2〈αm-1〉-1+〈υm2〉-1.
For noise (observation) model, consider
(60)q^(β)=𝒢(β∣Aβ,Bβ),(61)〈β〉=AβBβ,(62)q^(δn)=𝒢(δn∣Aδ,Bδn),(63)〈δn〉=AδBδn,(64)q^(d)=𝒢(d∣Ad,Bd),(65)〈d〉=AdBd.
When initial expectations are assumed, the iterative updating can be carried out in closed form.


## 5. Relations to ARD Update Procedure

Now let us inspect the relations between the proposed Double Lomax updating procedure and the ARD updating procedure.Considering (57) and (59) that are induced by the Double Lomax priors, it can be shown that if just (57) and (59) are iterated to convergence, then the final result is equivalent to directly setting
(66)〈αm-1〉=1〈xm2〉,
which is the exact ARD update [26]. To prove this, construct a function as ℱ(〈αm-1〉,〈υm2〉)=〈xm2〉〈αm-1〉+〈υm2〉/〈αm-1〉-ln(〈υm2〉). It can be shown that (57) minimizes ℱ partially with respect to 〈αm-1〉, and (59) minimizes ℱ partially with respect to 〈υm2〉. If we do the complete minimization over both 〈αm-1〉 and 〈υm2〉, then we directly get (66). Therefore, when noninformative hyperprior is assumed, the proposed Double Lomax updates only have one more latent variable 〈υm2〉 than the canonical ARD updates; this latent variable can be further vanished if we only consider the updating equations of (57) and (59).However, if we consider the full iteration procedure involving all the updating equations, (54), (55), (57), (59), (61), (63), and (65), then the difference is shown: the solution trajectories of Double Lomax are not necessarily consistent with those of ARD. To see this, we need to denote ARD’s Equation (66), and Double Lomax’s Equations (57) and (59) from an iteration perspective, so we have
(67)〈αm-1〉(l+1)=1〈xm2〉(l)
for ARD and
(68)〈αm-1〉(l+1)=〈υm2〉(l)〈xm2〉(l),(69)〈υm2〉(l+1)=2〈αm-1〉(l)-1+〈υm2〉(l)-1
for Double Lomax, where l denotes the lth iteration.It can be seen from (69) that the latent variable 〈υm2〉(l) of Double Lomax can be expressed as a function of 〈αm-1〉, which involves all the updates of 〈αm-1〉 from former iterations. That is,
(70)〈υm2〉(l)=Ω(l)(〈αm-1〉(l-1),〈αm-1〉(l-2),…,〈αm-1〉(0)).
Considering the above expression and comparing (67) with (68), the difference between ARD and Double Lomax (f→+0) is shown: in contrast to the ARD’s 〈αm-1〉(l+1) that depends only on 〈xm2〉(l), the Double Lomax’s 〈αm-1〉(l+1) depends not only on 〈xm2〉(l) but also on all the former updates of 〈αm-1〉, including 〈αm-1〉(l-1),〈αm-1〉(l-2),…,〈αm-1〉(0). Consequently, it makes the Double Lomax method being able to smooth the solution trajectory of 〈υm2〉(l) and also the trajectories of all other variables, as the results of smoothness propagation effect caused by the alternating-update strategy of VB inference. With these smooth trajectories, the algorithm may be prevented from a large deviation of the right direction when there are some disturbances or unreliable updates. As will be shown experimentally, comparing to ARD, Double Lomax tends to move more carefully around local minima and provides a better convergence performance.Following, we use linear regression as an example, the coefficient vector isx=[1,-1,0,0,1,-1,1,0,0,-1]T, the basis matrix Φ is a 12 × 10 random matrix whose entries are independently chosen from 𝒩(0,1), and e is zero-mean Gaussian noise with a standard deviation of 0.01. The VB-Double Lomax (f→+0) formulation established in the above section (denoted with DL) and an equivalent formulation with ARD prior (denoted with ARD) are used here for comparison. The maximum number of iterations are set to 50 for both of the algorithms. Figure 4 shows the estimation results of 50 trials, and the estimation mean (denoted with mean) and standard deviation (denoted with std) are also reported. Neither DL nor ARD can recover the exact coefficients due to the noise. In general, DL has shown more cases of convergence to global minima, while ARD has shown more cases of convergence to local minima, which has a relatively large deviation of global minima though with a small probability since that does not change the mean values much. Figure 5(a) reports the coefficient update traces of a typical case when ARD converges to a local minima but DL does not. Figure 5(b) reports the coefficient update traces of a typical case when both ARD and DL converge to the global minima. In both situations, the DL appears to have smoother coefficient update traces, which is in agreement with the above analysis. Furthermore, the update traces of latent variables 〈αm-1〉 and 〈υm2〉 and noise parameter β when both ARD and DL converge to the global minima (the case in Figure 5(b)) are also reported in Figure 6. Comparing to the traces of ARD’s 〈αm-1〉, smoother traces of DL’s are observed, even if they are plotted on a log scale. It is noted that the trace shapes of DL’s 〈αm-1〉 are similar to those of DL’s 〈υm2〉. But the differences can be seen in Figures 6(d) and 6(e) (nonzero coefficient’s traces are shown in Figure 6(d) while zero coefficients are shown in Figure 6(e) for better visualization). If we force 〈υm2〉 to equal to be 〈αm-1〉, an instance of ARD update is obtained. It is worth noting that since an extra level of iterations is introduced by 〈υm2〉, a slower convergence speed of DL can be expected in general.Figure 4
Regression coefficient estimation results using the Double Lomax formulation (DL) and the ARD formulation (ARD), reported with means and standard deviations.The update traces of coefficients from the regression coefficient estimation experiment shown in Figure4. (A) The update traces by ARD. (B) The update traces by DL. (C) The update traces of the first three coefficients [1,-1,0]T.
(a)
The update traces when ARD converges to a local minima but DL does not
(b)
The update traces when both ARD and DL converge to the global minimaThe update traces of latent variables and noise parameter when both ARD and DL converge to the global minima (shown in Figure5(b)).
(a)
log
〈
α
m
-
1
〉 in ARD
(b)
log
〈
α
m
-
1
〉 in DL
(c)
log
〈
υ
m
2
〉 in DL
(d)
〈
α
m
-
1
〉
-
〈
υ
m
2
〉
(e)
〈
α
m
-
1
〉
-
〈
υ
m
2
〉
(f)
〈
β
〉
-
0.5


## 6. Computer Simulations

In this section, the proposed VB-Double Lomax (f→+0) formulation (denoted with DL) is compared with the equivalent formulation with ARD prior (denoted with ARD). Furthermore, a hybrid formulation (denoted by hybrid), with first K iterations employing DL and the rest employing ARD until convergence, is also considered here for a tradeoff between performance and speed. This hybrid approach is based on the observation that ARD tends to converge to a local minima at early iterations. Two kinds of SLM applications with different Φ matrices are considered here. One is the AR model identification where the rows of Φ matrix are strongly correlated, the other is the standard compressive sensing problem where the rows of Φ matrix are uncorrelated. In the experiments reported below, the hybrid number K=20, and an experimental criterion, |∥x(k)∥-∥x(k-1)∥|<10-3 and lasting for more than 50 iterations, is used to terminate the iterative procedure.Autoregressive model is widely used to model time series data, which attempts to predict the future values based on previous values. An AR model of orderp is defined as(71)yn=∑i=1p‍xiyn-i+en,
where yn is the nth observation (measurement) in the series, xi is the AR coefficients, and en is the excitation noise. Using matrix notation, we have
(72)y=ΦARx+e,
where ΦAR is the N×p matrix whose nth row contains the lags for elements yn, that is, [yn-1,…,yn-p] and x is the AR coefficient vector to be estimated. If the model order p is unknown, the strategy used here is choosing an M where M is large enough to ensure M>p then computing the coefficient vector of M×1, and the extra coefficients are supposed to have estimates of zero.Figure7 shows an example of AR identification using DL with N=300, p=8, and M=20. The excitation noise follows a Student-t distribution with d=3 and β=1. It is shown in the left of the figure that for coefficients M≤p, the algorithm gives relatively satisfactory estimates, and for the extra coefficients M>p, the estimates have been set closely to zero. The center is plotted for the expected values of one-step-ahead predictions compared with the observations, and the right gives the comparison of the estimated excitation noise distribution with that used to generate the observations.Figure 7
AR model identification using DL. Left: coefficients and order estimation. Center: expected values of one-step-ahead predictions versus the observations. Right: excitation noise estimation.To compare DL, ARD, and hybrid, we use three AR models withp=5, 10, and 15. In each case, M=20 and the excitation noise is Student-t distributed with d=3 and β=1. The coefficient estimation error is calculated as ∥x^-x∥2/∥x∥2, where x^ and x are the estimated and true coefficient vectors. The prediction mean square errors (MSEs) are calculated as E[(y^-y)2], where y^ and y are the one-step-ahead predictions and actual observations. Each experiment is carried out 200 times and the average results for AR(5), AR(10), and AR(15) are shown in Figures 8(a), 8(b), and 8(c), respectively. The figures on the left report the coefficient estimation error with error bars (one standard deviation). It is shown that DL results in smaller errors than ARD, and the improvement is more significant for AR(10) and AR(15). It is also observed that hybrid gives almost the same performance as DL. As the number of observations increases, the coefficient estimation errors of DL, ARD, and Hybrid decrease and their differences become less significant. Similar observations are found in the prediction MSEs, which are reported in the center of Figure 8. However, the prediction MSEs do not decrease monotonously with the number of observations. This may be due to the reason that the current SLM formulation is based on the overall coefficient estimation error but does not put the partial autocorrelation of the AR model in consideration. Finally, the right of the figure reports the number of iterative steps required for convergence versus the number of observations. It is shown that for the same number of observations, DL requires more iterations than ARD for convergence. However, hybrid can achieve a comparable convergence rate as ARD. In general, as the number of observations increases, the convergence speeds of DL, ARD, and Hybrid also increase.Comparison of DL, ARD, and hybrid for AR(5), AR(10), and AR(15). Left: coefficient estimation error with error bars (one standard deviation). Center: prediction MSE. Right: the number of iterative step required for convergence.
(a)
AR(5): Actual coefficients [1.31, −1.53, 1.53, −1.23, 0.86]
(b)
AR(10): Actual coefficients [1.65, −2.59, 2.94, −3.78, 3.53, −3.26, 2.15, −1.61, 0.88, −0.47]
(c)
AR(15): Actual coefficients [1.82, −2.61, 3.04, −3.02, 3.33, −3.42, 3.55, −3.45, 2.95, −2.28, 1.59, −1.05, 0.63, −0.32, 0.09]CS is a technique for acquiring and recovering sparse signals with underdetermined linear systems that can be represented as(73)y=ΦCSx+e,
where y is N×1 observations, ΦCS is N×M measurement matrix, N<M, e represents the acquisition noise, and x is the M×1 unknown sparse signal, that is, most of its coefficients are zero.We use the following default CS setup in our experiments.ΦCS is induced by sampling i.i.d. columns from a unit sphere in RN. T coefficients at random locations of the signal are drawn from four different probability distributions, and the rest (M-T) of the coefficients are set to zero. The nonzero coefficients of the sparse signals are realizations of the following four distributions: (1) uniform ±1 random spikes, (2) nonnegative unit spikes, (3) zero-mean unit variance Gaussian, and (4) zero-mean Student-t with unit scale parameter and shape parameter equal 3. We fix M=256 and T=80, and we vary the number of measurements N. Moreover, both noiseless and noisy acquisitions are considered, where a zero-mean white Gaussian noise with a standard deviation of 0.005 is used in the noisy experiment. The reconstruction error is calculated as ∥x^-x∥2/∥x∥2, where x^ and x are the estimated and true signals. Each experiment is carried out 200 times and the average results are reported. Figure 9 reports the number of measurements versus reconstruction error with error ranges (one standard deviation) for the noise-free case. It is shown that DL provides improved overall reconstruction performance over ARD for a reasonable number of measurements, the improvement is more significant for signal types (1) and (2) and less significant for types (3) and (4) (but the differences can still be seen from the error ranges). And the hybrid can provide a performance close to DL. Figure 10 reports the number of iterative steps required for convergence versus the number of measurements for the noise-free case. DL shows a slower convergence rate than ARD, and hybrid has a comparable rate as of the ARD. Figure 11 shows the number of estimated nonzero coefficients for the noise-free case, whose values are larger than 10-3. Overall, DL gives the sparsest estimation results, then hybrid followed by ARD. Figure 12 reports the results for noisy acquisitions. Similar observations to the noise-free case (shown in Figure 9) are found in this figure, except that this time none of the algorithms can recover the coefficients exactly due to the noise.Number of measurements versus reconstruction error with error ranges (one standard deviation) for the noise-free case.
(a)
Uniform±1 random spikes
(b)
Nonnegative unit spikes
(c)
Zero-mean Gaussian
(d)
Zero-mean Student-tThe number of iterative step required for convergence for the noise-free case.
(a)
Uniform±1 random spikes
(b)
Nonnegative unit spikes
(c)
Zero-mean Gaussian
(d)
Zero-mean Student-tThe number of estimated nonzero coefficients for the noise-free case.
(a)
Uniform±1 random spikes
(b)
Nonnegative unit spikes
(c)
Zero-mean Gaussian
(d)
Zero-mean Student-tNumber of measurements versus reconstruction error with error ranges (one standard deviation) for noisy observations.
(a)
Uniform±1 random spikes
(b)
Non-negative unit spikes
(c)
Zero-mean Gaussian
(d)
Zero-mean Student-t


## 7. Conclusion

In this paper, we propose a new non-log-concave sparsity prior, referred to as Double Lomax Prior, which corresponds to a three-level hierarchical Bayes model. A VB inference procedure is introduced to solve for the SLM using Double Lomax priors. When noninformative hyperprior is assumed (f→+0), the proposed update procedure includes the canonical ARD update procedure as a special case. However, comparing to the canonical ARD update, the proposed update procedure has one more latent variable which has the smoothness effect that can result in an improved performance. A hybrid formulation of Double Lomax and ARD is also considered as a tradeoff between performance and computational speed. Experiments on both correlated and uncorrelated SLM simulations with applications to AR model identification and compressive sensing are carried out to demonstrate the effectiveness of the proposed approach.



---

*Source: 176249-2013-08-27.xml*