Algorytm dla zadania wygładzania

Odbić Hauseholdera można użyć do rozkładu macierzy A ∈ R^m×n na iloczyn ortogonalno-trójkątny.

Niech A = (~a1, ~a2, . . . , ~an), gdzie ~aj są wektorami-kolumnami macie-rzy A. Wybierzmy pierwsze odbicie Hauseholdera H1 = Im− ~u1~u^T₁/γ1

tak, aby przekształcało pierwszy wektor-kolumnę macierzy A na kieru-nek ~e1. Efektem pomnożenia macierzy A z lewej strony przez H1 będzie wtedy macierz

A⁽¹⁾ = (~a⁽¹⁾₁ , . . . , ~a⁽¹⁾_n ) = (H1~a1, . . . , H1~an),

w której pierwsza kolumna ~a⁽¹⁾₁ ma niezerową tylko pierwszą współ-rzędną. W następnym kroku wybieramy drugie przekształcenie Hause-holdera ¯H2 = Im−1− ~v2~v₂^T/γ2 wymiaru m − 1 tak, aby przeprowadzało wektor (a⁽¹⁾_i,2)^m_i=2 na kierunek pierwszego wersora w R^m−1. Rozszerzając

~v2 ∈ R^m−1 do wektora ~u2 ∈ R^m przez dodanie zera jako pierwszej współrzędnej, ~u2 = (0, ~v2)^T, otrzymujemy przekształcenie (macierz) Hauseholdera H2 = Im− ~u2~u^T₂/γ2 w R^m postaci

H₂ = 1 ~0^T

~0 ¯H2

!

.

Pomnożenie macierzy A⁽¹⁾ z lewej strony przez H2 spowoduje teraz wyzerowanie drugiej kolumny macierzy pod elementem a⁽¹⁾_2,2, przy czym pierwszy wiersz i pierwsza kolumna pozostaną niezmienione. Postępując tak dalej n razy (albo n − 1 razy gdy m = n) otrzymujemy

HnHn−1· · · H2H1A = R,

gdzie R ∈ R^m×n jest uogólnioną macierzą trójkątną górną, tzn. ri,j = 0 dla i > j. Stąd, podstawiając Q = H1H2· · · Hⁿ, dostajemy rozkład macierzy na iloczyn ortogonalno-trójkątny

A = Q · R. (5.3)

Rzeczywiście, macierz Q ∈ R^m×m jest ortogonalna, bo Q⁻¹ = (H1H2· · · Hⁿ)⁻¹ = H_n⁻¹· · · H2⁻¹H₁⁻¹

= H_n^T · · · H2^TH₁^T = (H1H2· · · Hⁿ)^T = Q^T.

5.3. ALGORYTM DLA ZADANIA WYGŁADZANIA 57 Dyspunując rozkładem (5.3) zadanie wygładzania liniowego można rozwiązać następująco. Ponieważ mnożenie przez macierz ortogonalną nie zmienia normy drugiej wektora, mamy

k~rk2 = k~b − A~xk² = k~b − QR~xk²

= kQ(Q^T~b − R~x)k2 = k~c − R~xk²,

gdzie ~c = Q^T~b = Hn· · · H²H1~b. Rozbijając wektor ~c na ~c = (~cI, ~cII)^T, gdzie ~cI ∈ Rⁿ i ~cII ∈ R^m−n, oraz macierz R na

R = R_I 0

!

,

gdzie RI ∈ R^n×njest macierzą trójkątną górną, a 0 jest macierzą zerową wymiaru (m − n) × n, otrzymujemy

k~rk²2 = k~c^I− R^I~xk²2 + k~c^IIk²2.

Rozwiązanie ~x^∗ zadania wygładzania jest więc rozwiązaniem układu liniowego trójkątnego,

~x^∗ = R⁻¹_I ~cI, oraz k~r^∗k2 = k~b − A~x^∗k2 = k~cIIk2.

Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde z kolejnych przekształceń Hauseholdera Hk wyznaczamy przez obliczenie γk oraz współrzędnych wektora ~uk. Wektor ten ma tylko m − k + 1 współrzędnych niezerowych, a ponadto uk,i= a^(k−1)_i,k dla k + 1 ≤ i ≤ m.

