Tuning Models

MLJ provides several built-in and third-party options for optimizing a model's hyper-parameters. The quick-reference table below omits some advanced keyword options.

Below we illustrate hyperparameter optimization using the Grid, RandomSearch, LatinHypercube and Explicit tuning strategies.

Overview

In MLJ model tuning is implemented as a model wrapper. After wrapping a model in a tuning strategy and binding the wrapped model to data in a machine called mach, calling fit!(mach) instigates a search for optimal model hyperparameters, within a specified range, and then uses all supplied data to train the best model. To predict using that model, one then calls predict(mach, Xnew). In this way, the wrapped model may be viewed as a "self-tuning" version of the unwrapped model. That is, wrapping the model simply transforms certain hyper-parameters into learned parameters.

A corollary of the tuning-as-wrapper approach is that the evaluation of the performance of a TunedModel instance using evaluate! implies nested resampling. This approach is inspired by MLR. See also below.

In MLJ, tuning is an iterative procedure, with an iteration parameter n, the total number of model instances to be evaluated. Accordingly, tuning can be controlled using MLJ's IteratedModel wrapper. After familiarizing oneself with the TunedModel wrapper described below, see Controlling model tuning for more on this advanced feature.

For a more in-depth overview of tuning in MLJ, or for implementation details, see the MLJTuning documentation. For a complete list of options see the TunedModel doc-string below.

Tuning a single hyperparameter using a grid search (regression example)

using MLJ
X = MLJ.table(rand(100, 10));
y = 2X.x1 - X.x2 + 0.05*rand(100);
Tree = @load DecisionTreeRegressor pkg=DecisionTree verbosity=0;
tree = Tree()

DecisionTreeRegressor(
  max_depth = -1, 
  min_samples_leaf = 5, 
  min_samples_split = 2, 
  min_purity_increase = 0.0, 
  n_subfeatures = 0, 
  post_prune = false, 
  merge_purity_threshold = 1.0, 
  feature_importance = :impurity, 
  rng = Random._GLOBAL_RNG())

Let's tune min_purity_increase in the model above, using a grid-search. To do so we will use the simplest range object, a one-dimensional range object constructed using the range method:

r = range(tree, :min_purity_increase, lower=0.001, upper=1.0, scale=:log);
self_tuning_tree = TunedModel(
    model=tree,
    resampling=CV(nfolds=3),
    tuning=Grid(resolution=10),
    range=r,
    measure=rms
);

DeterministicTunedModel(
  model = DecisionTreeRegressor(
        max_depth = -1, 
        min_samples_leaf = 5, 
        min_samples_split = 2, 
        min_purity_increase = 0.0, 
        n_subfeatures = 0, 
        post_prune = false, 
        merge_purity_threshold = 1.0, 
        feature_importance = :impurity, 
        rng = Random._GLOBAL_RNG()), 
  tuning = Grid(
        goal = nothing, 
        resolution = 10, 
        shuffle = true, 
        rng = Random._GLOBAL_RNG()), 
  resampling = CV(
        nfolds = 3, 
        shuffle = false, 
        rng = Random._GLOBAL_RNG()), 
  measure = RootMeanSquaredError(), 
  weights = nothing, 
  class_weights = nothing, 
  operation = nothing, 
  range = NumericRange(0.001 ≤ min_purity_increase ≤ 1.0; origin=0.5005, unit=0.4995; on log scale), 
  selection_heuristic = MLJTuning.NaiveSelection(nothing), 
  train_best = true, 
  repeats = 1, 
  n = nothing, 
  acceleration = CPU1{Nothing}(nothing), 
  acceleration_resampling = CPU1{Nothing}(nothing), 
  check_measure = true, 
  cache = true)

Incidentally, a grid is generated internally "over the range" by calling the iterator method with an appropriate resolution:

iterator(r, 5)

5-element Vector{Float64}:
 0.0010000000000000002
 0.005623413251903492
 0.0316227766016838
 0.1778279410038923
 1.0

Non-numeric hyperparameters are handled a little differently:

selector = FeatureSelector();
r2 = range(selector, :features, values = [[:x1,], [:x1, :x2]]);
iterator(r2)

