12  Przeszukiwanie iteracyjne

W poprzednich rozdziałach pokazano, jak wyszukiwanie oparte o siatkę przyjmuje wstępnie zdefiniowany zestaw wartości kandydujących, ocenia je, a następnie wybiera najlepsze ustawienia. Iteracyjne metody wyszukiwania realizują inną strategię. Podczas procesu wyszukiwania przewidują one, które wartości należy przetestować w następnej kolejności.

W tym rozdziale przedstawiono dwie metody wyszukiwania. Najpierw omówimy optymalizację bayesowską, która wykorzystuje model statystyczny do przewidywania lepszych ustawień parametrów. Następnie rozdział opisuje metodę globalnego wyszukiwania zwaną symulowanym wyżarzaniem (ang. simulated annealing).

Do ilustracji wykorzystujemy te same dane dotyczące charakterystyki komórek, co w poprzednim rozdziale, ale zmieniamy model. W tym rozdziale użyjemy modelu maszyny wektorów nośnych (ang. Support Vector Machine), ponieważ zapewnia on ładne dwuwymiarowe wizualizacje procesów wyszukiwania. Dwa dostrajane parametry w optymalizacji to wartość kosztu SVM i parametr jądra funkcji radialnej \(\sigma\). Oba parametry mogą mieć znaczący wpływ na złożoność i wydajność modelu.

Model SVM wykorzystuje iloczyn skalarny i z tego powodu konieczne jest wyśrodkowanie i przeskalowanie predyktorów. Obiekty tidymodels: svm_rec, svm_spec i svm_wflow definiują proces tworzenia modelu:

library(tidymodels)
library(finetune)
tidymodels_prefer()

data(cells)
cells <- cells %>% select(-case)

set.seed(1304)
cell_folds <- vfold_cv(cells)

svm_rec <- 
  recipe(class ~ ., data = cells) %>%
  step_YeoJohnson(all_numeric_predictors()) %>%
  step_normalize(all_numeric_predictors())

svm_spec <- 
  svm_rbf(cost = tune(), rbf_sigma = tune()) %>% 
  set_engine("kernlab") %>% 
  set_mode("classification")

svm_wflow <- 
  workflow() %>% 
  add_model(svm_spec) %>% 
  add_recipe(svm_rec)

Oto domyślne zakresy tuningowanych parametrów:

cost()
rbf_sigma()

Dla ilustracji, zmieńmy nieco zakres parametrów jądra, aby poprawić wizualizacje wyszukiwania:

svm_param <- 
  svm_wflow %>% 
  extract_parameter_set_dials() %>% 
  update(rbf_sigma = rbf_sigma(c(-7, -1)))

Zanim omówimy szczegóły dotyczące wyszukiwania iteracyjnego i jego działania, zbadajmy związek między dwoma parametrami dostrajania SVM a obszarem pod krzywą ROC dla tego konkretnego zestawu danych. Skonstruowaliśmy bardzo dużą regularną siatkę, składającą się z 2500 wartości kandydujących, i oceniliśmy siatkę przy użyciu resamplingu. Jest to oczywiście niepraktyczne w regularnej analizie danych i ogromnie nieefektywne. Jednakże, wyjaśnia ścieżkę, którą powinien podążać proces wyszukiwania i gdzie występuje numerycznie optymalna wartość (wartości).

Rys. 12.1: Mapa ciepła średniego obszaru pod krzywą ROC dla siatki o dużej gęstości wartości dostrajania parametrów. Najlepszy punkt jest oznaczony kropką w prawym górnym rogu.

Rys. 12.1 pokazuje wyniki oceny tej siatki, przy czym jaśniejszy kolor odpowiada wyższej (lepszej) wydajności modelu. W dolnej przekątnej przestrzeni parametrów znajduje się duży obszar, który jest stosunkowo płaski i charakteryzuje się słabą wydajnością. Grzbiet najlepszej wydajności występuje w prawej górnej części przestrzeni. Czarna kropka wskazuje najlepsze ustawienia. Przejście od płaskiego obszaru słabych wyników do grzbietu najlepszej wydajności jest bardzo ostre. Występuje również gwałtowny spadek obszaru pod krzywą ROC tuż po prawej stronie grzbietu.

Procedury wyszukiwania, które będziemy przestawiać, wymagają przynajmniej kilku statystyk obliczonych na podstawie resamplingu. W tym celu poniższy kod tworzy małą regularną siatkę, która znajduje się w płaskiej części przestrzeni parametrów. Funkcja tune_grid() dokonuje próbkowania tej siatki:

roc_res <- metric_set(roc_auc)

set.seed(1401)
start_grid <- 
  svm_param %>% 
  update(
    cost = cost(c(-6, 1)),
    rbf_sigma = rbf_sigma(c(-6, -4))
  ) %>% 
  grid_regular(levels = 2)

set.seed(1402)
svm_initial <- 
  svm_wflow %>% 
  tune_grid(resamples = cell_folds, grid = start_grid, metrics = roc_res)

collect_metrics(svm_initial)

# A tibble: 4 × 8
    cost rbf_sigma .metric .estimator  mean     n std_err .config             
   <dbl>     <dbl> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1 0.0156  0.000001 roc_auc binary     0.864    10 0.00864 Preprocessor1_Model1
