Supervised learning

This is the ordinary case that most applied machine learning still is: you have examples where the answer is known, and you want something that gives the answer for cases where it is not. The section covers the two shapes that problem takes, predicting a number and predicting a category, the handful of algorithms worth knowing for each, and why gradient-boosted trees remain the first thing to reach for on data that lives in a table.

What you have                                   Reach for first
a table, predict a number                       gradient boosting regressor
a table, predict a category                     gradient boosting classifier
a table, and someone must sign off the model    logistic or linear regression
under a few hundred rows                        regularised linear model
text, images, audio                             a pretrained network, fine-tuned
you need a number to beat before anything else  predict the mean, or the majority class

Two problems#

Supervised learning takes two shapes. Regression predicts a number: a price, a duration, a temperature. Classification predicts a category: spam or not, which of five classes. The distinction decides your loss function and your metrics, and nearly everything follows from it.

Regression is scored by how far off the number is, usually mean squared error during training and mean absolute error when reporting to anyone who has to act on it. Classification is scored on whether the category was right, and trained on cross-entropy, which is the loss that punishes being confidently wrong. Picking the wrong pair is not a subtle error: a classifier trained with squared error learns to hedge towards the middle instead of committing, because squared error treats a confident mistake far too gently.

Linear regression#

Linear regression assumes the answer is a weighted sum of the features plus a constant, and finds the weights that minimise squared error. There is a closed-form solution, so there is no training loop and no learning rate; it fits in milliseconds on data that would take a network an hour.

The reason to reach for it is that every coefficient has a unit and a meaning. If size is measured in square metres and the answer in thousands of pounds, a coefficient of 3.1 says one more square metre is worth three thousand one hundred pounds, holding everything else fixed. That last clause is the whole of the difficulty, and the next block is what happens when it stops being true.

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 2000
size = rng.normal(90, 25, n)                  # square metres
age = rng.integers(0, 60, n).astype(float)    # years
price = 3.1 * size - 1.4 * age + 40 + rng.normal(0, 8, n)   # thousands

lin = LinearRegression().fit(np.c_[size, age], price)
print(lin.coef_.round(2))        # [ 3.11 -1.38]  -> recovers 3.1 and -1.4
print(round(lin.intercept_, 1))  # 39.0          -> recovers 40
print(round(lin.score(np.c_[size, age], price), 3))   # 0.991

Now add a third feature that is the first one in different units. The model is just as accurate and the coefficients become nonsense.

size_sqft = size * 10.764 + rng.normal(0, 0.5, n)    # the same thing, in square feet
Xc = np.c_[size, age, size_sqft]

for seed in (0, 1, 2):
    idx = np.random.default_rng(seed).choice(n, n, replace=True)   # resample the rows
    m = LinearRegression().fit(Xc[idx], price[idx])
    print(f"size {m.coef_[0]:+8.2f}   size_sqft {m.coef_[2]:+7.2f}   R2 {m.score(Xc, price):.3f}")

# size   +11.37   size_sqft   -0.77   R2 0.991
# size    +5.49   size_sqft   -0.22   R2 0.991
# size   +11.14   size_sqft   -0.75   R2 0.991

Linear regression fails in four ways, in the order you will meet them. Correlated features, as above: predictions fine, coefficients meaningless, and you only notice if you refit. A genuinely curved relationship, which it fits with a straight line and is wrong at both ends. Outliers, because squared error weights a point ten times too far by a hundred. And extrapolation: it will confidently price a 400 square metre house from data that stops at 200, because nothing in the model knows the data ran out.

Logistic regression#

The same weighted sum, passed through a function that squashes any number onto the range zero to one so it can be read as a probability. The sum itself is the log-odds, and that is what makes the coefficients readable: a coefficient of minus 0.075 means one more unit multiplies the odds by exp(-0.075), which is 0.93, a seven per cent reduction in the odds per unit. Ten units multiply the odds by 0.47, roughly halving them.

from sklearn.linear_model import LogisticRegression

income = rng.normal(50, 15, n)
odds = np.exp(-0.08 * income + 3.0)
default = (rng.random(n) < odds / (1 + odds)).astype(int)

log = LogisticRegression().fit(income.reshape(-1, 1), default)
b = log.coef_[0, 0]
print(round(b, 4))                  # -0.0748   recovers the -0.08 it was built from
print(round(np.exp(b), 4))          #  0.9279   odds multiplier per extra unit of income
print(round(np.exp(10 * b), 4))     #  0.4733   odds roughly halve per 10 units

