Bayesian $l_0$-regularized Least Squares

TL;DR

This paper introduces a scalable method for Bayesian variable selection using $l_0$ regularization, linking Bayesian regularization with proximal updates, and demonstrates its effectiveness through simulations and real data examples.

Contribution

It establishes a connection between Bayesian regularization and proximal updates, enabling efficient spike-and-slab variable selection with the SBR algorithm.

Findings

SBR is faster and more scalable than traditional posterior sampling.

The method achieves comparable statistical properties to spike-and-slab priors.

Simulation and real data results validate the approach's efficiency and accuracy.

Abstract

Bayesian $l_0$-regularized least squares is a variable selection technique for high dimensional predictors. The challenge is optimizing a non-convex objective function via search over model space consisting of all possible predictor combinations. Spike-and-slab (a.k.a. Bernoulli-Gaussian) priors are the gold standard for Bayesian variable selection, with a caveat of computational speed and scalability. Single Best Replacement (SBR) provides a fast scalable alternative. We provide a link between Bayesian regularization and proximal updating, which provides an equivalence between finding a posterior mode and a posterior mean with a different regularization prior. This allows us to use SBR to find the spike-and-slab estimator. To illustrate our methodology, we provide simulation evidence and a real data example on the statistical properties and computational efficiency of SBR versus direct posterior sampling using spike-and-slab priors. Finally, we conclude with directions for future research.

Tables

Table 1. Table 1: Variable selection for the diabetes data. Out of all 64 64 64 variable, only those selected by at least one method are shown, chosen variables marked by ✓ ✓ \checkmark , and — otherwise. Lasso and elastic net select the most variables, whereas spike-and-slab the least. Compared with these two extremes, SBR selects all variables chosen by spike-and-slab, but no variables not chosen by Lasso and elastic net. The result indicates that SBR performs a reasonable variable selection task.

Variable	Lasso	Elastic Net	Spike & Slab	SBR
sex	$✓$	$✓$	$✓$	$✓$
bmi	$✓$	$✓$	$✓$	$✓$
map	$✓$	$✓$	$✓$	$✓$
hdl	$✓$	$✓$	$✓$	$✓$
ltg	$✓$	$✓$	$✓$	$✓$
glu	$✓$	$✓$	—	—
${\text{age}}^2$	$✓$	$✓$	—	—
${\text{bmi}}^2$	$✓$	$✓$	—	—
${\text{glu}}^2$	$✓$	$✓$	—	$✓$
$\text{age}\cdot \text{sex}$	$✓$	$✓$	$✓$	$✓$
$\text{age}\cdot \text{map}$	$✓$	$✓$	—	—
$\text{age}\cdot \text{ltg}$	$✓$	$✓$	—	—
$\text{age}\cdot \text{glu}$	$✓$	$✓$	—	—
$\text{sex}\cdot \text{map}$	—	$✓$	—	—
$\text{bmi}\cdot \text{map}$	$✓$	$✓$	—	$✓$

Equations

$$\underset{\beta \in {\mathbb{R}}^p}{\mathrm{minimize}}$$

$$p(y\mid \beta )\propto \exp \{-l(\beta )\},p(\beta )\propto \exp \{-{\text{pen}}_{\lambda }\left(\beta \right)\}.$$

$$y=X\beta +e,\text{where }e\sim N(0,{\sigma }_e^2{I}_n).$$

$${\beta }_i|\theta ,{\sigma }_{\beta }^2\sim (1-\theta ){\delta }_0+\theta N\left(0,{\sigma }_{\beta }^2\right).$$

$$\begin{array}[]{rcl}\gamma_{i}|\theta&\sim&\text{Bernoulli}(\theta)\ ;\\ \alpha_{i}|\sigma_{\beta}^{2}&\sim&N\left(0,\sigma_{\beta}^{2}\right)\ .\\ \end{array}$$

$$p\left({\gamma }_i,{\alpha }_i\mid \theta ,{\sigma }_{\beta }^2\right)={\theta }^{{\gamma }_i}{\left(1-\theta \right)}^{1-{\gamma }_i}\frac{1}{\sqrt{2\pi }{\sigma }_{\beta }}\exp \left\{-\frac{{\alpha }_i^2}{2{\sigma }_{\beta }^2}\right\},\text{for }1\le i\le p.$$

$$\begin{array}[]{rcl}p\left(\gamma,\alpha\mid\theta,\sigma_{\beta}^{2}\right)&=&\prod\limits_{i=1}^{p}p\left(\alpha_{i},\gamma_{i}\mid\theta,\sigma_{\beta}^{2}\right)\\ &=&\theta^{\|\gamma\|_{0}}\left(1-\theta\right)^{p-\|\gamma\|_{0}}\left(2\pi\sigma_{\beta}^{2}\right)^{-\frac{p}{2}}\exp\left\{-\frac{1}{2\sigma_{\beta}^{2}}\sum\limits_{i=1}^{p}\alpha_{i}^{2}\right\}\ .\end{array}$$

$$p\left(y\mid \gamma ,\alpha ,\theta ,{\sigma }_e^2\right)={\left(2\pi {\sigma }_e^2\right)}^{-\frac{n}{2}}\exp \left\{-\frac{1}{2{\sigma }_e^2}{\Vert y-X_{\gamma }{\alpha }_{\gamma }\Vert }_2^2\right\}.$$

$$\begin{array}[]{rcl}p\left(\gamma,\alpha\mid\theta,\sigma_{\beta}^{2},\sigma_{e}^{2},y\right)&\propto&p\left(\gamma,\alpha\mid\theta,\sigma_{\beta}^{2}\right)p\left(y\mid\gamma,\alpha,\theta,\sigma_{e}^{2}\right)\\ &\propto&\exp\left\{-\frac{1}{2\sigma_{e}^{2}}\left\|y-X_{\gamma}\alpha_{\gamma}\right\|_{2}^{2}-\frac{1}{2\sigma_{\beta}^{2}}\left\|\alpha\right\|_{2}^{2}-\log\left(\frac{1-\theta}{\theta}\right)\left\|\gamma\right\|_{0}\right\}\ .\end{array}$$