Dzięki takiej reprezentacji Hk, mnożenia Hk~x możemy dla dowolnego

~x realizować według wzoru

(Hk~x)i = xi − s uk,i, gdzie s = ~u^T_k~x/γk.

Uwzględnizjąc obecność zerowych elementów w ~uk, przejście od ma-cierzy A^(k−1)do A^(k)kosztuje rzędu 4(m−k+1)(n−k) operacji arytme-tycznych i obliczenie jednego pierwiastka kwadratowego. Cały rozkład A= QR kosztuje więc rzędu (dla dużych m i n)

Xn k=1

4(m − k + 1)(n − k) ≈ 4

3n³+ 2n²(m − n) = 2n²(m − n/3)

58 ROZDZIAŁ 5. ZADANIE WYGŁADZANIA LINIOWEGO operacji arytmetycznych i n pierwiastków kwadratowych. Zauważmy, że w przypadku m = n, a więc dla kwadratowego układu równań, koszt ten wynosi (4/3)n³ i jest dwa razy większy od kosztu eliminacji Gaussa.

Uwagi i uzupełnienia

U. 5.1 Pokazaliśmy, że rozwiązaniem zadania wygładzania liniowego jest wektor ~x^∗ = A⁺~b, gdzie

A⁺ = (A^TA)⁻¹A^T.

Macierz A⁺ ∈ R^n×m nazywa się macierzą pseudoodwrotną do A ∈ R^m×n, rank(A) = n. Dla nieosobliwych macierzy kwadratowych mamy oczywiście A⁺= A⁻¹, ponieważ wtedy ~x^∗ = A⁻¹.

U. 5.2 Pokażemy, że każdą macierz A^m×n, rank(A) = n, można rozłożyć na iloczyn

A = U Σ V^T, (5.4)

gdzie U ∈ R^m×m i V ∈ R^n×n są macierzami ortogonalnymi, a Σ ∈ R^m×n jest macierzą diagonalną (tzn. σ_i,j = 0 dla i 6= j) o dodatnich wyrazach σi,i. Ponieważ macierz A^TA jest symetryczna i dodatnio określona, znane twierdzenie z algebry liniowej mówi, że istnieje w Rⁿ baza ortonormalna (~ξ_j)ⁿ_j=1wektorów własnych tej macierzy, a odpowiadające im wartości własne są rzeczywiste i dodatnie, tzn.

h~ξi, ~ξ_ji2 =