It fails in the same four ways as linear regression, plus one of its own: the decision boundary is still a straight line through feature space, so any problem where the classes are separated by a curve needs the curve supplied as a feature. It is also the only model here that gives well-calibrated probabilities without extra work, which is why it survives in credit scoring and medicine long after it stopped being the most accurate option.

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

model = LogisticRegression(max_iter=1000).fit(X_train, y_train)
model.score(X_test, y_test)
model.coef_                 # which features mattered, and in which direction

Decision trees#

A tree asks one yes-or-no question about one feature at a time, and each answer sends the example down a branch to another question. It handles curves, thresholds and interactions between features without being told they exist, needs no scaling, and can be read aloud as a set of rules.

A split is chosen by brute force. At each node the tree tries every feature and every threshold in the data, scores each candidate by how pure the two resulting groups are, and keeps the best. Purity for classification is usually the Gini impurity: the chance of being wrong if you guessed a label at random from the group's own distribution. A group that is all one class scores zero; an even split of two classes scores 0.5.

def gini(y):
    if len(y) == 0:
        return 0.0
    p = np.bincount(y, minlength=2) / len(y)
    return 1.0 - (p ** 2).sum()

def best_split(x, y):
    parent = gini(y)
    best_gain, best_t = 0.0, None
    for t in np.unique(x)[:-1]:                 # every value is a candidate threshold
        left, right = y[x <= t], y[x > t]
        child = (len(left) * gini(left) + len(right) * gini(right)) / len(y)
        if parent - child > best_gain:
            best_gain, best_t = parent - child, t
    return best_gain, best_t

x = np.array([22, 25, 29, 31, 35, 38, 42, 47, 51, 60])
y = np.array([ 0,  0,  0,  1,  0,  1,  1,  1,  1,  1])

gain, t = best_split(x, y)
print(round(gini(y), 3))                 # 0.48  -> the parent node, five of each class
print(round(gain, 3), int(t))            # 0.32 35  -> split at age <= 35

threshold   left labels          right labels     weighted gini   gain
age <= 25   0 0                  0 1 0 1 1 1 1 1      0.300        0.180
age <= 29   0 0 0                1 0 1 1 1 1 1        0.171        0.309
age <= 35   0 0 0 1 0            1 1 1 1 1            0.160        0.320   <- chosen
age <= 38   0 0 0 1 0 1          1 1 1 1              0.267        0.213
age <= 47   0 0 0 1 0 1 1 1      1 1                  0.400        0.080

A tree left to run stops only when each leaf is pure, which on real data means one leaf per training row and a model that has memorised rather than learned. The controls are max_depth, min_samples_leaf and max_leaf_nodes, and they are the entire defence. A single tree constrained enough not to overfit is usually too weak to be interesting, which is why nobody ships one and why the next two sections exist.

Random forests#

Grow hundreds of deep trees, each on a bootstrap resample of the rows, and at each split let the tree consider only a random subset of the features. Average their predictions. Each tree individually overfits; the errors they make are different errors, and averaging cancels them.

The arithmetic of why is worth knowing, because it also says where the limit is. Averaging n predictions with variance v and correlation rho between them gives variance rho times v plus (1 minus rho) times v over n. The second term vanishes as you add trees. The first does not.

def ensemble_var(n_trees, rho, var=1.0):
    return rho * var + (1 - rho) * var / n_trees

for n in (1, 10, 100, 1000):
    print(n, round(ensemble_var(n, 0.0), 3), round(ensemble_var(n, 0.6), 3))

# 1     1.000   1.000
# 10    0.100   0.640
# 100   0.010   0.604
# 1000  0.001   0.600   <- adding trees past this point buys nothing

So the number of trees is not a parameter to tune, only a budget: more is never worse, and past a few hundred it stops helping. The parameter that matters is max_features, the size of that random subset, because it sets rho directly. Forests are hard to overfit, need almost no tuning, and land a little behind well-tuned boosting, which is a good trade when nobody has time to tune anything.

Gradient boosting#

Train a small tree, look at what it got wrong, train another tree on those errors, repeat. Each one corrects its predecessors. This is the strongest general approach for tabular data and remains so despite everything that has happened in deep learning.