2 2       0.000001 roc_auc binary     0.863    10 0.00867 Preprocessor1_Model2
3 0.0156  0.0001   roc_auc binary     0.863    10 0.00862 Preprocessor1_Model3
4 2       0.0001   roc_auc binary     0.866    10 0.00855 Preprocessor1_Model4

Ta początkowa siatka pokazuje dość równoważne wyniki, przy czym żaden pojedynczy punkt nie jest znacznie lepszy od pozostałych. Wyniki te mogą być wykorzystane przez funkcje iteracyjnego dostrajania jako wartości początkowe.

12.1 Optymalizacja bayesowska

Techniki optymalizacji bayesowskiej analizują bieżące wyniki próbkowania i tworzą model predykcyjny, aby zasugerować wartości parametrów dostrajania, które nie zostały jeszcze ocenione. Sugerowana kombinacja parametrów jest następnie ponownie próbkowana. Wyniki te są następnie wykorzystywane w innym modelu predykcyjnym, który rekomenduje więcej wartości kandydatów do testowania, i tak dalej. Proces ten przebiega przez ustaloną liczbę iteracji lub do momentu, gdy nie pojawią się dalsze poprawy. Shahriari i in. (2016) i Frazier (2018) są dobrymi wprowadzeniami do optymalizacji bayesowskiej.

Podczas korzystania z optymalizacji bayesowskiej, podstawowe problemy to sposób tworzenia modelu i wybór parametrów rekomendowanych przez ten model. Najpierw rozważmy technikę najczęściej stosowaną w optymalizacji bayesowskiej, czyli model procesu gaussowskiego.

12.1.1 Model procesu gaussowskiego

Modele procesu gaussowskiego (GP) (Schulz, Speekenbrink, i Krause 2016), to techniki statystyczne, które mają swoją historię w statystyce przestrzennej (pod nazwą metod krigingu). Mogą być wyprowadzone na wiele sposobów, w tym jako model Bayesowski.

Matematycznie, GP jest zbiorem zmiennych losowych, których wspólny rozkład prawdopodobieństwa jest wielowymiarowy normalny. W kontekście naszych zastosowań jest to zbiór metryk wydajności dla wartości kandydujących parametrów dostrajania. Dla początkowej siatki czterech próbek, realizacje tych czterech zmiennych losowych wynosiły 0.8639, 0.8625, 0.8627 i 0.8659. Zakłada się, że mają rozkład wielowymiarowy normalny. Wejściami definiującymi zmienne niezależne dla modelu GP są odpowiednie wartości dostrajania parametrów (przedstawione w Tab. 12.1).

Tab. 12.1: Wartości startowe w procedurze poszukiwania optymalnych parametrów

ROC	cost	rbf_sigma
0.8639	0.01562	1e-06
0.8625	2.00000	1e-06
0.8627	0.01562	1e-04
0.8659	2.00000	1e-04

Modele procesów gaussowskich są określone przez ich funkcje średniej i kowariancji, choć to ta ostatnia ma większy wpływ na charakter modelu GP. Funkcja kowariancji jest często parametryzowana w kategoriach wartości wejściowych (oznaczanych jako \(x\)). Przykładowo, powszechnie stosowaną funkcją kowariancji jest funkcja wykładnicza kwadratowa:

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

gdzie \(\sigma_{i,j}^2\) jest wariancją błędu modelu równą zero jeśli \(i=j\). Możemy to interpretować jako, że wraz ze wzrostem odległości pomiędzy dwoma kombinacjami parametrów, kowariancja pomiędzy metrykami wydajności rośnie wykładniczo. Z równania wynika również, że zmienność metryki wynikowej jest minimalizowana w punktach, które już zostały zaobserwowane (tzn. gdy \(|x_i - x_j|^2\) wynosi zero). Charakter tej funkcji kowariancji pozwala procesowi gaussowskiemu reprezentować wysoce nieliniowe zależności między wydajnością modelu a dostrajaniem parametrów, nawet jeśli istnieje tylko niewielka ilość danych.

Ważną zaletą tego modelu jest to, że ponieważ określony jest pełny model prawdopodobieństwa, przewidywania dla nowych wejść mogą odzwierciedlać cały rozkład wyniku. Innymi słowy, nowe statystyki wydajności mogą być przewidywane zarówno pod względem średniej jak i wariancji.

Na podstawie początkowej siatki czterech wyników, model GP jest dopasowywany, następnie są obliczane z modelu predykcje dla kandydatów, a piąta kombinacja parametrów dostrajania jest wybierana. Obliczamy szacunkową wydajność dla nowej konfiguracji, GP jest ponownie dopasowywany do pięciu istniejących wyników (i tak dalej).

12.1.2 Funkcja akwizycji

Jak wykorzystać model procesu gaussowskiego po dopasowaniu do aktualnych danych? Naszym celem jest wybranie następnej kombinacji dostrajania parametrów, która najprawdopodobniej da “lepsze wyniki” niż obecna najlepsza. Jednym z podejść do tego jest stworzenie dużego zbioru kandydatów, a następnie wykonanie prognoz średniej i wariancji dla każdego z nich. Korzystając z tych informacji, wybieramy najkorzystniejszą wartość parametru dostrajania.

