Extras - Iterative search

Practical Machine Learning with tidymodels

Startup! 

library(tidymodels)
library(important)
library(probably)

tidymodels_prefer()
theme_set(theme_bw())
options(pillar.advice = FALSE, pillar.min_title_chars = Inf)
mirai::daemons(parallel::detectCores())

More startup! 

# Make our example data for this section
set.seed(429)
class_data <- 
  # Nonlinear classification data
  sim_classification(5000) |> 
  # Add some noise predictors that are correlated with one another
  bind_cols(sim_noise(5000, num_var = 15, cov_type = "toeplitz", cov_param = 0.5))

set.seed(429)
sim_split <- initial_split(class_data, prop = 0.75, strata = class)
sim_train <- training(sim_split)
sim_test  <- testing(sim_split)

set.seed(523)
sim_rs <- vfold_cv(sim_train, v = 10, strata = class)

Neural network

From the seventh advanced slide deck:

nnet_spec <-
  mlp(hidden_units = tune(), penalty = tune(), learn_rate = tune(),
      epochs = 100, activation = tune()
  ) |>
  set_engine("brulee", stop_iter = 10) |>
  set_mode("classification")

rec <- 
  recipe(class ~ ., data = sim_train) |> 
  step_normalize(all_numeric_predictors())
  
thrsh_tlr <-
  tailor() |>
  adjust_probability_threshold(threshold = tune()) 
  
nnet_wflow <- workflow(rec, nnet_spec, thrsh_tlr)

nnet_param <-
  nnet_wflow |> 
  extract_parameter_set_dials() |> 
  update(threshold = threshold(c(0.0001, 0.1)))
  
cls_mtr <- metric_set(brier_class, roc_auc, sensitivity, specificity)

Iterative Search

A linear model probably isn’t the best choice though (more in a minute).

To illustrate the process, we resampled a large grid of learning rate values for our data to show what the relationship is between error and learning rate.

Now suppose that we used a grid of three points in the parameter range for learning rate…

A Large Grid

A Three Point Grid

Gaussian Processes and Optimization

We can make a “meta-model” with a small set of historical performance results.

Gaussian Processes (GP) models are a good choice to model performance.

  • It is a Bayesian model so we are using Bayesian Optimization (BO).
  • For regression, we can assume that our data are multivariate normal.
  • We also define a covariance function for the variance relationship between data points. A common one is:

\[\operatorname{cov}(\boldsymbol{x}_i, \boldsymbol{x}_j) = \exp\left(-\frac{1}{2}|\boldsymbol{x}_i - \boldsymbol{x}_j|^2\right) + \sigma^2_{ij}\]

Predicting Candidates

The GP model can take candidate tuning parameter combinations as inputs and make predictions for performance (e.g. Brier, ROC AUC, RMSE, etc.)

  • The mean performance
  • The variance of performance

The variance is mostly driven by spatial variability (the previous equation).

The predicted variance is zero at locations of actual data points and becomes very high when far away from any observed data.

Your turn

Your GP makes predictions on two new candidate tuning parameters.

We want to minimize error.

Which should we choose?

GP Fit (ribbon is mean +/- 1SD)

Choosing New Candidates

This isn’t a very good fit but we can still use it.

How can we use the outputs to choose the next point to measure?


Acquisition functions take the predicted mean and variance and use them to balance:

  • exploration: new candidates should explore new areas.
  • exploitation: new candidates must stay near existing values.

Exploration focuses on the variance, exploitation is about the mean.

Acquisition Functions

We’ll use an acquisition function to select a new candidate.

The most popular method appears to be expected improvement (EI) above the current best results.

  • Zero at existing data points.
  • The expected improvement is integrated over all possible improvement (“expected” in the probability sense).

We would probably pick the point with the largest EI as the next point.

(There are other functions beyond EI.)

Expected Improvement

Iteration

Once we pick the candidate point, we measure performance for it (e.g. resampling).


Another GP is fit, EI is recomputed, and so on.


We stop when we have completed the allowed number of iterations or if we don’t see any improvement after a pre-set number of attempts.

GP Fit with four points

Expected Improvement

AML4TD, TMwR

BO in tidymodels

We’ll use a function called tune_bayes() that has very similar syntax to tune_grid().


It has an additional initial argument for the initial set of performance estimates and parameter combinations for the GP model.

Initial grid points

initial can be the results of another tune_*() function or an integer (in which case tune_grid() is used under to hood to make such an initial set of results).

  • We’ll run the optimization more than once, so let’s make an initial grid of results to serve as the substrate for the BO.

  • I suggest at least the number of tuning parameters plus two as the initial grid for BO.