The whole method fits in a dozen lines. Start from the best constant, which for squared error is the mean. At each round compute the residual, which for squared error is exactly the negative gradient of the loss. Fit a shallow tree to that residual, and add a small fraction of its prediction. That fraction is the learning rate.

from sklearn.tree import DecisionTreeRegressor

rng = np.random.default_rng(1)
X = rng.uniform(-3, 3, (400, 1))
y = np.sin(X[:, 0] * 2) + 0.3 * X[:, 0] ** 2 + rng.normal(0, 0.1, 400)

pred = np.full(len(y), y.mean())          # the best constant, before any tree
trees, lr = [], 0.1
for round_ in range(200):
    residual = y - pred                   # for squared error, this IS the negative gradient
    stump = DecisionTreeRegressor(max_depth=2).fit(X, residual)
    pred += lr * stump.predict(X)         # a small step in the direction it points
    trees.append(stump)
    if round_ in (0, 50, 100, 199):
        print(round_, round(np.mean((y - pred) ** 2), 4))

# round   0  mse 1.0587
# round  50  mse 0.0186
# round 100  mse 0.0076
# round 199  mse 0.0058

Because each step is small and deliberate, a boosted model with too many rounds overfits, unlike a forest. Early stopping on a validation split is the standard defence and is built in.

from sklearn.ensemble import HistGradientBoostingClassifier

model = HistGradientBoostingClassifier(
    max_iter=300,
    learning_rate=0.1,
    early_stopping=True,
).fit(X_train, y_train)

Hyperparameter       What it does                              Where to start
learning_rate        size of each correction                    0.05 to 0.1
max_iter             number of trees                            as high as you can afford
max_depth            interaction order each tree can capture    3 to 6, or leave unset
min_samples_leaf     refuses splits on too few rows             20, raise if overfitting
l2_regularization    shrinks leaf values                        0, raise if overfitting
max_features         fraction of columns per split              1.0, lower on wide data

from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import train_test_split

Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0)
for lr_, iters in [(1.0, 20), (0.3, 60), (0.1, 200), (0.03, 700)]:
    m = HistGradientBoostingRegressor(learning_rate=lr_, max_iter=iters,
                                      early_stopping=False, random_state=0).fit(Xtr, ytr)
    print(lr_, iters, round(np.mean((yte - m.predict(Xte)) ** 2), 4))

# 1.0    20 trees   test mse 0.0321
# 0.3    60 trees   test mse 0.0293
# 0.1   200 trees   test mse 0.0289   <- best
# 0.03  700 trees   test mse 0.0291   <- 3.5x the compute, very slightly worse

As for which boosting library: HistGradientBoosting ships with scikit-learn and needs nothing installed, which is why it is used throughout here. LightGBM is usually the fastest on wide data, XGBoost is the one most job descriptions name, and CatBoost handles high-cardinality categorical columns with the least work. They are within a point or two of each other on most datasets; the choice rarely decides a project and never decides an interview.

An end-to-end fit#

The order matters more than any single model. Baseline first, so every later number means something; a linear model next, because it is fast and tells you whether the problem is even hard; then the ensembles. Everything compared by cross-validation on the training set, and the test set touched once.

from sklearn.datasets import load_breast_cancer
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True, as_frame=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=0)      # stratify keeps the class balance

models = {
    "majority class":    DummyClassifier(strategy="most_frequent"),
    "logistic":          make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)),
    "random forest":     RandomForestClassifier(n_estimators=300, random_state=0),
    "gradient boosting": HistGradientBoostingClassifier(random_state=0),
}
for name, m in models.items():
    print(name, round(cross_val_score(m, X_train, y_train, cv=5, scoring="roc_auc").mean(), 4))

# majority class      0.5000
# logistic            0.9945
# random forest       0.9956
# gradient boosting   0.9923   <- boosting loses here, and that is the point of the table

Boosting coming third is not a mistake in the code. This dataset has thirty smooth, correlated measurements and 569 rows, which is close to the ideal case for a linear model and close to the worst case for trees. Reaching for boosting by reflex and reporting 0.992 would have left a simpler, faster, more explainable model that scores higher on the table, and nobody would ever have known.

logit = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)).fit(X_train, y_train)
coef = logit[-1].coef_[0]
for i in np.argsort(-np.abs(coef))[:3]:
    print(f"{X.columns[i]:22} {coef[i]:+.2f}  odds x{np.exp(coef[i]):.2f} per SD")

