hdfe-stream

Regression with several high-dimensional fixed effects, on data that does not fit in memory: linear (feols_stream), and Poisson, logit and probit (fepois_stream, feglm_stream), with the Kline–Saggio–Sølvsten bias correction for AKM variance decompositions (leave_out_kss).

from hdfe_stream import feols_stream

fit = feols_stream(
    "log_earn ~ age_squared + age_cubed | worker_id + firm_id + year",
    "data/*.parquet",              # a Parquet path or glob, never fully loaded
    workdir="scratch",             # where intermediates go: disk, not memory
    vcov={"CRV1": "worker_id"},
)
fit.summary()

Formulas follow pyfixest/fixest syntax and are parsed by pyfixest's own parser, so the syntax, the coefficient names and the results all match. The difference is where the work happens: rows are streamed off disk in batches and never held as a design matrix.

The data can be a Parquet path or glob, a Polars LazyFrame, or a Polars DataFrame. Using file paths or LazyFrames keeps memory usage to a minimum, but using DataFrames works well for smaller data and for trying things out.

Documentation: leetucker.github.io/hdfe-stream, with examples, the KSS bias correction, benchmarks and the function reference.

Features

• Any number of fixed effects, including one or none, and interactions (firm_id^year)
• pyfixest formula syntax: transformations, categoricals, interactions, event studies
• Several models in one call, sharing passes over the data
• Standard errors: iid, robust, multi-way clustered, and one-way CRV3
• Weighted least squares and 2SLS
• Poisson, logit and probit
• Varying slopes in fixed effects (worker_id[t]), on the streamed dimension only
• Wide designs, with hundreds of covariates
• Connected components, the largest connected set, and min_degree pruning
• AKM variance decompositions, in a streaming pass
• KSS leave-out variance components, with standard errors and weak-identification diagnostics
• The sample, row by row: which rows were dropped and why, and residuals joined back to the source
• Regression tables and plots through pyfixest
• Memory and disk controls, a choice of solvers, and logging
• Experimental GPU options: the solve on an NVIDIA GPU (device="cuda") and the Polars queries on cudf-polars (polars_engine="gpu")
• Simulated panels with known effects, for testing

When to use this, and when not to

When your data do not fit in memory, use this. It was built for that case: a matched employer–employee panel — tens of millions of rows, millions of worker effects, an AKM variance decomposition — on a machine or a secure enclave where the data are larger than the memory you are allowed. Every other library compared here needs the whole panel in memory before it starts; this one holds a bounded slice of it at a time, and uses scratch disk for the rest.

When they do fit, it is still a sensible default for medium and large regressions. From about 400,000 workers (3.4 million rows) up, it was the fastest of the libraries benchmarked: against pyfixest, xhdfe, R's fixest and Julia's FixedEffectModels.jl on the AKM regression and on wide designs, against pyfixest on Poisson and logit, and against xhdfe on the KSS decomposition, and at the largest sizes it used a fraction of their memory. At 5 million workers it fitted the AKM regression in 23 seconds on a 16-vCPU machine, about as fast as xhdfe on an NVIDIA H100 (24 seconds), without a GPU (see should I use a GPU?). The benchmarks are simulated employer–employee panels, a shape that suits its design (how it works); one whose fixed-effect cells are no smaller than the data may do less well.

It is also a reasonable general-purpose regression engine. OLS and weighted least squares, 2SLS, and binary logit and probit fit with no fixed effects or one, from a formula, with robust and clustered standard errors. On simulated data with 10 covariates, from a few million rows up it was within a small factor of the fastest libraries, ahead of some, and far lighter. At 20 million rows an OLS took 6.5 s and 2.5 GB, against 5.4 s and 7.5 GB for scikit-learn, 8.5 s and 5.7 GB for polars-ols, and 9.1 s and 15.7 GB for pyfixest; a logit took 9.4 s and 2.4 GB, against 12.8 s and 7.7 GB for statsmodels and 15.9 s and 20.6 GB for pyfixest. It is not the fastest: fixest was about 2.2 times faster at OLS and glum about 2.2 times faster at logit, scikit-learn, polars-ols and statsmodels were faster from 1 to 5 million rows, and a hand-written solve in NumPy is faster still. Those packages return coefficients or, at most, classical standard errors; here you also get clustering, weights, IV and formulas, and a memory use that stays near 2.5 GB where theirs reaches 5 to 20 GB. See the OLS, IV, logit and probit benchmarks.

On small data it gains you nothing. Below about 400,000 workers the in-memory libraries were faster — by a second or less on the AKM regression, and by six seconds on the smallest KSS decomposition — and lighter: an hdfe_stream process starts at about 0.8 GB. The same holds for plain regressions below about a million rows: a fit takes about a second, against tenths of a second for the fast libraries.

Use pyfixest when you need what a streaming design does not offer. It is more mature and more widely used, and has inference and clustering options this library lacks or limits: CRV3 with fixed effects that are not nested in the clusters, separation checks on the covariates and bias corrections for Poisson and logit, and more (see limitations). Use xhdfe when GPU speed is the priority: its CUDA backend made the most of the GPU in the GPU benchmark.

This library is at its strongest when:

