01 — Feature Engineering

A machine-learning model cannot read a database row. It reads a fixed-length list of numbers. Feature engineering is the work of turning the row into that list.

Some of those transformations are required: without them, a particular model’s arithmetic is provably broken. Others are cosmetic: the model returns the identical answer with or without them. Telling the two apart is most of the skill.

Every rule below is derived from a line of the model’s own math. Each technique comes with the assumption it rests on. Three questions decide any transform: which models need it, what has to be true of your data for it to work, and how it fails silently when that assumption does not hold.

The four words used on every page

Four terms recur throughout. In plain English:

• A feature is one input column the model reads — age, income, zip_code.
• A label (also called the target, and written y) is the answer you want predicted: did this customer cancel, what did this house sell for.
• Training is the offline process that fits the model to historical rows whose label is already known.
• Serving (or inference) is the online process that runs the fitted model on a new row whose label is not yet known — and may not be known for weeks.

What goes in, and what comes out

The whole job has one concrete signature: one raw row goes in, one numeric vector comes out.

The block below is a before-and-after. The top half is a row as your warehouse stores it: mixed types, a date string, a null. The bottom half is what the model receives: seven numbers. The comments name the transform that produced each number, and each one is a decision this chapter works through.

INPUT   one raw row, as your warehouse stores it
        {user_id: 8814, age: 34, income: 82000, zip: "94107",
         last_login_ts: "2024-02-27", plan: "pro", refunds_30d: null}

OUTPUT  one fixed-length vector of numbers, as the model consumes it
        [ -0.42,   # age, standardized
           0.31,   # log(income), then standardized
           0.061,  # zip, TARGET-ENCODED (the category replaced by the average
                   #   label among its rows) on OUT-OF-FOLD data (computed only
                   #   from rows in other cross-validation folds, so this row's
                   #   own label cannot reach its own feature)
           3.0,    # days since last login, as of the prediction date
           1, 0,   # plan = "pro", ONE-HOT (one 0/1 column per level) with
                   #   "basic" as the dropped baseline
           1 ]     # refunds_30d was missing -- a missingness indicator

Four questions follow from that output:

• Why age is standardized, and why a tree would not care (Scaling which models care and exactly why).
• Why income gets a log before anything else (Log box cox and when a transform is doing real work).
• Why the zip encoding is computed on rows that exclude this one, and what happens to your production metrics if it is not (Categorical encoding).
• Why the null becomes a 1 in a new column rather than being replaced by an average (Missing values three mechanisms three different correct answers).

The property that matters more than any single transform

The same function must produce the same vector during training and during serving.

That sentence explains the chapter’s shape. Sections The pipeline and where the bugs live through Text and image features briefly build the vector. Feature selection cuts it back down — it is the Select stage of the pipeline in §1.

Feature stores and trainingserving skew is not about building the vector. It is about the “same function” clause above, which is where production models most often fail.

The organizing question is never “should I scale this?” It is “what in this model’s math touches the scale?” That question, answered once per model, replaces the checklist.

Chapter 02 covers the models themselves: how each one fits and what else it assumes.

1. The pipeline, and where the bugs live

Every transform in this chapter occupies one box on the same pipeline. The diagram runs left to right: six processing stages, then a store, then a fork into two consumers. Two facts about that shape cause most of the bugs ahead — the boxes are fit in a fixed order, and the last arrow splits in two — and both are covered below, along with leakage, the failure that recurs in every later section.

flowchart LR
    R[("Raw source<br/>events · tables · logs")] --> T["Transform<br/>scale · log · bin"]
    T --> E["Encode<br/>categoricals"]
    E --> I["Impute<br/>missing"]
    I --> X["Cross<br/>interactions"]
    X --> S["Select"]
    S --> FS[("Feature store")]
    FS --> TR["Training<br/>offline · batch"]
    FS --> SV["Serving<br/>online · per request"]

    style R fill:#1d3557,color:#fff
    style FS fill:#7209b7,color:#fff
    style TR fill:#2d6a4f,color:#fff
    style SV fill:#bc6c25,color:#fff

Box by box:

• Raw source — whatever your systems already record: events, tables, logs.
• Transform — the first decision on a numeric column: rescale it, reshape it, or bucket it (Scaling which models care and exactly why, Log box cox and when a transform is doing real work, Binning mostly a trap occasionally the point).
• Encode — turn text-valued columns such as zip into numbers. A model multiplies its inputs, and you cannot multiply the string "94107" (Categorical encoding).
• Impute — decide what to put where a value was never recorded. Imputation is just the act of filling in a missing value with a guess (Missing values three mechanisms three different correct answers).
• Cross — build combinations of features, so the model can react to two columns jointly rather than one at a time (Crosses and interactions).
• Select — drop the columns that do not earn their keep (Feature selection).

What survives lands in a feature store: a database of already-computed feature values.

The store then forks into two consumers. Training is offline — one batch job reads the whole history at once. Serving is online — one lookup answers one live request in milliseconds. Same feature name, two very different execution environments.

Two structural facts that generate most of this chapter

Fact one: every box is fit on data.

To fit a transform is to compute and store numbers from a dataset before the transform can run. A scaler stores a mean. An imputer stores a median. A target encoder stores one statistic per category. Those stored numbers are as much a part of the model as its weights, and they must ship with it.

That is where the recurring hazard comes from. Leakage is when a feature carries information that would not have been available at the moment of prediction. Most often it is the label itself, smuggled in through a statistic computed over rows you are about to score.

The rule is short: if any box is fit on rows that will later be scored, you have leaked.

The symptom is always the same. The model scores beautifully offline and disappoints in production, because the information it leaned on does not exist when a real request arrives (Categorical encoding, Missing values three mechanisms three different correct answers, Feature selection).

Fact two: the pipeline forks at the end.

Training and serving are two different executions of what is supposed to be the same function. When they diverge, nothing raises an exception and no monitor fires (Feature stores and trainingserving skew).

2. Scaling: which models care, and exactly why

Whether a model needs its inputs rescaled is not a matter of convention. It falls straight out of the model’s own math — some families are genuinely broken by unscaled features, others are completely indifferent — and standardization itself makes assumptions about your data that are worth stating out loud.

Scaling means putting every feature on a comparable numeric range.

The default form is standardization: subtract the column’s mean, then divide by its standard deviation (sd, the typical distance of a value from the mean). The result is a column centered at 0 with a spread of 1.

Suppose the training set’s mean age is 40.3 with an sd of 15. Then the 34-year-old from the row at the top of this chapter becomes:

(34 - 40.3) / 15  =  -6.3 / 15  =  -0.42

That is the first entry of the output vector. “This person is 0.42 standard deviations younger than average” — the units are gone, and every other standardized column now speaks the same language.

The word normalization is used loosely in practice, but it usually means min-max scaling, which squeezes a column into [0, 1].

Three model families give three different answers about whether any of this is needed, and each answer follows from one line of that model’s own math.

Gradient-based models: required, and here is the arithmetic

Most models are fit by gradient descent: guess the weights, compute the direction in which the error falls fastest, take a small step in that direction, repeat until the error stops improving.

Least squares is the fitting rule that picks the weights minimizing the sum of squared prediction errors. Picture its error surface as a bowl. The bottom of the bowl is the best set of weights, and gradient descent is a ball rolling toward it.

How fast the ball reaches the bottom depends entirely on the bowl’s shape. Four terms describe that shape, and each one is needed for the argument:

• The Hessian H is the matrix of second derivatives of the error. Read it as the curvature: how sharply the error rises as you move away from the best answer, in each direction.
• An eigenvalue of H is the steepness along one special direction of that bowl. Write L for the largest (the steepest direction, the narrow walls of the bowl) and mu for the smallest (the flattest direction, the long shallow trough).
• The condition number kappa = L/mu is the ratio of the two. In words: how much steeper the steepest direction is than the flattest. kappa = 1 is a perfectly round bowl. A large kappa is a long, narrow canyon.
• w_t is the weight vector after t steps, and w* is the best weight vector — the bottom of the bowl.

Why the eigenvalues are just the variances

For least squares on centered features (each column’s mean already subtracted), the Hessian is H = (2/n) X^T X, where X is the matrix of feature values and n is the row count.

Take two features that are uncorrelated with each other: age (sd 15) and income (sd 45,000). Uncorrelated and centered means every off-diagonal entry of X^T X / n is zero, so that matrix is diag(var_age, var_income). And the eigenvalues of a diagonal matrix are exactly its diagonal entries.

So H = 2 · diag(var_age, var_income), and its eigenvalues are the two variances scaled by the same constant 2. That constant divides out of the ratio L/mu.

So kappa here is simply the ratio of the two column variances. Which means it is set entirely by your choice of units.

Substituting the numbers

The first line below is the standard convergence bound for gradient descent; the rest is arithmetic on it.

||w_t - w*||  <=  ((kappa-1)/(kappa+1))^t · ||w_0 - w*||

var(age) = 15^2 = 225 ;  var(income) = 45,000^2 = 2,025,000,000
kappa = 2,025,000,000 / 225 = 9,000,000
per-step factor = (kappa-1)/(kappa+1) ≈ 1 - 2/kappa = 1 - 2.22e-7 = 0.999999778
steps to cut the error by 1e-6 = ln(1e6)/2.22e-7 = 13.8/2.22e-7 ≈ 62,000,000

The bound reads: the distance from the current weights to the best weights shrinks by the same factor on every step, and that factor is set entirely by kappa.

