Hyperparameter Tuning and Model Selection

Hyperparameter tuning uses data to choose settings such as tree depth or regularization strength. Its winning validation score can be optimistic because the search selected that score from noisy estimates. We need both an efficient search and an evaluation kept outside the selection process. Throughout the examples, higher scores are better; the classifier searches use accuracy. Run the Python blocks in order in an environment with NumPy, SciPy, and scikit-learn installed.

How grid and random search spend a budget

A grid with \(v\) values for each of \(d\) hyperparameters evaluates \(v^d\) combinations, but only \(v\) distinct values per parameter. Two parameters with three values each already require nine combinations. If only one parameter strongly affects the score, many grid evaluations reuse the same value of that parameter. Random sampling can explore more distinct values at the same candidate budget. This example uses a synthetic score, not a model’s accuracy: \(x_0\) controls a narrow peak, while the other four coordinates each add a smaller effect. Its maximum is 1.2, so scores above 1 are valid.

import numpy as np

def f(x):
    return np.exp(-((x[0] - 0.37) ** 2) / 0.005) + 0.05 * np.sum(np.cos(6 * x[1:]))

n_dim = 5
for per_axis in (3, 4):
    budget = per_axis ** n_dim
    axes = np.linspace(0, 1, per_axis)
    grid = np.stack(np.meshgrid(*[axes] * n_dim, indexing="ij"), -1).reshape(-1, n_dim)
    grid_best = max(f(c) for c in grid)
    rnd = []
    for rep in range(40):
        g = np.random.default_rng(rep)
        rnd.append(max(f(g.random(n_dim)) for _ in range(budget)))
    print(f"{per_axis} values/axis  budget {budget:5d}   grid {grid_best:.4f}"
          f"   random {np.mean(rnd):.4f}   random wins {np.mean(np.array(rnd) > grid_best):.2f}")
    print(f"    distinct x0 values tried:  grid {per_axis}   random {budget}")
# 3 values/axis  budget   243   grid 0.2340   random 1.0683   random wins 1.00
#     distinct x0 values tried:  grid 3   random 243
# 4 values/axis  budget  1024   grid 0.9642   random 1.1230   random wins 1.00
#     distinct x0 values tried:  grid 4   random 1024

At 243 evaluations the grid scores 0.2340, while the mean of 40 random-search best scores is 1.0683. Random search wins all 40 trials here. It samples 243 distinct values of the dominant coordinate, whereas the grid reuses only three. The remaining coordinates still matter: together their cosine terms range from −0.2 to 0.2. They were not removed from the objective.

The four-value grid includes \(1/3\), closer to the peak at 0.37 than the three-value grid’s nearest point, 0.5. Its score rises to 0.9642. Within that 1,024-evaluation budget, random search averages 1.1230, a difference of 0.1588; the 243-evaluation difference is 0.8343. Both rows compare equal candidate counts. These two grids are not nested. Adding points while retaining all old grid points cannot lower the best score of a deterministic objective, although it increases the cost.

This illustrates the motivation in Bergstra and Bengio’s Random Search for Hyper-Parameter Optimization: random trials can use a budget efficiently when only some dimensions strongly affect performance. A fixed random sampler does not learn which dimensions matter. The search ranges, distribution of good configurations, interactions, and evaluation noise still affect the outcome. Small grids remain useful when there are a few deliberate settings to compare.

Spending the budget adaptively

A fixed-budget random search evaluates all sampled configurations at the chosen resource level. Successive halving allocates less resource at first, retains the strongest fraction, and reevaluates survivors with more. In this example 27 candidates become 9, then 3, then 1. A resource level is the nominal number of rows used across a CV split; each fit uses its training portion. We fix every forest at 40 trees to make the resource comparison easier to follow, and keep 800 test rows outside all three searches.

import time
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import RandomizedSearchCV, train_test_split
from sklearn.experimental import enable_halving_search_cv
from sklearn.model_selection import HalvingRandomSearchCV
from scipy.stats import loguniform, randint

X, y = make_classification(n_samples=3200, n_features=25, n_informative=10,
                           random_state=0)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=800, stratify=y, random_state=0)
