Uwarunkowanie macierzy, a błąd w fl ν

zauważyliśmy wcześniej, mamy też kU1k∞ ¬ n2ⁿ⁻¹kAk∞. Przyjmując K = K1 + (K₂ + K₃)n²2ⁿ⁻¹ ostatecznie otrzymujemy

kEk∞ ¬ K ν kAk∞,

czyli numeryczną poprawność, ponieważ K nie zależy od ν i A. 

4.4 Uwarunkowanie macierzy, a błąd w fl

Pokazaliśmy, że eliminacja Gaussa jest numerycznie poprawna w klasie macierzy A ∈ R^n×n takich, że cond(A) ¬ M,

gdzie M < ∞ jest dowolna. Okazuje się, że wielkość uwarunkowania macierzy, cond(A), ma też zasadniczy wpływ na uwarunkowanie zadania rozwiązywania układu równań, a tym samym także na błąd wytworzony w flν przy realizacji algorytmu eliminacji Gaussa. Rzeczywiście, mamy bowiem następujące twierdzenie. (Poniżej norma wektorowa jest dowolna, ale ustalona, a norma macierzowa jest przez nią indukowana, zob. U. 1.3.)

Twierdzenie 4.3 Niech E i ~e będą zaburzeniami odpowiednio macierzy A i wektora ~b takimi, że kEk ¬ K₁ν kAk i k~e k ¬ K₂ν k~bk,

Jeśli

K₁ν cond(A) < 1

to układ zaburzony (A + E)~x = (~b + ~e ) ma jednoznaczne rozwiązanie ~z^∗ spełniające k~z^∗− ~x^∗k ¬ (K₁+ K2) ν cond(A)

1 − K₁ν cond(A)k~x^∗k.

Dowód Zauważmy najpierw, że jeśli F jest macierzą taką, że kF k < 1 to macierz (I − F ) (gdzie I jest macierzą identycznościową) jest nieosobliwa oraz

k(I − F )⁻¹k ¬ 1

1 − kF k. (4.2)

Rzeczywiście, gdyby (I − F ) była osobliwa to istniałby niezerowy wektor ~x taki, że (I − F )~x = 0, co implikuje kF~xk/k~xk = 1 i w konsekwencji kF k ­ 1. Aby pokazać (4.2) zauważmy, że

1 = kIk = k(I − F )(I − F )⁻¹k ­ k(I − F )⁻¹k − kF k k(I − F )⁻¹k = (1 − kF k) k(I − F )⁻¹k, skąd bezpośrednio wynika (4.2).

Po podstawieniu F = −A⁻¹E mamy teraz

kF k ¬ kA⁻¹k kEk ¬ K₁νkAk kA⁻¹k < 1,

co wobec równości A+E = A(I+A⁻¹E)daje, że macierz (A+E) jest nieosobliwa i układ zaburzony ma jednoznaczne rozwiązanie ~z^∗. Przedstawmy to rozwiązanie w postaci ~z^∗ = ~x^∗+ (~z^∗− ~x^∗). Rozpisując

34 ROZDZIAŁ 4. ANALIZA BŁĘDÓW W ELIMINACJI GAUSSA

Podsumowując Twierdzenia 4.2 i 4.3 otrzymujemy wniosek, który jest końcowym wynikiem tego rozdziału.

Wniosek 4.2 Niech A będzie macierzą nieosobliwą. Dla dostatecznie silnej arytmetyki,

ν K(n) cond(A)  1

(gdzie K(n) jest pewnym czynnikiem niezależnym od A i ν), eliminacja Gaussa z wyborem elementu głównego w kolumnie zastosowana do rozwiązania układu A ~x = ~b jest w fl_ν wykonalna i daje wynik fl_ν(~x^∗) spełniający nierówność

kfl_ν(~x^∗) − ~x^∗k / K(n) ν cond(A) k~x^∗k. 

Uwagi i uzupełnienia

U. 4.1 Jak zauważyliśmy w dowodzie Twierdzenia 4.1 stała kumulacji K(n) zależy zasadniczo od wzrostu maksymalnego elementu w macierzach A^(k) powstających w kolejnych krokach eliminacji. Okazuje się, że uzyskane, pesymistyczne oszacowanie 2^kmax_1¬i,j¬n|a_i,j| jest praktycznie niespotykane, chociaż teoretycznie może być osiągnięte. Przykładem jest macierz

A =

U. 4.2 Uważna analiza dowodu Twierdzenia 4.1 pokazuje, że rozkład trójkątno-trójkątny jest numerycznie poprawny bez przestawień wierszy gdy element maksymalny w kolejnych macierzach A^(k) wzrasta w sposób niezależny od A. Ma to miejsce np. wtedy gdy A jest symetryczna i dodatnio określona, albo ma dominującą główną przekątną, zob. U. 4.3 oraz Ćw. 4.5 i 4.6.

