Tuning hyperparameters

What this chapter covers

Every model has two kinds of numbers inside it. The first kind it learns from the data — the weights of a linear fit, the split points of a tree, the support vectors of an SVM. The second kind you have to hand it before it can learn anything at all: how deep the tree may grow, how many neighbors to poll, how hard an SVM should push its boundary. Those are the hyperparameters, and the model can't discover them on its own, because they're the settings that decide what "learning" even means for it. This chapter is about choosing them honestly.

Honestly is the whole point, and it's why this chapter sits where it does. The previous two chapters built the tools it leans on: cross-validation, which gives you a stable score for a fixed model, and the bias-variance U, which is the shape every hyperparameter is really steering you along. Tuning is the loop that wraps around cross-validation — try a configuration, score it with k-fold, keep the winner — plus the discipline to not fool yourself while you do it. We build that loop from scratch in NumPy: a k-fold splitter reused straight from the cross-validation chapter, a cross-validated score for one configuration, then grid search that scores every combination on a grid, random search that samples instead, and the validation curve that shows one hyperparameter's whole story at once. Then we hand the same job to scikit-learn's GridSearchCV, RandomizedSearchCV, and validation_curve and check that our search selects the identical configuration at the identical score.

The vehicle is an RBF-kernel SVM on the Pima diabetes data — two hyperparameters, C and gamma, that trade off in a way you can see. The centerpiece is an animation of the search itself: watch grid search crawl across a two-parameter grid, filling in each cell's cross-validated score, while random search scatters the same budget of fits and finds a good spot in a fraction of the evaluations. The number that ends the chapter is the honest one: what tuning bought us on data the search never touched.

A bit of history

Grid search has no inventor, which tells you something. It's what everyone did because it's the obvious thing to do: list some values for each knob, try every combination, keep the best. It shows up in statistics and engineering long before machine learning had the name, under headings like "parameter sweep," and for decades it was less a method than a reflex. If you had two hyperparameters and an afternoon, you made a grid. The reflex survived into the scikit-learn era because it's trivial to write, trivial to parallelize, and easy to explain to a reviewer — every point on the grid got its fair, identical trial.

The reflex was also quietly wasteful, and the paper that said so out loud is James Bergstra and Yoshua Bengio's 2012 "Random Search for Hyper-Parameter Optimization" in the Journal of Machine Learning Research. Their argument is sharp and a little counterintuitive: for the same budget of trials, random search usually beats grid search, and the gap widens as you add hyperparameters. The reason is that most hyperparameters don't matter equally. If only two of your ten knobs actually move the score, a grid wastes almost all its trials re-testing the eight that don't — it evaluates the two that matter at only a handful of distinct values each, because the grid's resolution got spent on the useless dimensions. Random search, sampling all knobs independently every trial, tries a different value of the important ones on every single fit, so it explores the dimensions that matter far more finely for free. Grid search's neat lattice, it turns out, is exactly the wrong shape for a problem where you don't know in advance which axes count. That paper is why random search is the sensible default today, and why the field kept going — to Bayesian optimization, which we'll get to at the end, where the search learns from its own trials instead of sampling blind.

The intuition

Here's the surface we're searching. Two hyperparameters, C and gamma, laid out on a grid — five values each, so twenty-five combinations — and every cell colored by its five-fold cross-validated accuracy. This is the real thing our code computes, not a sketch. Read it before we build the search, because the search is just a way of discovering this picture one cell at a time without being able to see it first.

The two knobs do different jobs, and the surface shows both. Gamma sets how far a single training point's influence reaches: too small and every point smears across the whole space so the boundary can't bend — the dark, flat underfit band at the bottom, stuck near 0.65. Too large and each point's influence collapses to a tiny bubble around itself, so the model memorizes the training set and generalizes to nothing — the darkening at the top. C is the softness of the margin: how much the SVM is willing to misclassify a training point to keep a simpler boundary. The bright ridge running through the middle is where both are in balance, and the single brightest cell — C is 10, gamma is 0.01 — is what the search is hunting. You can see it here in one glance. The search can't; it has to earn it.

