Bayesian Lasso : Concentration and MCMC Diagnosis
D.Ounaissi1 , N.Rahmania1

arXiv:1802.08572v1 [] 22 Feb 2018

Laboratoire LAMA, UFR de Mathématiques, Université de Paris 12, France.
Laboratoire Paul Painlevé, UFR de Mathématiques, Université de Lille, France.

Abstract
Using posterior distribution of Bayesian LASSO we construct a semi-norm
on the parameter space. We show that the partition function depends on the
ratio of the l1 and l2 norms and present three regimes. We derive the concentration of Bayesian LASSO, and present MCMC convergence diagnosis.
Keywords: LASSO, Bayes, MCMC, log-concave, geometry, incomplete
Gamma function
1. Introduction
Let p ≥ n be two positive integers, y ∈ Rn and A be an n × p matrix
with real numbers entries. Bayesian LASSO

 kAx − yk2
1
2
(1)
c(x) = exp −
− kxk1
Z
2
is a typically posterior distribution used in the linear regression
y = Ax + w.
Here
 kAx − yk2

2
exp −
− kxk1 dx
2
Rp

Z
Z=

(2)

is the partition function, k · k2 and k · k1 are respectively the Euclidean
and the l1 norms. The vector y ∈ Rn are the observations, x ∈ Rp is the
Email addresses: daoud.ounaissi@u-pec.fr (D.Ounaissi),
nadji.rahmania@univ-lille1.fr (N.Rahmania)
1
Corresponding Author
Preprint submitted to Elsevier

February 26, 2018

unknown signal to recover, w ∈ Rn is the standard Gaussian noise, and A
is a known matrix which maps the signal domain Rp into the observation
domain Rn . If we suppose that x is drawn from Laplace distribution i.e. the
distribution proportional to
exp(−kxk1 ),

(3)

then the posterior of x known y is drawn from the distribution c (1). The
mode
o
n kAx − yk2
2
+ kxk1 : x ∈ Rp
(4)
arg min
2
of c was first introduced in [18] and called LASSO. It is also called Basis
Pursuit De-Noising method [4]. In our work we select the term LASSO and
keep it for the rest of the article.
In general LASSO is not a singleton, i.e. the mode of the distribution c
is not unique. In this case LASSO is a set and we will denote by lasso any
element of this set. A large number of theoretical results has been provided
for LASSO. See [5], [6], [9], [12], [15] and the references herein. The most
popular algorithms to find LASSO are LARS algorithm [8], ISTA and FISTA
algorithms see e.g. [2] and the review article [14].
The aim of this work is to study geometry of bayesian LASSO and to
derive MCMC convergence diagnosis.
2. Polar integration
Using polar coordinates s =

x
kxk

∈ S, r = kxk, the partition function (2)

Z
Zp =

Zp (s)ds,

(5)

S

where k·k denotes one of l2 or l1 norms in Rp , ds denotes the surface measure
on the unit sphere S of the norm k · k, and
Z +∞
Zp (s) =
exp{−f (rs)}rp−1 dr,
(6)
0
kAx−yk22

where f (x) :=
+ kxk1 , x ∈ Rp .
2
We express the partition function (6) using the parabolic cylinder function. We also give an inequality of concentration and a geometric interpretation of the partition function Zp .
2

3. Parabolic cylinder function and partition function
We extend the function s ∈ S → Zp (s)
Z +∞
p
exp{−f (rx)}rp−1 dr.
x ∈ R → Zp (x) =

(7)

0

