A MOREAU-YOSIDA APPROXIMATION SCHEME FOR A CLASS
OF HIGH-DIMENSIONAL POSTERIOR DISTRIBUTIONS

arXiv:1505.07072v2 [] 26 Jun 2016

YVES F. ATCHADÉ

(Feb. 2016, first version May 2015)
Abstract. Exact-sparsity inducing prior distributions in Bayesian analysis typically lead to posterior distributions that are very challenging to handle by standard
Markov Chain Monte Carlo (MCMC) methods, particular in high-dimensional models with large number of parameters. We propose a methodology to derive smooth
approximations of such posterior distributions that are, in some cases, easier to
handle by standard MCMC methods. The approximation is obtained from the
forward-backward approximation of the Moreau-Yosida regularization of the neg√
ative log-density. We show that the derived approximation is within O( γ) of
the true posterior distribution in the β-metric, where γ > 0 is a user-controlled
parameter that defines the approximation. We illustrate the method with a highdimensional linear regression model.

1. Introduction
Successful handling of statistical models with large number of parameters from
limited data hinges on the ability to solve efficiently and simultaneously two problems: (a) weeding out non-significant variables, and (b) estimating the effect of the
significant variables. The concept of sparsity has come to play a fundamental role
in this endeavor. In the Bayesian framework, sparsity is naturally built in the prior
distribution using spike-and-slab priors (Mitchell and Beauchamp (1988); George and
McCulloch (1997)), which are mixtures of a point mass at the origin (the spike) and a
continuous density (the slab). We will refer to such priors as exact-sparsity inducing
priors. A number of recent works have established that these priors, with carefully
chosen slab densities, produce posterior distributions with optimal posterior contraction rates (Castillo et al. (2015); Atchadé (2015)). However, the flip side of such stellar
2000 Mathematics Subject Classification. 60F15, 60G42.
Key words and phrases. Moreau-Yosida approximation, Markov Chain Monte Carlo, Spike-andslab priors, Variable selection, High-dimensional models.
This work is partially supported by the NSF grant DMS 1228164 and 1513040.
Y. F. Atchadé: University of Michigan, 1085 South University, Ann Arbor, 48109, MI, United
States. E-mail address: yvesa@umich.edu.
1

2

YVES F. ATCHADÉ

statistical properties is the fact that these posterior distributions are computationally
difficult to handle, particularly in high-dimensional applications. Deriving tractable
and scalable approximations for such distributions is therefore a problem of practical
importance.
The most commonly used approach for dealing with posterior distributions from
exact-sparsity inducing priors consists in integrating out the regression coefficients
(George and McCulloch (1997); Bottolo and Richardson (2010)). Recent results by
Yang et al. (2015) have shown that such an approach indeed scales well with the
dimension of the parameter space. However, it holds the limitation that it does not
solve both the variable selection and sparse estimation problems jointly. Furthermore,
it does not easily extend to non-Gaussian slabs1, or to non-Gaussian models. Another
solution to dealing with these posterior distributions is to design specialized MCMC
samplers, typically using trans-dimensional MCMC techniques such as reversible jump
MCMC (Chen et al. (2011)), or the newly proposed shrinkage-thresholding Metropolis
adjusted Langevin algorithm (STMaLa; Schreck et al. (2013)). See also Ge et al.
(2011) for a specialized sampler for blind-deconvolution models. However, these transdimensional MCMC samplers currently do not scale well to large problems (Schreck
et al. (2013)).
The discussion above suggests that when dealing with exact-sparsity inducing priors in high-dimensional regression problems, scalable approximation of the posterior
distribution would be useful. Notice that the Laplace approximation (Tierney and
Kadane (1986)), one of the most standard approximation tool in Bayesian computation, cannot be straightforwardly applied when the dimension of the space is as big
as the sample size (Shun and McCullagh (1995)). Variational Bayes approximations
recently explored by Ormerod et al. (2014) form a promising solution, but remained
to be fully explored in high-dimensional settings.
1.1. Main contribution. We propose a methodology to approximate posterior distributions derived from exact-sparsity inducing priors. An interesting feature of the
approximation is that the approximation error is easily controlled by the user. Furthermore, in several important cases, the approximation thus obtained is easily explored by standard Markov Chain Monte Carlo (MCMC) algorithms. The approximation is obtained by taking the forward-backward approximation (closely related to
the Moreau-Yosida approximation) of the negative log-density. The Moreau-Yosida
approximation is a well-established regularization method in optimization for dealing
with non-smooth and constrained problems (Moreau (1965); Bauschke and Combettes
1

For optimal posterior contraction rate, currently available results suggest that one needs slab
densities with tails heavier than Gaussian.

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

3

(2011)). Several recent works have recognized the usefulness of the Moreau-Yosida
regularization for Bayesian computation. Pereyra (2015) noted that a log-concave
density can be well approximated by its Moreau-Yosida approximation. However, the
framework developed by Pereyra (2015) does not handle the class of posterior distributions considered here. Another related work is the STMaLa of Schreck et al.
(2013) mentioned above, which implicitly uses the Moreau approximation to design Metropolis-Hastings proposals to sample from posterior distributions with exactsparsity inducing prior distributions.
If Π̌(·|z) denotes the posterior distribution of interest on Rd × {0, 1}d given data z,
we write Π̌γ (·|z) to denote the proposed Moreau-Yosida approximation, where γ > 0
is a user-controlled parameter that defines the quality of the approximation. We
derive a general result (Theorem 9) that shows, under some regularity conditions,
that

√
(1)
dβ Π̌(·|z), Π̌γ (·|z) = O( γ),
where dβ is the β-metric that metricizes weak convergence (see Section 1.2 for precise
definition). One challenge with using the proposed approximation is to find values of
γ for which Π̌γ is close to Π̌, but not too close so that Markov Chain Monte Carlo
samplers with good mixing properties can be easily developed for Π̌γ . In Theorem 8 we
propose a choice of γ that strikes the aforementioned balance, as we show empirically
in the simulation examples. Furthermore, with this particular choice of γ, we show
that the constant in the big O in (1) degrades with the dimension d at most linearly.
We illustrate the method using a linear regression model with a spike-and-slab
prior, where the slab is the elastic net density (Li and Lin (2010)). The example
has relevance because the posterior distribution thus defined (actually a special case
thereof) is known to contract at the optimal rate (Castillo et al. (2015)). Our proposed
methodology produces an approximation Π̌γ of this posterior distribution, and we
develop an efficient Markov Chain Monte Carlo algorithm to sample from Π̌γ . In this
particular example, we show that the Moreau-Yosida approximation actually scales
very well with the dimension. More precisely, we derive an upper bound similar to
(1) that degrades at most like log(d) as d → ∞ (see Theorem 10). We illustrate
these results in a simulation study which shows that the method performs well, and
outperforms STMaLa for high-dimensional problems. A Matlab implementation can
be obtained from http://dept.stat.lsa.umich.edu/∼ yvesa/Research.html.
The remainder of the paper is organized as follows. We close the introduction with
some notation that will be used throughout the paper. In Section 2, we first introduce
the class of posterior distributions of interest, followed in Section 3 by the basic idea
of the Moreau-Yosida approximation. In Section 4, we develop how the idea can be

4

YVES F. ATCHADÉ

applied to approximate the posterior distributions of interest. Section 5 details an
application to linear regression models. We close the paper with further discussion in
Section 6. All the proofs are gathered in the Appendix, placed in a supplemental file.
1.2. Notation. Throughout the paper, d ≥ 1 is a given integer and Rd denotes the
d-dimensional Euclidean space equipped with its Borel sigma-algebra, its Euclidean
def Pd
norm k · k, and inner product h·, ·i. We also use the norms kθk1 =
j=1 |θj |, and
kθk0 defined as the number of non-zero components of θ. The Lebesgue measure on
Rd is written as dx when there is no confusion.
def
We set ∆ = {0, 1}d . For δ ∈ ∆, µδ denote the product measure on Rd defined
def Q
as µδ (dθ) = dj=1 νδj (dθj ), where ν0 (dz) is the Dirac mass at 0, and ν1 (dz) is the
Lebesgue measure on R. Hence integration with respect to µδ sets to zero all the
components for which δj = 0, and integrates the remaining components using the
standard Lebesgue measure.
For θ, ϑ ∈ Rd , θ · ϑ denotes the component-wise product: (θ · ϑ)j = θj ϑj , 1 ≤ j ≤ d.
For δ ∈ ∆, we shall write θδ to denote θ · δ, and we set
def

Rdδ = {θδ , θ ∈ Rd } = {θ ∈ Rd : θj = 0 for δj = 0, j = 1, . . . , d}.
We will need ways to evaluate the distance between two probability measures. Let
(X, dX ) be some arbitrary separable complete metric space equipped with its Borel
sigma-algebra. For any two probability measures µ1 , µ2 on X, the β-distance between
µ1 , µ2 is defined as
Z
Z
def
f (x)µ1 (dx) −
f (x)µ2 (dx) ,
dβ (µ1 , µ2 ) = sup
(2)
kf kBL ≤1

X

X

where the supremum is taken over all measurable functions f : X → R such that
def
kf kBL = kf k∞ + kf kL ≤ 1, where


|f (x1 ) − f (x2 )|
def
def
kf k∞ = sup |f (x)|, and kf kL = sup
, x1 , x2 ∈ X, x1 =
6 x2 .
dX (x1 , x2 )
x∈X
It is well-known that this metric metricizes weak convergence (see e.g. Dudley
(2002) Theorem 11.3.3). If the supremum in (2) is replaced by a supremum over
all measurable functions f : X → R such that kf k∞ ≤ 1 (resp. kf kL ≤ 1) one obtains
the total variation metric dtv (resp. the Wasserstein metric dw ).
2. High-dimensional posterior distributions with sparse priors
Let z be a realization of some random variable Z with conditional distribution fθ ,
given a parameter θ ∈ Rd . With a prior distribution Π on θ, the posterior distribution

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

5

for learning θ is
fθ (z)Π(dθ)
.
Rd fθ (z)Π(dθ)
Although the prior distribution Π can be constructed in a variety of ways, we focus
on exact-sparsity inducing priors (spike-and-slab priors). Such prior distributions
have been recently shown to produce posterior distributions with optimal contraction
properties (Castillo et al. (2015); Atchadé (2015)). More specifically, we consider a
prior distribution Π on ∆ × Rd of the form
Π̌(dθ|z) = R