The math

A configuration is a choice of values for every hyperparameter — call it $\lambda$. For our SVM, $\lambda = (C, \gamma)$. The set of configurations the search considers is $\Lambda$: for grid search it's the Cartesian product of the per-knob value lists; for random search it's whatever gets sampled. The whole task is a single optimization over that set.

The objective is the cross-validated score from the previous chapter, now read as a function of the configuration. Split the training rows into k folds. Let $\hat{f}_{\lambda}^{(-j)}$ be the model trained with settings $\lambda$ on every fold except fold j, and let $F_j$ be the held-out rows of fold j. The cross-validation score of a configuration is the mean accuracy of those held-out predictions across the k folds:

Read the inner sum as "accuracy on fold j" and the outer as "average over folds" — it's exactly the $\mathrm{CV}_k$ from the cross-validation chapter, with the model's settings now named explicitly because they're the thing we're varying. Tuning is then one line. Pick the configuration on the grid that maximizes it:

That's the entire objective. Grid search evaluates $\mathrm{CV}(\lambda)$ at every $\lambda$ in the lattice and returns the argmax by brute force. Random search draws $\lambda_1, \ldots, \lambda_n$ independently from a distribution over the space and returns the best of those:

Same objective, different $\Lambda$. The only real modeling choice buried in here is the distribution $p(\lambda)$, and for scale parameters like C and gamma it should be log-uniform — sample the exponent uniformly — because what matters is the order of magnitude, not the raw value. A search that draws gamma uniformly on $[0.0001, 1]$ spends 90% of its draws above 0.1 and never really tries the small values where the sweet spot lives.

What it's good at, what it isn't

The strength of both searches is that they need nothing from you but honesty and patience. There's no gradient, no assumption about the model, no requirement that the score be differentiable or even continuous — the objective is a black box that eats a configuration and returns a number, and the search only has to compare numbers. That's why the exact same loop tunes an SVM, a random forest, and a neural network without changing a line: it never looks inside the model. Grid search adds a second virtue, reproducibility. The lattice is fixed, so the search is deterministic and every reviewer who runs it gets the same winner — which is why, for two or three hyperparameters over a modest range, it's still a perfectly good default and the one I reach for when I want the result to be boring and auditable.

The weakness is cost, and it's exponential in the worst way. Every hyperparameter you add multiplies the grid: five values of one knob is five fits, but five values each of four knobs is 625, times k for the cross-validation, and that's before the model itself is slow. This is the curse of dimensionality wearing a different hat, and it's exactly the regime where Bergstra and Bengio's point bites — the grid spends its budget on combinations that vary knobs that don't matter. Random search sidesteps the multiplication by decoupling the budget from the dimensionality: you choose how many fits you can afford and it spends them, whether there are two hyperparameters or twenty. Neither search is clever, though. Both evaluate blind, learning nothing from the trials they've already run, which is the opening that Bayesian optimization walks through.

The data

Same patients as the cross-validation chapter: the Pima diabetes set, 768 people, eight measurements each — glucose, BMI, age, blood pressure and the rest — and a binary label, diabetes or not, positive in about 35% of them. The model this time is an RBF-kernel SVM, which cares about scale, so the features are standardized; and because standardization is a thing learned from data, it's refit inside every fold on that fold's training rows only, never on the held-out ones. That's a small honesty that's easy to skip and quietly poisons a score if you do.

