Bezpieczne stosowanie metod iteracyjnych na matrycach diagonalnie dominujących

9

Załóżmy, że podano następujący układ liniowy

(1)Lx=c,
gdzie L jest ważonym Laplacianem znanym jako dodatni semiokreślony z jednowymiarową przestrzenią zerową rozciągniętą przez 1n=(1,,1)Rnoraz wariancja tłumaczenia xRntzn. x+a1n nie zmienia wartości funkcji (której pochodną jest (1)). Jedyne pozytywne wpisy zL znajdują się na jego przekątnej, która jest sumą wartości bezwzględnych ujemnych wpisów poza przekątną.

Znalazłem to w jednym z bardzo cytowanych prac naukowych w tej dziedzinie L jest not strictly dominujące po przekątnej metody takie jak Conjugate Gradient, Gauss-Seidl, Jacobi, nadal mogą być bezpiecznie stosowane do rozwiązywania (1). Uzasadnieniem jest to, że ze względu na niezmienność tłumaczenia można bezpiecznie naprawić jeden punkt (np. Usunąć pierwszy wiersz i kolumnęL i pierwszy wpis z c ), w ten sposób przekształcając L do strictlymacierz diagonalnie dominująca. W każdym razie oryginalny system został rozwiązany w pełnej formie(1), z LRn×n.

Czy to założenie jest prawidłowe, a jeśli tak, jakie są alternatywne uzasadnienie? Próbuję zrozumieć, jak nadal zachowuje się zbieżność metod.