2-element Vector{Vector{Symbol}}:
 [:x1]
 [:x1, :x2]

Unbounded ranges are also permitted. See the range and iterator docstrings below for details, and the sampler docstring for generating random samples from one-dimensional ranges (used internally by the RandomSearch strategy).

Returning to the wrapped tree model:

mach = machine(self_tuning_tree, X, y);
fit!(mach, verbosity=0)

trained Machine; does not cache data
  model: DeterministicTunedModel(model = DecisionTreeRegressor(max_depth = -1, …), …)
  args: 
    1:	Source @026 ⏎ Table{AbstractVector{Continuous}}
    2:	Source @744 ⏎ AbstractVector{Continuous}

We can inspect the detailed results of the grid search with report(mach) or just retrieve the optimal model, as here:

fitted_params(mach).best_model

DecisionTreeRegressor(
  max_depth = -1, 
  min_samples_leaf = 5, 
  min_samples_split = 2, 
  min_purity_increase = 0.0021544346900318843, 
  n_subfeatures = 0, 
  post_prune = false, 
  merge_purity_threshold = 1.0, 
  feature_importance = :impurity, 
  rng = Random._GLOBAL_RNG())

For more detailed information, we can look at report(mach), for example:

entry = report(mach).best_history_entry

(model = DecisionTreeRegressor(max_depth = -1, …),
 measure = StatisticalMeasuresBase.RobustMeasure{StatisticalMeasuresBase.FussyMeasure{StatisticalMeasuresBase.RobustMeasure{StatisticalMeasuresBase.Multimeasure{StatisticalMeasuresBase.SupportsMissingsMeasure{StatisticalMeasures.RootMeanSquaredErrorOnScalars}, Nothing, StatisticalMeasuresBase.RootMean{Int64}, typeof(identity)}}, Nothing}}[RootMeanSquaredError()],
 measurement = [0.24315885564675355],
 per_fold = [[0.23497325986087308, 0.1907463204704172, 0.2930882043027645]],)

Predicting on new input observations using the optimal model, trained on all the data bound to mach:

Xnew  = MLJ.table(rand(3, 10));
predict(mach, Xnew)

3-element Vector{Float64}:
  0.5876731958281213
 -0.11262620502317926
  1.5001919496386296

Or predicting on some subset of the observations bound to mach:

test = 1:3
predict(mach, rows=test)

3-element Vector{Float64}:
 -0.0436042437015236
  1.1495316705361969
  1.1495316705361969

For tuning using only a subset train of all observation indices, specify rows=train in the above fit! call. In that case, the above predict calls would be based on training the optimal model on all train rows.

A probabilistic classifier example

Tuning a classifier is not essentially different from tuning a regressor. A common gotcha however is to overlook the distinction between supervised models that make point predictions (subtypes of Deterministic) and those that make probabilistic predictions (subtypes of Probabilistic). The DecisionTreeRegressor model in the preceding illustration was deterministic, so this example will consider a probabilistic classifier:

info("KNNClassifier").prediction_type

:probabilistic

X, y = @load_iris
KNN = @load KNNClassifier verbosity=0
knn = KNN()

KNNClassifier(
  K = 5, 
  algorithm = :kdtree, 
  metric = Distances.Euclidean(0.0), 
  leafsize = 10, 
  reorder = true, 
  weights = NearestNeighborModels.Uniform())

We'll tune the hyperparameter K in the model above, using a grid-search once more:

K_range = range(knn, :K, lower=5, upper=20);

NumericRange(5 ≤ K ≤ 20; origin=12.5, unit=7.5)

Since the model is probabilistic, we can choose either: (i) a probabilistic measure, such as brier_loss; or (ii) use a deterministic measure, such as misclassification_rate (which means predict_mean is called instead of predict under the hood).

Case (i) - probabilistic measure:

self_tuning_knn = TunedModel(
    model=knn,
    resampling = CV(nfolds=4, rng=1234),
    tuning = Grid(resolution=5),
    range = K_range,
    measure=BrierLoss()
);

mach = machine(self_tuning_knn, X, y);
fit!(mach, verbosity=0);

