Kryteria zatrzymania iteracyjnych solverów liniowych stosowanych w układach prawie pojedynczych

16

Ax=bAλ0Arn:=bAxnn v λ 0 A v = λ 0 v r 0r n/r 0< t o l x n - x x n - x n - 1rn/r0<tolnvλ0Av=λ0v. Załóżmy, że początkowa wartość rezydualna jest duża, to może się zdarzyć, że zatrzymamy się na ale błąd jest nadal duży. Jaki jest lepszy wskaźnik błędu w tym przypadku? Czydobry kandydat?r0rn/r0<tolxnxxnxn1

Hui Zhang
źródło
3
Możesz pomyśleć o swojej definicji „prawie pojedynczej”. Macierz (z i to matryca tożsamości) ma bardzo małą wartość własną, ale jest tak daleka od liczby pojedynczej, jak to tylko możliwe. ϵ 1 IIϵϵ1I
David Ketcheson
1
Równieżwygląda na niewłaściwy zapis. jest bardziej typowy, nie? ||rn/r0||||rn||/||r0||
Bill Barth
Tak, masz rację, Bill! Naprawię ten błąd.
Hui Zhang
1
Co z? a jaki dokładnie jest twój algorytm? bAx/b
shuhalo
2
Dodatek: Myślę, że poniższy artykuł dotyczy w dużej mierze źle kondycjonowanych systemów, o które się martwisz, przynajmniej jeśli używasz CG: Axelson, Kaporin: Szacowanie normy błędu i kryteria zatrzymania w wstępnie przygotowanych iteracjach gradientu sprzężonego. DOI: 10.1002 / nla.244
shuhalo

Odpowiedzi:

13

Proszę nie używać różnicę pomiędzy kolejnymi iteruje aby zdefiniować kryteria zatrzymania. To błędnie diagnozuje stagnację dla konwergencji. Większość niesymetrycznych iteracji macierzy nie jest monotoniczna, a nawet dokładna arytmetyka GMRES bez restartów może zastygać w przypadku dowolnej liczby iteracji (aż do wymiaru macierzy) przed nagłym zbiegiem. Zobacz przykłady w Nachtigal, Reddy i Trefethen (1993) .

Lepszy sposób na zdefiniowanie zbieżności

Zazwyczaj interesuje nas dokładność naszego rozwiązania, a nie wielkość reszty. W szczególności chcielibyśmy zagwarantować, że różnica między przybliżonym rozwiązaniem a dokładnym rozwiązaniem x jest wystarczająca | x n - x | < c dla niektórych określonych przez użytkownika c . Okazuje się, że można to osiągnąć przez znalezienie x n takiego, że | A x n - b | < c ϵ gdzie ϵ jest najmniejszą liczbą pojedynczą A ze względu naxnx

|xnx|<c
cxn
|Axnb|<cϵ
ϵA

|xnx|=|A1A(xnx)|1ϵ|AxnAx|=1ϵ|Axnb|<1ϵcϵ=c

gdzie zastosowaliśmy, że jest największą liczbą pojedynczą A - 1 (druga linia) i że x dokładnie rozwiązuje A x = b (trzecia linia).1/ϵA1xAx=b

Oszacowanie najmniejszej liczby pojedynczej ϵ

Dokładne oszacowanie najmniejszej pojedynczej wartości zwykle nie jest bezpośrednio dostępne na podstawie problemu, ale można je oszacować jako produkt uboczny gradientu sprzężonego lub iteracji GMRES. Należy zauważyć, że chociaż szacunki największych wartości własnych i wartości osobliwych jest zazwyczaj dość dobre po kilku iteracjach, dokładne oszacowanie najmniejszego EIGEN / wartość pojedynczej jest zwykle tylko uzyskane po osiągnięciu zbieżności. Przed konwergencją szacunek będzie zasadniczo znacznie większy niż wartość rzeczywista. To sugeruje, że musisz rozwiązać równania, zanim będziesz mógł zdefiniować prawidłową tolerancję c ϵ . Automatyczna tolerancja zbieżności, która wymaga dokładności podanej przez użytkownika cϵcϵcdla rozwiązania i oszacowań najmniejsza liczba pojedyncza z obecnym stanem metody Kryłowa może zbiegać się zbyt wcześnie, ponieważ oszacowanie ϵ było znacznie większe niż wartość rzeczywista.ϵϵ

Notatki

  1. Powyższa dyskusja działa również z zastąpionym przez wstępnie kondycjonowany operator P - 1 A i wstępnie kondycjonowany pozostały P - 1 ( A x n - b ) lub z odpowiednio kondycjonowanym operatorem A P - 1 i błędem P ( x n - x ) . Jeśli P - 1AP1AP1(Axnb)AP1P(xnx)P1jest dobrym kondycjonerem, kondycjonowany operator będzie dobrze kondycjonowany. W przypadku wstępnego kondycjonowania z lewej strony oznacza to, że wstępnie kondycjonowana reszta może być mała, ale prawdziwa reszta może nie być. Na prawej kondycjonowania jest łatwo zrobić mały, ale prawdziwy błąd | x n - x | Nie może być. To wyjaśnia, dlaczego wstępne kondycjonowanie z lewej strony jest lepsze do zmniejszania błędów, podczas gdy wstępne kondycjonowanie z prawej strony jest lepsze do zmniejszania pozostałych (i do debugowania niestabilnych warunków wstępnych).|P(xnx)||xnx|
  2. Zobacz tę odpowiedź, aby uzyskać więcej informacji na temat norm zminimalizowanych przez GMRES i CG.
  3. Szacunki ekstremalnych wartości pojedynczych można monitorować za -ksp_monitor_singular_valuepomocą dowolnego programu PETSc. Zobacz KSPComputeExtremeSingularValues ​​(), aby obliczyć pojedyncze wartości z kodu.
  4. Podczas korzystania z GMRES do szacowania pojedynczych wartości, bardzo ważne jest, aby nie uruchamiać ponownie (np. -ksp_gmres_restart 1000W PETSc).
