Podręcznik

Kierunek minimalizacji $\textbf{p}$ w procesie optymalizacji wyznacza taki wektor kierunkowy w przestrzeni wielowymiarowej, w którym następuje zmniejszanie się wartości funkcji celu w następnym kroku. Metody wyznaczania takiego kierunku opierają się na rozwinięciu w szereg Taylora funkcji celu $E$ w najbliższym sąsiedztwie (wzdłuż kierunku $\textbf{p}$) znanego, aktualnego punktu rozwiązania określonego wektorem $\textbf{w}$ [14,18,45]

W rozwinięciu tym $\textbf{w}$ oznacza $n$-wymiarowy wektor wszystkich wag sieci, $\textbf{p}$ – wektor kierunku poszukiwań, $\textbf{g}$ – wektor gradientu a $\textbf{H}$ - macierz drugich pochodnych funkcji celu (hesjan) wyznaczane w aktualnym punkcie pracy $\textbf{w}$

W praktyce wykorzystuje się co najwyżej trzy pierwsze składniki wzoru (2.7) pomijając pozostałe. Zależność (2.7) jest wówczas kwadratowym przybliżeniem funkcji $E( \textbf{w} )$ w najbliższym sąsiedztwie znanego punktu rozwiązania $\textbf{w}$. Dla uproszczenia zapisu wartości zmiennych uzyskane w $k$-tym cyklu oznaczane będą symbolem ze wskaźnikiem dolnym $k$. Punkt rozwiązania $\textbf{w} = \textbf{w}_k$ jest punktem odpowiadającym minimum wartości funkcji celu $E$, jeśli pierwsza pochodna funkcji celu (gradient) jest równy zeru, a macierz hesjanu $\textbf{H}$ jest dodatnio określona (warunki Kuhna-Tuckera [14]). Przy spełnieniu tych warunków funkcja w dowolnym punkcie należącym do sąsiedztwa $\textbf{w}$ ma wartość większą niż w samym punkcie $\textbf{w}$, a zatem punkt rozwiązania aktualnego $\textbf{w}$ jest rozwiązaniem odpowiadającym minimalnej wartości funkcji celu.

Na $k$-tym etapie poszukiwania punktu minimalnego wartości funkcji celu kierunek poszukiwań $\textbf{p} = \textbf{p}_k$ oraz krok $\eta$ powinny być dobrane w taki sposób, aby dla nowego rozwiązania