Π(δ, dθ) = πδ Π(dθ|δ),
for a discrete distribution {πδ , δ ∈ ∆} on ∆, and a prior Π(·|δ) that is built as follows.
Given δ, the components of θ are independent, and for 1 ≤ j ≤ d,
(
Dirac(0) if δj = 0
,
(3)
θj |δ ∼
p(·)
if δj = 1
where Dirac(0) is the Dirac measure on R with full mass at 0, and p(·) is a positive
density on R. By the standard data-augmentation trick, we will take the variable δ
as part of the posterior distribution. As defined, the support of Π(·|δ) is Rdδ = {θ ∈
Rd : θj = 0 for δj = 0, 1 ≤ j ≤ d}, and Π(·|δ) has a density with respect to the
measure µδ defined in Section 1.2:
Π(dθ|δ) = e−P (θ|δ) µδ (dθ), where
( P
− j: δj =1 log p(θj ) if θ ∈ Rdδ
def
P (θ|δ) =
+∞
otherwise .
In the above formula, and throughout the paper, we convene that e−∞ = 0, and
0 × ∞ = 0. We also define
def

`(θ) = − log fθ (z),

def

and h(θ|δ) = `(θ) + P (θ|δ), θ ∈ Rd ,

so that the posterior distribution writes
Π̌(δ, dθ|z) ∝ πδ e−h(θ|δ) µδ (dθ).

(4)

Monte Carlo simulation from this posterior distribution can be challenging. The
issue is related to the discrete-continuous mixture form of the spike-and-slab prior on
θ, which has the effect that any two distributions Π̌(δ, ·|z) and Π̌(δ 0 , ·|z) are mutually
singular for δ 6= δ 0 . As a result, if direct sampling from the conditional distribution
of θ|δ, z is not possible, then sampling from (4) requires the use of trans-dimensional
MCMC methods such as reversible jump (Chen et al. (2011)), or STMaLa (Schreck
et al. (2013)) which is shown to perform better than reversible jump. However,
one issue with STMaLa is that the algorithm has several tuning parameters that

6

YVES F. ATCHADÉ

are currently poorly understood. Furthermore, as we shall see in the simulations,
the mixing of the algorithm degrades significantly for high-dimensional problems,
particularly when the signal is weak.
3. The Moreau-Yosida approximation
Our goal in this work is to develop a more tractable approximation to the posterior distribution Π̌ in (4) using the Moreau-Yosida approximation. However to make
the ideas easy to follow, we start with some general discussion of the Moreau-Yosida
approximation. Let h : Rd → (−∞, +∞] be a convex, lower semi-continuous function that is not identically +∞, and let µ be a sigma-finite measure on Rd . In the
applications, µ will naturally be taken as the Lebesgue measure on the domain of h
(the domain of h is the set of points x ∈ Rd such that h(x) < ∞). Assuming that
def R
Z = Rd e−h(x) µ(dx) < ∞, we consider the probability measure
1 −h(x)
e
µ(dx).
(5)
Z
To fix the ideas, the reader may think of the case where h is finite everywhere and µ
is the Lebesgue measure on Rd . In that case ν is the probability distribution on Rd
with density (1/Z)e−h(x) . However our main interest is in the posterior distribution
(4) for which the slightly more general setting is needed.
Suppose that we are interested in drawing samples from ν. The lack of smoothness
of h, and the possibly complicated geometry of the support of ν can create difficulties
for standard MCMC algorithms. A smooth approximation of ν can be formed from
the Moreau-Yosida approximation of h defined for γ > 0 as


1
2
h̃γ (x) = min h(u) +
ku − xk , x ∈ Rd .
2γ
u∈Rd
ν(dx) =

Under the assumptions imposed on h above, the function h̃γ is known to be welldefined and finite everywhere. It is also convex, continuously differentiable with a
Lipschitz gradient, and h̃γ (x) ↑ h(x), as γ ↓ 0, for all x ∈ Rd . All these properties
are well-known and can be found in Bauschke and Combettes (2011) (Chapter 12).
def R
Assuming that Z̃γ = Rd e−h̃γ (x) dx < ∞, it seems natural to consider the probability
measure
1 −h̃γ (x)
ν̃γ (dx) =
e
dx,
Z̃γ
as an approximation of ν. To the best of our knowledge, the approximation ν̃γ was
first considered by Pereyra (2015), for a probability distribution ν for which h is finite
everywhere and µ is the Lebesgue measure on Rd . And we refer the reader to that
paper for a good discussion of the basic properties of ν̃γ , and how well it approximates

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

7

ν. In particular Pereyra (2015) showed that the smoothness of h̃γ can be exploited to
derive efficient gradient-based MCMC samplers for ν. An important limitation of the
Moreau-Yosida approximation is that it is typically not available in closed form, and
its computation leads to a d-dimensional, possibly complicated optimization problem.
In many problems the function h takes the particular form
h(x) = `(x) + P (x), x ∈ Rd
where ` is convex, finite everywhere and twice continuously differentiable, and P is
convex, not identically +∞ and lower semi-continuous. In such cases, one can approximate ` around a given point x by its linear approximation u 7→ `(x) + h∇`(x), u − xi,
where ∇`(x) denote the gradient of ` at x. This approximation leads to the so-called
forward-backward approximation of h, defined for γ > 0 as


1
def
2
hγ (x) = min `(x) + h∇`(x), u − xi + P (u) +
ku − xk , x ∈ Rd
2γ
u∈Rd


γ
1
2
2
= `(x) + − k∇`(x)k + min P (u) +
ku − x + γ∇`(x)k . (6)
2
2γ
u∈Rd
Under the assumptions imposed on ` and P above, the function hγ is finite everywhere,
continuously differentiable, and hγ ≤ h. These properties can be found in Patrinos
et al. (2014) Theorem 2.2, but are easy to derive. For instance, the differentiability
follows from the expression (6), the twice differentiability of `, and the differentiability
of the Moreau-Yosida approximation of P . Notice however that hγ is no longer convex
def R
in general. Assuming that Zγ = Rd e−hγ (x) dx < ∞ it seems also natural to consider
the resulting approximation of ν defined as
νγ (dx) =

1 −hγ (x)
e
dx.
Zγ

The main advantage of hγ over h̃γ is that in many problems of interest hγ is available
in closed form, whereas h̃γ is not. Furthermore, if the function P is separable, then
the computation of hγ leads to d separate one-dimensional optimization problems.
However, the price to pay for the computational convenience is that hγ may not be
convex, and it is a less accurate approximation of h. Indeed, by the convexity of `,
we have `(u) ≥ `(x) + h∇`(x), u − xi for all u ∈ Rd . Hence hγ (x) ≤ h̃γ (x) ≤ h(x) for
all x ∈ Rd . As we will see, the convergence hγ (x) ↑ h(x), as γ ↓ 0, for all x ∈ Rd ,
still holds. Figure 1 gives an illustrative example of the differences between hγ and
h̃γ and how both functions approximate h.
Since hγ converges pointwise to h as γ ↓ 0, it seems natural to expect that νγ
approaches ν for small γ. If the function h is finite everywhere, one can easily show
(see Proposition 1 below) that indeed, νγ converges to ν in the total variation metric,

8

YVES F. ATCHADÉ
γ=5

γ=1

γ = 0.1

2

2

2

1.8

1.8

1.8

1.6

1.6

1.6

1.4

1.4

1.4

1.2

1.2

1.2

1

1

1

0.8

0.8

0.8

−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

1

Figure 1. Figure showing the function h(x) = −ax + log(1 + eax ) + b|x|
for a = 0.8, b = 0.5 (blue/solid line), and the approximations hγ and h̃γ
(hγ ≤ h̃γ ), for γ ∈ {5, 1, 0.1}. For γ = 0.1, the curves of hγ and h̃γ are
almost undistinguishable on the figure.

as γ ↓ 0. However this result is no longer true when the domain of h has zero Rd Lebesgue measure. In this latter case, we will show that the convergence of νγ occurs
only weakly, or in the Wasserstein metric.
Proposition 1. Suppose µ in (5) is the Lebesgue measure on Rd , h = ` + P is
convex, finite everywhere, and hγ (x) ↑ h(x) for all x ∈ Rd . Suppose also that there
exists γ0 > 0 such that Zγ0 < ∞. Then for all γ ∈ (0, γ0 ], νγ is well-defined, and


Z
dtv (νγ , ν) ≤ 2 1 −
↓ 0, as γ ↓ 0.
Zγ
Proof. See the Appendix.



Remark 2.
(1) Notice that Proposition 1 can also be applied to ν̃γ by taking
` ≡ 0.
(2) We show in Lemma 2 in the Appendix that if ` is finite everywhere and differentiable, and P is finite everywhere, convex with a nonempty subdifferential
at x for all x ∈ Rd , then hγ ↑ h, as required in the proposition.
If the domain of h has Rd -Lebesgue measure 0, then ν and νγ are then automatically mutually singular and Proposition 1 cannot hold. The following toy example
illustrates this case.
Example 3. Suppose that we take Rd = R, ` ≡ 0, and we take P (x) = 0 is x =
and P (x) = +∞ if x 6= 0. In that case e−h(x) = 1 if x = 0, and e−h(x) = 0 if x 6=
Let µ = δ0 be the point mass probability measure at 0. Hence ν = δ0 . For γ >
hγ (x) = h̃γ (x) = x2 /(2γ), x ∈ R. Hence νγ is the normal distribution N(0, γ).

0,
0.
0,
It

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

9

follows that dtv (νγ , ν) = 2, for all γ > 0. But for any Lipschitz function f : R → R
with Lipschitz constant 1,
r
2γ
|νγ (f ) − ν(f )| = |νγ (f ) − f (0)| ≤ E(|Zγ |) =
,
π
q
where Zγ ∼ N(0, γ). By taking f = | · |, it can be easily seen that dw (νγ , ν) = 2γ
π .
Hence νγ converges in the Wasserstein metric to ν, but not in total variation. And
√
the convergence rate is O( γ).
Remark 4. The fact that we only have convergence in the Wasserstein metric has
practical implications. It implies that one needs to be cautious about what probability
ν(A) can be well approximated by νγ (A). For instance, in Example 3, if A is of the
form (a, 0] or [0, b], then ν(A) = 1, whereas limγ↓0 νγ (A) = 0.
In the next section we will use the approximating measure νγ introduced above to
approximate the posterior distribution (4). We will see that the situation is similar
to the one in Example 3, and as in that example we will show that the approximation
converges weakly to the posterior distribution Π̌.
4. The Moreau-Yosida approximation of the posterior distribution (4)
In this section, we return to the posterior distribution (4) defined in Section 2. And
we make the following assumptions on the functions ` and P .
(1) The function θ 7→ `(θ) is finite everywhere, convex, and twice continuously differentiable.
(2) For all δ ∈ ∆, the function θ 7→ P (θ|δ) is convex, lower semi-continuous, not
identically +∞, and admits a sub-gradient g(θ|δ) at θ, for all θ ∈ Rdδ .

H 1.

Remark 5.
(1) The convexity assumption on ` is fundamental and delineates
the type of problems to which the proposed approximation could be easily
applied. Extension beyond this set up is possible, but will require fundamentally different techniques.
(2) The convexity of P (·|δ) boils down to the log-concavity of the density p in the
prior (3). Most of the sparsity promoting prior densities used in practice are
log-concave.
Given δ ∈ ∆, we consider the forward-backward approximation of h(·|δ) defined as


1
def
hγ (θ|δ) = min `(θ) + h∇`(θ), u − θi + P (u|δ) +
ku − θk2 , θ ∈ Rd ,
(7)
2γ
u∈Rd

10

YVES F. ATCHADÉ

for some parameter γ > 0. Using hγ , we propose to approximate the posterior
distribution Π̌ in (4) by
Π̌γ (δ, dθ|z) ∝ πδ (2πγ)

kδk0
2

e−hγ (θ|δ) dθ,

(8)

that we call the Moreau-Yosida approximation of Π̌, although as we have seen above,
(7) is only the forward-backward approximation of h(·|δ). In the expression (8),
π denotes the irrational number. The function hγ (·|δ) is available in closed form
whenever the Moreau-Yosida approximation of P (·|δ) has a closed form expression.
More specifically, for δ ∈ ∆, and for γ > 0, we define


1
def
2
ku − θk , θ ∈ Rd
(9)
Pγ (θ|δ) = min P (u|δ) +
2γ
u∈Rd
the Moreau-Yosida approximation of P , and its associated proximal map


1
def
2
Proxγ (θ|δ) = Argmin u∈Rd P (u|δ) +
ku − θk , θ ∈ Rd .
2γ
From the definition of Pγ and Proxγ , we see that hγ can be alternatively written as
γ
hγ (θ|δ) = `(θ) − k∇`(θ)k2 + Pγ (θ − γ∇`(θ)|δ)
(10)
2
= `(θ) + h∇`(θ), Jγ (θ|δ) − θi + P (Jγ (θ|δ)|δ)
1
(11)
+ kJγ (θ|δ) − θk2 ,
2γ
where
def

Jγ (θ| δ) = Proxγ (θ − γ∇`(θ)|δ) .
For γ > 0, θ ∈ Rd , let sγ (θ) ∈ Rd be such that


1
def
2
(sγ (θ))j = Argmin u∈R − log p(u) +
(u − θj ) , 1 ≤ j ≤ d.
2γ
Then it is easy to check that Proxγ (θ|δ) = δ · sγ (θ). Hence by Equation (11), we
see that hγ (·|δ) is computationally tractable if the map sγ (the proximal map of the
negative log-prior) is easy to compute. Although this limits the applicability of the
method, there a several priors commonly used for which this holds, including the
Laplace prior and the elastic-net prior given respectively by


u2
−λ|u|
p(u) ∝ e
, and p(u) ∝ exp −αλ1 |u| − (1 − α)λ2
,
2
as well as the generalized double Pareto of Armagan et al. (2013), and the (improper)
prior distribution that arises from the MCP of Zhang (2010), given respectively by
!



Z |u| 
1
|u| −(α+1)
t
p(u) =
1+
, and p(u) = exp −λ
1−
dt .
2λ
αλ
αλ +
0

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

11

4.1. Connection with spike-and-slab priors. The proposed approximation Π̌γ is
closely related to the distribution obtained by replacing all the Dirac mass in (3) by
independent Gaussian distributions N(0, γ), γ > 0. More precisely, let Π̃γ denote the
posterior distribution of (δ, θ) in the following model.
(
N(0, γ) if δj = 0
δ ∼ {πδ }, θj |δ ∼
, 1 ≤ j ≤ d, and Z|δ, θ ∼ fθδ .
(12)
p(·)
if δj = 1
Notice that in (12) given (δ, θ), we draw Z from fθδ , with a sparse parameter θδ . The
posterior distribution thus defined is

Π̃γ (δ, dθ|z) ∝ πδ

1
2πγ

 d−kδk1
2

e

1
kθ−θδ k2 −h(θδ |δ)
− 2γ

e

dθ.

(13)

The distribution Π̃γ in turn, is closely related to another posterior distribution commonly used in practice and obtained from the following model:
(
N(0, γ) if δj = 0
δ ∼ {πδ }, θj |δ ∼
, 1 ≤ j ≤ d, and Z|δ, θ ∼ fθ ,
(14)
p(·)
if δj = 1
for some constant γ > 0. Here given (δ, θ), we draw Z from fθ . Model (14) is
widely used in practice as a more tractable alternative to the point-mass spike-andslab (George and McCulloch (1997); Ishwaran and Rao (2005); Rockova and George
(2014); Narisetty and He (2014)). Clearly, to the extend that the point-mass spikeand-slab is the gold-standard, the model in (12) is preferable to the one in (14).
The Moreau-Yosida approximation proposed in this paper can be viewed as a very
close approximation to Π̃γ , as we show that dtv (Π̌γ , Π̃γ ) = O(γ) (see Lemma 3 in
the Appendix, and (20)), where dtv denotes the total variation metric. The interest
of our method then comes from the fact that sampling from Π̌γ is much easier than
sampling from Π̃γ . Indeed, notice that in Π̃γ , the parameter δ appears also in the
likelihood function `(θδ ). As a result, both conditional distributions δ|z, θ and θ|z, δ
are typically intractable and require MCMC algorithms. Whereas in Π̌γ , given (z, θ),
the components of δ are independent Bernoulli random variables.
4.2. Approximation bounds. We will now derive a result that bounds the βdistance between Π̌γ and Π̌. We define
Z
def
%γ (z) = log erγ (δ,θ) Π̃γ (dδ, dθ|z),
(15)
where
def

rγ (δ, θ) = h∇`(θ) − ∇`(θδ ), θ − θδ )i +

γ
kδ · ∇`(θ) + δ · g (θδ |δ) k2 .
2

12

YVES F. ATCHADÉ

For simplicity, we shall omit the dependence of rγ (δ, θ) on z (same with `(θ) and
∇`(θ)). We note that by the convexity of `, rγ (δ, θ) ≥ 0. Hence %γ (z) ≥ 0.
Theorem 6. Assume H1, for some fixed data z. Suppose that there exists γ0 > 0
such that Π̌γ0 (·|z) is well-defined. Then for all γ ∈ (0, γ0 ], Π̌γ (·|z) is well-defined and


 p
dβ Π̌γ (·|z), Π̌(·|z) ≤ γd + 2 1 − e−%γ (z) .
(16)
Proof. See the Appendix.



Notice that 1 − e−x ≤ x for all x ≥ 0. Therefore, the convergence to zero of

dβ Π̌γ (·|z), Π̌(·|z) would follow if the term %γ (z) converges to 0 as γ → 0. In
the next result we impose some additional assumptions, which, together with H1
guarantee that %γ (z) converges to zero. In the process we derive an explicit bound
on the convergence rate which can be used to develop guidelines for choosing γ. We
make the following assumption.
H 2.

(1) There exists L1 < ∞ such that,
k∇`(θ1 ) − ∇`(θ2 )k ≤ L1 kθ1 − θ2 k, θ1 , θ2 ∈ Rd .

(17)

(2) There exists L2 < ∞ such that
kδ · ∇`(θ)k2 ≤ 2L2 `(θ), δ ∈ ∆, θ ∈ Rdδ .

(18)

(3) For all δ ∈ ∆, there exists c(δ) < ∞, such that
kδ · g(θ|δ)k2 ≤ c(δ) + 2L2 P (θ|δ), θ ∈ Rdδ .

(19)

Remark 7. H2-(1) is a standard Lipschitz assumption. H2-(2) and H2-(3) essentially
requires both functions ` and P (·|δ) to grow like O(kθk2 ), or o(kθk2 ), as kθk → ∞.
Theorem 8. Assume H1-H2, for some fixed data z, and suppose γ > 0 is such that
4γ max(L1 , L2 ) ≤ 1. Then Π̌γ (·|z) is well-defined and


 p
dβ Π̌γ (·|z), Π̌(·|z) ≤ γd + 2 1 − e−%γ (z) ,
where



1
%γ (z) ≤ 3γ
max c(δ) + d(L1 + 2L2 ) + L2 R(z) ,
2 δ∈∆

(20)

def

where R(z) = maxδ∈∆ inf θ∈Rd [`(θδ ) + P (θδ |δ)] ≤ maxδ [`(0) + P (0|δ)].
Proof. See the Appendix.



ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

13

Since 1 − e−x ≤ x, for all x ≥ 0, Theorem 8 shows that as γ → 0, the Moreau√
Yosida approximation Π̌γ (·|z) approaches Π̌(·|z) at the rate of γ, under H1 and H2.
The rate is optimal as Example 3 shows. Theorem 8 also provides some guidelines
for choosing γ, as it suggests that one can choose γ as
γ=

γ0
,
max(L1 , L2 )

1
with 0 < γ0 ≤ .
4

(21)

The bound in (20) seems to suggest that the quality of the approximation resulting
from choosing γ as in (21) degrades only linearly with the dimension d, as d increases2.
However, it is important to realize that the bound in (20) is most likely not tight, and
the dependence of %γ (z) on d could be even better than O(d) (see the linear regression
example below). In general we cautious against the use of too small values of γ, since
choosing Π̌γ very close to Π̌ limits the ability to construct good MCMC sampler to
explore Π̌γ .
5. Application to Bayesian linear regression with sparse priors
As an application we consider a high-dimensional linear regression problem, with
dependent variable z ∈ Rn , and design matrix X ∈ Rn×d . The variance term σ 2 is
assumed known. The negative-log-likelihood function ` for this problem can be taken
as
1
`(θ) = 2 kz − Xθk2 , θ ∈ Rd .
2σ
We will set up the prior distribution of θ using δ ∈ ∆, and using an auxiliary variable
def
φ = (q, λ1 , λ2 ), where q ∈ (0, 1) is a sparsity parameter, and λ1 > 0, λ2 > 0 are regularization parameters. Given φ, we assume that the components of δ are independent
and identically distributed, with distribution Ber(q). Hence πδ = qkδk0 (1 − q)d−kδk0 .
Given φ and δ, the components of θ are independent, and for 1 ≤ j ≤ d,
(
Dirac(0)

 if δj = 0 ,
θj |δ, φ ∼
λ1 λ2
EN σ2 , σ2
if δj = 1
where Dirac(0) is the Dirac measure on R with full mass at 0, and EN(λ1 /σ 2 , λ2 /σ 2 )
is the (elastic-net) distribution with density given by


1
λ1
λ2 2
exp −α 2 |x| − (1 − α) 2 x , x ∈ R,
(22)
Z(φ)
σ
2σ
2Indeed, the term R(z) always satisfies R(z) ≤ max [`(0) + P (0|δ)], and this latter expression
δ