This extension is homogeneous of order −p.
kyk2
If Ax = 0, then f (rx) = 2 2 + rkxk1 , and more if x 6= 0, then
kyk22
).
2
If Ax 6= 0, then we will express Zp (x) using the parabolic cylinder function.
We recall that for b ∈ R, z ≥ 0 the parabolic cylinder function is given by
Zp (x) = (p − 1)!kxk−p
1 exp(−

z2
)[1 + O(z −2 )],
(8)
4
when z → +∞ [16]. We also recall the integral representation of Erdlyi [7]
for the parabolic cylinder function
Z +∞
1
z2
exp(− r2 − zr)rν−1 dr, ν > 0,
exp( )Γ(ν)D−ν (z) =
4
2
0
R +∞
where Γ(ν) = 0 exp(−t)tν−1 dt is the Γ function.
Db (z) = z b exp(−

Proposition 3.1. The variable
ωlasso :=

kxk1 − hAx, yi
kAxk2

(9)

x
will play an important role. It depends only on s = kxk
∈ Sp−1,1 and the
1
1
function s ∈ Sp−1,1 → ωlasso (s) is bounded below by λ1,2 := min{ kAsk
−
2
hAs,yi
kAsk2

:

s ∈ Sp−1,1 }.

Now we can announce the following result.
Proposition 3.2. We have for Ax 6= 0
2
kyk22
ωlasso
−p
Zp (x) = (p − 1)! exp(−
)kAxk2 exp(
)D−p (ωlasso ).
2
4
If Ax → 0, then ωlasso → +∞ and

Zp (x) = (p − 1)! exp(−

kyk22
−2
)kxk−p
1 [1 + O(ωlasso )].
2
3

1
1
is bounded below by λ1,2
, where
Corollary 3.3. If y = 0 then ωlasso = kAsk
2
p
λ1,2 = max(kAsk2 : s ∈ Sp−1,1 ) is the norm of the operator A : (R , k · k1 ) →
(Rn , k · k2 ). The partition fucntion
p
Zp (s) = (p − 1)!ωlasso
exp(

2
ωlasso
)D−p (ωlasso )
4

is kAsk22 convex and decreasing.
Proof 3.4. It suffices to remark that Zp (s) =

R +∞
0

2

exp{− kAsk
r2 −r}rp−1 dr.
2

4. Geometric interpretation of the partition function
First we represent f (rx) for Ax 6= 0 in the form
f (rx) =

2
kyk22 ωlasso
(rkAxk2 + ωlasso )2
−
+
.
2
2
2

(10)

The function exp{−f (x)}, ∀x ∈ Rp is log-concav and integrable in Rp . Ob−1

serve that Zp p is a norm on the null space Ker(A) of A. A general result
[1] tells us that
−1

x ∈ Rp → Zp p (x) := kxkc
is a quasi-norm on Rp . The unit ball of k · kc is defined by
B(A, y) := {x ∈ Rp : kxkc ≤ 1}
= {x ∈ Rp : Zp (x) ≥ 1}
= {x = rs ∈ Rp : Zp (s) ≥ rp }
= {x = rs ∈ Rp : (p − 1)! exp(−

2
kyk22
ωlasso
)kAsk−p
exp(
)D−p (ωlasso ) ≥ rp }.
2
2
4

Its contour is equal to
C(A, y) := {x ∈ Rp : kxkc = 1}
= {x ∈ Rp : Zp (x) = 1}
= {x = rs ∈ Rp : Zp (s) = rp }
= {x = rs ∈ Rp : (p − 1)! exp(−

2
kyk22
ωlasso
)kAsk−p
exp(
)D−p (ωlasso ) = rp }.
2
2
4

We summarize our results in the following proposition.
4

Proposition 4.1. 1) For each s ∈ Sp−1,1 , the longest segment [0, r]s contained in B(A, y) holds for r = rmax (s) is solution of the equation
rp = (p − 1)! exp(−

2
ωlasso
kyk22
)kAsk−p
exp(
)D−p (ωlasso ).
2
2
4

2) The ball
B(A, y) =

[

[0, rmax (s)]s,

s∈Sp−1,1

and its contour is equal to
C(A, y) = {rmax (s)s :
3) The volume B(A, y) is
Z
Sp−1,1

s ∈ Sp−1,1 }.

p
Zp
rmax
(s)
ds =
.
p
p