Follow it through. With raw dollars against raw years the per-step factor is 0.999999778 — each step removes about two ten-millionths of the remaining error. Cutting that error by a millionfold therefore takes roughly 62 million iterations.

Now standardize both columns. Both variances become 1, so kappa = 1, so the per-step factor is (1-1)/(1+1) = 0, and gradient descent lands on the answer in one step.

That is the whole argument: 6.2e7 iterations versus 1. Scaling is not hygiene for gradient-based models. It is the difference between a solvable and an unsolvable optimization at any realistic iteration budget.

The same cause, seen as a learning-rate problem

The learning rate eta is the size of the step taken each iteration, and you get one value of it for all features.

That single value is squeezed from both sides. It must satisfy eta < 2/L or the ball overshoots the narrow walls and the run diverges. But the same eta then moves along the flat direction kappa times more slowly. With kappa = 9,000,000 there is no choice that serves both.

Per-coordinate optimizers paper over this. Adam maintains a separate effective step size for every feature, which is why deep nets trained with Adam are less scale-sensitive than ones trained with plain SGD (stochastic gradient descent — the same loop, with the direction estimated from a small random batch of rows rather than all of them). Less sensitive, not immune.

What standardization itself assumes

Standardization assumes two things. First, that the mean and standard deviation are meaningful summaries of the column, which needs a roughly symmetric distribution with finite variance and no dominant outliers. Second, that the statistics computed during training still describe the data arriving at serving time.

It breaks in two matching ways.

One heavy tail or data-entry error drags the mean and inflates the sd, which compresses every ordinary value into a narrow band near zero. That is the case for the robust alternative described at the end of this section.

And under distribution shift — where the pattern of live data drifts away from the training data — a mean computed six months ago no longer centers the column, so a value that should standardize to 0 arrives as 2.

Distance-based models: required, and the failure is silent

A second family of models measures similarity between rows, which makes it sensitive to units in a more direct way. These are the members:

• kNN (k-nearest neighbours) predicts a row’s label by finding the k most similar training rows and averaging their labels.
• k-means groups rows into k clusters by repeatedly assigning each row to the nearest cluster center.
• An SVM (support vector machine) with an RBF (radial basis function) kernel draws a boundary whose shape is a function of the distance between points.
• PCA (principal component analysis) finds the directions along which the data varies most — and variance is measured in the features’ own units.

Anything built on cosine similarity (the angle between two rows treated as arrows from the origin) or Euclidean distance (straight-line distance: the square root of the summed squared per-feature differences) inherits the same problem.

Here is the problem, on two customers who differ by 30 years of age and $500 of income. Note which gap dominates the distance before and after scaling.

A: age 25, income $50,000        B: age 55, income $50,500

raw:          d^2 = 30^2 + 500^2   = 900 + 250,000    age is 0.36% of the distance
standardized: d^2 = 2.00^2 + 0.011^2 = 4.000 + 0.00012  income is 0.003%

Take the raw line first. The 30-year age gap contributes 30^2 = 900 to the squared distance. The $500 income gap contributes 500^2 = 250,000. The income term is 250,000 / 900 = 278 times bigger — so a $500 income difference, which is noise at this income level, outweighs a 30-year age difference by 278x. Age accounts for 900 / 250,900 = 0.36% of the total.

Now the standardized line, using the same sds as before (age 15, income 45,000). The age gap becomes 30/15 = 2.00 standard deviations and the income gap becomes 500/45,000 = 0.011 standard deviations. Squared, that is 4.000 against 0.00012, and the ordering has completely reversed: income is now 0.003% of the distance.

On raw features, kNN is not a 2-feature model. It is a 1-feature model on whichever column happens to have the largest units.

Nothing errors. Accuracy is simply lower than it should be, and no diagnostic points at the cause.

Assumptions. Distance-based methods assume every feature is on a common scale and that each one deserves equal weight in the similarity judgement.

Scaling delivers the first. The second is an assertion you are making about the problem, and it fails when several features are irrelevant, because each irrelevant column adds its full share of noise to every distance. That is also why these models degrade as the number of features grows — the curse of dimensionality, worked through in Knn and the curse of dimensionality.

Trees: cosmetic, and here is the proof

A decision tree predicts by asking a nested series of yes/no questions about single features, then reading the answer off the leaf a row lands in.

Each question is a split of the form x <= t: pick a feature x and a threshold t, and send rows left or right. The tree chooses t using a criterion — gini, entropy or variance — that scores how mixed the labels are inside the two resulting groups, and keeps the split that leaves them least mixed.

That form is what makes trees immune to scaling. For any strictly increasing function f — one that always goes up, so it preserves order — x <= t holds if and only if f(x) <= f(t).

Concretely, with four ages and a split at 35:

raw:      25   30  | 35 |  40    60          left = {25, 30}
log:    3.219 3.401 | 3.555 | 3.689 4.094    left = {25, 30}   -- same rows

The threshold moved from 35 to ln(35) = 3.555, and not one row changed sides.

That generalizes. The set of achievable row partitions is identical before and after any such monotone transform, and gini, entropy and variance are all computed on row partitions. So the fitted tree is bit-for-bit the same; only the printed threshold changes.

This covers standardization, min-max, log, and rank transforms applied to a feature. It also carries over to the ensembles built from trees: random forests (many trees fit to random subsets of the data, then averaged) and GBDT (gradient-boosted decision trees, where each new tree is fit to the errors of the ones before it).

The exception people miss: a monotone transform of the target is not cosmetic for a regression tree. The split criterion is variance reduction on y, and Var(log y) != Var(y) under any reweighting of rows. Log-transforming y genuinely changes which splits win.

Assumptions. The invariance holds only when the transform is strictly increasing and applied to a feature. Apply something non-monotone — a residual, an absolute deviation from a mean, a hand-built bucket that reorders values — and the achievable partitions really do change. Apply it to the target and the criterion itself changes, as the note above says.

The verdict for every common model

The table collects the three arguments above and applies the same reasoning to the rest of the common models. The column that matters is the third one: in every row, the verdict falls out of a specific piece of that model’s math, never out of a convention.

Four rows in that table need a gloss.

• OLS (ordinary least squares) is least squares solved by a formula in one shot rather than by iteration. There is no convergence rate to ruin, so nothing to fix.
• Ridge, Lasso and ElasticNet are linear models with regularization — an extra penalty term added to the loss that charges the model for large coefficients, to stop it from fitting noise. The penalty is charged on the coefficient exactly as written, so a feature measured in dollars gets a different effective penalty from the same feature measured in thousands of dollars. Same information, different model.
• Saturating activations are the squashing functions inside a neural network, such as the sigmoid, which flatten out for large inputs. Feed one a raw income and it returns a value pinned at 1 with a gradient of essentially zero, so that unit stops learning entirely.
• Naive Bayes estimates one distribution per feature independently and never compares features to each other, so there is no shared metric for scale to distort.

Choosing a scaling method

Once you have decided to scale, three methods are in common use, and they differ in what they assume about your data rather than in what they compute. Read the last column first — it is where each one fails.

The last column is the assumption, stated as its violation.

Standardize assumes no single value dominates the mean and sd.

Min-max assumes the training minimum and maximum bound everything you will ever see. That is a promise about the future, and live data is free to break it.

Robust replaces both statistics with order-based ones. The median is the middle value when the data is sorted, and the IQR (interquartile range) is the distance between the 25th and 75th percentiles — the width of the middle half of the data. Neither one moves when a single point flies off to infinity, which is why its failure column is nearly empty.

3. Log, Box-Cox, and when a transform is doing real work

Scaling leaves a column’s shape alone; the log transform exists to change it. But “fix the skew” undersells what a log actually does — it performs three distinct jobs, each with its own arithmetic and its own assumption, and in one common setting it does no work at all.

Skew is the asymmetry of a distribution — a long tail running out on one side, which is what income, revenue and session length all look like. Only two of the three jobs below are about skew, which is why “my data is skewed so I logged it” is not a complete answer in an interview.

Job 1 — turning a multiplicative truth into an additive model

A linear model asserts y = w1·x1 + w2·x2. It adds contributions. That is all it can do.

Now suppose the world is actually multiplicative: y = a · x1^b · x2^c. No choice of w1 and w2 fits that, because the truth is a product and the model is a sum.

Take logs of both sides:

log y = log a + b·log x1 + c·log x2

That is exactly the form the model can represent — a constant plus weighted inputs. Feed the model log x1 and log x2, predict log y, and the weights it learns are b and c.

This is not cosmetic. It is a model-capacity fix: before the transform the relationship was unrepresentable, and after it, it is not.

Assumption: the underlying relationship really is a product of powers. That is common for quantities that compound — revenue, area, dose response — and wrong for quantities that add.

Job 2 — variance stabilization

Sometimes the noise around the prediction grows with the prediction itself. A model predicting $10 orders is off by cents; the same model predicting $10,000 orders is off by hundreds. The percentage error is roughly constant while the absolute error is not.

Written formally, that condition is sd(y | x) = c·E[y | x]: the spread is a constant fraction c of the level mu.

The delta method is the rule for pushing a noisy quantity through a smooth function g: the output variance is approximately g'(mu)^2 times the input variance, where g' is the derivative. Apply it with g = log, whose derivative is 1/y:

Var(log y) ≈ (d/dy log y)^2 · Var(y) = (1/mu)^2 · (c·mu)^2 = c^2

