In high-dimensional statistics, variable selection is an optimization problem
aiming to recover the latent sparse pattern from all possible covariate
combinations. In this paper, we propose a novel optimization method to solve
the exact L0​-regularized regression problem (a.k.a. best subset selection).
We reformulate the optimization problem from a discrete space to a continuous
one via probabilistic reparameterization. Within the framework of stochastic
gradient descent, we propose a family of unbiased gradient estimators to
optimize the L0​-regularized objective and a variational lower bound. Within
this family, we identify the estimator with a non-vanishing signal-to-noise
ratio and uniformly minimum variance. Theoretically, we study the general
conditions under which the method is guaranteed to converge to the ground truth
in expectation. In a wide variety of synthetic and semi-synthetic data sets,
the proposed method outperforms existing variable selection methods that are
based on penalized regression and mixed-integer optimization, in both sparse
pattern recovery and out-of-sample prediction. Our method can find the true
regression model from thousands of covariates in a couple of seconds.