Dlaczego źle uwarunkowane układy liniowe można precyzyjnie rozwiązać?

13

Zgodnie z odpowiedzią tutaj duża liczba warunków (dla liniowego rozwiązywania układu) zmniejsza gwarantowaną liczbę poprawnych cyfr w rozwiązaniu zmiennoprzecinkowym. Matryce różnicowania wyższego rzędu w metodach pseudospektralnych są zazwyczaj bardzo źle uwarunkowane. Dlaczego zatem wciąż są to bardzo dokładne metody?

Rozumiem, że niska precyzja pochodząca ze źle uwarunkowanych matryc jest jedynie wartością gwarantowaną , ale wciąż zastanawiam się, dlaczego źle uwarunkowane matryce są dokładnie rozwiązywane metodami bezpośrednimi w praktyce - np. LCOLKolumny tabeli 3.1 na stronie 11 Wang i wsp., DOBRZE STANOWANA METODA KOLOKACJI Z WYKORZYSTANIEM MATERIAŁU INTEGRACJI PSEUDOSPEKTRALNEJ , SIAM J. Sci. Comput., 36 (3) .

Zoltán Csáti
źródło
2
Moją intuicją jest to, że rozpuszczalność / dokładność układu Ax = b jest związana z wektorem wymuszającym b, a nie tylko z matrycą A. Być może, jeśli b nie „bada” ani „nie wzbudza” źle uwarunkowanych trybów A, wówczas dokładne rozwiązanie pozostaje możliwe. Jako ograniczający przykład, A może być dokładnie pojedynczą (nieskończoną liczbą warunkową), ale Ax = b może nadal posiadać rozwiązanie, które można dokładnie obliczyć, jeśli dane wymuszające b znajdują się w zakresie A. Przyznaję, że to ładna ręka -fale, dlatego komentuję tylko zamiast odpowiedzi.
rchilton1980
@ rchilton1980 „jeszcze Ax = b może nadal posiadać rozwiązanie” Ale to rozwiązanie nie jest unikalne. A przykłady, o których mówię, mają unikalne rozwiązanie.
Zoltán Csáti
To uczciwy kontrapunkt - być może artefakt wybierania nieskończonej liczby warunków (dokładnie zero wartości własnej). Myślę jednak, że można zastąpić tę zerową wartość własną maszyną epsilon i mój punkt widzenia jest nadal słuszny. (To znaczy, że system ma bardzo dużą liczbę stanów, system nie jest jednostkowy z unikalnym rozwiązaniem, które możemy bardzo dokładnie obliczyć pod warunkiem, że b nie ma elementu wzdłuż tej małej pary własnej).
rchilton1980
1
Mówiąc ściślej, mój eksperyment myślowy jest taki jak A = diag ([1 1 1 1 1 eps]), b = [b1 b2 b3 b4 b5 0]. Jest wymyślony, ale myślę, że uzasadnienie pierwotnego twierdzenia jest wystarczające: „czasami źle uwarunkowane A można dokładnie rozwiązać dla poszczególnych wyborów b”
rchilton1980,
1
Podaj kolejny przykład z blogu Moler blogs.mathworks.com/cleve/2015/02/16/…
percusse

Odpowiedzi:

7

Dodano po mojej pierwszej odpowiedzi:

Wydaje mi się teraz, że autor cytowanej pracy podaje numery warunków (pozornie 2-normowe numery warunków, ale być może numery warunków nieskończoności) w tabeli, jednocześnie podając maksymalne błędy bezwzględne zamiast normalnych błędów względnych lub maksymalnych błędów względnych elementarnych ( są to różne miary.) Zauważ, że maksymalny błąd względny elementu nie jest tym samym, co błąd względny normy nieskończoności. Ponadto błędy w tabeli odnoszą się raczej do dokładnego rozwiązania pierwotnego problemu wartości granicy równania różniczkowego, a nie dyskretnego liniowego układu równań. Tak więc informacje zawarte w artykule naprawdę nie są odpowiednie do użycia z błędem związanym z numerem warunku.

Jednak w mojej replikacji obliczeń widzę sytuacje, w których błąd względnej normy nieskończoności (lub błąd względny dwu-normalny) jest znacznie mniejszy niż wartość graniczna ustawiona przez numer warunku nieskończoności (odpowiednio numer warunku 2-normowego). Czasami po prostu masz szczęście.