5. Necessary and sufficient condition to have lasso equal zero
Now we can give the necessary and sufficient condition to have
lasso = {0}
Proposition 5.1. The following assertions are equivalent.
1) 0 = lasso.
2) ωlasso (s) ≥ 0 pour tout s ∈ Sp−1,1 .
3) kA> yk∞ ≤ 1.
6. Concentration around the lasso
6.1. The case lasso null
The polar coordinate formula tells us that, we can draw a vector x from
c(x)dx by drawing its angle s uniformly, and then simulate its distance r to
the origin from
c(r, s)dr =

1
exp{−f (rs)}rp−1 dr
Zp (s)

5

(11)

Now let’s estimate for r > 0 the probability
Z Z +∞
ds
c(r, s)dr ,
P(kxk > r) =
|S|
S r
where |S| denotes the Lebesgue measure of S. We introduce for each pair
a ≥ 0, b ∈ Rp the function
ga,b,p (r) := ga,b (r) − (p − 1) ln(r),

r > 0.

(12)

In the following
a = kAsk2 ,

b = ωlasso .

The function r ≥ 0 → ga,b (r) is increasing (because b := ωlasso ≥ 0). The
function r → ga,b,p (r) f est convexe et atteint son minimum au point r(s)
solution de l’quation
p−1
.
r

kAsk2 (rkAsk2 + ωlasso ) =
The positive root is given by
r(s) =

−ωlasso +

p
2
ωlasso
+ 4(p − 1)
2kAsk2

On one hand
Z +∞
Z
exp{−ga,b,p (r)}dr ≥ exp{−ga,b (r(s))}
0

(13)

r(s)

rp−1 dr = exp{−ga,b,p (r(s))}

0

On the other hand by using the convexity of r → ga,b (r), we have for all
r > 0,
ga,b (r) ≥ ga,b (r(s)) +

(p − 1)(r − r(s))
,
r(s)

p−1
because ∂r ga,b (r(s)) = r(s)
. We deduce for q > 0,
Z +∞
Z
exp{−ga,b,p (r)}dr ≤ exp{−ga,b (r) + p − 1}
qr(s)

Z

+∞

exp{−

qr(s)
+∞

≤ exp{−ga,b (r(s)) + p − 1}

exp(−r)rp−1 dr

q(p−1)

≤ exp{−ga,b (r(s)) + p − 1}
6

p−1
r}rp−1 dr
r(s)
r(s)
(p − 1)p

r(s)
Γ(p, q(p − 1)),
(p − 1)p

r(s)
.
p

R +∞
where Γ(ν, x) = x exp(−t)tν−1 dt is the upper incomplete gamma function.
Finally we get the following result.
Proposition 6.1. We have for all q > 0,
P(kxk ≥ qr(s)) ≤

p exp(p − 1)
Γ(p, q(p − 1)) := P (q, p).
(p − 1)p

(14)

Using the following estimate [? ]
xa−1 exp(−x) < Γ(a, x) < Bxa−1 exp(−x),

∀ a > 1, B > 1, x >

B
(a − 1),
B−1

we get for q > 1,
Γ(p, q(p − 1)) ≤ 2q p−1 (p − 1)p−1 exp(−q(p − 1)).
Therefore the quantity
2pq p−1
P (q, p) ≤
exp{−(q − 1)(p − 1)}.
(p − 1)
x
balance sheet . If x is drawn from the density c, alors r(θ)
∈ B2 (0, q)
with a probability at least equal to 1 − P (q, p).

In the figure(1) we plot for p = 2, n = 1, A = 1 1 and y = 0 the
density c(r, s)dr for a fixed value of s.

We notice that the mode c(r, s) = 0.6200 is very close to the value of
r( s) = 0.6290 (13) for the same fixed s.
6.2. The general case
We take the vector l ∈ lasso. We will study the concentration of c around
l. The variable of interest is u = x − l. The law of u has for density
c(u + l)du =

1
exp{−f (u + l, A, y)}du.
Zp

The change of variables formula gives for each norm k · k
c(u + l)du =

1
exp{−f (rθ + l)}rp−1 drdθ,
Zp
7

r > 0, θ ∈ S.

Figure 1: pour r ∈ [0.1; 10] → c(r, s).

