Rozmiar kroku opadającego gradientu adaptacyjnego, gdy nie można przeprowadzić wyszukiwania linii

Mam funkcję celu $E$ zależy od wartości $\phi(x, t = 1.0)$ , gdzie $\phi(x, t)$ jest rozwiązaniem dla PDE. Optymalizuję $E$ przez opadanie gradientu w początkowym stanie PDE: $\phi(x, t = 0.0)$ . To znaczy aktualizuję $\phi(x, t = 0.0)$ a następnie muszę zintegrować PDE, aby obliczyć resztę. Oznacza to, że gdybym szukał linii dla wielkości kroku spadku gradientu (nazwij to $\alpha$ ), dla każdej potencjalnej wartości $\alpha$ Musiałbym ponownie zintegrować PDE.

W moim przypadku byłoby to zbyt drogie. Czy istnieje inna opcja adaptacyjnego rozmiaru kroku spadku gradientu?

Nie szukam tu tylko schematów matematycznych (choć oczywiście lepiej, jeśli coś istnieje), ale byłbym zadowolony ze wszystkiego, co jest ogólnie lepsze niż statyczny rozmiar kroku.

Dzięki!

optimization pde conjugate-gradient

— NLi10Me
źródło

Nie sądzę, żebym chciał w tej chwili zmodyfikować sposób integracji PDE, ponieważ dla mnie byłoby to poważne przepisanie kodu. Poza tym PDE nie jest trudne, ponieważ muszę go rozwiązać na bardzo gęstej siatce w czasoprzestrzeni, ponieważ wymagam bardzo dużej dokładności numerycznej.

— NLi10Me

Z drugiej strony metoda BB (której nie znałem) wydaje się całkiem dobra; wszystko, co muszę zrobić, to śledzić stan i gradient poprzedniej iteracji i otrzymuję przybliżenie drugiego rzędu ... to wydaje się bardzo miłe. Jednak wyprowadzenie zakłada kwadrat wypukły, a mój problem prawie na pewno nie jest. Chociaż z pewnością znajduję (i jestem zadowolony) z minimów lokalnych zamiast globalnych. Czy wiesz, jak dobrze BB radził sobie z problemami o bardzo dużych wymiarach?

— NLi10Me

Chyba chodziło mi o lokalne minima, że w sąsiedztwie lokalnego minimum żadna funkcja nie jest w przybliżeniu kwadratowa? Myślę, że mój stan początkowy

ϕ^{(0)} (x, t = 0.0)

$\phi^{(0)}(x, t = 0.0)$ jest wystarczająco bliski minimum, ponieważ w wielu przypadkach uzyskuję płynną zbieżność nawet przy statycznym rozmiarze kroku. Tak więc, mimo że ma bardzo duże wymiary, i ogólnie biorąc pod uwagę całą przestrzeń poszukiwań, problem nie jest wypukły / niekwadratowy, czy BB nadal może być dobrym wyborem bez przeszukiwania linii?

— NLi10Me

Inne „składniki” do

E

$E$ są eksperymentalnymi danymi obrazu.

ϕ (x, t = 1.0)

$\phi(x, t = 1.0)$ próbuje dopasować jeden obraz, aby „dopasować” do drugiego (mierzony za pomocą funkcji dopasowania, takiej jak norma L2 zintegrowana z wokselami). W przypadku niektórych par obrazów uzyskuję płynną zbieżność z (moim obecnym wyborem) statycznym rozmiarem kroku. W przypadku innych par obrazów mam dużo oscylacji. System musi być w pełni zautomatyzowany, więc nie mogę cofnąć się i ręcznie edytować rozmiaru kroku dla kłopotliwych par obrazów.

— NLi10Me

Racja, muszę rozwiązać system przylegania, aby uzyskać gradient (który jest paskudniejszy i zajmuje więcej czasu). Ok, myślę, że spróbuję BB z wyszukiwaniem linii powrotnej. Dziękuję bardzo za radę; moi doradcy są często trudni do zdobycia i wielu z nich nie jest zainteresowanych wdrożeniem tak bardzo jak sam model. Uważam, że metody numeryczne są kluczowym elementem do wykazania, czy model jest dobry, czy nie, więc dzięki jeszcze raz bardzo to doceniam.

