mlr3MBOを用いたベイズ最適化によるハイパーパラメータチューニング

# 1. はじめに
機械学習では、lightGBMなどの決定木の勾配ブースティングモデルは、チューニングすべきハイパーパラメータが多いため、ベイズ最適化が良く用いられます。
以前に、Tidymodelを用いたベイズ最適化によるlightGBMのハイパーパラメータチューニングをしました。
[【R】tydymodelsを用いてlightgbmをやってみる（その３）](https://qiita.com/kenkenvw/items/d605050a88636cd677e7)

しかしながら、現在、上記はエラーが発生しています。原因は,

'reshape' argument is no longer supported.

とでるので、どうやらパッケージ側の問題みたいです。

そこで、今回はMLR3を用いたいと思います。
MLR3とは、Tidymodelと同じくRによる機械学習フレームワークですが、文法が異なるので、好みが分かれるところですが、実は私はこちらの文法が気に入っています。

しかしながら、マイナーな存在であるので、この機会に紹介してみたいと思います。
[mlr3を使ってみる（その１）](https://qiita.com/kenkenvw/items/c2c27d917a95b4a8b8a5)
MLR3によるベイズ最適化は、以前に投稿したところです。
[mlr3MBOを用いたベイズ最適化について（Bayesian Optimization）](https://qiita.com/kenkenvw/items/409d9f9a24d99f2f1cde)
上記は、ハイパーパラメータのチューニングではないため、今回、mlr3MBOをハイパーパラメータチューニングに適用してみます。

環境
```
R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
RStudio 2024.12.0 Build 467
OS:Ubuntu 22.04
```



# 2.DATA
ライブラリ⁠
```r
library(AmesHousing)
library(janitor)
library(tidyverse)
library(mlr3)
library(mlr3learners)
library(mlr3extralearners)  #lightGBMが含まれる
library(mlr3tuning)
library(mlr3pipelines)
library(mlr3mbo)
library(mlr3tuning)
library(bbotk)
library(lightgbm)
```
AmesHoisingパッケージで入手可能な住宅価格のデータを使用します。
```r
# set the random seed so we can reproduce any simulated results.
set.seed(123)
# load the housing data and clean names
ames_data = make_ames() %>% clean_names()
```
janitorのclean_names()関数は、自動でデータ列名をクリーニングしてくれる。どのような作業が実行されたというと、
Before
```
"MS_SubClass"        "MS_Zoning"          "Lot_Frontage"     "Kitchen_AbvGr"
```
After
```
 "ms_sub_class"       "ms_zoning"          "lot_frontage"     "kitchen_abv_gr"  
```
大文字を小文字で統一や、アンダーバーで区切りを自動で入れてくれたりします。地味に便利。

# 2.タスクの設定
多くのモデルでは、正確な予測を行うために、慎重かつ広範な変数の前処理を必要とします。
XGBoost、lightgbm、catboostのような決定木の勾配ブースティングツリーモデルは、高度に歪んだデータや相関性のあるデータに対して非常に頑健であるため、必要な前処理の量は最小限で済みます。XGBoostとは対照的に、lightgbmとcatboostはカテゴリー変数（因子）を扱うことができるので、変数をワンホットエンコードする必要はありません。
データをトレーニングデータとテストデータにわけます。
まず、データは、MLR3では最初にタスク化します。
```r
tsk_ames = TaskRegr$new(id="ames",backend = ames_data,target = "sale_price")
tsk_ames

#> <TaskRegr:ames> (2930 x 81)
#> * Target: sale_price
#> * Properties: -
#> * Features (80):
#>   - fct (46): alley, bldg_type, bsmt_cond, bsmt_exposure,
#>     bsmt_fin_type_1, bsmt_fin_type_2, bsmt_qual, central_air,
#>     condition_1, condition_2, electrical, exter_cond, exter_qual,
#>     exterior_1st, exterior_2nd, fence, fireplace_qu, foundation,
#>     functional, garage_cond, garage_finish, garage_qual,
#>     garage_type, heating, heating_qc, house_style, kitchen_qual,
#>     land_contour, land_slope, lot_config, lot_shape, mas_vnr_type,
#>     misc_feature, ms_sub_class, ms_zoning, neighborhood,
#>     overall_cond, overall_qual, paved_drive, pool_qc, roof_matl,
#>     roof_style, sale_condition, sale_type, street, utilities
#>   - int (22): bedroom_abv_gr, enclosed_porch, fireplaces,
#>     first_flr_sf, full_bath, gr_liv_area, half_bath, kitchen_abv_gr,
#>     lot_area, low_qual_fin_sf, misc_val, mo_sold, open_porch_sf,
#>     pool_area, screen_porch, second_flr_sf, three_season_porch,
#>     tot_rms_abv_grd, wood_deck_sf, year_built, year_remod_add,
#>     year_sold
#>   - dbl (12): bsmt_fin_sf_1, bsmt_fin_sf_2, bsmt_full_bath,
#>     bsmt_half_bath, bsmt_unf_sf, garage_area, garage_cars, latitude,
#>     longitude, lot_frontage, mas_vnr_area, total_bsmt_sf
```
# 2.学習器の設定
lightGBMの学習器を設定するとともに、ハイパーパラメータのチューニング設定をします。
lightGBMでは主に次のハイパーパラメータがあります。
+ num_leaves
+ feature_fraction 
+ num_iterations 
+ min_data_in_leaf 
+ max_depth 
+ learning_rate 
+ min_gain_to_split 
+ bagging_fraction 

今回の最適化をねらうターゲットは、num_leaves、min_data_in_leaf 、feature_fraction の３パラメータとし、learning_rate を0.1に、num_iterations を500で固定します。
```r
lgm_lrn = lrn("regr.lightgbm",
             learning_rate = 0.1, 
             num_leaves = to_tune(10,200),
             num_iterations = 500,
             feature_fraction = to_tune(0.5,1),
             min_data_in_leaf = to_tune(1, 100)
)
```
# 3.前処理との組み合わせ
ここでは、ハイパーパラメータのチューニングが主目的なので、前処理は最小限とします。
前処理として、factor列は全てラベルエンコーディングし、ターゲットをlog変換します。
ファクター化すると自動的に水準が振られるため、ラベルエンコーディングはas.numeric（数値化）するだけで良いです。pipline処理を用いて実行します。

po("colapply")は、列に対して関数を適用するパイプラインで、as.numeric関数を、factor列にだけに適用させます。学習器とパイプラインを組み合わせて新たな学習器を作成します。
```r
ple_lrn = as_learner(po("colapply", affect_columns=selector_type("factor"), applicator = as.numeric)　%>>% ppl("targettrafo",
  graph = lgm_lrn,
  targetmutate.trafo = function(x) log(x+1),
  targetmutate.inverter = function(x) list(response = exp(x$response)-1)) )
ple_lrn$plot()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/d62e6b45-8d8a-ee79-256c-07d0452714a1.png)

# 4.ベイズ最適化の実行
サロゲートモデル（代理モデル）にガウス過程回帰を設定し、獲得関数を”EI”とし、獲得関数のオプティマイザーにDIRECT法を指定します。
```r
bayesopt_ego = mlr_loop_functions$get("bayesopt_ego")
surrogate = srlrn(lrn("regr.km", covtype = "matern5_2",
  optim.method = "BFGS", control = list(trace = FALSE)))
acq_function = acqf("ei")
acq_optimizer = acqo(opt("nloptr", algorithm = "NLOPT_GN_ORIG_DIRECT"),
  terminator = trm("stagnation", iters = 100, threshold = 1e-5))
```
ベイズ最適化は上記の設定でほとんどの場合、行けそうです。
tunerにtnr("mbo")で、上記の設定でループさせます。
なお、ハイパーパラメータの場合は、オプティマイザー関数[opt("mbo")]ではなく、チューナー関数[ tnr("mbo"] を使います。tune（）関数でtunerに学習器、タスク、リサンプリング、目的関数、最適化数をセットし、ベイズ最適化のループを回します。
```r
rsmp_cv = rsmp("cv", folds = 5)
msr_rsme = msr("regr.rmse")

tuner = tnr("mbo",
  loop_function = bayesopt_ego ,
  surrogate = surrogate,
  acq_function = acq_function,
  acq_optimizer = acq_optimizer)

instance = tune(
  tuner,
  tsk_ames,
  ple_lrn,
  rsmp_cv,
  msr_rsme,
  20
)

#> ....
#> ....
#> INFO  [20:28:47.947] [mlr3] Finished benchmark
#> INFO  [20:28:47.980] [bbotk] Result of batch 9:
#> INFO  [20:28:47.990] [bbotk]
#> INFO  [20:28:48.029] [bbotk] Finished optimizing after 20 evaluation(s)
#> INFO  [20:28:48.031] [bbotk] Result:
```
```r
library(mlr3viz)
autoplot(instance)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/cea6abbb-f4ed-76a5-707c-9c7758124e31.png)

num_leavesは明らかに小さい方が良い傾向が出ています。
では、最適なパラメータ設定を見ます。
```r
instance$result_learner_param_vals

#> $colapply.applicator
#> function (x, ...)  .Primitive("as.double")
#> 
#> $colapply.affect_columns
#> selector_type("factor")
#> 
#> $regr.lightgbm.objective
#> [1] "regression"
#> 
#> $regr.lightgbm.verbose
#> [1] -1
#> 
#> $regr.lightgbm.learning_rate
#> [1] 0.1
#> 
#> $regr.lightgbm.num_threads
#> [1] 1
#> 
#> $regr.lightgbm.num_iterations
#> [1] 500
#> 
#> $targetmutate.trafo
#> function(x) log(x+1)
#> 
#> $targetmutate.inverter
#> function(x) list(response = exp(x$response)-1)
#> 
#> $regr.lightgbm.num_leaves
#> [1] 15
#> 
#> $regr.lightgbm.min_data_in_leaf
#> [1] 2
#> 
#> $regr.lightgbm.feature_fraction
#> [1] 0.5243773
```

またnum_iterations 数を500としていましたが、最適化させます。early_stopping（早期停止）を使います。
```r
lgm_lrn$param_set$set_values(
  early_stopping_rounds = 10,
  learning_rate =0.1,
  num_iterations =500,
  num_leaves=15,
  min_data_in_leaf=2,
  feature_fraction=0.5243773
)
```
```r
set_validate(lgm_lrn ,0.2)
lgm_lrn$train(tsk_ames)
lgm_lrn$internal_tuned_values[[1]]

#> [1] 116
```
バリデーションを20%に設定し、早期停止を使って、イテレーション数が500回から116回まで下がりました。これで必要以上な学習回数を防いで、時間の短縮とモデルの過学習を防止します。

以下、全パラメータを設定します。ここで注意する点は、早期設定時に用いた early_stoppingやバリデーション設定を元に戻（NULL)しておかないとエラーがでます。
```r
lgm_lrn$param_set$set_values(
  learning_rate =0.1,
  num_iterations =116,
  num_leaves=15,
  min_data_in_leaf=2,
  feature_fraction=0.5243773,
  early_stopping_rounds=NULL
)
lgm_lrn$validate=NULL
```
前処理と組み合わせます。
```r
ple_lrn = as_learner(po("colapply", affect_columns=selector_type("factor"), applicator = as.numeric)　%>>% ppl("targettrafo",
  graph = lgm_lrn,
  targetmutate.trafo = function(x) log(x+1),
  targetmutate.inverter = function(x) list(response = exp(x$response)-1)) )
ple_lrn 

#> <GraphLearner:colapply.targetmutate.regr.lightgbm.targetinvert>
#> * Model: -
#> * Parameters: colapply.applicator=<function>,
#>   colapply.affect_columns=<Selector>,
#>   regr.lightgbm.objective=regression, regr.lightgbm.verbose=-1,
#>   regr.lightgbm.learning_rate=0.1, regr.lightgbm.num_leaves=15,
#>   regr.lightgbm.num_threads=1, regr.lightgbm.min_data_in_leaf=2,
#>   regr.lightgbm.feature_fraction=0.5244,
#>   regr.lightgbm.num_iterations=116,
#>   regr.lightgbm.early_stopping_rounds=10,
#>   targetmutate.trafo=<function>, targetmutate.inverter=<function>
#> * Validate: NULL
#> * Packages: mlr3, mlr3pipelines, mlr3extralearners, lightgbm
#> * Predict Types:  [response], se, quantiles, distr
#> * Feature Types: logical, integer, numeric, character, factor,
#>   ordered, POSIXct
#> * Properties: featureless, hotstart_backward, hotstart_forward,
#>   importance, internal_tuning, marshal, missings,
#>   selected_features, validation, weights
```

# 5.モデルの評価
ベイズ最適化したパラメータとデフォルト設定と比較してモデルの評価をします。
```r
lgm=lrn("regr.lightgbm")
ple_base_lrn = as_learner(po("colapply", affect_columns=selector_type("factor"), applicator = as.numeric)　%>>% ppl("targettrafo",
  graph = lgm,
  targetmutate.trafo = function(x) log(x+1),
  targetmutate.inverter = function(x) list(response = exp(x$response)-1)) )
ple_lrn$id="bays-tune"
ple_base_lrn$id="original"
```
デフォルトをoriginalとしました。
評価は、全データの5分割交差検証法により実施します。
ベンチマークのグリッド設計をします。学習器だけを比較しています。
```r
design = benchmark_grid(tsk_ames,c(ple_base_lrn,ple_lrn), rsmp_cv)
```
ベンチマークを実行します。
```r
bmr = benchmark(design)
```
目的関数により結果を出します。
```r
bmr$aggregate(msr("regr.rmse"))[, .(task_id, learner_id, regr.rmse)]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/10b5df3e-7735-584b-bc95-1f1f7b9c59b4.png)

リサンプリング結果分布を見ます。
```r
autoplot(bmr)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/74a0dcf5-c685-4390-1ab1-b5dfe6ef2566.png)
ベイズ最適化でチューニングすることで、RSMEが小さくなりました。

# 6.Enjoy!