typically does not grow with d.

14

YVES F. ATCHADÉ

for a parameter α ∈ [0, 1], assumed known. The normalizing constant Z(φ) can be
written as


 q
 σ
αλ
2π
1
√
if α ∈ [0, 1)
(1−α)λ2 erfcx σ 2(1−α)λ2
Z(φ) =
,
 2σ2
if
α
=
1
λ1
where erfcx(x) is the scaled complementary error function, which can be written as
√
2
erfcx(x) = 2ex Φ(− 2x), where Φ is the cdf of standard normal distribution. The
prior density (22) is a reparametrization of the elastic-net (Zou and Hastie (2005))
prior used by Li and Lin (2010). Notice that α = 1 makes λ2 inactive, and setting
α = 0 makes λ1 inactive.
Given the prior specified above, the function P becomes
P (θ|δ) = kδk1 log Z(φ) +

(1 − α)λ2
αλ1
kθδ k1 +
kθδ k2 + ιRd (θ), θ ∈ Rd ,
2
δ
σ
2σ 2

where ιRd (θ) = 0 if θ ∈ Rdδ , and ιRd (θ) = +∞ otherwise. We recall that Rdδ = {θ ∈
δ

δ

Rd : θj = 0, if δj = 0, j = 1, . . . , d}. With h(θ|δ) = `(θ) + P (θ|δ), the posterior
distribution of (δ, θ) is
Π̌ (δ, dθ|z) ∝ πδ e−h(θ|δ) µδ (dθ).
With the elastic net prior (22), the proximal function Proxγ (θ|δ) is easy to compute.
For x ∈ R, define sign(x) as 1 is x > 0, −1 if x < 0 and 0 if x = 0. For γ > 0, let
sγ (θ) ∈ Rd denotes the vector whose j-th component is given by


sign(θj ) |θj | − αγ σλ12
+
(sγ (θ))j =
.
(23)
λ2
1 + γ σ2 (1 − α)
It is easy to show that
Proxγ (θ|δ) = δ · sγ (θ),
where θ1 · θ2 denotes the component-wise product. From (8), it follows that the
Moreau-Yosida approximation Π̌γ of Π̌ has a density π̌γ given by
π̌γ (δ, θ|z) ∝ πδ (2πγ)

kδk0
2

e−hγ (θ|δ) ,

(24)

where hγ (·|δ) is given by (11). In the next result, we show that H1 and H2 hold for
this problem, and Theorem 8 applies. For a matrix A, let λmax (A) denote its largest
eigenvalue.
Corollary 9. Suppose that (1 − α)λ2 ≤ λmax (X 0 X), and suppose that γ > 0 satisfies
4γ
λmax (X 0 X) ≤ 1.
σ2

(25)

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

15

Then for all z ∈ Rn , Π̌γ (·|z) is a well-defined probability measure on ∆ × Rd , and


 p
dβ Π̌γ (·|z), Π̌(·|z) ≤ γd + 2 1 − e−%γ (z) ,
where %γ (z) satisfies
3γ
%γ (z) ≤
2



αλ1
σ2

2



3γ
kzk2
0
d + 2 λmax (X X) 3d +
.
σ
2σ 2

Proof. See the Appendix.

(26)


As in the general case above, (25) suggests choosing
γ=

γ0 σ 2
,
λmax (X 0 X)

γ0 ∈ (0, 1/4].

(27)

And with this choice, the bound in (26) deteriorates only linearly in d. In fact, (26)
is a worst case analysis and better bounds can be derived if one takes into account
the sampling distribution of the data z. We prove one such result below.
We shall take the frequentist viewpoint and assume that the observed data z is a
realization of Z where
Z = Xθ? + ,

where

 ∼ N(0, σ 2 In ),

(28)

for a sparse unknown vector θ? ∈ Rd . In the Bayesian as in the frequentist framework,
the recovery of θ? when d > n depends on the positiveness of some restricted and
sparse eigenvalues of X 0 X. We define these quantities next. Let δ? ∈ ∆, denote the
def
sparsity structure of θ? . That is, δ?,j = 1 if and only if |θ?,j | > 0. We set s? = kθ? k0 ,
the number of non-zero components of θ? . We define
 0 0

θ (X X)θ
def
κ = inf
: θ 6= 0, kθ − θδ? k1 ≤ 7kθδ? k1 ,
nkθk2
and s ∈ {1, . . . , d}, we define


def

κ̄(s) = sup


θ0 (X 0 X)θ
: 1 ≤ kθk0 ≤ s .
nkθk2

Finally, we define
def

E =





n

z ∈ R : max | hXk , z − Xθ? i | ≤ λ1 /2 .
1≤k≤d

Theorem 10. Assume (28), with a design matrix X that satisfies κ > 0. Choose
p
α = 1, λ1 = 4σ nκ̄(1) log(d), and πδ = qkδk0 (1 − q)d−kδk0 , where q = d−u , for some
constant u > 1. Suppose that γ > 0 is small enough so that
4γ
λmax (X 0 X) ≤ 1,
σ2

and

24nκ̄(1)γ ≤ σ 2 (u − 1).

16

YVES F. ATCHADÉ

Then if d > 2, P(Z ∈ E) ≥ 1 − 2/d, and for Z ∈ E,




2
2
κ̄(s? )
E [%γ (Z)|Z ∈ E] ≤ − log 1 −
+ u−1 + 4 log(4) + s? log 1 +
d
d
16 log(d)



3γ
128 96γ
+ us? log(d) +
+ 2 n κ̄(1) log(d) + 2 nλmax (X 0 X) + Tr(X 0 X) . (29)
κ
σ
σ
Proof. See the Appendix in the Supplement.



Remark 11. To give some context, the theorem considers the posterior distribution
p
Π̌(·|Z), with α = 1 (Laplace prior), and λ1 = 4σ nκ̄(1) log(d) in (22), and with
πδ = qkδk0 (1 − q)d−kδk0 , where q = d−u , for some constant u > 1. It has recently
been shown by Castillo et al. (2015) that with the above choices, the θ-marginal
of this posterior distribution contracts to a point-mass at θ? at the optimal rate
p
O( s? log(d)/n). Theorem 10 gives a bound on the average error of the MoreauYosida approximation of this posterior distribution via

 p

E dβ Π̌γ (·|Z), Π̌(·|Z) |Z ∈ E ≤ γd + 2E [%γ (Z)|Z ∈ E] ,
for Z ∈ E.
The right-side of (29) does not converge to zero as γ → 0. However, it provides
some useful insights. We see that if we choose γ > 0 as in (27), then the rightside of (29) grows with d at most like log(d) + Tr(X 0 X)/λmax (X 0 X). From random
matrix theory it is known for several classes of random matrices that for fixed n,
Tr(X 0 X)/λmax (X 0 X) is typically O(1) as d → ∞. An inspection of the proof of
Theorem 10 suggests that the dependence of the right-side of (29) on the sample size
n can be improved. We leave this for possible future work.
5.1. Dealing with the hyper-parameter φ. We use a fully Bayesian for selecting
the hyper-parameter φ = (q, λ1 , λ2 ). We assume independent priors such that q ∼
Beta(1, du ) for some constant u > 1, λ1 ∼ U(a, M ), and λ2 ∼ U(a, M ) for some
small positive constant a (we use a = 10−5 in the simulations), and for a large positive
constant M such that (1 − α)M ≤ λmax (X 0 X). If γ > 0 is such that (25) holds then
the β-distance between the resulting posterior distribution and its Moreau-Yosida
approximation satisfies the same bound as in Theorem 9.
5.2. Markov Chain Monte Carlo. The density π̌γ in (24) is a “standard” density, and various MCMC schemes can be used to sample from it. We propose a
Metropolized-Gibbs strategy.

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

17

5.2.1. Updating δ. Given θ and φ, it is easy to see that hγ (θ|δ) depends on δj only
through the expression
"
#
αλ1 |dj | + 0.5(1 − α)λ2 d2j
d2j − 2θj dj
δj (∇`(θ))j dj + log Z(φ) +
,
+
σ2
2γ
where dj is the j-th component of sγ (θ − γ∇`(θ); λ1 /σ 2 , λ2 /σ 2 ). Hence, we update
jointly and independently the δj by setting δj = 1 with probability er /(1 + er ), where
r = log

q
1
+ log(2πγ)
1−q 2
"

#
αλ1 |dj | + 0.5(1 − α)λ2 d2j
d2j − 2θj dj
− (∇`(θ))j dj + log Z(φ) +
+
.
σ2
2γ