Qualitative parameters

  • What about non-numeric tuning parameters such as activation?

  • Currently, tidymodels converts these to dummy indicators and uses those in the GP. This is not unusual but also not great.

    • Our initial grid should include more points; one for each level of the qualitative parameter.
    • In our case, the activation function is preset to use 5 possible values.

  • An upcoming version of tune will use a different R package to fit the GP that uses factor or Gower kernels. This will avoid making indicators and require fewer initial points.

An Initial Grid

set.seed(12)
init_res <-
  nnet_wflow |>
  tune_grid(
    resamples = sim_rs,
    grid = nrow(nnet_param) + 6, # for activation values + 1 extra
    param_info = nnet_param,
    metrics = cls_mtr
  )

show_best(init_res, metric = "brier_class", n = 3) |> select(-.metric, -.estimator)
#> # A tibble: 3 × 9
#>   hidden_units      penalty activation learn_rate threshold   mean     n std_err
#>          <int>        <dbl> <chr>           <dbl>     <dbl>  <dbl> <int>   <dbl>
#> 1           30 0.001        log_sigmo…     0.331     0.0101 0.0743    10 0.00298
#> 2           40 0.0000000001 elu            0.0132    0.0600 0.0744    10 0.00285
#> 3            6 0.0001       elu            0.0251    0.0800 0.0762    10 0.00266
#> # ℹ 1 more variable: .config <chr>

BO using tidymodels

ctrl_bo <- control_bayes(verbose_iter = TRUE, no_improve = Inf) 

set.seed(125)
nnet_bayes_res <-
  nnet_wflow |>
  tune_bayes(
    resamples = sim_rs,
    initial = init_res,     # <- initial results
    iter = 25,
    control = ctrl_bo,
    param_info = nnet_param,
    metrics = cls_mtr
  )
