within

within provides high-performance solvers for projecting out high-dimensional fixed effects from regression problems.

By the Frisch-Waugh-Lovell theorem, estimating a regression of the form y = Xβ + Dα + ε reduces to a sequence of least-squares projections, one for y and one for each column of X, followed by a cheap regression fit on the resulting residuals. The projection step of solving the normal equations D'Dx = D'z is the computational bottleneck, which is the problem within is designed to solve.

within's solvers are tailored to the structure of fixed effects problems, which can be represented as a graph (as first noted by Correia, 2016). Concretely, within uses modified LSMR with a domain decomposition (Schwarz) preconditioner, backed by approximate Cholesky local solvers (Gao et al, 2025).

Scope

within is a low-level fixed-effects kernel. Callers pass categorical labels as uint32 F-order (column-major) arrays. Labels need not be contiguous: observed labels are compacted internally, while result layouts retain the caller-visible labels. Formula-level convenience — DataFrames, string/object categoricals, pandas.factorize, and formula parsing — is intentionally out of scope and belongs to a frontend layer built on top. The pyfixest-style workflow is served by such a frontend calling within underneath.

Installation

You can install Python bindings from PyPi by running

pip install within_py

Wheels ship for CPython 3.9+ (one abi3 wheel per platform), free-threaded CPython 3.14+ (cp314t; PyO3 dropped experimental 3.13t support), and PyPy 3.11 (PyO3's minimum supported PyPy version; earlier PyPy releases are rejected at build time). CPython 3.15's abi3t (PEP 803, PyO3/maturin#3064) will let the free-threaded build fold back into a single wheel per platform.

musllinux ships no PyPy wheel at all — the musllinux_1_2 build image bundles no PyPy interpreter.

Python Quickstart

within's main user-facing function is solve. Provide a 2-D uint32 array of category codes (one column per fixed-effect factor) and a response vector y. The solver finds x in the normal equations D'D x = D'y, where D is the sparse categorical design matrix.

import numpy as np
from within import solve, solve_batch, LsmrOptions, PreconditionerConfig

np.random.seed(1)
n = 100_000
fe = np.asfortranarray(np.column_stack([
    np.random.randint(0, 500, n).astype(np.uint32),
    np.random.randint(0, 200, n).astype(np.uint32),
]))
y = np.random.randn(n)

# Default: the adaptive diagonal→Schwarz ladder + LSMR
result = solve(fe, y)

# Custom tolerance / iteration cap
result = solve(fe, y, options=LsmrOptions(tol=1e-10, maxiter=2000))

# Weighted solve
result = solve(fe, y, weights=np.ones(n))

# Opt into diagonal/Jacobi preconditioning
result = solve(fe, y, preconditioner=PreconditionerConfig.Diagonal())

FWL regression example

beta_true = np.array([1.0, -2.0, 0.5])
X = np.random.randn(n, 3)
y = X @ beta_true + np.random.randn(n)

result = solve_batch(fe, np.column_stack([y, X]))
y_tilde, X_tilde = result.demeaned[:, 0], result.demeaned[:, 1:]
beta_hat = np.linalg.lstsq(X_tilde, y_tilde, rcond=None)[0]
print(np.round(beta_hat, 4))  # [ 0.9982 -2.006   0.5005]

Varying slopes

Pass a list of Effect terms instead of a categories array. Each term is a factor's level codes plus an optional intercept and zero or more slope covariates (per-level slopes, as in fixest's f[z] notation).

from within import solve, Effect

firm = np.random.randint(0, 500, n).astype(np.uint32)
year = np.random.randint(0, 20, n).astype(np.uint32)
x = np.random.randn(n)  # covariate whose slope varies by firm

result = solve(
    [
        Effect(firm, intercept=True, slopes=[x]),  # firm intercept + firm-specific x slope
        Effect(year, intercept=True),              # year intercept
    ],
    y,
)

# Read firm level 3's x-slope via the layout map (column 0 = intercept, 1 = first slope):
i = result.layout.index(0, 3, 1)
print(result.x[i])

Python API

High-level functions

design is a persistent Design, a 2-D uint32 array of shape (n_obs, n_factors), or a list of Effect terms (see Varying slopes). A UserWarning is emitted when a C-contiguous categories array is passed — use np.asfortranarray(design) for best performance.

Persistent solver

For repeated solves with the same design matrix, Design reuses the preprocessed observation storage and Solver builds the weight-dependent state and preconditioner once.

from within import Design, Solver, solve

design = Design(fe)
solver = Solver(design)
r = solver.solve(y)                            # reuses preconditioner
r = solver.solve_batch(np.column_stack([y, X]))

precond = solver.preconditioner                # picklable property
solver2 = Solver(design, preconditioner=precond)  # skip re-factorization

# Reuse design preprocessing while rebuilding weight-dependent solver state.
weighted = solve(design, y, weights=np.ones(n))

Solver configuration

Preconditioners

The preconditioner argument accepts any of:

PreconditionerConfig is a tagged union: each variant is a subclass, so instances support match/case and compare by value. A built preconditioner exposes the configuration used to construct it as .config and its original build duration in seconds as .build_duration_seconds; both are preserved across pickling and reuse.

Local solver configuration (advanced — within.config)

Result types

SolveResult: x (coefficients), unidentified (directions the data cannot identify, as UnidentifiedDirection(term, level, column) records), layout (a CoefficientLayout mapping a (term, level, column) address to its flat x index and back), demeaned (residuals), converged, iterations, residual, time_total, time_setup, time_solve.

BatchSolveResult: Same fields, with converged, iterations, residual, and time_solve as lists (one entry per RHS).

Coefficients for unidentified directions are pinned to the minimal-norm value 0 (never NaN). This is why x can differ from reference tools that instead drop a reference level; the identified fit — demeaned — is unaffected by the choice.

Rust API

use ndarray::Array2;
use within::{solve, LsmrOptions, PreconditionerConfig};
use within::config::{LocalSolverConfig, ReductionStrategy};

let categories = /* Array2<u32> of shape (n_obs, n_factors) */;
let y: &[f64] = /* response vector */;

// Default: LSMR + the adaptive ladder (None → library default)
let r = solve(categories.view(), &y, None, &LsmrOptions::default(), None)?;
assert!(r.converged);

// Tighter tolerance with an explicit additive preconditioner
let lsmr = LsmrOptions { tol: 1e-10, ..LsmrOptions::default() };
let precond = PreconditionerConfig::Additive {
    local_solver: LocalSolverConfig::default(),
    reduction: ReductionStrategy::default(),
};
let r = solve(categories.view(), &y, None, &lsmr, &precond)?;

// Opt into diagonal/Jacobi preconditioning
let diagonal = PreconditionerConfig::Diagonal;
let r = solve(categories.view(), &y, None, &lsmr, &diagonal)?;

Persistent solver — build once, solve many:

use ndarray::array;
use within::{solve, Channel, CoefficientAddress, Design, Solver};

// Caller labels are deliberately non-contiguous.
let categories = array![
    [10_u32, 100],
    [20, 100],
    [10, 900],
    [20, 900],
];
let y = [1.0, 2.0, 3.0, 4.0];

// Build and compact the design once, then share it across solvers.
let design = Design::from_categories(categories.view())?;
let unweighted = Solver::new(&design, None, None)?;
let weights = [1.0, 2.0, 1.0, 2.0];
let weighted = Solver::new(&design, Some(&weights), None)?;

let result = unweighted.solve(&y, None)?;
let weighted_result = weighted.solve(&y, None)?;
let one_shot_result = solve(&design, &y, None, None, None)?;
assert!(weighted_result.converged && one_shot_result.converged);

// Internal factor positions stay private; result lookup uses caller labels.
let address = CoefficientAddress {
    channel: Channel { term: 1, column: 0 },
    level: 900,
};
let slot = result.layout.index(address).expect("observed label");
assert_eq!(result.layout.address(slot), Some(address));
assert_eq!(result.x.len(), 4); // two observed levels in each of two factors

solve and Solver::new take the preconditioner as impl Into<PreconditionerInput>: None (library default), a &PreconditionerConfig or owned PreconditionerConfig (e.g. PreconditionerConfig::Off for the identity), or an owned/borrowed Preconditioner for reuse. LSMR options are impl Into<Option<&LsmrOptions>>, so None accepts the defaults and &opts overrides them.

Varying slopes

Pass a Vec<Effect> instead of a categories array. Each term is a factor's level codes plus an optional intercept and zero or more slope covariates (per-level slopes, as in fixest's f[z] notation).