5.2.2. Updating θ. Given δ and φ, we update the components of θ using a mix of
an independence Metropolis sampler, and a Metropolis Adjusted Langevin algorithm
(MaLa). The MaLa strategy needs some motivation. Although its definition might
perhaps suggest otherwise, the function Pγ in (9) is actually differential (Bauschke
and Combettes (2011) Proposition 12.29) and for all θ, H ∈ Rd ,
∇θ Pγ (θ|δ) · H =

1
hθ − Proxγ (θ|δ), Hi .
γ

And since ` is twice continuously differentiable in this example, the expression (10)
shows that hγ is in fact differential and for all θ, H ∈ Rd ,
D
E
∇θ hγ (θ|δ) · H = h∇`(θ), Hi − γ ∇`(θ), ∇(2) `(θ) · H




1
+
(θ − γ∇`(θ) − Jγ (θ|δ, φ)) , Id − γ∇(2) `(θ) · H
γ
D


E
1
θ − Jγ (θ|δ, φ), Id − γ∇(2) `(θ) · H .
=
γ
To avoid dealing with second order derivatives, and since γ is typically small, we make
the approximation Id − γ∇(2) `(θ) ≈ Id , and therefore, we approximate ∇θ hγ (θ|δ) by
def

Gγ (θ|δ) =

1
c
def
(θ − Jγ (θ|δ)) , and Ḡγ (θ|δ) =
Gγ (θ|δ),
γ
c ∨ kGγ (θ|δ)k

(30)

for a positive constant c. The function Ḡγ is introduced for further stability, in
the spirit of the truncated Metropolis adjusted Langevin algorithm (see e.g. Atchadé
(2006)). Hence, given δ and φ, one can update the components of θ using a MetropolizedLangevin-type algorithm where the drift function is given by the corresponding components of Ḡγ . This algorithm is similar to the proximal MaLa of Pereyra (2015).
However, when δj = 0, the corresponding component of Gγ (θ|δ) is θj /γ and is
typically very large and not very informative (particularly for γ small). To deal with

18

YVES F. ATCHADÉ

this, we use the following strategy. We update jointly the components θj for which
δj = 1 using the MaLa algorithm outlined above. Then, we group together all the
components for which δj = 0 and we update them jointly using an independence
Metropolis sampler. The proposal density of the Independence Metropolis sampler
is built by approximating Jγ (θ|δ) by Proxγ (θ|δ). This approximation makes sense
because, for γ ≈ 0, Jγ (θ|δ) = Proxγ (θ − γ∇`(θ)|δ) ≈ Proxγ (θ|δ).
To explain the detail of the independence sampler, let θδ = (θj , j : s.t. δj = 1),
u = (θj , j : s.t. δj = 0), and let us represent θ by the pair (θδ , u). Let h̃γ (θδ , u|δ) be
the function obtained by replacing Jγ (θδ , u|δ) by Proxγ (θδ , u|δ) in the expression of
hγ (θδ , u|δ). Because, Proxγ (θδ , u|δ) does not actually depend on u, we have
h̃γ (θδ , u|δ) = `(θ) + h∇`(θ), Proxγ (θδ , u|δ) − θi +
=

1
k Proxγ (θδ , u|δ) − θk2 + const.
2γ

1
kz − Xδ θδ − Xδc uk2
2σ 2
1
− 2 hz − Xδ θδ − Xδc u, X(Proxγ (θδ , u|δ) − δ · θ) + Xδc ui
σ
1
+ kuk2 + const.
2γ

It is then easy to see that u 7→ e−h̃γ (θδ ,u|δ) is proportional to the density of the
Gaussian distribution
−1


γ
γ
def
,
N 2 ΣXδ0 c X (Proxγ (θ|δ) − δ · θ) , γΣ , where Σ = Ikδc k − 2 Xδ0 c Xδc
σ
σ
where δ c is the vector 1 − δ, and for any δ ∈ ∆, Xδ ∈ Rn×kδk denote the sub-matrix
of X obtained by selecting the columns for which δj = 1. Notice that under the
2
assumption γ ≤ 4λmaxσ(X 0 X) , the matrix Σ is always positive definite. The acceptance
probability of this independence sampler is



exp hγ (θδ , u0 |δ) − h̃γ (θδ , u0 |δ)

 .
min 1,
exp hγ (θδ , u|δ) − h̃γ (θδ , u|δ)
We found this independence sampler to be extremely efficient, with an acceptance
probability typically above 90%.
5.2.3. Updating φ = (q, λ1 , λ2 ). We update q ∼ Beta(kδk1 + 1, d + du − kδk1 ), and
we update (λ1 , λ2 ) jointly using a Random Walk Metropolis algorithm with Gaussian proposal. For improved mixing, we adaptively tune the scale parameter of the
proposal density.

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

19

5.3. Simulation results and comparison with STMaLa. We illustrate the method
with a simulated data example. All the computations in this example were done using
Matlab 7.14 on a 2.8 GHz Quad-Core Mac Pro with 24 GB of 1066 DDR3 Ram.
We set n = 200, p = 500 and we generate the design matrix X by simulating the
rows of X independently from a Gaussian distribution with correlation ρ|j−i| between
components i and j. We set ρ = 0.9. Using X, we general the outcome z = Xθ? + σ,
with σ = 1 that we assume known. We build θ? by randomly selecting 10 components
that we fill with draws from the uniform distribution U(v/2, 3v/2), where  = ±1
with probability 1/2, all other components being set to zero. We consider two cases
p
for v: v = 1 (that we refer to below as SCENARIO 1), and v = log(d)/n ≈ 0.18
(SCENARIO 2). SCENARIO 2 is obviously more challenging since the average strength
of the signal is at the limit of what is detectable.
We set γ = γ0 σ 2 /λmax (X 0 X) as prescribed by (27) with two choices of γ0 : γ0 = 0.25,
and γ0 = 0.01.
We compare these two samplers to the STMaLa sampler of Schreck et al. (2013).
The comparison is slightly tricky because STMaLa uses a different prior, namely a
Gaussian “slab” prior. However, we expect both posterior distribution on (δ, θ) to be
close, and we expect (δ? , θ? ) to be close to the center of both distributions. For the
STMaLa, we use the Matlab code provided online by the authors, with the default
setting. Unlike our approach, this sampler requires the true value of the sparsity
parameter q, which we provide. We also edit their code to return the summary
statistics presented below.
We evaluate the mixing of these samplers by computing the following two metrics
along the MCMC iterations: the relative error and the F -score (to evaluate structure
recovery), defined respectively as
E (k) =

kθ(k) − θ? k
,
kθ? k

and F (k) =

2 × SEN(k) PREC(k)
SEN(k) + PREC(k)

,

where
Pd
SEN(k) =

j=1 1{|θ(k) |>0} 1{|θ?,j |>0}
j

Pd

j=1 1{|θ?,j |>0}

Pd
, PREC(k) =

j=1 1{|θ(k) |>0} 1{|θ?,j |>0}
j

Pd

j=1 1{|θ(k) |>0}

. (31)

j

In stationarity we expect values of E (k) (resp. F (k) ) to be close to zero (resp. one).
In the absence of a better metric, we will graphically access the mixing time of the
samplers by looking at how quickly the sequence E (k) (resp. F (k) ) converges towards
zero (resp. one). In order to account for the computing time, and for better comparison, we plot these metrics, not as function of the iterations k, but as function of the

20

YVES F. ATCHADÉ

Structure Recovery

0.6

gamma_0=0.25
gamma_0=0.01
stMaLa

0.0

0.2

Recovery

0.8

1.0

gamma_0=0.25
gamma_0=0.01
stMaLa

0.4

0.0 0.2 0.4 0.6 0.8 1.0 1.2

Relative error

Relative Error

0

20

40
Time (in seconds)

60

80

0

20

40

60

80

Time (in seconds)

Figure 2. Relative error and structure recovery as function of time
in SCENARIO 1. Based on 30 MCMC replications. The curves are
sub-sampled to improve the readability of the figure.
computing time needed to reach iteration k. For further stability in the comparison,
we repeat all the samplers 30 times and average the two metrics and the computing
times over these 30 replications.
All the chains are initialized by setting all components of θ(0) (and δ (0) ) to zero.
We run the samplers for a number of iterations that depends on θ? . In SCENARIO
1, we run the newly proposed sampler for 10, 000, and we run STMaLa for 120, 000
iterations. In SCENARIO 2, we run our proposed sampler for 40, 000, and we run
STMaLa for 250, 000 iterations.
Figure 2 and 3 present the results. First, we observe that that γ0 = 0.25 mixes
significantly better than γ0 = 0.01. We notice also that Π̌γ approximates (θ? , δ? )
only slightly better when γ0 = 0.01 compared to γ0 = 0.25. Overall, we found that
γ0 ∈ (0.1, 0.25) produces a very good approximation.
We also look at the usual sample path mixing of the proposed sampler by plotting
the trace plot, histogram, and the autocorrelation plot from a single run of the sampler
(Figure 4). Here, we consider only SCENARIO 1, and we set γ0 = 0.25. We look at
(k)
the MCMC output {θj , k ≥ 0}, for one component j for which δj = 0, and for one
component j for which δj = 1. From this sample path perspective, the plots suggest
that the proposed MCMC sampler has a good mixing.
5.4. Empirical Bayes implementation and further experimentation. A limitation of the methodology is that σ 2 is assumed known, which is rarely the case
in practice. We explore by simulation an empirical Bayes solution whereby σ 2 is

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

1.0
0.8

gamma_0=0.25
gamma_0=0.01
stMaLa

0.0

0.2

0.4
0.0

0.2

0.6

Recovery

0.8

gamma_0=0.25
gamma_0=0.01
stMaLa

0.6

Relative error

Structure Recovery

0.4

1.0

Relative Error

0

50

100

150

21

200

250

300

0

50

100

Time (in seconds)

150

200

250

300

Time (in seconds)

Figure 3. Relative error and structure recovery as function of time
in SCENARIO 2. Based on 30 MCMC replications. The curves are
sub-sampled to improve the readability of the figure.
ACF, non-selected coef.

Hist., non-selected coef.

ACF

0.0
-0.2

0

-0.15

20

0.2

0.4

0.6

80
60

Density

40

0.05
-0.05

V1[Rg]

0.8

1.0

Trace plot, non-selected coef.

0

5000

10000

15000

-0.15

-0.10

-0.05

0.00

0.05

0.10

0

200

400

Index

800

1000

ACF, non-selected coef.

Hist., selected coef.

0.6
ACF

0.4

10
0

0.0

-2.60

0.2

5

-2.55

Density

-2.50

0.8

15

-2.45

1.0

Trace plot, selected coef.

V2[Rg]

600
Lag

0

5000

10000
Index

15000

-2.60

-2.55

-2.50

-2.45

0

200

400

600

800

Lag

Figure 4. Trace plot, histogram, and autocorrelation plot, from one
MCMC run, using γ0 = 0.25. Top row is for a component j for which
the true value of δj is 0. Bottom row, true value of δj is 1.

1000

22

YVES F. ATCHADÉ

estimated from data. Following Reid et al. (2013) we estimate σ 2 by
σ̂n2 =

n 
2
X
1
yi − xi β̂λn ,
n − ŝλn
i=1

where β̂λ is the lasso estimate at regularization level λ, and λn is selected by 10fold cross-validation, and where ŝλn is the number of non-zeros components of β̂λn .
In the cross-validation, we choose λn as the value of λ that minimizes the MSE.
This leads to the empirical Bayes Moreau-Yosida posterior approximation that we
denote Π̌γ (·|z, σ̂n2 ). We do a simulation study using a semi-real dataset to compare
the distributions Π̌γ (·|z, σ̂n2 ) and Π̌γ (·|z) (with the true value of σ 2 set to one). We
use the colon dataset (Buhlmann and Mandozzi (2014)) downloaded from
http://stat.ethz.ch/~dettling/bagboost.html. The data gives microarray gene
expression levels for 2, 000 genes for n = 62 patients in a colon cancer study. We
randomly select a subset of p = 1, 000 variables to form a design matrix X ∈ R62×1,000 .
Following Buhlmann and Mandozzi (2014), we normalize each column of X to have
mean zero and variance unity. We simulate a sparse signal vector θ? ∈ Rp with
s = 5 non-zeros components, and where the non-zeros components are drawn from
U(−v − 1, −v) ∪ (v, v + 1). We consider two scenarios: v = 1 and v = 3. Using X and
θ? , we generate z = Xθ? + σ, with σ = 1, and  ∼ N(0, In ).
We set γ as in (27) with γ0 = 0.25. We evaluate the samplers along the same
metrics E and F. We average the results over 30 replications3 of the samplers, where
each sampler is run for 50, 000 iterations. The result is presented on Table 1. We
notice that the recoery of θ? is poor in both cases when v = 1. When the signal
is strong (v = 3), the empirical Bayes posterior distribution performs well, but as
expected, under-performs the posterior distribution with known variance.
Weak signal (v = 1) Strong signal (v = 3)
EB
True σ
EB
True σ
Relative error (in %) 97.3
91.7
12.4
9.4
F -score ( in %)
14.5
25.1
79.6
88.5

PN
(k)
,
k=B+1 E
P
N
−1
(k)
and (N − B)
F
,
averaged
over
30
MCMC
replications,
each
k=B+1
MCMC run is 5 × 104 iterations.

Table 1. Table showing the posterior estimates (N − B)−1

3here only X and θ are kept fixed. For each replication, the dataset z is re-simulated, and σ 2 is
?
n

re-estimated.

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

23

6. Further Discussion
In this work we have developed and analyzed a smooth approximation to highdimensional posterior distribution using the Moreau-Yosida envelop. The methodology can be readily extended to other high-dimensional statistical models (linear
and generalized linear regression models, graphical models, sparse PCA, and others).
Several theoretical issues remain. We have discussed some of these issues above. One
important problem that we did not directly address concerns the mixing properties
of the proposed MCMC algorithms, and the trade-off inherent to the methodology
between good approximation properties of Π̌γ , and good mixing of gradient-based
MCMC simulation from Π̌γ . Another potentially interesting direction of research
is the idea of treating Π̌γ itself as a quasi-posterior distribution, and investigating
directly its posterior contraction properties.
7. APPENDIX: Proof of the main results
def

For convenience, we introduce the product space Θ̄ = ∆ × Rd that we implicitly
def p
equip with the metric dΘ̄ (θ̄1 , θ̄2 ) = kδ1 − δ2 k20 + kθ1 − θ2 k2 , θ̄j = (δj , θj ), j = 1, 2.
7.1. Proof of Proposition 1. For all x ∈ Rd , and γ ∈ (0, γ0 ], e−h(x) ≤ e−hγ (x) ≤
e−hγ0 (x) . Hence Z ≤ Zγ ≤ Zγ0 < ∞. Since µ is the Lebesgue measure on Rd , we shall
write it as dx. For any bounded measurable function f : Rd → R, we have
Z


1
|Πγ (f ) − Π(f )| ≤
f (x) e−hγ (x) − e−h(x) dx
Zγ Rd
Z
(Zγ − Z)
+
|f (x)|e−h(x) dx,
Zγ Z
d
R
Z 

2kf k∞
≤
e−hγ (x) − e−h(x) dx
Zγ
Rd


Z
= 2kf k∞ 1 −
.
Zγ
The fact that Zγ → Z as γ ↓ 0, follows from Lebesgue’s monotone convergence applied
to e−hγ0 − e−hγ .
7.2. Proof of Theorem 6. We work on the product space Θ̄ = ∆ × Rd introduced
above. Throughout the proof, we assume that z is fixed, and at times we write Π̌(·|z)
simply as Π̌. Same for Π̃γ (·|z) and Π̌γ (·|z).
We prove the theorem in two steps. First in Lemma 12, we bound the Wasserstein
distance between the distributions Π̃γ and Π̌ by showing that for all γ > 0,
p
dw (Π̃γ , Π̌) ≤ γd.
(32)

24

YVES F. ATCHADÉ

Then in Lemma 14 we bound the total variation distance between Π̌γ and Π̃γ by
showing that for all γ ∈ (0, γ0 ],


dtv (Π̃γ , Π̌γ ) ≤ 2 1 − e−%γ (z) .
(33)
It is clear from their definitions that both the Wasserstein metric and the total variation metric are upper bounds for the metrix β, and the Theorem 6 follows by combining (32) and (33). The proof of Lemma 14 relies on a comparison result between
the functions h and hγ established in Lemma 13 that is also of independent interest.
Lemma 12. Let Π̃γ be the probability measure defined in (13).
r
r


p
1
2p
γd 1 − E(kηk0 ) ≤ dw (Π̃γ , Π̌) ≤ dγ 1 −
π
d

For all γ > 0,
1
E(kηk0 ),
d

(34)

where η is a random variable on ∆ with distribution given by the δ-marginal of Π̌,
R
that is P(η = δ) ∝ πδ Rd e−h(θ|δ) µδ (dθ), δ ∈ ∆.
Proof. For δ ∈ ∆, we set
def

Z

e−h(θ|δ) µδ (dθ),

C(δ) =

and C =

Rd

X

πδ C(δ).

(35)

δ∈∆

Using this notation, we can write
Π̌(δ, dθ|z) =

πδ C(δ)
1 −h(θ|δ)
def
Π̌(dθ|δ, z), where Π̌(dθ|δ, z) =
e
µδ (dθ).
C
C(δ)

For γ > 0, we notice that the normalizing constant of Π̃γ is
C

def

=

X


πδ

δ

=

X
δ

1
2πγ

Z
πδ

 d−kδk1 Z
2

1
− 2γ
kθ−θδ k2 −h(θδ |δ)

e

e

dθ

Rd

e−h(θ|δ) µδ (dθ),

(36)

Rd

which is the same as the normalizing constant of the posterior distribution Π̌. Hence
we get the following factorization of Π̃γ ,
Π̃γ (δ, dθ|z) =

πδ C(δ)
Π̃(dθ|δ, z),
C

def

where Π̃(dθ|δ, z) =

1 − 2γ1 kθ−θδ k2 −h(θδ |δ)
e
e
dθ.
C(δ)

We build following coupling of Π̌ and Π̃γ . First we generate η ∈ ∆ from the distribution δ 7→ πδ C(δ)
C , and we generate ϑ̌|η ∼ Π̌(dθ|z, η). Hence clearly, (η, ϑ̌) ∼ Π̌.
Given (η, ϑ̌), we generate ϑ̃ as follows. If ηj = 1, we set ϑ̃j = ϑ̌j . Otherwise we
√
generate independently Zj ∼ N(0, 1), and set ϑ̃j = γZj . It is also easy to check
that (η, ϑ̃) ∼ Π̃γ .

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

25

For any Lipschitz function on Θ̄ with Lipschitz constant less of equal to 1, we have
Z
Z
h
i
f (δ, θ)Π̃γ (dδ, dθ) − f (δ, θ)Π̌(dδ, dθ) = E f (η, ϑ̃) − f (η, ϑ̌)
h
i p r
1
≤ E kϑ̃ − ϑ̌k ≤ dγ 1 − E(kηk0 ).
d
P
Now consider the function f0 (δ, θ) = √1d dj=1 |θj |. It is Lipschitz with Lipschitz
constant 1. Hence


r
h
i
X
γ 
dw (Π̃γ , Π̌) ≥ E f0 (η, ϑ̃) − f0 (η, ϑ̌) =
E
|Zj |
d
j: ηj =0
r


2p
1
γd 1 − E(kηk0 ) ,
=
π
d
and the result is proved.



Lemma 13. Assume H1 and fix δ ∈ ∆. For all θ ∈ Rd ,
1
1
h(θδ |δ) +
kθ − θδ k2 ≥ hγ (θ|δ) ≥ h(θδ |δ) +
kθ − θδ k2 − rγ (θ|δ),
(37)
2γ
2γ
with
γ
def
rγ (δ, θ) = h∇`(θ) − ∇`(θδ ), θ − θδ )i + kδ · ∇`(θ) + δ · g (θδ |δ) k2 ,
2
and where g(θδ |δ) denotes a sub-gradient of P (·|δ) at θδ . It follows in particular that
for all θ ∈ Rd , hγ (θ|δ) ↑ h(θ|δ), as γ ↓ 0.
Proof. From the definition we have


1
2
hγ (θ|δ) = min `(θ) + h∇`(θ), u − θi + P (u|δ) +
ku − θk
2γ
u∈Rd
1
≤ `(θ) + h∇`(θ), θδ − θi + P (θδ |δ) +
kθ − θδ k2 .
2γ
By convexity of `, `(θ) + h∇`(θ), θδ − θi ≤ `(θδ ), which proves the first inequality in
(37). To prove the second inequality, we start by using again the convexity of ` to
write for all θ ∈ Rd ,
`(θ) ≥ `(θδ ) + h∇`(θδ ), θ − θδ i .
Hence for all θ ∈ Rd ,
`(θ) + h∇`(θ), Jγ (θ|δ) − θi ≥ `(θδ ) + h∇`(θδ ) − ∇`(θ), θ − θδ i
+ h∇`(θ), Jγ (θ|δ) − θδ i . (38)
By H1, P (·|δ) is convex, and if g(θδ |δ) denotes a sub-gradient of P (·|δ) at θδ , we have
P (Jγ (θ|δ)|δ) ≥ P (θδ |δ) + hg (θδ |δ) , Jγ (θ|δ) − θδ i .

(39)

26

YVES F. ATCHADÉ

(38)-(39) together with the expression (11) of hγ imply that
hγ (θ|δ) ≥ h(θδ |δ) − h∇`(θ) − ∇`(θδ ), θ − θδ i
+ h∇`(θ) + g (θδ |δ) , Jγ (θ|δ) − θδ i +

1
kθ − Jγ (θ|δ)k2 .
2γ

Since Jγ (θ|δ) ∈ Rdδ , we can split kθ − Jγ (θ|δ)k2 as kθ − θδ k2 + kθδ − Jγ (θ|δ)k2 . We use
this in the last inequality to conclude that
1
kθ − θδ k2 − h∇`(θ) − ∇`(θδ ), θ − θδ i
2γ
1
+ h∇`(θ) + g (θδ |δ) , Jγ (θ|δ) − θδ i +
kJγ (θ|δ) − θδ k2
2γ
1
≥ h(θδ |δ) +
kθ − θδ k2 − h∇`(θ) − ∇`(θδ ), θ − θδ i
2γ
γ
− kδ · ∇`(θ) + δ · g (θδ |δ) k2 ,
2

hγ (θ|δ) ≥ h(θδ |δ) +

as claimed. In the last inequality, the δ appearing in front of ∇`(θ) + g(θ|δ) comes
from the fact that Jγ (θ|δ) − θδ ∈ Rdδ .
It is obvious from its definition that hγ (θ|δ) is non-decreasing as γ ↓ 0. If θ ∈
/ Rdδ ,
then kθ − θ · δk > 0, and then both extreme sides of (37) converges to +∞ = h(θ|δ)
as γ ↓ 0. If θ ∈ Rdδ , then kθ − θ · δk = 0 and both extreme sides of (37) converges to
h(θ · δ|δ) = h(θ|δ) as γ ↓ 0.

Lemma 14. Assume H1. Suppose that there exists γ0 > 0 such that Π̌γ0 (·|z) is
well-defined. Then for all γ ∈ (0, γ0 ], Π̌γ (·|z) is well-defined and


dtv (Π̌γ , Π̃γ ) ≤ 2 1 − e−%γ (z) .

(40)

Proof. For all γ > 0, we define
def

Z

Cγ (δ) =

Rd

e−hγ (θ|δ) dθ,

and Cγ =

X

πδ (2πγ)

kδk0
2

Cγ (δ).

δ

The term Cγ is the normalizing constant of Π̌γ . The function hγ is nondecreasing as
γ ↓ 0. Hence, if Cγ0 < ∞, then Cγ < ∞ for all γ ∈ (0, γ0 ], which guarantees that Π̌γ
is well-defined for all γ ∈ (0, γ0 ]. For the remaining of the proof, we fix γ ∈ (0, γ0 ].
To derive the total variation majoration, we start with a bound on Cγ . Using the

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

27

second inequality of (37), we write
(2πγ)−d/2 Cγ

=

X


πδ

δ

≤

X


πδ

δ

1
2πγ
1
2πγ

 d−kδk0 Z
2

e−hγ (θ|δ) dθ

Rd

 d−kδk0 Z
2

erγ (δ,θ) e

1
− 2γ
kθ−θδ k2 −h(θδ |δ)

e

dθ.

Rd

where
rγ (δ, θ) = h∇`(θ) − ∇`(θδ ), θ − θδ )i +

γ
kδ · ∇`(θ) + δ · g (θδ |δ) k2 .
2

In view of this last inequality, and the definitions of Π̃γ , and %γ , we get
(2πγ)−d/2 Cγ
≤ e%γ (z) .
C

(41)

The total variation bound between Π̃γ (δ, dθ|z) and Π̌γ (δ, dθ|z) now follows from a
comparison of the two measures. Indeed, Using the first inequality of (37), and for
γ ∈ (0, γ0 ], we deduce that
1
πδ
Cγ

Π̌γ (δ, dθ|z) =



1
2πγ

− kδk1
2

e−hγ (θ|δ) dθθ

 kδk1

1
1 − 2 − 2γ1 kθ−θδ k2 −h(θδ |δ)
e
e
dθ
πδ
Cγ
2πγ
C
Π̃γ (δ, dθ|z),
d
(2πγ)− 2 Cγ

≥
=

≥ e−%γ (z) Π̃γ (δ, dθ|z),

(42)

using (41). By a standard coupling argument (see e.g. Lindvall (1992) Equation 5.1),
the minorization (42) implies (40).

7.3. Proof of Theorem 8. It suffices to establish the stated bound on %γ (z), and
apply Theorem 6. From its definition, we have
P
e%γ (z) =

δ∈∆ πδ

P



δ∈∆

1
2πγ



 d−kδk0 R
2

1
2πγ

Rd

 d−kδk0 R
2

1
− 2γ
kθ−θδ k2 −h(θδ |δ)
e
dθ

erγ (δ,θ) e

Rd

.
− 1 kθ−θδ k2 −h(θδ |δ)
e
dθ
e 2γ

28

YVES F. ATCHADÉ

It follows from H2 that
3γ
k∇`(θ) − ∇`(θδ )k2
2

rγ (δ, θ) ≤ h∇`(θ) − ∇`(θδ ), θ − θδ i +

3γ
3γ
kδ · ∇`(θδ )k2 +
kg(θδ |δ, φ)k2
2
2

3γ
3γ
≤ L1 1 +
L1 kθ − θδ k2 + 3γL2 `(θδ ) +
c(δ) + 3γL2 P (θδ |δ).(43)
2
2


def
def
We set hγ = 1 − 2γL1 1 + 3γ
L
2 1 , and a = 3L2 . Then (43) gives
+

Z

1
kθ−θδ k2 h(θδ |δ)
− 2γ

erδ (δ,θ) e

e

Rd

3γ

dθ ≤ e 2 c(δ)
Z
h
− γ kθ−θδ k2 −(1−γa)`(θδ )−(1−γa)P (θδ |δ)
×
e 2γ
e
dθ. (44)
Rd

Notice that the integral on the right-side of (44) can be factorized as the product
of two integrals, with one integral taken over the components for which δj = 0, and
the other taken over the components for which δj = 1. We introduce some notation
to do this rigorously. Fix δ ∈ ∆, and s = kδk0 . For a given function f : Rd → R, we
define f [s] : Rs → R as f [s] (u) = f (uδ ), where uδ ∈ Rd , and uδi = 0 if δi = 0, and
uδj = uPj δ if δj = 1. With this notation, and for 4γL1 ≤ 1 (which implies that
k=1 k
hγ > 0), the integral on the right-hand side of (44) is equal to


2πγ
hγ

 d−s Z
2

[s] (u)−(1−γa)P [s] (u|δ)

e−(1−γa)`

du.

Rs

A similar calculation on the denominator of e%γ (z) gives
Z
Z
d−s
[s]
[s]
− 1 kθ−θδ k2 −h(θδ |δ)
e 2γ
e
dθ = (2πγ) 2
e−` (u)−P (u|δ) du.
Rd

Rs

We conclude that
P
e%γ (z) ≤

δ∈∆ πδ e

3γ
c(δ)
2



P

1
hγ

 d−s R
2

Rs

R

δ∈∆ πδ Rs

e−(1−γa)`

[s] (u)−(1−γa)P [s] (u|δ)

[s] (u)−P [s] (u|δ)

e−`

du

du
,

(45)

For 4γL1 ≤ 1, and using the inequality log(1 − 2x − 3x2 ) ≥ −6x, valid for all
x ∈ [0, 1/4], we have


1
hγ

 d−s
2





d−s
2 2
= exp −
log 1 − 2γL1 − 3γ L1 ≤ e3dγL1 .
2

(46)

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

29

Fix u0 ∈ Rs , arbitrary. Since γ > 0 is taken such that 4γL2 ≤ 1, we see that
γa = 3γL2 ≤ 3/4. Then by the convexity of `[s] we have
(1 − γa)`[s] (u) = −γa`[s] (u0 ) + (1 − γa)`[s] (u) + γa`[s] (u0 )
≥ −γa`[s] (u0 ) + `[s] (γau0 + (1 − γa)u) .
Similarly, by the convexity of P [s] (·|δ),
(1 − γa)P [s] (u|δ) ≥ −γaP [s] (u0 |δ) + P [s] (γau0 + (1 − γa)u|δ) .
Using these last two inequalities, and the change of variable (1 − γa)u + γau0 = w,
we conclude that
Z
[s]
[s]
e−(1−γa)` (u)−(1−γa)P (u|δ) du
Rs
Z
[s]
[s]
[s]
[s]
e−` (u)−P (u|δ) du.
≤ eγa(` (u0 )+P (u0 |δ)) (1 − γa)−s
Rs