The bigger honesty is the split. Before any tuning happens, I set aside a quarter of the data — 192 patients — as a test set and do not touch it again until the final number. All the searching, all the cross-validation, all the argmax, happens on the other 576 patients. This is the golden rule of tuning, and it's worth stating flatly: the moment you choose a hyperparameter by a score, that score has been used for fitting and is no longer an unbiased estimate of anything. Cross- validation on the training set tells you which configuration to pick; it does not tell you how well the picked configuration will do, because you picked it for scoring well. The only clean estimate of that comes from data that played no part in the choice. So the structure is nested: an inner cross-validation loop selects the configuration, and an outer held-out set — evaluated exactly once, at the end — reports the number you'd actually quote.

Build it, one function at a time

The search is a small loop around a scoring function, and the scoring function is cross-validation, so we start by borrowing it. First the metric every fold needs:

def accuracy(y_true, y_pred):
    """Fraction of predictions that match the truth — unitless, 0 to 1."""
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    return float((y_true == y_pred).mean())

Then the folds. Tuning reuses cross-validation wholesale, so it reuses cross- validation's splitter — this is the stratified k-fold from the previous chapter, dropped in unchanged so the search's scores line up with scikit-learn's exactly. Stratified because the label is imbalanced and we want each fold to carry the same class balance:

def stratified_kfold_indices(y, k):
    """The week-25 splitter: k folds, each keeping the whole set's class balance.

    Tuning reuses cross-validation wholesale, so it reuses cross-validation's
    folds. This is the same StratifiedKFold(shuffle=False) allocation built in
    the cross-validation chapter — split each class across the folds separately
    so a rare label can't clump — reproduced here so the tuning loop is
    self-contained and its scores line up with scikit-learn's to the digit.
    """
    y = np.asarray(y)
    classes, y_enc = np.unique(y, return_inverse=True)
    ordered = np.sort(y_enc)
    per_fold = np.array([np.bincount(ordered[i::k], minlength=len(classes))
                         for i in range(k)])
    fold_id = np.empty(len(y), dtype=int)
    for c in range(len(classes)):
        fold_id[y_enc == c] = np.repeat(np.arange(k), per_fold[:, c])
    return [(np.where(fold_id != f)[0], np.where(fold_id == f)[0]) for f in range(k)]

Now the objective: the cross-validated score of one fixed configuration. This is the number the search maximizes. The estimator arrives as a fit callable that returns a predict function, so this loop is a black box to the model — it fits on each fold's training rows, scores the held-out rows, and averages:

def cv_score(fit, X, y, splits):
    """Cross-validated accuracy for one fixed configuration — the score to beat.

    `fit(X_train, y_train)` trains a fresh model and returns a `predict`
    callable. We rotate through the folds, fit on each fold's training rows,
    score the held-out rows, and return the mean across folds. That mean is the
    objective the search maximizes; the model is a black box to this function.
    """
    X, y = np.asarray(X), np.asarray(y)
    fold_acc = np.empty(len(splits))
    for j, (train, val) in enumerate(splits):
        predict = fit(X[train], y[train])
        fold_acc[j] = accuracy(y[val], predict(X[val]))
    return float(fold_acc.mean())

With one configuration scorable, grid search is a loop over the Cartesian product of the value lists. It records every combination's score — that's what the heatmap draws — and returns the argmax:

def grid_search(make_fit, grid, X, y, splits):
    """Exhaustive search: CV-score every combination on the grid, keep the best.

    `grid` maps each hyperparameter name to a list of candidate values; the
    search is the Cartesian product of those lists. `make_fit(params)` returns a
    `fit` callable wired with those settings. Returns every combination's score
    (for the heatmap) and the argmax — the single best configuration found.
    """
    names = list(grid)
    results = []
    for combo in itertools.product(*(grid[n] for n in names)):
        params = dict(zip(names, combo))
        score = cv_score(make_fit(params), X, y, splits)
        results.append({"params": params, "score": score})
    best = max(results, key=lambda r: r["score"])
    return results, best

Random search has the same skeleton with the loop body swapped: instead of walking a lattice, it draws each configuration from a sampler and spends a fixed budget of n_iter fits. The samplers are functions of the RNG — a log-uniform draw for C and gamma — and seeding the RNG makes the whole run repeatable:

def random_search(make_fit, samplers, n_iter, X, y, splits, rng):
    """Sample n_iter configurations, CV-score each, keep the best.

    `samplers` maps each hyperparameter to a function of the RNG that draws one
    value — a log-uniform draw for C and gamma, say. Unlike the grid, this never
    enumerates a lattice: it spends its whole budget on n_iter independent draws,
    which is why in high dimensions it finds a good point in far fewer fits.
    Seed the RNG and the run is exactly repeatable.
    """
    results = []
    for _ in range(n_iter):
        params = {name: sample(rng) for name, sample in samplers.items()}
        score = cv_score(make_fit(params), X, y, splits)
        results.append({"params": params, "score": score})
    best = max(results, key=lambda r: r["score"])
    return results, best

Last, the validation curve. This isn't a search — it's a diagnostic. Fix every hyperparameter but one, sweep that one across a range, and record two scores at each value: accuracy on the training folds and accuracy on the held-out folds. The training curve tells you how flexible the model has become; the validation curve tells you whether that flexibility is helping or hurting; and the gap between them is the overfitting, drawn as a distance:

def validation_curve(make_fit, param_name, values, X, y, splits, fixed=None):
    """Train and validation CV accuracy as one hyperparameter sweeps a range.

    Hold every other setting fixed, vary `param_name` across `values`, and at
    each value record two numbers: the mean accuracy on the training folds and
    the mean on the held-out folds. The training curve keeps rising as the model
    grows more flexible; the validation curve rises, peaks, then falls. That peak
    is the sweet spot, and the gap between the two curves is the overfitting.
    """
    fixed = dict(fixed or {})
    X, y = np.asarray(X), np.asarray(y)
    train = np.empty(len(values))
    val = np.empty(len(values))
    for i, v in enumerate(values):
        fit = make_fit({**fixed, param_name: v})
        tr_acc = np.empty(len(splits))
        va_acc = np.empty(len(splits))
        for j, (trn, va) in enumerate(splits):
            predict = fit(X[trn], y[trn])
            tr_acc[j] = accuracy(y[trn], predict(X[trn]))
            va_acc[j] = accuracy(y[va], predict(X[va]))
        train[i] = tr_acc.mean()
        val[i] = va_acc.mean()
    return train, val

Six small pieces: a metric, a splitter, a scorer, and three ways to point the scorer at the hyperparameter space. Nothing here knows it's tuning an SVM.

Watch it work

This is the search running for real. Two panels, one budget of fits each. On the left, grid search crawls the twenty-five-cell C-by-gamma grid, filling in one cell's cross-validated accuracy per frame, coloring it by score, with the running best marked by a star. On the right, random search spends the same fits — twelve of them — as log-uniform draws scattered across the same space, the faint squares showing where the grid would have looked. The caption names the current cell's configuration and cross-validated accuracy, the grid's best so far, and random search's best so far. The metric is accuracy — the fraction of held-out patients classified correctly, unitless, from 0 to 1.

Press play and watch the two strategies differ in character. Grid search is methodical and blind: it has to visit the dark underfit cells at the bottom and the dark overfit cells at the top with the same diligence as the bright ones, because it can't know they're wasted until it's scored them. It marches through all twenty-five, and its star lands on C is 10, gamma is 0.01 at a cross-validated accuracy of 0.781. Random search, meanwhile, ignores the lattice entirely. By its fourth draw it's already found a configuration scoring 0.778 — within 0.004 of the grid's eventual best — and it finishes its twelve fits at 0.778, having spent less than half the compute to land essentially the same place. On two knobs the gap is modest; the whole force of Bergstra and Bengio is that on ten knobs it's a chasm.

Now the validation curve, which is the other way to look at the same surface — one slice through it, in detail. Fix C at its best value and sweep gamma alone across a fine log range. The cyan line is the held-out cross-validation accuracy; the orange is the accuracy on the training folds; the marked point is the sweet spot where the validation curve peaks. Read it left to right and you're reading the bias-variance U from the last chapter, now indexed by a hyperparameter instead of a polynomial degree.

