Running in Colab? Install first: %pip install beamfeat.
See docs/installation.md for the other setups.
Search and scoring¶
This notebook covers how beamfeat explores the space of expressions, and how the choice of scorer changes what it finds.
The core problem is that expression space is enormous. With 5 input columns, 5 unary operators, and 4 binary operators, exhaustive expansion to depth 3 produces well over a million candidates, and each must be evaluated before it can be judged. Beam search avoids this by scoring candidates as they are generated and extending only the most promising.
import warnings
import numpy as np
warnings.filterwarnings("ignore")
from beamfeat import BeamSearch
rng = np.random.default_rng(0)
n = 500
data = {name: rng.uniform(1.0, 6.0, n) for name in "abcd"}
target = (data["a"] * data["b"]) / data["c"] + rng.normal(0, 0.05, n)
Running a search directly¶
BeamSearch is the layer beneath the estimators. Using it directly is useful
when you want to inspect the search rather than fit a model.
search = BeamSearch(max_depth=2, beam_width=30, random_state=0)
result = search.run(data, target)
print(result.summary())
beam search: 30 features from 3094 proposed / 3050 evaluated in 0.63s
depth 0: proposed 4, evaluated 4, kept 4, best 0.5941 (0.00s)
depth 1: proposed 56, evaluated 52, kept 30, best 0.8411 (0.02s)
depth 2: proposed 3038, evaluated 2994, kept 30, best 0.9999 (0.61s)
top features:
0.9999 ((a / c) * b)
0.9435 ((b / c) * (a - c))
0.9430 ((a / c) + sqrt(b))
0.9403 ((a / c) * (b - c))
0.9387 ((a * b) / sqrt(c))
0.9371 ((b / c) + log(a))
0.9350 ((a / c) * (b / c))
0.9302 ((a / c) + (b / c))
0.9201 ((a + d) * (b / c))
0.9109 ((b + d) * (a / c))
The trace¶
Every search records what happened at each depth. This is how you attribute cost, and how you notice a beam that collapsed or saturated.
print(f"{'depth':>6} {'proposed':>10} {'evaluated':>10} {'kept':>6} {'best':>8} {'seconds':>8}")
for record in result.trace:
print(
f"{record.depth:>6} {record.n_proposed:>10} {record.n_evaluated:>10} "
f"{record.n_kept:>6} {record.best_score:>8.4f} {record.elapsed:>8.3f}"
)
depth proposed evaluated kept best seconds
0 4 4 4 0.5941 0.001
1 56 52 30 0.8411 0.016
2 3038 2994 30 0.9999 0.611
How beam width controls cost¶
Beam width is the parameter that bounds work. Widening it explores more of the space at proportionally greater cost; narrowing it risks pruning a parent whose children would have been valuable.
print(f"{'beam':>6} {'proposed':>10} {'kept':>6} {'seconds':>9} top feature")
for width in (5, 10, 25, 50, 100):
outcome = BeamSearch(max_depth=2, beam_width=width, random_state=0).run(data, target)
top = outcome.names[0] if outcome.names else "-"
print(
f"{width:>6} {outcome.n_proposed_total:>10} {len(outcome):>6} "
f"{outcome.elapsed:>9.3f} {top}"
)
beam proposed kept seconds top feature
5 219 5 0.019 ((a / c) * b)
10 509 10 0.040 ((a / c) * b)
25 2265 25 0.214 ((a / c) * b)
50 3490 50 0.299 ((a / c) * b)
100 3490 100 0.302 ((a / c) * b)
Note that cost grows with beam width but the quality of the top feature plateaus quickly. In practice a moderate beam is usually enough, and the remaining budget is better spent on depth.
How depth controls expressiveness¶
Some relationships simply cannot be expressed below a certain depth. (a * b) / c needs depth 2: one step to build a * b, another to divide by c.
for depth in (0, 1, 2, 3):
outcome = BeamSearch(max_depth=depth, beam_width=25, random_state=0).run(data, target)
found = any(all(token in name for token in "abc") for name in outcome.names)
print(f"depth {depth}: {len(outcome):>3} features, "
f"three-way interaction found: {found}")
depth 0: 4 features, three-way interaction found: False depth 1: 25 features, three-way interaction found: False depth 2: 25 features, three-way interaction found: True
depth 3: 25 features, three-way interaction found: True
Choosing a scorer¶
Three scoring strategies are available, and they differ in what they can detect and what they cost.
"correlation"— absolute Pearson correlation against the residual. One matrix product per batch. The default."mutual_information"— nearest-neighbour mutual information. Detects dependence that correlation misses."gradient_boosting"— measured out-of-fold predictive improvement. Scores what is actually being optimised.
Where they agree¶
On a relationship whose nonlinearity is captured by the expression itself, all three find the same thing, and correlation is much cheaper.
import time
for name in ("correlation", "mutual_information"):
started = time.perf_counter()
outcome = BeamSearch(scorer=name, max_depth=2, beam_width=20, random_state=0).run(data, target)
elapsed = time.perf_counter() - started
print(f"{name:>20}: {elapsed:>6.2f}s top: {outcome.names[0]}")
correlation: 0.11s top: ((a / c) * b)
mutual_information: 7.73s top: ((a / c) * b)
Where they differ¶
Correlation is blind to symmetric relationships. If the target depends on the magnitude of a feature but not its sign, correlation sees nothing while mutual information sees the dependence clearly.
from beamfeat.scoring import CorrelationScorer, MutualInformationScorer
symmetric_feature = rng.uniform(-3, 3, n)
symmetric_target = symmetric_feature**2 + rng.normal(0, 0.3, n)
noise = rng.normal(size=n)
candidates = np.column_stack([symmetric_feature, noise])
corr = CorrelationScorer().score_batch(candidates, symmetric_target)
mutual = MutualInformationScorer().score_batch(candidates, symmetric_target)
print(f"{'':>12} {'true feature':>14} {'noise':>10}")
print(f"{'correlation':>12} {corr[0]:>14.4f} {corr[1]:>10.4f}")
print(f"{'mutual info':>12} {mutual[0]:>14.4f} {mutual[1]:>10.4f}")
true feature noise correlation 0.0149 0.0493 mutual info 1.7938 0.0000
Correlation cannot distinguish the true feature from noise here. Mutual information can.
In practice this matters less than it appears, because beamfeat applies the
scorer to transformed columns. Once square(x) has been constructed, a
linear scorer detects it easily — which is why correlation remains a sensible
default despite this weakness.
squared = symmetric_feature**2 # what the search would hand the scorer
print(f"correlation on the raw column: {CorrelationScorer().score(symmetric_feature, symmetric_target):.4f}")
print(f"correlation on the squared column: {CorrelationScorer().score(squared, symmetric_target):.4f}")
correlation on the raw column: 0.0149 correlation on the squared column: 0.9926
Redundancy control¶
A beam scored purely on individual merit fills with variants of the same signal. beamfeat scores candidates against the residual left by what is already selected, and additionally drops candidates too correlated with the current beam.
from beamfeat import Evaluator
for threshold in (0.5, 0.9, 0.999):
outcome = BeamSearch(
max_depth=2, beam_width=30, redundancy_threshold=threshold, random_state=0
).run(data, target)
evaluator = Evaluator(data)
_, matrix = evaluator.evaluate_many(outcome.nodes)
correlations = np.abs(np.corrcoef(matrix, rowvar=False))
np.fill_diagonal(correlations, 0.0)
print(
f"threshold {threshold:<6}: {len(outcome):>3} features, "
f"max pairwise correlation {np.nanmax(correlations):.3f}"
)
threshold 0.5 : 14 features, max pairwise correlation 0.497
threshold 0.9 : 30 features, max pairwise correlation 0.899
threshold 0.999 : 30 features, max pairwise correlation 0.997
Auditing what was rejected¶
Candidates excluded during evaluation are recorded rather than silently dropped. This matters when a search returns less than expected: the log tells you whether candidates were numerically invalid, degenerate, or simply never proposed.
signed_data = dict(data)
signed_data["e"] = rng.normal(0, 2, n) # negative values break log and sqrt
outcome = BeamSearch(max_depth=1, beam_width=20, random_state=0).run(signed_data, target)
print("rejections by reason:")
for reason, count in outcome.evaluation_log.counts().items():
print(f" {reason.value:>15}: {count}")
print("\nexamples:")
for record in list(outcome.evaluation_log)[:5]:
print(f" {record}")
rejections by reason:
domain_error: 6
examples:
abs(a) [domain_error]: input outside domain of 'abs'
abs(b) [domain_error]: input outside domain of 'abs'
abs(c) [domain_error]: input outside domain of 'abs'
abs(d) [domain_error]: input outside domain of 'abs'
log(e) [domain_error]: input outside domain of 'log'