• the data do not fit in memory, or fit so tightly that nothing else does;
• you have one very high-cardinality dimension (workers, patients, students) whose groups each touch only a few levels of the others;
• you want worker-specific slopes as well as intercepts (worker_id[t]), which multiplies the parameter count by the number of slopes;
• the design is wide: hundreds of covariates, such as a categorical expanded into fine indicators or many interactions. An in-memory library needs rows × covariates in memory before it starts; this needs one batch of rows × covariates at a time, with any number of fixed effects, including one or none (see many covariates);
• you need a hard, predictable memory ceiling — a shared cluster, or a job that must not be the one that gets killed.

You are trading memory for disk, and you need scratch space: about five times the size of your Parquet input while the fit runs, and only the result files once it finishes. See performance for what that buys. Point workdir= at it, or set the HDFE_STREAM_WORKDIR environment variable once (for example to a cluster's scratch filesystem) and leave workdir out; a workdir passed to a call takes precedence, and with neither, the system temporary directory is used.

How it works, in one paragraph

One fixed-effect dimension — normally the largest — is streamed: it is never represented as a vector in memory. Every other dimension is held as a vector sized by its number of levels, not rows. Rows are read in batches, hash-partitioned into buckets by the streamed dimension, sorted one bucket at a time, and reduced to a table of cells (one per combination of fixed-effect levels). The reduced normal equations for the non-streamed dimensions are then solved by conjugate gradient, and the streamed effects are recovered group by group in a final pass that writes residuals and fixed effects straight to Parquet. So memory scales with firms × years, not with workers × rows, and that asymmetry is the whole design. The covariates travel with the rows: what is held at once is one bucket or batch of rows × covariates, never all of them, which is why a wide design fits where an in-memory one does not. Full detail is in hdfe_stream/__init__.py.

Install

pip install hdfe-stream[formula]     # formula syntax needs pyfixest

polars, numpy, numba, pyarrow and scipy are required. The extras:

Without formula you still get StreamingHDFE, the lower-level interface that takes column names or Polars expressions instead of a formula string.

The experimental GPU options need CuPy and cudf-polars for your CUDA version (pip install cupy-cuda13x cudf-polars-cu13); there is no extra for them, since the packages differ by CUDA version. See CPU and GPU.

Performance

The standard three-way AKM regression, log_earn ~ age_squared + age_cubed | worker_id + firm_id + year, clustered by worker, on a simulated panel of 5 million workers and 42.5 million rows, against pyfixest, xhdfe, R's fixest and Julia's FixedEffectModels.jl. Every configuration agrees on the coefficients to within 2e-9.

On small data the in-memory libraries use less memory: an hdfe_stream process costs about 0.8 GB before it has done anything. The benchmarks cover every size from 212,000 rows up, designs with up to 503 covariates, Poisson and logit, and the KSS decomposition, with how each was measured. They also cover plain OLS, 2SLS, logit and probit without fixed effects or with one, against scikit-learn, polars-ols, statsmodels and others: at 20 million rows hdfe_stream is not the fastest, but is within a small factor of most and uses a fraction of their memory.

Should I use a GPU?

Probably not for this. hdfe_stream's GPU options are experimental, and on an NVIDIA H100 their gains were modest throughout and appeared only at large sizes: about 4% at 5 million workers with both options on, and no consistent gain below that. A GPU-oriented library such as xhdfe is the better choice if GPU performance is the priority: its CUDA backend took 24 s at 5 million workers, against 107 s on its CPU backend. But hdfe_stream on a 16-vCPU CPU machine took 23 s, and that machine costs about an eighth as much per hour as AWS's smallest H100 machine.

For plain OLS, 2SLS and logit with no fixed effects the answer is less clear-cut. At 20 million rows GPU libraries fitted them in 3 to 5 s (cuML, and a hand-written CuPy solve), 1.3 to 3.4 times faster than hdfe_stream on the CPU of the same machine, and not at all below a few million rows. But hdfe_stream's own GPU options did little for them: polars_engine="gpu" was slower at every size, and device="cuda", which runs a logit's IRLS steps on the GPU, saved 26% at 20 million rows, 10% at 5 million and nothing below that. And fixest, on the 16-vCPU CPU machine, fitted the same OLS in 2.9 s. See plain regressions on the GPU. CPU and GPU has the numbers.

Documentation

• Examples: one page per feature, on simulated data that ships with the package: formulas and standard errors, weights and IV, Poisson and logit, wide designs, varying slopes, connected sets, running out of core, and reporting.
• Leave-out variance components: the Kline–Saggio–Sølvsten bias correction for the AKM decomposition, its standard errors and weak-identification diagnostic, and every way it differs from the reference implementations.
• Resource limits and profiling: profile=True estimates a fit's memory and disk before you run it; the packages cap their threads and Polars' memory at a PBS/Slurm/cgroup allocation (the standalone stream-limits package), and their memory guards follow what is left of it.
• Benchmarks: time, memory and disk against pyfixest, xhdfe, fixest and FixedEffectModels.jl.
• Function reference: every option.
• Limitations: what it does not do.
• Changelog.

The scripts behind the examples are in examples/, and the benchmarks in benchmarks/.