Dlaczego w niektórych przypadkach dynamika hamiltonowska jest lepsza niż propozycja chodzenia losowego w MCMC?

W niektórych przypadkach dynamika hamiltonowska zawsze przewyższa osiągi w algorytmie Metropolis. Czy ktoś mógłby wyjaśnić przyczynę prostymi słowami bez zbyt dużej matematyki?

Po pierwsze, pozwól mi stwierdzić, że nie sądzę, aby współczynnik akceptacji dla HMC (Hamiltonian Monte Carlo) był zawsze wyższy niż dla algorytmu Metropolis. Jak zauważył @JuhoKokkala, współczynnik akceptacji Metropolis jest dostrajany, a wysoki wskaźnik akceptacji nie oznacza, że ​​Twój algorytm dobrze sprawdza się w badaniu rozkładu tylnego. Jeśli użyjesz bardzo wąskiej dystrybucji propozycji (na przykład$T(q|q')=\mathcal{N}(q',\sigma I)$ z bardzo małym $\sigma$), uzyskasz bardzo wysoki wskaźnik akceptacji, ale tylko dlatego, że zasadniczo zawsze pozostajesz w tym samym miejscu, bez badania pełnego rozkładu z tyłu.

To, co myślę, że naprawdę pytasz (a jeśli mam rację, to proszę odpowiednio zmodyfikuj swoje pytanie), dlaczego Hamiltonian Monte Carlo ma (w niektórych przypadkach) lepszą wydajność niż Metropolis. „Lepsza wydajność” oznacza, że ​​dla wielu aplikacji, jeśli porównasz łańcuch wygenerowany przez HMC z łańcuchem o tej samej długości (ta sama liczba próbek ) wygenerowanym przez algorytm Metropolis, łańcuch HMC osiągnie stan ustalony wcześniej niż Łańcuch Metropolis znajduje niższą wartość ujemnego prawdopodobieństwa logarytmu (lub podobną wartość, ale w mniejszej liczbie iteracji), efektywny rozmiar próbki jest mniejszy, autokorelacja próbek rozpada się szybciej z opóźnieniem itp.$N$

Spróbuję wyjaśnić, dlaczego tak się dzieje, nie wnikając zbytnio w matematyczne szczegóły. Przypomnijmy więc przede wszystkim, że algorytmy MCMC są ogólnie przydatne do obliczania całek wielowymiarowych (oczekiwań) funkcji (lub większej liczby funkcji) w odniesieniu do gęstości docelowej , zwłaszcza gdy nie mamy sposób bezpośredniego pobrania próbki z docelowej gęstości:$f$$\pi(q)$

$\mathbb{E}_{\pi}{[f]}=\int_{\mathcal{Q}} f(q)\pi(q)\text{d}q_1\dots\text{d}q_d$

gdzie $q$ jest wektorem $d$ parametry, na których $f$ i $\pi$ zależeć i $\mathcal{Q}$to przestrzeń parametrów. Teraz, w wysokich wymiarach, objętość przestrzeni parametrów, która najbardziej przyczynia się do powyższej całki, nie jest sąsiedztwem trybu$\pi(q)$ (tj. nie jest to wąska objętość wokół oszacowania MLE na $q$), ponieważ tutaj $\pi(q)$ jest duży, ale głośność jest bardzo mała.

Załóżmy na przykład, że chcesz obliczyć średnią odległość punktu $q$ od pochodzenia $\mathbb{R}^d$, gdy jego współrzędne są niezależnymi zmiennymi Gaussa o średniej zerowej i wariancji jednostkowej. Następnie powyższa całka staje się:

$\mathbb{E}_{\pi}{[X]}=\int_{\mathcal{Q}} ||q||(2\pi)^{-d/2}\exp{(-||q||^2/2)}\text{d}q_1\dots\text{d}q_d$

