A shrinkage-thresholding Metropolis adjusted Langevin algorithm for Bayesian variable selection

Introduction

We focus on variable selection in regression problems: the objective is to explain a response variable with a (possibly very) large number of explanatory variables, which can be either discrete or continuous. In many applications, it is known that only a small fraction of explanatory variables explains a large fraction of the observations, and using this information is crucial for inference. Variable selection is particularly challenging in high dimensional settings.

A variety of algorithms to explore the collection of models and criteria for selecting among competing models has been proposed. In the Bayesian framework, the variable selection problem is transformed into posterior inference: rather than searching a highly hypothetical ”best” model, Bayesian analysis attempts to estimate the joint posterior distribution of the collection of all subsets of parameters. In high dimension, this aim is often overly ambitious: estimating the marginal posterior probability that a variable should be included in the model is already challenging.

In the last three decades, Markov Chain Monte Carlo (MCMC) methods have been the most commonly used computational procedures to sample posterior distributions . An early attempt to perform variable selection is the Reversible Jump MCMC (RJMCMC) introduced in . RJMCMC is a trans-dimensional sampler which produces a Markov chain evolving between spaces of different dimensions. The dimension of the sample varies at each iteration as active (nonzero) parameters are added or discarded from the model. Each new sample is accepted or rejected using a Metropolis-Hastings step where the acceptance probability is adjusted to the trans-dimensional moves. RJMCMC requires ingenuity in designing appropriate jumping rules to produce computationally efficient and theoretically effective methods. Despite many attempts , this algorithm is prone to fail when the dimension of the parameter space is large (as illustrated in our numerical section).

considers another setting that encompasses all the models jointly: at each iteration, pseudo-prior distributions are used to jointly sample regression parameters associated with all models. For high dimensional statistical problems, sampling jointly all models is of course out of reach. A more efficient algorithm, the Metropolized Carlin and Chib (MCC), simultaneously proposed by and later improved by , does not require to sample from the whole collection of models and therefore can be implemented in practice. The mixing of this algorithm depends critically on the specification of pseudo-priors, which requires also a fair amount of tuning.

Other MCMC approaches for Bayesian variable selection define a posterior distribution on the model space, where a model is a binary vector locating the active (nonzero) components of the regression vector. The objective is to estimate probabilities of activation for each regression parameter. In for example, this exploration is performed with a Gibbs sampler. Variants and adaptive versions of the Gibbs sampler for this problem have been proposed in . Samples from the posterior distribution of the models are obtained in and in with particle filters.

In this paper, we introduce a novel algorithm, the Shrinkage-Thresholding Metropolis-Adjusted Langevin Algorithm (STMALA) to perform sparse regression in high dimensional models. This algorithm might be seen as a trans-dimensional MCMC method relying on the MALA algorithm (see ). The proposal distribution in the STMALA algorithm goes as follows:

compute a noisy gradient step of the logarithm of the smooth part of the target distribution;

apply a shrinkage-thresholding operator to ensure sparsity and to shrink values of the regression parameters toward zero;

use an accept-reject step to guarantee the convergence to the correct target distribution.

Each iteration of the STMALA algorithm may be seen as a randomized version of the Shrinkage-Thresholding algorithm (see ) to guide variable selection. The Shrinkage-Thresholding algorithm (and its accelerated version FISTA) is one of the most effective method to solve sparse inverse problems. Our intuition is that a single iteration of the Shrinkage-Thresholding algorithm (with some additional noise added to ensure irreducibility) is a sensible way to visit collection of models. This intuition is supported both by very promising experimental results obtained in a variety of challenging situations and by some theoretical results. In particular, we have established the geometric ergodicity of the STMALA algorithm for a large class of target distributions. To our best knowledge, it is the first result providing a rate of convergence for a trans-dimensional MCMC algorithm (like RJMCMC and MCC); usually, only Harris recurrence is proved, see .

Our algorithm is closely related to the proximal MCMC algorithm of ; the main difference stems from the fact that our algorithm is designed to sample jointly the models and their parameters, whereas is a method to sample from high-dimensional posterior distributions with sparsity inducing priors.

