Running in Colab? Install first: %pip install beamfeat (before the PyPI release: %pip install git+https://github.com/LD-Shell/beamFeat).

Getting started with beamfeat¶

beamfeat constructs interpretable features from tabular data. Given columns a, b, and c, it searches expressions like log(a), a * b, and (a * b) / c, keeps the ones that predict the target, and reports them as readable formulas.

Two things distinguish it from exhaustive feature engineering:

1. The search is guided. Candidates are scored against the target and only the best are extended, so cost does not compound with expression depth.
2. Selection is calibrated. The retained features carry a false discovery rate guarantee rather than surviving a heuristic threshold.

This notebook covers the basic workflow. Later notebooks cover the search and selection machinery in detail.

In [1]:

import warnings

import numpy as np
import pandas as pd

warnings.filterwarnings("ignore")

from beamfeat import BeamFeatRegressor

rng = np.random.default_rng(0)

import warnings import numpy as np import pandas as pd warnings.filterwarnings("ignore") from beamfeat import BeamFeatRegressor rng = np.random.default_rng(0)

A problem with a known answer¶

We generate data where the target is a known function of the inputs, so we can check whether beamfeat recovers it. The relationship is y = (a * b) / c, which no linear model on the raw columns can represent.

In [2]:

n = 500
a = rng.uniform(1.0, 6.0, n)
b = rng.uniform(1.0, 6.0, n)
c = rng.uniform(1.0, 6.0, n)
d = rng.uniform(1.0, 6.0, n)  # an irrelevant column

# A DataFrame rather than a bare array: beamfeat takes column names from it,
# so the formulas come back as `(a * b) / c` instead of `(x0 * x1) / x2`.
# A NumPy array works identically and reports columns as x0...xp.
X = pd.DataFrame({"a": a, "b": b, "c": c, "d": d})
y = (a * b) / c + rng.normal(0, 0.05, n)

print(f"{X.shape[0]} rows, {X.shape[1]} columns")
print(f"target range: {y.min():.2f} to {y.max():.2f}")

n = 500 a = rng.uniform(1.0, 6.0, n) b = rng.uniform(1.0, 6.0, n) c = rng.uniform(1.0, 6.0, n) d = rng.uniform(1.0, 6.0, n) # an irrelevant column # A DataFrame rather than a bare array: beamfeat takes column names from it, # so the formulas come back as `(a * b) / c` instead of `(x0 * x1) / x2`. # A NumPy array works identically and reports columns as x0...xp. X = pd.DataFrame({"a": a, "b": b, "c": c, "d": d}) y = (a * b) / c + rng.normal(0, 0.05, n) print(f"{X.shape[0]} rows, {X.shape[1]} columns") print(f"target range: {y.min():.2f} to {y.max():.2f}")

500 rows, 4 columns
target range: 0.35 to 25.99

Fitting¶

BeamFeatRegressor follows the usual scikit-learn interface. The two parameters that matter most are max_depth, which bounds expression complexity, and beam_width, which bounds how many expressions survive each depth.

In [3]:

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

print(f"R^2: {model.score(X, y):.4f}")
print(f"features constructed: {model.n_features_out_}")

model = BeamFeatRegressor(max_depth=2, beam_width=30, random_state=0) model.fit(X, y) print(f"R^2: {model.score(X, y):.4f}") print(f"features constructed: {model.n_features_out_}")

R^2: 0.9999
features constructed: 1

Watching it work¶

verbose is a level. 0 is silent, 1 prints a line per stage, and 2 adds per-depth search detail and the strongest certified candidates with their p- and q-values. It is the quickest way to see where a fit spent its effort, and to tell a search that found nothing from a screen that certified nothing.

In [4]:

BeamFeatRegressor(max_depth=2, beam_width=30, random_state=0, verbose=1).fit(X, y)

BeamFeatRegressor(max_depth=2, beam_width=30, random_state=0, verbose=1).fit(X, y)

[beamfeat] fit: 500 rows x 4 columns, split 250 search / 250 selection
[beamfeat] search: depth 2, beam 30 -> 3094 proposed, 30 candidates (0.20s)
[beamfeat] screen: 30 candidates, permutation (BY) at q=0.1 -> 30 certified
[beamfeat] parsimony: 30 certified -> 1 term (|S|/|S'| = 30.00)
[beamfeat] result set: 1 term in the fitted equation (the guarantee covers the certified set, not this subset)
[beamfeat] result: y = 0.9991*((b / c) * a) + 0.003388   [1 of 30 certified terms; parsimony=None prints all 30]  [FDR controlled over the screened set]

Out[4]:

BeamFeatRegressor(beam_width=30, verbose=1)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Reading what it found¶

This is the part a gradient-boosted model cannot give you. The fitted model exposes both the individual feature formulas and the full equation.

In [5]:

for formula in model.formulas()[:5]:
    print(" ", formula)

for formula in model.formulas()[:5]: print(" ", formula)

  ((b / c) * a)

In [6]:

print(model.equation())

print(model.equation())

y = 0.9991*((b / c) * a) + 0.003388   [1 of 30 certified terms; parsimony=None prints all 30]

The true generating expression should appear among the selected features. Note that d — the irrelevant column — should be largely absent.

The bracketed suffix is the equation telling you that it is a subset. Screening certified a larger set; the parsimony step kept the compact predictive part of it and fitted that, so the guarantee covers the set the terms came from rather than the terms themselves. Notebook 04 takes that apart.

Comparing against a linear baseline¶

The point of constructing features is to let a simple model fit a relationship it otherwise could not.

In [7]:

from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)

baseline = Ridge().fit(X_train, y_train)
beamfeat = BeamFeatRegressor(max_depth=2, beam_width=30, random_state=0).fit(X_train, y_train)

print(f"Ridge on raw columns:  R^2 = {baseline.score(X_test, y_test):.4f}")
print(f"beamfeat features:     R^2 = {beamfeat.score(X_test, y_test):.4f}")

from sklearn.linear_model import Ridge from sklearn.model_selection import train_test_split X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0) baseline = Ridge().fit(X_train, y_train) beamfeat = BeamFeatRegressor(max_depth=2, beam_width=30, random_state=0).fit(X_train, y_train) print(f"Ridge on raw columns: R^2 = {baseline.score(X_test, y_test):.4f}") print(f"beamfeat features: R^2 = {beamfeat.score(X_test, y_test):.4f}")