Follow the cancellation. Var(y) is (c·mu)^2 because the sd is c·mu. The derivative term contributes (1/mu)^2. The two mu^2 factors cancel, leaving c^2 — a constant, independent of the level.

Heteroscedasticity is the name for the original condition: noise that varies with the level, which breaks the equal-noise assumption behind least squares’ standard errors. After the log it is gone by construction, which restores the OLS efficiency and standard-error assumptions (ch 02).

Assumption: the noise is proportional to the level. If it is instead constant, the log creates heteroscedasticity where none existed.

Job 3 — killing leverage

Leverage measures how much a single row can pull the fitted line toward itself. It runs from 1/n (this row is no more influential than any other) to 1 (the line is obliged to pass exactly through this point).

The OLS leverage of row i is:

h_i = 1/n + (x_i - xbar)^2 / sum_j (x_j - xbar)^2

Read that as a share. The numerator is how far row i sits from the mean, squared. The denominator is the same quantity summed over every row. So a point far from the mean takes a large share of the total and therefore a large share of the influence.

Now a concrete case: 1,000 incomes, lognormally distributed around $60k, with one data-entry error at $10,000,000. Compare the two scales.

raw dollars:   999 normal points ≈ 999 · 28,500^2 = 8.1e11 ; outlier (1e7-7e4)^2 = 9.86e13
               h_outlier = 0.001 + 9.86e13/9.94e13 = 0.992

natural log:   sd(ln income) ≈ 0.45 ; ln(1e7) - 11.0 = 5.12
               999 · 0.2025 = 202.3 ; outlier 26.2 ; total 228.5
               h_outlier = 0.001 + 26.2/228.5 = 0.116

In dollars, the 999 ordinary points sit about $28,500 from the mean each, contributing 999 · 28,500^2 = 8.1e11 between them. The outlier sits 1e7 - 7e4 = 9,930,000 from the mean, contributing 9.86e13 on its own — 122 times the other 999 combined. Its share of the denominator is 9.86e13 / 9.94e13 = 0.991, and adding the baseline 1/n = 0.001 gives h = 0.992.

On the log scale the same point is 5.12 units from the mean against a spread of 0.45, so it contributes 5.12^2 = 26.2 against the other points’ 999 · 0.45^2 = 202.3. Its share is 26.2 / 228.5 = 0.115, giving h = 0.116.

A leverage of 0.99 means the fitted line passes through that one point and treats the other 999 as noise. The log takes it to 0.12. That is the transform earning its place, measurably.

When log is pure ceremony

Feeding a tree ensemble. Skew in a feature is invisible to a monotone-invariant splitter, for exactly the reason worked through in Scaling which models care and exactly why: the log moves the threshold and nothing else.

Three relatives, and what each one assumes

Box-Cox generalizes the log into a family: (y^lambda - 1)/lambda, where lambda = 0 recovers log y. It picks the exponent lambda by maximum likelihood — choosing the value under which the observed data is most probable — targeting approximate normality. It requires strictly positive input.

Yeo-Johnson is the variant of Box-Cox that tolerates zeros and negative values.

Both fit lambda on training data, so lambda is a parameter you must persist and ship with the model, exactly like a scaler’s stored mean.

The quantile (or rank) transform is the most aggressive of the three. A quantile is a value’s position in the sorted order, expressed as a fraction — the 300th of 1,000 sorted values sits at quantile 0.30. Replace each value by its quantile and you force a uniform or normal marginal distribution, destroying all spacing information and keeping only order.

That makes its assumption the strongest: that spacing carries no signal. Right for meaningless units such as composite scores. Wrong when the spacing is the signal, as with dollars and durations.

4. Binning: mostly a trap, occasionally the point

After rescaling and reshaping, the third thing you can do to a numeric column is chop it into pieces — and this one, unlike the other two, throws information away. It buys exactly one capability, and the price is paid on every row.

Binning (or discretization) replaces a numeric column with the label of the range it falls in. age becomes one of [18-25, 26-35, 36-50, 51+].

What binning buys

Exactly one thing: a linear model can now fit a non-monotone relationship — one that goes up and then back down.

Say default risk is high for the young, low in the middle, and high again for the old. A linear model has one coefficient for age, and one coefficient can only say “risk rises with age” or “risk falls with age.” It cannot say both.

Bin the column and the model gets three dummies — 0/1 indicator columns, one per bin beyond the first — plus the intercept, which is the model’s baseline value when every input is zero. Four free numbers instead of one, and each bin’s risk is now estimated independently of the others, so up-then-down is expressible.

What binning costs

The price is exact and it is paid on every row. The model can no longer distinguish 26 from 35, and it now treats 25 and 26 as maximally different. You deleted the ordering inside each bin and manufactured a discontinuity at an arbitrary cut point.

The three lines below are three ages fed through the same binned logistic model: what it does at the boundary, and what it does inside a bin.

age 25.9 -> bin 1 -> coefficient -0.31
age 26.1 -> bin 2 -> coefficient +0.44        a 0.75 logit jump across 0.2 years
age 35.0 -> bin 2 -> coefficient +0.44        identical to a 26-year-old

The jump is measured in logits, the log-odds scale a logistic regression works on. A move of 0.75 logits multiplies the odds of the outcome by e^0.75 = 2.1.

So two people about two and a half months apart in age get predictions whose odds differ by a factor of 2.1, while a 26-year-old and a 35-year-old get exactly the same prediction. Neither is a fact about the world; both are artifacts of where you put the cut.

When the trade is worth making

Binning is worth it only when the relationship is genuinely non-monotone and the model cannot represent non-monotonicity on its own.

That narrows the list to linear and logistic models, and essentially nothing else. A tree already learns cut points, chooses them by a supervised criterion rather than by quantiles, and is free to use more of them wherever the signal warrants.

Three ways to choose the cuts

If you do bin, the strategy determines what you are assuming. The third row is the only supervised one, which is also the only one that can leak.

Equal width assumes the data is spread evenly across its range, which skewed data violates by definition.

Equal frequency assumes the quantile boundaries are stable. They are not: refit on next month’s data and the cut points move, so the same customer can change bins without changing age.

WoE (weight of evidence) is supervised — it consults the label y when choosing where to cut, merging adjacent bins until the event rate rises or falls monotonically across them. That makes it a target-derived transform, so it carries the leakage risk of Categorical encoding and needs the same out-of-fold treatment described there.

5. Categorical encoding

So far every column has been numeric. A plan or zip column has to become numbers before any model can touch it. The right way to do that depends almost entirely on how many distinct values the column has, and one of the options is a common source of a fake offline score.

A categorical feature takes one of a finite set of unordered values: plan, zip_code, device_type.

Its cardinality is the number of distinct values it can take. This chapter writes it k. A plan tier has k = 3; US zip codes have k around 40,000.

Cardinality is the axis every decision below turns on. The flowchart routes to an encoding by asking about k first, then one follow-up question. The branch labels carry the thresholds.

flowchart TD
    C{"Cardinality k"} -->|"k <= ~15"| OH["One-hot"]
    C -->|"k in ~15..1000"| T{"Enough rows<br/>per category?"}
    C -->|"k > ~1000"| H{"Model type?"}
    T -->|"yes, 100+"| OH2["One-hot, or<br/>target encoding OOF"]
    T -->|"no, long tail"| TE["Target encoding<br/>+ smoothing, OOF"]
    H -->|"linear / streaming"| HA["Hashing"]
    H -->|"NN / factorization"| EM["Learned embedding"]
    H -->|"GBDT"| NA["Native categorical<br/>LightGBM / CatBoost"]

    style OH fill:#2d6a4f,color:#fff
    style HA fill:#2d6a4f,color:#fff
    style NA fill:#2d6a4f,color:#fff
    style TE fill:#bc6c25,color:#fff
    style EM fill:#7209b7,color:#fff

The chart reads as three bands.

Below about fifteen categories: one-hot, and stop thinking about it.

Between fifteen and a thousand: the question is whether you have enough rows per category. With a hundred or more rows behind each level, one-hot still works. With a long tail of levels seen two or three times, target encoding with smoothing is better.

Above a thousand: the question is model type. Hashing on the linear or streaming branch. A learned embedding for neural networks and factorization machines. Native categorical handling for gradient-boosted trees, which LightGBM and CatBoost both provide.

What the leaf colours mean

The colours grade one thing: how much of the label the encoding itself touches. That is the same axis as the “Leakage risk” column of the summary table at the end of this section.

• Green leaves never look at y at all, so there is nothing they can leak — one-hot, hashing, and a GBDT’s native categorical handling.
• Orange is target encoding, the one leaf computed directly from y, and therefore the one that will hand you a fake offline score if it is fitted even slightly wrong.
• Purple is the learned embedding. It is also fitted against y, but only inside the model’s own training loop, so it leaks exactly as much as the model does and no more.

One leaf is left unfilled: the middle band’s “one-hot, or target encoding OOF”. It names two answers from two different colours and cannot honestly take either.

One-hot: the cost is statistical, not just memory

One-hot encoding gives each category its own 0/1 column. k categories become k-1 dummy columns; one level is dropped as the baseline, and the intercept absorbs it.

The memory objection is the obvious one and usually the least important. Sparse matrices, which store only the non-zero entries, handle 50,000 columns fine.

The real cost is statistical, and it comes from one formula. For an orthogonal dummy (a column that does not overlap with any other, which one-hot columns satisfy by construction), the variance of the estimated coefficient is:

Var(beta_j) = sigma^2 / n_j

Two symbols there. n_j is the number of rows in category j. sigma is the residual noise — the part of y the model has not explained. That is not the same as sd(y), and confusing the two is a common slip.

The formula says the estimate for a category is built from that category’s rows and nothing else. Few rows means a wide, unreliable estimate.

Now put real numbers in. You have 50,000 zip codes and 100,000 rows, so the average category holds 100,000 / 50,000 = 2 rows. Then:

sd(beta_j) = sigma / sqrt(2) = 0.71 · sigma

The estimated effect of a zip code has a standard error equal to 71% of the residual standard deviation. Those coefficients are noise with a name on them.

Regularization will shrink them toward zero. That is correct behavior, and it is also a confession that the encoding gave the model nothing.

Assumptions. One-hot assumes each category has enough rows to estimate its own effect independently, and that categories share nothing with each other — every level is treated as equidistant from every other. It breaks exactly when cardinality is high and the distribution is long-tailed, which is the common case for identifiers.

Ordinal encoding: a false metric, and a smaller partition set

Ordinal encoding replaces each category with an integer.

Map {red: 0, green: 1, blue: 2} and you have told a linear model two things it will believe. First, that green sits exactly between red and blue. Second, that going red → green has the same effect on the prediction as going green → blue, because both are a step of 1 and the model has one coefficient.

Both are hard constraints on an unordered variable, and both are false.

For trees the damage is subtler but real, and it is a counting argument.

A tree splitting a true categorical can send any subset of levels left and the rest right. The number of ways to split k levels into two non-empty groups is 2^(k-1) - 1.

Given an arbitrary integer code, the tree can only cut the code line at some threshold x <= t, so it can only produce contiguous groups: {0} vs the rest, {0,1} vs the rest, and so on. That is k - 1 options.

At k = 10: 2^9 - 1 = 511 possible groupings collapse to 9. And which 9 you get depends on an arbitrary alphabetical ordering of the level names.

Assumptions. Ordinal encoding assumes the levels are genuinely ordered and roughly equally spaced. Reserve it for cases where that is true (small < medium < large), where the constraint is not a lie but a useful piece of prior knowledge.

Target encoding, and the exact leakage mechanism

Target encoding replaces category c with mean(y | c), the average label among the rows in that category. If 6.1% of customers in zip 94107 churned, then zip = "94107" becomes 0.061.

It compresses arbitrary cardinality into a single column and works well. It is also the most common source of a model that scores high offline and ships badly.

The mechanism, stated exactly

The encoded value for row i is computed from a set of rows that includes row i. So x_encoded[i] contains y[i].

Take the extreme case. For a category with exactly one row:

encoded = mean([y_i]) = y_i

The feature is a verbatim copy of the label. The model does not have to learn anything; it just reads the answer off the input.

With a long-tailed categorical, a large fraction of rows sit in small categories. So a large fraction of the column is literally the target, and gradient descent finds it in the first few steps.

What it looks like on the scoreboard

Three metric names first, because the block below uses all three.

• AUC (area under the ROC curve) is the ranking metric used here. It runs from 0.5 (a coin flip) to 1.0 (a perfect ranking), and it equals the probability that a randomly chosen positive row scores above a randomly chosen negative one.
• Cross-validation (CV) is the standard honesty check: split the training rows into K equal parts called folds, fit on K-1 of them, score on the one held out, rotate through all K, average the results.
• A held-out set is a block of rows set aside before any fitting and scored exactly once.

The same model, evaluated four ways, from most optimistic to reality:

train AUC          0.971      the model memorized y through the encoding
5-fold CV AUC      0.943      the encoder was fit before the split -- also leaked
held-out AUC       0.622      first honest number
production AUC     0.618

The second line is the one to stare at. Fitting the encoder on the full dataset and then cross-validating does not detect the leak. The CV test fold’s labels were already baked into the encoding it is being scored on, so “my CV said it was fine” is not a defense.

The only numbers reflecting reality are the last two, and they agree to within 0.004. Offline and online agreeing is the mark of an honest evaluation.

The fix, part one: out-of-fold encoding

Out-of-fold encoding (OOF) closes the loop directly. Split into K folds. The encoding applied to the rows in fold k is computed only from the folds that are not k.

Row i therefore gets an encoding that never saw y[i]. The copy-the-label path is cut.

The fix, part two: smoothing

Out-of-fold alone is not enough, because it still estimates a 2-row category from 2 labels. That estimate is honest but wildly noisy.

Smoothing shrinks each category’s mean toward the global mean:

enc(c) = (n_c·mean_c + m·global_mean) / (n_c + m)

Read it as a weighted average between the category’s own evidence (n_c rows saying mean_c) and a prior (m fictional rows saying global_mean).

That is the posterior mean under a conjugate prior. In plain terms: the prior is what you believe before seeing a category’s rows — here, the overall average label. The posterior is your belief after seeing them. Conjugate means the arithmetic of combining the two collapses neatly to that single weighted average instead of requiring an integral.

The parameter m has a real interpretation. Formally it is the ratio of within-category noise variance to between-category variance. Operationally it is the number of observations at which you trust a category as much as you trust the prior.

The table below sets global_mean = 0.05 and m = 20, so the prior is worth 20 fictional rows at the base rate. Each row shows a category with n_c rows and a raw average of mean_c, and what smoothing does to it.

Read the first and last rows together.

The singleton in row one had a raw encoding of 1.00 — a perfect copy of its label. Smoothing takes it to 0.095, against a base rate of 0.05. That is nearly no signal, which is exactly right, because one observation carries nearly no signal.

The 500-row category in row four keeps 0.290 of its raw 0.30. Five hundred observations do carry signal, and the prior’s 20 fictional rows barely move them.

Rows two and three show the crossover: at n_c = 5 the prior still dominates (0.60 pulled down to 0.16), and by n_c = 100 the category is mostly speaking for itself.

Both halves in code

The function below does OOF and smoothing together. The comment block inside the loop is the important part — it marks the subtle version of the same bug, where the prior rather than the category mean carries the leak.

import numpy as np
from sklearn.model_selection import KFold


def oof_target_encode(cat, y, m=20.0, n_splits=5, seed=0):
    """Out-of-fold smoothed target encoding. Returns (train_enc, fit_map)."""
    cat, y = np.asarray(cat), np.asarray(y, dtype=float)
    prior = y.mean()
    out = np.full(len(y), prior)

    for tr, va in KFold(n_splits, shuffle=True, random_state=seed).split(cat):
        # The prior must come from the TRAINING fold only. Using the global
        # y.mean() here leaks y[i] into out[i] through the smoothing term --
        # for a singleton category at m=20 that is 20/21 of the encoded value.
        # It is a 1/n-order leak, which is exactly the size that survives a
        # code review and then shows up as an unexplained CV-to-serving gap.
        prior_tr = y[tr].mean()
        stats = {}
        for c in np.unique(cat[tr]):
            sel = cat[tr] == c
            n_c, mean_c = sel.sum(), y[tr][sel].mean()
            stats[c] = (n_c * mean_c + m * prior_tr) / (n_c + m)
        out[va] = [stats.get(c, prior_tr) for c in cat[va]]   # unseen -> prior

    full = {}                                              # map for inference
    for c in np.unique(cat):
        sel = cat == c
        full[c] = (sel.sum() * y[sel].mean() + m * prior) / (sel.sum() + m)
    return out, (full, prior)

Three details in that code are the answers to the usual follow-up questions.

Unseen categories fall back to the prior, never to NaN or 0. Zero is a real predicted rate — “nobody in this category ever churned” — so filling with 0 is a confident wrong answer rather than an admission of ignorance.

The map used at inference is fit on all training rows. That looks like the leak you just removed, but it is not: at inference there is no row whose label could leak into itself, because the label does not exist yet. More data means a better estimate, so use all of it.

For time-ordered data the folds must be time-ordered too, not random (Temporal features and lookahead leakage). A random fold lets the future encode the past.

Assumptions. Three, and they get progressively less forgiving.

It assumes the relationship between a category and the label is stable over time. You are baking a historical rate into a feature, so if the rate moves, the feature is not merely noisy — it is stale, and confidently so.

It assumes enough rows per category that the shrunken estimate still carries some signal.

And it assumes the encoding is fit strictly out-of-fold. Violate that last one and every offline number you produce is fiction.

Hashing: why collisions do not hurt as much as they should

Hashing replaces the vocabulary with arithmetic. Map c -> hash(c) mod d: run the category name through a hash function and take the result modulo d, which lands it in one of d buckets.

You get four things for that. A fixed output width. No vocabulary to store or ship. Unseen categories handled for free, since the hash of any string is defined. And a single pass over the data, with no prior scan to build a dictionary.

The price is collisions: two different categories landing in the same bucket and therefore sharing one coefficient. Once k > d collisions are guaranteed by the pigeonhole principle.

Here is how bad it is for k = 50,000 categories hashed into d = 2^18 = 262,144 buckets.

P(a given category collides) = 1 - (1 - 1/d)^(k-1) ≈ 1 - e^(-k/d) = 1 - e^(-0.19) = 0.173

expected colliding pairs, ALL 50,000 categories = C(50000,2)/d = 1.25e9/262144 = 4,768
expected colliding pairs, the TOP 100           = C(100,2)/d   = 4,950/262144  = 0.019

