Czy można rozwiązywać ukośne i stałe symetryczne układy liniowe w czasie kwadratowym po wstępnym obliczeniu?

21

Czy istnieje metoda do rozwiązania układów liniowych formy gdzie jest stałą macierzą SPD, a są dodatnimi macierzami diagonalnymi? $O(n^3+n^2 k)$ $k$ $(D_i + A) x_i = b_i$ $A$ $D_i$

Na przykład, jeżeli każdy jest skalarem, wystarczy obliczyć SVD . Jest to jednak podział na ogólne powodu braku przemienności. $D_i$ $A$ $D$

Aktualizacja : Jak dotąd odpowiedzi są „nie”. Czy ktoś ma jakąś ciekawą intuicję, dlaczego? Brak odpowiedzi oznacza, że nie ma niebanalnego sposobu kompresji informacji między dwoma nieprzybywającymi operatorami. Nie jest to zaskakujące, ale byłoby wspaniale zrozumieć to lepiej.

— Geoffrey Irving
źródło

SPD = półpodatnie określone?

— autor

Tak, choć problem jest zasadniczo taki sam bez SPD. Dodałem to ograniczenie tylko po to, aby systemy nigdy nie były osobliwe.

— Geoffrey Irving

19

Najbliższe pozytywne odpowiedzi na twoje pytanie, jakie mogłem znaleźć, dotyczą rzadkich przekątnych zaburzeń (patrz poniżej).

To powiedziawszy, nie znam żadnych algorytmów dla ogólnego przypadku, chociaż istnieje uogólnienie techniki, o której wspomniałeś, dla przesunięć skalarnych z macierzy SPD na wszystkie macierze kwadratowe:

Biorąc pod uwagę dowolną macierz kwadratową , istnieje rozkład Schura , gdzie jest jednostkowy, a jest górny trójkątny, a zapewnia rozkład Schura . Zatem Twój pomysł na wstępne obliczenia obejmuje wszystkie macierze kwadratowe za pomocą algorytmu: $A$ $A=U T U^H$ $U$ $T$ $A+\sigma I = U (T + \sigma I) U^H$ $A + \sigma I$

Oblicz co najwyżej . $[U,T]=\mathrm{schur}(A)$ $\mathcal{O}(n^3)$
Rozwiąż każdą przez in praca (środkowa inwersja jest po prostu powrót podstawienia). $(A+\sigma I) x = b$ $x := U (T +\sigma I)^{-1} U^H b$ $\mathcal{O}(n^2)$

Ta linia rozumowania ogranicza się do podejścia, o którym wspominałeś, gdy jest SPD, ponieważ rozkład Schura zmniejsza się do EVD dla normalnych macierzy, a EVD pokrywa się z SVD dla Hermitianowych macierzy dodatnich. $A$

Odpowiedź na aktualizację: dopóki nie mam dowodu, którego nie mam, odmawiam twierdzenia, że odpowiedź brzmi „nie”. Mogę jednak podać wgląd w to, dlaczego jest to trudne, a także inną podkladę, w której odpowiedź brzmi „tak”.

Podstawowa trudność polega na tym, że chociaż aktualizacja jest diagonalna, to nadal jest na ogół w pełnej randze, więc podstawowe narzędzie do aktualizacji odwrotności, formuła Shermana-Morrisona-Woodbury'ego , nie wydaje się pomóc. Chociaż skrzynia przesunięcia skalarnego ma również pełną rangę, jest to wyjątkowo wyjątkowy przypadek, ponieważ, jak wspomniałeś, dojeżdża do każdej macierzy.

Powiedziawszy to, jeśli każde było rzadkie, tj. Każdy miał matthcal nonzeros, wówczas formuła Shermana-Morrisona-Woodbury daje rozwiązanie z każdą parą . Na przykład z pojedynczą niezerową na tej pozycji po przekątnej, tak że : $D$ $\mathcal{O}(1)$ $\mathcal{O}(n^2)$ $\{D,b\}$ $j$ $D=\delta e_j e_j^H$