— NLi10Me

Zacznę od ogólnej uwagi: informacje pierwszego rzędu (tj. Użycie tylko gradientów, które kodują nachylenie) mogą dać tylko informacje kierunkowe: mogą powiedzieć, że wartość funkcji maleje w kierunku wyszukiwania, ale nie na jak długo . Aby zdecydować, jak daleko iść w kierunku wyszukiwania, potrzebujesz dodatkowych informacji (opadanie gradientu ze stałymi długościami kroków może się nie powieść nawet w przypadku wypukłych problemów kwadratowych). W tym celu masz zasadniczo dwie możliwości:

Użyj informacji drugiego rzędu (która koduje krzywiznę), na przykład stosując metodę Newtona zamiast spadku gradientu (dla którego zawsze możesz użyć długości kroku $1$ wystarczająco blisko minimizatora).
Próba i błąd (przez co oczywiście mam na myśli właściwe wyszukiwanie linii, takie jak Armijo).

Jeśli, jak piszesz, nie masz dostępu do drugich pochodnych, a ocena funkcji obejctive jest bardzo droga, jedyną nadzieją jest kompromis: użyj wystarczającej przybliżonej informacji drugiego rzędu, aby uzyskać odpowiednią długość kroku kandydata, na przykład linię szukaj tylko $\mathcal{O}(1)$ oceny (tj. co najwyżej (mała) stała wielokrotność wysiłku potrzebnego do oceny gradientu).

Jedną z możliwości jest zastosowanie długości kroku Barzilai - Borwein (patrz np. Fletcher: W metodzie Barzilai-Borwein . Optymalizacja i kontrola za pomocą aplikacji, 235–256, Appl. Optim., 96, Springer, New York, 2005 ). Chodzi o to, aby użyć przybliżenia skończonej różnicy krzywizny wzdłuż kierunku wyszukiwania, aby uzyskać oszacowanie wielkości kroku. W szczególności wybierz $\alpha_0>0$ dowolne, ustawione $g^0:=\nabla f(x^0)$ a następnie dla $k=0,...$ :

Zestaw $s^k = -\alpha_k^{-1} g^k$ i $x^{k+1}=x^k+s^k$
Oceniać $g^{k+1}=\nabla f(x^{k+1})$ i nastaw $y^k = g^{k+1}-g^{k}$
Zestaw $\alpha_{k+1} = \frac{(y^k)^Ty^k}{(y^k)^Ts^k}$

Ten wybór można wykazać jako zbieżny (w praktyce bardzo szybko) dla funkcji kwadratowych, ale zbieżność nie jest monotoniczna (tj. Wartość funkcji $f(x^{k+1})$ może być większy niż $f(x^k)$ , ale tylko raz na jakiś czas; patrz wykres na stronie 10 w pracy Fletchera). W przypadku funkcji niekwadratowych należy połączyć to z wyszukiwaniem linii, które należy zmodyfikować, aby poradzić sobie z niemonotonicznością. Jedną z możliwości jest wybór $\sigma_k \in (0,\alpha_k^{-1})$ (np. przez cofanie), takie jak

f (x^{k} - σ_{k} g^{k}) \leq max_{max (k - M, 1) \leq j \leq k} f (x^{j}) - γ σ_{k} (g^{k})^{T} g^{k},

$f(x^k - \sigma_k g^k) \leq \max_{\max(k-M,1)\leq j\leq k} f(x^j) - \gamma \sigma_k (g^k)^Tg^k,$ gdzie jest typowym parametrem Armijo, a kontroluje stopień monotoniczności (np. ). Istnieje również wariant, który używa wartości gradientu zamiast wartości funkcji, ale w twoim przypadku gradient jest nawet droższy do oceny niż funkcja, więc nie ma to sensu. (Uwaga: Możesz oczywiście spróbować ślepo zaakceptować długość kroku BB i zaufać swojemu szczęściu, ale jeśli potrzebujesz jakiejkolwiek solidności - jak napisałeś w komentarzach - byłby to naprawdę zły pomysł).

γ \in (0, 1)

$\gamma\in(0,1)$

