pybiglasso

Lasso and elastic-net regression for data larger than memory.

If your penalized regression won't fit in memory, pybiglasso fits it where the data already is: on disk. A 12 GB design matrix fits under a 3 GB memory limit, on an ordinary desktop, in under two minutes (see Performance).

Save the matrix as a numpy memmap and pybiglasso fits the whole regularization path, with cross-validation to choose the penalty, without loading it into memory.

To try it without installing anything, open the playground. It runs the same Rust solver, compiled to WebAssembly, in your browser, on example datasets or your own CSV or Parquet file, which never leaves your machine. It fits lasso and elastic-net paths for Gaussian or binomial responses, and can cross-validate them.

It's also fast when the data does fit in memory. The timings quoted here use binomial (logistic) fits. On a full 100-step lasso path over 5,000 samples by 50,000 features, it is 2.4x faster than R's glmnet, 2.5x faster than R biglasso, and 5.2x faster than skglm, the fastest of the other Python solvers we tested. Across the four benchmark problems, with 1,000 to 50,000 features, it runs 2.1–4.7x faster than glmnet, 2.5–9x faster than R biglasso, and 5–9x faster than skglm.

It is an unofficial Python port of the R package biglasso, and reproduces its results to floating-point precision (see Validation).

Install

pip install pybiglasso

The core depends only on numpy and scipy; the solver is compiled Rust, in the wheel. Extras add optional integrations:

Quick start

import numpy as np
import pybiglasso as bl

# open the matrix on disk without loading it (a read-only memmap); store it
# column-major (Fortran order)
X = np.load("X.npy", mmap_mode="r")
y = np.load("y.npy")

path = bl.biglasso(X, y, alpha=1.0, nlambda=100)   # alpha=1: lasso
path
# BiglassoPath(500000 features, 100 lambdas from 0.112 to 0.000112, alpha=1.0)

path.coef                  # scipy.sparse (features x lambdas)
path.n_nonzero             # model size along the path
path.predict(X_new)        # (rows x lambdas); reads only the columns in use
path.score(X_val, y_val, metric="mse")

path.to_polars()           # lambda_, feature, coefficient (nonzero only)

To write such a file, open it as a column-major memmap and fill it a block of columns at a time, so the whole matrix never has to be in memory:

X = np.lib.format.open_memmap(
    "X.npy", mode="w+", dtype=np.float32,
    shape=(n_samples, n_features), fortran_order=True,
)
for start, block in read_blocks():
    X[:, start : start + block.shape[1]] = block
X.flush()

X can be float64, float32, any integer type, or bool. Use the narrowest type that holds your values exactly: the fit is mostly reading X, so a narrower type is faster as well as smaller (see storage).

Unpenalized covariates (say age, sex, and principal components) take a penalty factor of 0:

pf = np.ones(X.shape[1])
pf[:12] = 0
path = bl.biglasso(X, y, penalty_factor=pf)

Gaussian regression

Gaussian is the default family in 0.3.0, matching R biglasso. This uses the shared Rust solver in Python and WASM; the deprecated numba backend supports binomial fits only. Gaussian paths predict continuous responses and report residual sum of squares as deviance; scoring also accepts "mse" and "r2". BiglassoRegressor and BiglassoRegressorCV provide sklearn's regression API. For logistic regression, explicitly pass family="binomial":

classification = bl.biglasso(X, y_binary, family="binomial")
classification.predict_proba(X_new)

This is a breaking default change from 0.2.x: existing logistic callers must select binomial explicitly. Classifier wrappers already do so. Numba remains available for explicit binomial fits through 0.3.x and is removed in 0.4.0.

scikit-learn

pybiglasso also offers a scikit-learn interface: BiglassoClassifier fits at one penalty strength, and BiglassoClassifierCV chooses it by cross-validation. BiglassoRegressor and BiglassoRegressorCV provide the corresponding Gaussian estimators. All four pass scikit-learn's estimator checks and work in pipelines and grid searches (install pybiglasso[sklearn]).

from pybiglasso import BiglassoClassifierCV

clf = BiglassoClassifierCV(l1_ratio=1.0, cv=5, rule="1se").fit(X, y_binary)
clf.alpha_, clf.coef_

The estimators use scikit-learn's parameter names, which differ from R's:

Why this exists