Klasa funkcji celu, zwanych funkcjami akwizycji, ułatwia kompromis pomiędzy średnią a wariancją. Przypomnijmy, że przewidywana wariancja modeli GP zależy głównie od tego, jak bardzo są one oddalone od istniejących danych. Kompromis pomiędzy przewidywaną średnią i wariancją dla nowych kandydatów jest często postrzegany przez pryzmat eksploracji i eksploatacji:

• Eksploracja (ang. exploration) - powoduje wybór tych regionów, w których jest mniej obserwowanych modeli kandydujących. W ten sposób nadaje się większą wagę kandydatom o wyższej wariancji i koncentruje się na poszukiwaniu nowych wyników.
• Eksploatacja (ang. exploitation) - zasadniczo opiera się istniejących wynikach, w celu odnalezienia najlepszej wartości średniej.

Aby zademonstrować, spójrzmy na przykład z jednym parametrem, który ma wartości pomiędzy [0, 1], a metryką wydajności jest \(R^2\). Prawdziwa funkcja jest pokazana na Rys. 12.2 wraz z pięcioma wartościami kandydującymi, które mają istniejące wyniki jako punkty.

Rys. 12.2: Hipotetyczny rzeczywisty profil wydajności dla arbitralnie wybranego parametru dostrajania, z pięcioma szacowanymi punktami

Dla tych danych dopasowanie modelu GP przedstawiono na Rys. 12.3. Zacieniowany obszar wskazuje średnią \(\pm\) 1 błąd standardowy. Dwie pionowe linie wskazują dwa punkty kandydujące, które są później bardziej szczegółowo badane. Zacieniowany obszar ufności pokazuje funkcję wykładniczej wariancji kwadratowej; staje się ona bardzo duża między punktami i zbiega do zera w istniejących punktach danych.

Rys. 12.3: Przykładowy przebieg procesu gaussowskiego z zaznaczoną funkcją średniej i wariancji

Ta nieliniowa funkcja przechodzi przez każdy obserwowany punkt, ale model nie jest doskonały. Nie ma obserwowanych punktów w pobliżu prawdziwego optimum ustawienia, a w tym regionie dopasowanie mogłoby być znacznie lepsze. Pomimo tego, model GP może skutecznie wskazać nam właściwy kierunek.

Z punktu widzenia czystej eksploatacji, najlepszym wyborem byłoby wybranie wartości parametru, który ma najlepszą średnią predykcję. W tym przypadku byłaby to wartość 0,106, tuż na prawo od istniejącego najlepszego zaobserwowanego punktu na poziomie 0,09.

Jako sposób na zachęcenie do eksploracji, prostym (ale nie często stosowanym) podejściem jest znalezienie parametru dostrajania związanego z największym przedziałem ufności.

Jedną z najczęściej stosowanych funkcji akwizycji jest oczekiwana poprawa. Pojęcie poprawy wymaga wartości dla bieżących najlepszych wyników (w przeciwieństwie do podejścia opartego na przedziałach ufności). Ponieważ GP może opisać nowy punkt kandydacki za pomocą rozkładu, możemy ważyć fragmenty rozkładu, które wykazują poprawę, używając prawdopodobieństwa wystąpienia poprawy.

Na przykład, rozważmy dwie wartości parametrów kandydujących 0,10 i 0,25 (wskazane przez pionowe linie na Rys. 12.3). Używając dopasowanego modelu GP, ich przewidywane \(R^2\) są pokazane na Rys. 12.3 wraz z linią odniesienia dla aktualnych najlepszych wyników.

Rys. 12.4: Przewidywane rozkłady wydajności dla dwóch próbkowanych wartości parametrów dostrajania

Rozpatrując tylko średnią \(R^2\) lepszym wyborem jest wartość parametru 0,10 (patrz Tab. 12.1). Rekomendacja parametru dostrajania dla 0,25 ma gorsze przewidywanie średnie niż aktualny najlepszy kandydat. Jednakże, ponieważ ma wyższą wariancję, ma większy ogólny obszar prawdopodobieństwa powyżej aktualnego najlepszego. W rezultacie ma większą oczekiwaną poprawę:

tab2 <- tibble::tribble(
  ~Parameter.Value,  ~Mean,  ~Std.Dev, ~Expected.Improvment,
               0.1, 0.8679, 0.0004317,              0.00019,
              0.25, 0.8671, 0.0039301,             0.001216
  )
gt(tab2)

Tab. 12.2: Oczekiwana poprawa dla dwóch kandydujących parametrów dostrajania.

Parameter.Value	Mean	Std.Dev	Expected.Improvment
0.10	0.8679	0.0004317	0.000190
0.25	0.8671	0.0039301	0.001216

Kiedy oczekiwana poprawa jest obliczana w całym zakresie dostrajania parametrów, zalecany punkt do próbkowania jest znacznie bliższy 0,25 niż 0,10, jak pokazano na Rys. 12.5.

Rys. 12.5: Szacowany profil wydajności wygenerowany przez model procesu gaussowskiego (górny panel) oraz oczekiwana poprawa (dolny panel). Pionowa linia wskazuje punkt maksymalnej poprawy.

Aby zaimplementować wyszukiwanie iteracyjne poprzez optymalizację bayesowską, należy użyć funkcji tune_bayes(). Jej składnia jest bardzo podobna do tune_grid(), z kilkoma dodatkowymi argumentami:

