Running in Colab? Install first: %pip install beamfeat.

See docs/installation.md for the other setups.

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

warnings.filterwarnings("ignore")

from beamfeat import BeamFeatRegressor

rng = np.random.default_rng(0)

import warnings import numpy as np 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

X = np.column_stack([a, b, c, 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 X = np.column_stack([a, b, c, 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

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 [4]:

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

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

  ((x1 / x2) * x0)

In [5]:

print(model.equation(max_terms=3))

print(model.equation(max_terms=3))

y = 0.9971*((x1 / x2) * x0) + 0.0129

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

Comparing against a linear baseline¶

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

In [6]:

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.9998

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 [7]:

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 [8]:

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(x0) * log(x1))
  ((x0 + x2) + (x0 * x1))
  ((x0 + x2) + (x1 - x2))
  ((x0 * x1) - (x0 + x2))
  (log(x1) * (x0)^2)
  (log(x0) * (x1)^2)
  ((x0 * x1))^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.