trained Machine; does not cache data
  model: ProbabilisticTunedModel(model = KNNClassifier(K = 5, …), …)
  args: 
    1:	Source @021 ⏎ Table{AbstractVector{Continuous}}
    2:	Source @788 ⏎ AbstractVector{Multiclass{3}}

Case (ii) - deterministic measure:

self_tuning_knn = TunedModel(
    model=knn,
    resampling = CV(nfolds=4, rng=1234),
    tuning = Grid(resolution=5),
    range = K_range,
    measure=MisclassificationRate()
)

mach = machine(self_tuning_knn, X, y);
fit!(mach, verbosity=0);

trained Machine; does not cache data
  model: ProbabilisticTunedModel(model = KNNClassifier(K = 5, …), …)
  args: 
    1:	Source @020 ⏎ Table{AbstractVector{Continuous}}
    2:	Source @333 ⏎ AbstractVector{Multiclass{3}}

Let's inspect the best model and corresponding evaluation of the metric in case (ii):

entry = report(mach).best_history_entry

(model = KNNClassifier(K = 5, …),
 measure = StatisticalMeasuresBase.RobustMeasure{StatisticalMeasuresBase.FussyMeasure{StatisticalMeasuresBase.RobustMeasure{StatisticalMeasuresBase.Multimeasure{StatisticalMeasuresBase.SupportsMissingsMeasure{StatisticalMeasures.MisclassificationRateOnScalars}, Nothing, StatisticalMeasuresBase.Mean, typeof(identity)}}, Nothing}}[MisclassificationRate()],
 measurement = [0.02666666666666667],
 per_fold = [[0.0, 0.02631578947368421, 0.0, 0.08108108108108109]],)

entry.model.K

5

Recall that fitting mach also retrains the optimal model on all available data. The following is therefore an optimal model prediction based on all available data:

predict(mach, rows=148:150)

3-element UnivariateFiniteVector{Multiclass{3}, String, UInt32, Float64}:
 UnivariateFinite{Multiclass{3}}(setosa=>0.0, versicolor=>0.0, virginica=>1.0)
 UnivariateFinite{Multiclass{3}}(setosa=>0.0, versicolor=>0.0, virginica=>1.0)
 UnivariateFinite{Multiclass{3}}(setosa=>0.0, versicolor=>0.0, virginica=>1.0)

Specifying a custom measure

Users may specify a custom loss or scoring function, so long as it complies with the StatisticalMeasuresBase.jl API and implements the appropriate orientation trait (Score() or Loss()) from that package. For example, we suppose define a "new" scoring function custom_accuracy by

custom_accuracy(yhat, y) = mean(y .== yhat); # yhat - prediction, y - ground truth

custom_accuracy (generic function with 1 method)

In tuning, scores are maximised, while losses are minimised. So here we declare

import StatisticalMeasuresBase as SMB
SMB.orientation(::typeof(custom_accuracy)) = SMB.Score()

For full details on constructing custom measures, see StatisticalMeasuresBase.jl.

self_tuning_knn = TunedModel(
    model=knn,
    resampling = CV(nfolds=4),
    tuning = Grid(resolution=5),
    range = K_range,
    measure = [custom_accuracy, MulticlassFScore()],
    operation = predict_mode
);

mach = machine(self_tuning_knn, X, y)
fit!(mach, verbosity=0)
entry = report(mach).best_history_entry

(model = KNNClassifier(K = 5, …),
 measure = StatisticalMeasuresBase.RobustMeasure[Main.custom_accuracy, MulticlassFScore(beta = 1.0, …)],
 measurement = [0.8866666666666667, 0.4660896415303711],
 per_fold = [[1.0, 0.9210526315789472, 0.918918918918919, 0.7027027027027027], [0.33333333333333337, 0.6462585034013605, 0.6083530338849489, 0.2751322751322752]],)

entry.model.K

5

Tuning multiple nested hyperparameters

The forest model below has another model, namely a DecisionTreeRegressor, as a hyperparameter:

tree = Tree() # defined above
forest = EnsembleModel(model=tree)