$$\arg \underset{\gamma ,\alpha }{\max }p\left(\gamma ,\alpha \mid \theta ,{\sigma }_{\beta }^2,{\sigma }_e^2,y\right).$$

$$\underset{\gamma ,\alpha }{\min }{\Vert y-X_{\gamma }{\alpha }_{\gamma }\Vert }_2^2+\frac{{\sigma }_e^2}{{\sigma }_{\beta }^2}{\Vert \alpha \Vert }_2^2+2{\sigma }_e^2\log \left(\frac{1-\theta }{\theta }\right){\Vert \gamma \Vert }_0.$$

$$\underset{\beta }{\min }\frac{1}{2}{\Vert y-X\beta \Vert }_2^2+\lambda {\Vert \beta \Vert }_0.$$

$$\theta <\frac{1}{2},{\sigma }_{\beta }^2\gg {\sigma }_e^2,\frac{{\sigma }_e^2}{{\sigma }_{\beta }^2}{\Vert \alpha \Vert }_2^2\to 0,$$

$$\underset{\gamma ,\alpha }{\min }\frac{1}{2}{\Vert y-X_{\gamma }{\alpha }_{\gamma }\Vert }_2^2+\lambda {\Vert \gamma \Vert }_0,\text{where }\lambda :={\sigma }_e^2\log \left(\left(1-\theta \right)/\theta \right)>0.$$

$$\begin{array}[]{rcl}y|\beta&\sim&N(\beta,\sigma_{e}^{2})\ ;\\ \beta&\sim&p\left(\beta\right)\ .\end{array}$$

$$\hat{\beta }=E\left[\beta \mid y\right]=y+{\sigma }_e^2\frac{d}{dy}\log m(y),$$

$$m(y)=\int f\left(y\mid \beta \right)p(\beta )d\beta .$$

$${\sigma }_e^2\frac{d}{dy}\log m(y),$$

$$\hat{\beta }=y+{\hat{\sigma }}_e^2\frac{d}{dy}\log \hat{m}(y).$$

$$R(\hat{\beta })={\sigma }_e^2\left(1-{\sigma }_e^{2}I\left(m\right)\right),$$

$$y=X\beta +e,\text{where }e\sim N\left(0,\Sigma \right).$$

$$E\left[\beta \mid y\right]={\left(X^T{\Sigma }^{-1}X\right)}^{-1}{X}^T\left({\Sigma }^{-1}y+{\nabla }_y\log m\left(y\right)\right).$$

$$p\left(\beta \mid y\right)=\frac{p\left(y\mid \beta \right)p\left(\beta \right)}{m\left(y\right)}=\frac{1}{m\left(y\right)}N\left(y;X\beta ,\Sigma \right)p\left(\beta \right).$$

$$\begin{array}[]{rcl}E\left[\Sigma^{-1}\left(y-X\beta\right)\mid y\right]&=&\int_{\beta}\Sigma^{-1}\left(y-X\beta\right)p\left(\beta\mid y\right)d\beta\\ &=&\frac{1}{m\left(y\right)}\int_{\beta}\Sigma^{-1}\left(y-X\beta\right)N\left(y;X\beta,\Sigma\right)p\left(\beta\right)d\beta\ .\end{array}$$

$${\Sigma }^{-1}\left(y-X\beta \right)N\left(y;X\beta ,\Sigma \right)=-{\nabla }_{y}N\left(y;X\beta ,\Sigma \right),$$

$$E\left[{\Sigma }^{-1}\left(y-X\beta \right)\mid y\right]=\frac{1}{m\left(y\right)}{\int }_{\beta }-{\nabla }_{y}N\left(y;X\beta ,\Sigma \right)p\left(\beta \right)d\beta =-{\nabla }_y\log m\left(y\right).$$

$$\hat{\beta }=E\left[\beta \mid y\right]=\arg \underset{\beta }{\min }\left\{\frac{1}{2{\sigma }_e^2}(y-\beta {)}^2+\varphi (\beta )\right\}.$$

$$\begin{array}[]{rcl}E_{\gamma f}(x)&=&\inf_{z}\left\{f(z)+\frac{1}{2\gamma}\|z-x\|_{2}^{2}\right\}\leq f(x)\ ;\\ \mathop{\mathrm{prox}}_{\gamma f}(x)&=&\arg\min_{z}\left\{f(z)+\frac{1}{2\gamma}\|z-x\|_{2}^{2}\right\}\,.\end{array}$$

$$\nabla \left\{f(\hat{z})+\frac{1}{2\gamma }\Vert \hat{z}-x{\Vert }_2^2\right\}=\nabla f(\hat{z})+\frac{1}{\gamma }(\hat{z}-x)=0,$$

$$\nabla E_{\gamma f}(x)=\nabla \underset{z}{\inf }\left\{f(z)+\frac{1}{2\gamma }\Vert z-x{\Vert }_2^2\right\}=\frac{1}{\gamma }[x-\hat{z}(x)],$$

Full text

Keywords: Spike-and-Slab Prior, $l_{0}$ Regularization, Proximal Updating, Single Best Replacement, Lasso, Variable Selection, Sparsity, Bayes, Model Choice

1 Introduction

Bayesian $l_{0}$ regularization is an attractive solution for high dimensional variable selection as it directly penalizes the number of predictors. The caveat is the need to search over all possible model combinations, as a full solution requires enumeration over all possible models which is NP-hard. The gold standard for Bayesian variable selection are spike-and-slab priors, or Bernoulli-Gaussian mixtures (Mitchell and Beauchamp, 1988; George and McCulloch, 1993; Scott and Varian, 2014). Whilst spike-and-slab priors provide full model uncertainty quantification, they can be hard to scale to very high dimensional problems, and can have poor sparsity properties (Amini et al., 2012). On the other hand, techniques like proximal algorithms (Polson et al., 2015; Polson and Sun, 2017) can solve non-convex optimization problems which are fast and scalable, but they generally don’t provide a full assessment of model uncertainty (Jeffreys, 1961; Hans, 2009; Scott and Berger, 2010; Li and Lin, 2010; Marjanovic and Solo, 2013).