def
Setting R(z) = maxδ∈∆ inf u∈Rs `[s] (u) + P [s] (u|δ) , and using the inequality log(1 −
3x) ≥ −6x, x ∈ [0, 1/4] we obtain,
Z
Z
(s)
¯
¯(s)
−(1−γa1 )`(u)−(1−γa
γaR(z) 6dγL2
2 )P̄ (u)
e
du ≤ e
e
e−` (u)−P̄ (u|δ) du.
Rs

Rs

It follows from this last inequality, (46) and (45) that
%γ (z) ≤

3γ
max c(δ) + 3γd(L1 + 2L2 ) + 3γL2 R(z),
2 δ∈∆

as claimed.

7.4. Proof of Corollary 9. We show that H1-H2 hold and apply Theorem 8.
The function ` is clearly convex and ∇`(θ) = − σ12 X 0 (z − Xθ). Hence H1(1) holds.
The elastic-net density in (22) is log-concave and continuous, which implies that
P (·|δ) is convex and lower semi-continuous for any given δ. Furthermore, For θ ∈ Rdδ ,
def αλ1
sign(θ)
σ2

sign(θ) is a subgradient of x 7→ kxk1 at θ. Hence g(θ|δ) =
subgradient of P (·|δ) at θ ∈ Rdδ . Hence H1 holds.
From the expression of ∇`, we have