At 23andMe, I used R's biglasso package to fit penalized models on large genotype/phenotype datasets without loading them into memory. Over a holiday break, I ported it to Python, partly to gain a better understanding of how it worked and partly because, assuming I could get it working, it would be convenient to stay within the Python ecosystem. I validated it against R biglasso, matching its coefficients, lambda grid, and iteration counts to floating-point precision, and benchmarked it on data both smaller and larger than memory.

In 2026 I revived the port for release, working with coding agents (Claude Code and Codex). That work added the scikit-learn estimators and cross-validation, extended the validation (against R across many more configurations, and with checks independent of R, which turned up the differences listed below), and carried the performance work and benchmarks further.

This is an independent port, not affiliated with or endorsed by the biglasso authors: Yaohui Zeng, Chuyi Wang, Tabitha Peter, and Patrick Breheny. If you use it, please cite their paper (below).

Validation

pybiglasso is checked two ways: against R biglasso, which it ports, and against checks that don't depend on R, which also cover the few places where the two differ (see Differences).

Against R. Every push runs pybiglasso against stored R biglasso 1.7.2 output for 192 configurations (152 binomial and 40 Gaussian), comparing coefficients, deviance, the lambda grid, centering and scaling, strong-rule rejections, and the iteration count at every lambda, with rtol=1e-9. The grid covers wide and tall problems, genotype-like dosage data (including int8 and uint8 storage), memmaps, float32, row subsets, penalty factors (including zeros), user-supplied and linear lambda grids, dfmax truncation, and saturation. The comparison runs pybiglasso with r_compatible=True; for cases with penalty factors it passes R pybiglasso's lambda grid, which sets aside the zero-penalty-factor case under Differences. A weekly job reruns the grid live against the current CRAN release, to catch upstream changes.

Independent of R.

• Optimality. Converged test paths satisfy the elastic-net optimality (KKT) conditions to 1e-6, on wide, correlated, partly unpenalized, and genotype-like designs.
• Other solvers. R's glmnet, which minimizes the same objective, and scikit-learn's saga solver, both tightly converged, find the same binomial coefficients: each pair passes the optimality conditions, and then their distance is certified, to about 1e-7 on tall data and, through the duality gap, to about 1e-5 on wide data (all but the wide lasso, whose optimum need not be unique). glmnet's path is stored from a pinned release and checked at every lambda where no fitted log-odds exceeds 10, beyond which pybiglasso, like R biglasso, clips them. In the benchmarks, glum reaches the same elastic-net objective as well.
• Gaussian checks. Six pinned glmnet paths use its generic gaussian() family, avoiding the response scaling of its fast Gaussian path. Independent checks cover KKT conditions, curvature, exact-arithmetic duality gaps, analytic solutions, and invariances.
• Where R differs. R's elastic-net update (reproduced by r_compatible=True) doesn't satisfy the optimality conditions, and its coefficients differ from saga's by 2–50% of their size on the test problem. This is the basis for the elastic-net difference under Differences.
• Shifted data. Adding 1e8 to every column of X leaves the solution unchanged and still optimal; this is the case described under Differences as data far from zero.

uv run pytest               # against stored R output; no R needed
make test-r                 # live, in a container (CONTAINER=podman works too)
make test-r R_VERSION=latest
make soak                   # larger memory-mapped problems; for a big machine
make soak-hypothesis        # property tests, wide random search (part of soak)

Performance

The binomial benchmarks below measure speed only. Correctness is established separately, against R biglasso and by checks independent of it (see Validation). Each solver fits the same problem along the same 100-value lambda path, at its own default tolerance, and its accuracy is reported alongside its time. Machine: AMD Ryzen 5 5600X (6 cores, 12 threads) and 64 GB of RAM, with Python 3.14 and R 4.6.1. Each time is the median of 3 runs, or a single run for a solver slow enough to take over a third of its time limit. Every multithreaded solver is shown at its fastest thread count: 1, 6, or 12. The scripts and raw results are in benchmarks/.

Which solvers. R biglasso, the original; R's glmnet, whose coordinate descent biglasso builds on; and every maintained Python package, installable with pip, that fits L1 or elastic-net logistic regression on a CPU: scikit-learn (its liblinear and saga solvers), glum, skglm, and statsmodels.