Our goal is to build on the single best replacement (SBR) algorithm of Soussen et al. (2011) as one approach to Bayesian $l_{0}$ regularization, in contrast with other methods based on proximal algorithms (Parikh and Boyd, 2014; Polson et al., 2015; Polson and Sun, 2017). Our approach also builds on other Bayesian regularization methods including, for example, the Bayesian bridge (Polson et al., 2014), horseshoe regularization (Carvalho et al., 2010; Bhadra et al., 2017a), SVMs (Polson and Scott, 2011), Bayesian lasso (Hans, 2009; Park and Casella, 2008; Carlin and Polson, 1991), Bayesian elastic net (Li and Lin, 2010; Hans, 2011), spike-and-slab lasso (Rockova and George, 2017), and global-local shrinkage priors (Bhadra et al., 2017a; Griffin and Brown, 2010).

To fix notation, statistical regularization requires the specification of a measure of fit, denoted by $l\left(\beta\right)$ and a penalty function, denoted by $\text{pen}_{\lambda}\left(\beta\right)$, where $\lambda$ is a global regularization parameter. From a Bayesian perspective, $l\left(\beta\right)$ and $\text{pen}_{\lambda}\left(\beta\right)$ correspond to the negative logarithms of the likelihood and prior distribution, respectively. Regularization leads to an optimization problem of the form

[TABLE]

Taking a probabilistic approach leads to a Bayesian hierarchical model

[TABLE]

The solution to the minimization problem estimated by regularization corresponds to the posterior mode, $\hat{\beta}={\rm arg\;max}_{\beta}\;p(\beta|y)$, where $p(\beta|y)$ denotes the posterior distribution. For example, regression with a least squares log-likelihood subject to a penalty such as an $l_{2}$-norm (ridge) Gaussian probability model or $l_{1}$-norm (lasso) double exponential probability model.

The rest of the paper is outlined as follows. Section 2 defines the Bayesian $l_{0}$ regularization problem and explores the connections with spike-and-slab priors. Section 3 introduces a novel connection between regularization and Bayesian inference. Section 4 surveys recent developments on $l_{0}$ regularization and best subset selection problems. Section 5 describes the SBR algorithm for $l_{0}$ regularization. Section 6 provides a comparison between SBR and other variable selection methods and models including Lasso and elastic net (Tibshirani, 1994; Zou and Hastie, 2005), Bayesian bridge (Polson et al., 2014), and spike-and-slab (George and McCulloch, 1993). Finally, Section 7 concludes with directions for future research.

2 Bayesian $l_{0}$ regularization

Consider a standard Gaussian linear regression model, where $X=[X_{1},\ldots,X_{p}]\in\mathbb{R}^{n\times p}$ is a design matrix, $\beta=(\beta_{1},\ldots,\beta_{p})^{T}\in\mathbb{R}^{p}$ is a $p$-dimensional coefficient vector, and $e$ is an $n$-dimensional independent Gaussian noise. After centralizing $y$ and all columns of $X$, we ignore the intercept term in the design matrix $X$ as well as $\beta$, and we can write

[TABLE]

To specify a prior distribution $p\left(\beta\right)$, we impose a sparsity assumption on $\beta$, where only a small portion of all $\beta_{i}$’s are non-zero. In other words, $\|\beta\|_{0}=k\ll p$, where $\|\beta\|_{0}\mathrel{\mathop{:}}=\#\{i:\beta_{i}\neq 0\}$, the cardinality of the support of $\beta$, also known as the $l_{0}$ pseudo-norm of $\beta$. A multivariate Gaussian prior ($l_{2}$ norm) leads to poor sparsity properties in this situation (see, e.g., Polson and Scott, 2010).

Sparsity-inducing prior distributions for $\beta$ can be constructed to impose sparsity. The gold standard is a spike-and-slab priors (Jeffreys, 1961; Mitchell and Beauchamp, 1988; George and McCulloch, 1993). Under these assumptions, each $\beta_{i}$ exchangeably follows a mixture prior consisting of $\delta_{0}$, a point mass at [math], and a Gaussian distribution centered at zero. Hence we write,

[TABLE]

Here $\theta\in\left(0,1\right)$ controls the overall sparsity in $\beta$ and $\sigma_{\beta}^{2}$ accommodates non-zero signals. This family is termed as the Bernoulli-Gaussian mixture model in the signal processing community.

A useful re-parameterization, the parameters $\beta$ is given by two independent random variable vectors $\gamma=\left(\gamma_{1},\ldots,\gamma_{p}\right)^{\prime}$ and $\alpha=\left(\alpha_{1},\ldots,\alpha_{p}\right)^{\prime}$ such that $\beta_{i}=\gamma_{i}\alpha_{i}$, with probabilistic structure

[TABLE]

Since $\gamma_{i}$ and $\alpha_{i}$ are independent, the joint prior density becomes

[TABLE]

The indicator $\gamma_{i}\in\{0,1\}$ can be viewed as a dummy variable to indicate whether $\beta_{i}$ is included in the model. (Soussen et al., 2011)

Let $S=\{i:\gamma_{i}=1\}\subseteq\{1,\ldots,p\}$ be the “active set” of $\gamma$, and $\|\gamma\|_{0}=\sum\limits_{i=1}^{p}\gamma_{i}$ be its cardinality. The joint prior on the vector $\{\gamma,\alpha\}$ then factorizes as

[TABLE]

Let $X_{\gamma}\mathrel{\mathop{:}}=\left[X_{i}\right]_{i\in S}$ be the set of “active explanatory variables” and $\alpha_{\gamma}\mathrel{\mathop{:}}=\left(\alpha_{i}\right)^{\prime}_{i\in S}$ be their corresponding coefficients. We can write $X\beta=X_{\gamma}\alpha_{\gamma}$. The likelihood can be expressed in terms of $\gamma$, $\alpha$ as

[TABLE]

Under this re-parameterization by $\left\{\gamma,\alpha\right\}$, the posterior is given by