• iter to maksymalna liczba iteracji wyszukiwania.
• initial może być liczbą całkowitą, obiektem utworzonym przy pomocy tune_grid(), albo jedną z funkcji wyścigowych. Użycie liczby całkowitej określa rozmiar konstrukcji wypełniającej przestrzeń, która jest próbkowana przed pierwszym modelem GP.
• objective jest argumentem, dla którego należy użyć funkcji akwizycji. Pakiet tune zawiera funkcje takie jak exp_improve() lub conf_bound().
• Argument param_info, w tym przypadku określa zakres parametrów, jak również wszelkie transformacje.

Argument control funkcji tune_bayes() ustawia się za pomocą control_bayes(). Niektóre istotne argumenty to:

• no_improve to liczba całkowita, która zatrzyma wyszukiwanie, jeśli ulepszone parametry nie zostaną odkryte w ciągu iteracji no_improve.
• uncertain jest również liczbą całkowitą (lub Inf), używaną do ustalenia liczby przejść algorytmu wg reguły eksploatacji bez poprawy, aby następnie wybrać próbę z zakresu z wysoką wariancją, po to żeby dokonać eksploracji.
• verbose jest parametrem, który decyduje co będzie się wyświetlało podczas przebiegu algorytmu.

Użyjmy pierwszych wyników SVM jako początkowego podłoża dla modelu GP. Przypomnijmy, że w tym zastosowaniu chcemy zmaksymalizować obszar pod krzywą ROC.

ctrl <- control_bayes(verbose = TRUE)

set.seed(1403)
svm_bo <-
  svm_wflow %>%
  tune_bayes(
    resamples = cell_folds,
    metrics = roc_res,
    initial = svm_initial,
    param_info = svm_param,
    iter = 25,
    control = ctrl
  )

collect_metrics(svm_bo) |> 
  gt()

Tab. 12.3: Przebieg optymalizacji baysowskiej

cost	rbf_sigma	.metric	.estimator	mean	n	std_err	.config	.iter
0.015625000	1.000000e-06	roc_auc	binary	0.8638724	10	0.008637894	Preprocessor1_Model1	0
2.000000000	1.000000e-06	roc_auc	binary	0.8625326	10	0.008672545	Preprocessor1_Model2	0
0.015625000	1.000000e-04	roc_auc	binary	0.8627495	10	0.008624554	Preprocessor1_Model3	0
2.000000000	1.000000e-04	roc_auc	binary	0.8659439	10	0.008545691	Preprocessor1_Model4	0
0.110831031	2.665887e-02	roc_auc	binary	0.8882389	10	0.008697493	Iter1	1
0.002570159	1.498759e-02	roc_auc	binary	0.8733510	10	0.008784422	Iter2	2
0.076743095	3.257686e-02	roc_auc	binary	0.8820147	10	0.008962769	Iter3	3
0.158777868	8.699292e-04	roc_auc	binary	0.8646797	10	0.008717232	Iter4	4
0.028725163	2.614512e-02	roc_auc	binary	0.8748373	10	0.009077821	Iter5	5
0.110836919	4.063308e-07	roc_auc	binary	0.8640364	10	0.008588054	Iter6	6
0.161280423	2.376079e-02	roc_auc	binary	0.8912983	10	0.008546152	Iter7	7
0.191548445	1.753100e-02	roc_auc	binary	0.8924207	10	0.008604149	Iter8	8
0.322506923	2.024837e-02	roc_auc	binary	0.8951338	10	0.008488570	Iter9	9
0.616915747	2.547641e-02	roc_auc	binary	0.8953431	10	0.008324862	Iter10	10
0.528787926	2.052531e-02	roc_auc	binary	0.8954957	10	0.008263982	Iter11	11
0.552359295	2.381201e-02	roc_auc	binary	0.8954230	10	0.008235443	Iter12	12
0.518642133	2.172084e-02	roc_auc	binary	0.8957039	10	0.008241457	Iter13	13
0.699416391	2.187984e-02	roc_auc	binary	0.8954909	10	0.008317203	Iter14	14
22.086197445	3.848670e-02	roc_auc	binary	0.8647523	10	0.009452636	Iter15	15
0.415152612	2.230321e-02	roc_auc	binary	0.8951613	10	0.008306736	Iter16	16
0.001157284	2.435740e-03	roc_auc	binary	0.8656808	10	0.008630693	Iter17	17
0.924338505	7.885258e-02	roc_auc	binary	0.8858168	10	0.008167638	Iter18	18
29.083907487	1.192898e-03	roc_auc	binary	0.8979424	10	0.008064529	Iter19	19
31.924080961	6.212242e-04	roc_auc	binary	0.8942442	10	0.008618062	Iter20	20
10.599233064	1.275900e-03	roc_auc	binary	0.8952500	10	0.008543058	Iter21	21
30.740481343	1.907822e-03	roc_auc	binary	0.8985736	10	0.007908688	Iter22	22
24.204296906	1.712143e-03	roc_auc	binary	0.8983898	10	0.007782671	Iter23	23
31.939709359	1.703217e-03	roc_auc	binary	0.8988150	10	0.007769977	Iter24	24
0.001100257	1.000011e-07	roc_auc	binary	0.8639531	10	0.008702518	Iter25	25