DeterministicEnsembleModel(
  model = DecisionTreeRegressor(
        max_depth = -1, 
        min_samples_leaf = 5, 
        min_samples_split = 2, 
        min_purity_increase = 0.0, 
        n_subfeatures = 0, 
        post_prune = false, 
        merge_purity_threshold = 1.0, 
        feature_importance = :impurity, 
        rng = Random._GLOBAL_RNG()), 
  atomic_weights = Float64[], 
  bagging_fraction = 0.8, 
  rng = Random._GLOBAL_RNG(), 
  n = 100, 
  acceleration = CPU1{Nothing}(nothing), 
  out_of_bag_measure = Any[])

Ranges for nested hyperparameters are specified using dot syntax. In this case, we will specify a goal for the total number of grid points:

r1 = range(forest, :(model.n_subfeatures), lower=1, upper=9);
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0);
self_tuning_forest = TunedModel(
    model=forest,
    tuning=Grid(goal=30),
    resampling=CV(nfolds=6),
    range=[r1, r2],
    measure=rms);

X = MLJ.table(rand(100, 10));
y = 2X.x1 - X.x2 + 0.05*rand(100);

mach = machine(self_tuning_forest, X, y);
fit!(mach, verbosity=0);

trained Machine; does not cache data
  model: DeterministicTunedModel(model = DeterministicEnsembleModel(model = DecisionTreeRegressor(max_depth = -1, …), …), …)
  args: 
    1:	Source @409 ⏎ Table{AbstractVector{Continuous}}
    2:	Source @566 ⏎ AbstractVector{Continuous}

We can plot the grid search results:

using Plots
plot(mach)

Instead of specifying a goal, we can declare a global resolution, which is overridden for a particular parameter by pairing its range with the resolution desired. In the next example, the default resolution=100 is applied to the r2 field, but a resolution of 3 is applied to the r1 field. Additionally, we ask that the grid points be randomly traversed and the total number of evaluations be limited to 25.

tuning = Grid(resolution=100, shuffle=true, rng=1234)
self_tuning_forest = TunedModel(
    model=forest,
    tuning=tuning,
    resampling=CV(nfolds=6),
    range=[(r1, 3), r2],
    measure=rms,
    n=25
);
fit!(machine(self_tuning_forest, X, y), verbosity=0);

trained Machine; does not cache data
  model: DeterministicTunedModel(model = DeterministicEnsembleModel(model = DecisionTreeRegressor(max_depth = -1, …), …), …)
  args: 
    1:	Source @277 ⏎ Table{AbstractVector{Continuous}}
    2:	Source @657 ⏎ AbstractVector{Continuous}

For more options for a grid search, see Grid below.

Tuning using a random search

Let's attempt to tune the same hyperparameters using a RandomSearch tuning strategy. By default, bounded numeric ranges like r1 and r2 are sampled uniformly (before rounding, in the case of the integer range r1). Positive unbounded ranges are sampled using a Gamma distribution by default, and all others using a (truncated) normal distribution.

self_tuning_forest = TunedModel(
    model=forest,
    tuning=RandomSearch(),
    resampling=CV(nfolds=6),
    range=[r1, r2],
    measure=rms,
    n=25
);
X = MLJ.table(rand(100, 10));
y = 2X.x1 - X.x2 + 0.05*rand(100);
mach = machine(self_tuning_forest, X, y);
fit!(mach, verbosity=0)

trained Machine; does not cache data
  model: DeterministicTunedModel(model = DeterministicEnsembleModel(model = DecisionTreeRegressor(max_depth = -1, …), …), …)
  args: 
    1:	Source @956 ⏎ Table{AbstractVector{Continuous}}
    2:	Source @549 ⏎ AbstractVector{Continuous}

using Plots
plot(mach)

The prior distributions used for sampling each hyperparameter can be customized, as can the global fallbacks. See the RandomSearch doc-string below for details.

Tuning using Latin hypercube sampling

One can also tune the hyperparameters using the LatinHypercube tuning strategy. This method uses a genetic-based optimization algorithm based on the inverse of the Audze-Eglais function, using the library LatinHypercubeSampling.jl.

We'll work with the data X, y and ranges r1 and r2 defined above and instantiate a Latin hypercube resampling strategy:

latin = LatinHypercube(gens=2, popsize=120)

LatinHypercube(
  gens = 2, 
  popsize = 120, 
  ntour = 2, 
  ptour = 0.8, 
  interSampleWeight = 1.0, 
  ae_power = 2, 
  periodic_ae = false, 
  rng = Random._GLOBAL_RNG())

Here gens is the number of generations to run the optimisation for and popsize is the population size in the genetic algorithm. For more on these and other LatinHypercube parameters refer to the LatinHypercubeSampling.jl documentation. Pay attention that gens and popsize are not to be confused with the iteration parameter n in the construction of a corresponding TunedModel instance, which specifies the total number of models to be evaluated, independent of the tuning strategy.

For this illustration we'll add a third, nominal, hyper-parameter:

r3 = range(forest, :(model.post_prune), values=[true, false]);
self_tuning_forest = TunedModel(
    model=forest,
    tuning=latin,
    resampling=CV(nfolds=6),
    range=[r1, r2, r3],
    measure=rms,
    n=25
);
mach = machine(self_tuning_forest, X, y);
fit!(mach, verbosity=0)

trained Machine; does not cache data
  model: DeterministicTunedModel(model = DeterministicEnsembleModel(model = DecisionTreeRegressor(max_depth = -1, …), …), …)
  args: 
    1:	Source @241 ⏎ Table{AbstractVector{Continuous}}
    2:	Source @508 ⏎ AbstractVector{Continuous}

using Plots
plot(mach)

Comparing models of different type and nested cross-validation

Instead of mutating hyperparameters of a fixed model, one can instead optimise over an explicit list of models, whose types are allowed to vary. As with other tuning strategies, evaluating the resulting TunedModel itself implies nested resampling (e.g., nested cross-validation) which we now examine in a bit more detail.

tree = (@load DecisionTreeClassifier pkg=DecisionTree verbosity=0)()
knn = (@load KNNClassifier pkg=NearestNeighborModels verbosity=0)()
models = [tree, knn]

The following model is equivalent to the best in models by using 3-fold cross-validation:

multi_model = TunedModel(
    models=models,
    resampling=CV(nfolds=3),
    measure=log_loss,
    check_measure=false
)

Note that there is no need to specify a tuning strategy or range but we do specify models (plural) instead of model. Evaluating multi_model implies nested cross-validation (each model gets evaluated 2 x 3 times):

X, y = make_blobs()

e = evaluate(multi_model, X, y, resampling=CV(nfolds=2), measure=log_loss, verbosity=6)

PerformanceEvaluation object with these fields:
  model, measure, operation, measurement, per_fold,
  per_observation, fitted_params_per_fold,
  report_per_fold, train_test_rows, resampling, repeats
Extract:
┌──────────────────────┬───────────┬─────────────┬─────────┬────────────────────
│ measure              │ operation │ measurement │ 1.96*SE │ per_fold          ⋯
├──────────────────────┼───────────┼─────────────┼─────────┼────────────────────
│ LogLoss(             │ predict   │ 0.36        │ 0.999   │ [0.721, 2.22e-16] ⋯
│   tol = 2.22045e-16) │           │             │         │                   ⋯
└──────────────────────┴───────────┴─────────────┴─────────┴────────────────────

Now, for example, we can get the best model for the first fold out of the two folds:

e.report_per_fold[1].best_model

KNNClassifier(
  K = 5, 
  algorithm = :kdtree, 
  metric = Distances.Euclidean(0.0), 
  leafsize = 10, 
  reorder = true, 
  weights = NearestNeighborModels.Uniform())

And the losses in the outer loop (these still have to be matched to the best performing model):

e.per_fold

1-element Vector{Vector{Float64}}:
 [0.7208730677823432, 2.220446049250313e-16]

It is also possible to get the results for the nested evaluations. For example, for the first fold of the outer loop and the second model:

e.report_per_fold[2].history[1]