space = {"max_depth": randint(2, 20), "max_features": loguniform(0.05, 1.0),
         "min_samples_leaf": randint(1, 40)}
rf = lambda: RandomForestClassifier(n_estimators=40, random_state=0)
searches = {
    "random": RandomizedSearchCV(rf(), space, n_iter=27, cv=3, random_state=0),
    "halving smallest": HalvingRandomSearchCV(
        rf(), space, n_candidates=27, cv=3, random_state=0, factor=3),
    "halving exhaust": HalvingRandomSearchCV(
        rf(), space, n_candidates=27, cv=3, random_state=0, factor=3,
        min_resources="exhaust"),
}
for name, search in searches.items():
    t = time.perf_counter()
    search.fit(X_train, y_train)
    elapsed = time.perf_counter() - t
    resources = getattr(search, "n_resources_", [len(X_train)])
    candidates = getattr(search, "n_candidates_", [27])
    fits = 3 * sum(candidates) + 1
    print(f"{name}: selection CV {search.best_score_:.4f}"
          f"  held-out accuracy {search.score(X_test, y_test):.4f}")
    print(f"  resources {resources}  candidates {candidates}  fits {fits}")
    print(f"{name} elapsed {elapsed:.2f} s")
# random: selection CV 0.8375  held-out accuracy 0.8600
#   resources [2400]  candidates [27]  fits 82
# random elapsed 12.25 s # varies by machine
# halving smallest: selection CV 0.7130  held-out accuracy 0.7662
#   resources [12, 36, 108, 324]  candidates [27, 9, 3, 1]  fits 121
# halving smallest elapsed 6.48 s # varies by machine
# halving exhaust: selection CV 0.8346  held-out accuracy 0.8600
#   resources [88, 264, 792, 2376]  candidates [27, 9, 3, 1]  fits 121
# halving exhaust elapsed 7.79 s # varies by machine

The output separates the score used for selection from the accuracy of the selected, refitted model on the common test set. Times include the final refit and vary with the machine and load; one timing run does not establish a stable speed ratio. Scores can also vary with library versions. The three fixed candidates here are search procedures, and using their test results to choose a procedure would require a new final test set.

With HalvingRandomSearchCV, the default min_resources="smallest" starts this three-fold binary example at \(3\times2\times2=12\) nominal rows. Training and validation folds are both subsampled: at the first rung each fit has eight training rows and four validation rows. Such small validation sets make rankings coarse and unstable. Inspect the printed resource and candidate lists instead of assuming the final rung reaches the full training size. After selection, the default refit=True still trains the winner on all 2,400 training rows.

The exhaust setting increases the starting resource so the last rung approaches the maximum permitted by the schedule. Integer rounding means it need not use every row. The final-rung CV scores can therefore use different sample sizes; the common held-out scores compare the final refits more directly. Here random search and the exhaust schedule both reach held-out accuracy 0.8600, while the smallest-resource schedule reaches 0.7662. Random search performs 82 forest fits and each halving schedule performs 121; many of the latter use far fewer rows. Fit count alone therefore does not measure total work. Matching one test accuracy neither guarantees the same winner nor establishes equivalence on other data.

Halving depends on low-resource rankings being informative about high-resource performance. Early noise, strong regularization, or a model that benefits disproportionately from more rows can eliminate a promising candidate. A larger first rung or a less aggressive elimination factor can reduce that risk at extra cost. The experiment measures one schedule and dataset, not how often rankings remain stable in general.

Bayesian optimization fits a surrogate, a model that predicts evaluation scores from hyperparameter settings. It then uses an acquisition rule to balance promising regions against uncertainty. It can be attractive when evaluations are expensive enough to justify this extra modeling and optimization work. Sequential decisions depend on earlier results, but batch acquisition methods can propose several evaluations together. Random search remains a useful baseline; there is no universal trial count at which one approach becomes preferable.

The winner’s score is not the winner’s performance

Taking the largest of 25 validation estimates can favor candidates whose scores contain positive estimation errors. The following no-signal experiment repeats this selection on 30 independently generated datasets per sample size. The SVC is a classifier with two settings: C controls inverse regularization and gamma the locality of its RBF kernel. Understanding that kernel is not needed for the comparison; the grid simply supplies 25 candidate configurations.