By definition for any vector x, the convex function r ≥ 0 → f (rs + l) reaches
its minimum at the point r = 0. Therefore r ≥ 0 → f (rs + l) is increasing.
The function
f (rs + l, p) := f (rs + l) − (p − 1) ln(r),

r > 0,

(15)

is strictly convex. Its critical point rl (s) is solution of the equation
p−1
∂r f (rs + l) =
.
r
By a similar proof to that of propostion (6.1) we have the following result;
x−l
Proposition 6.2. If x is drawn from the density c, and s = kx−lk
, then for
all q > 0,
p exp(p − 1)
P(kx − lk ≥ qrl (s)) ≤
Γ(p, q(p − 1)) := P (q, p).
(16)
(p − 1)p

7. Applications
7.1. The contour in the case p = 2, n = 1
Let A := (a1 , a2 ) a matrix of order 1 × 2. Its null-space Ker(A) =
{(x1 , x2 ) : a1 x1 + a2 x2 = 0}. We have that B(a1 , a2 , y) contains
\
Ker(A) B2,1 .
8

This intersection is a symmetric segment noted [(x1 (a1 , a2 ), x2 (a1 , a2 )), −(x1 (a1 , a2 ), x2 (a1 , a2 ))].
To determine the other points of the set B(a1 , a2 , y), we will directly
calculate Z2 (s). A simple calculation gives
Z +∞
2
(|As|r + ωlasso )2
y 2 ωlasso
exp{−
Z2 (s) = exp(− +
)
}rdr,
2
2
2
0
et
Z
|As|
0

+∞

(|As|r + ωlasso )2
}rdr + ωlasso
exp{−
2

Z

+∞

exp{−
0

(|As|r + ωlasso )2
}dr = 1.
2

Finally we have the following proposition.
Proposition 7.1. 1) If As 6= 0, then
Z2 (s) = exp(−

2
y 2 ωlasso
ωlasso √
+
)|As|−1 {1 −
2π(1 − F (ωlasso ))},
2
2
|Aω|

where F is the distibution function of the normal law.
2) If As 6= 0 and y = 0, then
Z2 (s) = ωlasso exp(

2
√
ωlasso
2
){1 − ωlasso
2π(1 − F (ωlasso ))}.
2

3) Ifi s ∈ S1,1 , As 6= 0 and y = 0, then the function z2
z2 (b2 ) =

1
1√
1
1
exp( 2 ){1 − 2 2π(1 − F ( ))}
b
2b
b
b

defined on (0, λ21,2 ] is convex and decreasing, where λ1,2 = maxs∈S1,1 |As|.
4) We have for s ∈ S1,1
Z2 (s) = z2 (

1
2
ωlasso

),

∀ ω ∈ S1,1 .

the ball
Z2 (s) ≥ r2 }.