show_best(svm_bo) |> 
  gt()

Tab. 12.4: Zestaw optymalnych parametrów na podstawie OB

cost	rbf_sigma	.metric	.estimator	mean	n	std_err	.config	.iter
31.9397094	0.001703217	roc_auc	binary	0.8988150	10	0.007769977	Iter24	24
30.7404813	0.001907822	roc_auc	binary	0.8985736	10	0.007908688	Iter22	22
24.2042969	0.001712143	roc_auc	binary	0.8983898	10	0.007782671	Iter23	23
29.0839075	0.001192898	roc_auc	binary	0.8979424	10	0.008064529	Iter19	19
0.5186421	0.021720838	roc_auc	binary	0.8957039	10	0.008241457	Iter13	13

autoplot(svm_bo, type = "performance")

Rys. 12.6: Jakość dopasowania modelu w poszczególnych iteracjach OB

Poniższa animacja wizualizuje wyniki wyszukiwania. Czarne “x” pokazują wartości początkowe zawarte w svm_initial. Niebieski panel u góry po lewej stronie pokazuje przewidywaną średnią wartość obszaru pod krzywą ROC. Czerwony panel na górze po prawej stronie pokazuje przewidywaną zmienność wartości ROC, podczas gdy dolny wykres wizualizuje oczekiwaną poprawę. W każdym panelu ciemniejsze kolory wskazują mniej atrakcyjne wartości (np. małe wartości średnie, duże zróżnicowanie i małe ulepszenia).

Video

Rys. 12.7: Przykład działania optymalizacji bayesowskiej

Powierzchnia przewidywanej średniej jest bardzo niedokładna w pierwszych kilku iteracjach wyszukiwania. Pomimo tego, pomaga ona poprowadzić proces w rejon dobrej wydajności. W ciągu pierwszych dziesięciu iteracji, wyszukiwanie jest dokonywane w pobliżu optymalnego miejsca.

Podczas gdy najlepsza kombinacja dostrajania parametrów znajduje się na granicy przestrzeni parametrów, optymalizacja bayesowska często wybierze nowe punkty poza granicami przestrzeni parametrów.

Powyższy przykład startował z punktów startowych wybranych nieco arbitralnie ale nieco lepsze wyniki można osiągnąć stosując losowe wypełnienie przestrzeni parametrów.

12.2 Symulowane wyżarzanie

Symulowane wyżarzanie (ang. simulated annealing) (Kirkpatrick, Gelatt, i Vecchi 1983; Laarhoven i Aarts 1987); jest nieliniową procedurą wyszukiwania zainspirowaną procesem stygnięcia metalu. Jest to metoda globalnego wyszukiwania, która może efektywnie poruszać się po wielu różnych obszarach poszukiwań, w tym po funkcjach nieciągłych. W przeciwieństwie do większości procedur optymalizacji opartych na gradiencie, symulowane wyżarzanie może ponownie ocenić poprzednie rozwiązania.

Proces użycia symulowanego wyżarzania rozpoczyna się od wartości początkowej i rozpoczyna kontrolowany losowy spacer przez przestrzeń parametrów. Każda nowa wartość parametru-kandydata jest niewielką perturbacją poprzedniej wartości, która utrzymuje nowy punkt w lokalnym sąsiedztwie.

Punkt kandydujący jest oceniana przy zastosowaniu resamplingu, aby uzyskać odpowiadającą mu wartość wydajności. Jeśli osiąga ona lepsze wyniki niż poprzednie parametry, jest akceptowana jako nowa najlepsza i proces jest kontynuowany. Jeśli wyniki są gorsze niż poprzednia wartość, procedura wyszukiwania może nadal używać tego parametru do określenia dalszych kroków. Zależy to od dwóch czynników. Po pierwsze, prawdopodobieństwo zatrzymania złego kandydata maleje wraz z pogorszeniem się wyników. Innymi słowy, tylko nieco gorszy wynik od obecnie najlepszego ma większą szansę na akceptację niż ten z dużym spadkiem wydajności. Drugim czynnikiem jest liczba iteracji wyszukiwania. Symulowane wyżarzanie próbuje zaakceptować mniej suboptymalnych wartości w miarę postępu wyszukiwania. Z tych dwóch czynników prawdopodobieństwo akceptacji złego wyniku można sformalizować jako:

\[ \operatorname{Pr}[\text{accept suboptimal parameters at iteration } i] = \exp(c\times D_i \times i) \tag{12.2}\]

gdzie \(i\) jest numerem iteracji, \(c\) jest stałą określoną przez użytkownika, \(D_i\) jest procentową różnicą pomiędzy starą i nową wartością (gdzie wartości ujemne oznaczają gorsze wyniki). Dla złego wyniku określamy prawdopodobieństwo akceptacji i porównujemy je z liczbą wylosowaną z rozkładu jednostajnego. Jeśli liczba ta jest większa od wartości prawdopodobieństwa, wyszukiwanie odrzuca bieżące parametry i następna iteracja tworzy swoją wartość kandydata w sąsiedztwie poprzedniej wartości. W przeciwnym razie następna iteracja tworzy kolejny zestaw parametrów na podstawie bieżących (suboptymalnych) wartości.

