Running in Colab? Install first: %pip install "beamfeat[units]".

The units extra pulls in pint, which the dimensional-analysis section below needs. See docs/installation.md for the other setups.

Selection and units¶

Feature construction proposes far more candidates than any dataset can support. Selection is where spurious features are either excluded or silently admitted, and it is where most automated feature engineering tools are weakest: an importance threshold gives no guarantee about how many reported features are noise, and the expected number of false positives grows with the size of the candidate pool.

beamfeat controls the false discovery rate — the expected proportion of selected features that are spurious.

In [1]:

import warnings

import numpy as np

warnings.filterwarnings("ignore")

from beamfeat.selection import KnockoffSelector, PermutationSelector

rng = np.random.default_rng(0)

import warnings import numpy as np warnings.filterwarnings("ignore") from beamfeat.selection import KnockoffSelector, PermutationSelector rng = np.random.default_rng(0)

The problem, demonstrated¶

With a pure-noise target, every selected feature is by definition a false positive. A well-calibrated selector should return almost nothing.

In [2]:

n_trials = 15
selections = []

for trial in range(n_trials):
    trial_rng = np.random.default_rng(100 + trial)
    features = trial_rng.standard_normal((200, 30))
    noise_target = trial_rng.standard_normal(200)

    result = PermutationSelector(target_fdr=0.1, n_permutations=20, random_state=trial).select(
        features, noise_target
    )
    selections.append(result.n_selected)

print(f"pure-noise target, {n_trials} trials, 30 candidate features each")
print(f"mean features selected: {np.mean(selections):.2f}")
print(f"trials selecting nothing: {sum(s == 0 for s in selections)}/{n_trials}")

n_trials = 15 selections = [] for trial in range(n_trials): trial_rng = np.random.default_rng(100 + trial) features = trial_rng.standard_normal((200, 30)) noise_target = trial_rng.standard_normal(200) result = PermutationSelector(target_fdr=0.1, n_permutations=20, random_state=trial).select( features, noise_target ) selections.append(result.n_selected) print(f"pure-noise target, {n_trials} trials, 30 candidate features each") print(f"mean features selected: {np.mean(selections):.2f}") print(f"trials selecting nothing: {sum(s == 0 for s in selections)}/{n_trials}")

pure-noise target, 15 trials, 30 candidate features each
mean features selected: 0.07
trials selecting nothing: 14/15

Measured calibration¶

More usefully: with real signal present, what fraction of selections are false? This is the quantity the library claims to control.

In [3]:

def gaussian_design(seed, n=300, n_signal=5, n_noise=20, effect=3.0):
    trial_rng = np.random.default_rng(seed)
    n_features = n_signal + n_noise
    features = trial_rng.standard_normal((n, n_features))
    coefficients = np.zeros(n_features)
    coefficients[:n_signal] = effect
    target = features @ coefficients + trial_rng.standard_normal(n)
    truth = np.zeros(n_features, dtype=bool)
    truth[:n_signal] = True
    return features, target, truth


def realised_fdr(selector_factory, nominal, n_trials=15):
    false_rates, recalls = [], []
    for trial in range(n_trials):
        features, target, truth = gaussian_design(3000 + trial)
        result = selector_factory(nominal, trial).select(features, target)
        if result.n_selected:
            false_rates.append(np.sum(~truth[result.selected]) / result.n_selected)
        else:
            false_rates.append(0.0)
        recalls.append(np.sum(truth[result.selected]) / truth.sum())
    return np.mean(false_rates), np.mean(recalls)


print(f"{'nominal':>8} {'realised FDR':>14} {'power':>8}")
for nominal in (0.05, 0.1, 0.2):
    fdr, power = realised_fdr(
        lambda f, s: PermutationSelector(target_fdr=f, random_state=s, n_permutations=25), nominal
    )
    print(f"{nominal:>8.2f} {fdr:>14.3f} {power:>8.2f}")