# radius error           -1.32  odds x0.27 per SD
# worst radius           -0.98  odds x0.37 per SD
# worst texture          -0.97  odds x0.38 per SD

from sklearn.metrics import roc_auc_score
best = HistGradientBoostingClassifier(random_state=0).fit(X_train, y_train)
print(round(roc_auc_score(y_test, best.predict_proba(X_test)[:, 1]), 4))    # 0.9868

Why trees beat neural networks on tables#

• Tabular features are heterogeneous and often uninformative individually. Trees select; networks must learn to ignore.
• Trees are invariant to monotonic transformations of a feature, so scaling and skew do not matter.
• They handle missing values natively rather than requiring imputation.
• Tabular datasets are usually small by deep learning standards, and networks need data.

You can already fit these models and read their scores. This section is about what each one does internally: the assumptions a linear model makes, how a tree picks a split, why fitting new trees to the errors of previous ones works at all, and where support vector machines and nearest neighbours still fit. The point is to choose one for a stated reason rather than by habit, and to know which of them you can read an explanation out of, where the answer for linear coefficients is often no.

Model                  Fitted by              Capacity set by        Read an explanation out of it?
linear / logistic      closed form or convex  the features you give  yes, if features are uncorrelated
decision tree          greedy search          max_depth, min_leaf    yes, as rules, if shallow
random forest          many trees, averaged   max_features           only as feature importances
gradient boosting      many trees, sequential learning_rate x trees  only as SHAP or importances
SVM                    convex optimisation    C and the kernel       not really
k-nearest neighbours   no fitting at all      k                      yes: here are the neighbours

Linear models#

Linear regression finds the weights minimising squared error. It assumes the relationship is linear and the errors are roughly independent with constant variance. When those assumptions hold it is fast, stable and interpretable, and each coefficient says how the prediction moves with that feature.

Squared error is not an arbitrary choice. It makes the problem convex with a single solution reachable in closed form, and it is what you get if you assume the noise is Gaussian and ask for the most likely weights. The same derivation with a different noise assumption gives a different loss, which is where absolute error and the rest come from.

Regularisation is the fix for that instability and for overfitting generally: add a penalty on the size of the weights, so the model prefers a simpler explanation unless the data insists otherwise. The two penalties behave differently in a way worth seeing rather than being told. L2, called ridge, penalises squared weights and shrinks everything smoothly towards zero without reaching it. L1, called lasso, penalises absolute weights, and its solution sets coefficients exactly to zero, which makes it a feature selector.

import numpy as np
from sklearn.linear_model import Ridge, Lasso

rng = np.random.default_rng(0)
X = rng.normal(0, 1, (200, 30))
true = np.zeros(30); true[:4] = [3.0, -2.0, 1.5, 1.0]     # only 4 of the 30 matter
y = X @ true + rng.normal(0, 0.5, 200)

for name, m in [("plain", Ridge(alpha=1e-9)),
                ("ridge a=10", Ridge(alpha=10.0)),
                ("lasso a=0.1", Lasso(alpha=0.1))]:
    m.fit(X, y)
    print(name, (np.abs(m.coef_) > 1e-6).sum(), np.abs(m.coef_[4:]).max().round(3))

# plain         nonzero 30/30    largest junk coefficient 0.101
# ridge a=10    nonzero 30/30    largest junk coefficient 0.091
# lasso a=0.1   nonzero  4/30    largest junk coefficient 0.000   <- exactly zero

So the choice is not about which regularises better. Use L2 when you believe most features contribute a little, which is the usual case and the safer default. Use L1 when you believe most do not and you want the model to tell you which, accepting that with correlated features it picks one of each group arbitrarily. ElasticNet mixes both and exists mainly to soften that arbitrariness.

How a tree splits#

At each node the tree tries every feature and every threshold and picks the split that makes the two resulting groups purest. For classification purity is usually measured by Gini impurity or entropy; for regression, by reduction in variance. It repeats until a stopping rule.

Node: 100 samples, 50 positive  (impurity high)

Split on age < 30:
  left   40 samples, 35 positive   (purer)
  right  60 samples, 15 positive   (purer)

Chosen because it reduced impurity more than any other split.