The two curves tell the whole story of gamma. At the far left, small gamma, both curves sit low together near 0.65 — the model is too stiff to bend, high bias, bad on training and validation alike. As gamma grows the validation curve climbs to its peak, cross-validated accuracy 0.781 around gamma 0.003, and this is the sweet spot the search is looking for. Then the curves split and it turns ugly fast: the training curve keeps rising toward 1.0, a perfect fit to the data it can see, while the validation curve turns and falls. By the right edge the training accuracy is exactly 1.0 and the validation accuracy has collapsed all the way back to 0.65 — the same place it started, which is just the majority-class baseline you'd get by predicting "no diabetes" for everyone. The model has memorized its training points and learned nothing that transfers. That widening gap is overfitting made visible, and the peak of the cyan line is exactly the value of gamma you'd want the search to return.

The full implementation

The whole file, top to bottom — the metric, the reused splitter, the scorer, and the three searches. No scikit-learn anywhere in it; the model is handed in from outside. This is the code both panels of the animation and the validation curve actually ran:

"""Hyperparameter tuning, built from scratch in pure NumPy.

The problem: a model has settings it cannot learn from the data on its own — an
SVM's C and gamma, a tree's depth, a neighbor count. You have to choose them,
and you want to choose them honestly: by how well the model does on data it did
not train on. Cross-validation (week 25) is the scoring engine that answers
"how good is *this* configuration"; tuning is the loop around it that asks the
question for many configurations and keeps the winner.

This file is only the search machinery — the k-fold splitter, the CV score for
one configuration, exhaustive grid search, random search, and the validation
curve. The estimator is handed in as a `fit` callable that returns a `predict`
function, so nothing here imports scikit-learn; the same loop tunes any model.
Every function appears in the chapter one step at a time (the `# region:`
markers are what the book's include directives pull in).
"""

import itertools

import numpy as np
import pandas as pd


# region: accuracy
def accuracy(y_true, y_pred):
    """Fraction of predictions that match the truth — unitless, 0 to 1."""
    y_true = np.asarray(y_true)
    y_pred = np.asarray(y_pred)
    return float((y_true == y_pred).mean())
# endregion


# region: stratified_kfold_indices
def stratified_kfold_indices(y, k):
    """The week-25 splitter: k folds, each keeping the whole set's class balance.

    Tuning reuses cross-validation wholesale, so it reuses cross-validation's
    folds. This is the same StratifiedKFold(shuffle=False) allocation built in
    the cross-validation chapter — split each class across the folds separately
    so a rare label can't clump — reproduced here so the tuning loop is
    self-contained and its scores line up with scikit-learn's to the digit.
    """
    y = np.asarray(y)
    classes, y_enc = np.unique(y, return_inverse=True)
    ordered = np.sort(y_enc)
    per_fold = np.array([np.bincount(ordered[i::k], minlength=len(classes))
                         for i in range(k)])
    fold_id = np.empty(len(y), dtype=int)
    for c in range(len(classes)):
        fold_id[y_enc == c] = np.repeat(np.arange(k), per_fold[:, c])
    return [(np.where(fold_id != f)[0], np.where(fold_id == f)[0]) for f in range(k)]
# endregion


# region: cv_score
def cv_score(fit, X, y, splits):
    """Cross-validated accuracy for one fixed configuration — the score to beat.

    `fit(X_train, y_train)` trains a fresh model and returns a `predict`
    callable. We rotate through the folds, fit on each fold's training rows,
    score the held-out rows, and return the mean across folds. That mean is the
    objective the search maximizes; the model is a black box to this function.
    """
    X, y = np.asarray(X), np.asarray(y)
    fold_acc = np.empty(len(splits))
    for j, (train, val) in enumerate(splits):
        predict = fit(X[train], y[train])
        fold_acc[j] = accuracy(y[val], predict(X[val]))
    return float(fold_acc.mean())