def gaussian_design(seed, n=300, n_signal=5, n_noise=20, effect=3.0): trial_rng = np.random.default_rng(seed) n_features = n_signal + n_noise features = trial_rng.standard_normal((n, n_features)) coefficients = np.zeros(n_features) coefficients[:n_signal] = effect target = features @ coefficients + trial_rng.standard_normal(n) truth = np.zeros(n_features, dtype=bool) truth[:n_signal] = True return features, target, truth def realised_fdr(selector_factory, nominal, n_trials=15): false_rates, recalls = [], [] for trial in range(n_trials): features, target, truth = gaussian_design(3000 + trial) result = selector_factory(nominal, trial).select(features, target) if result.n_selected: false_rates.append(np.sum(~truth[result.selected]) / result.n_selected) else: false_rates.append(0.0) recalls.append(np.sum(truth[result.selected]) / truth.sum()) return np.mean(false_rates), np.mean(recalls) print(f"{'nominal':>8} {'realised FDR':>14} {'power':>8}") for nominal in (0.05, 0.1, 0.2): fdr, power = realised_fdr( lambda f, s: PermutationSelector(target_fdr=f, random_state=s, n_permutations=25), nominal ) print(f"{nominal:>8.2f} {fdr:>14.3f} {power:>8.2f}")

 nominal   realised FDR    power
    0.05          0.067     1.00

    0.10          0.130     1.00
    0.20          0.223     1.00

Realised FDR tracks the nominal level, and power is full. Two reading notes, both scientifically important. First, FDR is an expectation: individual trials may exceed the level, and a mean over 15 trials carries visible Monte Carlo noise, so values a little above nominal are consistent with control. Second, an exact test spends its error budget — realised FDR near nominal is the signature of a calibrated procedure, where a realised FDR pinned at zero would signal wasted power.

Knockoffs: two constructions, two sets of assumptions¶

The knockoff filter comes in two forms, and it matters which one you mean.

Fixed-X knockoffs (Barber & Candès, 2015) treat the design as fixed and make no assumption about the distribution of the features — deterministic engineered columns are admissible. The guarantee needs n >= 2p and Gaussian noise in y = Xβ + ε. Its weakness on engineered designs is power, not validity: near-duplicate columns drive the construction's s toward zero, the knockoffs become nearly identical to the originals, and the filter loses the ability to tell them apart.

Model-X knockoffs (Candès et al., 2018) instead assume the features are jointly Gaussian — which engineered features are not: log(a) and a * b are deterministic functions of shared parents. beamfeat only uses this construction when n < 2p, where fixed-X does not exist, and warns.

beamfeat routes between the two automatically and reports which assumptions are under strain, rather than pretending otherwise.

In [4]:

def engineered_design(seed, n=300):
    trial_rng = np.random.default_rng(seed)
    a = trial_rng.uniform(1.0, 5.0, n)
    b = trial_rng.uniform(1.0, 5.0, n)
    c = trial_rng.uniform(1.0, 5.0, n)
    columns = [
        a, b, c, np.log(a), np.log(b), np.sqrt(a), 1.0 / a,
        a * b, a / b, b / a, a * c, a - b, a + b, a**2, b**2,
    ]
    return np.column_stack(columns), a * b + trial_rng.normal(0, 0.1, n)


features, target = engineered_design(5000)

knockoff_result = KnockoffSelector(target_fdr=0.1).select(features, target)
print("knockoff selector on an engineered design:")
for message in knockoff_result.warnings_raised:
    print(f"  warning: {message}")

permutation_result = PermutationSelector(target_fdr=0.1, n_permutations=25).select(features, target)
print(f"\npermutation selector warnings: {permutation_result.warnings_raised or 'none'}")