Teraz gęstość docelowa $\pi(q)=(2\pi)^{-d/2}\exp{(-||q||^2/2)}$ ma oczywiście maksimum na 0. Jednak przez zmianę na współrzędne sferyczne i wprowadzenie $r=||q||$, widać, że całka staje się proporcjonalna do $r^{d-1}\exp{(-r^2/2)} \text{d}r$. Ta funkcja ma oczywiście maksimum w pewnej odległości od źródła. Region w środku$\mathcal{Q}$który w największym stopniu przyczynia się do wartości całki nazywany jest zestawem typowym , a dla tej całki zestaw typowy jest kulistą powłoką o promieniu$R\propto\sqrt{d}$.

Teraz można pokazać, że w idealnych warunkach łańcuch Markowa wygenerowany przez MCMC najpierw zbiega się do punktu w typowym zestawie, a następnie zaczyna eksplorować cały zestaw, a na koniec kontynuuje eksplorację jego szczegółów. W ten sposób oszacowanie MCMC oczekiwań staje się coraz bardziej dokładne, z tendencyjnością i wariancją, która zmniejsza się wraz ze wzrostem liczby kroków.

Jednak, gdy geometria typowego zestawu jest skomplikowana (na przykład, jeśli ma guzek w dwóch wymiarach), wówczas standardowy algorytm Metropolis z losowym przejściem ma wiele trudności w badaniu „patologicznych” szczegółów zestawu. Ma tendencję do przypadkowego przeskakiwania „wokół” tych regionów, bez ich eksploracji. W praktyce oznacza to, że oszacowana wartość całki ma tendencję do oscylowania wokół prawidłowej wartości, a przerwanie łańcucha przy skończonej liczbie kroków spowoduje źle tendencyjne oszacowanie.

Hamiltonian Monte Carlo próbuje przezwyciężyć ten problem, wykorzystując informacje zawarte w rozkładzie docelowym (w jego gradiencie) w celu poinformowania propozycji nowego punktu próbnego, zamiast po prostu użyć rozkładu propozycji niezwiązanego z docelowym. Dlatego mówimy, że HMC wykorzystuje pochodne rozkładu docelowego do bardziej efektywnego badania przestrzeni parametrów. Sam gradient rozkładu docelowego nie jest jednak wystarczający do poinformowania o etapie wniosku. Jak w przykładzie średniej odległości losowego punktu od początku$\mathbb{R}^d$gradient gradientu docelowego sam w sobie kieruje nas do trybu rozkładu, ale obszar wokół tego trybu niekoniecznie jest regionem, który najbardziej przyczynia się do całki powyżej, tj. nie jest to typowy zbiór.

Aby uzyskać właściwy kierunek, w konsoli HMC wprowadzamy pomocniczy zestaw zmiennych, zwany zmiennymi pędu . Fizyczny analog może tu pomóc. Satelita krążący wokół planety pozostanie na stabilnej orbicie tylko wtedy, gdy jej pęd ma wartość „właściwą”, w przeciwnym razie albo odpłynie na otwartą przestrzeń, albo zostanie przyciągnięty do planety przez przyciąganie grawitacyjne (tutaj odgrywa rolę gradientu docelowej gęstości, który „przyciąga” do trybu). W ten sam sposób parametry pędu mają na celu utrzymanie nowych próbek w typowym zestawie, zamiast dryfowania w kierunku ogonów lub w kierunku trybu.

To jest krótkie streszczenie bardzo interesującego artykułu Michaela Betancourta na temat wyjaśniania Hamiltonian Monte Carlo bez nadmiernej matematyki. Papier, który jest bardziej szczegółowy, znajduje się tutaj .

Jedną z rzeczy, której artykuł nie opisuje wystarczająco szczegółowo, IMO, jest to, kiedy i dlaczego HMC może zrobić coś gorszego niż przypadkowa metropolia. Nie zdarza się to często (z mojego ograniczonego doświadczenia), ale może się zdarzyć. W końcu wprowadzasz gradienty, które pomagają znaleźć drogę w przestrzeni parametrów wielowymiarowych, ale podwajasz także wymiar problemu. Teoretycznie może się zdarzyć, że spowolnienie spowodowane wzrostem wymiarowości przezwycięży przyspieszenie wynikające z wykorzystania gradientów. Również (i jest to omówione w artykule), jeśli typowy zestaw ma regiony o wysokiej krzywiźnie, HMC może „przeregulować”, tj. Może rozpocząć próbkowanie bezużytecznych punktów bardzo daleko w ogonach, co nie ma nic wspólnego z oczekiwaniami. Jednak, powoduje to niestabilność integratora symplektycznego, który jest stosowany w praktyce do implementacji numerycznej HMC. Tak więc tego rodzaju problem można łatwo zdiagnozować.

