Benchmark Methodology

Transparent and reproducible performance testing

Overview

This page documents our benchmark methodology to ensure transparency and reproducibility. All benchmarks can be reproduced using the provided scripts.

YOCO Compression + Sparse Matrices

leanfe uses the YOCO (You Only Compress Once) strategy from Wong et al. (2021) combined with sparse matrices for all standard error types (IID, HC1, and clustered):

For IID/HC1 Standard Errors

  1. Compress: GROUP BY (regressors + fixed effects) → compute n, sum(y), sum(y²) per group
  2. Build sparse design matrix: FE dummies as scipy.sparse (Python) or Matrix (R)
  3. Solve sparse WLS: Weighted least squares with sparse operations
  4. Standard errors: Compute from grouped residual sum of squares

For Clustered Standard Errors (Section 5.3.1)

Clustered SEs require additional steps to aggregate scores within clusters:

  1. Within-cluster compression: GROUP BY (regressors + fixed effects + cluster_id) — each compressed record belongs to exactly one cluster
  2. Build sparse design matrix: Same as above
  3. Solve sparse WLS: Same as above
  4. Compute ẽ⁰: Sum of residuals per group = sum_y - n × ŷ
  5. **Build sparse cluster matrix W̃_C**: Indicator matrix (G × C) mapping groups to clusters
  6. Aggregate scores: cluster_scores = W̃_Cᵀ @ (X × ẽ⁰) — vectorized, no loops!
  7. Meat matrix: meat = cluster_scoresᵀ @ cluster_scores
  8. Sandwich: V(β̂) = (X'X)⁻¹ × meat × (X'X)⁻¹ × adjustment

The sparse cluster matrix W̃_C enables vectorized aggregation across all clusters simultaneously, avoiding the traditional loop over C clusters.

This is lossless - coefficients and standard errors are identical to the full computation, but up to 23x faster for discrete regressors at scale.

Why Sparse Matrices? (leanfe’s Contribution)

The original YOCO paper by Wong et al. (2021) describes the compression strategy and the mathematical framework for computing coefficients and standard errors from sufficient statistics. However, the paper focuses on the statistical methodology and does not address the practical implementation challenge of how to efficiently represent the design matrix when fixed effects are included.

Specifically, the paper:

  • ✅ Describes compression via GROUP BY on features (M̃, ỹ⁰, ỹ⁰⁰, ñ)
  • ✅ Shows how to compute β̂ and V(β̂) from compressed data
  • ✅ Covers IID, HC1, and clustered standard errors
  • ❌ Does not discuss how to handle FE dummy variables efficiently
  • ❌ Does not mention sparse matrices

leanfe extends YOCO with sparse matrices for the FE dummy variables. This is a natural and mathematically equivalent extension that provides significant performance benefits.

The Problem: Dense FE Dummies Are Wasteful

After YOCO compression, we have G groups (compressed rows) and need to estimate:

\[\bar{y}_g = X_g \beta + D_g \gamma + \epsilon_g\]

where \(D_g\) is a matrix of FE dummy variables. With F total FE levels across all FEs, \(D\) has dimensions G × (F-1). For example:

  • G = 50,000 compressed groups
  • F = 600 FE levels (500 + 100)
  • Dense \(D\) matrix: 50,000 × 599 = 30 million entries

But each row of \(D\) has only one 1 per FE (and zeros elsewhere). With 2 FEs, each row has exactly 2 non-zero entries out of 599 columns — that’s 99.7% zeros.

The Solution: Sparse Matrix Representation

Sparse matrices store only non-zero entries as (row, column, value) triplets:

Representation Storage for G=50K, F=600
Dense matrix 30 million floats (~240 MB)
Sparse matrix 100K triplets (~2.4 MB)
Savings ~100x less memory

Why This Works Mathematically

The key insight is that sparse matrix operations are mathematically identical to dense operations — they just skip the zero multiplications:

# Dense: X'X computes all G×F×F multiplications
XtX_dense = X.T @ X  # O(G × F²)

# Sparse: X'X only computes non-zero contributions  
XtX_sparse = X_sparse.T @ X_sparse  # O(G × nnz_per_row × F)

Since nnz_per_row (non-zeros per row) equals the number of FEs (typically 2-3), sparse operations are O(F) faster than dense operations.

Correctness Guarantee

