sparsevb 0.1.2

Spike-and-Slab Variational Bayes for Linear and Logistic Regression

sparsevb

A Python package for Spike-and-Slab Variational Bayes for sparse variable selection in linear and logistic regression.

Description

This package implements the mean-field variational Bayes algorithm of Ray and Szabo (2022) for scalable Bayesian variable selection. Given a potentially high-dimensional design matrix X and response y, the algorithm approximates the posterior distribution under a spike-and-slab prior using coordinate-ascent variational inference (CAVI).

The spike-and-slab prior on each coefficient takes the form:

theta_i ~ gamma_i * Slab(lambda) + (1 - gamma_i) * delta_0

where gamma_i ~ Bernoulli(r) is the inclusion indicator and the slab is a Laplace (default) or Gaussian distribution. The variational approximation returns posterior inclusion probabilities gamma_i for each variable, along with posterior mean (mu) and standard deviation (sigma) parameters.

It is a direct port of the R sparsevb package, using the same underlying C++ implementation with a Pythonic API.

Installation

pip install sparsevb

System requirements: LAPACK, BLAS, and a C++14 compiler (usually available by default on Linux and macOS).

To install from source:

pip install .

Quick start

Linear regression

import numpy as np
from sparsevb import svb_fit_linear

# Synthetic data: 200 samples, 20 features, 3 active
np.random.seed(42)
n, p = 200, 20
X = np.random.randn(n, p)
beta_true = np.zeros(p)
beta_true[:3] = [1.5, -2.0, 3.0]
y = X @ beta_true + 0.5 * np.random.randn(n)

# Fit sparse model
res = svb_fit_linear(X, y)

# Point estimate of coefficients
print("Estimated coefficients:", res['mu'] * res['gamma'])

# Variable selection: gamma_i > 0.5 indicates inclusion
selected = res['gamma'] > 0.5
print("Selected variables:", np.where(selected)[0])

Logistic regression

from sparsevb import svb_fit_logit

# Binary response
logits = X @ beta_true
probs = 1 / (1 + np.exp(-logits))
y_binary = (np.random.rand(n) < probs).astype(float)

res = svb_fit_logit(X, y_binary)
print("Estimated coefficients:", res['mu'] * res['gamma'])

With intercept

y_with_intercept = X @ beta_true + 5.0 + 0.5 * np.random.randn(n)
res = svb_fit_linear(X, y_with_intercept, intercept=True)
print("Intercept:", res['intercept'])
print("Coefficients:", res['mu'] * res['gamma'])

API Reference

svb_fit_linear(X, y, **kwargs)

Fit a sparse linear regression model Y = X @ theta + noise.

Parameters:

Parameter	Type	Default	Description
X	array (n, p)	required	Design matrix
y	array (n,)	required	Response vector
max_iter	int	1000	Maximum CAVI iterations
tol	float	1e-5	Convergence tolerance on the max change in binary entropy of gamma
slab	str	"laplace"	Slab distribution: "laplace" or "gaussian". Laplace is recommended
intercept	bool	False	Include an intercept term
mu	array (p,)	None	Initial mean parameters. If None, initialized via Ridge regression
sigma	array (p,)	None	Initial std parameters. If None, set to ones
gamma	array (p,)	None	Initial inclusion probabilities. If None, initialized via Lasso
noise_sd	float	None	Known noise std. If None, estimated from data

Returns a dict with keys:

Key	Type	Description
mu	ndarray (p,)	Posterior mean of the slab component
sigma	ndarray (p,)	Posterior std of the slab component
gamma	ndarray (p,)	Posterior inclusion probabilities in [0, 1]
noise_sd	float	Noise std used (estimated or provided)
intercept	float or None	Estimated intercept value

svb_fit_logit(X, y, **kwargs)

Fit a sparse logistic regression model P(y=1 | X) = logistic(X @ theta).

Same parameters as svb_fit_linear except:

• y must be binary (0 or 1)
• noise_sd is not available (no noise parameter in logistic regression)

Returns the same dict, with noise_sd always set to None.

Interpreting the output

• Point estimate: the posterior mean of theta under the variational approximation is mu * gamma.
• Variable selection: a variable i is typically selected if gamma[i] > 0.5.
• Posterior uncertainty: sigma[i] gives the posterior standard deviation of the slab component for coefficient i. Note that variational Bayes tends to underestimate posterior variance.

Algorithm details

The algorithm uses coordinate-ascent variational inference (CAVI) to minimize the KL divergence between the mean-field variational family and the true posterior. Each iteration updates mu_i, sigma_i, and gamma_i for every coordinate (see Algorithm 1 in Ray & Szabo, 2022).

Key implementation choices:

• Laplace slab prior: unlike many VB approaches that use Gaussian slabs, this implementation defaults to Laplace slabs. Laplace priors avoid the excessive shrinkage of non-zero coefficients that Gaussian slabs can cause, leading to better signal recovery.
• Prioritized updating order: the CAVI algorithm is sensitive to the order in which coordinates are updated. This implementation sorts coordinates by decreasing |mu_init| (from a preliminary Ridge estimate), updating the most important coefficients first. This avoids poor local minima that can arise with lexicographic or random orderings.
• Noise estimation: for linear regression with unknown noise variance, the data is rescaled by an estimate of the noise standard deviation (empirical Bayes approach). When n > p, Ridge residuals are used; otherwise, the MAD of Lasso residuals is used.
• Convergence: the algorithm stops when the maximum change in binary entropy of the inclusion probabilities falls below tol.

Dependencies

• Python >= 3.8
• NumPy >= 1.20.0
• SciPy >= 1.6.0
• scikit-learn >= 1.0.0

Build dependencies: Cython >= 0.29.0, a C++14 compiler, LAPACK, BLAS.

License

GPL-3.0-or-later

References

• Ray, K. and Szabo, B. (2022). Variational Bayes for high-dimensional linear regression with sparse priors. Journal of the American Statistical Association, 117(539), 1270-1281. DOI:10.1080/01621459.2020.1847121