Rozwiązanie rzadkiego i bardzo źle uwarunkowanego systemu

Zamierzam rozwiązać Ax = b, gdzie A jest złożona, rzadka, niesymetryczna i bardzo źle uwarunkowana (liczba warunku ~ 1E + 20) macierz kwadratowa lub prostokątna. Udało mi się dokładnie rozwiązać system za pomocą ZGELSS w LAPACK. Jednak wraz ze wzrostem stopni swobody w moim systemie rozwiązanie systemu na PC za pomocą ZGELSS zajmuje dużo czasu, ponieważ rzadkość nie jest wykorzystywana. Ostatnio wypróbowałem SuperLU (używając pamięci Harwell-Boeinga) dla tego samego systemu, ale wyniki były niedokładne dla stanu warunku> 1E + 12 (nie jestem pewien, czy jest to problem numeryczny z obrotem).

Jestem bardziej skłonny do używania już opracowanych rozwiązań. Czy istnieje solidny solver, który może rozwiązać system, o którym wspomniałem szybko (tj. Wykorzystując rzadkość) i niezawodnie (biorąc pod uwagę liczby warunków)?

linear-solver sparse-matrix lapack condition-number użytkownik1234
źródło

Czy potrafisz to zrobić? Jeśli tak, metody podprzestrzeni Kryłowa mogą być skuteczne. Nawet jeśli nalegasz na bezpośrednie metody, przygotowanie wstępne pomoże kontrolować błędy numeryczne.

Geoff Oxberry

Mam również dobre doświadczenia z przygotowaniem wstępnym w sposób opisany tutaj: en.wikipedia.org/wiki/… Przygotowanie wstępne można wykonać z dokładną arytmetyką. Moje matryce są jednak gęste, więc nie mogę wskazać tutaj bardziej szczegółowych metod / procedur.

AlexE

Dlaczego liczba stanów jest tak duża? Być może formułę można poprawić, aby system był lepiej warunkowany? Ogólnie rzecz biorąc, nie możesz oczekiwać, że będziesz w stanie oszacować resztę dokładniej niż , co sprawia, że Kryłow ma małą wartość, gdy zabraknie Ci bitów. Jeśli liczba warunków jest naprawdę , powinieneś użyć poczwórnej precyzji ( z GCC, obsługiwanym przez kilka pakietów, w tym PETSc).

(machine precision) \cdot (condition number)

$(\text{machine precision})\cdot(\text{condition number})$

10^{20}

$10^{20}$ __float128

Jed Brown

Skąd bierzesz ten szacunkowy numer warunku? Jeśli poprosisz Matlaba o oszacowanie liczby warunków macierzy z zerową spacją, może to dać ci nieskończoność, a czasem może po prostu dać naprawdę ogromną liczbę (jak to, co masz). Jeśli system, na który patrzysz, ma puste miejsce i wiesz, co to jest, możesz go wyświetlić, a to, co zostało, może mieć lepszy numer stanu. Następnie możesz użyć PETSc lub Trilinos lub co masz.

Daniel Shapero

Daniel - skrócona metoda SVD stosowana przez ZGELSS określa przestrzeń zerową (wektory osobliwe powiązane z małymi wartościami osobliwymi w SVD są podstawą dla N (A)) i znajduje rozwiązanie najmniejszych kwadratów dlaponad .

min ‖ A x - b ‖

$\min \| Ax-b \|$

p e r p (N (A))

$perp(N(A))$

Brian Borchers,

Odpowiedzi:

Kiedy używasz ZGELSS do rozwiązania tego problemu, używasz okrojonego rozkładu pojedynczej wartości, aby uregulować ten niezwykle źle uwarunkowany problem. ważne jest, aby zrozumieć, że ta procedura biblioteczna nie próbuje znaleźć rozwiązania najmniejszych kwadratów dla , ale raczej próbuje znaleźć równowagę między znalezieniem rozwiązania, które minimalizujeprzed minimalizowaniem. $Ax=b$ $\| x \|$ $\| Ax-b \|$

Zauważ, że parametr RCOND przekazany do ZGELSS może być użyty do określenia, które wartości osobliwe powinny zostać uwzględnione i wykluczone z obliczeń rozwiązania. Każda liczba pojedyncza mniejsza niż RCOND * S (1) (S (1) jest największą wartością pojedynczą) zostanie zignorowana. Nie powiedziałeś nam, jak ustawiłeś parametr RCOND w ZGELSS, i nie mamy nic na temat poziomu hałasu współczynników w macierzy lub po prawej stronie , więc trudno powiedzieć, czy użyłeś odpowiednia ilość regularyzacji. $A$ $b$