[TABLE]

Our goal then is to find the regularized maximum a posterior (MAP) estimator

[TABLE]

By construction, the $\gamma$ $\in\left\{0,1\right\}^{p}$ will directly perform variable selection. Spike-and-slab priors, on the other hand, will sample the full posterior and calculate the posterior probability of variable inclusion.

Finding the MAP estimator is equivalent to minimizing over $\left\{\gamma,\alpha\right\}$ the regularized least squares objective function (Soussen et al., 2011).

[TABLE]

This objective possesses several interesting properties:

The first term is essentially the least squares loss function. 2. 2.

The second term looks like a ridge regression penalty and has connection with the signal-to-noise ratio (SNR) $\sigma_{\beta}^{2}/\sigma_{e}^{2}$. Smaller SNR will be more likely to shrink the estimate of $\alpha$ towards [math]. If $\sigma_{\beta}^{2}\gg\sigma_{e}^{2}$, the prior uncertainty on the size of non-zero coefficients is much larger than the noise level, that is, the SNR is sufficiently large, this term can be ignored. This is a common assumption in spike-and-slab framework in that people usually want $\sigma_{\beta}\to\infty$ or to be “sufficiently large” in order to avoid imposing harsh shrinkage to non-zero signals. 3. 3.

If we further assume that $\theta<\frac{1}{2}$, meaning that the coefficients are known to be sparse a priori, then $\log\left(\left(1-\theta\right)/\theta\right)>0$, and the third term can be seen as an $l_{0}$ regularization.

Therefore, our Bayesian objective inference is connected to $l_{0}$-regularized least squares, which we summarize in the following proposition.

Proposition 1**.**

(Spike-and-slab MAP & $l_{0}$ regularization)

For some $\lambda>0$, assuming $\theta<\frac{1}{2}$, $\sigma_{\beta}^{2}\gg\sigma_{e}^{2}$, the Bayesian MAP estimate defined by (2.4) is equivalent to the $l_{0}$ regularized least squares objective, for some $\lambda>0$,

[TABLE]

Proof.

First, assuming that

[TABLE]

gives us an objective function of the form

[TABLE]

Equation (2.6) can be seen as a variable selection version of equation (2.5). The interesting fact is that (2.5) and (2.6) are equivalent. To show this, we need only to check that the optimal solution to (2.5) corresponds to a feasible solution to (2.6) and vice versa. This is explained as follows.

On the one hand, assuming $\hat{\beta}$ is an optimal solution to (2.5), then we can correspondingly define $\hat{\gamma}_{i}\mathrel{\mathop{:}}=I\left\{\hat{\beta}_{i}\neq 0\right\}$, $\hat{\alpha}_{i}\mathrel{\mathop{:}}=\hat{\beta}_{i}$, such that $\left\{\hat{\gamma},\hat{\alpha}\right\}$ is feasible to (2.6) and gives the same objective value as $\hat{\beta}$ gives (2.5).

On the other hand, assuming $\left\{\hat{\gamma},\hat{\alpha}\right\}$ is optimal to (2.6), implies that we must have all of the elements in $\hat{\alpha}_{\gamma}$ should be non-zero, otherwise a new $\tilde{\gamma}_{i}\mathrel{\mathop{:}}=I\left\{\hat{\alpha}_{i}\neq 0\right\}$ will give a lower objective value of (2.6). As a result, if we define $\hat{\beta}_{i}\mathrel{\mathop{:}}=\hat{\gamma}_{i}\hat{\alpha}_{i}$, $\hat{\beta}$ will be feasible to (2.5) and gives the same objective value as $\left\{\hat{\gamma},\hat{\alpha}\right\}$ gives (2.6).

Combining both arguments shows that the two problems (2.5) and (2.6) are equivalent. Hence we can use results from non-convex optimization literature to find Bayes MAP estimators. ∎

3 Bayesian regularization and Proximal Updating

Section 2 provides a connection between spike-and-slab and $l_{0}$ regularization, in the sense that $l_{0}$ regularization can be viewed as the MAP estimator of Bayesian spike-and-slab regression. From a Bayesian perspective, the posterior mean is preferred to the MAP estimator due to its superior risk minimization properties under mean squared error. Starck et al. (2013) also discusses the inadequacy of interpreting sparse regularization as Bayesian MAP estimator. We now show that the posterior mean estimation is also connected to regularization and introduce the connection with the help of proximal operators and Tweedie’s formula (Efron, 2011).

Gribonval (2011) considers the relationship between the MAP estimator and the posterior mean. Polson and Scott (2016) considers the relationship between posterior modes and envelopes. These approaches try to uncover the implicit prior that maps a posterior mean to a mode. As the posterior mean is designed to minimize the mean squared error and Bayes risk, this is the appropriate calculation.

3.1 Posterior mean optimality

Consider a normal mean problem in Bayesian setting,

[TABLE]

Then the optimal estimator of $\beta$ with respect to the quadratic loss is the posterior mean. To calculate this posterior mean, $\hat{\beta}$, Efron (2011) introduced Tweedie’s formula which, for the normal mean problem, gives

[TABLE]

where the marginal density of $y$ is

[TABLE]

Here $f\left(y\mid\beta\right)=\frac{1}{\sqrt{2\pi}\sigma_{e}}\exp\left\{-\frac{(y-\beta)^{2}}{2\sigma_{e}^{2}}\right\}$ is the probability density function (pdf) of $N(\beta,\sigma_{e}^{2})$. The interpretation of the Bayesian correction term,

[TABLE]

is to “regularize” the unbiased maximum likelihood estimator $y$, which provides the optimal bias-variance trade-off for prediction, see Pericchi and Smith (1992) for further discussion. When both $m$ and $\sigma_{e}^{2}$ are unknown, Donoho and Reeves (2013) proposes a procedure to estimate each term, which achieves the optimal Bayes risk, resulting in the plug-in estimator

[TABLE]

Another useful result applies Stein’s risk function to Tweedie’s formula. Then we can derive this optimal Bayes risk (Robbins, 1956),