This paper is organized as follows. STMALA and its application to Bayesian variable selection is described in Section 2. The geometric ergodicity of the STMALA algorithm is addressed in Section 3. Numerical experiments on simulated and real data sets to assess the performance of STMALA are given in Section 4. All the proofs are postponed to Section 6.

The Shrinkage-Thresholding MALA algorithm

Lemma 2.1 shows that the soft thresholding operator with vanishing shrinkage $\Psi_{2}$, known as the empirical Wiener operator (see ), compromises between minimizing a (non-convex) function $h$ and being close to $u$.

For any $x$, $z\mapsto q_{\Psi}(x,z)$ consists in

sampling each component of a new model $m^{\prime}\in\mathcal{M}$ as independent $\{0,1\}$-Bernoulli random variable with success parameter $\rho(\mu_{i}(x))$, $1\leq i\leq p$;

The proposal (6) is then accepted and $X^{n+1}=Z$ with probability $\alpha_{\Psi}(X^{n},Z)$ given by

otherwise, $X^{n+1}=X^{n}$. In high dimensional settings, STMALA may encounter some difficulties to accept the proposed moves. Following , we introduce a variant of the algorithm in which only a fixed number $\eta$ of components of $X^{n}$ is updated at each iteration $n$. This is achieved by combining STMALA and a Gibbs sampler in a STMALA-within-Gibbs algorithm.

In this section, we address the $V$-geometric ergodicity of the STMALA chain $(X^{n})_{n\geq 0}$ under the following assumptions: for any $m\in\mathcal{M}$,

$\omega_{m}>0$ and $\pi_{m}>0$ on $S_{m}$.

$\pi_{m}(x)\mathds{1}_{S_{m}}(x)\to 0$ when $\|x\|\to\infty$.

Let $b,\epsilon>0$ and $u\in(0,b)$. For any $m\in\mathcal{M}$ and $x\in S_{m}$, define

$W_{m}(x)$ is the cone of $S_{m}$ with apex $x-u\ n(x)$ and aperture $2\epsilon$. We will prove (see Lemma 6.6) that AA3 guarantees that, the probability to accept a move from $x$ to any point of $W_{m}(x)$ converges to one as $\|x\|$ goes to infinity.

There exist $b,R,\epsilon>0$ and $u\in(0,b)$ such that for any $m\in\mathcal{M}$, for any $x\in S_{m}\cap\{\|x\|\geq R\}$, for all $y\in S_{m}\cap W_{m}(x)$: $\pi_{m}(x-u\ n(x))\leq\pi_{m}(y)$.

When for any $m\in\mathcal{M}$, $\pi_{m}$ is differentiable on $S_{m}$, AA2 and AA3 are satisfied if (see for details), for all $m\in\mathcal{M}$,

(see [20, Section 4 and the proof of Theorem 4.3] for details).

Let $P_{\Psi}$ denote the transition kernel associated to the Hastings-Metropolis move with proposal (6) and acceptance-rejection ratio (10).

Numerical illustrations

In the examples below, $\pi$ is the posterior distribution of a regression vector in a logistic regression model; $\pi_{m}$ is the conditional distribution of the regression vector conditionally to the observations and to the model $m$.

Let $G$ be a known $N\times p$ design matrix. We have $N$ independent observations $Y=(Y_{1},\ldots,Y_{N})$ such that for all $i$, $Y_{i}$ is a Bernoulli random variable with parameter $\exp(G_{i,\cdot}X)/(1+\exp(G_{i,\cdot}X))$. Conditionally to a model $m\in\mathcal{M}$, the prior on the nonzero components of the regression vector $X\in S_{m}$ is $\mathcal{N}(0,c(G^{\prime}_{m}G_{m})^{-1})$, where $c$ is a known scaling parameter, and $G_{m}$ denotes the matrix with columns $\{G_{\cdot,i},i\in I_{m}\}$. The prior on the models $\omega_{m}$ is equal to $\theta_{\star}^{|m|}(1-\theta_{\star})^{p-|m|}$ for $\theta_{\star}\in(0,1)$. In this experiment, we choose $p=50$ and $N=100$ to assess the performance of STMALA in a simple framework; the components of $G$ are i.i.d. $\mathcal{N}(0,1)$ and $\theta_{\star}=0.05$. The algorithm is run with $c=100$, $\sigma=0.3$ and $\eta=5$. The choice of the threshold $\gamma$ in $\Psi_{2}$ is crucial (if $\gamma$ is too large, few nonzero samples are proposed and the algorithm converges slowly and if $\gamma$ is too small, the algorithm proposes non-sparse solutions that are not likely to be accepted): $\gamma$ is set to $0.4$ to get a mean acceptance rate of around $20\%$.