+

(1−α)λ2
θ
σ2

k∇`(θ2 ) − ∇`(θ1 )k ≤ L1 kθ − θ2 k.
def

with L1 = λmax (X 0 X)/σ 2 . Furthermore, for all δ ∈ ∆ and θ ∈ Rdδ ,
kδ · ∇`(θ)k2 =

1
1
(z − Xθ)0 Xδ Xδ0 (z − Xθ) ≤ 2L1 2 kz − Xθk2 .
4
σ
2σ

is a

30

YVES F. ATCHADÉ

From the expression of g(·|δ), we have
kg(θ|δ)k



2



λ1
2(1 − α)λ2
λ2
2
≤
α 2 kθk1 + (1 − α) 2 kθk
kδk0 +
σ2
σ
2σ
2(1 − α)λ2
= c(δ) +
P (θ|δ),
σ2
≤ c(δ) + 2L1 P (θ|δ), θ ∈ Rdδ

2

αλ1
σ2


2
def
1
kδk0 , and using the assumption (1 − α)λ2 ≤ λmax (X 0 X). These
where c(δ) = αλ
σ2
inequalities show that H2 holds. The corollary then follows from Theorem 8 by noting

2
2
1
that maxδ c(δ) ≤ αλ
d, and R(z) ≤ `(0) = kzk
.
2
σ
2σ 2