B(A, 0) = {rs :

is contained in the unit disk kxk1 ≤ 1 for the norm l1 . The contour is defined
by the equation
Z2 (s) = r2 .
9

The norm of the linear operator A : (R2 , k · k1 ) → (R, k · k2 ) is defined by
λ1,2 = sup kAsk2 .
s:ksk1 =1

the function s → Z2 (s) = z2 (λ21,2 ) is constant on est constante sur
Ω1,2 = {s :

ksk1 = 1,

kAsk = λ1,2 }.

If A = (1, 1) then
Ω1,2 = {s :

ksk1 = 1, kAsk = λ1,2 }
[
= [(1, 0), (0, 1)] [(−1, 0), (0, −1)].

If A = (a1 , a2 ) with |a1 | < |a2 |, then
Ω1,2 = {s : ksk1 = 1, kAsk2 = λ1,2 }
= {(0, sgn(a2 )), (0, −sgn(a2 ))}.
In both case
1

{z2 (λ21,2 )} 2 Ω1,2
is part of the contour. The other points of the contour are deduced from the
equation
z2 (b2 ) = a2 ,

b ∈ (0, λ1,2 ).

Each pair (a, b) generate four points of B((a1 , a2 ), 0) of the form as where
|s1 | + |s2 | = 1,

|a1 s1 + a2 s2 | = b.

We plot in the figure 2 the contour of B(a1 , a2 , 0) for different choices of
the matrix (a1 , a2 ). We notice that the surface of B((a1 , a2 ), 0) is decreasing
function relatively the norm λ1,2 of the matrix A.
Remark 7.2. The numerics show thatZ(ωlasso ) exploses for the large values
of ωlasso , it means that ω is closes to the null-space of A. to eleminate that
explosion we need to estimate the tail of the gaussian density. Using the
Gordon estimation [10]
2

2

exp(− x2 )
exp(− x2 ) √
≤
2π(1
−
F
(x))
≤
,
x
x + x1
we have the following approximation
1
1
1
1
− ≤ z2 (b2 ) ≤ 2 −
,
2
b
b
b
1 + b3
10

x > 0,

near to 0.

(17)

Figure 2: Contours of B(a1 , a2 , 0) for different matrices A, p = 2 and n = 1.

8. MCMC diagnosis
Here we take p = 7, n = 4, A ∼ B(± √1n ) and for simplicity we consider
y = 0. We sample from the distribution c (1) using Hastings-Metropolis
algorithm (x(t) ) and propose the test kx(t) k2 ≤ qr(θ(t) ) as a criterion for
(t)
the convergence. Here θ(t) := kxx(t) k2 . We recall that if x is drawn from the
target distribution c, then kxk2 ≤ qr(θ) with the probability at least equal
to P (q, p). Table 2 gives the values of the probability P (q, p). Note that for
q ≥ 2.5 the criterion kx(t) k2 ≤ qr(θ(t) ) is satisfied with a large probability.
q
2
2.5
3
3.5
4
4.5
5
P (q, p) 0.6672 0.9446 0.9924 0.9991 0.9999 1.0000 1.0000
Table 1: Values of the probability P (q, p) for p = 7.

8.1. Independent sampler (IS)
The proposal distribution
Q(x2 , x1 ) = p(x2 ) =

1
exp(−kx2 k1 ),
2p
11

∀ x1 , x2 .

The ratio
c(x)
2p
≤ ,
p(x)
Z

∀ x.

It’s known that MCMC (x(t) ) with the target distribution c and the proposal
distribution p is uniformly ergodic [13]:
Z
Z
(t)
(0)
c(x)dx| ≤ (1 − p )t .
sup |P(x ∈ A | x ) −
2
A⊂B(Rp )
A
Here Z ≈ 2.2142 and then (1 − 2Zp ) = 0.9827. Figure 4(a) shows respectively
the plot of t → 5r(θ(t) ) and t → kx(t) k2 .
8.2. Random-walk (RW) Metropolis algorithm
We do not know if the target distribution c satisfies the curvature condition in [17] Section 6. Here we propose to analyse the convergence of the Random walk Metropolis algorithm (x(t) ) using the criterion kx(t) k2 ≤ qr(θ(t) ).
Figure 4(b) shows respectively the plot of t → 5r(θ(t) ) and t → kx(t) k2 .
Figures 4 show that contrary to independent sampler algorithm, the random walk (RW) algorithm satisfies early the criterion kx(t) k2 ≤ 5r(θ). More
precisely
1) the independent sampler (IS) algorithm begins to satisfy the criterion
kx(t) k2 ≤ 5r(θ(t) ) at t = 8 × 105 iteration.
2) The RW algorithm begins to satisfy the criterion kx(t) k2 ≤ 3.5r(θ(t) ) at
t = 939065 iteration, but the IS algorithm never satisfies the criterion
kx(t) k2 ≤ 3.5r(θ(t) ).
R
We finally compare IS and RW algorithms using the fact that Rp xc(x)dx =
0.
R of the integral
R The best algorithm will furnish the best1 approximation
PN (t)
x
≈
xc(x)dx and
xc(x)dx.
Table
3
gives
the
estimators
t=1 IS
N
Rp
Rp
R
P
P
(t)
(t)
N
N
1
1
t=1 xRW ≈ Rp xc(x)dx. It follows that k N
t=1 xIS k2 = 0.0187 and
N
P
(t)
N
1
k N t=1 xRW k2 = 0.0041. We conclude that the random walk algorithm
wins for both criteria against independent sampler algorithm.

