GEMSS: A Variational Bayesian Method for Discovering Multiple Sparse Solutions in Classification and Regression Problems

Selecting interpretable feature sets in underdetermined ($n \ll p$) and highly correlated regimes constitutes a fundamental challenge in data science, particularly when analyzing physical measurements. In such settings, multiple distinct sparse subsets may explain the response equally well. Identifying these alternatives is crucial for generating domain-specific insights into the underlying mechanisms, yet conventional methods typically isolate a single solution, obscuring the full spectrum of plausible explanations. We present GEMSS (Gaussian Ensemble for Multiple Sparse Solutions), a variational Bayesian framework specifically designed to simultaneously discover multiple, diverse sparse feature combinations. The method employs a structured spike-and-slab prior for sparsity, a mixture of Gaussians to approximate the intractable multimodal posterior, and a Jaccard-based penalty to further control solution diversity. Unlike sequential greedy approaches, GEMSS optimizes the entire ensemble of solutions within a single objective function via stochastic gradient descent. The method is validated on a comprehensive benchmark comprising 128 synthetic experiments across classification and regression tasks. Results demonstrate that GEMSS scales effectively to high-dimensional settings ($p=5000$) with sample size as small as $n = 50$, generalizes seamlessly to continuous targets, handles missing data natively, and exhibits remarkable robustness to class imbalance and Gaussian noise. GEMSS is available as a Python package ‘gemss’ at PyPI. The full GitHub repository at https://github.com/kat-er-ina/gemss/ also includes a free, easy-to-use application suitable for non-coders.

💡 Research Summary

The paper introduces GEMSS (Gaussian Ensemble for Multiple Sparse Solutions), a variational Bayesian framework designed to uncover several distinct sparse feature subsets simultaneously in underdetermined, highly correlated settings where the number of predictors far exceeds the number of observations (n ≪ p). Traditional sparse methods such as Lasso, standard spike‑and‑slab Bayesian inference, or stability selection typically return a single solution or rely on sequential masking strategies that either enforce disjointness or incur high computational cost. GEMSS overcomes these limitations by modeling the posterior distribution of the regression coefficients as a mixture of m diagonal Gaussian components, each intended to capture a different mode corresponding to a plausible sparse solution.

Key methodological components include: (1) a structured spike‑and‑slab (sss) prior that explicitly controls the sparsity level D by defining a support set of exactly D non‑zero entries; for small p all supports are enumerated, while for large p they are sampled. (2) Variational approximation of the multimodal posterior via a mixture of Gaussians, with mixture weights α_k learned jointly. (3) Maximization of the Evidence Lower Bound (ELBO) using stochastic gradient descent (Adam) and the implicit reparameterization trick to obtain unbiased gradients for the mixture parameters (means μ_k, variances σ_k, weights α_k). (4) A Jaccard‑based diversity regularizer, λ_J · J_avg(q), that penalizes overlap between the supports of different mixture components, thereby encouraging genuinely distinct solutions. (5) Native handling of missing data through masking and variational expectation, avoiding separate imputation steps.

The authors evaluate GEMSS on a comprehensive synthetic benchmark comprising 128 experiments organized into seven tiers: baseline (n < p), high‑dimensional (p up to 5 000), sample‑rich (n ≥ p), noisy observations, missing values, class imbalance, and regression tasks. Performance metrics include support recovery rate (fraction of true variables recovered), Jaccard diversity between recovered supports, and predictive accuracy (AUC for classification, RMSE for regression). Across all tiers, GEMSS consistently recovers a larger number of true sparse supports than competing methods such as Sequential Lasso, ALFESE, SES, and evolutionary niching algorithms. In the most challenging high‑dimensional, low‑sample regime (p = 5 000, n = 50), GEMSS identifies 5–10 distinct solutions with an average support recovery of >85 % while maintaining predictive performance comparable to or better than single‑solution baselines. The diversity regularizer proves essential for preserving distinctness; without it, mixture components tend to collapse onto the same mode.

A notable practical contribution is the release of a well‑documented Python package (gemss) on PyPI and a GUI‑based “GEMSS Explorer” application, allowing non‑programmers to load data, set hyperparameters (D, λ_J, m), run the algorithm, and visualize the resulting feature sets. The paper also provides guidance on selecting hyperparameters, discusses computational complexity (linear in p for each gradient step, with additional cost for support sampling), and acknowledges limitations: the need to pre‑specify the number of solutions m and sparsity level D, potential scaling issues when p is extremely large, and the current focus on linear models.

Future work outlined includes automatic hyperparameter tuning, extensions to non‑linear models (e.g., Bayesian neural networks), and distributed implementations to handle massive datasets. Overall, GEMSS offers a principled, scalable, and user‑friendly approach for generating multiple interpretable sparse models, addressing the “Rashomon effect” in high‑dimensional scientific data analysis and providing domain experts with a richer set of hypotheses for downstream validation.