def

7.5. Proof of Theorem 10. Since  = Z −Xθ? ∼ N(0, σ 2 In ), by standard Gaussian
tail bound and union bound inequalities, we have




λ1
λ21
2
P [Z ∈
/ E] = P max | hXk , i | >
≤ 2 exp log(d) −
= ,
(47)
2
1≤k≤d
2
8nσ κ̄(1)
d
p
given the choice λ1 = 4σ κ̄(1)n log(d). By Jensen’s inequality,
h
i
E [%γ (Z)|Z ∈ E] ≤ log E e%γ (Z) |Z ∈ E


h
i
2
≤ − log 1 −
+ log E e%γ (Z) 1E (Z) , on Z ∈ E.
d
In the particular case of the linear model, we have
P
e%γ (Z) =

δ πδ



λ1
2σ 2

kδk0 R

 d−kδk0 − 1 kθ−θδ k2 − 1 kZ−Xθδ k2 − λ1 kθδ k1
1
2
erγ (δ,θ) 2π
e 2γ
e 2σ2
e σ2
dθ


λ
kδk0 R
P
− 12 kZ−Xδ uk2 − 12 kuk1
λ1
2σ
e σ
du
δ πδ 2σ 2
Rkδk0 e
Rd

Using Lemma 17 of Atchadé (2015), the denominator of e%γ (Z) satisfies the lower
bound
e

1
kZ−Xθ? k2
2σ 2

e

λ1
kθ? k1
σ2

X
δ


πδ

λ1
2σ 2

kδk0 Z

1

2

λ1

e− 2σ2 kZ−Xδ uk e− σ2 kuk1 du

Rkδk0


