Jaka jest zasada zbieżności metod podprzestrzeni Kryłowa do rozwiązywania liniowych układów równań?

24

Jak rozumiem, istnieją dwie główne kategorie iteracyjnych metod rozwiązywania liniowych układów równań:

  1. Metody stacjonarne (Jacobi, Gauss-Seidel, SOR, Multigrid)
  2. Metody podprzestrzeni Kryłowa (Gradient sprzężony, GMRES itp.)

Rozumiem, że większość metod stacjonarnych działa poprzez iteracyjne relaksowanie (wygładzanie) trybów błędu Fouriera. Jak rozumiem, metoda gradientu sprzężonego (metoda podprzestrzeni Kryłowa) działa poprzez „przechodzenie” przez optymalny zestaw kierunków wyszukiwania z mocy macierzy zastosowanej do n tego reszty. Czy ta zasada jest wspólna dla wszystkich metod podprzestrzeni Kryłowa? Jeśli nie, to jak ogólnie scharakteryzujemy zasadę konwergencji metod podprzestrzeni Kryłowa?

Paweł
źródło
2
Twoja analiza metod stacjonarnych jest obciążona przez proste problemy modelowe, ponieważ można je analizować pod kątem trybów Fouriera. Ignoruje także ukryty kierunek naprzemienny (ADI) i wiele innych metod. Celem większości „metod stacjonarnych” jest połączenie wielu prostych „przybliżonych częściowych” solverów w jeden iteracyjny solver. Celem metod Kryłowa jest przyspieszenie (lub nawet wymuszenie) zbieżności danej stacjonarnej iteracji liniowej.
Thomas Klimpel
4
Artykuł, który, jak sądzę, został napisany, aby odpowiedzieć na twoje pytania, to Ipsen i Meyer. Pomysł krystalicznych metod, Amer. Matematyka Monthly 105 (1998) s. 889–899. To cudownie dobrze napisany i wyjaśniający artykuł, dostępny tutaj .
Andrew T. Barker,
@ AndrewT.Barker: Awesome! Dzięki Andrew! :)
Paweł

Odpowiedzi:

21

Ogólnie rzecz biorąc, wszystkie metody Kryłowa zasadniczo szukają wielomianu, który jest mały, gdy jest oceniany na widmie macierzy. W szczególności ta reszta metody Kryłowa (z zerowym początkowym domysłem) może być zapisana w formien

rn=P.n(ZA)b

gdzie to jakiś monomiczny wielomian stopnia n .P.nn

Jeśli jest diagonalizowalne, przy A = V Λ V - 1 , mamyZAZA=V.ΛV.-1

rnV.P.n(Λ)V.-1b=κ(V.)P.n(Λ)b.

W przypadku, gdy jest normalny (np. Symetryczny lub jednostkowy) wiemy, że κ ( V ) = 1. GMRES konstruuje taki wielomian poprzez iterację Arnoldiego, podczas gdy CG konstruuje wielomian przy użyciu innego produktu wewnętrznego (szczegóły w tej odpowiedzi ) . Podobnie, BiCG konstruuje swój wielomian poprzez niesymetryczny proces Lanczosa, podczas gdy iteracja Czebyszewa wykorzystuje wcześniejsze informacje o widmie (zwykle szacunki największych i najmniejszych wartości własnych dla symetrycznych określonych macierzy).ZAκ(V.)=1.

Jako fajny przykład (motywowany przez Trefethen + Bau), rozważ macierz o spektrum:

Widmo macierzy

W MATLAB zbudowałem to z:

A = rand(200,200);
[Q R] = qr(A);
A = (1/2)*Q + eye(200,200);

Jeśli weźmiemy pod uwagę GMRES, który konstruuje wielomianów, które faktycznie minimalizują resztki na wszystkich wielomianach monicznych stopnia , możemy łatwo przewidzieć historię resztkową, patrząc na kandydujący wielomiann

P.n(z)=(1-z)n

co w naszym przypadku daje

|P.n(z)|=12)n

do w widmie A .zZA

Teraz, jeśli uruchomimy GMRES na losowym RHS i porównamy resztkową historię z tym wielomianem, powinny one być dość podobne (potencjalne wartości wielomianowe są mniejsze niż resztkowe GMRES, ponieważ ):b2)>1