Gini and entropy almost always choose the same split, and the difference is not worth an argument: Gini is the default because it avoids a logarithm. The consequence worth knowing is a bias in the search itself. Trying every threshold means a feature with a thousand distinct values gets a thousand chances to look good by luck, while a binary feature gets one. Tree feature importances therefore favour high-cardinality columns, which is why an id column left in the data will appear near the top of the importance list having taught the model nothing.

Why boosting works#

Each new tree is fitted to the errors the ensemble currently makes, specifically to the gradient of the loss. The ensemble moves step by step towards lower loss, which is gradient descent in the space of functions. A small learning rate and many trees beats a large rate and few, at the cost of time.

That explains why boosting and bagging fix opposite problems, which is the thing to be able to say. A deep tree has low bias and high variance: it can represent anything and is wildly sensitive to which rows it saw. Averaging hundreds of them, as a forest does, leaves the bias alone and cancels the variance. A stump has the reverse problem, high bias and low variance, and boosting drives that bias down one correction at a time. Both arrive somewhere good from opposite directions.

from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor, HistGradientBoostingRegressor
from sklearn.model_selection import train_test_split

X2 = rng.uniform(-3, 3, (600, 1))
y2 = np.sin(X2[:, 0] * 2) + 0.3 * X2[:, 0] ** 2 + rng.normal(0, 0.2, 600)
Xtr, Xte, ytr, yte = train_test_split(X2, y2, test_size=0.3, random_state=0)

for name, m in [
    ("stump (depth 1)",      DecisionTreeRegressor(max_depth=1)),
    ("deep tree",            DecisionTreeRegressor()),
    ("forest of deep trees", RandomForestRegressor(n_estimators=300, random_state=0)),
    ("boosted stumps",       HistGradientBoostingRegressor(
                                 max_depth=1, max_iter=300,
                                 early_stopping=False, random_state=0)),
]:
    m.fit(Xtr, ytr)
    print(name, round(np.mean((ytr - m.predict(Xtr)) ** 2), 4),
                round(np.mean((yte - m.predict(Xte)) ** 2), 4))

# stump (depth 1)        train 0.5619  test 0.5512   high bias: wrong on train too
# deep tree              train 0.0000  test 0.0810   high variance: perfect on train, not on test
# forest of deep trees   train 0.0089  test 0.0588   averaging killed the variance
# boosted stumps         train 0.0507  test 0.0674   correcting killed the bias

Support vector machines#

An SVM draws the boundary that sits as far as possible from the nearest point of either class, and that widest-gap criterion is what distinguishes it from logistic regression, which will accept any boundary that separates the classes. Only the points on the edge of the gap, the support vectors, affect the answer; move a point deep inside its own class and the boundary does not shift.

The kernel trick is what made them famous. Instead of building the non-linear features explicitly, an SVM only ever needs the dot product between pairs of points, so you can substitute a function that returns what the dot product would have been in a much larger space without ever going there. An RBF kernel corresponds to an infinite-dimensional space, computed in one line. C controls how much misclassification is tolerated for a wider margin, gamma how local the RBF kernel's influence is, and getting both wrong is the usual reason an SVM disappoints.

The reason they are no longer a first choice is cost: fitting is between quadratic and cubic in the number of rows, so past roughly fifty thousand examples they stop being practical while boosting barely notices. They remain a sensible choice on a few hundred rows with many features, such as a small text classification problem, where trees struggle and there is not enough data for anything deeper.

K-nearest neighbours#

Store the training data. To predict, find the k closest examples and take their majority vote or their mean. There is no fitting step at all, which makes it the fastest model to train and the slowest to use: every prediction is a search over the whole training set, so the cost moves from a one-off to every request you ever serve.

It fails in high dimensions, and not gently. As dimensions grow, the distances between random points converge on each other, so the nearest neighbour stops being meaningfully nearer than the farthest and the vote becomes noise. This is the concrete form of the curse of dimensionality, and it is worth measuring once rather than believing on authority.

from scipy.spatial.distance import pdist      # scipy ships with scikit-learn

for d in (2, 10, 100, 1000):
    P = np.random.default_rng(0).normal(0, 1, (500, d))
    dists = pdist(P)
    print(d, round(dists.min() / dists.max(), 3))

# d=2      nearest / farthest = 0.001    neighbours mean something
# d=10     nearest / farthest = 0.108
# d=100    nearest / farthest = 0.504
# d=1000   nearest / farthest = 0.833    the nearest point is nearly as far as the farthest