[TABLE]

where $I(m)=E_{y}\left[\left(\frac{d}{dy}\log m(y)\right)^{2}\right]$ is the Fisher Information for $m$. This risk is optimal given the use of the posterior mean estimator.

Tweedie’s formula can be generalized to Gaussian linear regression in the following theorem.

Theorem 3.1**.**

Suppose in the Gaussian linear regression model,

[TABLE]

Let $p\left(\beta\right)$ denote the prior density of $\beta$, and $m\left(y\right)=\int_{\beta}p\left(y\mid\beta\right)p\left(\beta\right)d\beta$ the marginal (prior predictive) density of $y$. Then, assuming $\left(X^{T}\Sigma^{-1}X\right)^{-1}$ exists, the posterior mean of $\beta$ given $y$ is

[TABLE]

Proof.

Let $N\left(y;X\beta,\Sigma\right)$ denote the multivariate normal density of $y|\beta\sim N\left(X\beta,\Sigma\right)$, then the posterior density of $\beta$ given $y$,

[TABLE]

Therefore, the posterior mean of the quantity $\Sigma^{-1}\left(y-X\beta\right)$,

[TABLE]

Note that by the property of the multivariate normal density,

[TABLE]

and so (3.4) becomes

[TABLE]

Multiplying both sides by $X^{T}$ and assuming $\left(X^{T}\Sigma^{-1}X\right)^{-1}$ exists, the posterior mean of $\beta$ given $y$ becomes

$\displaystyle E\left[\beta\mid y\right]=\left(X^{T}\Sigma^{-1}X\right)^{-1}X^{T}\left(\Sigma^{-1}y+\nabla_{y}\log m\left(y\right)\right)\ .$ ∎

It’s easy to see that, similar to Tweedie’s formula for the normal means problem, the posterior mean in the Gaussian linear regression (3.3) consists of two parts. One is the usual weighted least squares solution, and the other is a Bayesian correction by the gradient of the prior predictive score, $\nabla_{y}\log m\left(y\right)$. Griffin and Brown (2010) gives the equivalent result in the form when the least squares estimator $\hat{\beta}$ instead of $y$ is conditioned on. Masreliez (1975) discusses the posterior mean under Gaussian prior but non-Gaussian likelihood. (Pericchi and Smith, 1992)

3.2 Regularized linear regression

Proximal operators and Tweedie’s formula also provide a way to connect the Bayesian posterior mean and the regularized least squares. Specifically, we want to find a $\phi$, such that the Bayesian posterior estimator $\hat{\beta}=E\left[\beta\mid y\right]$ with the prior $p$ is the same as the solution to the $\phi$-regularized least squares. That is,

[TABLE]

We now use the theory of proximal mappings to re-write this estimator. First, here are some definitions. The Moreau envelope $E_{\gamma f}(x)$ and proximal mapping $\mathop{\mathrm{prox}}_{\gamma f}(x)$ of a convex function are defined as

[TABLE]

The Moreau envelope is a regularized version of $f$ and approximates $f$ from below, and has the same set of minimizing values as $f$. The proximal mapping returns the value that solves the minimization problem defined by the Moreau envelope. It balances two goals: minimizing $f$, and staying near $x$.

Now, observe that if $\hat{z}(x)=\mathop{\mathrm{prox}}_{\gamma f}(x)$ is the value that achieves the minimum,

[TABLE]

which leads to $\hat{z}=x-\gamma\nabla f(\hat{z})\,.$ By construction of the envelope,

[TABLE]

This leads to the fundamental proximal relation $\hat{z}=x-\gamma\nabla E_{\gamma f}(x)\,.$ Therefore, write

[TABLE]

Meanwhile, the definition of proximal mapping (3.6) and its property (3.7) give us

[TABLE]

Combining with Tweedie’s formula (3.2), gives

[TABLE]

Hence, if we want to match a regularized least squares with a posterior mean, we can “solve” for the penalty $\phi$, given a marginal distribution $m(y)$, via the equation for the proximal mapping (3.8). If $\phi$ is non-differentiable at a point, we replace $\nabla$ by $\partial$.

Gribonval (2011) provides the following answer. Given any $z$, find $\hat{y}$ such that $E\left[\beta\mid\hat{y}\right]=z$. Then the penalty

[TABLE]

with the constant $c$ to ensure that $\phi(0)=0$. To see why this construction makes sense, simply take derivatives with respect to $\hat{y}$ on both sides, and get

[TABLE]

To summarize, the solution to the regularized least squares problem

[TABLE]

is the posterior mode with the prior $p\left(\beta\right)\propto\exp\left(-\phi\left(\beta\right)\right)$. Our proximal operators and Tweedie’s formula discussion shows that the regularized least squares solution can also be viewed as the posterior mean under an implied prior $p\left(\beta\right)$, see Strawderman et al. (2013).

To illustrate our result, When $p$ is sparsity-inducing, such as the spike-and-slab, we can construct the associated penalty $\phi$, which is typically non-convex for both Gaussian and Laplace cases.

3.3 Example: Spike-and-slab Gaussian & Laplace prior

For the normal mean problem (3.1), assuming $p\left(\beta\right)$ is the aforementioned spike-and-slab (Bernoulli Gaussian) prior (2.2), the marginal distribution of $y$ is a mixture of two mean zero normals,

[TABLE]

The posterior mean $E\left[\beta|y\right]$ is given by

[TABLE]

Thus, $\forall z\in\mathbb{R}$, we can find $\hat{y}$ such that $w\left(\hat{y}\right)\hat{y}=z$; then the penalty $\phi^{BG}$ associated with the Bayesian posterior mean with the Bernoulli-Gaussian prior can be obtained by (3.9). $\phi^{BG}$ doesn’t have an analytical form, but can be computed numerically.

Amini et al. (2012) argued that Bernoulli-Gaussian priors are usually not applicable to real-world signals, and proposed Bernoulli-Laplace priors which are infinitely divisible and more appropriate for sparse signal processing. The Bernoulli-Laplace priors are very similar to the Bernoulli-Gaussian ones, their only difference being that the “slab” parts are replaced by Laplace distributions. Therefore, the prior

