$\epsilon_i$ are zero-mean independent random noises (but not necessarily identically distributed). Moreover, we assume that the target vector $\bar{{{\mathbf w}}}$ is sparse. That is, $\bar{k}=\|\bar{{{\mathbf w}}}\|_0$ is small. Here we use the standard notation $${{\mathrm{supp}}}({{\mathbf w}})= \{j: {{\mathbf w}}_j \neq 0\}  \qquad
\|{{\mathbf w}}\|_0 = |{{\mathrm{supp}}}({{\mathbf w}})|$$ for any vector ${{\mathbf w}}\in R^p$.

This paper focuses on the feature selection problem, where we are interested in estimating the set of nonzero coefficients ${{\mathrm{supp}}}(\bar{{{\mathbf w}}})$ (also called support set). Let ${{\mathbf y}}$ denote the vector of $[y_i]$ and $X$ be the $n \times d$ matrix with each row a vector ${{\mathbf x}}_i$. The standard statistical method is subset selection ($L_0$ regularization), which computes the following estimator $$\hat{{{\mathbf w}}}_{L_0} = \arg\min_{{{\mathbf w}}\in R^p} \|X {{\mathbf w}}- {{\mathbf y}}\|_2^2 
\qquad \text{subject to } \|{{\mathbf w}}\|_0 \leq k , \label{eq:L0}$$ where $k$ is a tuning parameter. This method is arguably a natural method for feature selection because if noise $\epsilon_i$ are iid Gaussian random variables, then (\[eq:L0\]) can be regarded as a Bayes procedure with an appropriately defined sparse prior over ${{\mathbf w}}$. However, because the optimization problem in (\[eq:L0\]) is nonconvex, the global solution of this problem cannot be efficiently computed. In practice, 