Użyłem pakietu DMSUITE MATLAB i rozwiązałem przykładowy problem z tego artykułu, stosując metodę pseudospektralną z wielomianami Czebyszewa. Moje liczby warunków i maksymalne błędy bezwzględne były podobne do podanych w pracy.

Widziałem także błędy względne normy, które były nieco lepsze niż można by się spodziewać na podstawie liczby warunków. Na przykład na przykładowym problemie z , używając , otrzymujęN = 1024ϵ=0.01N=1024

cond (A, 2) = 7,9e + 8

cond (A, inf) = 7,8e + 8

norma (u-uexact, 2) / norma (uexact, 2) = 3,1e-12

norm (u-uact, inf) / norma (uexact, inf) = 2,7e-12

Wydaje się, że rozwiązanie jest dobre na około 11-12 cyfr, a numer warunku jest rzędu 1e8.

Jednak sytuacja z błędami elementarnymi jest bardziej interesująca.

max (abs (u-uact)) = 2,7e-12

To nadal wygląda dobrze.

max (abs ((u-uact) ./ uexact) = 6,1e + 9

Wow - w co najmniej jednym elemencie rozwiązania występuje bardzo duży błąd względny.

Co się stało? Dokładne rozwiązanie tego równania ma małe elementy (np. 1,9e-22), podczas gdy przybliżone rozwiązanie osiąga znacznie większą wartość 9e-14. Jest to ukryte przez pomiar błędu względnego normy (niezależnie od tego, czy jest to norma 2, czy norma nieskończoności) i staje się widoczny tylko wtedy, gdy spojrzysz na błędy względne elementarne i weźmiesz maksimum.

Moja oryginalna odpowiedź poniżej wyjaśnia, dlaczego można uzyskać błąd względny normy w rozwiązaniu, który jest mniejszy niż granica podana przez numer warunku.


Jak zauważyłeś w pytaniu, liczba warunków macierzy niepodzielnej daje najgorszy możliwy błąd względny związany z rozwiązaniem zaburzonego układu równań. Oznacza to, że jeśli rozwiążemy A ( x + Δ x ) = b + Δ b dokładnie i rozwiążemy A x = b dokładnie, toκ(A)A(x+Δx)=b+ΔbAx=b

Δxxκ(A)Δbb

Numery stanów można obliczyć w odniesieniu do różnych norm, ale często stosuje się dwuznakowy numer warunku, a jest to numer warunku użyty w dokumencie, do którego się odwołujesz.

Najgorszy przypadek występuje wtedy, gdy błąd jest lewa pojedynczej wektor odpowiada najmniejszej liczbie pojedynczej wartości . Najlepszy przypadek występuje wtedy, gdy jest lewy pojedynczej wektora odpowiada wielkości pojedynczej wartości . Gdy jest losowa, musisz spojrzeć na rzuty na wszystkie lewe wektory w liczbie pojedynczej i odpowiadające im wartości w liczbie pojedynczej. W zależności od spektrum , sprawy mogą pójść bardzo źle lub bardzo dobrze. A A hemibursztynianu b A A hemibursztynianu b hemibursztynianu b A AΔbAAΔbAAΔbΔbAA

Rozważ dwie macierze , obie o 2-normowym stanie warunku . Pierwsza macierz ma liczby pojedyncze , , , . Druga macierz ma pojedyncze wartości , , , , . 1,0 × 10 10 1 1 x 10 - 10 ... 1 x 10 - 10 1 1 ... 1 1 x 10 - 10A1.0×101011×10101×10101111×1010

W pierwszym przypadku jest mało prawdopodobne, aby przypadkowe zaburzenie było w kierunku pierwszego lewego pojedynczego wektora, a bardziej prawdopodobne, aby znajdowało się w pobliżu jednego z pojedynczych wektorów o wartości pojedynczej . Zatem względna zmiana w rozwiązaniu będzie prawdopodobnie bardzo duża. W drugim przypadku prawie każde zaburzenie będzie zbliżone do pojedynczego wektora o pojedynczej wartości , a względna zmiana w roztworze będzie niewielka. 1×10101

PS (dodane później po powrocie z zajęć jogi ...)

Wzór na rozwiązanie toAΔx=Δb

Δx=VΣ1UTΔb=i=1nUiTΔbσiVi

Według twierdzenia Pitagorasa

Δx22=i=1n(UiTΔbσi)2

Jeśli zatrzymamy , wtedy ta suma jest zmaksymalizowana, gdy i zminimalizowana, gdy .Δb2=1Δb=UnΔb=U1

W rozważanej tutaj sytuacji jest wynikiem losowych błędów zaokrąglania, więc wartości powinny być w przybliżeniu tej samej wielkości. Warunki z mniejszymi wartościami przyczynią się dużo do błędu, podczas gdy warunki z większymi wartościami nie przyczynią się wiele. W zależności od spektrum może to być znacznie mniejsze niż w najgorszym przypadku. ΔbUiTΔbσiσi

Brian Borchers
źródło
Czy ten argument nie sugerowałby, że możliwe jest (nawet jeśli jest to mało prawdopodobne) osiągnięcie najgorszego przypadku dla macierzy w przykładzie? AFAIU, na podstawie mojej odpowiedzi i dokumentacji nie powinno być to możliwe. κ(A)?getrs
Kirill
@BrianBorchers Czy mógłbyś wyjaśnić, dlaczego „Najgorszy przypadek występuje, gdy jest lewym osobliwym wektorem odpowiadającym najmniejszej liczbie pojedynczej Najlepszy przypadek występuje, gdy jest lewym pojedynczym wektorem odpowiadającym największa pojedyncza wartość ” trzyma? Z poniższego przykładu jest to logiczne, ale potrzebowałbym kilku formuł. Niech SVD jest . W pierwszym przypadku, . Jak postępować? ΔbAAΔbAAAA=UΣVTA=Δbσ1v1T+i=2NuiσiviT
Zoltán Csáti
Nie omawiałem błędów zaokrąglania w macierzy , ale ogólny efekt jest podobny - chyba że naprawdę nie dopisze ci szczęście w zaokrąglaniu błędów, zazwyczaj robisz coś lepiej niż pesymistyczne ograniczenie w najgorszym przypadku. A
Brian Borchers
(-1) Omówienie względnych błędów w danych wyjściowych jest bardzo mylące.
Kirill
1

tl; dr Zgłosili się liczbę warunek, niekoniecznie właściwą liczbę Warunkiem matrycy, ponieważ istnieje różnica.

Jest to specyficzne dla macierzy i wektora po prawej stronie. Jeśli spojrzysz na dokumentację*getrs , to zobaczysz , że granica błędu przekazywania to Tutaj nie jest całkiem zwykłym numerem warunku , ale raczej (Wewnątrz normy są to bezwzględne wartości składowe.) Zobacz na przykład iteracyjne udoskonalenie układów liniowych i LAPACK firmy Higham lub Dokładność i stabilność algorytmów numerycznych Highama (7.2).

xx0xcond(A,x)ucond(A)u.
cond(A,x)κ(A)
cond(A,x)=|A1||A||x|x,cond(A)=|A1||A|.

Dla twojego przykładu wziąłem pseudospektralny operator różnicowy dla podobnego problemu z , i faktycznie istnieje duża różnica międzyi , obliczyłem i , co wystarcza do wyjaśnienia obserwacji, że dzieje się to dla wszystkich prawych stron, ponieważ rzędy wielkości z grubsza odpowiadają temu, co jest patrz Tabela 3.1 (lepsze błędy o 3-4 zamówienia). To nie działa, gdy próbuję to samo tylko na losowej źle uwarunkowanym matrycy, więc to musi być właściwością .n=128|A1||A|κ(A)7×1032.6×107A

Wyraźnym przykładem, dla którego dwie liczby nie pasują do siebie, którą wziąłem od Highama (7.17, s. 124), z powodu Kahana, jest Innym przykładem, który znalazłem, jest zwykła macierz Vandermonde z losowym . Przeszedłem i niektóre inne źle uwarunkowane matryce również dają tego typu wynik, jak i .

(2111ϵϵ1ϵϵ),(2+2ϵϵϵ).
[1:10]bMatrixDepot.jltriwmoler

Zasadniczo chodzi o to, że analizując stabilność rozwiązywania układów liniowych w odniesieniu do zaburzeń, najpierw musisz określić, które zaburzenia rozważasz. Podczas rozwiązywania układów liniowych za pomocą LAPACK, to ograniczenie błędu uwzględnia perturbacje składowe w , ale brak perturbacji w . To różni się od zwykłego, która uwzględnia zaburzenia normalne zarówno w jak i .Abκ(A)=A1AAb

Zastanów się (jako kontrprzykład) również, co by się stało, gdybyś nie rozróżniał. Wiemy, że stosując iteracyjne udoskonalanie z podwójną precyzją (patrz link powyżej) możemy uzyskać najlepszy możliwy błąd względny w przód dla tych macierzy z . Więc jeśli weźmiemy pod uwagę pomysł, że układów liniowych nie da się rozwiązać z dokładnością lepszą niż , to jak udałoby się dopracować rozwiązania?O(u)κ(A)1/uκ(A)u

PS Ważne jest, że ?getrsobliczone rozwiązanie jest prawdziwym rozwiązaniem (A + E)x = bz perturbacją w , ale bez perturbacji w . Sytuacja wyglądałaby inaczej, gdyby perturbacje były dozwolone w .EAbb

Edytuj Aby pokazać temu działającemu bardziej bezpośrednio, w kodzie, że nie jest to przypadek ani szczęście, ale raczej (nietypowa) konsekwencja dwóch liczb warunków bardzo różniących się dla niektórych określonych macierzy, np.

cond(A,x)cond(A)κ(A).
function main2(m=128)
    A = matrixdepot("chebspec", m)^2
    A[1,:] = A[end,:] = 0
    A[1,1] = A[end,end] = 1
    best, worst = Inf, -Inf
    for k=1:2^5
        b = randn(m)
        x = A \ b
        x_exact = Float64.(big.(A) \ big.(b))
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        best, worst = min(best, err), max(worst, err)
    end
    @printf "Best relative error:       %.3e\n" best
    @printf "Worst relative error:      %.3e\n" worst
    @printf "Predicted error κ(A)*ε:    %.3e\n" cond(A, Inf)*eps()
    @printf "Predicted error cond(A)*ε: %.3e\n" norm(abs.(inv(A))*abs.(A), Inf)*eps()
end

julia> main2()
Best relative error:       2.156e-14
Worst relative error:      2.414e-12
Predicted error κ(A)*ε:    8.780e-09
Predicted error cond(A)*ε: 2.482e-12

Edycja 2 Oto kolejny przykład tego samego zjawiska, w którym różne liczby warunków nieoczekiwanie różnią się znacznie. Tym razem Tutaj jest macierzą Vandermonde 10 × 10 na , a gdy jest wybierane losowo, jest zauważalnie mniejsze niż , a najgorszy przypadek jest podany przez dla niektórych .

cond(A,x)cond(A)κ(A).
A1:10xcond(A,x)κ(A)xxi=iaa
function main4(m=10)
    A = matrixdepot("vand", m)
    lu = lufact(A)
    lu_big = lufact(big.(A))
    AA = abs.(inv(A))*abs.(A)
    for k=1:12
        # b = randn(m) # good case
        b = (1:m).^(k-1) # worst case
        x, x_exact = lu \ b, lu_big \ big.(b)
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        predicted = norm(AA*abs.(x), Inf)/norm(x, Inf)*eps()
        @printf "relative error[%2d]    = %.3e (predicted cond(A,x)*ε = %.3e)\n" k err predicted
    end
    @printf "predicted κ(A)*ε      = %.3e\n" cond(A)*eps()
    @printf "predicted cond(A)*ε   = %.3e\n" norm(AA, Inf)*eps()
end

Średnia wielkość liter (prawie 9 rzędów lepszych błędów):

julia> T.main4()
relative error[1]     = 6.690e-11 (predicted cond(A,x)*ε = 2.213e-10)
relative error[2]     = 6.202e-11 (predicted cond(A,x)*ε = 2.081e-10)
relative error[3]     = 2.975e-11 (predicted cond(A,x)*ε = 1.113e-10)
relative error[4]     = 1.245e-11 (predicted cond(A,x)*ε = 6.126e-11)
relative error[5]     = 4.820e-12 (predicted cond(A,x)*ε = 3.489e-11)
relative error[6]     = 1.537e-12 (predicted cond(A,x)*ε = 1.729e-11)
relative error[7]     = 4.885e-13 (predicted cond(A,x)*ε = 8.696e-12)
relative error[8]     = 1.565e-13 (predicted cond(A,x)*ε = 4.446e-12)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Najgorszy przypadek ( ):a=1,,12

julia> T.main4()
relative error[ 1]    = 0.000e+00 (predicted cond(A,x)*ε = 6.608e-13)
relative error[ 2]    = 1.265e-13 (predicted cond(A,x)*ε = 3.382e-12)
relative error[ 3]    = 5.647e-13 (predicted cond(A,x)*ε = 1.887e-11)
relative error[ 4]    = 8.895e-74 (predicted cond(A,x)*ε = 1.127e-10)
relative error[ 5]    = 4.199e-10 (predicted cond(A,x)*ε = 7.111e-10)
relative error[ 6]    = 7.815e-10 (predicted cond(A,x)*ε = 4.703e-09)
relative error[ 7]    = 8.358e-09 (predicted cond(A,x)*ε = 3.239e-08)
relative error[ 8]    = 1.174e-07 (predicted cond(A,x)*ε = 2.310e-07)
relative error[ 9]    = 3.083e-06 (predicted cond(A,x)*ε = 1.700e-06)
relative error[10]    = 1.287e-05 (predicted cond(A,x)*ε = 1.286e-05)
relative error[11]    = 3.760e-10 (predicted cond(A,x)*ε = 1.580e-09)
relative error[12]    = 3.903e-10 (predicted cond(A,x)*ε = 1.406e-09)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Edycja 3 Kolejnym przykładem jest macierz Forsythe, która jest zaburzonym blokiem Jordana o dowolnym rozmiarze postaci Ma to , , więc , ale , więc . I jak można zweryfikować ręcznie, rozwiązywanie układów równań liniowych, takich jak jest niezwykle dokładne, pomimo potencjalnie nieograniczonego . Więc ta matryca przyniesie nieoczekiwanie precyzyjne rozwiązania.

A=(010000100001ϵ000).
A=1A1=ϵ1κ(A)=ϵ1|A1|=A1=|A|1A x = b κ ( A )cond(A)=1Ax=bκ(A)

Edycja 4 macierzy Kahana jest również taka, z :cond(A)κ(A)

A = matrixdepot("kahan", 48)
κ, c = cond(A, Inf), norm(abs.(inv(A))*abs.(A), Inf)
@printf "κ=%.3e c=%.3e ratio=%g\n" κ c (c/κ)

κ=8.504e+08 c=4.099e+06 ratio=0.00482027
Cyryl
źródło
Numery warunków w dokumencie, do którego odnosi się OP, są numerami warunków o dwóch normach. Jeśli wrócisz do odniesienia [17] ElBarbary, zobaczysz, że we wcześniejszej pracy były to dwuznakowe numery warunków. Poza tym ustawiłem przykłady z tego artykułu za pomocą DMsuite i uzyskałem prawie dokładnie takie same 2-normowe numery warunków, jakie podano w artykule.
Brian Borchers
Numery norm warunków nieskończoności dla tych przykładów, które otrzymałem przy użyciu interpolacji dmsuite i Chebysheva, były podobne pod względem wielkości do dwuwymiarowych liczb warunków. Nie sądzę, aby różnica między 2-normą w liczbach warunków normy nieskończoności była tak ważna dla tego konkretnego przykładu.
Brian Borchers
Uważam, że błędy podane w artykule są błędami absolutnymi, a nie względnymi (nie robi to zbyt dużej różnicy, z wyjątkiem , gdzie rozwiązanie spada w dół do zeraϵ=0.01
Brian Borchers,
Dla i błędy względne dla części rozwiązania zbliżonych do 0 są ogromne, ale błędy bezwzględne są niewielkie. Zgadzam się, że artykuł był bardzo niejasny w kwestii tego, jaki numer warunku został użyty, i czym dokładnie były „błędy” (błędy względne lub bezwzględne).N = 1024ϵ=0.01N=1024
Brian Borchers,
@BrianBorchers Nie jestem pewien, co masz na myśli: to nie jest różnica między liczbami warunków 2-normalnych i infty-normalnych, ale raczej liczbami warunkowymi normalnie i składowymi (względne zaburzenia na wejściu, a nie składowe -względne błędy względne w danych wyjściowych jak w odpowiedzi).
Kirill