[TABLE]

Mitchell (1994) and Hans (2009) studied the marginal and posterior distribution with the Laplace prior. With their results, the marginal density of $y$ with the Bernoulli-Laplace prior (3.10) is given by

[TABLE]

where

[TABLE]

Here $\Phi$ is the cumulative distribution function (cdf) of the standard normal. The posterior mean $E\left[\beta|y\right]$ is then given by

[TABLE]

Similarly, we are also able to find the penalty $\phi^{BL}$ associated with this Bernoulli-Laplace prior numerically by (3.9). It’s worth noting that this prior is a special case of the spike-and-slab Lasso prior proposed by Rockova and George (2017). In that paper the authors use a mixture of two Laplace distributions, one of which is very close to $\delta_{0}$ as its variance goes to zero. Both priors are capable of striking a balance between hard-thresholding and soft-thresholding.

For comparison, $\hat{\beta}^{BG}$ and $\hat{\beta}^{BL}$ are plotted in Figure 1; $\phi^{BG}$ and $\phi^{BL}$ are plotted in Figure 2. Both priors shrink small observations towards zero. For large observations, when $\sigma_{\beta}$ is small, Bernoulli-Gaussian, like ridge regression, unnecessarily penalizes large observations too much, whereas Bernoulli-Laplace is more like Lasso. As $\sigma_{\beta}$ gets larger, both priors get closer to hard-thresholding, and their associated penalties $\phi$ closer to SCAD-like non-convex penalties (Fan and Li, 2001).

4 Computing the $l_{0}$-regularized regression solution

We now turn to the problem of computation. $l_{0}$-regularized least squares (2.5) is closely related to the best subset selection in linear regression as follows.

[TABLE]

The $l_{0}$-regularized least squares (2.5) can be seen as (4.1)’s Lagrangian form. However, due to high non-convexity of the $l_{0}$-norm, (2.5) and (4.1) are connected but not equivalent. In particular, for any given $\lambda\geq 0$, there exists an integer $k\geq 0$, such that (2.5) and (4.1) have the same global minimizer $\hat{\beta}$. However, it’s not true the other way around. It’s possible, even common, that for a given $k$, we cannot find a $\lambda\geq 0$, such that the solutions to (4.1) and (2.5) are the same.

Indeed, for $k\in\left\{1,2,\ldots,p\right\}$, let $\hat{\beta}_{k}$ be respective optimal solutions to (4.1) and $f_{k}$ respective optimal objective values, and so $f_{1}\geq f_{2}\geq\cdots\geq f_{p}$. If we want a solution $\hat{\beta}_{\lambda}$ to (2.5) has $\left\|\hat{\beta}_{\lambda}\right\|_{0}=k$, we need to find a $\lambda$ such that

[TABLE]

with the caveat that such $\lambda$ needs not exist.

Both problems involve discrete optimization and have thus been seen as intractable for large-scale data sets. As a result, in the past, $l_{0}$ norm is usually replaced by its convex relaxation $l_{1}$ norm to facilitate computation. However, it’s widely known that the solutions of $l_{0}$ norm problems provide superior variable selection and prediction performance compared with their $l_{1}$ convex relaxation such as Lasso. Zhang et al. (2014) studies the statistical properties of the theoretical solution to (2.5), and points out that the solution to the $l_{0}$-regularized least squares should be better than Lasso in terms of variable selection especially when we have a design matrix $X$ that has high collinearity among its columns.

Bertsimas et al. (2016) introduced a first-order algorithm to provide a stationary solution $\beta^{*}$ to a class of generalized $l_{0}$-constrained optimization problem, with convex $g$,

[TABLE]

Let $L$ be the Lipschitz constant for $\nabla g$ such that $\forall\beta_{1},\beta_{2}$, $\|\nabla g(\beta_{1})-\nabla g(\beta_{2})\|\leq L\|\beta_{1}-\beta_{2}\|$. Their “Algorithm 1” is as follows.

Initialize $\beta^{0}$ such that $\left\|\beta^{0}\right\|_{0}\leq k$. 2. 2.

For $t\geq 0$, obtain $\beta^{t+1}$ as

[TABLE]

until convergence to $\beta^{*}$.

where the operator $H_{k}\left(\cdot\right)$ is to keep the largest $k$ elements of a vector as the same, whilst to set all else to zero. It can also be called the hard thresholding at the $k^{\text{th}}$ largest element. In the least squares setting when $g(\beta)=\frac{1}{2}\|y-X\beta\|_{2}^{2}$, $\nabla g$ and $L$ are easy to compute. Bertsimas et al. (2016) then uses the stationary solution $\beta^{*}$ obtained by the aforementioned algorithm (4.3) as a warm start for their mixed integer optimization (MIO) scheme to produce a “provably optimal solution” to the best subset selection problem (4.1).

It’s worth pointing out that the key iteration step (4.3) is connected to the proximal gradient descent (PGD) algorithm many have used to solve the $l_{0}$-regularized least squares (2.5), as well as other non-convex regularization problems. PGD methods solve a general class of problems such as

[TABLE]

where $g$ is the same as in (4.2), and $\phi$, usually non-convex, is a regularization term. In this framework, in order to obtain a stationary solution $\beta^{*}$, the key iteration step is

[TABLE]

where $\beta^{t}-\frac{1}{L}\nabla g(\beta^{t})$ can be seen as a gradient descent step for $g$ and $\mathop{\mathrm{prox}}_{\lambda\phi}$ is the proximal operator for $\lambda\phi$. In $l_{0}$-regularized least squares, $\lambda\phi\left(\cdot\right)=\lambda\left\|\cdot\right\|_{0}$, and its proximal operator $\mathop{\mathrm{prox}}_{\lambda\|\cdot\|_{0}}$ is just the hard thresholding at $\lambda$. That is, $\mathop{\mathrm{prox}}_{\lambda\|\cdot\|_{0}}$ is to keep the same all elements no less than $\lambda$, whilst to set all else to zero. As a result, the similarity between (4.3) and (4.5) are quite obvious.