import numpy as np
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV, cross_val_score, StratifiedKFold, KFold

for n in (60, 120, 300):
    flat, nested, held = [], [], []
    for rep in range(30):
        g = np.random.default_rng(rep)
        X = g.normal(size=(n, 20))
        y = (g.random(n) < 0.5).astype(int)      # no signal at all: truth is 0.5
        grid = {"C": [0.01, 0.1, 1, 10, 100], "gamma": [1e-4, 1e-3, 1e-2, 1e-1, 1]}
        gs = GridSearchCV(SVC(), grid, cv=StratifiedKFold(3, shuffle=True, random_state=0))
        gs.fit(X, y)
        flat.append(gs.best_score_)              # the winner's own CV score
        nested.append(cross_val_score(
            gs, X, y, cv=KFold(3, shuffle=True, random_state=1)).mean())
        X_test = g.normal(size=(2000, 20))
        y_test = g.integers(0, 2, 2000)
        held.append(gs.score(X_test, y_test))
    mc_se = np.std(nested, ddof=1) / np.sqrt(len(nested))
    print(f"n={n:4d}  selected inner CV {np.mean(flat):.4f}"
          f"  nested {np.mean(nested):.4f}  nested MC SE {mc_se:.4f}"
          f"  independent test {np.mean(held):.4f}")
# n=  60  selected inner CV 0.5733  nested 0.4800  nested MC SE 0.0123  independent test 0.4973
# n= 120  selected inner CV 0.5572  nested 0.5039  nested MC SE 0.0091  independent test 0.5002
# n= 300  selected inner CV 0.5474  nested 0.5064  nested MC SE 0.0074  independent test 0.4993

The labels are independent fair coins, so any classifier fitted without access to a new label has expected accuracy 0.5 on that label. A measured score need not equal 0.5000. At 60 rows, the selected inner-CV mean is 0.5733, about 7.33 percentage points above 0.5. The nested mean is 0.4800 with Monte Carlo standard error 0.0123, while the separate test mean is 0.4973. The selected inner-CV mean is compared with both nested CV and a new 2,000-row test set for the final fitted winner. This separates the reference performance from finite-run estimates of it.

Nested CV repeats the complete search inside each outer training fold, then scores its selected model on the untouched outer fold. Here the outer KFold split is independent of the labels; stratification is used only inside each training fold. Each outer fit uses two thirds of the available rows. The output’s Monte Carlo standard error, nested MC SE, is the standard deviation of the 30 dataset-level nested estimates divided by \(\sqrt{30}\). It measures uncertainty in their mean, not the uncertainty of one CV run. A departure from 0.5 in a finite set of repetitions does not show that the inner search transmitted optimism to the outer scores.

For independent held-out accuracy, the outer test labels cannot influence their own model’s search or fit. Stratifying a split by labels conditions which rows are held out, so exact finite-sample expectation arguments require care; the label-independent outer split makes that distinction explicit here. Stratification remains useful for keeping class counts stable; this simulation uses ordinary outer folds to isolate the independence argument. In real problems, nested CV evaluates the tuning procedure at the outer training size. Preprocessing, feature selection, and search-space decisions based on data must also remain inside that outer training process; grouped or temporal data requires suitable inner and outer splits.

For a fixed fitted classifier and independent test rows from the target distribution, mean test accuracy is an unbiased estimate of that classifier’s accuracy. This statement does not automatically extend to every nonlinear metric or a shifted deployment distribution. A separate test set can evaluate the final fit without nested CV. Nested CV evaluates repeated executions of the search-and-fit procedure on smaller training subsets; a final search on all development data is still needed to produce the deployment model.

Exercises

1. How selection optimism grows. Give candidates the same true accuracy and vary how many independent noisy validation scores are searched.

Compare the increase with the square-root-of-logarithm approximation, and identify where its assumptions differ from an actual model search.

import numpy as np