Zagrożenie

Prawdopodobieństwa akceptacji symulowanego wyżarzania pozwalają na postępowanie w złym kierunku, przynajmniej na krótką metę, z potencjałem znalezienia znacznie lepszego regionu przestrzeni parametrów w dłuższej perspektywie.

Mapa ciepła na Rys. 12.8 pokazuje, jak prawdopodobieństwo akceptacji może się zmieniać w zależności od iteracji, wydajności i współczynnika określonego przez użytkownika.

Rys. 12.8: Mapa ciepła prawdopodobieństwa akceptacji symulowanego wyżarzania dla różnych wartości współczynnika

Użytkownik może dostosować współczynniki \(c\), aby znaleźć profil prawdopodobieństwa, który odpowiada jego potrzebom. W finetune::control_sim_anneal(), domyślnym dla tego argumentu cooling_coef jest 0.02. Zmniejszenie tego współczynnika zachęci wyszukiwanie do bycia bardziej wyrozumiałym dla słabych wyników.

Proces ten trwa przez określoną ilość iteracji, ale może zostać zatrzymany, jeśli w ciągu określonej liczby iteracji nie pojawią się globalnie najlepsze wyniki. Bardzo pomocne może być ustawienie progu restartu. Jeśli wystąpi ciąg niepowodzeń, funkcja ta powraca do ostatnich globalnie najlepszych ustawień parametrów i zaczyna od nowa.

Najważniejszym szczegółem jest określenie sposobu perturbacji parametrów dostrajania z iteracji na iterację. W literaturze można znaleźć wiele metod na to. My stosujemy metodę podaną przez Bohachevsky, Johnson, i Stein (1986) zwaną uogólnionym symulowanym wyżarzaniem. Dla ciągłych parametrów dostrajania definiujemy mały promień, aby określić lokalne “sąsiedztwo”. Na przykład załóżmy, że są dwa parametry dostrajania i każdy z nich jest ograniczony przez zero i jeden. Proces symulowanego wyżarzania generuje losowe wartości na otaczającym promieniu i losowo wybiera jedną z nich jako aktualną wartość kandydacką.

Wielkość promienia kontroluje, jak szybko wyszukiwanie bada przestrzeń parametrów. Im większy jest promień tym szybciej przestrzeń parametrów będzie przeszukiwana przez algorytm ale jednocześnie mniej dokładnie.

Dla zilustrowania użyjemy dwóch głównych parametrów dostrajania glmnet:

• Wielkość kary (penalty). Domyślny zakres tego parametru to od \(10^{-10}\) do \(10^0\). Najczęściej jest podawany w skali logarytmicznej o podstawie 10.
• Proporcja kary lasso (mixture) - wartość pomiędzy 0 i 1 określająca balans pomiędzy karą L1 i L2.

Proces rozpoczyna się od wartości początkowych penalty = 0,025 i mixture = 0,050. Używając promienia, który losowo waha się między 0,050 a 0,015, dane są odpowiednio skalowane, losowe wartości są generowane na promieniach wokół punktu początkowego, a następnie jedna jest losowo wybierana jako kandydat. Dla ilustracji przyjmiemy, że wszystkie wartości kandydujące faktycznie poprawiają wydajność modelu. Korzystając z nowej wartości, generowany jest zestaw nowych losowych sąsiadów, wybierany jest jeden, i tak dalej. Rys. 12.9 przedstawia sześć iteracji sukcesywnie podążających do lewego górnego rogu.

Rys. 12.9: Kilka iteracji procesu symulowanego wyżarzania

Zauważmy, że podczas niektórych iteracji zestawy kandydatów wzdłuż promienia wykluczają punkty poza granicami parametrów. Ponadto, nasza implementacja przesuwa wybór kolejnych konfiguracji dostrajania parametrów od nowych wartości, które są bardzo podobne do poprzednich konfiguracji.

Dla parametrów kategorycznych i całkowitych, każdy z nich jest generowany z określonym wcześniej prawdopodobieństwem. Argument flip w control_sim_anneal() może być użyty do określenia tego prawdopodobieństwa. Dla parametrów całkowitoliczbowych, używana jest najbliższa wartość całkowita.

Przeszukiwanie metodą symulowanego wyżarzania może nie być optymalna w przypadku, gdy wiele parametrów jest nienumerycznych lub całkowitoliczbowych o niewielu unikalnych wartościach. W tych przypadkach jest prawdopodobne, że ten sam zestaw kandydatów może być testowany więcej niż raz.

Aby zaimplementować wyszukiwanie iteracyjne poprzez symulowane wyżarzanie, użyj funkcji tune_sim_anneal(). Składnia tej funkcji jest niemal identyczna jak tune_bayes(). Nie ma żadnych opcji dotyczących funkcji akwizycji czy próbkowania niepewności. Funkcja control_sim_anneal() posiada pewne szczegóły, które definiują lokalne sąsiedztwo i harmonogram chłodzenia:

• no_improve jest liczbą całkowitą, która zatrzyma wyszukiwanie, jeśli w ciągu iteracji określonym przez no_improve nie zostaną odkryte globalnie najlepsze wyniki. Przyjęte suboptymalne lub odrzucone parametry liczą się jako “brak poprawy”.
• restart to liczba iteracji, która musi minąć bez poprawy jakości modelu, po której następuje przejście do najlepszego poprzednio ustalonego zestawu parametrów.
• radius to wektor liczbowy w przedziale (0, 1), który określa minimalny i maksymalny promień lokalnego sąsiedztwa wokół punktu początkowego.
• flip to wartość prawdopodobieństwa określająca szanse zmiany wartości parametrów kategorycznych lub całkowitych.
• cooling_coef jest współczynnikiem \(c\) w \(\exp(c\times D_i \times i)\), który moduluje jak szybko prawdopodobieństwo akceptacji maleje w trakcie iteracji. Większe wartości cooling_coef zmniejszają prawdopodobieństwo akceptacji suboptymalnego ustawienia parametrów.

Dla danych dotyczących segmentacji komórek składnia jest bardzo spójna z poprzednio stosowanymi funkcjami:

ctrl_sa <- control_sim_anneal(verbose = TRUE, no_improve = 10L)

set.seed(1404)
svm_sa <-
  svm_wflow %>%
  tune_sim_anneal(
    resamples = cell_folds,
    metrics = roc_res,
    initial = svm_initial,
    param_info = svm_param,
    iter = 50,
    control = ctrl_sa
  )

Dane wyjściowe dla poszczególnych iteracji:

collect_metrics(svm_sa) |> 
  gt()

Tab. 12.5: Podsumowanie symulowanego wyżarzania

cost	rbf_sigma	.metric	.estimator	mean	n	std_err	.config	.iter
0.0156250	1.000000e-06	roc_auc	binary	0.8638724	10	0.008637894	initial_Preprocessor1_Model1	0
2.0000000	1.000000e-06	roc_auc	binary	0.8625326	10	0.008672545	initial_Preprocessor1_Model2	0
0.0156250	1.000000e-04	roc_auc	binary	0.8627495	10	0.008624554	initial_Preprocessor1_Model3	0
2.0000000	1.000000e-04	roc_auc	binary	0.8659439	10	0.008545691	initial_Preprocessor1_Model4	0
2.7792215	4.231877e-05	roc_auc	binary	0.8635104	10	0.008642153	Iter1	1
8.9751798	1.317101e-04	roc_auc	binary	0.8733437	10	0.008401879	Iter2	2
8.0489544	3.113457e-05	roc_auc	binary	0.8672404	10	0.008380301	Iter3	3
7.4631879	6.401369e-05	roc_auc	binary	0.8704169	10	0.008147365	Iter4	4
25.4276256	6.964316e-05	roc_auc	binary	0.8744436	10	0.008459084	Iter5	5
26.0442746	1.620597e-05	roc_auc	binary	0.8696570	10	0.008110480	Iter6	6
27.8247939	6.830484e-06	roc_auc	binary	0.8658910	10	0.008447841	Iter7	7
11.5246937	1.459848e-06	roc_auc	binary	0.8623315	10	0.008657469	Iter8	8
11.5911446	4.938754e-06	roc_auc	binary	0.8622997	10	0.008647646	Iter9	9
16.0739132	2.533161e-06	roc_auc	binary	0.8623219	10	0.008658226	Iter10	10
22.9676231	1.507608e-06	roc_auc	binary	0.8623751	10	0.008670690	Iter11	11
14.6174389	3.625603e-07	roc_auc	binary	0.8624605	10	0.008645675	Iter12	12
30.2128582	1.354001e-07	roc_auc	binary	0.8625366	10	0.008612533	Iter13	13
14.1812532	7.415027e-05	roc_auc	binary	0.8726092	10	0.008351119	Iter14	14
5.1974958	1.892985e-04	roc_auc	binary	0.8727956	10	0.008301061	Iter15	15
6.2799813	1.101844e-03	roc_auc	binary	0.8897887	10	0.008859945	Iter16	16
1.7551679	2.895039e-03	roc_auc	binary	0.8918849	10	0.008922869	Iter17	17
2.7450163	7.152477e-03	roc_auc	binary	0.8958891	9	0.009100454	Iter18	18
0.6580640	7.237199e-03	roc_auc	binary	0.8944414	10	0.008826226	Iter19	19
0.2894265	2.800490e-03	roc_auc	binary	0.8781751	10	0.008519248	Iter20	20
0.8648804	9.061603e-04	roc_auc	binary	0.8744488	10	0.008274188	Iter21	21
2.6900002	4.370896e-04	roc_auc	binary	0.8748006	10	0.008467929	Iter22	22
1.9738868	8.518383e-04	roc_auc	binary	0.8787410	10	0.008594452	Iter23	23
2.8390760	1.343734e-04	roc_auc	binary	0.8694466	10	0.008255440	Iter24	24
0.6443050	2.818297e-04	roc_auc	binary	0.8660129	10	0.008642925	Iter25	25
0.4332128	8.494378e-04	roc_auc	binary	0.8702872	10	0.008407344	Iter26	26
1.5062749	4.117203e-02	roc_auc	binary	0.8868902	10	0.008384564	Iter27	27
1.4313131	1.549218e-02	roc_auc	binary	0.8963429	10	0.008367098	Iter28	28
2.7699626	1.362418e-02	roc_auc	binary	0.8929540	10	0.008489741	Iter29	29
4.9266582	4.583077e-02	roc_auc	binary	0.8698810	10	0.009097231	Iter30	30
7.7758900	9.490760e-03	roc_auc	binary	0.8874485	10	0.008767151	Iter31	31
4.8976312	1.580979e-02	roc_auc	binary	0.8825035	10	0.009000858	Iter32	32
12.4918485	8.493479e-03	roc_auc	binary	0.8829449	10	0.008897035	Iter33	33
9.4138354	2.056779e-03	roc_auc	binary	0.8971685	10	0.008342496	Iter34	34
16.7499203	2.398614e-03	roc_auc	binary	0.8985179	10	0.007916078	Iter35	35
18.8139030	4.403197e-04	roc_auc	binary	0.8839150	9	0.009020320	Iter36	36
21.9409340	1.216970e-03	roc_auc	binary	0.8973481	10	0.008280488	Iter37	37
8.7296351	5.154540e-03	roc_auc	binary	0.8961246	10	0.008263055	Iter38	38
13.5493211	2.655684e-02	roc_auc	binary	0.8593977	10	0.009612964	Iter39	39
3.4729934	1.312798e-03	roc_auc	binary	0.8847509	9	0.009126495	Iter40	40
3.1103793	4.869128e-03	roc_auc	binary	0.8969694	10	0.008168916	Iter41	41
1.2292254	6.540634e-03	roc_auc	binary	0.8965152	10	0.008550536	Iter42	42
0.5769315	2.846348e-03	roc_auc	binary	0.8842725	10	0.008540335	Iter43	43
22.9584171	1.241411e-02	roc_auc	binary	0.8618903	10	0.009740857	Iter44	44
10.6297370	4.059190e-03	roc_auc	binary	0.8975137	10	0.008203966	Iter45	45
10.6164560	2.519507e-02	roc_auc	binary	0.8617737	10	0.009517440	Iter46	46
24.5183146	7.016539e-03	roc_auc	binary	0.8775665	10	0.008934404	Iter47	47
10.2189037	1.330918e-02	roc_auc	binary	0.8733744	10	0.009258351	Iter48	48
7.1021781	1.966405e-03	roc_auc	binary	0.8966379	10	0.008539778	Iter49	49
9.2002334	5.308068e-04	roc_auc	binary	0.8841932	10	0.008839438	Iter50	50

