Revisiting Bayesian Variable Selection via Optimization

1. Problem Statement

The paper addresses a fundamental enigma in Bayesian variable selection for linear regression: while methods using discrete or continuous shrinkage priors (e.g., spike-and-slab, horseshoe) are empirically successful, the underlying marginal likelihood (obtained by integrating out regression coefficients $\theta$ and noise variance $\sigma^2$) is not log-concave.

• The Challenge: Because the negative log-marginal likelihood is non-convex, standard optimization techniques offer no guarantee of finding the global optimum. This raises doubts about the reliability of Type-II Maximum Likelihood (MLE) or Maximum Marginal Posterior estimates for the precision parameters ($D$), which control variable selection.
• The Gap: Existing literature relies heavily on Markov Chain Monte Carlo (MCMC) for posterior sampling. While MCMC convergence theory exists, it does not guarantee finding the mode of the posterior efficiently. Coordinate-wise optimization methods (like those in Relevance Vector Machines) ensure local optimality but not global optimality.

2. Methodology

The author reframes the Bayesian variable selection problem as a non-convex optimization task, specifically exploiting a Difference-of-Convex (DC) structure.

A. Model Formulation

The paper utilizes a hierarchical normal model:
$$y | \theta, \sigma^2 \sim N(X\theta, \sigma^2 I_n)$$
$$\theta | \sigma^2, D \sim N(0, \sigma^2 D^{-1})$$
$$\sigma^2 \sim IG(a_0, b_0)$$
where $D = \text{diag}(d_1, \dots, d_p)$ represents precision parameters. Larger $d_j$ implies stronger shrinkage of $\theta_j$ toward zero.

B. The DC Decomposition

The negative log-marginal likelihood is transformed into a loss function $f(D) = g(D) - h(D)$, where both $g$ and $h$ are convex functions:

• $g(D)$: Contains the term $-\frac{1}{2}\log|D|$ (convex) and indicator functions for constraints.
• $h(D)$: Contains terms involving $\log|X^\top X + D|$ and the residual sum of squares. The convexity of $h$ is derived from the concavity of $-y^\top X(X^\top X + D)^{-1}X^\top y$ and the properties of the log function.

C. The DC Algorithm

The proposed algorithm iteratively minimizes the loss:

1. Linearization: At iteration $t$, replace the convex function $h(D)$ with its first-order affine approximation at $D^{(t)}$.
2. Subproblem: Solve the convex subproblem:
   $$D^{(t+1)} \leftarrow \arg\min_{D \in \mathcal{D}} \left\{ g(D) - \langle \nabla h(D^{(t)}), D - D^{(t)} \rangle \right\}$$
3. Update Rule: For box constraints $a \le d_j \le b$, the update has a closed-form solution involving the projection of the gradient of $h$.
4. Acceleration: To handle boundary convergence (where $d_j \to \infty$ implies $\theta_j = 0$), the algorithm includes a heuristic to explicitly test boundary values if the current iterate is close to them.

D. Handling Structured Sparsity

The framework extends to structured sparsity (groups, overlaps, monotonicity) by defining the constraint set $\mathcal{D}$ as a closed convex set. The projection step $P_\mathcal{D}$ is adapted using:

• Group constraints: Averaging within groups.
• Overlapping groups: Isotonic regression on decomposed regions.
• Radial monotonicity: Weighted Pool-Adjacent-Violators Algorithm (PAVA) for spatial decay constraints.

3. Key Contributions

Theoretical Insights

• Global Optimality Guarantee: Under mild compact set conditions, the DC algorithm is proven to converge to the global optimum at a linear rate.
• Benign Landscape: The paper demonstrates that despite the lack of log-concavity, the optimization landscape is "benign." Non-optimal stationary points (saddle points or local minima) are either ruled out by the DC update rule or correspond to boundary solutions that are easily identifiable.
• Extension to Posterior Modes: The theory extends from Type-II MLE to finding the mode of the marginal posterior, provided the prior can be expressed as a difference of convex functions.

Algorithmic Advantages

• Tuning-Free: The algorithm requires no hyperparameter tuning (unlike gradient descent step sizes).
• Efficiency: It is computationally efficient, with complexity $O(\min(n, p)^3)$ per iteration using the Woodbury matrix identity.
• Warm-Start Capability: It serves as an efficient alternative or a high-quality warm-start for MCMC methods, significantly reducing burn-in periods.

4. Experimental Results

A. Simulation Studies

• Setup: Compared the DC algorithm against Automatic Relevance Determination (ARD/EM) and Projected Gradient Descent (PGD) across dimensions $p \in \{500, \dots, 5000\}$.
• Performance:
  • Accuracy: All three methods achieved perfect variable selection (TPR=1, FPR=0) in simulations, validating the theoretical global optimality.
  • Speed: The DC algorithm was the fastest, converging in 79–292 iterations with negligible variance. It completed scenarios in under 2.3 seconds on a standard CPU.
  • Stability: Unlike ARD (which showed high variance in iterations due to coordinate-wise sensitivity) and PGD (which required more iterations due to non-adaptive step sizes), DC was robust and stable.

B. Real-World Applications

1. Dorothea Dataset (Drug Discovery):
   • Scale: $n=1,150$, $p=91,598$ features.
   • Result: The algorithm ran in ~9 minutes on a GPU. It drove 99.97% of features to the upper boundary (indicating zero coefficients), selecting only 30 features. This selection size was comparable to Lasso and Group Lasso but incorporated uncertainty in both coefficients and noise variance.
2. Ridgecrest Earthquake Aftershocks:
   • Task: Quantifying aftershock risk with radial monotonicity constraints (shrinkage increases with distance from the fault).
   • Result: The algorithm converged in 35 iterations (0.2 seconds). It successfully identified the fault rupture zone while smoothly shrinking activity to zero at the periphery, demonstrating the utility of structured sparsity constraints.

5. Significance and Conclusion

This paper fundamentally shifts the perspective on Bayesian variable selection from a purely sampling-based paradigm to an optimization-based paradigm.

• Resolution of the "Enigma": It proves that the non-convexity of the marginal likelihood does not hinder the ability to find the global mode reliably.
• Practical Impact: The proposed DC algorithm is a "plug-and-play" solution that is faster, more stable, and easier to implement than existing MCMC or coordinate-descent approaches.
• Future Directions: While currently focused on linear models, the author notes potential extensions to Generalized Linear Models (GLMs) via data augmentation, though this introduces theoretical challenges regarding tractability.

In summary, the work provides a rigorous theoretical foundation and a highly efficient algorithmic tool that makes Bayesian variable selection more accessible, reliable, and scalable for high-dimensional data.