Podręcznik
2. Gradientowe algorytmy rozwiązywania zadań optymalizacji bez ograniczeń
2.2. Algorytm gradientu sprzężonego
Zajmiemy się najpierw przypadkiem pierwszym – poprawką addytywną.
Przyjmijmy na moment, że chcemy minimalizować metodą gradientu prostego funkcję kwadratową o dodatnio określonej macierzy Q:
Rysunek zrobimy dla funkcji dwu zmiennych.
Algorytm wystartował z punktu i w kierunku antygradientu dokładnie minimalizując w kierunku znalazł punkt . Następny krok zaprowadzi go do punktu leżącego na prostej wyznaczonej przez antygradient , i tak dalej. Zauważmy, że tak wyznaczane kolejne kierunki poprawy (antygradienty) są do siebie prostopadłe
nie tylko dla funkcji kwadratowych, ale dla dowolnych funkcji różniczkowalnych, ponieważ kolejne punkty są punktami minimalizującymi funkcję
co oznacza, że spełniony jest warunek (różniczkowanie funkcji złożonej)
Rysunek pokazuje, że w punkcie istnieje kierunek lepszy od antygradientu, oznaczmy go q, mający tą własność, że dokładna minimalizacja na tym kierunku doprowadzi do rozwiązania . Wyznaczyć ten kierunek jest łatwo. Przy naszych założeniach rozwiązanie zadania optymalizacji musi spełniać równanie (Twierdzenie 3.6):
Zatem . Jak na razie, te obliczenia nic ciekawego nie dają, ale policzmy czemu równa się iloczyn . Najpierw zauważmy, że
Teraz możemy policzyć wartość interesującego nas iloczynu
czyli
Wprowadźmy następujące definicje:
Niezerowe, wektory oraz takie, że istnieje dodatnio określona macierz , dla której
nazywamy sprzężonymi względem tej macierzy.
elementowy zbiór niezerowych wektorów przestrzeni n-wymiarowej, , nazywamy układem wektorów wzajemnie sprzężonych względem dodatnio określonej macierzy Q (układem Q-sprzężonym) jeżeli
Prawdzie jest następujące stwierdzenie
Układ r wektorów Q-sprzężonych jest liniowo niezależny
Najprostszym przykładem dwu wektorów sprzężonych są wektory prostopadłe, które są sprzężone względem macierzy jednostkowej. Trzy wektory w przestrzeni trójwymiarowej są sprzężone względem macierzy jednostkowej jeżeli każdy z nich jest prostopadły do płaszczyzny wyznaczonej przez dwa pozostałe. Zatem już w przestrzeni trójwymiarowej sprzężenie względem macierzy jednostkowej to inna własność niż prostopadłość.
Wzór (4.3): , pokazuje, że wektor antygradientu obliczony w punkcie początkowym i wektor ustalający kierunek na którym leży punkt minimalizujący są sprzężone względem macierzy określającej minimalizowaną funkcję kwadratową.
Wektor jest dany wzorem , zatem aby go obliczyć trzeba odwrócić macierz .
Zastanowimy się teraz nad odpowiedzią na pytanie: czy można go wyznaczyć bez odwracania macierzy?
Przyjmiemy zatem, że
gdzie współczynnik trzeba wyznaczyć. Do jego wyliczenia posłużmy się warunkiem Q-sprzężoności wektorów i czyli równaniem . Zauważmy najpierw, że
a także
Mamy zatem
co daje
Czyli
Wzory (4.4) i (4.5) pozwalają więc wyznaczyć dla zadania optymalizacji o kwadratowej funkcji celu z macierzą Q, posługując się tylko gradientami, kierunek sprzężony względem Q z początkowym kierunkiem antygradientu.
Dla dwuwymiarowego przykładu tak określony kierunek „zaprowadzi” do rozwiązania (jest wektorem q). Jednak w przypadku wielowymiarowym tak być nie musi, bo wektorów sprzężonych do antygradientu jest nieskończenie wiele (jak wiemy, już w przestrzeni trójwymiarowej, jest nieskończenie wiele wektorów prostopadłych do danego wektora – leżą one na płaszczyźnie prostopadłej do tego wektora) a tylko jeden „prowadzi” do punktu minimalizującego. Nie znaleziono metody pozwalającej od razu wybrać ten właściwy wektor, udowodniono natomiast twierdzenie następujące.
Dla danej funkcji kwadratowej o dodatnio określonej macierzy (a więc kwadratowej funkcji ścisłe wypukłej) i dowolnego punktu początkowego określamy ciąg n punktów następująco:
gdzie
a) Określony powyżej n-elementowy zbiór wektorów jest zbiorem kierunków poprawy, a ponadto jest układem wektorów Q-sprzężonych.
b) , co przy poczynionych założeniach, że funkcja f jest kwadratowa i macierz Q jest dodatnio określona, oznacza (Twierdzenie 3.6) że .
Teza b) powyższego twierdzenia stwierdza więc, że za pomocą opisanego algorytmu można znaleźć rozwiązanie zadania minimalizacji funkcji kwadratowej w co najwyżej n krokach, inaczej mówiąc – dla funkcji kwadratowej algorytm ma skończoną zbieżność (pod warunkiem, że minimalizacja w kierunku jest dokładna, wzór (4.6)). Napisaliśmy „w co najwyżej...”, ponieważ teza a) twierdzenia mówi o tym, że kolejne kierunki są sprzężone z pierwszym – wektorem antygradientu, zatem zgodnie z naszymi wcześniejszymi rozważaniami kierunek prowadzący do minimum może zostać wyznaczony szybciej niż dopiero w n-tym kroku. Zauważmy też, że wzór (4.8) określający współczynnik ustalający „wkład” poprzedniego kierunku do nowego kierunku poprawy jest uogólnieniem wyprowadzonego poprzednio wzoru (4.5).
W twierdzeniu 4.2 został sformułowany algorytm o skończonej zbieżności pozwalający znaleźć rozwiązanie zadania minimalizacji ściśle wypukłej funkcji kwadratowej. Około 1963 R. Fletcher i C.M. Reeves zastosowali ten algorytm do minimalizacji funkcji wypukłej. Byli niepoprawnymi optymistami, bo z dowodu twierdzenia wynika, że skończona zbieżność algorytmu jest konsekwencją tego, że kolejne kierunki powiększają zbiór wektorów Q-sprzężonych, a przecież dla nie kwadratowej funkcji wypukłej nie można mówić o kierunkach sprzężonych względem macierzy Q, bo nie ma macierzy Q! Co więcej ta „kolejna sprzężoność” może nie być zachowana nawet dla funkcji kwadratowej gdy minimalizacja w kierunku nie jest dokładna. Ale udało się – algorytm zachowywał się nadspodziewanie dobrze, potrafił na przykład bez zygzakowania, w rozsądnej liczbie kroków przy interpolacyjnej metodzie minimalizacji w kierunku znaleźć rozwiązanie dla bananowej funkcji Rosenbrocka. W artykule opublikowanym w 1964 r. Fletcher i Reeves nazwali swój algorytm Conjugate Gradient Method – metoda gradientu sprzężonego, i tak zostało, pomimo tego że kierunki generowane przez ich algorytm, choć oparte o gradient (patrz wzór (4.7)) dla funkcji nie kwadratowych nie są sprzężone. Intensywne badania metody gradientu sprzężonego prowadzone w drugiej połowie lat sześćdziesiątych doprowadziły do opublikowania w 1969 r. następującego algorytmu, w którym współczynnik jest wyliczany według innej zależności, zwanej wzorem Poljaka, Polaka, Ribière’a
Algorytm gradientu sprzężonego (wersja Poljaka, Polaka, Ribière’a z odnową)
Inicjalizacja
Wybierz punkt początkowy , podstaw .
Kroki algorytmu
- (Odnowa) Oblicz , podstaw , jako kierunek poprawy wybierz kierunek antygradientu podstawiając . Podstaw , idź do 4.
- Podstaw , (stosujemy wzór Poljaka, Polaka, Ribière’a (4.9)). Jako kandydata na kierunek poprawy wybierz wektor , podstawiając (wzór (4.7)).
- Jeżeli ( jest kierunkiem poprawy), podstaw oraz i przejdź do 4. W przeciwnym przypadku podstaw i przejdź do 1 (dokonaj odnowy).
- Wyznacz długość kroku w kierunku poprawy podstaw.
- Wylicz następny punkt podstawiając .
- Oblicz , oraz , podstaw , podstaw .
- Jeżeli , to przyjmij za rozwiązanie, podstaw i stop.
- W przeciwnym przypadku idź do 7.
- Podstaw podstaw Jeżeli , gdzie jest wymiarem przestrzeni poszukiwań, to idź do 2. W przeciwnym przypadku podstaw i idź do 1 (odnowa po krokach od ostatniej odnowy).
Przeanalizujmy istotne aspekty działania przedstawionego wyżej algorytmu gradientu sprzężonego.
Jeżeli w kroku 4 algorytmu wykonamy dokładną minimalizację w kierunku, to kierunki wyznaczone w kroku 2 będą kierunkami poprawy, co poniżej pokażemy.
Jak pamiętamy, przy naszych założeniach, kierunek jest kierunkiem poprawy gdy . Policzmy zatem ten iloczyn
Pokażemy, że . Przy dokładnej minimalizacji w kierunku, dla punktu zachodzą równości
co oznacza, że spełniony jest warunek
Zatem z równości (4.10) wynika, że
Oznacza to, że przy dokładnej minimalizacji w kierunku krok 3, sprawdzający czy kierunki wyznaczone w kroku 2 są kierunkami poprawy, nie jest potrzebny.
Wiemy jednak, że każda praktyczna realizacja minimalizacji w kierunku jest obarczona pewnym błędem. Zatem konstruując algorytm, który ma naprawdę służyć do rozwiązywania zadań optymalizacji musimy pamiętać o stosownych zabezpieczeniach pozwalających pokonać trudności związane z różnego rodzaju niedokładnościami obliczeniowymi. Warunek kroku 3 sprawdza czy wyliczony kierunek jest kierunkiem poprawy. Jeżeli nie jest, dokonuje się tak zwanej odnowy – jako kierunek poszukiwań wybierany jest antygradient, który (gdy nie jest zerowy) jest zawsze kierunkiem poprawy. Jeżeli jednak zbyt często będziemy zmuszeni stosować odnowę, to działanie algorytmu będzie niewiele różnić się od działania algorytmu największego spadku – w szczególności wystąpi „zygzakowanie”. Dlatego przy implementacjach metody gradientu sprzężonego do znajdowania długości kroku w kierunku najczęściej wybiera się jedną z metod dających dobrą aproksymację minimum w kierunku.
- Analizując zachowanie algorytmu największego spadku (gradientu prostego) upatrywaliśmy przyczyn jego złego zachowania – zygzakowania – w braku pamięci. Wzór (4.7) określający nowy kierunek poprawy jest realizacją pierwszego z zaproponowanych podejść mających temu zaradzić – dodaniu poprawki zależnej od historii, ponieważ wiąże nowy kierunek poprawy z poprzednim kierunkiem poprawy, a także poprzez współczynnik (k) wyliczany wg wzoru (4.9) z gradientem obliczonym w kroku poprzednim. Ponieważ w kroku 4 nie możemy zapewnić dokładnej minimalizacji w kierunku, to okazało się, że kierunki wyznaczone w ten sposób mogą „prowadzić na manowce” – w tę stronę gdzie funkcja celu maleje powoli. Wobec tego aby zabezpieczyć się przed taką sytuacją, lub ją przerwać („oczyścić” pamięć) gdy już wystąpiła, przyjęto że odnowa następuje gdy od poprzedniej „upłynęło” n kroków. Wybór długości okresu, po którym następuje odnowa jako równego ilości zmiennych jest w zasadzie wyborem heurystycznym i odnosi się do zadań o dużej wymiarowości (kilkanaście zmiennych). Dla zadań o kilku zmiennych zręczniej jest dokonywać odnowy co 3n lub co 2n kroków, ponieważ z odnową co n kroków działanie algorytmu dla takich zadań może różnić się niewiele od działania algorytmu największego spadku.
- Zauważmy, że przy dokładnej minimalizacji w kierunku wzory Fletchera – Reevesa (4.8) i Poljaka – Polaka – Ribière’a (4.9) dają tą samą wartość współczynnika , ponieważ wtedy
a więc odjemnik w liczniku wzoru (4.9) jest równy zeru:
Dlaczego wybraliśmy wersję Poljaka – Polaka – Ribière’a (Rosjanina B.T. Poljaka, Amerykanina E. Polaka i Francuza G. Ribière) (4.9) a nie Fletchera – Reevesa (4.8)?
Ponieważ testy różnych wersji wzoru określającego wpływ poprzedniego kierunku poprawy na kierunek bieżący (zaproponowano ich więcej niż tylko przedstawione dwa wyrażenia) pokazały, że algorytm posługujący się wzorem Poljaka – Polaka – Ribière’a częściej znajduje rozwiązanie w mniejszej liczbie iteracji aniżeli algorytm ze wzorem Fletchera – Reevesa. Ponadto częściej potrafi znaleźć rozwiązanie w sensownej liczbie iteracji przy obniżonej dokładności minimalizacji w kierunku – jest więc bardziej odporny błędy wyznaczania kolejnych punktów. Formalnego wyjaśnienia tego faktu nie ma.
M.J.D. Powell przedstawił następującą hipotezę. Ponieważ minimalizowana funkcja nie jest kwadratowa i minimalizacja w kierunku nie jest dokładna, to przy minimalizacji wielu funkcji celu mamy do czynienia z sytuacją, w której kolejne kroki są coraz krótsze – leży coraz bliżej , algorytm rusza się coraz wolniej, jak się często mówi – grzęźnie. Bliskość kolejnych punktów oznacza, że gradienty w nich wyliczane są prawie takie same: .Dla formuły Fletchera – Reevesa (4.8) mamy więc Natomiast dla formuły Poljaka – Polaka – Ribière’a (4.9) dostajemy , co daje Algorytm ma zdolność do samo-odnowy w sytuacji ugrzęźnięcia. Takiej własności dla implementacji algorytmu ze wzorem Fletchera – Reevesa nie zaobserwowano.
Ponieważ odnowa – ruch w kierunku antygradientu – następuje wtedy gdy wyznaczony według wzoru (4.7) kierunek przestaje być kierunkiem poprawy, a także gdy od ostatniej odnowy „upłynęło” n kroków, to w nieskończonym ciągu generowanym przez algorytm zawiera się nieskończony podciąg punktów wyznaczonych na kierunkach określonych przez stosowne antygradienty . Zatem twierdzenie określające warunki dające zbieżność algorytmu gradientu sprzężonego jest podobne do twierdzenia dla algorytmu gradientu prostego.
Jeżeli
funkcja f jest różniczkowalna oraz ograniczona od dołu a jej gradient spełnia warunek Lipschitza na ,
to
nieskończony ciąg generowany przez algorytm gradientu sprzężonego, w którym długość kroku w kierunku antygradientu . jest ustalona:
– przez rozwiązanie zadania minimalizacji w kierunku:
– albo przez algorytm oparty na regule Armijo,
jest taki, że dla związanego z nim ciągu normy gradientu prawdziwa jest równość
.
Przy odpowiednio mocniejszych założeniach o funkcji wyboru oraz wymaganiu dokładnej minimalizacji w kierunku, można pokazać, że algorytm ten zbiega superliniowo do punktu minimalizującego.
Powiedzieliśmy, że algorytm gradientu sprzężonego jest lepszy od algorytmu największego spadku. Aby nie być gołosłownym musimy teraz to pokazać.
Wyniki testów muszą być porównywalne, rozważmy więc tą samą funkcję dwu zmiennych co w Przykładzie 8.2:
Algorytm wykonał siedem kroków i zatrzymał się w punkcie w którym , wartość funkcji celu , a odległość proponowanego rozwiązania od punktu optymalnego jest równa (w Przykładzie 4.1 była równa 0.3).
Krok pierwszy z do był oczywiście taki sam jak dla metody gradientu prostego. Przejście z kroku szóstego do siódmego jest na rysunku niezauważalne, ponieważ a . Zadanie jest dwuwymiarowe, zatem odnowa (ruch w kierunku antygradientu) nastąpiła w kroku trzecim, piątym i ostatnim – siódmym. Ponieważ punkty wyznaczone w drugim kroku przez oba algorytmy leżą jeszcze blisko siebie to i antygradienty są prawie równe, więc ruch obu algorytmów w trzecim kroku musi dać punkty , które też będą leżały blisko siebie. Lecz antygradient policzony w punkcie wyliczonym przez algorytm najszybszego spadku jest równy , a kierunek poprawy wyznaczony w swoim punkcie przez algorytm gradientu sprzężonego jest równy . Ta różnica spowodowała, że algorytm gradientu sprężonego wyznaczył jako punkt punkt o współrzędnych , co dało wartość funkcji celu , a algorytm gradientu prostego – punkt o współrzędnych z wartością funkcji celu . W piątym kroku algorytm gradientu sprzężonego wyznaczył punkt dający wartość funkcji celu i wartość normy gradientu , a więc punkt lepszy niż ostatni (ósmy) punkt uzyskany przez algorytm najszybszego spadku (przypominamy był to punkt w którym , oraz .
Algorytm gradientu sprzężonego (z za częstą odnową) zatrzymał się w punkcie oddalonym od rozwiązania o nieco ponad 0.02 z wartością funkcji celu większą od optymalnej o 0.0012 (algorytm prostego gradientu osiągnął, odległość – 0.09 i wartość większa od minimalnej o 0.0062). Okazał się wiec lepszy od współzawodnika, ale też nie bardzo poradził sobie z praktyczną utratą ścisłej wypukłości przez funkcję celu.
A jak algorytm gradientu sprzężonego poradził sobie ze znalezieniem minimum funkcji Rosenbrocka ?
Rysunek pokazuje kolejne kroki algorytmu gradientu sprzężonego bez odnowy (bo funkcja Rosenbrocka jest dwuwymiarowa) i z minimalizacją w kierunku wykorzystującą wielomian interpolacyjny. Widać, że Algorytm wykonał 24 iteracje poruszając się po dnie doliny bez zygzakowania (8 iteracji), dobrze przebył zakręt (8 iteracji) a następnie dość szybko „pomaszerował” w kierunku punktu minimalizującego (8 iteracji). Zastosowana modyfikacja gradientu w tym teście się sprawdziła. Widać jednak, że dla twórców algorytmów trudniej jest nauczyć Algorytm jak ma sobie poradzić z „rozpłaszczeniem się rzeźby terenu” niż z zaproponowaniem efektywnej metody poruszania się „w krętych wąskich dolinach”.
Jak już powiedzieliśmy algorytm gradientu sprzężonego jest realizacją pierwszego z zaproponowanych podejść wzbogacenia algorytmu gradientu prostego w pamięć – dodaniu poprawki zależnej od historii. Zastanówmy się teraz nieco głębiej nad tym jak ta zależność od historii wygląda.
Przepiszmy wzór z piątego kroku określający kolejne punkty generowane przez algorytm gradientu sprzężonego (zapominając o potencjalnej odnowie wynikającej z nie wygenerowania kierunku poprawy rozpoczynanej w trzecim kroku algorytmu), a także wyrażenia określające jego poszczególne składniki (przyjmujemy, że
Widać, że algorytm ma dwa rodzaje pamięci zawarte w składniku wzoru
- Pierwszą, potrzebną do obliczenia bieżącej wartości współczynnika skalującego
Jej stan nie zależy od wartości tego współczynnika w kroku poprzednim, a z wielkości wyznaczonych uprzednio wykorzystuje tylko gradient. Taką pamięć można nazwać krótkookresową, bo nie ma w niej efektu kumulacji informacji zawartej w stanach poprzednich współczynnika skalującego.
- Drugą – określoną przez sposób obliczania kierunków poprawy
Po pierwsze, jej stan d jest w naturalny sposób opóźniony o takt (kierunek poprawy w takcie k możemy obliczyć tylko na podstawie kierunku wyznaczonego w takcie poprzednim). Bieżący stan tej pamięci zależy od jej stanu poprzedniego, co wobec tego że wzór jest sumą oznacza, że w kolejnym wektorze kierunku poprawy jest w swoisty sposób „podsumowywana” cała dotychczasowa historia działania algorytmu. Jest to więc pamięć długookresowa. Kumulacja niedokładności obliczeniowych w tej pamięci może powodować, wspomniane wyżej „prowadzenie na manowce”. Odnowa, to oczyszczenie pamięci długookresowej ze złej informacji i wprowadzenie do niej kierunku, o którym wiemy, że na pewno jest kierunkiem poprawy.