( 0 i6= j, 1 i = j, oraz (A^TA)~ξ_j = λjξ~_j, gdzie

λ1 ≥ λ2≥ · · · ≥ λn> 0.

Zauważmy, że wektory ~η_i= λ^−1/2_i (A~ξ_i) są ortonormalne w R^m, bowiem h~ηi, ~η_ji2 = (λ_iλ_j)^−1/2hA~ξi, A~ξ_ji2

= (λ_iλ_j)^−1/2h(A^TA)~ξ_i, ~ξ_ji2

= (λ_iλ_j)^−1/2λ_ih~ξi, ~ξ_ji2 =

( 0 i6= j, 1 i = j.

5.3. ALGORYTM DLA ZADANIA WYGŁADZANIA 59 Wektory ~ηi, 1 ≤ i ≤ n, można uzupełnić m − n dodatkowymi wektorami tak, aby cały układ (~ηi)^m_i=1był bazą ortonormalną w R^m. Zdefiniujmy teraz macierze ortogonalne U o kolumnach ~η_i, 1 ≤ i ≤ m, oraz V o kolumnach

~ξ_j, 1 ≤ j ≤ n. Bezpośredni rachunek pokazuje, że macierz U^TAV jest dia-gonalna z wyrazami na przekątnej ^pλj. Stąd A = UΣV^T z σj,j = ^pλj, 1 ≤ j ≤ n.

Liczby^pλ_j nazywa się wartościami szczególnymi macierzy A, a rozkład (5.4) rozkładem macierzy według wartości szczególnych.

Zauważmy jeszcze, że jeśli A = UΣV^T to A⁺ = V Σ⁺U^T oraz Σ⁺ ∈ R^n×m jest macierzą diagonalną z wyrazami na diagonali równymi 1/^pλ_j. Rozwiązanie zadania wygładzania liniowego można więc zapisać jako

~x^∗ = V Σ⁺U^T~b.

U. 5.3 Przedstawimy teraz jedną z możliwych implementacji algorytmu Hauseholdera rozkładu macierzy na iloczyn A = QR. Po wykonaniu poniż-szego programu wyrazy r_i,j macierzy R dla 1 ≤ i ≤ j zostaną zapamiętane na miejscach a[i, j], współrzędne uj,i wektora ~uj dla j + 1 ≤ i ≤ m na miej-scach a[i, k], a współrzędna uj,j i liczba γj odpowiednio na u[j] i gam[j], 1 ≤ j ≤ n.

for k := 1 to n do begin

{ obliczanie kolejnego odbicia Hauseholdera } norm2 := 0.0;

for l := k to m do begin

aa := a[l, k];

norm2 := norm2 + aa∗ aa end;

norm := sqrt(norm2);

aa := a[k, k];

if (aa ≥ 0.0) then begin

uu := aa + norm;

akk := −norm end else

begin

uu := aa − norm;

akk := norm

60 ROZDZIAŁ 5. ZADANIE WYGŁADZANIA LINIOWEGO end;

gamma := norm2 + abs(aa)∗ norm;

u[k] := uu;

gam[k] := gamma;

{ modyfikacja kolumn macierzy } a[k, k] := akk;

for j := k + 1 to n do begin

s := uu∗ a[k, j];

for l := k + 1 to m do s := s + a[l.k]∗ a[l, j];

s := s/gamma;

a[k, j] := a[k, j] − s ∗ uu;

for l := k + 1 to m do a[l, j] := a[l, j] − s ∗ a[l, k]

end;

end.

U. 5.4 Można pokazać, że przedstawiony algorytm Hauseholdera rozkładu macierzy na iloczyn ortogonalno-trójkątny jest numerycznie poprawny. To znaczy, otrzymane w fl_ν macierz trójkątna górna R^ν i ortogonalna Q^ν (re-prezentowana przez n wektorów ~u_k i liczb γ_k) spełniają (A + E) = Q^νR^ν, gdzie kEk ≤ K(m, n)νkAk.

U. 5.5 Zadanie wygładzania liniowego można również rozwiązać stosując innego rodzaju rozkład ortogonalno-trójkątny macierzy A ∈ R^m×n; miano-wicie przez zastosowanie ortogonalizacji Grama-Schmidta do wektorów ko-lumn ~a_j macierzy A. W wyniku otrzymujemy ciąg wektorów ~q_j, 1 ≤ j ≤ n, tworzący układ ortonormalny w R^m, oraz spełniający

span{~q1, ~q₂, . . . , ~q_j} = span{~a1, ~a₂, . . . , ~a_j}, 1 ≤ j ≤ n. (5.5) Zauważmy, że jeśli z wektorów ~qj stworzymy macierz Q = (~q1, ~q₂, . . . , ~q_n) ∈ R^m×n, to (5.5) implikują istnienie kwadratowej macierzy trójkątnej górnej R∈ R^n×n takiej, że

A = Q· R.

(W odróżnieniu od efektu działania algorytmu Hauseholdera, macierz Q nie jest tu kwadratowa, ale za to kwadratowa jest macierz R.)

5.3. ALGORYTM DLA ZADANIA WYGŁADZANIA 61 Wektory ortonormalne ~qj oraz współczynniki ri,j macierzy R można wy-znaczyć z układu równań

~a_j = Xj s=1

r_s,j~q_s, 1 ≤ j ≤ n.

Mnożąc j-te równanie skalarnie przez ~qi, 1 ≤ i ≤ j, otrzymujemy ri,j = h~aj, ~q_ii2= ~q^T_i~a_j, a stąd wzory rekurencyjne

for j := 1 to n do begin

for i := 1 to j − 1 do ri,j := h~qi, ~a_ji2;

~

p_j := ~a_j − ^Pj−1_s=1r_s,j~q_s; r_j,j := k~pjk2;

~q_j := ~pj/r_j,j end.

Dysponując rozkładem A = QR, rozwiązanie ~x^∗ zadania wygładzania wyznaczamy z równania R~x = Q^T~b. Rzeczywiście, układ (~qj)ⁿ_j=1 można for-malnie uzupełnić do układu ortonormalnego w R^m dodając do niego pewne wektory ~qj dla n+1 ≤ j ≤ m. Tworząc macierz ortogonalną Q = (~q1, . . . , ~q_m) i rozszerzając macierz R ∈ R^n×n do macierzy R ∈ R^m×n przez dodanie do niej (m−n) zerowych wierszy, otrzymujemy A = Q R, a więc rozkład taki jak w algorytmie Hauseholdera. Postępując dalej tak, jak w analizie algorytmu Hauseholdera, otrzymujemy żądany wynik.

U. 5.6 Okazuje się, że algorytm ortogonalizacyjny Grama-Schmidta z U.

5.5 ma niedobre własności numeryczne, gdy wektory-kolumny macierzy wyj-ściowej A są “słabo liniowo niezależne”; tzn. otrzymane w fl_ν wektory ~q_j⁽¹⁾ mogą być wtedy “słabo” ortogonalne. W takim wypadku należy stosować podwójną ortogonalizację, najpierw do wektorów ~a_j, a potem do ~q_j⁽¹⁾, 1 ≤ j≤ n. (Zob. Ćw. 5.7.)

Ćwiczenia

Ćw. 5.1 Pokazać, że dla macierzy pseudoodwrotnej A⁺ do danej macierzy A ∈ R^m×n, rank(A) = n, macierz A⁺A = In jest identycznością w R^n×n, natomiast

AA⁺ = I_n 0

0 0

!

∈ R^m×m,

62 ROZDZIAŁ 5. ZADANIE WYGŁADZANIA LINIOWEGO czyli AA⁺ jest rzutem prostopadłym na podprzestrzeń rozpiętą na pierw-szych n wersorach w R^m.

Ćw. 5.2 Niech A ∈ R^m×n. Uzasadnić, że jądro macierzy A^T, ker(A^T) = { ~y ∈ R^m: A^T~y = ~0},

jest podprzestrzenią prostopadłą do obrazu A,

im(A) = { A~x ∈ R^m: ~x ∈ Rⁿ}.

Ćw. 5.3 Pokazać, że macierze A^TA i AA^T mają takie same niezerowe war-tości własne, a podprzestrzenie własne im odpowiadające mają ten sam wy-miar.

Ćw. 5.4 Uzasadnić, że dana macierz kwadratowa Q ∈ R^m×m jest ortogo-nalna wtedy i tylko wtedy gdy jej kolumny (wiersze) tworzą bazę ortonor-malną w R^m.

Ćw. 5.5 Niech ~u będzie niezerowym wektorem w R^m oraz γ = k~uk²2/2.

Uzasadnić, że algorytm obliczania

H~x = ~x − s~u, s = (~u^T~x)/γ,

według powyższego wzoru, jest numerycznie poprawny ze względu na dany wektor ~x ∈ R^m.

Ćw. 5.6 Policzyć złożoność algorytmu ortogonalizacyjnego opisanego w U.

5.5 rozwiązania zadania wygładzania liniowego.

Ćw. 5.7 Zastosujmy ortogonalizację Grama-Schmidta do dwóch wektorów liniowo niezależnych ~a i ~b o normach k~ak2 = 1 = k~bk2. Załóżmy dla uprosz-czenia, że w obliczeniach jedynie iloczyn skalarny tych wektorów s = h~a,~bi2

liczy się z błędem, fl_ν(s) = s(1 + ε), gdzie |ε| jest dodatnie i na poziomie ν. Pokazać, że wtedy dla otrzymanej w wyniku ortogonalizacji unormowanej pary wektorów ~a i ~c mamy

h~a,~ci2 = −εs p1 − s²(1 − ε²).

Wywnioskować stąd, że gdy ~a i ~b są “prawie liniowo zależne”, to ~a i ~c są dalekie od ortogonalnych.

5.3. ALGORYTM DLA ZADANIA WYGŁADZANIA 63 Ćw. 5.8 Zaprogramować algorytm Grama-Schmidta z U. 5.5 (z podwójną ortogonalizacją) rozwiązujący zadanie wygładzania.

Ćw. 5.9 Zaprogramować algorytm obliczający macierz odwrotną do danej nieosobliwej macierzy A ∈ R^n×n wykorzystujący rozkład Hauseholdera A = QR.

æ

64 ROZDZIAŁ 5. ZADANIE WYGŁADZANIA LINIOWEGO

Rozdział 6

Interpolacja wielomianowa

Dotychczas rozpatrywaliśmy zadania, w których danymi są ciągi liczb rzeczywistych. Teraz zajmiemy się zadaniami, w których danymi są funkcje o wartościach rzeczywistych. Pierwszym z nich jest zadanie in-terpolacji wielomianowej.

6.1 Sformułowanie zadania interpolacji

Niech D ⊂ R i niech F będzie pewnym zbiorem funkcji f : D → R.

Niech x0, x1, . . . , xnbędzie ustalonym zbiorem parami różnych punktów z D, zwanych później węzłami.

Powiemy, że wielomian w interpoluje funkcję f ∈ F w węzłach x^j, gdy

w(x_j) = f(x_j), 0 ≤ j ≤ n.

Oznaczmy przez Π_nprzestrzeń liniową wielomianów stopnia co naj-wyżej n o współczynnikach rzeczywistych,

Πn = { w(x) = aⁿxⁿ+ an−1xⁿ⁻¹+ · · · + a¹x+ a0 : aj ∈ R, 0 ≤ j ≤ n }.

Lemat 6.1 Dla dowolnej funkcji f : D → R istnieje dokładnie jeden wielomian w_f ∈ Πn interpolujący f w węzłach x_j, 0≤ j ≤ n.

Dowód Wybierzmy w Πn dowolną bazę wielomianów ϕj, 0 ≤ j ≤ n, Πn = span{ ϕ⁰, ϕ1, . . . , ϕn}.

65

66 ROZDZIAŁ 6. INTERPOLACJA WIELOMIANOWA Wtedy każdy wielomian z Πn można jednoznacznie przedstawić w po-staci rozwinięcia względem wybranej bazy. Warunkiem koniecznym i dostatecznym na to, aby wielomian wf(·) =^Pnj=0c_jϕ_j(·) interpolował f jest spełnienie układu n + 1 równań liniowych

Xn j=0

cjϕj(xi) = f(xi), 0 ≤ i ≤ n,

z n + 1 niewiadomymi cj, który w postaci macierzowej wygląda nastę-pująco:

Aby wykazać, że układ ten ma jednoznaczne rozwiązanie wystarczy, aby wektor zerowy był jedynym rozwiązaniem układu jednorodnego. Rze-czywiście, układ jednorodny odpowiada interpolacji danych zerowych, f(xi) = 0, ∀i. Istnienie niezerowego rozwiązania byłoby więc równo-ważne istnieniu niezerowego wielomianu stopnia nie większego od n, który miałby n + 1 różnych zer xi, co jest niemożliwe. 2

Zadanie znalezienia dla danej funkcji f jej wielomianu interpolacyj-nego stopnia co najwyżej n jest więc dobrze zdefiniowane, tzn. rozwią-zanie istnieje i jest wyznaczone jednoznacznie. Zauważmy, że wielomian interpolacyjny wf jako taki nie może być wynikiem obliczeń w naszym modelu obliczeniowym, możemy natomiast wyznaczyć jego współczyn-niki c_j w wybranej bazie.

Definicja 6.1 Niech (ϕj)ⁿ_j=0 będzie bazą w przestrzeni Πnwielomianów stopnia co najwyżej n. Zadanie (obliczeniowe) interpolacji wielomino-wej polega na obliczeniu dla danej funkcji f współczynników cj takich, że wielomian

wf(·) = ^Xn

j=0

cjϕj(·) (6.2)

interpoluje f w punktach xj, 0≤ j ≤ n.

W dokumencie Trzynaście Wykładów z Matematyki Obliczeniowej (Stron 61-72)