n_val, true_acc = 500, 0.70
print(f"{'candidates':>11} {'reported best':>14} {'true':>7} {'optimism':>9}")
for k in (1, 2, 5, 10, 50, 200, 1000):
    best = []
    for rep in range(2000):
        g = np.random.default_rng(rep)
        scores = g.binomial(n_val, true_acc, k) / n_val
        best.append(scores.max())
    print(f"{k:11d} {np.mean(best):14.4f} {true_acc:7.4f} {np.mean(best)-true_acc:+9.4f}")
#  candidates  reported best    true  optimism
#           1         0.7002  0.7000   +0.0002
#           2         0.7115  0.7000   +0.0115
#           5         0.7237  0.7000   +0.0237
#          10         0.7317  0.7000   +0.0317
#          50         0.7461  0.7000   +0.0461
#         200         0.7556  0.7000   +0.0556
#        1000         0.7651  0.7000   +0.0651

2. Sampling on the right scale. Search a regularization strength with a uniform and a log-uniform distribution over the same range, and compare what each actually explores.

Count how often each sampler explores strong regularization, then compare the best observed CV scores without treating either sampled winner as the true optimum.

import numpy as np
from scipy.stats import loguniform, uniform
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import RandomizedSearchCV

X, y = make_classification(n_samples=1500, n_features=60, n_informative=8,
                           n_redundant=20, random_state=0)

for label, dist in (("uniform(1e-4, 100)", uniform(loc=1e-4, scale=100-1e-4)),
                    ("loguniform(1e-4, 100)", loguniform(1e-4, 100))):
    s = RandomizedSearchCV(LogisticRegression(max_iter=5000), {"C": dist},
                           n_iter=30, cv=5, random_state=0).fit(X, y)
    draws = np.asarray(s.cv_results_["param_C"], dtype=float)
    print(f"{label:24s} best C {s.best_params_['C']:9.4f}"
          f"  CV {s.best_score_:.4f}   draws below 1.0: {(draws < 1).sum():2d}/30")
# uniform(1e-4, 100)       best C   54.8814  CV 0.9080   draws below 1.0:  0/30
# loguniform(1e-4, 100)    best C    0.0200  CV 0.9113   draws below 1.0: 18/30

3. When is nested CV worth its cost. Compare a nested-CV estimate against a simple held-out test set at several dataset sizes, counting the fits each requires.

Distinguish fit counts from runtime, and distinguish one observed score gap from evidence about estimator variance.

import numpy as np, time
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import (GridSearchCV, cross_val_score,
                                     train_test_split, StratifiedKFold)

grid = {"max_depth": [3, 6, 12], "min_samples_leaf": [1, 5, 20]}
for n in (300, 1200, 5000):
    X, y = make_classification(n_samples=n, n_features=20, n_informative=8, random_state=0)
    mk = lambda: GridSearchCV(RandomForestClassifier(n_estimators=100, random_state=0),
                              grid, cv=StratifiedKFold(4, shuffle=True, random_state=0))

    t = time.perf_counter()
    nested = cross_val_score(mk(), X, y,
                             cv=StratifiedKFold(4, shuffle=True, random_state=1)).mean()
    t_nested = time.perf_counter() - t

    t = time.perf_counter()
    Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=1, stratify=y)
    holdout = mk().fit(Xtr, ytr).score(Xte, yte)
    t_hold = time.perf_counter() - t

    print(f"n={n:5d}  nested {nested:.4f} ({4*(9*4+1)} fits)"
          f"   holdout {holdout:.4f} ({9*4+1} fits)"
          f"   gap {nested - holdout:+.4f}")
    print(f"n={n} elapsed nested {t_nested:.2f} s, holdout {t_hold:.2f} s")
# n=  300  nested 0.8833 (148 fits)   holdout 0.8533 (37 fits)   gap +0.0300
# n=300 elapsed nested 18.72 s, holdout 4.61 s # varies by machine
# n= 1200  nested 0.9300 (148 fits)   holdout 0.9133 (37 fits)   gap +0.0167
# n=1200 elapsed nested 24.92 s, holdout 6.25 s # varies by machine
# n= 5000  nested 0.9474 (148 fits)   holdout 0.9392 (37 fits)   gap +0.0082
# n=5000 elapsed nested 59.07 s, holdout 14.89 s # varies by machine