beamfeat 0.1.0

Beam-search feature construction with false-discovery-rate controlled selection

beamfeat

Beam-search feature construction with false-discovery-rate (FDR) controlled selection, packaged as scikit-learn estimators.

beamfeat builds interpretable mathematical expressions from your columns — (x0 / x2) * x1, log(a) * log(b) — searches for the ones that explain the target, and then screens them with a selection procedure whose statistical guarantee is stated, tested, and honestly reported when it cannot be given.

from beamfeat import BeamFeatRegressor

model = BeamFeatRegressor(max_depth=2, beam_width=30).fit(X, y)

model.formulas()        # ['((x0 / x2) * x1)', ...]
model.equation()        # 'y = 0.9847*((x0 / x2) * x1) + ...'
model.fdr_controlled_   # True: the features carry the stated FDR guarantee
model.predict(X_new)

Install

pip install beamfeat            # core: numpy + scikit-learn only
# before the first PyPI release, install from source:
#   pip install git+https://github.com/LD-Shell/beamfeat
pip install "beamfeat[units]"   # + pint, for unit-aware search (units as
                                #   strings like "kg" or as pint quantities)
pip install "beamfeat[all]"     # + sympy display support and test deps

Python ≥ 3.10. There are no upper version pins: the test suite runs against the oldest supported stack (numpy 1.26 / scikit-learn 1.6) and through the current releases (verified up to scikit-learn 1.9 / numpy 2.5).

What it does

1. Search. A beam search over an expression DAG proposes and scores compositions of unary (log, sqrt, reciprocal, square, abs) and binary (*, /, +, -) operators, with structural deduplication, redundancy pruning, and cost O(depth × beam²) rather than the exhaustive O(operators^depth).
2. Selection. Candidates are tested for association with the target under procedures that control the false discovery rate (details below), on a held-out split of the training rows that the search never saw.
3. Fit. A parsimony step (greedy forward selection within the screened set, on by default) keeps the compact subset the equation uses — the full screened set with per-candidate p- and q-values stays auditable in selection_report_ — then a linear model fits it, readable via .equation().

Numerical failures (overflow, domain errors) are recorded and excluded, never silently masked. Optional unit propagation — units given as pint quantities or as plain strings like "kg" or "m / s**2" — rejects dimensionally invalid expressions at construction time — on a mechanics example this pruned 62% of proposals before any data work.

The guarantees, stated precisely

• Permutation selector (default). Tests marginal association with an exact permutation test: the statistic (|Pearson r| for regression, eta-squared for classification) is a fixed function of the data, and p-values use the add-one estimator of Phipson & Smyth (2010). Multiplicity is corrected with Benjamini–Yekutieli by default (valid under arbitrary dependence); Benjamini–Hochberg is available and is valid under positive dependence (PRDS).
• Knockoff selector. Fixed-X knockoffs (Barber & Candès, 2015) when n ≥ 2p: no distributional assumption on the features — engineered columns are fine — with the guarantee requiring Gaussian noise in the linear model. Model-X Gaussian knockoffs (Candès et al., 2018) only as an n < 2p fallback, with a warning, since engineered features violate its Gaussianity assumption.
• Selective inference. Search retains candidates because they correlate with the target in its sample, so testing them on the same rows would bias p-values optimistically. By default the estimators search on one half of the training rows and select on the other (selection_holdout=0.5).
• Association ≠ generalisation. A post-fit check on the selection holdout raises a visible DegenerateFitWarning if the assembled model fails to generalise (negative held-out R², or accuracy below the majority class) — FDR certifies the features' association, and this check covers the separate claim that the equation built on them is sound.
• Honest failure. If nothing passes selection, the default (on_no_discoveries="empty") keeps no constructed features: the model degrades to an intercept-only fit and raises a visible NoDiscoveriesWarning. "fallback" (keep the unfiltered search output, flagged fdr_controlled_ = False) and "raise" are available. It never silently claims a guarantee it does not have.

Three readings of these guarantees matter. FDR is an expectation over repetitions, not a per-run bound. The marginal null counts a near-duplicate of a real signal as a true discovery — de-duplication is the redundancy threshold's job, not the error-control procedure's. And every number below carries Monte Carlo error from finite trials.

Measured, not assumed

All of the following are measured by the test suite or the committed benchmark scripts, not asserted:

• Selector calibration (Gaussian designs, 25 trials): realised FDR 0.065/0.120/0.210 at nominal 0.05/0.10/0.20 for BH, 0.020/0.033/0.065 for BY, power 1.00 throughout. Fixed-X knockoff+ realised 0.201 at nominal 0.2.
• End-to-end pipeline (search + holdout selection) at nominal FDR 0.10: empirical FDR 0.0000 with power 1.000 over 200 replicates, zero fallbacks; zero selections over 60 global-null replicates.
• Feynman-equation panel (12 physics laws, 0.1% noise): 9–10/12 solved (a borderline surrogate flips across numeric stacks) at the SRBench criterion and a stable 8/12 in exact symbolic form, ~0.9 s per equation; the misses are named boundaries (a depth-3 rational, a literal constant, depth-4 nesting, the Gaussian's exponential).
• Pure-noise stress: BH breached in 3/25 trials (mean FDP 0.12 vs nominal 0.10, the PRDS caveat materialising on correlated candidates); BY in 1/25 (mean FDP 0.04). This measurement is why BY is the estimator default.
• Stress suite (signal among ten distractor columns): beamfeat's false-feature rate — returned formulas touching only irrelevant columns — is 0.000 at noise levels from 5% to 50% of the signal scale and at n as small as 120, with recovery and the guarantee intact throughout; autofeat recovers the formula equally often but returns distractor-only features at a mean rate of 0.29 (worst 0.75 at 50% noise).
• Real data (six datasets, numeric columns): best on mpg (0.871 vs LightGBM 0.853) and tips — where gradient boosting overfits below ridge while the FDR gate holds beamfeat at the ceiling — tied with ridge on diamonds from one interpretable feature, behind LightGBM on penguins and breast cancer, and behind ridge on diabetes (0.376; autofeat 0.303, OpenFE 0.235): no nonlinear structure exists there to find, and beamfeat does not invent any.
• OpenFE on the same stress suite (its features + LightGBM, after shimming a removed sklearn parameter): false-feature rate also 0.000, but 0/5 formulas recovered and lower R² throughout — it optimises tree-model lift, not equations.
• Heavy-tailed noise (Student t(2), infinite variance): exact formula recovered at R² 0.956 vs LightGBM 0.866 — the permutation test is distribution-free, so heavy tails cost it nothing.
• Known boundaries, measured: on piecewise targets trees win outright (LightGBM 0.998 vs 0.808 — if a tree model beats beamfeat by a wide margin, the signal is likely piecewise, not algebraic), and on Friedman #1 the 0.744-vs-0.964-oracle gap splits into a measured 0.875 screening-admissible ceiling (the marginally-quiet quadratic is priced by the guarantee itself) and a search shortfall below it that is insensitive to beam width and diversity settings — stated rather than hidden.
• Test suite: 370 tests at 95% statement coverage (a 90% floor is enforced in CI), including scikit-learn's full estimator-conformance checks.
• Selection-only baseline (knockpy on raw columns, modified-FDR offset — the only satisfiable configuration at these dimensionalities): 0/15 formulas recovered, core-suite R² at the ridge level (0.845), and a 0.11 mean false-feature rate on the stress suite — the construction step and the holdout are the delta.
• Independent 315-fit study across nine datasets and seven methods (benchmarks/independent/): mean R² 0.803 against random forest 0.810, LightGBM 0.798, autofeat 0.751, OpenFE 0.736, featuretools −2.478 — and a worst single fit of 0.355 where autofeat reached −2.28 and featuretools −57.3, at ~48× autofeat's speed, with the FDR flag delivered on 45/45 fits.
• Benchmarks (10 synthetic problems with known formulas, 70/30 splits): beamfeat mean R² 0.9993, ~0.4 s fits, 9/9 formulas recovered; autofeat 0.9974, ~9 s, 8/9; LightGBM 0.9798; ridge 0.8442. On the purely linear problem — where feature construction should not win — ridge ties beamfeat.

Relationship to autofeat

beamfeat's design began from a code review of autofeat (Horn et al., 2019) and keeps its best idea — compact symbolic features feeding a linear model — while replacing exhaustive expansion with beam search, heuristic noise-threshold selection with FDR-controlled procedures, and upper-pinned dependencies with tested floors. The generate-then-FDR-filter design has a direct precedent in tsfresh (Christ et al., 2018), which filters mass-generated time-series features under Benjamini–Yekutieli control; beamfeat applies that design to symbolic expressions on tabular columns.

Independent comparison

benchmarks/independent/ holds a 315-fit study (nine datasets, seven methods, five splits, average ranks plus a Friedman test) with its harness, pinned environment, data, and a notebook that reproduces the analysis. Mean held-out R²: beamfeat 0.803, random forest 0.810, LightGBM 0.798, autofeat 0.751, OpenFE 0.736, ridge 0.704, featuretools −2.478 — and worst single fit of 45: beamfeat 0.355, autofeat −2.28, featuretools −57.3, with the FDR guarantee delivered on 45/45. Read PROVENANCE.md first: it records what is reproducible (beamfeat, bit-identically) and what is not (autofeat seeds nothing internally and swung from +0.955 to −77.86 on one identical split).

Documentation

Tutorial notebooks live in notebooks/ (getting started; search and scoring; selection and units). The statistical guarantees are documented in detail in docs/guarantees.md and in the module docstrings of beamfeat.selection. Users coming from autofeat will find a full API mapping and version-compatibility notes in docs/migrating-from-autofeat.md.

Citing

See CITATION.cff. If you use the selection procedures, please also cite the underlying statistics: Benjamini & Hochberg (1995), Benjamini & Yekutieli (2001), Barber & Candès (2015), Candès et al. (2018), and Phipson & Smyth (2010).

License

MIT.