Jed Brown
źródło
1
„” działa również z A zastąpionym przez kondycjonowany operator ”” - Jednak wówczas stosuje się tylko do wstępnie kondycjonowanej pozostałości jeśli zastosowano P - 1 A , odpowiednio. do warunku wstępnego błędu P - 1 δ x, jeśli zastosowano A P - 1 . P1rP1AP1δxAP1
Arnold Neumaier,
1
Dobrze, zredagowałem swoją odpowiedź. Zauważ, że prawidłowa skrzynka kondycjonowana daje ci kontrolę nad , odwijanie kondycjonera wstępnego (zastosowanie P - 1 ) zwykle wzmacnia błędy w trybach niskoenergetycznych. PδxP1
Jed Brown
6

Innym sposobem spojrzenia na ten problem jest rozważenie narzędzi z dyskretnych problemów odwrotnych, to znaczy problemów obejmujących rozwiązanie lub min | | A x - b | | 2 gdzie A jest bardzo źle uwarunkowane (tzn. Stosunek pierwszej i ostatniej liczby pojedynczej σ 1 / σ n jest duży).Ax=bmin||Axb||2Aσ1/σn

Tutaj mamy kilka metod wyboru kryterium zatrzymania, a dla metody iteracyjnej zaleciłbym kryterium krzywej L, ponieważ dotyczy ono tylko ilości, które są już dostępne (OŚWIADCZENIE: Mój doradca wprowadził tę metodę, więc jestem zdecydowanie stronniczy to). Z powodzeniem wykorzystałem to w metodzie iteracyjnej.

Chodzi o monitorowanie normy resztkowej i norma rozwiązania η k = | | x k | | 2 , gdzie x k jest iteracją k . Podczas iteracji zaczyna to rysować kształt litery L na wykresie dziennika (rho, eta), a punkt na rogu tej litery L jest optymalnym wyborem.ρk=||Axkb||2ηk=||xk||2xkk

(ρk,ηk)

abs(log(ηk)log(ηk1)log(ρk)log(ρk1))

Istnieją również bardziej szczegółowe metody znajdowania narożnika, które działają lepiej, ale wymagają przechowywania znacznej liczby iteracji. Baw się z tym trochę. Jeśli korzystasz z Matlaba, możesz skorzystać z przybornika Narzędzia do regularyzacji, które implementują niektóre z tych funkcji (w szczególności obowiązuje funkcja „narożnika”).

Należy pamiętać, że takie podejście jest szczególnie odpowiednie w przypadku problemów na dużą skalę, ponieważ dodatkowy czas obliczeniowy jest niewielki.

OscarB
źródło
1
Wielkie dzięki! Czyli na wykresie loglog (rho, eta) zaczynamy od prawej strony krzywej L i kończymy na górze L, prawda? Po prostu nie znam zasady tego kryterium. Czy możesz wyjaśnić, dlaczego zawsze zachowuje się jak krzywa L i dlaczego wybieramy róg?
Hui Zhang
||Axb||2=||e||2ebexact=b+e. Aby uzyskać więcej analiz, zobacz Hansen, PC i O'Leary, DP (1993). Zastosowanie krzywej L w regularyzacji dyskretnych źle postawionych problemów. SIAM Journal on Scientific Computing, 14. Zwróć uwagę, że właśnie zaktualizowałem post.
OscarB
4
@HuiZhang: nie zawsze jest to L. Jeśli regularyzacja jest niejednoznaczna, może to być podwójna litera L, co prowadzi do dwóch kandydatów do rozwiązania, jednego z rażącą cechą lepiej rozwiązać, a drugiego z pewnymi szczegółami lepiej. (I oczywiście mogą pojawić się inne złożone kształty.)
Arnold Neumaier
Czy krzywa L ma zastosowanie do problemów źle uwarunkowanych, w których powinno istnieć unikalne rozwiązanie? To znaczy, interesują mnie problemy Ax = b, gdzie b jest znane „dokładnie”, a A jest prawie pojedyncza, ale wciąż technicznie odwracalna. Wydaje mi się, że jeśli użyjesz czegoś takiego jak GMRES, norma twojego obecnego przypuszczenia x nie zmienia się zbytnio w czasie, szczególnie po pierwszych, jakkolwiek wielu iteracjach. Wydaje mi się, że pionowa część krzywej L występuje, ponieważ nie ma unikalnego / ważnego rozwiązania źle postawionego problemu; czy ta pionowa cecha byłaby obecna we wszystkich źle uwarunkowanych problemach?
nukeguy
W pewnym momencie osiągniesz taką pionową linię, zwykle dlatego, że błędy numeryczne w twojej metodzie rozwiązania powodują, że || Ax-b || nie maleje. Jednak masz rację, że w takich bezszumowych problemach krzywa nie zawsze wygląda jak L, co oznacza, że ​​zwykle masz kilka narożników do wyboru, a wybranie jednego z drugiego może być trudne. Uważam, że w artykule, do którego nawiązałem w moim komentarzu powyżej, krótko omówiono scenariusze pozbawione hałasu.
OscarB