In a recent work, Jewell and Witten (2017) proposes an exact algorithm to obtain the global minimum of $l_{0}$-regularized optimization in a computational neuroscience context. Consider the optimization problem

[TABLE]

By exploiting the sequential time series nature of the problem, one can recast the problem as a changepoint detection problem and use the available results in that literature to design a dynamic programming algorithm. Instead of trying to simultaneously find all $i$’s where $c_{i}-\gamma c_{i-1}\neq 0$, the algorithm aims to find them sequentially, from $i=1$ to $i=n$ at each step. In this sense, it can reach a global minimum within $\mathbb{O}\left(n^{2}\right)$. The authors further speed up the algorithm by pruning the set of possible changepoints to at each step of the sequential search, and reduce the expected time cost to $\mathbb{O}\left(n\right)$.

5 Single best replacement (SBR) algorithm

The single best replacement (SBR) algorithm, originally developed by Soussen et al. (2011), provides solution to the variable selection regularization (2.6). Since (2.6) and the $l_{0}$-regularized least squares (2.5) are equivalent, SBR also provides a practical way to give a sufficiently good local optimal solution to the NP-hard $l_{0}$ regularization.

Take a look at the objective (2.6). For any given variable selection indicator $\gamma$, we have an active set $S=\left\{i:\gamma_{i}=1\right\}$, based on which the minimizer $\hat{\alpha}_{\gamma}$ of (2.6) has a closed form. $\hat{\alpha}_{\gamma}$ will set every coefficients outside $S$ to zero, and regress $y$ on $X_{\gamma}$, the variables inside $S$. Therefore, the minimization of the objective function can be determined by $\gamma$ or $S$ alone. Accordingly, the objective function (2.6) can be rewritten as follows.

[TABLE]

The SBR algorithm thus tries to minimize $f_{SBR}(S)$ via choosing the optimal $\hat{S}$.

The algorithm works as follows. Suppose we start as an initial $S$, usually the empty set. At each iteration, SBR aimes to find a “single change of $S$”, that is, a single removal from or adding to $S$ of one element, such that this single change decreases $f_{SBR}(S)$ the most. SBR stops when no such change is available, or in other words, any single change of $\gamma$ or $S$ will only give the same or larger objective value. Therefore, intuitively SBR stops at a local optima of $f_{SBR}(S)$.

SBR is essentially a stepwise greedy variable selection algorithm. At each iteration, both adding and removal are allowed, so this algorithm is one example of the “forward-backward” stepwise procedures. It’s provable that with this feature the algorithm “can escape from some [undesirable] local minimizers” of $f_{SBR}(S)$ (Soussen et al., 2015). Therefore, SBR can solve the $l_{0}$-regularized least squares in a sub-optimal way, providing a satisfactory balance between efficiency and accuracy.

We are now writing out the algorithm more formally. For any currently chosen active set $S$, define a single replacement $S\cdot i,i\in\left\{1,\ldots,p\right\}$ as $S$ adding or removing a single element $i$,

[TABLE]

Then we compare the objective value at current $S$ with all of its single replacements $S\cdot i$, and choose the best one. SBR proceeds as follows.

Step 0: Initialize $S_{0}$. Usually, $S_{0}=\emptyset$. Compute $f_{SBR}(S_{0})$. Set $k=1$.

Step $k$: For every $i\in\{1,\ldots,p\}$, compute $f_{SBR}(S_{k-1}\cdot i)$. Obtain the single best replacement $j\mathrel{\mathop{:}}=\arg\min\limits_{i}f_{\text{SBR}}(S_{k-1}\cdot i)$.

1. 1.

If $f_{SBR}(S_{k-1}\cdot j)\geq f_{SBR}(S_{k-1})$, stop. Report $\hat{S}=S_{k-1}$ as the solution. 2. 2.

Otherwise, set $S_{k}=S_{k-1}\cdot j$, $k=k+1$, and repeat step $k$.

Soussen et al. (2011) shows that SBR always stops within finite steps. With the output $\hat{S}$, the locally optimal solution to the $l_{0}$-regularized least squares $\hat{\beta}$ is just the coefficients of $y$ regressed on $X_{\hat{S}}$ and zero elsewhere.

In order to include both forward and backward steps in the variable selection process, the algorithm needs to compute $f_{SBR}(S_{k-1}\cdot i)$ for every $i$ at every step. Because it involves a one-column update of current design matrix $X_{S_{k-1}}$, this computation can be made very efficient by using the Cholesky decomposition, without explicitly calculating $p$ linear regressions at each step (Soussen et al., 2011). An R package implementation of the algorithm is available upon request.

6 Applications

6.1 Statistical properties of SBR and $l_{0}$ regularization

The design matrix $X$ in this experimentation has $n=120$ rows and $p=100$ columns. In order to impose high collinearity in the columns of $X$, we construct it in the following way.

Construct a $p\times d$ matrix $L$ consisting of $N\left(0,1\right)$ random samples and obtain $\Sigma_{X}=LL^{T}+I_{p}$. If $d\ll p$, $\Sigma_{X}$ will have a low-rank structure. Here we use $d=5$. 2. 2.

Sample each of the $n$ rows of $X$ from the multivariate normal distribution $N_{p}\left(0,\Sigma_{X}\right)$. 3. 3.

Centralize and normalize the columns of $X$ such that each column sums to zero and has unit $l_{2}$ norm.

A design matrix $X$ constructed this way has highly collinear columns. $\beta$ is a highly sparse coefficient vector with $100$ elements, $90$ of which are zero, and the rest $10$ randomly chosen to be $\left\{-5,-4,-3,-2,-1,1,2,3,4,5\right\}$. The noise vector $e$ is sampled from $N\left(0,\sigma_{e}^{2}\right)$. In this setting, the signal-to-noise ratio was defined as

[TABLE]

and $\sigma_{e}$ is determined such that $\text{SNR}=20\text{ dB}$. Finally, let $y=X\beta+e$ be the observation.