Line by line.

The first line asks: for one particular category, what is the chance that at least one of the other k-1 = 49,999 categories lands on top of it? Each of them misses with probability 1 - 1/d, so all of them miss with probability (1 - 1/d)^49999, which is about e^(-50000/262144) = e^(-0.19) = 0.827. So 17.3% of categories collide with something.

The second and third lines count expected pairs. C(n,2) is the number of ways to pick 2 items from n, and each pair collides with probability 1/d. Across all 50,000 categories there are about 1.25 billion pairs, giving 4,768 expected collisions. Across just the 100 most frequent categories there are only 4,950 pairs, giving 0.019 — under a 2% chance that any two of the top 100 share a bucket.

17% of categories colliding sounds fatal. Here is why it is not:

Categorical frequency is Zipf-distributed, so almost every collision is between two rare categories whose coefficients were noise anyway — that is the one-hot variance argument from earlier in this section — while the chance that any two of the top 100 collide is under 2%.

A Zipf distribution is the steep popularity curve you see in real-world categories: a handful of values cover most of the rows, and a very long tail covers the rest. Because buckets are assigned uniformly at random, and because the tail is where nearly all the categories are, collisions overwhelmingly pair tail with tail.

Hashing degrades exactly the part of the feature that had no signal. That is a much better interview answer than “collisions are rare,” because collisions are not rare.

Assumptions. The argument depends on the frequency curve being Zipf-like. On a flat, uniform categorical, collisions land on informative levels as readily as uninformative ones, and the damage is real.

Hashing also assumes you never need to invert it. There is no asking “what is the coefficient for zip 94107” afterward, because you cannot recover which categories share a bucket.

Embeddings

An embedding is a learned dense vector per category: a short list of numbers — say 77 of them — standing in for zip = "94107". Those 77 numbers start random and are trained jointly with the model, the same way its weights are. Nobody sets them by hand.

The common heuristic for choosing the length is min(600, round(1.6 · k^0.56)). At k = 1,000 that is 1.6 · 1000^0.56 = 1.6 · 47.9 = 76.6, rounding to 77 dimensions.

That formula is the fast.ai rule of thumb — it ships as emb_sz_rule in fastai.tabular — and it is an empirical fit to what worked across their tabular benchmarks, not a quantity derived from theory. Treat it as a starting point to tune from, not a result.

Why embeddings beat one-hot comes down to sharing. One-hot forces every category to be equidistant from every other, so nothing learned about one level transfers to any other. An embedding places similar categories near each other in the numeric space, so a rare category borrows statistical strength from the frequent ones it resembles. The two-row zip code inherits most of its behaviour from the neighbouring zip codes that have thousands of rows.

That is the same lossy-compression-for-similarity property described in Embeddings and why dense search misses err_4021, and it comes with the same failure mode: an embedding is optimized for similarity, so it will not preserve the identity of a rare literal. Do not expect an embedding of sku_id (stock-keeping unit, a product identifier) to encode “this exact SKU” for a SKU seen 3 times.

Assumptions. Embeddings assume similarity between categories is meaningful and learnable.

They also assume you have enough labeled rows to fit k · dim extra parameters on top of the model itself. At k = 1,000 and dim = 77 that is 1,000 × 77 = 77,000 new parameters, all of which have to be estimated from your labels.

And they assume you have reserved an OOV (out-of-vocabulary) slot for categories never seen in training, because an unknown key has no vector to look up.

All six encodings side by side

One row per encoding. “Output width” is how many columns the encoding adds. “Handles unseen” is what happens at serving time when a category shows up that training never saw — the one column people forget until production.

6. Missing values: three mechanisms, three different correct answers

“What should I fill in?” is the wrong first question about a missing value. The right first question is “why is this value missing?” — there are three mechanisms by which a value can be absent, they take three different correct responses, and mean imputation, the default in most tutorials, damages a model in three separate, quantifiable ways.

The chart below branches on that one question — what does the missingness depend on? — and the three answers give the three mechanisms, each with its own recommended response in the box below it. The names are unhelpful, so plain-English definitions follow.

flowchart TD
    M{"Does missingness depend on...?"} -->|"nothing"| MCAR["MCAR<br/>sensor dropped a packet"]
    M -->|"other OBSERVED columns"| MAR["MAR<br/>income missing more<br/>often for the young"]
    M -->|"the MISSING VALUE itself"| MNAR["MNAR<br/>high earners decline<br/>to state income"]
    MCAR --> F1["Any imputation is unbiased.<br/>Listwise deletion also valid,<br/>just wasteful"]
    MAR --> F2["Conditional imputation<br/>+ missingness indicator"]
    MNAR --> F3["No imputation recovers it.<br/>Model the missingness<br/>as its own signal"]

    style MCAR fill:#2d6a4f,color:#fff
    style MAR fill:#bc6c25,color:#fff
    style MNAR fill:#9d0208,color:#fff

MCAR stands for missing completely at random. The absence has nothing to do with any value, observed or not — a sensor dropped a packet. Under MCAR, any imputation is unbiased in the sense that it does not systematically shift the column’s average. Listwise deletion, meaning you drop any row that has a missing value anywhere, is also valid here. Just wasteful.

MAR stands for missing at random, which is a badly chosen name: it does not mean random at all. It means the absence depends on other observed columns — income is missing more often for the young. Because the columns it depends on are ones you have, conditioning on them recovers what is missing.

MNAR stands for missing not at random. The absence depends on the missing value itself, as when high earners decline to state their income. No imputation recovers this, because the information you would need is precisely the information that is absent. The correct response is to stop trying and model the missingness as its own signal.

Why “impute with the mean” is usually wrong

Three separate damages, each with arithmetic behind it.

Damage 1: variance collapse

Take n = 1000 rows with 30% missing, a true mean of 50 and a true sd of 10. You impute all 300 missing values with 50, the mean.

Every imputed value now sits at exactly the mean, so it contributes zero to the squared deviations:

Var_imputed = (700·100 + 300·0) / 1000 = 70,000 / 1000 = 70
sd_imputed  = sqrt(70) = 8.37        (the truth is 10)

In general Var_imputed = (1-f)·Var_true, where f is the missing fraction. Here (1 - 0.3) · 100 = 70, which matches.

For a mean the damage is worse than that factor suggests, because two errors compound. The step that makes them compound is the definition of the standard error — the expected size of the gap between your estimate and the truth. For a mean, SE = sd / sqrt(n).

Mean imputation gets both pieces of that formula wrong, and in the same direction:

sd:  understated by sqrt(1-f)                                  = 0.837
n :  the formula divides by sqrt(n), but only n(1-f) rows carry
     information, so the divisor is too big by
     sqrt(n) / sqrt(n(1-f))  =  1/sqrt(1-f)                    = 1.195

SE = sd/sqrt(n)   ->   reported/true  =  sqrt(1-f) · sqrt(1-f)  =  (1-f)  =  0.70

Check it against the actual numbers. What you report is 8.37 / sqrt(1000) = 0.265. What is true, since only 700 rows carry information, is 10 / sqrt(700) = 0.378. The ratio is 0.265 / 0.378 = 0.70, exactly (1-f).

Both errors point the same way and each contributes one factor of sqrt(1-f). So your reported standard error is 70% of the truth. Confidence intervals come out 30% too narrow, and t-statistics — the estimate divided by its standard error, the number a significance test compares against a threshold — come out inflated by 1/0.70 = 1.43, or 43% too large.

Mean imputation manufactures statistical significance.

An aside: is a fitted regression spared?

Partly, and it is worth being precise about which part, because this is a common interview follow-up.

With a single mean-imputed predictor under MCAR, the coefficient is exactly unbiased. The cancellation is easy to see. Filling the gaps with the observed mean scales both Cov(x*, y) and Var(x*) by the same factor (1-f). And beta = Cov(x*, y)/Var(x*) is a ratio, so the factor divides out.

That is the only claim that survives.

Its standard error does not cancel. Simulating n = 1000, f = 0.3, y = x1 + noise over 40,000 replicates, the reported SE comes out at 1.139x the true sampling spread of beta_1 — so the interval is 14% too wide, not 30% too narrow.

The sign flipped relative to the mean case, and for a specific reason: in a regression the sqrt(n) term is not what broke. What broke is the residual variance, which got inflated by the rows whose x1 is now a constant and therefore explains nothing.

And the unbiasedness itself is fragile. It holds for one predictor and fails the moment a second correlated one exists — which is the next damage.

Damage 2: coefficient bias, even under MCAR

The common belief is that MCAR makes mean imputation safe. It does not, the moment you have two correlated predictors.

Set up the truth: y = 1.0·x1 + 1.0·x2 + noise, with x1 and x2 both standard normal and correlated at 0.6. Half of x1 is missing completely at random and imputed with its mean, which is 0 since the variables are centered. Write x1* for the imputed column.

Replacing half a centered column with 0 halves its variance and halves every covariance it takes part in. x2 is untouched. That asymmetry is the whole bug.

Var(x1*)     = 0.5 · 1.0 = 0.5        Cov(x1*, x2) = 0.5 · 0.6 = 0.3
Cov(x1, y)   = 1.0 + 0.6 = 1.6   ->   Cov(x1*, y)  = 0.5 · 1.6 = 0.8
Cov(x2, y)   = 0.6 + 1.0 = 1.6