def engineered_design(seed, n=300): trial_rng = np.random.default_rng(seed) a = trial_rng.uniform(1.0, 5.0, n) b = trial_rng.uniform(1.0, 5.0, n) c = trial_rng.uniform(1.0, 5.0, n) columns = [ a, b, c, np.log(a), np.log(b), np.sqrt(a), 1.0 / a, a * b, a / b, b / a, a * c, a - b, a + b, a**2, b**2, ] return np.column_stack(columns), a * b + trial_rng.normal(0, 0.1, n) features, target = engineered_design(5000) knockoff_result = KnockoffSelector(target_fdr=0.1).select(features, target) print("knockoff selector on an engineered design:") for message in knockoff_result.warnings_raised: print(f" warning: {message}") permutation_result = PermutationSelector(target_fdr=0.1, n_permutations=25).select(features, target) print(f"\npermutation selector warnings: {permutation_result.warnings_raised or 'none'}")

[beamfeat.selection] design covariance is near-singular (min eigenvalue 0.00e+00); the knockoff construction remains valid for fixed-X but its power degrades toward zero, since knockoffs become nearly identical to the originals

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 3 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

knockoff selector on an engineered design:
  warning: design covariance is near-singular (min eigenvalue 0.00e+00); the knockoff construction remains valid for fixed-X but its power degrades toward zero, since knockoffs become nearly identical to the originals
  warning: knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 3 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

permutation selector warnings: none

So the practical division of labour: the permutation selector is the default because its exactness does not depend on the design at all; fixed-X knockoffs are a strong choice for wide-enough problems with roughly Gaussian noise (the features can be anything); model-X is a last resort for n < 2p.