−s?
nσ 2 κ̄(s? )
≥ π δ? 1 +
. (48)
λ21

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS
def 1
X 0X
σ2

As in the proof of Theorem 6, and setting Mγ =

rδ (δ, θ) ≤ h∇`(θ) − ∇`(θδ ), θ − θδ i +



Id +

3γ
X 0X
2σ 2



31

, we have

3γ
k∇`(θ) − ∇`(θδ )k2
2

3γ
3γ
kδ · ∇`(θδ )k2 +
kg(θδ |δ, φ)k2
2
2
≤ (θ − θδ )0 Mγ (θ − θδ )
+

+

3γ
3γ
λmax (Xδ Xδ0 )kZ − Xθδ k2 +
2σ 4
2



λ1
σ2

2
kδk0 .

def

def

Set aδ = (3/σ 2 )λmax (Xδ Xδ0 ), and Hγ = I − 2γMγ . Using the last inequality, we get
the bound
Z
λ
− 1 kθ−θδ k2 − 12 kZ−Xθδ k2 − 12 kθδ k1
erγ (δ,θ) e 2γ
e σ
dθ
e 2σ
Rd

≤e

3γ
2



λ1
σ2

2

kδk0

Z
e

1
− 2γ
(θ−θδ )0 Hγ (θ−θδ ) −

e

1−γaδ
2σ 2

kZ−Xθδ k2 −

λ1
kθδ k1
σ2

dθ. (49)

Rd

And if we call Jδ (Z) the integral on the right-side of (49), then by Fubini’s theorem,

E 1E (Z)e

1
kZ−Xθ? k2
2σ 2

λ1
kθ? k1
σ2



Z

−

1

(θ−θ )0 H (θ−θ )

λ1

γ
δ
δ
Jδ (Z) =
e 2γ
e− σ2 [kθδ k1 −kθ? k1 ]
d
R



1
1 − γaδ
2
2
× E 1E (Z) exp
kZ − Xθ? k −
kZ − Xθδ k
dθ. (50)
2σ 2
2σ 2

e

We write
1
1
kZ − Xθδ k2 = 2 kZ − Xθ? k2 +
2
2σ
2σ




1 0
X (Z − Xθ? ), θδ − θ?
σ2
1
+ 2 (θδ − θ? )0 (X 0 X)(θδ − θ? ),
2σ

and for Z ∈ E, | σ12 X 0 (Z − Xθ? ), θδ − θ? | ≤ (λ1 /2σ 2 )kθδ − θ? k1 . Therefore the
expectation on the right-side of (50) is upper bounded by
λ1

e 2σ2 kθδ −θ? k1 e−

1−γaδ
2σ 2

(θδ −θ? )0 (X 0 X)(θδ −θ? )

=e

h γaδ
i
2
E e 2σ2 kZ−Xθ? k

λ1
kθδ −θ? k1
2σ 2

−

e

1−γaδ
2σ 2

(θδ −θ? )0 (X 0 X)(θδ −θ? )



1
1 − γaδ

n/2

Letting
def

B(θ) = −

λ1
λ1
1 − γaδ
[kθk1 − kθ? k1 ] + 2 kθ − θ? k1 −
(θ − θ? )0 (X 0 X)(θ − θ? ),
2
σ
2σ
2σ 2

.

32

YVES F. ATCHADÉ

for θ ∈ Rd , we conclude that

E 1E (Z)e

1
kZ−Xθ? k2
2σ 2

e

λ1
kθ? k1
σ2


Jδ (Z)


≤

1
1 − γaδ

n/2 Z
e

1
− 2γ
(θ−θδ )0 Hγ (θ−θδ ) B(θδ )

e

dθ. (51)

Rd

Using an argument that can be found in Castillo et al. (2015) (proof of Theorem 10),
and also in Atchadé (2015) (proof of Lemma 5), it can be shown that the function B
satisfies
B(θ) ≤

λ1
8λ21 s?
− 2 kθ − θ? k1 , θ ∈ Rd .
2
nσ κ
4σ

(52)

Combining (52), (51), (49), and (48), we conclude that

s
 2
2s
i
h
3γ λ1
1 ? X
1
nσ 2 κ̄(s? ) ? 8λ
kδk0 − n log(1−γaδ )
2
nσ 2 κ
σ2
π
e
e 2
E 1E (Z)e%γ (Z) ≤
1+
e
δ
2
πδ?
λ1
δ


×

λ1
2σ 2

kδk0 

1
2π

 d−kδk0 Z
2
e

1
− 2γ
(θ−θδ )0 Hγ (θ−θδ ) −

e

λ1
kθδ −θ? k1
4σ 2

dθ. (53)

Rd

Since (γ/σ 2 )λmax (X 0 X) ≤ 1/4, and using the inequality log(1 − 2x − 3x2 ) ≥ −6x,
valid for all x ∈ [0, 1/4], it can be shown that the integral on the right-side of (53) is
upper bound by
 2 kδk0
d−kδk0 3γ
0
8σ
(2π) 2 e σ2 Tr(X X) .
λ1
We conclude that
h

%γ (Z)

E e

s

2s
i
1 ? 3γ
0
1
nσ 2 κ̄(s? ) ? 8λ
2κ
nσ
e
1E (Z) ≤
1+
e σ2 Tr(X X)
2
πδ?
λ1
 2
X
3γ λ1
kδk0 − n log(1−γaδ ) log(4)kδk0
2
σ2
e
.
πδ e
e 2
δ

Since − n2 log(1 − γaδ ) ≤ 3γnλmax (X 0 X)/σ 2 , and with λ1 = 4σ
h

%γ (Z)

log E e

p

κ̄(1)n log(d),




κ̄(s? )
κ̄(1)
+ 128
log(d)
16 log(d)
κ
 2
X
3γ λ1

3γ
kδk0 log(4)kδk0
+ 2 nλmax (X 0 X) + Tr(X 0 X) + log
πδ e 2 σ2
e
.
σ


1E (Z) ≤ − log(πδ? ) + s? log 1 +
i

δ

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS
def

Set A = log(4) +

3γ
2



log

λ1
σ2

2

X

for

24γ
κ̄(1)
σ2

. Then
Akδk0

πδ e

δ

33

d  
X
d s
= log
q (1 − q)d−s eAs
s
s=0

= d log 1 + qeA (1 − e−A ) ,
 2
λ1
≤ 4 log(4) + 6γ
,
σ2

≤ u − 1. Also,
− log(πδ? ) ≤ us? log(d) +

2
du−1

.

The theorem is proved.


References
Armagan, A., Dunson, D. B. and Lee, J. (2013). Generalized double Pareto
shrinkage. Statist. Sinica 23 119–143.
Atchadé, Y. F. (2006). An adaptive version for the metropolis adjusted langevin
algorithm with a truncated drift. Methodol Comput Appl Probab 8 235–254.
Atchadé, Y. F. (2015). On the contraction properties of some high-dimensional
quasi-posterior distributions. ArXiv e-prints .
Bauschke, H. H. and Combettes, P. L. (2011). Convex analysis and monotone operator theory in Hilbert spaces. CMS Books in Mathematics/Ouvrages de
Mathématiques de la SMC, Springer, New York.
Bottolo, L. and Richardson, S. (2010). Evolutionary stochastic search for
bayesian model exploration. Bayesian Anal. 5 583–618.
Buhlmann, P. and Mandozzi, J. (2014). High-dimensional variable screening and
bias in subsequent inference, with an empirical comparison. Computational Statistics 29 407–430.
Castillo, I., Schmidt-Hieber, J. and van der Vaart, A. (2015). Bayesian linear
regression with sparse priors. Ann. Statist. 43 1986–2018.
Chen, X., Wang, Z. J. and McKeown, M. J. (2011). A bayesian lasso via
reversible-jump {MCMC}. Signal Processing 91 1920 – 1932.

34

YVES F. ATCHADÉ

Dudley, R. (2002). Real Analysis and Probability. Cambridge Series in advanced
mathematics, Cambridge University Press, NY.
Ge, D., Idier, J. and Carpentier, E. L. (2011). Enhanced sampling schemes for
MCMC based blind bernoulli-gaussian deconvolution. Signal Processing 91 759 –
772.
George, E. I. and McCulloch, R. E. (1997). Approaches to bayesian variable
selection. Statist. Sinica 7 339–373.
Ishwaran, H. and Rao, J. S. (2005). Spike and slab variable selection: Frequentist
and bayesian strategies. Ann. Statist. 33 730–773.
Li, Q. and Lin, N. (2010). The bayesian elastic net. Bayesian Anal. 5 151–170.
Lindvall, T. (1992). Lectures on the coupling method. John Wiley & Sons, Inc.,
New York.
Mitchell, T. J. and Beauchamp, J. J. (1988). Bayesian variable selection in linear
regression. JASA 83 1023–1032.
Moreau, J.-J. (1965). Proximité et dualité dans un espace hilbertien. Bull. Soc.
Math. France 93 273–299.
Narisetty, N. and He, X. (2014). Bayesian variable selection with shrinking and
diffusing priors. Ann. Statist. 42 789–817.
Ormerod, J. T., You, C. and Muller, S. (2014). A variational Bayes approach
to variable selection. Tech. rep., Preprint.
Patrinos, P., Stella, L. and Bemporad, A. (2014). Forward-backward truncated
Newton methods for convex composite optimization. ArXiv e-prints .
Pereyra, M. (2015). Proximal Markov chain Monte Carlo algorithms. Statistics
and Computing (To Appear) http://dx.doi.org/10.1007/s11222–0159567–4.
Reid, S., Tibshirani, R. and Friedman, J. (2013). A Study of Error Variance
Estimation in Lasso Regression. ArXiv e-prints .
Rockova, V. and George, E. I. (2014). Emvs: The em approach to bayesian
variable selection. Journal of the American Statistical Association 109 828–846.
Schreck, A., Fort, G., Le Corff, S. and Moulines, E. (2013). A shrinkagethresholding Metropolis adjusted Langevin algorithm for Bayesian variable selection. ArXiv e-prints .
Shun, Z. and McCullagh, P. (1995). Laplace approximation of high-dimensional
integrals. J. Roy. Statist. Soc. Ser. B 57 749–760.
Tierney, L. and Kadane, J. B. (1986). Accurate approximations for posterior
moments and marginal densities. J. Amer. Statist. Assoc. 81 82–86.
Yang, Y., Wainwright, M. J. and Jordan, M. I. (2015). On the Computational
Complexity of High-Dimensional Bayesian Variable Selection. ArXiv e-prints .

ON THE MOREAU-YOSIDA APPROXIMATION OF POSTERIOR DISTRIBUTIONS

35

Zhang, C.-H. (2010). Nearly unbiased variable selection under minimax concave
penalty. The Annals of Statistics 38 894–942.
Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic
net. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67
301–320.