M

$M$

M = 10

$M=10$

Alternatywnym (i moim zdaniem znacznie lepszym) podejściem byłoby zastosowanie tego przybliżenia różnic skończonych już w obliczeniach kierunku wyszukiwania; nazywa się to metodą quasi-Newtona . Chodzi o stopniowe budowanie przybliżenia Hesji przy użyciu różnic gradientów. Na przykład możesz wziąć (macierz tożsamości), a dla rozwiązują i ustaw pomocą jak wyżej i . (To się nazywa aktualizacja Broyden $\nabla^2 f(x^k)$ $H_0=\mathrm{Id}$ $k=0,\dots$

\begin{matrix} (1) & H_{k} s^{k} = - g^{k}, \end{matrix}

$H_{k}s^{k} = -g^{k},\label{cc1}\tag{1}$

H_{k + 1} = H_{k} + \frac{(y^{k} - H_{k} s^{k})^{T} (s^{k})^{T}}{(s^{k})^{T} s^{k}}

$H_{k+1} = H_k + \frac{(y^k-H_ks^k)^T(s^k)^T}{(s^k)^Ts^k}$

y^{k}

$y^k$

x^{k + 1} = x^{k} + s^{k}

$x^{k+1} = x^k +s^k$ i jest rzadko stosowany w praktyce; lepszą, ale nieco bardziej skomplikowaną aktualizacją jest aktualizacja BFGS , do której - i więcej informacji - odnoszę się do książki Nocedal i Wrighta pt . Optymalizacja numeryczna .) Minusem jest to, że a) wymagałoby to rozwiązania systemu liniowego na każdym etapie (ale tylko wielkości nieznanej, co w twoim przypadku jest warunkiem początkowym, dlatego wysiłek powinien być zdominowany przez rozwiązywanie PDE w celu uzyskania gradientu; istnieją też zasady aktualizacji przybliżeń odwrotnego Hesji, które wymagają obliczenia tylko jednej macierzy - produkt wektorowy) ib) nadal potrzebujesz wyszukiwania linii, aby zagwarantować zbieżność ...

Na szczęście w tym kontekście istnieje alternatywne podejście, które wykorzystuje każdą ocenę funkcji. Chodzi o to, że dla symetryczny i dodatni określony (co jest gwarantowane w przypadku aktualizacji BFGS) rozwiązanie jest równoważne zminimalizowaniu modelu kwadratowego W metodzie regionu zaufania zrobiłbyś to z dodatkowym ograniczeniem, że , gdzie jest odpowiednio wybranym promieniem regionu zaufania (który odgrywa rolę długości kroku ). Kluczową ideą jest teraz, aby wybrać ten promień adaptacyjnie, w oparciu o obliczony krok. W szczególności patrzysz na stosunek $H_k$ $\eqref{cc1}$

q_{k} (s) = \frac{1}{2} s^{T} H_{k} s + s^{T} g^{k} .

$q_k(s) = \frac12 s^T H_k s + s^T g^k.$

‖ s ‖ \leq Δ_{k}

$\|s\| \leq \Delta_k$

Δ_{k}

$\Delta_k$

σ_{k}

$\sigma_k$

ρ_{k} := \frac{f (x^{k}) - f (x^{k} + s^{k})}{f (x^{k}) - q_{k} (s^{k})}

$\rho_k := \frac{f(x^k)-f(x^k+s^k)}{f(x^k)-q_k(s^k)}$ rzeczywistej i przewidywanej redukcji wartości funkcji. Jeśli jest bardzo mały, twój model był zły, a ty i próbujesz ponownie z . Jeśli jest bliskie , twój model jest dobry i ustawiasz i zwiększasz . W przeciwnym razie po prostu ustaw i pozostaw spokoju. Aby obliczyć rzeczywisty minimalizator z

ρ_{k}

$\rho_k$

s^{k}

$s^k$

Δ_{k + 1} < Δ_{k}

$\Delta_{k+1}<\Delta_k$

ρ_{k}

$\rho_k$

1

$1$

x^{k + 1} = x^{k} + s^{k}

$x^{k+1}=x^k+s^k$

Δ_{k + 1} > Δ_{k}