Two things follow. Distances must be scaled, or the feature measured in metres silently outvotes the one measured in kilometres. And kNN is the conceptual ancestor of vector search: an embedding model's job is precisely to produce a space where nearness means similarity, so that this algorithm starts working again in a thousand dimensions.

Most machine learning at most companies is still a table of rows and a column you want to predict. This section works one through end to end on a real public dataset: loading it and looking at it properly, getting a crude baseline in place before anything clever, building the preprocessing into a pipeline so it cannot leak, comparing a linear model against gradient boosting, then examining the errors rather than only the score. It is the workflow every scikit-learn project reuses.

Step                     Why it is in this order                      Skipping it costs you
1 look at the data       you cannot choose a metric without it        a metric nobody wanted
2 hold out a test set    everything after this can overfit            a number you cannot trust
3 dummy baseline         gives every later number a meaning           shipping 0.87 that was free
4 pipeline               preprocessing refits per fold                a score that will not reproduce
5 linear model           tells you whether the problem is hard        no reference point
6 boosting               the actual candidate                         nothing, if 5 was good enough
7 test once              the only unbiased estimate you get           a slowly rotting test set
8 error analysis         says what to do next                         a model you cannot improve

The data#

from sklearn.datasets import fetch_openml

adult = fetch_openml("adult", version=2, as_frame=True)
X, y = adult.data, (adult.target == ">50K").astype(int)
print(X.shape, round(y.mean(), 4))            # (48842, 14)  0.2393
print((X.dtypes == "category").sum())         # 8 categorical, 6 numeric
print(X.isna().sum()[X.isna().sum() > 0].to_dict())
# {'workclass': 2799, 'occupation': 2809, 'native-country': 857}

Three facts from those four lines decide the rest of the project. The positive class is 24 per cent, so accuracy will be misleading and the metric should be ROC AUC or average precision rather than accuracy. Eight of the fourteen columns are categorical, which means either one-hot encoding for a linear model or a booster that accepts categories directly. And three columns have missing values, which is not a data error but an artefact of how the census was collected; an imputer plus an indicator will do, and boosting can take them as they are.

Split once, and hold out a test set#

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=0)

This happens before anything is inspected in detail, and stratify keeps the 24 per cent positive rate identical in both halves. The test set is now closed. Every decision from here is made on cross-validation over the training set.

The baseline#

from sklearn.dummy import DummyClassifier
from sklearn.model_selection import cross_val_score

print(cross_val_score(DummyClassifier(), X_train, y_train, scoring="roc_auc").mean())   # 0.5

A dummy scores exactly 0.5 ROC AUC by construction, which makes it look like a pointless step. It is not: the same dummy scores 0.76 on accuracy, and that is the number a stakeholder will hear first. Knowing that 0.76 is free is what stops anyone celebrating a model that scores 0.80.

A pipeline#

from sklearn.compose import make_column_transformer, make_column_selector
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.impute import SimpleImputer

prep = make_column_transformer(
    (make_pipeline(SimpleImputer(strategy="median"), StandardScaler()),
     make_column_selector(dtype_include="number")),
    (make_pipeline(SimpleImputer(strategy="most_frequent"), OneHotEncoder(handle_unknown="ignore")),
     make_column_selector(dtype_include="category")),
)
logreg = make_pipeline(prep, LogisticRegression(max_iter=2000))
print(cross_val_score(logreg, X_train, y_train, cv=5, scoring="roc_auc").mean())   # 0.9055

Everything that learns from data, the imputer, the scaler and the encoder, lives inside the pipeline, so cross-validation fits it on each training fold only. That one habit prevents the most common leakage in tabular work.

It is worth seeing how large that effect can be, because described in words it sounds like pedantry. Here is a dataset of pure noise with labels assigned at random, where the honest score is 50 per cent and nothing can beat it.

import numpy as np
from sklearn.feature_selection import SelectKBest, f_classif

rng = np.random.default_rng(0)
Xn = rng.normal(0, 1, (200, 5000))     # pure noise
yn = rng.integers(0, 2, 200)           # labels unrelated to it

# WRONG: pick the best 20 columns using every row, then cross-validate
best = SelectKBest(f_classif, k=20).fit(Xn, yn).transform(Xn)
print(cross_val_score(LogisticRegression(max_iter=2000), best, yn, cv=5).mean())   # 0.795