use within::{solve, Channel, CoefficientAddress, Effect};

let firm: &[u32] = /* level codes */;
let year: &[u32] = /* level codes */;
let x: &[f64] = /* covariate whose slope varies by firm */;

let terms = vec![
    Effect::new(firm, true, [x])?,  // firm intercept + firm-specific x slope
    Effect::new(year, true, [])?,   // year intercept
];
let r = solve(terms, &y, None, None, None)?;

// Read firm level 3's x-slope via the layout map (column 0 = intercept, 1 = first slope):
let at = CoefficientAddress { channel: Channel { term: 0, column: 1 }, level: 3 };
println!("{}", r.x[r.layout.index(at).unwrap()]);

Solver::new takes the same Vec<Effect>, so a slope design can be reused across solves like any other.

Lower-level access

Project structure

crates/
  schwarz-precond/   Generic domain decomposition library (traits, solvers, Schwarz preconditioners)
  within/            Core fixed-effects solver (designs, domains, operators, solver)
  within-py/         PyO3 bridge (cdylib → within._within)
python/within/       Python package re-exporting the Rust extension
benchmarks/          Python benchmark framework

Development

Uses pixi as the task runner.

pixi run develop                         # Build Rust extension (release mode)
pixi run test                            # Rebuild + pytest
cargo test --workspace                   # Rust tests only
cargo bench -p within                    # Criterion benchmarks
pixi run python -m benchmarks run all    # Python benchmarks

Rust changes require rebuilding before running Python code (pixi run develop).

License

MIT

References

• Correia, Sergio. "A feasible estimator for linear models with multi-way fixed effects." Preprint at http://scorreia.com/research/hdfe.pdf (2016).
• Gao, Y., Kyng, R. & Spielman, D. A. (2025). AC(k): Robust Solution of Laplacian Equations by Randomized Approximate Cholesky Factorization. SIAM Journal on Scientific Computing.
• Toselli & Widlund (2005). Domain Decomposition Methods — Algorithms and Theory. Springer.
• Xu, J. (1992). Iterative Methods by Space Decomposition and Subspace Correction. SIAM Review, 34(4), 581--613.