# endregion


# region: grid_search
def grid_search(make_fit, grid, X, y, splits):
    """Exhaustive search: CV-score every combination on the grid, keep the best.

    `grid` maps each hyperparameter name to a list of candidate values; the
    search is the Cartesian product of those lists. `make_fit(params)` returns a
    `fit` callable wired with those settings. Returns every combination's score
    (for the heatmap) and the argmax — the single best configuration found.
    """
    names = list(grid)
    results = []
    for combo in itertools.product(*(grid[n] for n in names)):
        params = dict(zip(names, combo))
        score = cv_score(make_fit(params), X, y, splits)
        results.append({"params": params, "score": score})
    best = max(results, key=lambda r: r["score"])
    return results, best
# endregion


# region: random_search
def random_search(make_fit, samplers, n_iter, X, y, splits, rng):
    """Sample n_iter configurations, CV-score each, keep the best.

    `samplers` maps each hyperparameter to a function of the RNG that draws one
    value — a log-uniform draw for C and gamma, say. Unlike the grid, this never
    enumerates a lattice: it spends its whole budget on n_iter independent draws,
    which is why in high dimensions it finds a good point in far fewer fits.
    Seed the RNG and the run is exactly repeatable.
    """
    results = []
    for _ in range(n_iter):
        params = {name: sample(rng) for name, sample in samplers.items()}
        score = cv_score(make_fit(params), X, y, splits)
        results.append({"params": params, "score": score})
    best = max(results, key=lambda r: r["score"])
    return results, best
# endregion


# region: validation_curve
def validation_curve(make_fit, param_name, values, X, y, splits, fixed=None):
    """Train and validation CV accuracy as one hyperparameter sweeps a range.

    Hold every other setting fixed, vary `param_name` across `values`, and at
    each value record two numbers: the mean accuracy on the training folds and
    the mean on the held-out folds. The training curve keeps rising as the model
    grows more flexible; the validation curve rises, peaks, then falls. That peak
    is the sweet spot, and the gap between the two curves is the overfitting.
    """
    fixed = dict(fixed or {})
    X, y = np.asarray(X), np.asarray(y)
    train = np.empty(len(values))
    val = np.empty(len(values))
    for i, v in enumerate(values):
        fit = make_fit({**fixed, param_name: v})
        tr_acc = np.empty(len(splits))
        va_acc = np.empty(len(splits))
        for j, (trn, va) in enumerate(splits):
            predict = fit(X[trn], y[trn])
            tr_acc[j] = accuracy(y[trn], predict(X[trn]))
            va_acc[j] = accuracy(y[va], predict(X[va]))
        train[i] = tr_acc.mean()
        val[i] = va_acc.mean()
    return train, val
# endregion


FEATURES = [
    "Pregnancies", "Glucose", "BloodPressure", "SkinThickness",
    "Insulin", "BMI", "DiabetesPedigreeFunction", "Age",
]


def load_xy(path="../data/diabetes.csv"):
    """Pima Indians Diabetes: 768 patients, 8 features, binary Outcome."""
    df = pd.read_csv(path)
    return df[FEATURES].to_numpy(float), df["Outcome"].to_numpy(int)

The library version

Nobody hand-rolls a tuning loop in production once they've written it, and scikit-learn's model_selection module has all three of ours as drop-in tools. The estimator is a Pipeline so the standardization is refit per fold exactly as we did by hand — the default C=1.0, gamma="scale" is the untuned baseline:

def make_pipeline(C=1.0, gamma="scale"):
    """Standardize-then-SVM as one estimator, so scaling is refit per fold.

    StandardScaler subtracts the mean and divides by the (population, ddof=0)
    standard deviation — the same transform impl.py's `fit` applies by hand.
    Wrapping it with the SVM in a Pipeline means cross-validation refits the
    scaler on each fold's training rows, so nothing about the held-out rows
    leaks into the fit. The default C=1.0, gamma="scale" is the untuned model.
    """
    return Pipeline([
        ("scale", StandardScaler()),
        ("svc", SVC(kernel="rbf", C=C, gamma=gamma)),
    ])

GridSearchCV is the counterpart of grid_search. Hand it the estimator, a grid as a dict, a scorer, and the same stratified folds, and it does the identical sweep, exposing every cell's mean cross-validation score in cv_results_ and the winner in best_params_:

def sklearn_grid(C_values, gamma_values, X, y, k):
    """GridSearchCV over the same C x gamma lattice, same stratified folds.

    Scores accuracy on every combination and reports the best. `cv_results_`
    holds the per-combination mean CV score we check our grid against; the
    integer `cv=k` on a classifier would already stratify, but we pass the
    splitter object so the folds are provably the ones impl.py builds.
    """
    grid = {"svc__C": C_values, "svc__gamma": gamma_values}
    search = GridSearchCV(
        make_pipeline(), grid, scoring="accuracy",
        cv=StratifiedKFold(n_splits=k, shuffle=False),
    )
    search.fit(X, y)
    return search

RandomizedSearchCV is the counterpart of random_search, and here scikit-learn's scipy.stats.loguniform does the log-scale sampling we wrote by hand. It uses a different RNG than ours, so its sampled points and its winner differ run to run — the point was never that the two random searches match, only that both spend a fixed budget and land near the grid's best:

def sklearn_random(C_lo, C_hi, g_lo, g_hi, n_iter, X, y, k, seed):
    """RandomizedSearchCV: n_iter log-uniform draws over the same box.

    Same method as our random_search — draw configurations from a log-uniform
    distribution over C and gamma and CV-score each — but a different RNG, so
    the sampled points and the winner differ run to run. The point isn't that
    the two random searches match; it's that both spend a fixed budget of fits
    and land near the grid's best.
    """
    dists = {"svc__C": loguniform(C_lo, C_hi),
             "svc__gamma": loguniform(g_lo, g_hi)}
    search = RandomizedSearchCV(
        make_pipeline(), dists, n_iter=n_iter, scoring="accuracy",
        cv=StratifiedKFold(n_splits=k, shuffle=False), random_state=seed,
    )
    search.fit(X, y)
    return search

And validation_curve is the counterpart of ours, returning the train and held-out scores for every value of a swept hyperparameter, one column per fold:

def sklearn_validation_curve(param_name, values, C, X, y, k):
    """validation_curve for one hyperparameter, other settings fixed.

    Returns (train_scores, test_scores), each shaped (len(values), k). Averaging
    the last axis gives the two curves impl.py's validation_curve returns.
    """
    est = make_pipeline(C=C)
    train, test = validation_curve(
        est, X, y, param_name=param_name, param_range=values,
        scoring="accuracy", cv=StratifiedKFold(n_splits=k, shuffle=False),
    )
    return np.asarray(train), np.asarray(test)

The one thing worth flagging is a convenience that hides the lesson: passing a bare integer, cv=5, to any of these on a classifier silently gives you stratified folds. That's a good default, but it means the stratification we were careful about is invisible unless you look for it. We pass the splitter object explicitly on both sides so the comparison is provably the same folds.

Scratch versus library

The claim is that our grid search and scikit-learn's are the same search, so they should select the same configuration at the same score — not close, identical. They are. Here's every one of the twenty-five cells, our cross-validated score against GridSearchCV's for the same cell; each point sits exactly on the diagonal because the underlying numbers match to the digit:

The generator asserts this equality on every run — all twenty-five cells and the argmax, plus the validation curve train and held-out scores, checked against scikit-learn to a maximum absolute difference of zero — so if the two ever drifted, the data wouldn't build. Matching the library isn't the achievement; it's the proof that the from-scratch search is the real thing and not a lookalike.