[A^{- 1} + δ e_{j} e_{j}^{H}]^{- 1} = A^{- 1} - \frac{δ A^{- 1} e_{j} e_{j}^{H} A^{- 1}}{1 + δ (e_{j}^{H} A^{- 1} e_{j})},

$[A^{-1}+\delta e_j e_j^H]^{-1} = A^{-1} - \frac{\delta A^{-1} e_j e_j^H A^{-1}}{1+\delta (e_j^H A^{-1} e_j)},$

gdzie jest tym standardowym wektorem podstawowym . $e_j$ $j$

Kolejna aktualizacja: powinienem wspomnieć, że wypróbowałem warunek wstępny który @GeoffOxberry zasugerował na kilku losowych matrycach SPD za pomocą PCG i, co może dziwić, wydaje się znacznie zmniejszać liczbę iteracji, gdy jest mały, ale nie, gdy jest to lub większy. $A^{-1}$ $1000 \times 1000$ $||D||_2/||A||_2$ $\mathcal{O}(1)$

— Jack Poulson
źródło

12

Jeśli jest dominujący po przekątnej dla każdego , to do rozwiązania każdego systemu w można zastosować najnowsze prace Koutisa, Millera i Penga (patrz strona Koutisa dotycząca pracy nad symetrycznymi macierzami diagonalnie dominującymi) w czasie (w rzeczywistości czas, w którym jest maksymalna liczba niezerowych wpisów w w stosunku do wszystkich , więc możesz również skorzystać ze sparaliżacji). Wówczas całkowity czas działania wyniósłby , co jest lepsze niż $(D_{i} + A)$ $i$ $\mathcal{O}(n^2 \log(n))$ $\mathcal{O}(m\log(n))$ $m$ $(D_{i} + A)$ $i$ $\mathcal{O}(n^2 \log(n) k)$ $\mathcal{O}(n^3 k)$ podejście do rozwiązywania każdego układu naiwnie za pomocą gęstej algebry liniowej, ale nieco gorszej niż kwadratowy czas pracy, o który prosisz.

Znaczące sparsity w dla wszystkich może być wykorzystywana przez rzadkie rozwiązują otrzymując algorytm, ale zgaduję że jeśli miały znaczący, to sparsity wspomniałby o tym. $(D_{i} + A)$ $i$ $\mathcal{O}(n^2 k)$

Możesz także użyć jako warunku wstępnego rozwiązania każdego systemu za pomocą iteracyjnych metod i zobacz, jak to działa. $A^{-1}$

Odpowiedź na aktualizację : @JackPaulson ma świetny punkt widzenia z punktu widzenia numerycznej algebry liniowej i algorytmów. Zamiast tego skupię się na argumentach złożoności obliczeniowej.

Złożoność obliczeniowa rozwiązania układów liniowych i złożoność obliczeniowa mnożenia macierzy są zasadniczo równe. (Zobacz Teoria złożoności algebraicznej ). Jeśli potrafisz znaleźć algorytm, który mógłby kompresować informacje między dwoma nieprzystosowanymi do pracy operatorami (ignorując pozytywną część półfinałową) i bezpośrednio rozwiązać zbiór proponowanych systemów w czasie kwadratowym w , to jest to prawdopodobnie możesz użyć takiego algorytmu do wyciągania wniosków na temat szybszego mnożenia macierzy. Trudno dostrzec, w jaki sposób można zastosować pozytywną strukturę półfinałową w gęstej, bezpośredniej metodzie dla systemów liniowych w celu zmniejszenia jej złożoności obliczeniowej. $n$