Other Python solvers. Each competitor was stopped once it took 10x as long as pybiglasso, which would otherwise have stretched the run to many hours. (The competitors ran against 0.1.x, whose times set the cap; 0.2.0 is faster, so against it the caps fall at 11–25x, and glum's one finish, the 50,000 x 2,000 elastic net on 12 threads, is at 18x.) Only skglm finishes within that on every dataset, at 5–7x for the lasso and 6–9x for the elastic net. statsmodels (lasso only) ran only on the tall datasets: its L1 fit uses a general-purpose optimizer with dense matrices over twice the columns, which would need about 12 GB and 76 GB on the wide ones. To show how large the gaps get, the smallest problem (5,000 x 1,000) also ran with no cap:

R biglasso 1.7.2. Times are in seconds; "best" is each package at its fastest thread count. For the lasso, both packages fit identical paths. For the elastic net, each ran its own update (see Differences); pybiglasso with r_compatible=True, which matches R, takes about as long as its default.

R's adaptive screening rule is no faster than SSR on these problems.

glmnet 5.1. glmnet standardized X itself (standardize = TRUE, its default), as pybiglasso does, so both times include standardizing; the Python solvers, which don't standardize, were given X already standardized. glmnet fit with an intercept, along the same lambda grid, at its default tolerance; it runs on one thread. Its objective is the benchmark's. It fit every path in full, at 2.1–4.7x pybiglasso's best time (1.8–3.1x pybiglasso's on one thread), and about as accurately: a worst relative suboptimality of 2e-8 to 3e-6, against pybiglasso's 6e-8 to 3e-6.

Out of core. This test fits a 12 GB float64 matrix (20,000 x 75,000) with memory capped at 3 GB. pybiglasso read a numpy memmap under a cgroup memory limit; R read a file-backed big.matrix in a container with the same limit, and both were measured from their cgroup's accounting. Before each run, the matrix was evicted from the page cache. Both fit the same 20-lambda lasso path and agree to 1e-9, with identical iteration counts. pybiglasso ran twice, and both runs are shown; R ran once.

Here the disk sets the pace, and times vary by a third from run to run. Extra threads don't reliably help: they evict pages each other still need.

Not compared, by the rule above:

• Other languages: Julia's GLMNet.jl, a wrapper of glmnet's earlier Fortran solver, and Lasso.jl, written in Julia.
• Other hardware: RAPIDS cuML fits on a GPU, and only data that fits in GPU memory.
• Unmaintained: python-glmnet, archived, with its last release in 2020.
• A different problem: celer's logistic regression fits the lasso only, and no intercept. skglm, from the same authors, fits an intercept and the elastic net, and is compared above.
• Not a solver: benchopt is a benchmarking framework; its logistic-regression benchmark fits without an intercept or standardization, which pybiglasso doesn't support yet.

Gaussian path benchmarks

Gaussian timings are separate from the binomial results above. The runner benchmarks/run-gaussian.sh compares Rust, R biglasso's SSR solver and glmnet's native family="gaussian" path on the same 100 supplied lambdas. The grid comes from the standardized score formula, not a pybiglasso fit. All receive raw X and standardize it during the fit; data generation, imports, container startup and R big.matrix construction are outside the timers. Gaussian Python times biglasso() itself; the older binomial Python table also includes coefficient export. Each repetition uses a fresh process. The raw results record every time, input fingerprints, versions, warnings, fitted path length, first-principles KKT residuals and objective differences versus glmnet. R's BLAS is limited to one thread inside the container; biglasso's explicit thread count controls its own parallel work.

Glmnet internally rescales the response by its population standard deviation sᵧ and returns coefficients in the original response units. To match the same standardized elastic-net objective, the native path receives α′ = α / (α + (1−α)sᵧ) and λ′ = λ(α + (1−α)sᵧ). For the lasso this leaves α and λ unchanged. No further coefficient rescaling is needed. The generic gaussian() path remains a correctness reference in tests, not a timed competitor. See the glmnet response-scaling documentation. The solvers retain their default tolerances, so accuracy accompanies timing; the KKT residual is measured, not an eps-to-accuracy promise.

The table uses the same AMD Ryzen 5 5600X (6 cores, 12 threads) as the binomial benchmarks, with pybiglasso 0.3.0, Python 3.14.0, NumPy 2.5.4, R 4.6.1, biglasso 1.7.2 and glmnet 5.1. Times are seconds, the median of three runs; both thread counts are shown rather than choosing the faster one. The CPU was 96–97% idle across twelve ten-second samples before the run; load was sampled throughout and returned to 95–97% idle afterward.