# RIGHT: selection inside the pipeline, refitted on each training fold
pipe = make_pipeline(SelectKBest(f_classif, k=20), LogisticRegression(max_iter=2000))
print(cross_val_score(pipe, Xn, yn, cv=5).mean())                                  # 0.510

With 5,000 noise columns and 200 rows, some columns correlate with the labels by luck. Choosing them while looking at all the rows lets that luck into every fold, and the model then re-finds it. The same mechanism operates more quietly with a scaler or an imputer, where it inflates a score by a point or two rather than thirty, and is correspondingly harder to notice.

Gradient boosting#

from sklearn.ensemble import HistGradientBoostingClassifier

hgb = HistGradientBoostingClassifier(categorical_features="from_dtype", random_state=0)
print(cross_val_score(hgb, X_train, y_train, cv=5, scoring="roc_auc").mean())   # 0.9282

That gap, a little over two points of AUC, is the honest size of the improvement you should expect from boosting on a table like this one. It is real and worth having, and it is not the order-of-magnitude difference the reputation suggests. If the linear model had come within half a point, the right decision would have been to ship the linear model.

Test once, then look at the errors#

from sklearn.metrics import roc_auc_score, classification_report

hgb.fit(X_train, y_train)
proba = hgb.predict_proba(X_test)[:, 1]
print(roc_auc_score(y_test, proba))                                # 0.9283
print(classification_report(y_test, (proba > 0.5).astype(int), digits=3))

#               precision    recall  f1-score   support
#            0      0.895     0.943     0.918      7431
#            1      0.781     0.649     0.709      2338
#     accuracy                          0.873      9769

The test AUC of 0.9283 against a cross-validated 0.9282 is the result you want: cross-validation was not optimistic, so the process was sound. The classification report is where the useful detail is. Recall on the positive class is 0.649, meaning the model misses roughly a third of the high earners, while its accuracy reads 0.873. If the application were finding people to contact, that 0.649 is the number that matters and it would be improved by lowering the threshold from 0.5, not by improving the model.

errors = X_test.assign(y=y_test, p=proba).query("(p > 0.5) != (y == 1)")

total = X_test.groupby("education", observed=True).size()
wrong = errors.groupby("education", observed=True).size()
print(((wrong / total).dropna().sort_values(ascending=False) * 100).round(1).head(6))

# Doctorate     18.8      overall error rate: 12.7
# Bachelors     17.7
# Masters       17.7
# Assoc-voc     16.6
# Prof-school   14.4
# Assoc-acdm    14.0

The model is half again as wrong on doctorates as it is overall. That is the finding, and it is actionable in a way the single number 0.9283 is not: highly educated people are the group where education alone stops predicting income, so the model needs a feature that separates them, or the documentation needs to say the predictions are less reliable there. Reporting the AUC and stopping would have hidden this entirely.

The test set is used once, at the end. Then read the errors grouped by slice: the groups where the model fails most are where more data, better features, or a note in the model's documentation are needed.

Save it#

import joblib
joblib.dump(hgb, "income_model.joblib")     # load with the same library versions

Bias and variance, regularisation, missing data, why high-dimensional space behaves strangely: a small set of questions about classical machine learning comes up again and again, and the usual answers are memorised phrases that collapse under one follow-up. This section works through each of them properly, so what you leave with is an explanation you could defend rather than a definition you can recite.

How to use this section#

These questions come up in almost every ML interview loop. Turn on 'Test yourself' below, answer each out loud before revealing, and practise the follow-ups. The follow-ups are where interviews are actually decided.

First, though, the five ideas the questions are drawn from, each worked rather than defined. A memorised definition of the bias-variance trade-off survives one question. Having watched the two columns of numbers behave differently survives three.

Bias and variance, in two columns#

Bias is error from a model that cannot represent the pattern. Variance is error from a model that is too sensitive to the particular rows it was shown. You do not need a decomposition to tell them apart: fit the model, print the training error next to the validation error, and read them.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures

rng = np.random.default_rng(0)
X = rng.uniform(-1, 1, (60, 1))
y = np.sin(3 * X[:, 0]) + rng.normal(0, 0.15, 60)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.4, random_state=0)

for deg in (1, 3, 9, 20):
    m = make_pipeline(PolynomialFeatures(deg), LinearRegression()).fit(Xtr, ytr)
    print(deg, round(np.mean((ytr - m.predict(Xtr)) ** 2), 4),
               round(np.mean((yte - m.predict(Xte)) ** 2), 4))