The sparse implementation produces bit-identical results to a dense implementation because:

  1. Same linear algebra: \(X'X\), \(X'y\), and \((X'X)^{-1}\) are computed identically
  2. Same sufficient statistics: YOCO’s n, sum(y), sum(y²) are unchanged
  3. Same standard errors: The sandwich formula uses the same matrices

We validate this in our test suite by comparing sparse vs dense results.

Performance Impact

Operation Dense Sparse Speedup
Build design matrix O(G × F) O(G × n_fe) ~100x
Matrix multiply X’X O(G × F²) O(G × n_fe × F) ~100x
Memory allocation O(G × F) O(G × n_fe) ~100x

Where n_fe is the number of FE columns (typically 2-3) and F is total FE levels (can be 1000s).

Implementation Details

Python uses scipy.sparse.csr_matrix:

from scipy import sparse

# Build in COO format (fast construction)
X_fe = sparse.coo_matrix((data, (rows, cols)), shape=(G, F-1))

# Convert to CSR for fast arithmetic
X_fe = X_fe.tocsr()

# Combine with regressors
X = sparse.hstack([X_reg_sparse, X_fe], format='csr')

R uses Matrix::sparseMatrix:

library(Matrix)

X_fe <- sparseMatrix(i = rows, j = cols, x = data, dims = c(G, F-1))
X <- cbind(Matrix(X_reg, sparse = TRUE), X_fe)

When Sparse Matrices Help Most

Sparse matrices provide the biggest benefit when:

  • Many FE levels: F > 100 (more columns to skip)
  • Good compression: G << N (fewer rows to process)
  • Few FEs: n_fe = 2-3 (very sparse rows)

This is exactly the typical panel data scenario (firm + year FEs with discrete treatments).

For clustered SEs, the cluster matrix W̃_C is also sparse (one 1 per row), enabling efficient score aggregation even with thousands of clusters.

Data Generation

Synthetic Panel Data

All benchmarks use the same standardized data generation for consistency:

import numpy as np
import polars as pl

def generate_benchmark_data(n_obs, n_fe1=500, n_fe2=100, n_clusters=None, seed=42):
    """
    Generate synthetic benchmark data.
    
    Parameters
    ----------
    n_obs : int
        Number of observations
    n_fe1 : int
        Number of FE1 levels (default 500, e.g., firms)
    n_fe2 : int
        Number of FE2 levels (default 100, e.g., products)
    n_clusters : int, optional
        Number of clusters for clustered SE benchmarks
    seed : int
        Random seed for reproducibility
    """
    np.random.seed(seed)
    
    # Fixed effect IDs
    fe1 = np.random.randint(1, n_fe1 + 1, n_obs)
    fe2 = np.random.randint(1, n_fe2 + 1, n_obs)
    
    # Fixed effect values
    fe1_effects = np.random.normal(0, 1, n_fe1 + 1)[fe1]
    fe2_effects = np.random.normal(0, 0.5, n_fe2 + 1)[fe2]
    
    # Regressors - discrete for YOCO compression benefits
    treatment = np.random.binomial(1, 0.3, n_obs).astype(float)
    x1 = np.random.choice([0.0, 1.0, 2.0], n_obs)  # Discrete
    
    # Outcome: y = 2.5*treatment + 1.5*x1 + fe1 + fe2 + noise
    y = (2.5 * treatment + 1.5 * x1 + 
         fe1_effects + fe2_effects + 
         np.random.normal(0, 1, n_obs))
    
    data = {
        "y": y,
        "treatment": treatment,
        "x1": x1,
        "fe1": fe1,
        "fe2": fe2,
    }
    
    # Add cluster column if needed
    if n_clusters is not None:
        cluster_id = np.random.randint(0, n_clusters, n_obs)
        cluster_effects = np.random.normal(0, 0.5, n_clusters)[cluster_id]
        data["cluster_id"] = cluster_id
        data["y"] = y + cluster_effects  # Add cluster-level variation
    
    return pl.DataFrame(data)

Data Characteristics

Parameter Value Rationale
True treatment effect 2.5 Known ground truth for validation
True x1 effect 1.5 Known ground truth for validation
FE1 cardinality 500 High-dimensional (e.g., firms)
FE2 cardinality 100 Medium-dimensional (e.g., products)
Treatment Binary (0/1) 30% probability, discrete for compression
x1 Discrete (0, 1, 2) Enables YOCO compression
Error variance 1.0 Standard normal errors
Seed 42 Reproducibility

Dataset Sizes

All benchmarks use consistent sizes for comparability:

Size Observations Use Case
Small 100,000 Quick tests
Medium 1,000,000 Standard benchmarks
Large 10,000,000 Stress tests
XL 100,000,000 Memory efficiency / scale tests

Hardware and Software

Test Environment

Benchmarks were run on:

  • Machine: Apple MacBook Pro
  • Processor: Apple M-series
  • RAM: 16 GB
  • OS: macOS

Software Versions

import sys
import polars
import duckdb
import numpy

print("Software Versions:")
print(f"  Python: {sys.version.split()[0]}")
print(f"  Polars: {polars.__version__}")
print(f"  DuckDB: {duckdb.__version__}")
print(f"  NumPy: {numpy.__version__}")
Software Versions:
  Python: 3.11.14
  Polars: 1.36.1
  DuckDB: 1.4.3
  NumPy: 2.3.5

Measurement Methodology

Time Measurement

import time

start = time.time()
result = leanfe(...)
elapsed = time.time() - start
  • Metric: Wall-clock time in seconds
  • Iterations: 3 runs per configuration
  • Reported: Mean time across iterations

Memory Measurement

import tracemalloc

tracemalloc.start()
result = leanfe(...)
current, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()
  • Metric: Peak memory usage in MB
  • Note: Measures Python heap only, not total process memory

Validation

Coefficient Matching

All benchmarks verify that packages produce matching coefficients:

tolerance = 1e-6

polars_coef = result_polars['coefficients']['treatment']
duckdb_coef = result_duckdb['coefficients']['treatment']

assert abs(polars_coef - duckdb_coef) < tolerance

Ground Truth Comparison

With known true coefficients (treatment=2.5, x1=1.5), we verify estimates are close:

import numpy as np
import polars as pl
from leanfe import leanfe

# Generate data with known coefficients
np.random.seed(42)
n = 100_000

fe1 = np.random.randint(1, 501, n)
fe2 = np.random.randint(1, 101, n)
fe1_effects = np.random.normal(0, 1, 501)[fe1]
fe2_effects = np.random.normal(0, 0.5, 101)[fe2]

treatment = np.random.binomial(1, 0.3, n).astype(float)
x1 = np.random.normal(0, 1, n)

# True coefficients: treatment=2.5, x1=1.5
y = 2.5 * treatment + 1.5 * x1 + fe1_effects + fe2_effects + np.random.normal(0, 1, n)

df = pl.DataFrame({"y": y, "treatment": treatment, "x1": x1, "fe1": fe1, "fe2": fe2})

result = leanfe(formula="y ~ treatment + x1 | fe1 + fe2", data=df, vcov="iid")

print("Ground Truth Validation:")
print(f"  True treatment: 2.5, Estimated: {result['coefficients']['treatment']:.4f}")
print(f"  True x1: 1.5, Estimated: {result['coefficients']['x1']:.4f}")
Ground Truth Validation:
  True treatment: 2.5, Estimated: 2.4996
  True x1: 1.5, Estimated: 1.4924

Limitations and Caveats

What These Benchmarks Show

✅ Relative performance between leanfe backends ✅ Scaling behavior with dataset size ✅ Memory efficiency of DuckDB backend

What These Benchmarks Don’t Show

❌ Comparison with PyFixest/fixest (different feature sets) ❌ Performance on real-world data (may vary) ❌ GPU acceleration (not implemented)

Known Limitations

  1. Synthetic data: Real data may have different characteristics
  2. Single machine: Results may vary on different hardware
  3. Python overhead: Timing includes Python interpreter overhead
  4. Memory measurement: tracemalloc may underestimate total memory

Reproducible Scripts

Download Benchmark Code

The complete benchmark suite is available in the repository:

# Clone repository
git clone https://github.com/diegogentilepassaro/leanfe.git

# Navigate to benchmarks
cd leanfe/package/docs/benchmarks

# Run benchmarks
python benchmark_runner.py

Quick Benchmark

from leanfe import leanfe
import numpy as np
import polars as pl
import time

# Generate data
np.random.seed(42)
n = 1_000_000
df = pl.DataFrame({
    "y": np.random.normal(0, 1, n),
    "x": np.random.normal(0, 1, n),
    "fe1": np.random.randint(1, 1001, n),
    "fe2": np.random.randint(1, 501, n),
})

# Benchmark Polars
start = time.time()
leanfe(data=df, formula="y ~ x | fe1 + fe2", backend="polars")
print(f"Polars: {time.time() - start:.2f}s")

# Benchmark DuckDB
start = time.time()
leanfe(data=df, formula="y ~ x | fe1 + fe2", backend="duckdb")
print(f"DuckDB: {time.time() - start:.2f}s")

See Also