$\Delta_{k+1}>\Delta_k$

x^{k + 1} = x^{k} + s^{k}

$x^{k+1}=x^k+s^k$

Δ_{k}

$\Delta_k$

s^{k}

$s^k$

min_{‖ s ‖ \leq Δ_{k}} q_{k} (s)

$\min_{\|s\|\leq \Delta_k} q_k(s)$ , istnieje kilka strategii pozwalających uniknąć konieczności rozwiązania pełnego problemu optymalizacji; moją ulubioną jest skrócona metoda CG Steihauga . Aby uzyskać więcej informacji, ponownie odnoszę się do Nocedal i Wright.

— Christian Clason
źródło

Właśnie teraz patrzę na to ponownie i zdaję sobie sprawę, że mam pytanie. W kroku trzecim dla metody BB masz ; gdzie i . Licznik i mianownik w wyrażeniu dla wyglądają jak produkty wewnętrzne. W moim przypadku , gdzie jest przestrzenią wektorową z nietrywialną metryką Riemanniana: K. To znaczy, . Czy to wpływa na definicję ?

α_{k + 1} = \frac{(y^{k})^{T} y^{k}}{(y^{k})^{T} s^{k}}

$\alpha_{k+1} = \frac{(y^k)^Ty^k}{(y^k)^Ts^k}$

y^{k} = g^{k + 1} - g^{k}

$y^{k} = g^{k+1} - g^k$

s^{k} = - α_{k}^{- 1} g^{k}

$s^k = -\alpha_k^{-1}g^k$

α_{k + 1}

$\alpha_{k+1}$

g^{k} \in V^{*}

$g^k \in V^*$

V^{*}

$V^*$

⟨ g^{k}, g^{k} ⟩_{V^{*}} = ⟨ g^{k}, K g^{k} ⟩_{L_{2}}

$\langle g^k, g^k \rangle _{V^*} = \langle g^k, Kg^k \rangle_{L_2}$

α_{k + 1}

$\alpha_{k+1}$

— NLi10Me

Tak, jeśli masz nietrywialną strukturę przestrzeni wektorowej, powinieneś ją uszanować w algorytmach. W szczególności należy rozróżnić między iloczynami wewnętrznymi dwóch funkcji w tej samej przestrzeni (np. i ) i iloczynami dualności między funkcją w przestrzeni a jedną w przestrzeni podwójnej (np. i ) - w tym drugim przypadku musisz dołączyć mapowanie Riesz, aby najpierw przekształcić go w produkt wewnętrzny. (Można to interpretować jako warunek wstępny.)

y^{k}

$y^k$

y^{k}

$y^k$

s^{k}

$s^k$

y^{k}

$y^k$

— Christian Clason,

Dr Clason, prześlę artykuł na ISBI 2017, szczegółowo opisujący niektóre eksperymenty, które przeprowadziłem przy użyciu metody wyszukiwania linii BB + dla zadania rejestrowania obrazu diffeomorficznego. Czy chciałbyś zostać uwzględniony jako autor manuskryptu? Nie napisałem tego jeszcze, ale większość eksperymentów mam ukończonych lub trwających. Proszę daj mi znać.

— NLi10Me

@ NLi10Me Dziękuję za uprzejmą ofertę, ale nie zrobiłem nic, co zasługiwałoby na współautorstwo - wszystko, co napisałem, to standardowy podręcznik. Jeśli czujesz do tego silne zdanie, możesz podziękować za „pomocne uwagi na temat (cokolwiek to pomogło)”, ale nawet to nie byłoby wymagane. Wystarczy wiedzieć, że to, co napisałem, było pomocne!

— Christian Clason

Przepraszamy, masz rację, to literówka - naprawione! (Warunek Armijo jest często zapisywany jako , gdzie jest kierunkiem wyszukiwania - niekoniecznie ujemnym gradient - i rozmiar kroku, co powinno wyjaśnić, co się dzieje.)

f (x + σ s) - f (x) \leq γ \nabla f (x)^{T} (σ s)

$f(x+\sigma s) - f(x) \leq \gamma\nabla f(x)^T(\sigma s)$

s

$s$

σ

$\sigma$

— Christian Clason,