Jak wspomniano w komentarzach @JuhoKokkala, wysoki wskaźnik akceptacji niekoniecznie daje dobre wyniki. Wskaźnik akceptacji Metropolis Hastings można zwiększyć, zmniejszając rozkład propozycji. Spowoduje to jednak podjęcie mniejszych kroków, co sprawi, że zbadanie dystrybucji docelowej zajmie więcej czasu. W praktyce istnieje kompromis między wielkością kroku a stopą akceptacji, a do uzyskania dobrej wydajności potrzebna jest odpowiednia równowaga.

Hamiltonian Monte Carlo ma tendencję do przewyższania Metropolis Hastings, ponieważ może dotrzeć do bardziej odległych punktów z większym prawdopodobieństwem akceptacji. Pytanie zatem brzmi: dlaczego HMC ma większe prawdopodobieństwo akceptacji niż MH dla bardziej odległych punktów ?

MH ma problem z dotarciem do odległych punktów, ponieważ jego propozycje są składane bez wykorzystywania informacji o rozkładzie docelowym. Rozkład propozycji jest zazwyczaj izotropowy (np. Symetryczny gaussowski). Tak więc w każdym punkcie algorytm próbuje przesunąć losową odległość w losowym kierunku. Jeśli odległość jest niewielka w stosunku do tego, jak szybko rozkład celu zmienia się w tym kierunku, istnieje duża szansa, że ​​gęstość w obecnych i nowych punktach będzie podobna, dając co najmniej rozsądną szansę na akceptację. Na większych odległościach rozkład docelowy mógł się nieco zmienić w stosunku do bieżącego punktu. Tak więc szansa na losowe znalezienie punktu o podobnej lub (miejmy nadzieję) wyższej gęstości może być niska, szczególnie gdy wzrasta wymiarowość. Na przykład jeśli bieżący punkt leży na wąskim grzbiecie, wówczas „

Natomiast konsola HMC wykorzystuje strukturę dystrybucji docelowej. Jego mechanizm propozycji można pomyśleć o zastosowaniu fizycznej analogii, jak opisano w Neal (2012). Wyobraź sobie krążek ślizgający się po pagórkowatej, pozbawionej tarcia powierzchni. Lokalizacja krążka reprezentuje bieżący punkt, a wysokość powierzchni reprezentuje log ujemny rozkładu celu. Aby uzyskać nowy proponowany punkt, krążek otrzymuje pęd z losowym kierunkiem i wielkością, a jego dynamika jest następnie symulowana podczas przesuwania się po powierzchni. Krążek będzie przyspieszał w kierunku z góry i zwalniał w kierunku pod górę (być może nawet zatrzyma się i ześlizgnie z powrotem w dół). Trajektorie poruszające się w bok wzdłuż ściany doliny będą zakręcać w dół. Sam krajobraz wpływa więc na trajektorię i ciągnie ją w kierunku regionów o wyższym prawdopodobieństwie. Pęd może pozwolić, aby krążek przeszedł przez małe wzgórza, a także przekroczył małe baseny. Lokalizacja krążka po kilku krokach czasowych daje nowy proponowany punkt, który jest akceptowany lub odrzucany przy użyciu standardowej reguły Metropolis. Wykorzystanie rozkładu docelowego (i jego gradientu) pozwala konsoli HMC dotrzeć do odległych punktów przy wysokich wskaźnikach akceptacji.

Oto dobra recenzja:

Neal (2012) . MCMC z wykorzystaniem dynamiki hamiltonowskiej.

Jako luźna odpowiedź (która wydaje się być tym, czego szukasz) metody hamiltonowskie uwzględniają pochodną prawdopodobieństwa logarytmicznego, podczas gdy standardowy algorytm MH nie.