Podobnie jak @JackPaulson, nie chcę powiedzieć, że odpowiedź jest „nie” bez dowodu, ale biorąc pod uwagę powyższe powiązania, problem jest bardzo trudny i jest przedmiotem zainteresowania badawczego. Z asymptotycznego punktu widzenia najlepsze, co możesz zrobić bez użycia specjalnej struktury, to ulepszenie algorytmu Coppersmith i Winograd, uzyskując algorytm , gdzie . Algorytm ten byłby trudny do kodowania i prawdopodobnie byłby powolny w przypadku małych macierzy, ponieważ stały czynnik poprzedzający asymptotyczną ocenę jest prawdopodobnie ogromny w stosunku do eliminacji Gaussa. $\mathcal{O}(n^{\alpha}k)$ $\alpha \approx 2.375$

— Geoff Oxberry
źródło

3

Nie widziałem jeszcze konkretnego stwierdzenia, gdzie może być crossover, ale kilka renomowanych źródeł stwierdziło, że (pomijając kwestie implementacyjne), Coppersmith-Winograd nie jest w stanie pokonać standardowych metod dla rozmiarów matryc, które będą mogły zmieścić się w pamięci w najbliższej przyszłości (kilka dekad). Biorąc pod uwagę, że test porównawczy Linpack zajmuje więcej niż jeden dzień na obecnych najlepszych maszynach, nie wydaje się prawdopodobne, aby Coppersmith-Winograd kiedykolwiek był stosowany w praktyce. Strassen jest w rzeczywistości praktyczny w przypadku dużych problemów, choć jest nieco mniej stabilny numerycznie.

— Jed Brown,

To mnie nie zaskakuje. +1 za szczegóły implementacji.

— Geoff Oxberry,

6

Rozszerzenia Taylora pierwszego rzędu można użyć do poprawy konwergencji w stosunku do zwykłego opóźniania. Załóżmy, że mamy środek przygotowujący (lub czynniki bezpośrednie rozwiązać) dostępny dla , i chcemy ją wykorzystać do wstępnego . Możemy obliczyć $A+D$ $A$

\begin{aligned} A^{- 1} & = (A + D - D)^{- 1} (A + D) (A + D)^{- 1} \\ = [(A + D)^{- 1} (A + D - D)]^{- 1} (A + D)^{- 1} \\ = [I - (A + D)^{- 1} D]^{- 1} (A + D)^{- 1} \\ \approx [I + (A + D)^{- 1} D] (A + D)^{- 1} \end{aligned}

$\begin{align} A^{-1} &= (A+D-D)^{-1} (A+D) (A+D)^{-1} \\ &= [(A+D)^{-1} (A+D-D)]^{-1} (A+D)^{-1} \\ &= [I - (A+D)^{-1} D]^{-1} (A+D)^{-1} \\ &\approx [I + (A+D)^{-1} D] (A+D)^{-1} \end{align}$

gdzie do napisania ostatniego wiersza użyto rozszerzenia Taylora. Zastosowanie tego prekondycjonerze wymaga dwóch rozwiązuje z . $A+D$

Działa dość dobrze, gdy warunek wstępny jest przesunięty od 0 o wartość podobną lub większą niż operator, z którym próbujemy rozwiązać (np. ). Jeśli przesunięcie w warunku wstępnym jest mniejsze ( ), operator warunków wstępnych staje się nieokreślony. $D\gtrsim 0$ $D \lesssim \min \sigma(A)$

Jeśli przesunięcie w kondycjonerze jest znacznie większe niż u operatora, metoda ta ma tendencję do generowania liczby warunków o połowę mniejszej niż w przypadku kondycjonowania przez operatora opóźnionego (w losowych testach, które przeprowadziłem, może być lepsza lub gorsza dla konkretnej klasy matryce). Ten współczynnik 2 w liczbie warunków daje współczynnik w liczbie iteracji. Jeśli koszt iteracji jest zdominowany przez rozwiązania z , nie jest to wystarczający czynnik, aby uzasadnić rozszerzenie Taylora pierwszego rzędu. Jeśli aplikacja matrycowa jest proporcjonalnie droga (np. Masz niedrogi w użyciu preparat wstępny dla ), wówczas ta metoda pierwszego rzędu może mieć sens. $\sqrt 2$ $A+D$ $A+D$

— Jed Brown
źródło