STMALA is used to estimate the posterior probabilities of activation of the components of $X$, defined for all $1\leq i\leq p$ as the conditional probability of the event $\{X_{i}\neq 0\}$ given the observations $Y$. The estimation is given by $\sum_{n=B}^{N_{it}+B}\mathds{1}_{\{X^{n}_{i}\neq 0\}}/N_{it}$ where $N_{it}$ is the number of iterations of the algorithm and $B$ denotes the number of iterations discarded as a burn-in period. We choose $N_{it}=50.000$ and $B=10.000$. Figure 2 (top) provides the true regression vector, the posterior mean of the regression vector given by STMALA and the estimated activation probabilities over $100$ independent Monte Carlo runs. This experiment highlights the ability of STMALA to choose the good model (the $3$ nonzero components of $X$ are recovered) and to get high posterior probabilities of activation for the selected components of $X$.

2 Linear regression

where $G$ is a $N\times p$ (known) design matrix, $E$ is a Gaussian random vector with i.i.d. standard entries and $\tau$ is the (known) precision. The prior on the models is $\omega_{m}=\theta_{\star}^{|m|}(1-\theta_{\star})^{p-|m|}$ for some (known) $\theta_{\star}\in(0,1)$. The conditional distribution of $X$ given the observations $Y$ and the model $m$ is given by

Such a posterior distribution can be obtained from the following hierarchical model: (i) given $m\in\mathcal{M}$ and positive precisions $(\vartheta_{1},\ldots,\vartheta_{p})$, the entries $X=(X_{1},\ldots,X_{p})$ are independent with distribution

The standard deviation of the RJMCMC proposal is chosen so that STMALA and RJMCMC have similar acceptance rates (between $15\%$ and $20\%$).

Figure 3 shows the true regression vector $X$ and its estimates obtained by STMALA and RJMCMC; these estimates $\hat{X}$ are defined as the posterior mean along a trajectory of length $10^{6}$ (the first $10\%$ samples are discarded). It shows that STMALA provides a sparse estimation while RJMCMC needs a lot of components to explain the observations. This is probably because RJMCMC is more or less equivalent to test each model in turn, which yields slow convergence in high dimensional settings. This slow convergence is also illustrated in Figure 4. $50$ independent trajectories of length $10^{6}$ are run; Figure 4 (top) shows the evolution of the mean number (over the $50$ runs) of active components $|m|$. RJMCMC has not converged after the $300.000$ iterations while the mean number of active components of STMALA is stable after few iterations. Figure 4 (bottom) displays the boxplots of the estimation of the first component $X_{1}$ estimated by STMALA and RJMCMC as a function of the number of iterations.

The evolution of the mean test error $\mathcal{E}_{\rm test}$ over 100 independent runs, is displayed in Figure 5 (bottom). Both figures show that RJMCMC is subject to some over fitting, which is not the case of STMALA.

3 Regression for spectroscopy data

We use the biscuits data set composed of near infrared absorbance spectra of 70 cookies with different water, fat, flour and sugar contents studied in and . The data are divided into a training data set containing measurements for $N=39$ cookies, and a test data set containing measurements for $31$ cookies. The observation model is given by

where $G$ is the design matrix, $X$ is the unknown regression vector and $E\sim\mathcal{N}(0,I)$ is the measurement noise. Each row of the design matrix consists of absorbance measurements for $p=300$ different wavelengths from $1202$ nm to $2400$ nm with gaps of $4$ nm. We compare the results obtained by STMALA with those obtained by RJMCMC for the prediction of fat content. To improve the stability of the algorithm, the columns of the matrix $G$ containing the measurements are centered and a column with each entry being equal to one is added.