#> Optimizing brier_class using the expected improvement
#> 
#> ── Iteration 1 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 14.1%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=8, penalty=0.000173, activation=tanh, learn_rate=0.154,
#>   threshold=0.0936
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.08848 (+/-0.00524)
#> 
#> ── Iteration 2 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 93.3%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=26, penalty=0.000276, activation=elu, learn_rate=0.106,
#>   threshold=0.0588
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07849 (+/-0.00322)
#> 
#> ── Iteration 3 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 34.4%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=32, penalty=1.99e-10, activation=elu, learn_rate=0.00856,
#>   threshold=0.00222
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.0757 (+/-0.00258)
#> 
#> ── Iteration 4 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 78.5%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=50, penalty=4.27e-08, activation=log_sigmoid,
#>   learn_rate=0.261, threshold=0.0609
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.08127 (+/-0.00709)
#> 
#> ── Iteration 5 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 61.7%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=18, penalty=2.42e-07, activation=relu, learn_rate=0.00106,
#>   threshold=0.00355
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07863 (+/-0.00248)
#> 
#> ── Iteration 6 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 43%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=24, penalty=5.05e-05, activation=log_sigmoid,
#>   learn_rate=0.00112, threshold=0.000702
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.1182 (+/-0.0137)
#> 
#> ── Iteration 7 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 84.6%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=50, penalty=1e-10, activation=relu, learn_rate=0.00204,
#>   threshold=0.00494
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07925 (+/-0.00316)
#> 
#> ── Iteration 8 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 19.3%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=32, penalty=2.41e-07, activation=elu, learn_rate=0.544,
#>   threshold=0.097
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07549 (+/-0.003)
#> 
#> ── Iteration 9 ───────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 36%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=42, penalty=0.143, activation=log_sigmoid, learn_rate=0.598,
#>   threshold=0.0701
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.2153 (+/-0.0101)
#> 
#> ── Iteration 10 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 83.6%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=3, penalty=0.0025, activation=elu, learn_rate=0.0029,
#>   threshold=0.0649
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.08331 (+/-0.00785)
#> 
#> ── Iteration 11 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07427 (@iter 0)
#> ✔ Gaussian process model (LOO R²: 97.2%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=32, penalty=8.51e-07, activation=log_sigmoid,
#>   learn_rate=0.175, threshold=0.0303
#> i Estimating performance
#> ✓ Estimating performance
#> ♥ Newest results:    brier_class=0.07422 (+/-0.00294)
#> 
#> ── Iteration 12 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 78.3%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=30, penalty=3.36e-07, activation=elu, learn_rate=0.217,
#>   threshold=0.029
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07808 (+/-0.00355)
#> 
#> ── Iteration 13 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 79.9%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=31, penalty=6.31e-07, activation=relu, learn_rate=0.341,
#>   threshold=0.0429
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07925 (+/-0.00386)
#> 
#> ── Iteration 14 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 84.3%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=12, penalty=0.00907, activation=elu, learn_rate=0.0882,
#>   threshold=0.00584
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07646 (+/-0.00241)
#> 
#> ── Iteration 15 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 60.1%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=13, penalty=0.00122, activation=relu, learn_rate=0.606,
#>   threshold=0.00433
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07853 (+/-0.00315)
#> 
#> ── Iteration 16 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 80.9%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=6, penalty=1.83e-07, activation=log_sigmoid, learn_rate=0.502,
#>   threshold=0.014
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07651 (+/-0.00293)
#> 
#> ── Iteration 17 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 84.8%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=17, penalty=3.39e-06, activation=elu, learn_rate=0.144,
#>   threshold=0.095
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07625 (+/-0.00341)
#> 
#> ── Iteration 18 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 80.6%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=49, penalty=0.00177, activation=elu, learn_rate=0.103,
#>   threshold=0.087
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07462 (+/-0.0027)
#> 
#> ── Iteration 19 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 81.6%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=49, penalty=3.05e-09, activation=elu, learn_rate=0.0332,
#>   threshold=0.0977
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07447 (+/-0.00299)
#> 
#> ── Iteration 20 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 81.1%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=44, penalty=1.52e-06, activation=tanh, learn_rate=0.168,
#>   threshold=0.00779
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.08905 (+/-0.00812)
#> 
#> ── Iteration 21 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 85.8%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=42, penalty=1.54e-08, activation=elu, learn_rate=0.0031,
#>   threshold=0.0551
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07514 (+/-0.00274)
#> 
#> ── Iteration 22 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 81.4%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=13, penalty=2.1e-10, activation=elu, learn_rate=0.00104,
#>   threshold=0.0978
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07869 (+/-0.00345)
#> 
#> ── Iteration 23 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 77.8%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=47, penalty=5.56e-08, activation=elu, learn_rate=0.00272,
#>   threshold=0.00438
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07482 (+/-0.00276)
#> 
#> ── Iteration 24 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 73.7%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=28, penalty=6.07e-10, activation=relu, learn_rate=0.401,
#>   threshold=0.0952
#> i Estimating performance
#> ✓ Estimating performance
#> ⓧ Newest results:    brier_class=0.07944 (+/-0.00283)
#> 
#> ── Iteration 25 ──────────────────────────────────────────────────────
#> 
#> i Current best:      brier_class=0.07422 (@iter 11)
#> ✔ Gaussian process model (LOO R²: 88.8%)
#> ℹ Generating 5000 candidates.
#> i hidden_units=47, penalty=1.53e-05, activation=elu, learn_rate=0.61,
#>   threshold=0.0373
#> i Estimating performance
#> ✓ Estimating performance
#> ♥ Newest results:    brier_class=0.07404 (+/-0.00307)

Best results

show_best(nnet_bayes_res, metric = "brier_class") |> select(-.metric, -.estimator)
#> # A tibble: 5 × 10
#>   hidden_units      penalty activation learn_rate threshold   mean     n std_err
#>          <int>        <dbl> <chr>           <dbl>     <dbl>  <dbl> <int>   <dbl>
#> 1           47     1.53e- 5 elu            0.610     0.0373 0.0740    10 0.00307
#> 2           32     8.51e- 7 log_sigmo…     0.175     0.0303 0.0742    10 0.00294
#> 3           30     1   e- 3 log_sigmo…     0.331     0.0101 0.0743    10 0.00298
#> 4           40     1   e-10 elu            0.0132    0.0600 0.0744    10 0.00285
#> 5           49     3.05e- 9 elu            0.0332    0.0977 0.0745    10 0.00299
#> # ℹ 2 more variables: .config <chr>, .iter <int>

Plotting BO Results

autoplot(nnet_bayes_res, metric = "brier_class")

Plotting BO Results

autoplot(nnet_bayes_res, metric = "brier_class", type = "parameters")

Plotting BO Results

autoplot(nnet_bayes_res, metric = "brier_class", type = "performance")

Your turn

Let’s try a different acquisition function: conf_bound(kappa).

We’ll use the objective argument to set it.

Choose your own kappa value:

  • Larger values will explore the space more.
  • “Large” values are usually less than one.

Bonus points: Before the optimization is done, press <esc> and see what happens.

Notes

  • Stopping tune_bayes() will return the current results.

  • Parallel processing can still be used to more efficiently measure each candidate point.

  • There are a lot of other iterative methods that you can use.

  • The finetune package also has functions for simulated annealing search.