beta = Sigma^-1 c,  Sigma = [[0.5, 0.3], [0.3, 1.0]],  det = 0.41
beta_1 = ( 1.0·0.8 - 0.3·1.6) / 0.41 = 0.32/0.41 = 0.780      (true 1.0)
beta_2 = (-0.3·0.8 + 0.5·1.6) / 0.41 = 0.56/0.41 = 1.366      (true 1.0)

Where each line comes from. Cov(x1, y) = 1.0 + 0.6 because x1 contributes its own variance of 1.0 through its coefficient, plus 0.6 through its correlation with x2. The 0.5 · factors are the (1-f) = 0.5 shrinkage from imputing half the column.

The last two lines solve the ordinary least-squares formula beta = Sigma^-1 c, where the coefficient vector equals the inverse of the predictors’ covariance matrix Sigma times the vector c of their covariances with y. For a symmetric 2x2 matrix [[a, b], [b, d]] the inverse is [[d, -b], [-b, a]] / det, which is where the cross-multiplied terms come from. The determinant here is 0.5 · 1.0 - 0.3 · 0.3 = 0.41.

x1’s effect is understated by 22% and x2 absorbs the difference, inflating by 37% — under the most benign missingness mechanism there is.

The reason is one line above: Cov(x1*, y) and Cov(x1*, x2) both shrink by (1-f) while Cov(x2, y) does not, so the regression reallocates the credit to the predictor whose numbers are intact. Have this example ready when someone tells you MCAR is fine.

Damage 3: under MNAR, the mean is the worst possible guess

If income is missing precisely because it is high, imputing the mean places every high earner at the center of the distribution.

That is maximally wrong, in a coordinated direction, on exactly the rows that matter most.

What to do instead

Each row below pairs a method with the condition that makes it correct. The first row is the only unconditional recommendation in the table — do it regardless of which mechanism you are facing.

Four terms in that table need unpacking.

Identifiability means the data contains enough information to pin down the answer uniquely. Adding the missingness indicator restores it because the model can now fit a separate effect for “we do not know,” instead of blending that case into the imputed value and pretending it is an observation.

MICE stands for multivariate imputation by chained equations. Predict each column that has missing values from the other columns, sweep through all of them repeatedly, and stop when the fills stop changing.

Multiple imputation runs that process m times with different random draws, fits the model on each completed dataset, and pools the results. That way the uncertainty about the fills survives into the reported standard errors instead of being silently discarded — which is exactly what single imputation throws away.

In the fourth row, p is the number of features.

Assumptions, stated plainly.

Mean imputation assumes MCAR and uncorrelated predictors. Damage 2 above shows the second condition is not optional.

MICE assumes MAR — that the observed columns contain what is needed to predict the missing one. Under MNAR it is confidently wrong, which is worse than being uncertain.

The missingness indicator assumes almost nothing. That is why it is the one universal recommendation here: at worst it is a column of near-constant bits that the model ignores.

One operational rule follows straight from the pipeline diagram in The pipeline and where the bugs live: fit the imputer on training folds only, after the split and after deduplication.

A median computed over train and test together is a test-set value entering the training features. It is a small leak, but a real one, and both duplicate rows and test rows silently move the statistic.

What interviewers probe: “You have 40% missing on your best feature — what do you do?” The wrong answer is a method name. The right answer starts with “first I’d find out whether it’s MCAR, MAR, or MNAR, by regressing the missingness indicator on the observed columns” — if that model has any predictive power at all, it is not MCAR, and the indicator is a feature.

7. Crosses and interactions

Everything so far has treated columns one at a time. Some signal only exists in combinations of columns — and an additive model, provably, cannot see a combination at all.

An interaction is when the effect of one feature depends on the value of another. A discount drives conversion for new users and does nothing for loyal ones — that is an interaction between discount and tenure.

A cross is the explicit feature you build to expose an interaction: the product of two numeric columns, or the pair of two categorical values treated as a single new category.

Why a linear model cannot see an interaction at all

A linear model is additive by construction. It adds up one contribution per feature and does nothing else.

That makes it unable to represent XOR, the pattern “exactly one of the two inputs is true.” The proof is two lines.

Suppose some f(x1, x2) = a·x1 + b·x2 + c classifies XOR correctly: negative when the two inputs match, positive when they differ. Then evaluate f at all four input combinations and add them in pairs.

f(0,0) < 0  and  f(1,1) < 0   ->   sum:  a + b + 2c < 0
f(0,1) > 0  and  f(1,0) > 0   ->   sum:  a + b + 2c > 0        contradiction

Walk the first line. f(0,0) = c and f(1,1) = a + b + c. Both are required to be negative, so their sum a + b + 2c is negative.

Now the second. f(0,1) = b + c and f(1,0) = a + c. Both are required to be positive, so their sum a + b + 2c is positive.

Both lines compute the same expression, and it cannot be negative and positive at once. So no a, b, c exists. This is not “linear models are bad at XOR” — it is impossible, for any data and any fitting procedure.

Add the cross term x1·x2 as a third input and the problem disappears:

f = x1 + x2 - 2·(x1·x2) - 0.5
   (0,0) -> -0.5   (0,1) -> +0.5   (1,0) -> +0.5   (1,1) -> -0.5     correct

Check the last one: 1 + 1 - 2·(1·1) - 0.5 = 2 - 2 - 0.5 = -0.5, negative, which is right because the inputs match. The cross term is what fires only when both inputs are 1, and that is the whole fix.

Trees find interactions, but they pay in depth — and depth costs data

Trees find interactions on their own, but pay for them in depth, and depth costs data.

Representing a k-way interaction requires k levels of splits, one per feature involved. Each level roughly halves the number of rows reaching a node, so the rows available to estimate the leaf value fall as n / 2^k.

n = 10,000, balanced splits
2-way interaction -> leaf sees ~2,500 rows     estimate is stable
4-way interaction -> leaf sees ~625 rows       usable
7-way interaction -> leaf sees ~78 rows        noise

Those are 10,000/4, 10,000/16 and 10,000/128. By the seventh level you are averaging 78 rows and calling it a prediction.

That is the honest version of “trees handle interactions automatically.” They handle low-order interactions automatically, and high-order ones only with a great deal of data.

It is also why explicit crosses still earn their place in ads and recommendation systems. There the important interactions are already known from the domain — user_country × item_category — so you hand them to the model directly and skip the depth the tree would otherwise spend rediscovering them.

The cardinality problem

Crossing two categoricals multiplies their level counts, and the products get large fast:

• zip (40,000) × device (5) = 200,000 levels
• zip (40,000) × hour (24) = 960,000 levels

Most of those levels have zero support, meaning no training row exhibits that combination at all, so nothing can be estimated for them. You have added 960,000 columns and real information for perhaps a few thousand.

This is exactly the setting hashing was built for (Categorical encoding): hash the cross, do not enumerate it.

Numeric crosses are cheaper. All pairwise products of p = 100 features is only C(100,2) = 4,950 columns. But they explode the same way once you go three-way.

Which crosses to try first

In order of hit rate:

1. Pairs the domain says interact. Cheapest to evaluate and most likely to work.
2. A high-importance numeric feature crossed with a low-cardinality segment. Small blast radius, and it directly asks “does this feature behave differently for this group?”
3. Ratios and differences of same-unit features. price / category_median_price beats having price and category_median_price as two separate columns, because the ratio is the interaction the model would otherwise have to learn.

Ratios rank high for a specific reason: a single division encodes a relationship that a tree needs several splits to approximate and a linear model cannot express at all.

Assumptions. A cross assumes the interaction is real and that you have enough rows in each combination to estimate it. Otherwise you have added a wide, sparse block of columns that regularization will spend its budget shrinking.

A ratio additionally assumes the denominator is never zero and is on a meaningful scale. price / category_median_price is only interpretable because both sides are prices.

8. Temporal features and lookahead leakage

Features built over time — recency, rolling counts, trends — all answer to a single rule, and four common traps break it. Each trap produces the same symptom: an offline score far better than anything production will ever show you.

The diagram shows the same three-stage timeline twice: a stretch of history, the features computed from it, and the label window whose outcome you are predicting. The only difference between the two rows is where the history stops. Compare the top-left box with the bottom-left box.

flowchart LR
    subgraph OK["Correct — as-of the prediction time"]
        direction LR
        H1["history<br/>t-90d .. t"] --> P1["features<br/>computed at t"]
        P1 --> L1["label window<br/>t .. t+30d"]
    end
    subgraph BAD["Leaked — the window crosses t"]
        direction LR
        H2["history<br/>t-90d .. t+30d"] --> P2["features<br/>see the label period"]
        P2 --> L2["label window<br/>t .. t+30d"]
    end

    style OK fill:#2d6a4f,color:#fff
    style BAD fill:#9d0208,color:#fff

In the correct version, the history ends at the prediction time t. The features computed at t see only that history, and the label window runs forward from t to t+30d. No overlap.

In the leaked version, the history runs to t+30d. The window crosses t, so the features see the very period whose outcome they are supposed to predict.

The rule: every feature value must be computable from data that was both (a) timestamped before the prediction time t and (b) actually available at t.

Those are two different conditions, and the second is the one that gets missed. A credit score can be about March 1 and still not have existed on March 1. Trap 3 below is that case.

Trap 1 — the unbounded aggregate

Consider a churn model whose labels are snapshotted at 2024-03-01, with a feature days_since_last_login. The query below looks completely ordinary. The comment marks the flaw.