Shapes are requested sizes; dropping generated constant columns leaves 973 and 19,199 features respectively. Accuracy and warnings describe each row's final repetition. All 20 paths contain 100 lambdas with no warnings. The KKT diagnostic is an absolute subgradient residual in response units (the columns have unit variance). The largest residuals are about 1.0e-4 for Rust/R biglasso and 1.1e-4 for native glmnet. On tall data the smallest L1 penalty λ_min α is 1.50e-3: the largest path-wide residual is therefore about 7% of that threshold at default tolerances. On wide data λ_min α is 0.0762, and the corresponding ratio is about 0.01%. Raw records report this conservative path-wide ratio; it is not a relative coefficient error. The largest relative objective difference versus native glmnet is 1.9e-6; this pairwise difference is not a certified suboptimality bound. Regenerate the table with:

python benchmarks/gaussian-table.py benchmarks/results/path-gaussian.jsonl
make r-image R_VERSION=4.6.1  # once, before the timing run
benchmarks/run-gaussian.sh
# Quick adapter/timing smoke (not a reportable benchmark):
DATASETS=smoke-500-100 THREADS=1 REPEATS=1 benchmarks/run-gaussian.sh /tmp/gaussian-smoke.jsonl

Differences from R biglasso

pybiglasso follows R biglasso 1.7.2 closely, matching it coefficient for coefficient and down to its iteration and screening counts. It departs in a few places, mostly where we think R's behavior may not be what was intended; if we've erred we'd be glad to be corrected. Each departure is marked in the code with a comment beginning "Differs from R biglasso".

• Elastic net. For alpha < 1, R's logistic coordinate update divides by v * (1 + l2), where the elastic-net objective in its documentation calls for v + l2, so the two solve somewhat different problems. pybiglasso uses v + l2, as glmnet and scikit-learn do; r_compatible=True reproduces R exactly. The lasso is the same either way.
• Zero penalty factors. With a penalty.factor of 0 and the default lambda grid, R's grid comes out as lambda = 0 at every step. pybiglasso instead fits the unpenalized columns first, and starts the path at the smallest lambda that keeps every penalized column at zero.
• Default lambda_min. R chooses 0.001 or 0.05 by comparing nrow(X) with ncol(X); pybiglasso compares the number of rows actually fit, which differs when row_indices is given.
• Data far from zero. For columns whose mean is large next to their spread (say, values near 1e8 that vary by 1), standardizing from raw sums loses precision in floating point. R's scales, and so its fit, can be affected without a warning; in our tests some scales came out as NaN. pybiglasso computes these columns' moments from centered values, so a constant shift of X leaves the solution unchanged.
• Invalid input. NaN or infinite values in the rows being fit, and non-finite penalty factors or lambdas, raise ValueError; R instead drops the affected column or returns the null model.
• Binomial deviance. deviance is the deviance, -2 log-likelihood, of the returned model. R reports loss, half the deviance, computed before the final coordinate sweep and with probabilities clipped at 0 and 1, so it can occasionally be infinite. r_compatible=True keeps R's version.
• Binomial saturation. stop_on_saturation=False, an option R doesn't have, fits every lambda, even past R's saturation cutoff.
• Gaussian constant responses. R's strict convergence comparison against a zero null loss uses the full iteration budget for a constant response. pybiglasso stops once the coefficients no longer change; r_compatible=True retains R's comparison.
• Gaussian zero maximum penalty. R's logarithmic grid calculation can produce non-finite values when the maximum penalty is zero. pybiglasso returns a finite all-zero grid, whose fits solve the unpenalized problem.
• Unrepresentable grids. If an automatic grid overflows (for example, with extremely small positive alpha), R returns non-finite penalties, and pybiglasso raises ValueError, for both families and in both modes. Finite supplied grids remain supported, including in compatibility mode.
• Gaussian loss after an iteration limit. R can leave loss at zero when the budget ends with outstanding KKT violations. pybiglasso reports the returned model's residual sum of squares; r_compatible=True retains R's bookkeeping.

Threads, cross-validation, and storage

pybiglasso reads the matrix one column at a time, through the memmap, and never loads or copies it; strong-rule screening skips most columns at most steps. Cross-validation folds are row indices into the same matrix on disk.

