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

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) .

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 = $\epsilon=0.01$$N=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$\kappa(A)$$A(x+\Delta x)=b+\Delta b$$Ax=b$

$\frac{\| \Delta x \|}{\| x \|} \leq \kappa(A) \frac{\| \Delta b \|}{\| b \|}$

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 $\Delta b$$A$$A$$\Delta b$$A$$A$$\Delta b$$\Delta b$$A$$A$

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 - $A$$1.0 \times 10^{10}$$1$$1 \times 10^{-10}$$\ldots$$1 \times 10^{-10}$$1$$1$$\ldots$$1$$1 \times 10^{-10}$

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 \times 10^{-10}$$1$

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

Wzór na rozwiązanie to$A\Delta x = \Delta b$

$\Delta x = V \Sigma^{-1} U^{T} \Delta b=\sum_{i=1}^{n} \frac{U_{i}^{T}\Delta b}{\sigma_{i}} V_{i}$

Według twierdzenia Pitagorasa

$\| \Delta x \|_{2}^{2}= \sum_{i=1}^{n} \left( \frac{U_{i}^{T}\Delta b}{\sigma_{i}} \right)^{2}$

Jeśli zatrzymamy , wtedy ta suma jest zmaksymalizowana, gdy i zminimalizowana, gdy .$\| \Delta b \|_{2}=1$$\Delta b=U_{n}$$\Delta b=U_{1}$

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. $\Delta b$$U_{i}^{T}\Delta b$$\sigma_{i}$$\sigma_{i}$

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).

$$\frac{\|x-x_0\|_\infty}{\|x\|_\infty} \lesssim \mathrm{cond}(A,x)u \leq \mathrm{cond}(A)u.$$ $\mathrm{cond}(A,x)$$\kappa_\infty(A)$ $$\mathrm{cond}(A,x) = \frac{\||A^{-1}||A||x|\|_\infty}{\|x\|_\infty},\\ \mathrm{cond}(A) = \||A^{-1}||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$$\||A^{-1}||A|\|$$\kappa_\infty(A)$$7\times 10^3$$2.6\times 10^7$$A$

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 .

$$\begin{pmatrix}2&-1&1\\-1&\epsilon&\epsilon\\1&\epsilon&\epsilon\end{pmatrix}, \qquad \begin{pmatrix}2+2\epsilon\\-\epsilon\\\epsilon\end{pmatrix}.$$ [1:10]$b$MatrixDepot.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 .$A$$b$$\kappa(A) = \|A^{-1}\|\|A\|$$A$$b$

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)$$\kappa(A)\ll 1/u$$\kappa(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 .$E$$A$$b$$b$

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.

$$\mathrm{cond}(A,x) \approx \mathrm{cond}(A) \ll \kappa(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 .

$$\mathrm{cond}(A, x) \ll \mathrm{cond}(A) \approx \kappa(A).$$ $A$$1:10$$x$$\mathrm{cond}(A,x)$$\kappa(A)$$x$$x_i = i^a$$a$

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,\ldots,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 = \begin{pmatrix} 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \\ \epsilon&0&0&0 \end{pmatrix}.$$ $\|A\|=1$$\|A^{-1}\|=\epsilon^{-1}$$\kappa_\infty(A) = \epsilon^{-1}$$|A^{-1}| = A^{-1} = |A|^{-1}$A x = b κ ∞ ( A $\mathrm{cond}(A) = 1$$Ax=b$$\kappa_\infty(A)$

Edycja 4 macierzy Kahana jest również taka, z :$\mathrm{cond}(A)\ll \kappa(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