była spełniona zależność  $E (\textbf{w}_{k+1}) < E (\textbf{w}_{k})$. Poszukiwania minimum funkcji celu trwają dopóty, dopóki norma gradientu nie spadnie poniżej założonej, dopuszczalnej wartości błędu bądź też nie będzie przekroczony określony czas obliczeń (liczba iteracji). Ogólny algorytm adaptacji wag sieci przedstawić można w następującej postaci (zakłada się, że dana jest wartość początkowa wektora optymalizowanego $\textbf{w}_{k} = \textbf{w}_{0}$.

1. Test na zbieżność i optymalność aktualnego rozwiązania $\textbf{w}_{k}$. Jeżeli aktualny punkt $\textbf{w}_{k}$ spełnia warunki zatrzymania procesu ze względu na gradient - koniec obliczeń. W przeciwnym razie przejście do punktu 2.

2. Wyznaczenie wektora kierunku minimalizacji $\textbf{p}_{k}$ w punkcie rozwiązania $\textbf{w}_{k}$.

3. Wyznaczenie takiego kroku $\eta$ na kierunku $\textbf{p}_{k}$, że spełniony jest warunek  $E(\textbf{w}_{k} +\eta \textbf{p}_k) < E (\textbf{w}_{k})$.

4. Określenie nowego punktu rozwiązania $\textbf{w}_{k+1} = \textbf{w}_{k} + \eta \textbf{p}_{k}$ oraz odpowiadającej mu wartości $E (\textbf{w}_{k+1}) < E (\textbf{w}_{k})$ i powrót do punktu 1.

2.5.1 Algorytm największego spadku

Przy ograniczeniu się w rozwinięciu Taylora (2.7) do liniowego przybliżenia funkcji celu w najbliższym sąsiedztwie aktualnego rozwiązania $\textbf{w}$ otrzymuje się algorytm największego spadku. Aby była spełniona zależność $E (\textbf{w}_{k+1}) < E (\textbf{w}_{k})$ , wystarczy dobrać $\textbf{g}_k^T \textbf{p}_k < 0$.

Przyjęcie wektora kierunkowego w postaci [18]

spełnia wymagany warunek zmniejszenia wartości funkcji celu. Adaptacja wag przebiega wówczas według wzoru

określającego metodę największego spadku. Ograniczenie się w rozwinięciu Taylora do składnika pierwszego rzędu powoduje niewykorzystanie informacji o krzywiźnie funkcji zawartej w hesjanie. To sprawia, że metoda jest wolno zbieżna (zbieżność liniowa: oznacza ona, iż przy spełnieniu założeń metody, odległości pomiędzy kolejnymi przybliżeniami $\mathbf{x}_k$, a rzeczywistym minimum $\mathbf{x}^{*}$ funkcji maleją liniowo, czyli $\left\|\mathbf{x}^*-\mathbf{x}_{k+1}\right\| \leq c\left\|\mathbf{x}^*-\mathbf{x}_k\right\|$). Wada ta oraz brak postępów minimalizacji w okolicy punktu optymalnego, gdzie gradient przyjmuje bardzo małe wartości, spowodowały, że metoda ta jest mało efektywna. Niemniej jednak, ze względu na jej prostotę, małe wymagania co do wielkości pamięci i małą złożoność obliczeniową stanowiła ona przez lata i nadal pozostaje podstawową metodą stosowaną w uczeniu sieci wielowarstwowych.

Poprawy efektywności poszukuje się przez modyfikację (zwykle heurystyczną) wzoru określającego kierunek. Dobre wyniki przynosi zwykle zastosowanie uczenia z momentem rozpędowym (w skrócie uczenie z momentem) W przypadku zastosowania tej techniki zmiana wektora wag sieci $\Delta\mathbf{w}$ w $k$-tej iteracji uwzględnia nie tylko gradient, ale również zmianę wartości wag z poprzedniej $(k-1)$ iteracji. Adaptacja wag przebiega wówczas według wzoru [43,46]

$\mathbf{w}_{k+1}=\mathbf{w}_k-\eta \mathbf{g}_k+\alpha\left(\mathbf{w}_k-\mathbf{w}_{k-1}\right)$

(2.11)

w którym $\alpha$ jest współczynnikiem momentu, przyjmowanym zwykle w przedziale wartości [0, 1].

Składnik $\eta \mathbf{g}_k$tego wyrażenia odpowiada zwykłej metodzie uczenia największego spadku, natomiast drugi $\alpha\left(\mathbf{w}_k-\mathbf{w}_{k-1}\right)$uwzględnia ostatnią zmianę wag i jest niezależny od aktualnej wartości gradientu. Im większa jest wartość współczynnika $\alpha$, tym składnik wynikający z momentu rozpędowego ma większy wpływ na adaptację wag. Jego wpływ wzrasta w sposób istotny na płaskich odcinkach funkcji celu oraz w pobliżu minimum lokalnego, gdzie wartość gradientu jest bliska zeru. Na płaskich odcinkach funkcji celu przyrost wag $\Delta \mathbf{w}$, przy stałej wartości współczynnika uczenia $\eta$ jest w przybliżeniu równy, co oznacza, że $\Delta \mathbf{w}_k=-\eta \mathbf{g}_k+\alpha \Delta \mathbf{w}_k$. Wynika stąd, że efektywny przyrost wartości wag określa teraz relacja

Przy wartości współczynnika $\alpha=0,9$ oznacza to 10-krotny wzrost efektywnej wartości współczynnika uczenia, a więc również 10-krotne przyśpieszenie procesu uczenia (na płaskich odcinkach funkcji celu). W pobliżu minimum lokalnego składnik momentu jako nie związany z wartością gradientu może spowodować zmianę wag, prowadzącą do wzrostu wartości funkcji celu i opuszczenie strefy ,,przyciągania" tego minimum. Przy małej wartości gradientu czynnik momentu może stać się dominujący we wzorze (2.11) i spowodować przyrost wag odpowiadający wzrostowi wartości funkcji celu i umożliwić opuszczenie bieżącego minimum lokalnego.

Czynnik momentu nie może przy tym całkowicie zdominować procesu uczenia, gdyż prowadziłoby to do niestabilności algorytmu. Zwykle kontroluje się wartość funkcji celu $E$ w procesie uczenia, dopuszczając do jej wzrostu jedynie w ograniczonym zakresie, np. 4%. W takim przypadku, jeśli w kolejnych iteracjach $k$ oraz $k+1$ spełniona jest relacja $E(k+1) < 1,04E(k)$, krok jest akceptowany i następuje uaktualnienie wartości wag; jeżeli natomiast $E(k+1) > 1,04E(k)$, zmiany są ignorowane i przyjmuje się $\left(\mathbf{w}_k-\mathbf{w}_{k-1}\right)=\mathbf{0}$. W takim przypadku składnik gradientowy odzyskuje dominację nad składnikiem momentu i proces przebiega zgodnie z kierunkiem minimalizacji wyznaczonym przez wektor gradientu. Należy podkreślić, że dobór wartości współczynnika momentu jest zagadnieniem niełatwym i wymaga wielu eksperymentów, mających na celu dopasowanie jego wartości do specyfiki rozwiązywanego problemu.

2.5.2 Algorytm zmiennej metryki

Metoda zmiennej metryki wywodzi się z metody optymalizacyjnej Newtona, wykorzystującej kwadratowe przybliżenie funkcji celu w sąsiedztwie znanego rozwiązania $\mathbf{w}_k$. Ograniczając się we wzorze (2.7) do trzech pierwszych składników, minimum funkcji celu $E$ wyznacza się przyrównując do zera pochodną $E(\textbf{w}+\textbf{p})$ względem wektora $\textbf{p}$, czyli $\frac{\partial E(\mathbf{w}+p)}{\partial \mathbf{p}}=\mathbf{0}$. Po wykonaniu tej operacji mamy

Ze wzoru tego otrzymuje się wyrażenie określające kierunek newtonowski wyznaczający minimum funkcji celu w $k$-tej iteracji [45]

Wzór (2.14) wyznacza w sposób jednoznaczny taki kierunek $\mathbf{p}_k$, który zapewnia minimum funkcji celu w danym kroku i stanowi istotę newtonowskiego algorytmu optymalizacji, jest zależnością czysto teoretyczną, gdyż wymaga zapewnienia dodatniej określoności hesjanu w każdym kroku, co w przypadku ogólnym jest niemożliwe do spełnienia.

Z tego powodu w praktycznych implementacjach algorytmu rezygnuje się z dokładnego wyznaczania hesjanu $\mathbf{H}$, a w zamian stosuje się jego przybliżenie $\mathbf{G}$ spełniające warunek dodatniej określoności. Jedną z najpopularniejszych metod jest przy tym metoda zmiennej metryki, posiadająca wiele odmian. W metodzie tej, w każdym kroku modyfikuje się hesjan (lub jego odwrotność) z kroku poprzedniego o pewną poprawkę. Oznaczmy przez $\mathbf{s}_{\mathrm{k}}=\mathbf{w}_{\mathrm{k}}-\mathbf{w}_{\mathrm{k}-1}$ oraz $\mathbf{r}_{\mathrm{k}}=\mathbf{g}_{\mathrm{k}}-\mathbf{g}_{\mathrm{k}-1}$, odpowiednio przyrosty wektora $\mathbf{w}$ oraz gradientu $\mathbf{g}$ w dwu kolejnych iteracjach, a przez $\mathbf{V}$ - macierz odwrotną przybliżonego hesjanu $\mathbf{V} = \mathbf{G}^{-1}$, to zgodnie z najefektywniejszą formułą Broydena-Fletchera-Goldfarba-Shanno (BFGS) proces uaktualniania macierzy $\mathbf{V}$ opisuje się zależnością rekurencyjną [14]

Jako wartość startową przyjmuje się zwykle $\mathbf{V}_{0}=\mathbf{1}$, co oznacza, że pierwsza iteracja przeprowadzana jest zgodnie z algorytmem największego spadku. Zostało wykazane [14], że przy wartości startowej $\mathbf{V}_{0}=\mathbf{1}$ i zastosowaniu procedury minimalizacji kierunkowej w każdym kroku optymalizacyjnym można zawsze zapewnić dodatnią określoność aproksymowanej macierzy hesjanu.

Metoda zmiennej metryki charakteryzuje się szybszą zbieżnością (zbieżność superliniowa) niż metoda największego spadku (zbieżność liniowa). Poza tym fakt, że hesjan w każdym kroku spełnia warunek dodatniej określoności, daje pewność, że spełnienie warunku $\mathbf{g}=\mathbf{0}$ odpowiada rozwiązaniu problemu minimalizacji. Metoda ta należy obecnie do grupy najlepszych technik optymalizacji funkcji wielu zmiennych. Jej wadą jest znaczna złożoność obliczeniowa (konieczność wyznaczenia w każdym cyklu $n^2$ elementów hesjanu), a także duże wymagania w odniesieniu do pamięci przy przechowywaniu macierzy hesjanu, co w przypadku funkcji z dużą liczbą zmiennych może być pewnym problemem.

Najważniejszym ograniczeniem jest jednak pogarszanie się uwarunkowania problemu przy znacznym wzroście liczby wag. Wzrastający wymiar hesjanu pociąga za sobą rosnącą wartość współczynnika uwarunkowania $\text{cond}(\mathbf{H})$ i pogarszającą się dokładność obliczeń. Z tego względu stosuje się ją do niezbyt dużych sieci. W przypadku komputerów PC jej skuteczne działanie zostało sprawdzone dla sieci o liczbie wag rzędu kilkuset (poniżej tysiąca).

2.5.3 Algorytm Levenberga-Marquardta

Inną implementacją newtonowskiej strategii optymalizacji jest metoda Levenberga-Marquardta. W metodzie tej dokładną wartość hesjanu $\mathbf{H}$ we wzorze (2.14) zastępuje się wartością również aproksymowaną $\mathbf{G}$, ale określaną na podstawie informacji zawartej w gradiencie z uwzględnieniem czynnika regularyzacji. Dla wprowadzenia metody przyjęta będzie definicja funkcji celu w postaci odpowiadającej istnieniu jednego wzorca uczącego, co przy $M$ wyjściach sieci odpowiada

gdzie $e_i(\mathbf{x})=\left(y_i(\mathbf{x})-d_i\right)$ oznacza błąd na $i$-tym wyjściu sieci przy prezentacji wektora $\mathbf{x}$. Przy wprowadzeniu oznaczeń ($\mathbf{J}$ – macierz Jacobiego)

wektor gradientu $\mathbf{g}$ i aproksymowana macierz hesjanu $\mathbf{G}$ odpowiadające funkcji celu (2.16) określane są w postaci [14,46]

W równaniu tym współczynnik $v$ zwany parametrem regularyzacyjnym Levenberga-Marquardta jest wielkością skalarną, dążącą do zera w trakcie procesu optymalizacyjnego. Na starcie procesu uczenia, gdy aktualna wartość wektora wagowego $\mathbf{w}_k$  jest daleka od rozwiązania (duże wartości składowych wektora błędu $\mathbf{e}$), przyjmuje się wartość parametru $v$ bardzo dużą w porównaniu z największą wartością własną macierzy $\mathbf{J}^T \mathbf{J}$. W takim przypadku hesjan jest reprezentowany przede wszystkim przez czynnik regularyzacyjny $v \mathbf{1}$ a poszukiwanie kierunku odbywa się podobnie jak w metodzie największego spadku $\mathbf{p}_k=-\mathbf{g}_k / v$. W miarę zmniejszania się residuum błędu i zbliżania do rozwiązania, parametr $v$ ulega redukcji i pierwszy składnik sumy we wzorze (2.18) określającym hesjan odgrywa coraz większą rolę.

O skuteczności działania algorytmu decyduje odpowiedni dobór wartości $v$. W praktyce wartość ta podlega ciągłej adaptacji przy zastosowaniu metody Levenberga-Marquardta [43]. Przy aktualnym rozwiązaniu bliskim optymalnemu hesjan jest dostatecznie wiernie reprezentowany przez składnik  $\mathbf{J}^T \mathbf{J}$ wzoru aproksymacyjnego (2.18). W takim przypadku czynnik regularyzacji $v \mathbf{1}$  może być pominięty ($v=0$) i proces określania hesjanu opiera się bezpośrednio na aproksymacji pierwszego rzędu, a kierunek wyznaczony jest wzorem $\mathbf{p}_k=-\left[\mathbf{J}_k^T \mathbf{J}_k\right]^{-1} \mathbf{g}_k$. Algorytm Levenberga-Marquardta przechodzi wówczas w klasyczny algorytm Gaussa-Newtona [14], charakteryzujący się kwadratową zbieżnością do rozwiązania optymalnego. W praktyce algorytm Levenberga-Marquardta ma zbieżność super liniową i obok metody zmiennej metryki jest uważany za najlepszy algorytm uczący sieci MLP przy założeniu sieci o niezbyt dużej liczbie adaptowanych wag.

2.5.4 Algorytm gradientów sprzężonych

Metody newtonowskie bardzo dobrze sprawdzają się w zadaniach uczenia sieci MLP przy liczbie wag nie przekraczających kilkuset. Przy zbyt dużej liczbie zmiennych optymalizowanych pojawia się problem złego uwarunkowania macierzy $\mathbf{H}$. W takim przypadku lepsze wyniki uczenia uzyskuje się zwykle stosując metodę gradientów sprzężonych [14]. W metodzie tej podczas wyznaczania kierunku minimalizacji rezygnuje się z bezpośredniej informacji o hesjanie. W zamian za to kierunek poszukiwań $\mathbf{p}$ jest konstruowany w taki sposób, aby był ortogonalny oraz sprzężony ze wszystkimi poprzednimi kierunkami poszukiwań. W $k$-tej iteracji warunek ten jest spełniony przy doborze kierunku $\mathbf{p}$ według wzoru

Ze wzoru (2.19) wynika, że nowy kierunek minimalizacji zależy tylko od wartości gradientu w aktualnym punkcie rozwiązania oraz od poprzedniego kierunku poszukiwań $\mathbf{p}_{k-1}$, pomnożonego przez współczynnik sprzężenia $\beta$. Współczynnik sprzężenia odgrywa bardzo ważną rolę, kumulując w sobie informacje o poprzednich kierunkach poszukiwań i zapewniając ich ortogonalność względem kierunku aktualnego. Spośród wielu sposobów wyznaczania wartości tego współczynnika ograniczymy się do reguły Polaka-Ribiery [46]

Ze względu na kumulację błędów zaokrąglenia podczas kolejnych cykli obliczeniowych metoda gradientów sprzężonych w praktyce zatraca własność ortogonalności między wektorami kierunków minimalizacji. Dlatego też po wykonaniu $m$ iteracji ($m$ przyjmuje się jako funkcję liczby zmiennych $n$ podlegających optymalizacji) przeprowadza się jej restart, zakładając w pierwszym kroku po restarcie kierunek minimalizacji zgodnie z algorytmem największego spadku w aktualnie osiągniętym punkcie rozwiązania.

W praktyce obliczeniowej metoda gradientów sprzężonych wykazuje zbieżność super liniową, choć o mniejszej wartości wykładnika zbieżności niż metoda zmiennej metryki. Z tego powodu jest mniej skuteczna niż metody pseudo newtonowskie (zmiennej metryki i Levenberga-Marquardta), ale zdecydowanie szybsza niż metoda największego spadku. Stosuje się ją powszechnie jako skuteczny algorytm optymalizacji przy bardzo dużej liczbie zmiennych, sięgających nawet kilku tysięcy. Dzięki małym wymaganiom co do pamięci i niewielkiej złożoności obliczeniowej metoda ta pozwala na efektywne rozwiązanie stosunkowo dużych problemów optymalizacyjnych.

^1.