# degree  1   train 0.1694   test  0.3032    both bad        -> bias
# degree  3   train 0.0224   test  0.0317    close together  -> about right
# degree  9   train 0.0118   test  0.0485    gap opening     -> variance
# degree 20   train 0.0088   test 32.9775    best train, catastrophic test

training error   validation error   diagnosis   what actually helps
high             high, similar      bias        more features, more capacity, less regularisation
low              much higher        variance    more data, more regularisation, less capacity
low              low                neither     you are done; stop tuning
high             lower than train   a bug       leakage, a broken split, or mismatched sets

The follow-up that catches people is whether more data helps. It helps variance and does nothing for bias, and the difference is stark enough to measure in four lines.

for deg in (1, 15):
    for n in (30, 100, 1000, 10000):
        Xb = rng.uniform(-1, 1, (n, 1))
        yb = np.sin(3 * Xb[:, 0]) + rng.normal(0, 0.15, n)
        m = make_pipeline(PolynomialFeatures(deg), LinearRegression()).fit(Xb, yb)
        Xv = rng.uniform(-1, 1, (3000, 1))
        yv = np.sin(3 * Xv[:, 0]) + rng.normal(0, 0.15, 3000)
        print(deg, n, round(np.mean((yv - m.predict(Xv)) ** 2), 4))

# degree 1  (bias)      30: 0.1876   100: 0.1927   1000: 0.1877   10000: 0.1908
# degree 15 (variance)  30: 2559.14  100: 0.0307   1000: 0.0225   10000: 0.0223

Regularisation#

Regularisation adds a cost for complexity, so the model prefers a simpler explanation unless the data pays for the extra. L2 penalises squared weights and shrinks them smoothly; L1 penalises absolute weights and sets many exactly to zero, which makes it a feature selector. The strength is a hyperparameter and it is chosen by cross-validation, never by inspection.

The useful observation for an interview is that almost everything people do to stop overfitting is the same idea wearing different clothes. Early stopping limits how far the weights can travel from their small initial values, which is a budget on complexity. Dropout forces the network not to depend on any one unit. Data augmentation says the label should not change under transformations you specify. Max depth on a tree, and the number of boosting rounds, are the same lever again. If asked to name regularisation techniques, naming the mechanism they share is what separates a strong answer.

Choosing a model family#

This question is testing whether you have constraints in mind other than accuracy. Accuracy is where beginners stop, and it is rarely the binding constraint in a real deployment.

Constraint that actually binds        Pushes you towards
someone must sign off on decisions    logistic regression, a shallow tree
predictions in under a millisecond    linear model, or a small compiled tree ensemble
a few hundred rows                    regularised linear model
retraining daily on fresh data        anything that fits in minutes
features arrive with gaps             gradient boosting, which handles them natively
the answer must be a real probability logistic regression, or calibrate afterwards
accuracy above everything else        gradient boosting, tuned

Missing data#

The first question is never which imputation strategy. It is why the value is missing, because that determines whether imputing is honest. If a sensor fails at random, filling with the median loses almost nothing. If income is missing because the respondent declined to answer, and declining correlates with income, then the fact of missingness is one of the better features in the dataset and imputing it away destroys signal.

from sklearn.impute import SimpleImputer

imp = SimpleImputer(strategy="median", add_indicator=True)   # keeps the missingness

Two rules survive every variation of this question. Fit the imputer inside the cross-validation pipeline, or the median has been computed from rows the model is about to be tested on. And check whether the missingness pattern is the same in training and in production, because a column that is 2 per cent missing historically and 90 per cent missing in the live feed will not announce itself; the model will simply get quietly worse.

High dimensions#

Volume grows exponentially with dimensions, so data that felt dense in two dimensions is isolated points in two hundred. Covering the space to the same density needs exponentially more rows, which nobody has. The practical consequence is that distances stop discriminating: in a thousand dimensions the nearest of five hundred random points sits about 83 per cent as far away as the farthest, so any method built on neighbourhoods degrades towards a coin flip.

What follows is the answer to the follow-up about why anything works at all in high dimensions. Real data is not spread uniformly through its space; it lies near a much lower-dimensional surface within it, because the features are related. Feature selection, regularisation, PCA and embedding models are all ways of finding that surface and working on it instead. That is also the honest account of why vector search functions in 1,536 dimensions when the geometry above says it should not.