Historia szczątkowa

Reid.Atcheson
źródło
Czy możesz wyjaśnić, co rozumiesz przez „mały w spektrum macierzy”?
Paweł
2
Traktować jako złożoną wielomianu wielomian ma mały moduł sprężystości w obszarze płaszczyzny zespolonej, które obejmuje zakres A . Wyobraź sobie wykres konturowy nałożony na wykres rozproszenia wartości własnych. Jak mały jest mały? Zależy to od problemu, czy A jest normalne, a od prawej b . Podstawową ideą jest jednak to, że sekwencja wielomianów ( P n ) dąży do stopniowego zmniejszania się widma, tak że resztkowa wartość szacunkowa w mojej odpowiedzi wynosi 0P.nZAZAb.(P.n)0 .
Reid.Atcheson
@ Reid.Atcheson: Bardzo dobrze powiedziane. Czy mogę polecić napisanie jako κ ( V ) i wspomnienie, że jest to jeden dla normalnych matryc? V.V.-1κ(V.)
Jack Poulson
Laplacian wstępnie przygotowany przez optymalną SOR ma widmo bardzo podobne do tej przykładowej matrycy. Szczegóły tutaj: scicomp.stackexchange.com/a/852/119
Jed Brown
Ściśle mówiąc, CGNE jest niezależny od widma, ponieważ zależy tylko od pojedynczych wartości.
Jed Brown
17

W sprawie norm

Jako uzupełnienie odpowiedzi Reid.Atcheson chciałbym wyjaśnić niektóre kwestie dotyczące norm. W iteracji GMRES znajduje wielomianunthP.n2)

rn=ZAxn-b=(P.n(ZA)-1)b-b=P.n(ZA)b.

ZAZAZA-1

rnZA-1=rnT.ZA-1rn=(ZAmin)T.ZA-1ZAmin=minT.ZAmin=minZA

gdzie użyliśmy błędu

min=xn-x=xn-ZA-1b=ZA-1rn

ZAZA-1ZA2)ZAT.ZAZA

Ostrość granic konwergencji

Wreszcie, istnieje interesująca literatura na temat różnych metod Kryłowa i subtelności konwergencji GMRES, szczególnie dla nietypowych operatorów.

Jed Brown
źródło
Zrezygnowałeś
11

Metody iteracyjne w pigułce:

  1. Metody stacjonarne są w istocie iteracjami stałymi punktami : do rozwiązaniaZAx=b, wybierasz odwracalną macierz do i znajdź stały punkt

    x=x+dob-doZAx
    Jest to zbieżne z twierdzeniem Banacha o punkcie stałym if ja-doZA<1. Różne metody odpowiadają konkretnemu wyborowido (np. w przypadku iteracji Jacobi, do=re-1, gdzie re jest matrycą ukośną zawierającą elementy ukośne ZA).
  2. Metody Kryłowa Metody podprzestrzeni są w istocie metodami projekcji : wybierasz podprzestrzenieU,V.don i poszukaj x~U tak, że resztkowe b-ZAx~ jest prostopadły do V.. W przypadku metod KryłowaU oczywiście jest to przestrzeń rozpięta przez moce ZAzastosowane do początkowej pozostałości. Różne metody odpowiadają następnie konkretnym wyboromV. (na przykład, V.=U dla CG i V.=ZAU dla GMRES).

    Właściwości zbieżności tych metod (i ogólnie metod projekcji) wynikają z faktu, że ze względu na odpowiedni wybór V., x~ są optymalne U(np. minimalizują błąd w normie energetycznej dla CG lub resztkowy dla GMRES). Jeśli zwiększysz wymiarU w każdej iteracji masz gwarancję (dokładnie arytmetyki), aby znaleźć rozwiązanie po skończonej liczbie kroków.

    Jak zauważył Reid Atcheson, używając spacji Kryłowa dla U pozwala udowodnić wskaźniki konwergencji pod względem wartości własnych (a tym samym liczby warunków) ZA. Ponadto mają one kluczowe znaczenie dla uzyskania wydajnych algorytmów do obliczania projekcjix~.

    Jest to dobrze wyjaśnione w książce Youcefa Saada na temat metod iteracyjnych .

Christian Clason
źródło