The parameters of the algorithms are given by $\tau=0.5$, $\eta=15$, $\gamma=0.35$ for STMALA. The computations are made over $100$ independent trajectories of $N_{it}=2.10^{6}$ iterations, with a burn-in $B=10^{5}$. The design parameters of STMALA and RJMCMC are chosen so that the two algorithms have similar acceptance-rejection ratios (the final ratios are about $45\%$ for STMALA and $42\%$ for RJMCMC). Figure 6 shows the regression vectors $\hat{X}$ obtained by STMALA and RJMCMC, and computed as the posterior mean along one trajectory (left) and the mean regression vector estimated by STMALA and RJMCMC over $100$ independent trajectories (right).

The regression vector estimated by STMALA has a spike around $1726$ nm, which is known to be in a fat absorbance region (see ), in almost all the trajectories.

Figure 7 displays the boxplots of the $100$ independent values of the components of the regression vectors estimated by STMALA and RJMCMC associated to $9$ wavelengths close to $1726$ nm. It illustrates that the location of the spike retrieved by RJMCMC is not stable, while STMALA retrieves a spike centered at $1726$ nm in almost every trajectory.

Figure 8 (top) shows the estimated emitted signal $G\hat{X}$ obtained by STMALA and RJMCMC as a function of the observations $Y$. In this numerical experiment, STMALA provides better results than RJMCMC for both the training set and the test set. This is confirmed by Figure 8 (bottom) which displays the evolution of the mean square error (MSE) on the test dataset, defined by

as a function of the number of iterations (mean over 100 independent trajectories). The mean MSE after $2.10^{6}$ iterations is about $0.75$ for STMALA and about $1.6$ times greater for RJMCMC.

Conclusions

In this paper, we propose a new trans-dimensional MCMC algorithm to perform Bayesian variable selection in a high-dimensional regression setting. This algorithm is closely related to but is adapted to sample models which are exactly sparse in the sense that a certain number of components are equal to zero. In addition, under fairly weak assumptions, the STMALA algorithm is shown to be geometrically ergodic. In the high-dimensional setting, the STMALA algorithm outperforms the RJMCMC algorithm which is considered as the state of the art. The performance of the STMALA algorithm depends on the tuning of a set of parameters: an adaptive version is currently under investigation. Also, the algorithm has still to be adapted to the ultra large scale framework, which likely requires additional specific procedures.

Proofs

Consider first the case $p=1$. We first compute the derivative of $t$ on $\left(0,\infty\right)$ (note that $t$ is symmetric). For any $x\in\left(0,\infty\right)$,

Using straightforward computations, we get

2 Proof of Lemma 2.2

It is sufficient to compute integrals of the form

This concludes the proof for $\Psi_{1}$. The proof for $\Psi_{2}$ follows the same lines as the proof of Lemma 2.2, with the function $\psi$ replaced by $\widetilde{\psi}(z)=g\left(\gamma^{2}/\|z\|^{2}\right)\ z$.

3 Proof of Theorem 3.1

For ease of notations, we denote by $q$ the proposal distribution. Lemma 2.2 shows that for any $m\in\mathcal{M}$ and $y\in S_{m}$

where $\rho$ and $f$ are given by Lemma 2.2 and $\mu(x)=(\mu_{1}(x),\cdots,\mu_{p}(x))$ is given by (8). We start with a preliminary lemma which will be fundamental for the proofs since it allows to compare the proposal distribution $q$ to Gaussian proposals. Denote by $g_{\sigma}$ the one-dimensional centered Gaussian density with standard deviation $\sigma$.

which implies $\left|y_{i}+\gamma\ \textrm{sign}(y_{i})-\mu_{i}(x)\right|^{2}\geq\frac{1}{2}\left|y_{i}-x_{i}\right|^{2}-\left(\gamma+D\sigma^{2}/2\right)^{2}$. Similarly, $\left|y_{i}+\gamma\ \textrm{sign}(y_{i})-\mu_{i}(x)\right|^{2}\leq 2\left|y_{i}-x_{i}\right|^{2}+2\left(\gamma+D\sigma^{2}/2\right)^{2}$.

The proof of Theorem 3.1 also requires a lower bound on the probability that a component of the proposed point will be set to zero. Such a bound is given in Lemma 6.3.

Let $\rho$ and $\mu$ be given by Lemma 2.2 and (8). It holds

For $i\not\in I_{m}$, by (8), $|\mu_{i}(z)|\leq D\sigma^{2}/2$. Hence, there exists a constant $C>0$ such that

