Bezpieczne stosowanie metod iteracyjnych na matrycach diagonalnie dominujących

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

$$Lx=c,\tag1$$ gdzie $L$ jest ważonym Laplacianem znanym jako dodatni $semi-$określony z jednowymiarową przestrzenią zerową rozciągniętą przez $1_n=(1,\dots,1)\in\mathbb{R}^n$oraz wariancja tłumaczenia $x\in\mathbb{R}^{n}$tzn. $x+a1_n$ nie zmienia wartości funkcji (której pochodną jest $(1)$). Jedyne pozytywne wpisy z$L$ 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 $strictly$macierz diagonalnie dominująca. W każdym razie oryginalny system został rozwiązany w pełnej formie$(1)$, z $L\in\mathbb{R}^{n\times 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 $\rho$ macierzy iteracji $D^{-1}(D-L)$, gdzie $D$ to macierz diagonalna z wpisami $L$na przekątnej? Jest$\rho(D^{-1}(D-L)\leq1$, a zatem różni się od ogólnych gwarancji konwergencji dla $\rho(D^{-1}(D-L))<1$? Pytam o to, ponieważ wartości własne macierzy Laplaciana$D^{-1}L$z 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:

$$L · x(t + 1) = L(x(t),y(t)) · x(t) \\ L · y(t + 1) = L(x(t),y(t)) · y(t) \tag 8$$ 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.

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

Pierwszą rzeczą, którą powinieneś się upewnić, jest to $c^T \mathbf{1}_n = 0$, który jest warunkiem istnienia rozwiązania (zakładam $L=L^T$, w przeciwnym razie potrzebujesz $c\in (\mathrm{Ker} L^T)^\perp$), ponieważ powiedziałeś $V_0:=\mathrm{Ker} L = \mathrm{span}\{\mathbf{1}_n\}$. Wykorzystamy konwencję, która$V_0$jest również macierzą z kolumnami będącymi jej podstawą ortonormalną. W Twoim przypadku,$V_0:=\mathbf{1}_n/\sqrt{n}$.

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

$$e_1 = (I-D^{-1}L)e_0 = (I-D^{-1}L) (P e_0 + V_0a)=(I-D^{-1}L) P e_0 + V_0a,$$ gdzie $P:=I-V_0V_0'$ jest rzutem ortogonalnym na $V_1:=V_0^\perp$. Z powyższej iteracji wiemy to $$P e_1 = P (I-D^{-1}L) P e_0,$$
z którego mamy macierz iteracji $S$ w $V_1$, $$S: = P (I-D^{-1}L) P.$$ Nie to $S$ ma te same widma (oprócz zer) z następującą macierzą $$\tilde{S}:= (I-D^{-1}L) P P=(I-D^{-1}L) P=(I-D^{-1}L)(I-V_0V_0')\\ =I-D^{-1}L-V_0V_0'.$$ 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, $V_0V_0'=\frac{1}{n}\mathbf{1}_{n\times n}.$ I możesz to zweryfikować $D^{-1}L+V_0V_0'$ jest ściśle dominujący po przekątnej, przyjmując założenie, że wpisy z $L$są dodatnie na przekątnej, a w przeciwnym razie ujemne. Aby pokazać wartości własne $D^{-1}L+V_0V_0'$ są prawdziwe, zauważamy, że matryca jest samoprzylegająca pod produktem wewnętrznym $<x,y>:=y^TDx.$

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

Zauważ, że $V_0$ jest wektorem własnym odpowiadającym wartości własnej $1$ z $I-D^{-1}L$. 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 Niech$u$ i $v$ być dwojgiem $n$-wymiarowe wektory kolumnowe takie, że $u$ jest wektorem własnym $A$ związane z wartością własną $\lambda_1$. Następnie wartości własne$A+uv^T$ są $\{\lambda_1+u^Tv,\lambda_2,\ldots,\lambda_n\}$ zliczanie krotności algebraicznej.

Następnie wiemy, że widma $\tilde{S}$ jest taki sam jak $I-D^{-1}L$ z wyjątkiem tego, że wartość własna $1$ w tym drugim przesunięto o $-1$w zero wartości własnej w pierwszym. Od$\rho(I-D^{-1}L)\subset (-1,1]$, mamy $\rho(\tilde{S})\subset (-1,1)$.

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 $A x$ż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ł do$A x = 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ójny$b$ i początkowe przypuszczenie, że iteraty (dokładnie arytmetyka) nie kumulują składników w przestrzeni zerowej.