Ridge on raw columns:  R^2 = 0.7489
beamfeat features:     R^2 = 0.9999

Using it in a pipeline¶

BeamFeatTransformer constructs features without fitting a model, so it composes with any downstream estimator. Because construction happens inside fit, cross-validating the pipeline does not leak information across folds.

In [8]:

from sklearn.model_selection import cross_val_score
from sklearn.pipeline import Pipeline

from beamfeat import BeamFeatTransformer

pipeline = Pipeline(
    [
        ("features", BeamFeatTransformer(max_depth=2, beam_width=20, random_state=0)),
        ("model", Ridge()),
    ]
)

scores = cross_val_score(pipeline, X, y, cv=5)
print(f"cross-validated R^2: {scores.mean():.4f} (+/- {scores.std():.4f})")

from sklearn.model_selection import cross_val_score from sklearn.pipeline import Pipeline from beamfeat import BeamFeatTransformer pipeline = Pipeline( [ ("features", BeamFeatTransformer(max_depth=2, beam_width=20, random_state=0)), ("model", Ridge()), ] ) scores = cross_val_score(pipeline, X, y, cv=5) print(f"cross-validated R^2: {scores.mean():.4f} (+/- {scores.std():.4f})")

cross-validated R^2: 0.9999 (+/- 0.0000)

Classification¶

BeamFeatClassifier works the same way, with predict_proba and decision_function available as usual.

In [9]:

from beamfeat import BeamFeatClassifier

labels = ((a * b) > np.median(a * b)).astype(int)

classifier = BeamFeatClassifier(
    max_depth=2, beam_width=20, selector="permutation", target_fdr=0.1, random_state=0
)
classifier.fit(X, labels)

print(f"accuracy: {classifier.score(X, labels):.4f}")
print(f"features retained after FDR control: {classifier.n_features_out_}")
for formula in classifier.formulas():
    print(" ", formula)

from beamfeat import BeamFeatClassifier labels = ((a * b) > np.median(a * b)).astype(int) classifier = BeamFeatClassifier( max_depth=2, beam_width=20, selector="permutation", target_fdr=0.1, random_state=0 ) classifier.fit(X, labels) print(f"accuracy: {classifier.score(X, labels):.4f}") print(f"features retained after FDR control: {classifier.n_features_out_}") for formula in classifier.formulas(): print(" ", formula)

accuracy: 0.9920
features retained after FDR control: 7
  (log(a) * log(b))
  ((a + c) + (a * b))
  ((a + c) + (b - c))
  ((a * b) - (a + c))
  (log(b) * (a)^2)
  (log(a) * (b)^2)
  ((a * b))^2

What to read next¶

• 02: Search and scoring — how the beam search works, and how the three scoring strategies differ in what they detect and what they cost.
• 03: Selection and units — how false discovery rate control works, why the default is not knockoffs, and how dimensional analysis constrains the search.
• 04: Reading a fit — what to do when the equation contains columns you did not expect, and how to tell a search failure from a selection failure.