Now the number the whole chapter is for. We tuned on the 576 training patients and never touched the 192 test patients until this moment. Here's the default SVM against the grid-tuned one, scored two ways: the cross-validation score on the training set that the tuning optimized, and the held-out test accuracy that it didn't:

Tuning lifted the cross-validation score from 0.774 to 0.781 — a real gain of about seven tenths of a point, robust across the folds, and exactly what the search was built to maximize. On the held-out test set it moved from 0.734 to 0.776, a bigger jump. I want to be careful about that test number, though, because the test set is only 192 patients, so one patient is half a point and the test estimate has real scatter; the honest reading is that tuning bought a solid, trustworthy improvement in the cross-validated score, and the untouched test set agrees rather than disagrees. What you should not do is read the test bar as precise to the digit. It's one draw, which is exactly why you tune on cross-validation and only glance at the test set once, at the end, as a sanity check on the estimate you already trust.

Finally, grid against random, the Bergstra comparison in miniature. Both searched the same space; grid spent 25 fits, random spent 12:

Grid search's best cross-validated accuracy was 0.781 for 25 fits; random search's was 0.778 for 12. Random landed three thousandths of a point below the grid's best using less than half the evaluations — on this small two-knob problem, a modest trade. The reason to care is that the trade doesn't stay modest. Add hyperparameters and the grid's fit count explodes while random search's stays exactly whatever you budgeted, so the same three-thousandths-of-a-point gap comes at a tenth of the cost, then a hundredth. That's the whole argument for reaching for random search first, and it's why scikit-learn's own RandomizedSearchCV, run here with its different RNG, happened to land the grid's exact best cell at 0.781 in its twelve draws.

Takeaways

Start with random search, not grid. On one or two hyperparameters over a small range the two are interchangeable and grid's determinism is a mild plus, but that's the only case where grid wins, and it's rarer than it feels. The moment you have three or more knobs — and most real models do — the grid is spending most of its budget varying things that don't matter, and random search gets you a better result for the same compute, or the same result for a fraction of it. Fix your fit budget to what you can afford and let random search spend it; that's the default that scales.

Never tune on the test set, and mean it structurally, not just as a slogan. The score you optimize is contaminated the instant you optimize it — that's not a risk, it's a definition — so the configuration you pick by cross-validation needs a clean, untouched hold-out to be measured against, evaluated exactly once. Tune on the inner loop, report on the outer one. Every time I've seen a model look great in development and disappoint in production, the first place to check is whether the "held-out" number was quietly used to pick something, because a test set you looked at twice is a validation set wearing a disguise.

Search on a log scale for anything that's a scale — C, gamma, learning rate, regularization strength, the number of trees. What matters for these is the order of magnitude, and a linear grid or a uniform sampler wastes almost all its trials at the top of the range and never probes the small values where the answer usually lives. Sample the exponent, not the number.

And keep the whole enterprise in proportion, because this is the takeaway I'd fight for hardest. Tuning bought us seven tenths of a point of cross-validated accuracy here, which is real and worth having, but it's a polish, not a transformation. A better model, a better feature, more or cleaner data — any of those moves the needle more than perfect hyperparameters ever will, and they move it in a way that doesn't evaporate when the data shifts. The junior instinct is to grid-search a mediocre model into the ground; the senior move is to get a decent baseline working end to end, tune it once with random search to claim the easy gain, and spend the rest of your time on the data and the problem framing where the real accuracy hides. The next step past blind search is Bayesian optimization — tools like Optuna and Hyperopt that build a model of the score surface from the trials they've already run and sample where the payoff looks likeliest, turning tuning from exhaustive into adaptive. It's the right tool when each fit is expensive and the space is large. But it's the same objective we wrote in one line at the top of this chapter, just searched more cleverly, and it's still no substitute for a model that was worth tuning in the first place.