4.4. UWARUNKOWANIE MACIERZY, A BŁĄD W FL_ν 35 U. 4.3 Dla macierzy symetrycznych i dodatnio określonych, A = A^T > 0, eliminacja Gaussa jest wykonalna bez przestawień wierszy, zob. U. 3.7. W klasie tych macierzy eliminacja bez przestawień wierszy jest też numerycznie poprawna, a więc w szczególności algorytm Banachiewicza-Choleskiego jest numerycznie po-prawny. Wykażemy to pokazując, że maksymalny element w kolejnych macierzach A^(k) wzrasta co najwyżej dwukrotnie. W tym celu wykorzystamy fakt, że dla dowolnej symetrycznej i dodatnio określonwj macierzy B = (b_i,j) mamy b_i,i> 0, oraz spełniona jest nierówność

b²_i,j < b_i,ib_j,j, ∀i, j (4.3)

(zob. Ćw. 4.1). Jak wykazaliśmy w U. 3.7, każda z macierzy A^(k)jest symetryczna i dodatnio określona, Stąd

|a^(k)_i,j| = |a^(k−1)_i,j − l_i,ka^(k−1)_k,j | ¬ |a^(k−1)_i,j | + |a^(k−1)_i,k |

Ćw. 4.1 Wykazać, że macierz symetryczna 2 × 2 a c c b

!

jest dodatnio określona wtedy i tylko wtedy, gdy a > 0 i ab > c². Wywnioskować stąd nierówność (4.3) dla macierzy symetrycznych i dodatnio określonych dowolnego wymiaru. Ponadto, największy co do modułu element takiej macierzy leży na głównej diagonali.

Ćw. 4.2 Wykazać, że macierz A jest dodatnio określona wtedy i tylko wtedy gdy dla każdego ~x ∈ Rⁿ wektory A~x i ~x tworzą w Rⁿ kąt ostry.

Ćw. 4.3 Pokazać, że jeśli eliminację Gaussa z wyborem elementu głównego w kolumnie zastosujemy do macierzy trójdiagonalnej, to wzrost elementu maksymalnego macierzy nie będzie zależał od n. Dokładniej,

maxi,j,k |a^(k)_i,j| ¬ 2 max

i,j |a_i,j|.

Ćw. 4.4 Pokazać, że dla macierzy Hessenberga (a_i,j = 0 dla i ­ j + 2) eliminacja Gaussa z wyborem elementu głównego w kolumnie daje

max

i,j,k |a^(k)_i,j| ¬ (k + 1) max

i,j |a_i,j|.

Ćw. 4.5 Pokazać numeryczną poprawność algorytmu przeganiania z U. 3.6.

Ćw. 4.6 Wykazać numeryczną poprawność eliminacji Gaussa bez przestawień wierszy dla macierzy z domi-nującą przekątną, tzn. gdy

36 ROZDZIAŁ 4. ANALIZA BŁĘDÓW W ELIMINACJI GAUSSA

Ćw. 4.7 Jeśli

(A + E)~z = ~b, (4.4)

gdzie kEk_p ¬ KνkAk_p, to oczywiście dla residuum ~r = ~b − A~z mamy

k~rk_p ¬ KνkAk_pk~zkp. (4.5)

Pokazać, że dla p = 1, 2, ∞ zachodzi też twierdzenie odwrotne, tzn. jeśli spełniony jest warunek (4.5) to istnieje macierz pozornych zaburzeń E taka, że kEk_p¬ KνkAk_p oraz spełniona jest równość (4.4).

Wskazówka. Rozpatrzyć E = ~r (sgn z_i)^T_1¬i¬n/k~z k₁ dla p = 1, E = ~r ~z^T/k~z k²₂ dla p = 2, oraz E =

~

r (sgn z_k)~e_k^T/k~z k∞ dla p = ∞, gdzie k jest indeksem dla którego |z_k| = k~zk∞.

Rozdział 5

Zadanie wygładzania liniowego

W tym rozdziale zajmiemy się zadaniem wygładzania liniowego, nazywanym też często liniowym zadaniem najmniejszych kwadratów. Jest ono uogólnieniem zadania rozwiązywania kwadratowych układów równań liniowych do przypadku, gdy układ jest nadokreślony.

5.1 Układ normalny

Niech A będzie daną macierzą o m wierszach i n kolumnach, A ∈ R^m×n, taką że m ­ n = rank(A),

albo równoważnie, taką że jej wektory kolumnowe są liniowo niezależne. Niech także dany będzie wektor ~b ∈ R^m. Jasne jest, że wtedy układ równań A~x = ~b nie zawsze ma rozwiązanie - mówimy, że układ jest nadokreślony.