12

x̂IS
x̂RW

x1
x2
x3
x4
x5
x6
x7
-0.0005 -0.0037 0.0016 0.0164 0.0050 0.0021 -0.0058
0.0005 -0.0019 -0.0002 0.0012 -0.0005 0.0031 -0.0011

Table 2: IS and RW estimators using N = 106 iterations.

Figure 3: (a): Test of convergence of MCMC algorithm with proposal distribution p(x2 ).
(b): Test of convergence of MCMC algorithm with N (0, 0.5Ip ) proposal distribution. N =
106 iterations.

9. Conclusion
We studied the geometry of bayesian LASSO using polar coordinates and
calculated the partition function. We obtained a concentration inequality
and derived MCMC convergence diagnosis for the convergence of Hasting
Metropolis algorithm. We showed that the random walk MCMC with the
variance 0.5 wins again the independent sampler with Laplace proposal distribution.

13

References
[1] K. Ball, Logarithmically concave functions and sections of convex sets
in Rn , Stud. Math. T. LXXXVIII (1988) 69-84.
[2] A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm
for linear inverse problem, SIAM J. Imag. Sci. (2009) 183-202. .
[3] E.T. Copson, Asymptotic Expansions, Cambridge University Press,
1965.
[4] S. Chen, D. L. Donoho, M. Saunders, Atomic decomposition by basis
pursuit, SIAM J. Sci. Comput. 20 (1) (1998) 33-61.
[5] I. Daubechies, M. Defrise, C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Comm. Pure
Appl. Math. 57 (11) (2004) 14131457.
[6] A. Dermoune, D. Ounaissi, N. Rahmania, Oscillation of MetropolisHastings and simulated annealing algorithms around LASSO estimator,
Math. Comput. Simulation (2015).
[7] A. Erdlyi, Higher transcendental functions, California institute of technology, bateman manuscript project, vol. 2 p. 119 (1953).
[8] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, Least angle regression,
Ann. Statist. 37 (2004) 407-499.
[9] G. Fort, S. Le Corff, E. Moulines, A. Schreck, A shrinkage-thresholding
Metropolis adjusted Langevin algorithm for bayesian variable selection,
arXiv:1312.5658v3 [] (2015).
[10] R.D. Gordon, Values of Mills’ ratio of area to bounding ordinate of the
normal probability integral for large values of the argument, Annals of
Mathematical Statistics, vol. 12 364–366 (1941).
[11] B. Klartag, V.D. Milman, Geometry of log-concave functions and measures, Geom. Dedicata 112 (1) (2005) 169-182.
[12] M. Mendoza, A. Allegra, T.P. Coleman, Bayesian Lasso Posterior Sampling Via Parallelized Measure Transport, arXiv:1801.02106v1 [stat.CO]
(2018).
14

[13] K. Mengersen, R.L Tweedie, Rates of convergence of the Hastings and
Metropolis algorithms, Ann. Statist. 24 (1994) 101-121.
[14] N. Parikh, S. Boyed, Proximal algorithms, Found. Trends Optim. 1 (3)
(2003) 123-231.
[15] M. Pereyra, Proximal Markov chain Monte Carlo algorithms, Stat. Comput. (2015) 1-16.
[16] N.M. Temme, Numerical and asymptotic aspects of parabolic cylinder
functions, ournal of Computational and Applied Mathematics 121 221–
246 (2000).
[17] G.O. Roberts, A.L. Tweedie, Geometric convergence and central limit
theorems for multidimensional Hastings and Metropolis algorithms,
Biometrika 83 (1) (1996) 95-110.
[18] R. Tibshirani, Regression shrinkage and selection via Lasso, J. R. Stat.
Soc. Ser. B Stat. Methodol. 58 (1) (1996) 267-288.
[19] R. Tibshirani, The Lasso problem and uniqueness, Electron. J. Stat. 7
(2013) 1456-1490.

15