The Markov kernel $P_{\Psi}$ is psi-irreducible and aperiodic.

so that it is enough to establish a minorization on the kernel for any $x\in C\cap S_{m_{\star}}$ whatever $m_{\star}\in\mathcal{M}$. Let $m_{\star}\in\mathcal{M}$. By definition of $P$, $q$ (see (13)) and $\nu$

where, for any $x\in S_{m_{\star}}$ and $y\in S_{m}$, we have

where the last inequality follows from Lemma 6.1. For any $x\in S_{m_{\star}}$ and $y\in S_{m}$, we have

For any $m\in\mathcal{M}$, $\limsup\limits_{\|x\|\to\infty}T_{m}(x)=0$.

The proof is adapted from and . Let $m\in\mathcal{M}$ be fixed. Define

Moreover, by definition of $A_{m}(x)$, for any $z\in A_{m}(x)$ it holds

This yields there exists $C_{\star}$ such that for any $a,u,r>0$,

Let $z\in\mathcal{C}_{m}^{c}(x,u)\cap\{z:\pi(z^{[m]})<\pi(x)\}$. By AA1(ii), $h:s\mapsto\pi(z^{[m]}-s\ n(z^{[m]}))-\pi(x)$ is continuous, and by definition of $\mathcal{C}_{m}^{c}(x,u)$, $h(s)\neq 0$ for any $0\leq s\leq u$. Since $h(0)<0$ (we assumed that $\pi(z^{[m]})<\pi(x)$), this implies that $h(u)<0$ i.e. $\pi(z^{[m]}-sn(z^{[m]}))\leq\pi(x)$. Then,

If $z\in\mathcal{C}_{m}^{c}(x,u)\cap\{z:\pi(z^{[m]})\geq\pi(x)\}$, we obtain similarly that $\pi(x)/\pi(z^{[m]})\leq d_{r}(u)$. Hence, we established that

As a conclusion, there exists $C_{\star}>0$ and for any $\epsilon,a,u>0$, there exists $M>0$ such that $\sup_{\|x\|\geq M}T_{m,3}(x,a,u)\leq C_{\star}\epsilon$.

Following the same lines as for the control of $T_{m,3}(x,a,u)$, it can be shown that there exists $C_{\star}>0$ and for any $\epsilon,a,u>0$, there exists $M>0$ such that $\sup_{\|x\|\geq M}T_{m,4}(x,a,u)\leq C_{\star}\epsilon$. ∎

Let $u,b,\epsilon,R$ be given by AA3 and $W_{m}(x)$ be defined by (11). There exists $r>R$ such that for any $m\in\mathcal{M}$ and $x\in S_{m}\cap\{\|x\|\geq r\}$, $W_{m}(x)\subset\{y\in S_{m},\alpha_{\Psi}(x,y)=1\}$.

The proof is adapted from . Let $m\in\mathcal{M}$ and $x\in S_{m}$ such that $\|x\|\geq r$ for some $r>R$ to be fixed later (the constant $R$ is given by AA3). We first prove that there exists a positive constant $C_{b}$ such that

By (13), Lemma 6.1 and Lemma 6.3, there exist $C,C_{b}>0$ - independent of $x\in S_{m}$ - such that

By AA2, we can choose $r$ large enough so that for all $\|x\|\geq r$, $\pi(x)/\pi(x-un(x))\leq C_{b}$. This yields (18). Let $z\in W_{m}(x)$. Then, $\|z-x\|\leq b$ so that $z\in\mathcal{B}_{m}(x,b)$. Hence, by (18), $q(z^{[m]},x)/q(x,z^{[m]})\geq C_{b}$. In addition,

where in the last inequality we used AA3. Hence,

and $\alpha_{\Psi}(x,z^{[m]})=1$ thus showing the lemma. ∎

By Lemma 6.6, for any $x\in S_{m_{\star}}$ large enough,

so that the integrals in (19) depend on $x$ only through $m_{\star}$. Since $\mathcal{M}$ is finite, there exists a constant $C^{\prime}>0$ independent of $x$ such that for any $m\in\mathcal{M}$ and $x\in S_{m}$,

The result follows from (15) and Lemmas 6.5 and 6.7. ∎

References