Zadanie wygładzania liniowegopolega na znalezieniu wektora ~x^∗ ∈ Rⁿ, który najbardziej “pasuje”

do równania w tym sensie, że minimalizuje wektor residualny ~r = ~b − A~x w normie drugiej, tzn.

k~b − A~x^∗k₂ = min

~

x∈Rⁿk~b − A~xk₂.

Przykład 5.1 Przypuśćmy, że dla pewnej funkcji f : [a, b] → R obserwujemy jej wartości fi (do-kładne lub zaburzone) w punktach ti, 1 ¬ i ¬ m. Funkcję tą chcielibyśmy przybliżyć inną funkcją w należącą do pewnej n wymiarowej przestrzeni liniowej W , np. przestrzeni wielomianów stopnia mniejszego niż n. Jakość przybliżenia mierzymy wielkością

m

X

i=1

(f_i− w(t_i))². (5.1)

Wybierając pewną bazę (wj)ⁿ_j=1w W i rozwijając w w tej bazie, w =^Pn_j=1c_jw_j, sprowadzamy problem do minimalizacji wyrażenia

m

X

i=1

f_i−

n

X

j=1

c_jw_j(ti)

!2

względem cj, a więc do zadania wygładzania liniowego. Rzeczywiście, kładąc A = (ai,j) ∈ R^m×n z a_i,j = w_j(t_i), ~b = (fi)^m_i=1 i ~x = (cj)ⁿ_j=1, wielkość (5.1) jest równa k~b − A~xk²2.

37

38 ROZDZIAŁ 5. ZADANIE WYGŁADZANIA LINIOWEGO Lemat 5.1 Zadanie wygładzania liniowego ma jednoznaczne rozwiązanie ~x^∗, które spełnia układ rów-nań

A^TA ~x = A^T~b. (5.2)

Dowód Niech P ⊂ R^m będzie obrazem A jako odwzorowania liniowego z Rⁿ w R^m, P = { A~x : ~x ∈ Rⁿ}.

Ponieważ kolumny macierzy A są liniowo niezależne, tworzą one bazę w P . Stąd dim(P ) = n i odwzo-rowanie A : Rⁿ → P jest różnowartościowe. Ponadto przestrzeń R^m z normą drugą jest przestrzenią unitarną. Residuum jest więc minimalizowane dla wektora ~x^∗ ∈ Rⁿ takiego, że A~x^∗ jest rzutem pro-stopadłym wektora ~b na podprzestrzeń P . (Przypomnijmy, że skończony wymiar P zapewnia, że taki rzut istnieje i jest wyznaczony jednoznacznie.) Równoważnie można powiedzieć, że residuum ~b − A~x^∗ jest prostopadłe do P ,

∀~x ∈ Rⁿ (A~x)^T(~b − A~x^∗) = 0, albo

∀~x ∈ Rⁿ x^T(A^T~b − A^TA~x^∗) = 0.

Otrzymaliśmy więc, że wektor A^T~b − A^TA~x^∗ jest prostopadły w Rⁿ do każdego innego wektora.

Ponieważ jedynym wektorem o tej własności jest wektor zerowy, to A^TA~x^∗ = A^T~b. 

Zauważmy, że jeśli macierz A jest kwadratowa, m = n, to rozwiązaniem zadania jest ~x^∗ = A⁻¹~b i residuum jest zerowe. Zadanie wygładzania liniowego jest więc uogólnieniem rozwiązywania kwadra-towych układów równań liniowych.

Równanie (5.2) nazywa się układem normalnym. Może ono nam sugerować sposób rozwiązania zadania wygładzania liniowego. Wystarczy bowiem pomnożyć macierz A^T przez A i rozwiązać układ normalny. Zauważmy ponadto, że macierz A^TA jest symetryczna i dodatnio określona, bo (A^TA)^T = A^TAi dla ~x 6= 0 mamy ~x^T(A^TA)~x = (A~x)^T(A~x) = kA~xk2 > 0, przy czym ostatnia nierówność wynika z faktu, że kolumny macierzy A są liniowo niezależne i dlatego A~x 6= ~0. Przy mnożeniu A^T przez A wystarczy więc obliczyć tylko elementy na głównej przekątnej i pod nią, a do rozwiązania równania z macierzą A^TAmożna zastosować algorytm Banachiewicza-Choleskiego opisany w U. 3.7. Jak łatwo się przekonać, koszt takiego algorytmu wynosi n²(k +n/3), przy czym dominuje koszt obliczenia macierzy A^TA.

Ma on jednak pewne wady. Mnożenie macierzy powoduje w flν powstanie “po drodze” dodatko-wych błędów, które mogą nawet zmienić rząd macierzy. Na przykład, dla macierzy

A =

5.2. ODBICIA HOUSEHOLDERA 39