-- LEAKED
SELECT user_id, DATEDIFF('2024-03-01', MAX(login_ts)) AS days_since_last_login
FROM logins GROUP BY user_id;                 -- MAX over the WHOLE table

The query asks for the most recent login per user with no upper bound on time. So a user who logged in on 2024-04-10 — six weeks after the label date — has that row picked up by MAX, and DATEDIFF returns a negative days-since value.

The model duly learns that negative or tiny values mean “not churned.” That is true, and it is useless, because at scoring time no future login exists to produce a negative number.

The damage, measured:

offline AUC, random 80/20 split   0.941
production AUC, first week        0.683

The fix is one clause: WHERE login_ts < '2024-03-01'. It must appear in every feature query, not just the ones that look suspicious.

Trap 2 — rolling windows without a shift

A rolling window computes a statistic over the last n rows at each point in a time series. By default, that window includes the current row — which is the bug.

# LEAKED: today's value is inside today's mean
df["avg_7d"] = df["amount"].rolling(7).mean()

# CORRECT: the window ends yesterday
df["avg_7d"] = df["amount"].rolling(7).mean().shift(1)

The leaked version hands each row one seventh of its own target-correlated value. One seventh sounds small; it is more than enough for the model to reconstruct the answer, since the other six sevenths are also highly correlated with it.

On a daily-revenue forecast that alone can move R^2 — the fraction of the target’s variance the model explains, where 1.0 is perfect — from 0.40 to 0.92 offline, while changing nothing at all in production.

.shift(1) moves the whole window back one period, so the mean for a given day is computed from the seven days before it.

Trap 3 — restated data, and why you need two timestamps

A credit score for 2024-03-01 is revised on 2024-04-15, and your warehouse overwrites the old value, holding the revised number under the original date. Query it and you get a number nobody could have known on March 1.

Notice that the timestamp is not lying. The score really is about March 1. It just did not exist yet.

The fix is bitemporal storage, meaning every fact carries two dates rather than one:

• valid_time — when the fact was true in the world.
• transaction_time — when your system learned it.

A point-in-time join then filters on both: valid_time <= t AND transaction_time <= t. That retrieves what was both true and known at t.

A store that keeps only valid_time cannot be made correct by any amount of query care. The information about when you learned each value was never recorded, so no query can recover it. This is a schema decision, not a SQL decision.

The same trap turns up in a friendlier disguise. A nightly job that lands at 06:00 is not available to a 05:30 request, so serving quietly falls back to yesterday’s value while training used today’s. That is a skew bug (Feature stores and trainingserving skew) wearing a leakage bug’s clothes.

Trap 4 — the split itself

The train/test split is how you decide which rows the model fits on and which it is scored on.

Random k-fold cross-validation on temporal data trains on the future and tests on the past. No deployed model ever gets to do that.

Here is the same model under three evaluation schemes. RMSE (root mean squared error) is the typical size of a prediction error in the target’s own units, so lower is better.

random 5-fold CV        RMSE 0.31        wrong by construction
forward-chaining CV     RMSE 0.58        matches reality
production, month 1     RMSE 0.61

The random-CV number is 0.31 and reality is 0.61 — off by a factor of two, in the flattering direction. The forward-chaining number is 0.58, within 0.03 of production.

Forward chaining — train on [0, k], test on [k+1], then advance and repeat — is the only scheme whose number you can quote to a stakeholder. It is the only one that reproduces the shape of the real task: fit on the past, predict the future.

One addition is required. Put a purge gap between the training window and the test window, equal to the label horizon. If your label takes 30 days to resolve, the last 30 days of the training window overlap the test window’s label period, so those rows must be dropped.

Assumptions. Every temporal feature assumes your event timestamps mean what you think they mean — event time, not ingestion time — and that the arrival delay you saw historically is the delay you will see live.

Forward chaining additionally assumes the relationship is stable enough that the past predicts the future at all. Under a sharp regime change, even a correct evaluation overstates what the model will do next month.

Temporal features worth building

These five carry most of the signal in churn, fraud and demand problems. The seasonality row is the one that trips people up.

Feeding the raw hour as an integer tells the model that 23:00 and 00:00 are 23 units apart, when in fact they are one hour apart. On an integer axis, midnight is maximally far from 11pm — the worst possible arrangement for any pattern that runs across midnight.

Mapping the hour onto a circle fixes it. Take sin(2·pi·h/24) and cos(2·pi·h/24) as two features. At h = 23 and h = 0 those two points sit next to each other on the circle, which is what the clock actually does. You spend one extra column and get correct adjacency everywhere, including the wraparound.

9. Text and image features, briefly

Tabular pipelines regularly pick up a text column or an image column along the way, and in both cases the right first move is the cheap classical one, not the expensive modern one.

Text

TF-IDF (term frequency–inverse document frequency) is still the correct first move for classification over a fixed vocabulary:

tfidf(w, d) = tf(w, d) · log(N / df(w))

It scores a word w in a document d as two factors multiplied together. tf(w, d) is how often the word appears in that document. log(N / df(w)) is how rare the word is across the whole collection, where N is the number of documents and df(w) is the number of documents containing w.

The second factor is the interesting one. Put numbers in:

a term in 3 of 100,000 documents:      log(100,000/3) = log(33,333) = 10.4
a term in half of them:                log(2)                       =  0.69

A term appearing in 3 documents is weighted 15 times more heavily than one appearing in 50,000. That is the IDF term doing its job, and it does the same job as BM25, the standard keyword-ranking function in search engines described in Embeddings and why dense search misses err_4021.

Why chapter 00 shows 10.26 for the same example

Follow that link and you will meet 10.26 where this chapter says 10.4. Both are right; they are different formulas.

Plain TF-IDF, used here, is log(N/df), giving log(33,333) = 10.41.

BM25, used there, applies a smoothed variant: ln((N − df + 0.5)/(df + 0.5) + 1) = ln(28,572) = 10.26. The half-count corrections exist so the weight stays finite and non-negative when df reaches 0 or N.

At N = 100,000 and df = 3 the smoothing moves the answer by 1.5%. On the point both chapters are making — a rare term outweighs a common one by an order of magnitude — it changes nothing.

Rare literals are up-weighted, which is exactly the property dense embeddings lack. Same complementarity as in retrieval, one level down the stack.

Character n-grams are the other standard trick: every contiguous run of 3 to 5 characters, so error yields err, rro, ror, erro, rror, error. Because a typo only corrupts the few n-grams that touch it, most of the representation survives — which is why n-grams are the default for short, noisy text like product titles.

Which to use

Sentence embeddings when paraphrase matters and exact wording does not. TF-IDF when the discriminating tokens are literals: error codes, SKUs, model numbers.

Concatenating a 300-dimensional embedding with the top-5,000 TF-IDF columns beats either alone, for exactly the reason hybrid retrieval beats either alone — the two methods fail on different inputs.

Images

Raw pixels are a poor feature for a classical model. Shift an image one pixel to the right and every value changes while nothing about its meaning does.

The standard move is a pretrained backbone as a frozen feature extractor. Two words there: the backbone is a large network already trained on millions of images, and frozen means you do not update its weights at all.

Take the penultimate layer — the second-to-last layer, whose output is the network’s compressed description of the image, typically 512 to 2,048 numbers — of a ResNet or a ViT (vision transformer, the attention-based architecture that replaced convolutional networks for many vision tasks). Then fit a logistic regression or a GBDT on those numbers.

With 2,000 labeled images this beats fine-tuning the whole network, and the reason is a parameter count: you are estimating 2,048 parameters instead of 25 million, from the same 2,000 rows.

Assumptions. The frozen-backbone recipe assumes your images resemble the ones the backbone was pretrained on. On natural photographs that holds; on medical scans, satellite imagery or engineering diagrams the pretrained description may not separate your classes at all, and fine-tuning stops being the wasteful choice.

10. Feature stores and training/serving skew

The most common production ML bug, by a wide margin, lives not in any single transform but in the fork at the end of the pipeline — and it is the hardest bug to see, because nothing raises an exception and no monitor fires.

Training/serving skew is the condition where the feature vector computed during training differs from the one computed at serving, for the same entity at the same moment.

The diagram below shows why that happens architecturally. There are two vertical paths that never touch, and one dotted line between them. The dotted line is the bug.

flowchart TD
    S1[("Warehouse<br/>full history")] --> B["Batch job<br/>Spark SQL"]
    B --> OFF[("Offline store<br/>point-in-time correct")]
    OFF --> TRAIN["Training"]

    S2[("Event stream<br/>last 30d")] --> R["Streaming job<br/>Java / Flink"]
    R --> ON[("Online store<br/>Redis, p99 < 10ms")]
    ON --> SERVE["Serving"]

    B -. "same NAME<br/>different CODE" .-> R

    style OFF fill:#2d6a4f,color:#fff
    style ON fill:#bc6c25,color:#fff
    style TRAIN fill:#2d6a4f,color:#fff
    style SERVE fill:#bc6c25,color:#fff

The diagram is one feature definition realized twice.

The training path. A warehouse holding full history feeds a batch job — typically Spark SQL, a distributed engine for large scans. That job writes an offline store that is point-in-time correct, meaning each row records the feature values as of that row’s own prediction time.

The serving path. An event stream carrying only the last 30 days feeds a streaming job written in Java or Flink. That job writes an online store, usually Redis — an in-memory key-value database fast enough that its p99 (the latency 99% of requests come in under) is below 10 milliseconds.