One further caveat, measured rather than assumed: the offset parameter dominates knockoff power on narrow designs. offset=1 (knockoff+) requires (1 + #negatives) / #positives <= target_fdr, which cannot be satisfied when there are fewer than 1 / target_fdr features, regardless of signal strength.

In [5]:

print(f"{'offset':>8} {'power':>8}")
for offset in (0, 1):
    recalls = []
    for trial in range(10):
        features, target, truth = gaussian_design(6000 + trial)
        result = KnockoffSelector(target_fdr=0.1, random_state=trial, offset=offset).select(
            features, target
        )
        recalls.append(np.sum(truth[result.selected]) / truth.sum())
    print(f"{offset:>8} {np.mean(recalls):>8.2f}")

print(f"{'offset':>8} {'power':>8}") for offset in (0, 1): recalls = [] for trial in range(10): features, target, truth = gaussian_design(6000 + trial) result = KnockoffSelector(target_fdr=0.1, random_state=trial, offset=offset).select( features, target ) recalls.append(np.sum(truth[result.selected]) / truth.sum()) print(f"{offset:>8} {np.mean(recalls):>8.2f}")

  offset    power

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 5 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 5 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

       0     1.00

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 6 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 11 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 5 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 6 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

[beamfeat.selection] knockoff+ (offset=1) selected nothing at target FDR 0.1, but offset=0 — which controls only a modified FDR — would have selected 6 feature(s). If that trade-off is acceptable, pass KnockoffSelector(offset=0) as the selector; the permutation selector is the recommended default for engineered candidates

       1     0.30

End to end¶

Putting search and selection together: construct candidates, then keep only those that survive FDR control.

In [6]:

from beamfeat import BeamSearch, Evaluator

n = 500
data = {name: rng.uniform(1.0, 6.0, n) for name in "abcd"}
y = (data["a"] * data["b"]) / data["c"] + rng.normal(0, 0.05, n)

search_result = BeamSearch(max_depth=2, beam_width=40, random_state=0).run(data, y)
evaluator = Evaluator(data)
nodes, matrix = evaluator.evaluate_many(search_result.nodes)

selection = PermutationSelector(target_fdr=0.1, n_permutations=30).select(matrix, y)

print(f"proposed: {search_result.n_proposed_total}")
print(f"kept by search: {len(search_result)}")
print(f"kept by selection: {selection.n_selected}")
print("\nselected features:")
for index in selection.selected:
    print(f"  {nodes[index].name}")

from beamfeat import BeamSearch, Evaluator n = 500 data = {name: rng.uniform(1.0, 6.0, n) for name in "abcd"} y = (data["a"] * data["b"]) / data["c"] + rng.normal(0, 0.05, n) search_result = BeamSearch(max_depth=2, beam_width=40, random_state=0).run(data, y) evaluator = Evaluator(data) nodes, matrix = evaluator.evaluate_many(search_result.nodes) selection = PermutationSelector(target_fdr=0.1, n_permutations=30).select(matrix, y) print(f"proposed: {search_result.n_proposed_total}") print(f"kept by search: {len(search_result)}") print(f"kept by selection: {selection.n_selected}") print("\nselected features:") for index in selection.selected: print(f" {nodes[index].name}")

proposed: 3490
kept by search: 40
kept by selection: 40

selected features:
  ((a / c) * b)
  ((b / c) * (a - c))
  ((a / c) + sqrt(b))
  ((a / c) * (b - c))
  ((a * b) / sqrt(c))
  ((b / c) + log(a))
  ((a / c) * (b / c))
  ((a / c) + (b / c))
  ((a + d) * (b / c))
  ((b + d) * (a / c))
  (b + (a - c))
  ((b / c) + (a - c))
  ((a / c) + (b - c))
  ((a * b) * (a - c))
  ((a / c) - (c / b))
  (a / c)
  ((a / c) + b)
  ((a * b) / (c + d))
  ((b / c) + a)
  ((a * b) * (b - c))
  (b / c)
  ((a / c) - (d / b))
  ((a / c) - (a / b))
  ((b / c) - (b / a))
  ((a / c))^2
  ((d / c) * (a * b))
  ((b / c) - (d / a))
  ((a - c) - (d / b))
  ((b - c) - (d / a))
  ((a / c) + (b / d))
  ((b / c) / c)
  ((a + d) * (a / c))
  ((a / c) - sqrt(d))
  (b * (a - c))
  ((a / d) + (b / c))
  ((a / d) + (b - c))
  (a * (b - c))
  ((a - c) / (c + d))
  ((a / c) + sqrt(c))
  ((a / c) * sqrt(c))

Dimensional analysis¶

If your columns carry physical units, supplying them restricts the search to dimensionally valid expressions. Adding a mass to a length is rejected at construction time, before any numerical work happens.

This is both a correctness feature and a performance one: it prunes large parts of the search space for free.

In [7]:

import pint

ureg = pint.UnitRegistry()

n = 400
physical = {
    "mass": rng.uniform(1.0, 5.0, n),
    "length": rng.uniform(1.0, 5.0, n),
    "time": rng.uniform(1.0, 5.0, n),
}
units = {
    "mass": 1.0 * ureg.kilogram,
    "length": 1.0 * ureg.meter,
    "time": 1.0 * ureg.second,
}

momentum = physical["mass"] * physical["length"] / physical["time"]
physical_target = momentum + rng.normal(0, 0.05, n)

unconstrained = BeamSearch(max_depth=2, beam_width=25, random_state=0).run(physical, physical_target)
constrained = BeamSearch(max_depth=2, beam_width=25, random_state=0).run(
    physical, physical_target, units=units
)

print(f"without units: {unconstrained.n_proposed_total} candidates proposed")
print(f"with units:    {constrained.n_proposed_total} candidates proposed")
print(f"reduction:     {100 * (1 - constrained.n_proposed_total / unconstrained.n_proposed_total):.0f}%")

import pint ureg = pint.UnitRegistry() n = 400 physical = { "mass": rng.uniform(1.0, 5.0, n), "length": rng.uniform(1.0, 5.0, n), "time": rng.uniform(1.0, 5.0, n), } units = { "mass": 1.0 * ureg.kilogram, "length": 1.0 * ureg.meter, "time": 1.0 * ureg.second, } momentum = physical["mass"] * physical["length"] / physical["time"] physical_target = momentum + rng.normal(0, 0.05, n) unconstrained = BeamSearch(max_depth=2, beam_width=25, random_state=0).run(physical, physical_target) constrained = BeamSearch(max_depth=2, beam_width=25, random_state=0).run( physical, physical_target, units=units ) print(f"without units: {unconstrained.n_proposed_total} candidates proposed") print(f"with units: {constrained.n_proposed_total} candidates proposed") print(f"reduction: {100 * (1 - constrained.n_proposed_total / unconstrained.n_proposed_total):.0f}%")

without units: 1058 candidates proposed
with units:    366 candidates proposed
reduction:     65%

In [8]:

print("top features with units enforced:")
for name in constrained.names[:5]:
    print(f"  {name}")

invalid = [name for name in constrained.names if " + " in name or " - " in name]
print(f"\ndimensionally invalid sums or differences: {len(invalid)}")

print("top features with units enforced:") for name in constrained.names[:5]: print(f" {name}") invalid = [name for name in constrained.names if " + " in name or " - " in name] print(f"\ndimensionally invalid sums or differences: {len(invalid)}")

top features with units enforced:
  ((length / time) * mass)
  ((length / time) * (length * mass))
  ((mass / time) * (length * mass))
  ((length / time) * (mass / time))
  sqrt((length / time))

dimensionally invalid sums or differences: 0

Mass plus length never appears, because it was never constructible.

Units through the estimator API¶

The same constraint is available on the estimators via the units argument. The fitted equation is compact because, after FDR screening, the estimator keeps a parsimonious forward-selected subset of the screened features by default (parsimony="forward"); the full screened set with per-candidate p- and q-values is available in selection_report_.

In [9]:

from beamfeat import BeamFeatRegressor

# Rebuild from a fresh generator so this cell does not depend on how many
# draws earlier cells happened to consume.
cell_rng = np.random.default_rng(7)
physical = {
    "mass": cell_rng.uniform(1.0, 5.0, n),
    "length": cell_rng.uniform(1.0, 5.0, n),
    "time": cell_rng.uniform(1.0, 5.0, n),
}
physical_target = (
    physical["mass"] * physical["length"] / physical["time"] + cell_rng.normal(0, 0.05, n)
)

X_physical = np.column_stack([physical["mass"], physical["length"], physical["time"]])
# Units may be pint quantities or plain strings; strings are parsed with
# pint at fit time.
estimator_units = {"x0": "kg", "x1": "meter", "x2": "second"}

model = BeamFeatRegressor(
    max_depth=2, beam_width=25, units=estimator_units, selector="permutation",
    random_state=0,
).fit(X_physical, physical_target)

print(f"R^2: {model.score(X_physical, physical_target):.4f}")
print(f"features retained: {model.n_features_out_}")
print(model.equation())

from beamfeat import BeamFeatRegressor # Rebuild from a fresh generator so this cell does not depend on how many # draws earlier cells happened to consume. cell_rng = np.random.default_rng(7) physical = { "mass": cell_rng.uniform(1.0, 5.0, n), "length": cell_rng.uniform(1.0, 5.0, n), "time": cell_rng.uniform(1.0, 5.0, n), } physical_target = ( physical["mass"] * physical["length"] / physical["time"] + cell_rng.normal(0, 0.05, n) ) X_physical = np.column_stack([physical["mass"], physical["length"], physical["time"]]) # Units may be pint quantities or plain strings; strings are parsed with # pint at fit time. estimator_units = {"x0": "kg", "x1": "meter", "x2": "second"} model = BeamFeatRegressor( max_depth=2, beam_width=25, units=estimator_units, selector="permutation", random_state=0, ).fit(X_physical, physical_target) print(f"R^2: {model.score(X_physical, physical_target):.4f}") print(f"features retained: {model.n_features_out_}") print(model.equation())

R^2: 0.9997
features retained: 1
y = 0.9971*((x0 / x2) * x1) + 0.0138