Jeśli metoda Jacobi jest zbieżna (1), co można powiedzieć o promieniu spektralnym ρ macierzy iteracji D1(DL), gdzie D to macierz diagonalna z wpisami Lna przekątnej? Jestρ(D1(DL)1, a zatem różni się od ogólnych gwarancji konwergencji dla ρ(D1(DL))<1? Pytam o to, ponieważ wartości własne macierzy LaplacianaD1Lz tymi na przekątnej powinny być w zasięgu[0,2].

Z oryginalnej pracy:

......................................

Przy każdej iteracji obliczamy nowy układ (x (t +1), y (t + 1)), rozwiązując następujący układ liniowy:

(8)L·x(t+1)=L(x(t),y(t))·x(t)L·y(t+1)=L(x(t),y(t))·y(t)
Bez utraty ogólności możemy ustalić położenie jednego z czujników (wykorzystując stopień swobody translacji zlokalizowanego naprężenia) i uzyskać matrycę ściśle dominującą po przekątnej. Dlatego możemy bezpiecznie używać iteracji Jacobi do rozwiązywania (8)

.......................................

W powyższym pojęciu „iteracja” jest związana z podstawową procedurą minimalizacji i nie należy jej mylić z iteracją Jacobiego. Tak więc system jest rozwiązywany przez Jacobiego (iteracyjnie), a następnie rozwiązanie jest kupowane po prawej stronie (8), ale teraz w celu kolejnej iteracji podstawowej minimalizacji. Mam nadzieję, że to wyjaśnia sprawę.

Zauważ, że znalazłem Które iteracyjne solwery liniowe są zbieżne dla dodatnich macierzy półprzewodnikowych? , ale szukam bardziej szczegółowej odpowiedzi.

usero
źródło
Czy możesz opublikować link lub cytat do wysoko cytowanej pracy?
Geoff Oxberry
Można go odzyskać z: citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.164.1421 Ponieważ nie oczekuje się, że przeczytasz całą pracę, spójrz na str. 7 (na dole). Podejrzewam, że wybór solverów iteracyjnych jest uzasadniony, ale uważam, że potrzebne jest lepsze (a przynajmniej inne) uzasadnienie.
usero
Zastanawiam się, czy ci ludzie pochodzą z tej samej społeczności, co kombinatoryjne warunki wstępne.
shuhalo,

Odpowiedzi:

5

Iterację Jacobi można wykazać zbieżnością.

Pierwszą rzeczą, którą powinieneś się upewnić, jest to cT1n=0, który jest warunkiem istnienia rozwiązania (zakładam L=LT, w przeciwnym razie potrzebujesz c(KerLT)), ponieważ powiedziałeś V0:=KerL=span{1n}. Wykorzystamy konwencję, któraV0jest również macierzą z kolumnami będącymi jej podstawą ortonormalną. W Twoim przypadku,V0:=1n/n.

Następnie, w przypadku błędów iteracji Jacobi w oryginalnym systemie, masz

e1=(ID1L)e0=(ID1L)(Pe0+V0a)=(ID1L)Pe0+V0a,
gdzie P:=IV0V0 jest rzutem ortogonalnym na V1:=V0. Z powyższej iteracji wiemy to
Pe1=P(ID1L)Pe0,

z którego mamy macierz iteracji S w V1,
S:=P(ID1L)P.
Nie to S ma te same widma (oprócz zer) z następującą macierzą
S~:=(ID1L)PP=(ID1L)P=(ID1L)(IV0V0)=ID1LV0V0.
Chcemy promienia spektralnego S mniej niż jeden, aby udowodnić zbieżność.

Poniższy cytat jest stary i przechowywany wyłącznie w celach informacyjnych. Zobacz później nowy dowód.

W Twoim przypadku, V0V0=1n1n×n. I możesz to zweryfikować D1L+V0V0 jest ściśle dominujący po przekątnej, przyjmując założenie, że wpisy z Lsą dodatnie na przekątnej, a w przeciwnym razie ujemne. Aby pokazać wartości własne D1L+V0V0 są prawdziwe, zauważamy, że matryca jest samoprzylegająca pod produktem wewnętrznym <x,y>:=yTDx.

Gdyby V0nie jest w twojej konkretnej formie, nie znalazłem odpowiedzi na pytanie o zbieżność. Czy ktoś może to wyjaśnić?

Zauważ, że V0 jest wektorem własnym odpowiadającym wartości własnej 1 z ID1L. Na podstawie obserwacji nazywamy Twierdzenie 2.1 z wartości własnych zaktualizowanych macierzy rangi 1 z niektórymi aplikacjami Jiu Dinga i Ai-Hui Zhou.

Twierdzenie 2.1 Niechu i v być dwojgiem n-wymiarowe wektory kolumnowe takie, że u jest wektorem własnym A związane z wartością własną λ1. Następnie wartości własneA+uvT{λ1+uTv,λ2,,λn} zliczanie krotności algebraicznej.

Następnie wiemy, że widma S~ jest taki sam jak ID1L z wyjątkiem tego, że wartość własna 1 w tym drugim przesunięto o 1w zero wartości własnej w pierwszym. Odρ(ID1L)(1,1], mamy ρ(S~)(1,1).

Hui Zhang
źródło
Dzięki za odpowiedź. Coś podobnego było tym, co rozważałem: mianowicie, z ważonym Laplacianem skonstruowanym jakoD1L powyżej można wykazać, że jego wartości własne są wewnątrz [0,2), stąd z promieniem spektralnym w środku (0,2) (jedna wartość własna jest większa niż 0, a przynajmniej jeden to 0). Dlatego promień widmowy macierzy iteracjiID1L jest mniej niż 1, stąd z konwergentnym Jacobi. Być może powyższe założenie dotyczące promienia spektralnegoID1L (nie licząc 0) nie jest bezpieczny?
usero
Myślę, że widma D1L powinien być w [0,2], to jest zamknięte o godz 2. Nie wiem jak możesz dostać2wyłączony. Z mojego punktu widzenia (twierdzenie o okręgu Gershgorina) [ en.wikipedia.org/wiki/Gershgorin_circle_theorem] może podać tylko oszacowanie obejmujące2. W takim przypadku oszacuj promień widmowyID1L jest 1 z równością możliwą do osiągnięcia z wektorami w jądrze L. Myślę, że pożądana zbieżność polega na tym, że w ortogonalnej przestrzeni dopełniaczaV1jak zauważono w powyższej „odpowiedzi”.
Hui Zhang
Możesz rzucić okiem na Lemma 1.7 (v) math.ucsd.edu/~fan/research/cb/ch1.pdf Matryca D1L można by uznać za ważonego Laplaciana na pełnym wykresie, stąd z wyłączeniem 2. Wydaje mi się, że jest to wystarczający argument za potwierdzeniem zbieżności? ........... Czy twoje podejście wymaga innego wstępnego / końcowego przetwarzania iteracji poza centrowaniemc. Pytam, ponieważ się przedstawiłeśV0 A jeśli chodzi o widma ID1LV0V0: biorąc pod uwagę, że promień widmowy (sr) z ID1L jest (0,1], dodanie 1n, dałby sr<1. Czy to nie wystarczający argument?
usero
Cześć, dziękuję za wskazanie dobrej książki. Ale stwierdziłem, że nie mogę rzucić okiem. Jeśli chodzi o twój ostatni argument, pojawia się on prawie tak samo jak „odpowiedź” powyżej. Bądź ostrożny, nie dodajesz1n ale 1n1n×n, więc nie jest to prosty dodatek do sr z ID1L. Ogólnie rzecz biorąc,srsumy dwóch macierzy nie są prostą sumąsrposzczególnych matryc.
Hui Zhang,
Dobrze, że to wskazałeś. Czy twoje podejście wymaga innego wstępnego / końcowego przetwarzania iteracji poza centrowaniem c. Pytam, ponieważ się przedstawiłeśV0i pomyślałem, że mówisz o rzutowaniu na pustą przestrzeń. Jeśli tak, to czy projekcja jest pustą przestrzenią naprawdę niezbędną do konwergencji?
usero
5

Metody Kryłowa nigdy nie używają jawnie wymiarów przestrzeni, w której się iterują, dlatego można je uruchamiać w systemach pojedynczych, o ile iteraty są przechowywane w podprzestrzeni o wartości innej niż zero. Zwykle odbywa się to poprzez rzutowanie pustego miejsca przy każdej iteracji. Są dwie rzeczy, które mogą pójść nie tak, pierwsza jest znacznie częstsza niż druga.

  1. Kondycjoner wstępny jest niestabilny po zastosowaniu do pojedynczego operatora. Bezpośrednie solwery i niepełne rozkładanie na czynniki mogą mieć tę właściwość. Ze względów praktycznych wybieramy tylko różne warunki wstępne, ale istnieją bardziej zasadnicze sposoby projektowania warunków wstępnych dla pojedynczych systemów, np. Zhang (2010) .
  2. W pewnej iteracji x znajduje się w podprzestrzeni niepustej, ale Axżyje całkowicie w pustej przestrzeni. Jest to możliwe tylko w przypadku macierzy niesymetrycznych. Niezmodyfikowane GMRES rozkłada się w tym scenariuszu, ale zobacz Reichel i Ye (2005), aby uzyskać warianty wolne od awarii.

Aby rozwiązać pojedyncze systemy za pomocą PETSc, patrz KSPSetNullSpace(). Większość metod i warunków wstępnych może rozwiązać pojedyncze systemy. W praktyce mała zerowa przestrzeń dla PDE z warunkami brzegowymi Neumanna prawie nigdy nie stanowi problemu, o ile poinformujesz solver Krylova o pustej przestrzeni i wybierzesz rozsądny warunek wstępny.

Z komentarzy wynika, że ​​jesteś szczególnie zainteresowany Jacobi. (Dlaczego? Jacobi jest przydatny jako wygładzacz wielosieciowy, ale istnieją znacznie lepsze metody do wykorzystania jako solwery.) Jacobi zastosował doAx=b nie zbiega się, gdy wektor b ma komponent w pustym miejscu A, jednak część rozwiązania prostopadła do przestrzeni zerowej jest zbieżna, więc jeśli rzutujesz przestrzeń zerową z każdej iteracji, zbiega się. Alternatywnie, jeśli wybierzesz spójnyb i początkowe przypuszczenie, że iteraty (dokładnie arytmetyka) nie kumulują składników w przestrzeni zerowej.

Jed Brown
źródło
Możesz wykonać prostopadłą zmianę podstawy, aby na przekątnej było zero (znajdź dowolną macierz ortogonalną Qw którym pierwsza kolumna jest wektorem stałym). W ramach tej transformacjiA1=QTAQ, macierz A1jest nadal symetryczny dodatni półokreślony, ale pierwsza wartość po przekątnej wynosi 0, więc bezpośrednie zastosowanie Jacobiego zakończy się niepowodzeniem. OdA1jest gęsta, nie zrobiłbyś tego w praktyce, ale to pokazuje, że podstawa ma znaczenie. GdybyZ jest ortogonalną podstawą dla przestrzeni zerowej, przewidywany GMRES właśnie rozwiązuje (IZ)P1Ax=(IZ)P1b.
Jed Brown
Hmm, wygląda na to, że odpowiedziałem na komentarz, który został usunięty. Zostawię tutaj komentarz na wypadek, gdyby był użyteczny.
Jed Brown
Dzięki za odpowiedź, jest na znacznie wyższym poziomie specjalistycznym, niż się spodziewałem. Dlatego potrzebuję przewodników na temat: 1) jak rzutować puste miejsce przy każdej iteracji? 2) W moim rozumieniu stwierdziłeś, że aplikacja Jacobi do systemu, jak wspomniano przede wszystkim, może nie zbiegać się z dokładnym rozwiązaniem (tzn. Iterandy nie uzyskują lepszych oszacowań rozwiązania). Dlatego sugeruje się wybór różnych warunków wstępnych? Jeśli tak, to czy w praktyce oznacza to dynamiczną kontrolę zachowaniadiag(A)i zmienić, jeśli wystąpi problem (z powyższym przypadkiem układu liniowego)?
usero
Mój 1) z góry należy uznać za: biorąc pod uwagę iterację Jacobiego z pierwotnie opublikowanym systemem, czy konieczne jest wyświetlenie pustej przestrzeni, a jeśli tak, to w jaki sposób można ją włączyć do aktualizacjiXk+1=D1(b(AD)Xk)? Przetwarzanie iteratu po przetworzeniuXk+1i biorąc pod uwagę wersję przetworzoną dla Xk?
usero
1
W rozsądnym zakresie Jacobi powinien być stabilny. Może również użyć 1 na przekątnej, jeśli element matrycy po przekątnej ma wartość 0, rzut nadal usuwa pustą przestrzeń. Czy planujesz użyć metody Kryłowa, takiej jak CG lub GMRES? Jeśli nie, dlaczego nie? Jeśli tak, potrzebujesz po prostu ortogonalnej podstawy dla przestrzeni zerowej. Masz tylko tryb stały w swojej pustej przestrzeni, więc ortogonalny projektor w pustej przestrzeni jestN=ZZT gdzie Zto wektor kolumny. Tak więc jest projektor ortogonalny, który usuwa pustą przestrzeńIN. (Mój pierwszy komentarz miał błąd, jeśliZ jest podstawą N=IZZTto projektor.)
Jed Brown