Feature Selection Methods

You inherit a click-through-rate model with 300 candidate features, a laptop-class machine, and a deadline. Training a wrapper on all 300 is hours per fold; you need a defensible first cut before any expensive modeling. Filter methods do exactly that: they score features using only the data matrix (and optionally the target), with no model in the loop, which makes them fast, model-independent, and reusable across estimators — at the price of being blind to feature interactions.

The cheapest filter is variance itself. VarianceThreshold(0.01) drops features whose variance falls below a floor — a column that is 0 for 99.97% of rows (a "pixel" that never fired, a rare error code) carries no information for any split, distance, or coefficient. Scikit-learn's default drops exactly-zero-variance columns; with standardized features, set thresholds on the raw frame before scaling.

Next, score each feature's association with the target. Pearson correlation is the reflex choice:

import numpy as np
import pandas as pd

corr = df.corr(numeric_only=True)["clicked"].abs().sort_values(ascending=False)
top_linear = corr.head(20)

# redundancy pass: drop one of any pair correlated above 0.95
feat_corr = df.drop(columns="clicked").corr(numeric_only=True).abs()
upper = feat_corr.where(np.triu(np.ones(feat_corr.shape), k=1).astype(bool))
redundant = [c for c in upper.columns if (upper[c] > 0.95).any()]

Pearson has one blindness that matters: it only detects linear association. A feature where risk rises and then falls with the input (a U-shape: young and old drivers crash more) has correlation ~0 with the target and a filter based on it deletes a genuinely predictive feature. That gap is what mutual information fills:

I(X; Y) = SUM over x,y of  p(x,y) * log2[ p(x,y) / (p(x) p(y)) ]
        = KL divergence between the joint distribution and the product of marginals

Mutual information is zero if and only if X and Y are independent, it captures any functional relationship (U-shapes, thresholds, XOR), and it is symmetric — but unlike correlation it is not signed, so you cannot tell which way a feature pushes the target from the score alone. In scikit-learn:

from sklearn.feature_selection import SelectKBest, mutual_info_classif

sel = SelectKBest(mutual_info_classif, k=25)
X25 = sel.fit_transform(X_train, y_train)     # estimate MI on train, always
scores = sel.scores_                           # one MI estimate per feature

mutual_info_classif (and mutual_info_regression) estimates the entropies from k-nearest-neighbor distances, so continuous features work directly; documented gotcha: discrete-valued columns get a small random jitter added, so pass random_state for reproducibility. A second gotcha bites everyone once: finite-sample MI of pure noise is strictly positive. A random column scores maybe 0.003 while a weak-but-real feature scores 0.010; rank relative to a permutation baseline (shuffle y, re-run, use the noise ceiling) instead of believing the absolute numbers.

The filter family, compared:

And the structural failure you must be able to state in an interview: filters score features individually, so they cannot see interactions. The canonical counterexample is XOR: let y = x1 XOR x2 with independent Bernoulli x1, x2. Both features have exactly zero correlation and zero mutual information with y, yet together they predict it perfectly. A SelectKBest run on XOR data drops both features and destroys the signal. Nothing in the filter's math can rescue it, because the information lives entirely in the joint distribution. When interaction terms plausibly matter — and in most tabular problems they do — filters are a pre-filter, and the final say belongs to model-aware selection: wrappers, which the next chapter covers, and embedded methods like L1, which follow it. If your dimensionality problem is really "too many redundant measurements" rather than "too many features", the compression route from Principal Component Analysis remains the lighter tool.

Filters ask "how informative is each feature alone?". Wrapper methods ask the only question that ultimately matters: how well does a model perform on this subset of features? They treat feature selection as a search problem — evaluate candidate subsets by cross-validated model score, keep what wins. Model-awareness is the superpower: interactions like XOR are visible because the actual estimator sees the actual subset. The price is compute, because every candidate subset means retraining.

Exhaustive search over 2^p subsets dies immediately: p = 20 gives ~10^6 subsets, p = 40 gives ~10^12. Every practical wrapper is therefore greedy. Forward selection starts empty and adds the feature that most improves CV score; backward elimination starts full and removes the least useful; and the workhorse, recursive feature elimination (RFE), fits the model on the current set, ranks features by the model's own importance (coefficients or feature_importances_), removes the weakest, and repeats:

from sklearn.feature_selection import RFECV
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold

est = LogisticRegression(max_iter=2000, solver="liblinear")
rfecv = RFECV(estimator=est, step=1,
              cv=StratifiedKFold(5), scoring="accuracy")
rfecv.fit(X_train, y_train)

print(rfecv.n_features_)        # e.g. 37
print(rfecv.support_)           # boolean mask of kept features
print(rfecv.ranking_)           # 1 = kept; rank = elimination round + 1
X_train_sel = X_train[:, rfecv.support_]

RFECV wraps RFE with cross-validation and a feature-count sweep (it evaluates every k from 1 to p by default, using the ranking order), so it answers both "which features" and "how many" in one run. Swap in a tree ensemble and RFE switches to impurity importances automatically — the estimator's importance attribute is the entire interface.

Now the cost arithmetic, because "wrappers are expensive" is only useful when you can quantify it. RFE with step=1 and cv=5 fits 5 * (p + (p-1) + ... + 1) = 5p(p+1)/2 models: p = 50 gives 6,375 fits, p = 200 gives over 100,000. Three standard mitigations:

That last row is the pattern senior engineers use in practice: let a fast embedded method (next chapter) or a filter do the coarse ranking, then run RFE only on the surviving shortlist.

Two traps complete the picture. First, selection bias: you pick the subset with the best CV score, so that score is optimistically biased — you selected noise peaks among 50 candidates. The honest estimate comes from nested cross-validation (selection inside the outer folds) or, minimally, from a held-out test set that influenced no part of the selection. Second, instability: rerun RFE on bootstrap resamples and correlated features will trade ranks freely — two features carrying the same signal look interchangeable, so which one "wins" is a coin flip. Measure it (Jaccard overlap of selected sets across resamples), and when sets are unstable, break ties on engineering cost — latency, availability of the feature at serving time — rather than on a CV difference of 0.002.

import numpy as np
from sklearn.utils import resample
from sklearn.feature_selection import RFE

sets = []
for _ in range(10):
    Xb, yb = resample(X_train, y_train, random_state=None)
    est.fit(Xb, yb)
    rfe = RFE(est, n_features_to_select=37).fit(Xb, yb)
    sets.append(set(np.array(X_train.columns)[rfe.support_]))

overlap = np.mean([[len(a & b) / len(a | b) for b in sets] for a in sets])
print("mean pairwise Jaccard:", overlap)   # <0.5 means the subset choice is fragile

Wrappers give you the most accurate subset-per-model answer available; they are the right tool when p is already moderate (tens, not thousands) and the model is cheap enough to retrain dozens of times. When p is in the thousands, you now know the full escalation path: filters cut to tens, embedded methods (L1, tree importances) cut to the final set, and wrappers — if budget remains — verify. The roadmap next makes the case that the embedded tier is often where the search should end.

Embedded methods let the estimator select its own features as a side effect of training: model-aware like wrappers, single-fit cost like filters. In interviews, two embedded mechanisms dominate — L1 regularization, which produces exact zeros, and tree-based importances, which rank without thresholding — and you are expected to explain why each works, not just that it does.

L1 and the geometry of zeros. L1-regularized linear models optimize loss + alpha * sum(|w_j|). The question interviewers pose: why does L1 set coefficients to exactly zero when L2 never does? The answer is geometric. Both penalties restrict w to a region — L1's region is the diamond |w_1| + |w_2| <= t (a cross-polytope in higher dimensions), L2's is the circle w_1^2 + w_2^2 <= t. The optimum sits where a contour of the loss first touches that region. The L1 diamond has corners on the axes, and a contour touching a corner sets the perpendicular coordinate to exactly zero; the smooth L2 circle almost never touches an axis, so coefficients shrink toward zero but land near it, never on it. That one picture is the whole explanation.

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np

# Standardize first: the penalty applies to coefficient SCALE, which depends on feature units
l1 = Pipeline([
    ("scaler", StandardScaler()),
    ("clf", LogisticRegression(penalty="l1", solver="liblinear",
                               C=0.1, max_iter=2000)),
]).fit(X_train, y_train)

mask = l1.named_steps["clf"].coef_[0] != 0
print(f"{mask.sum()} of {mask.size()} features kept")
X_train_sel = X_train.loc[:, mask]

Sweep C (the inverse of alpha) and you trace the regularization path: features enter one at a time as C grows, roughly ordered by individual signal. Two standard caveats: correlated features are interchangeable under L1, which picks one arbitrarily (unstable selections, same resample problem as wrappers — elastic net, with its l1_ratio mixing, groups correlated features instead); and standardization is mandatory, because a feature measured in units 100x smaller gets a coefficient 100x larger and takes a 100x bigger L1 hit. SelectFromModel(l1_clf) wraps the mask into a pipeline step so selection and fitting stay inside cross-validation folds.

Trees rank without a penalty. Gradient boosting and random forests expose feature_importances_: the mean decrease in impurity (Gini or entropy) each feature contributed across all splits, weighted by sample counts. Convenient — one fit, full ranking — but two biases you must know. First, MDI importances are biased toward high-cardinality continuous features: a feature with 1,000 distinct values offers 1,000 candidate split points and more chances to find a spurious impurity drop, while a binary flag gets one; scikit-learn's own docs warn about this. Second, correlated features split credit: each looks individually weak because the other covers its splits.

The remedy is permutation importance, computed on data the model did not train on:

from sklearn.inspection import permutation_importance

r = permutation_importance(lgbm, X_val, y_val, n_repeats=10, random_state=0)
order = r.importances_mean.argsort()[::-1]
for i in order[:10]:
    print(X_val.columns[i], r.importances_mean[i].round(4))

Each feature is shuffled on the validation set and the score drop is measured: a feature's importance is how much the model's predictions suffer without it. Unbiased by cardinality, computed against real performance. Remaining caveat: with correlated features, shuffling one is masked by its twin — both look dispensable. Cluster correlated features first (or use grouped permutation) before trusting the ranking.

The full taxonomy, as an interview summary table:

The decision heuristic worth reciting: start with filters if p is in the thousands, finish with embedded methods, add a wrapper only when p is in the tens and the accuracy stakes justify thousands of fits — and always keep selection inside the CV pipeline, the same leakage rule that governs Principal Component Analysis. Selection keeps your original, interpretable features; the next topic trades interpretability for something selection can never give you — a two-dimensional map of the data, from t-SNE and UMAP Visualization.