(model = DecisionTreeClassifier(max_depth = -1, …),
 measure = StatisticalMeasuresBase.RobustMeasure{StatisticalMeasuresBase.FussyMeasure{StatisticalMeasuresBase.RobustMeasure{StatisticalMeasures._LogLossType{Float64}}, typeof(StatisticalMeasures.l2_check)}}[LogLoss(tol = 2.22045e-16)],
 measurement = [2.1626192033470293],
 per_fold = [[2.1202149052421855, 2.1202149052421855, 2.252728336819822]],)

Reference

r = range(model, :hyper; values=nothing)

r = range(model, :hyper; upper=nothing, lower=nothing,
          scale=nothing, values=nothing)

iterator([rng, ], r::NominalRange, [,n])
iterator([rng, ], r::NumericRange, n)

sampler(r::NominalRange, probs::AbstractVector{<:Real})
sampler(r::NominalRange)
sampler(r::NumericRange{T}, d)

rand(s)             # for one sample
rand(s, n)          # for n samples
rand(rng, s [, n])  # to specify an RNG

r = range(Char, :letter, values=collect("abc"))
s = sampler(r, [0.1, 0.2, 0.7])
samples =  rand(s, 1000);
StatsBase.countmap(samples)
Dict{Char,Int64} with 3 entries:
  'a' => 107
  'b' => 205
  'c' => 688

r = range(Int, :k, lower=2, upper=6) # numeric but discrete
s = sampler(r, Normal)
samples = rand(s, 1000);
UnicodePlots.histogram(samples)
           ┌                                        ┐
[2.0, 2.5) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 119
[2.5, 3.0) ┤ 0
[3.0, 3.5) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 296
[3.5, 4.0) ┤ 0
[4.0, 4.5) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 275
[4.5, 5.0) ┤ 0
[5.0, 5.5) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 221
[5.5, 6.0) ┤ 0
[6.0, 6.5) ┤▇▇▇▇▇▇▇▇▇▇▇ 89
           └                                        ┘

Distributions.fit(D, r::MLJBase.NumericRange)

tuned_model = TunedModel(; model=<model to be mutated>,
                         tuning=RandomSearch(),
                         resampling=Holdout(),
                         range=nothing,
                         measure=nothing,
                         n=default_n(tuning, range),
                         operation=nothing,
                         other_options...)

tuned_model = TunedModel(; models=<models to be compared>,
                         resampling=Holdout(),
                         measure=nothing,
                         n=length(models),
                         operation=nothing,
                         other_options...)

Grid(goal=nothing, resolution=10, rng=Random.GLOBAL_RNG, shuffle=true)

range(model, :hyper1, lower=1, origin=2, unit=1)

[(range(model, :hyper1, lower=1, upper=10), 15),
  range(model, :hyper2, lower=2, upper=4),
  range(model, :hyper3, values=[:ball, :tree])]

# a range generating the grid `[1, 2, 10, 20, 30]` for `:hyper1`:
[range(model, :hyper1, values=[1, 2]),
 (range(model, :hyper1, lower= 10, upper=30), 3)]

RandomSearch(bounded=Distributions.Uniform,
             positive_unbounded=Distributions.Gamma,
             other=Distributions.Normal,
             rng=Random.GLOBAL_RNG)

using Distributions

range1 = range(model, :hyper1, lower=0, upper=1)

range2 = [(range(model, :hyper1, lower=1, upper=10), Arcsine),
          range(model, :hyper2, lower=2, upper=Inf, unit=1, origin=3),
          (range(model, :hyper2, lower=2, upper=4), Normal(0, 3)),
          (range(model, :hyper3, values=[:ball, :tree]), [0.3, 0.7])]

# uniform sampling of :(atom.λ) from [0, 1] without defining a NumericRange:
struct MySampler end
Base.rand(rng::Random.AbstractRNG, ::MySampler) = rand(rng)
range3 = (:(atom.λ), MySampler())

LatinHypercube(gens = 1,
               popsize = 100,
               ntour = 2,
               ptour = 0.8.,
               interSampleWeight = 1.0,
               ae_power = 2,
               periodic_ae = false,
               rng=Random.GLOBAL_RNG)

tuned_model = TunedModel(model=...,
                         tuning=LatinHypercube(...),
                         range=...,
                         measures=...,
                         n=...)