Two independent implementations of one definition, written by different people in different languages against different data. The dotted line, where the same feature name is produced by different code, is where the money goes.

A concrete trace

Take the feature avg_order_value_30d. Its two implementations differ by one word, and they shipped six months apart.

Read the block in three parts: the two queries, then what they produce for two customer segments, then what that does to the metrics.

Spark (training):   SUM(amount)/COUNT(*)  WHERE status IN ('completed','refunded')
Flink (serving):    SUM(amount)/COUNT(*)  WHERE status = 'completed'

customers with no refunds (88%):  training  84.20   serving  84.20    identical
customers with refunds    (12%):  training  71.40   serving 133.60    +87%

overall AUC          0.812 -> 0.794     -0.018    "looks like normal drift"
AUC on refund slice  0.791 -> 0.712     -0.079    nobody was looking at this
approval rate         5.0% ->  5.4%     +8% relative

The mechanism is in the word 'refunded'. Refunded orders are small — a refunded $12 item drags the average down in training, and is simply absent at serving. So for the 12% of customers with refunds, serving computes 133.60 where training computed 71.40, an 87% inflation.

Those customers are disproportionately the high-risk ones the model exists to catch, and the model is being handed a rosier number for exactly them.

The aggregate metric absorbs a subpopulation failure, because 88% of the population is unaffected and drags the average back. Overall AUC moved 0.018 — noise, by most standards. The refund slice moved 0.079, four times as far.

The business cost follows. At 2M decisions a day, the 0.4-point rise in approval rate is 8,000 additional approvals daily on the worst-scoring segment. It was found six weeks later, by the fraud-loss report.

Every guard in a normal pipeline passes. The feature is present, non-null, numeric, inside its historical min/max, and stable over time.

That last one is the key. Drift monitoring cannot find this class of bug, because drift detectors compare today against yesterday, and this feature has been identically wrong since the day it launched. There is no change to detect.

The three kinds of skew

Skew has three sources, and each one is detected by a different check and closed by a different fix. A monitor built for any one of them is blind to the other two, which is why the table has separate Detection and Fix columns rather than a single recommendation.

An as-of join, in the last row, means retrieving the value a column held at a specific past moment rather than the value it holds now. Training does this by construction. Serving, reading a live key-value store, does not.

The fixes, in the order they actually pay off

The four fixes below are ranked by how much of the failure class each one removes, not by how easy each is to start.

The split that matters: the first two prevent skew, the last two only detect it. That is the argument for not beginning with the monitoring everyone reaches for first.

1. Log the served feature vector and train on it. This eliminates implementation skew by construction: there is exactly one producer, so the training data is by definition what serving computed. Cost: you cannot train until you have logged for a while, and you can only train on features you already serve. It is first because it removes the failure class rather than detecting it.
2. One definition, one engine. A feature store where the transformation is declared once and both the batch and streaming paths are generated from it. Weaker than (1), because the two generated paths can still differ in engine semantics — null handling, integer division, timezone. But it works from day one and it covers features you have not launched yet.
3. Skew monitoring. For a sample of entities, recompute the offline feature at the exact serving timestamp and diff it against what was logged online. Alert on the p99 absolute difference, not the mean. In the trace above, the mean difference across all customers is 0.12 · 62 = +7.4 — small enough to look like noise — while the p99 is +62. Averaging over a mixed population cancels exactly the signal you need.
4. Per-slice metrics. A slice is a subpopulation you evaluate separately. Overall AUC moved 0.018; refund-slice AUC moved 0.079. If your dashboard has one number on it, it will not show you this.

What interviewers probe: “Your offline AUC is 0.92 and online it’s 0.71 — walk me through it.” Ordered diagnosis, and the order is by base rate: (1) leakage in a feature or the split (Categorical encoding, Temporal features and lookahead leakage) — most common, and it makes the offline number wrong; (2) training/serving skew — second most common, makes the online number wrong; (3) distribution shift, meaning the live data no longer resembles the training period; (4) the eval itself. Check leakage first because it is the only one you can confirm without any production data, and skew second because a single recompute-and-diff confirms or clears it.

11. Feature selection

Building features is generative work; the pipeline’s last box is the opposite — getting from thousands of candidate columns down to the few dozen worth serving. The cheapest way to do it is blind to the most valuable features, and one particular way of doing it wrong produces the most convincing fake result in all of applied ML.

The methods divide into three families by how much model fitting each one requires. Throughout, p is the number of candidate features and k is the number you keep.

A filter scores each feature against the label on its own, with no model involved. Mutual information measures how much knowing one variable reduces uncertainty about the other, in bits. Chi-square and the ANOVA F statistic are significance tests, for a categorical and a numeric predictor respectively.

An embedded method gets selection for free from a single model fit. L1 regularization is the penalty that drives weak coefficients exactly to zero, so the survivors are the selection. Tree gain is the total improvement in the split criterion attributable to each feature.

A wrapper repeatedly refits the model on candidate subsets and keeps the subset that scores best. RFE (recursive feature elimination) is the standard version: fit, drop the weakest features, refit, repeat.

Why filters miss interactions, exactly

Take XOR again: y = x1 XOR x2, with both inputs uniform on {0,1}.

Look at x1 on its own. When x1 = 0, y is 1 half the time (whenever x2 = 1). When x1 = 1, y is also 1 half the time. So:

P(y=1 | x1=0) = P(y=1 | x1=1) = 0.5   ->   MI(x1 ; y) = 0
P(y=1 | x2=0) = P(y=1 | x2=1) = 0.5   ->   MI(x2 ; y) = 0
MI({x1, x2} ; y) = 1 bit

Both individual scores are exactly zero — not small, zero. Knowing x1 alone tells you literally nothing about y. But the pair determines y perfectly, so the joint mutual information is a full bit.

A mutual-information filter therefore deletes both of the only two real features and keeps whatever noise column happened to score 0.001.

Marginal relevance and joint relevance are different quantities, and a filter only measures the first.

Why wrappers are expensive, exactly

Forward selection from p = 200 candidates down to k = 20 keepers evaluates every remaining candidate at each of the 20 rounds. Round 1 tries 200 features, round 2 tries the 199 left, and so on down to 181:

200 + 199 + ... + 181 = 3,810 model fits
3,810 fits x 5 folds x 30 s per fit = 571,500 s = 159 hours

That is the derivation behind “wrappers do not scale,” and it also names the lever. RFE removes a block of features per round instead of one, cutting the fit count by the block size.

Two traps in embedded selection

Tree gain importance is biased toward high-cardinality and continuous features. More distinct values means more candidate split points, which means more chances to fit noise. That is the multiple-comparison mechanism explained in Decision trees: test enough candidate thresholds and one of them looks good by chance.

Permutation importance on a held-out set is the honest version. Shuffle one column’s values at random, rescore the already-fitted model, and see how much the metric drops. If the column mattered, breaking it hurts.

But permutation importance has its own trap: it splits credit between correlated features. Shuffle one and the other covers for it, so the metric barely moves, so both look unimportant — and you delete both. The fix is to permute correlated groups together rather than one column at a time.

The leakage trap that produces the most convincing fake result

Selecting features on the full dataset and then cross-validating.

Here is why it works so well. Take n = 100 rows and p = 5,000 columns of pure noise, with no relationship to y whatsoever.

sampling sd of a correlation at n = 100:   1/sqrt(n-3) = 1/sqrt(97) = 0.102
E[max |z|] over 5,000 standard normals:    3.85
best noise feature by chance:              3.85 x 0.102 = 0.39

Step by step. Even with zero real signal, a sample correlation is not exactly zero — it scatters around zero with a standard deviation of about 0.102 at this sample size. Draw 5,000 of them and the largest in absolute value lands around 3.85 standard deviations out (E[max |z|] over 5,000 standard normals, by numerical integration and confirmed by simulation).

So the best of your 5,000 noise columns shows |r| ≈ 0.39 purely by chance. That looks like a real feature.

Now select the top 10 using all the data, then run 5-fold CV. You will report an accuracy well above chance from a dataset containing zero signal — because the selection step already consulted the labels of every fold, including the ones CV is about to “hold out.”

Feature selection is model fitting, so it goes inside the CV loop, not before it.

The practical ordering, and why it is this order

Five steps, arranged so the expensive methods only ever see small feature sets. Cost per step rises by roughly 10x each time.

1. Drop constants and duplicates. Free, and no risk of dropping anything useful.
2. Drop features you cannot serve. A feature that needs a 2-second join is not a feature, regardless of its importance score.
3. Use L1 or tree importance to build a shortlist. One model fit takes you from 5,000 candidates to about 200.
4. Run permutation importance on held-out data for an honest ranking. This costs p re-scorings of the already fitted model and zero refits — which is exactly why it is affordable at 200 features and not at 5,000.
5. Run a wrapper search on that final shortlist only, where O(p·k) is survivable.

Assumptions. Every selection method assumes the rows you select on are representative of the rows you will serve, and that the feature set is stable enough to be worth freezing.

Under distribution shift, a feature dropped today for low importance may be the one that matters next quarter. That is the argument for re-running selection on a schedule rather than treating it as a one-time step.

Cheat sheet

Each row is a symptom you can observe, the mechanism that produces it, and the fix — the whole chapter compressed into the order you will actually meet these problems.

Next: 02 — Classical Models — the models these features feed, and why each one cares about the transforms above.