Essentially all kinds of regularization methods, including ridge regression, Lasso, and $l_{0}$ regularization, share a common difficulty: to find a suitable regularization parameter $\lambda$. A thorough theoretical treatment on the optimal $\lambda$ hasn’t been established, and in practice cross validation is often used.

Meanwhile, one of the advantages of $l_{0}$ regularization is that, compared with its $l_{1}$ counterpart, the estimates of coefficients are relatively insensitive to $\lambda$ when it’s large enough, as shown in Figure 3. Therefore, we don’t have to worry too much about choosing a particularly optimal $\lambda$.

In a large scale numerical experiment, we conducted $1000$ simulation trials with random $X,y$ generated as above. We compare SBR with the ordinary least squares (OLS), two sparsity regularization methods, Lasso and elastic net with the tuning parameter $\alpha=0.5$, and two Bayesian MCMC methods, Bayesian bridge (Polson et al., 2014) and Bayesian spike-and-slab shrinkage (Scott and Varian, 2014), each with their state-of-the-art implementations. The regularization parameter $\lambda$ in SBR is determined by a 10-fold cross validation. For Lasso and elastic net, the regularization parameter $\lambda$ is chosen by the built-in cross validation in glmnet (Friedman et al., 2010). For Bayesian bridge, we use the default prior in the package BayesBridge (Polson et al., 2012). For spike-and-slab, the prior inclusion probability is set to be $0.5$, and other hyperparameters are the same as in the default setting in the package BoomSpikeSlab (Scott, 2016).

Figure 4 shows the accuracy in estimating $\beta$ for the six methods. For OLS, Lasso, elastic net, and SBR, $\hat{\beta}$ are the minimizers of the corresponding optimization problems, and for Bayesian bridge and spike-and-slab, $\hat{\beta}$ are the posterior means. SBR performs as well as the gold standard Spike-and-slab and better than all else, including the two widely used convex regularization methods.

In terms of variable selection, we compare SBR with the other $3$ sparsity-inducing methods, Lasso, elastic net, and Bayesian spike-and-slab. For spike-and-slab, an coefficient is selected when its posterior inclusion probability is greater than $0.5$. Figure 5 shows that all methods are able to select all the true non-zero coefficients almost all the time. However, when compared in terms of preventing false selection, SBR and spike-and-slab are the best, whereas Lasso and elastic net tend to drastically over-select.

Recently, Hastie et al. (2017) compared Lasso with best subset selection (Bertsimas et al., 2016), forward stepwise selection, and found similar phenomenon. Namely, the performance of best subset selection and forward stepwise selection is overall similar, and both tend to outperform Lasso in the high signal-to-noise setting. Since SBR is a forward-backward stepwise selection algorithm, it’s no surprise then that SBR should give better results than Lasso in such setting.

6.2 Computational efficiency and scalability of SBR

In section 6.1 we’ve shown that the $l_{0}$ regularization has superior statistical properties in terms of both minimizing the estimation risk and selecting correct variables, especially its statistical performance improvement on convex regularization methods such as Lasso and elastic net.

Full Bayesian spike-and-slab performs very well statistically. In order to achieve this good performance, spike-and-slab needs to do a complete MCMC sampling, and this task could take a significant amount of time, especially in high-dimensional settings, whereas regularization methods are usually able to handle large-scale computation efficiently.

In order to compare the computational efficiency of different methods, two sets of experiments, one with $n=120,p=100$, the other $n=300,p=200$, are run, and the results are plotted in Figure 6. SBR, as well as Lasso and elastic net, is almost as efficient as OLS, and only changes proportionally when the size of the problem increases. On the contrary, the two full Bayesian methods, Bayesian bridge and especially spike-and-slab, are costly and scale badly with the problem size. Actually when $n=300,p=200$, it could take as much as $40$ minutes to run even one spike-and-slab MCMC, whereas SBR finishes all $200$ simulation trials under $10$ seconds. When $n$ and $p$ are in thousands, spike-and-slab is computationally intractable.

6.3 Diabetes data

Now we examine the performance of SBR on the classic diabetes data, available in the R package lars (Efron et al., 2004). The design matrix $X$ has $64$ columns, including all $10$ biochemical attributes and certain interactions. Each column of $X$ has been normalized to have zero mean and unit $l_{2}$ norm. The response $y$ is centralized to have zero mean. We compare SBR with sparsity-inducing methods including Lasso, elastic net, and spike-and-slab priors, with the same settings as in Section 6.1, and $\lambda$ determined by cross validation. The results shown in Table 1 indicate that SBR’s performance on variable selection is in line with popular sparse linear regression alternatives.

7 Discussion

Bayesian $l_{0}$ regularization can be solved using a fast and scalable single best replacement (SBR) algorithm. In variable selection, this estimator possesses much of the statistical properties of spike-and-slab priors. We provide theoretical links between the spike-and-slab MAP estimator and $l_{0}$ regularization.

We also explore the connection between regularized MAP estimators and posterior means (Strawderman et al., 2013). Tweedie–Masreliez construction of the posterior mean is re-interpreted as a proximal update rule. This proximal update identity shows how the sparse posterior mode can be viewed as a posterior mean under a suitably defined prior. Bernoulli-Gaussian (BG) and Bernoulli-Laplace (BL) priors are used for illustration. Our approach demonstrates how regularized estimators can have good out-of-sample mean squared error.

In simulated and real data applications, SBR performs favorably compared with popular convex regularization methods such as Lasso and elastic net, as well as full Bayesian sampling methods including Bayesian bridge and spike-and-slab priors.

Recently non-convex feature selection methods for sparse signals estimation have gained increasing attention from the statistical learning community, including the classic SCAD penalty (Fan and Li, 2001), $l_{q}$ penalty (Marjanovic and Solo, 2013), and horseshoe regularization (Bhadra et al., 2017b). There are a number of future directions for research, such as regularized logistic regression (Gramacy and Polson, 2012) and structural sparsity learning (Polson and Sun, 2017). A comprehensive theoretical treatment and empirical comparison on different non-convex regularizations on the trade-off between statistical accuracy and computational efficiency remains open.