Wygląda na to, że jesteś zadowolony ze standardowych rozwiązań, które otrzymujesz dzięki ZGELSS, więc wydaje się, że regularyzacja dokonana przez skróconą metodę SVD (która znajduje rozwiązanie minimalne wśród rozwiązań najmniejszych kwadratów, które minimalizują na przestrzeni roztworów obejmujących wektory osobliwe związane z wartościami osobliwymi większymi niż RCOND * S (1)) jest dla Ciebie zadowalający. $\| x \|$ $\| Ax-b \|$

Twoje pytanie można przeformułować w następujący sposób: „Jak skutecznie uzyskać uregulowane rozwiązania dotyczące najmniejszych kwadratów dla tego dużego, rzadkiego i bardzo źle uwarunkowanego liniowego problemu najmniejszych kwadratów?”

Zalecam użycie metody iteracyjnej (takiej jak CGLS lub LSQR) w celu zminimalizowania wyraźnie uregulowanego problemu najmniejszych kwadratów

$\min \| Ax-b \|^{2} + \alpha^{2} \| x \|^{2}$

gdzie parametr regulowania jest dostosowywany, aby problem tłumienia najmniejszych kwadratów był dobrze uwarunkowany i abyś był zadowolony z wynikowych rozwiązań regulowanych. $\alpha$

Brian Borchers
źródło

Przepraszam, że nie wspomniałem o tym na wstępie. Problemem do rozwiązania jest równanie akustyczne Helmholtza z wykorzystaniem MES. Układ jest źle kondycjonowany z powodu fali płaskiej zastosowanej do przybliżenia rozwiązania.

user1234 28.04.13

Skąd pochodzą współczynniki i ? Czy są to zmierzone dane? „dokładne” wartości z projektu jakiegoś obiektu (którego w praktyce nie można obrabiać z tolerancjami 15 cyfr ...)?

A

$A$

b

$b$

Brian Borchers,

Macierze A i b są tworzone przy użyciu słabego sformułowania Helmholtz PDE, patrz: asadl.org/jasa/resource/1/jasman/v119/i3/…

user1234

Jed Brown wskazał to już w komentarzach do pytania, ale tak naprawdę niewiele można zrobić w zwykłej podwójnej precyzji, jeśli liczba warunków jest duża: w większości przypadków prawdopodobnie nie uzyskasz ani jednej cyfry dokładności w swojego rozwiązania, a co gorsza, nie możesz nawet powiedzieć, ponieważ nie możesz dokładnie ocenić reszty odpowiadającej wektorowi rozwiązania. Innymi słowy: nie jest kwestią wyboru solwera liniowego - żaden liniowy solver nie może zrobić czegoś przydatnego dla takich matryc.

Tego rodzaju sytuacje zwykle zdarzają się, ponieważ wybierasz nieodpowiednią podstawę. Na przykład, macie takie źle uwarunkowane macierze, jeśli wybraliście funkcje jako podstawę metody Galerkina. (Prowadzi to do macierzy Hilberta, która jest bardzo źle uwarunkowana.) Rozwiązaniem w takich przypadkach nie jest pytanie, który solver może rozwiązać układ liniowy, ale pytanie, czy można zastosować lepsze zasady. Zachęcam do zrobienia tego samego: pomyśl o przeformułowaniu swojego problemu, aby nie skończyć z tego rodzaju matrycami. $1,x,x^2,x^3,...$

Wolfgang Bangerth
źródło

Gdy dyskretyzujemy źle postawiony problem dla PDE, np. Równanie wstecznego ciepła, na pewno otrzymamy źle uwarunkowane równanie macierzowe. Nie jest to przypadek, który możemy rozwiązać, zmieniając formułę równania lub wybierając efektywny solver macierzowy lub poprawiając precyzję liczby zmiennoprzecinkowej. W takim przypadku [tj. Problemy z odwrotnością akustyczną] wymagana jest metoda regularyzacji.

tqviet

Najprostszym / najszybszym sposobem rozwiązania źle uwarunkowanych problemów jest zwiększenie precyzji obliczeń (brutalną siłą). Innym (choć nie zawsze możliwym) sposobem jest przeformułowanie problemu.

Może być konieczne zastosowanie poczwórnej precyzji (34 cyfry dziesiętne). Mimo że 20 cyfr zostanie utraconych podczas kursu (z powodu numeru warunku), nadal otrzymasz 14 poprawnych cyfr.

Jeśli jest to interesujące, teraz w MATLAB-ie dostępne są teraz również quad-precyzyjne rzadkie solvery .

(Jestem autorem wspomnianego zestawu narzędzi).

Pavel Holoborodko
źródło