show_best(svm_sa) |> 
  gt()

Tab. 12.6: Najlepsze kombinacje parametrów na podstawie SA

cost	rbf_sigma	.metric	.estimator	mean	n	std_err	.config	.iter
16.749920	0.002398614	roc_auc	binary	0.8985179	10	0.007916078	Iter35	35
10.629737	0.004059190	roc_auc	binary	0.8975137	10	0.008203966	Iter45	45
21.940934	0.001216970	roc_auc	binary	0.8973481	10	0.008280488	Iter37	37
9.413835	0.002056779	roc_auc	binary	0.8971685	10	0.008342496	Iter34	34
3.110379	0.004869128	roc_auc	binary	0.8969694	10	0.008168916	Iter41	41

Podobnie jak w przypadku innych funkcji tune_*(), odpowiadająca im funkcja autoplot() tworzy wizualną ocenę wyników. Użycie autoplot(svm_sa, type = "performance") pokazuje wydajność w czasie iteracji (Rys. 12.10), podczas gdy autoplot(svm_sa, type = "parameters") przedstawia wydajność w zależności od konkretnych wartości dostrajania parametrów (Rys. 12.11).

autoplot(svm_sa, type = "performance")

Rys. 12.10: Przebieg procesu symulowanego wyżarzania

autoplot(svm_sa, type = "parameters")

Rys. 12.11: Wydajność modeli w kontekście hiperparametrów modelu

Wizualizacja ścieżki wyszukiwania pomaga zrozumieć, gdzie proces wyszukiwania poradził sobie dobrze, a gdzie pobłądził:

Video

Rys. 12.12: Wizualizacja przebiegu symulowanego wyżarzania

Bohachevsky, Ihor O., Mark E. Johnson, i Myron L. Stein. 1986. „Generalized Simulated Annealing for Function Optimization”. Technometrics 28 (3): 209–17. https://doi.org/10.1080/00401706.1986.10488128.

Frazier, Peter I. 2018. „A Tutorial on Bayesian Optimization”. https://doi.org/10.48550/ARXIV.1807.02811.

Kirkpatrick, S., C. D. Gelatt, i M. P. Vecchi. 1983. „Optimization by Simulated Annealing”. Science 220 (4598): 671–80. https://doi.org/10.1126/science.220.4598.671.

Laarhoven, Peter J. M. van, i Emile H. L. Aarts. 1987. „Simulated annealing”. W, 7–15. Springer Netherlands. https://doi.org/10.1007/978-94-015-7744-1_2.

Schulz, Eric, Maarten Speekenbrink, i Andreas Krause. 2016. „A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions”. http://dx.doi.org/10.1101/095190.

Shahriari, Bobak, Kevin Swersky, Ziyu Wang, Ryan P. Adams, i Nando de Freitas. 2016. „Taking the Human Out of the Loop: A Review of Bayesian Optimization”. Proceedings of the IEEE 104 (1): 148–75. https://doi.org/10.1109/jproc.2015.2494218.