• A single fit is effectively single-threaded, as is R biglasso's. Coordinate descent updates one coefficient at a time, and that sequential sweep is most of the work: profiling an elastic-net path on 5,000 x 1,000 puts over 80% of the time in it, against 1-2% in the strong-rule and KKT scans that run in parallel (n_threads=). make profile reproduces this with Python 3.15's sampling profiler. Out of core, threads can help when memory to spare holds the page cache, by keeping more disk reads in flight. Under a tight memory limit they hurt (see Performance), so use n_threads=1 there. Results are bit-identical at any thread count.
• Cross-validation scales with cores. BiglassoClassifierCV(n_jobs=...) fits folds in parallel threads that share one X (and, for a memmap, one page cache); results don't depend on n_jobs. With the deprecated numba backend this needs numba's tbb or omp threading layer; where numba falls back to workqueue, as on macOS by default, folds run one at a time.
• Cross-validation can stop each fold early. With BiglassoClassifierCV(early_stopping=10), a fold's path stops once its held-out score has gone 10 alphas without improving, as bigstatsr's big_spLogReg does, and the fit on all the rows stops at the chosen alpha. On wide data, where the end of the path costs most, this saves much of the work; the cost is that a better score past the stopping point goes unseen. It is off by default.
• Store X compactly. Since the work is mostly moving bytes, a narrower dtype is faster: small integers (counts, categories, 0/1 indicators) as int8 or uint8 take 8x less memory and disk than float64, and fit 1.3–1.6x faster in memory; float32 halves the memory and is about as fast. Each value is read as a float64, so values that a narrower dtype holds exactly give bit-identical results.
• Bounds checking is on by default. Every index into an array is checked, which guards against a bug in pybiglasso reading past the end of X (inputs, row_indices among them, are validated before a fit either way). It costs little where X's rows are read in place: a whole Fortran-ordered X, or a subset of at least half its rows without repeats, as in cross-validation, which is read through a mask. Where every read goes through an index or a stride (C-ordered X, or a smaller or repeating row_indices), it costs up to about 15%. PYBIGLASSO_BOUNDSCHECK=0 turns it off; the path is the same (lambdas, iterations, screening), with coefficients identical in the Rust backend and equal to about 1e-11 in numba's, whose vectorized sums the check reorders. numba reads the setting at import, and caches compiled kernels only with it on.

How the code is organized

R biglasso implements the logistic solver in two long C++ functions: cdfit_binomial_ssr standardizes, sets up the lambda grid, and runs the lambda loop with its strong rule and KKT checks, and cdfit_binomial_lambda_core runs the coordinate descent at each lambda. That keeps the algorithm in a few places, close to how it's written in the papers.

pybiglasso splits the same algorithm into one small function per step, each documenting the state it changes, plus a short driver, fit_lambda, that reads like the algorithm's outline. To my eye that's easier to read, test, and modify. It's a matter of taste, though; some excellent programmers prefer the long inlined function you can read straight through (John Carmack has a good essay on why). The solver takes exactly the same iterations as R's at every lambda; the R comparison checks this in every case.

The solver is written twice, in the same steps. The one that runs is Rust (rust/src/solver.rs): one call from Python fits the whole path. The readable one, tests/support/reference.py, is plain numpy, step by step, and imports nothing from pybiglasso; the tests hold the Rust solver to it, and both to R.

The Gaussian driver and coordinate sweeps live in rust/src/gaussian.rs, porting R's Gaussian SSR solver. They share column access, standardization, and fixed-order reductions with the binomial solver, but use residual updates with constant curvature, without IRLS. The independent NumPy reference also covers Gaussian paths.

Every sum over rows runs in an order fixed by the data's shape alone. Threads split the work by column, and a tall column's rows into 2,048-row chunks whose sums are added in order. Within a chunk, eight running sums are combined in a fixed order, so vector instructions of any width (AVX2, AVX-512, ARM's NEON) give the same result, and no multiply and add are fused into one instruction, which would round differently. Exponentials and logarithms come from a portable implementation rather than the operating system's math library, whose last bits differ between Linux, macOS, and Windows. So results are bit-identical at any thread count, on any CPU, and on every operating system with a wheel.

The deprecated numba backend (backend="numba") is in src/pybiglasso/binomial.py, one numba kernel per step; it is slated for removal in 0.4.0.

Scope

• Binomial (logistic) and Gaussian (least-squares) families; Cox is not ported.
• SSR (sequential strong rule) screening only. R uses its newer Hybrid and Adaptive rules only for the plain lasso without penalty factors, and they were no faster than SSR on our benchmarks.
• With the deprecated numba backend, the first call in a session pays a few seconds of compilation; after that, compiled kernels are cached on disk.

Citation

Zeng Y, Breheny P (2021). The biglasso Package: A Memory- and Computation-Efficient Solver for Lasso Model Fitting with Big Data in R. The R Journal, 12(2), 6–19. doi:10.32614/RJ-2021-001

License

GPL-3.0-or-later, as R biglasso is GPL-3.