Algorytm dla zadania wygładzania

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

Niech A = (~a1, ~a₂, . . . , ~a_n), gdzie ~aj są wektorami-kolumnami macierzy A. Wybierzmy pierwsze odbicie Householdera H1 = Im− ~u1~u₁^T/γ1 tak, aby przekształcało pierwszy wektor-kolumnę macierzy A na kierunek ~e1. Efektem pomnożenia macierzy A z lewej strony przez H1 będzie wtedy macierz

A⁽¹⁾ = (~a₁⁽¹⁾, . . . , ~a_n⁽¹⁾) = (H₁~a₁, . . . , H₁~a_n),

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

~v₂ ∈ R^m−1 do wektora ~u2 ∈ R^m przez dodanie zera jako jego pierwszej współrzędnej, ~u2 = (0, ~v₂)^T,

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 H_nH_n−1· · · H₂H₁A = 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 = H₁H₂· · · H_n, dostajemy rozkład macierzy na iloczyn ortogonalno-trójkątny

A = Q R. (5.3)

5.3. ALGORYTM DLA ZADANIA WYGŁADZANIA 41 Rzeczywiście, macierz Q ∈ R^m×m jest ortogonalna, bo

Q⁻¹ = (H₁H₂· · · H_n)⁻¹ = H_n⁻¹· · · H₂⁻¹H₁⁻¹

= H_n^T · · · H₂^TH₁^T = (H1H2· · · Hn)^T = Q^T.

Dyspunując rozkładem (5.3) zadanie wygładzania liniowego można rozwiązać następująco. Ponie-waż mnożenie przez macierz ortogonalną nie zmienia normy drugiej wektora, mamy

k~r k₂ = k~b − A~xk₂ = k~b − QR~xk₂ = kQ(Q^T~b − R~x)k₂ = k~c − R~xk₂,

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

R = R_I 0

!

,

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

k~r k²₂ = k~c_I− R_I~xk²₂ + k~c_IIk²₂.

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

~

x^∗ = R⁻¹_I ~c_I, 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ń Ho-useholdera 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~xmożemy dla dowolnego ~x realizować według wzoru

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

Uwzględnizjąc obecność zerowych elementów w ~uk, przejście od macierzy A^(k−1) do A^(k) kosztuje 4(m − k + 1)(n − k) operacji arytmetycznych i obliczenie jednego pierwiastka kwadratowego. Cały rozkład A = QR kosztuje więc

n

X

k=1

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

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

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×mnazywa się macierzą pseudoodwrotną do A ∈ R^m×n, o ile rank(A) = n. Dla nieosobliwych macierzy kwadratowych mamy oczywiście A⁺= A⁻¹, ponieważ wtedy ~x^∗= A⁻¹.

42 ROZDZIAŁ 5. ZADANIE WYGŁADZANIA LINIOWEGO

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=1 wektorów własnych tej macierzy, a odpowiadające im wartości własne są rzeczywiste i dodatnie, tzn.

h~ξi, ~ξji₂ =

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

λ₁ ­ λ₂ ­ · · · ­ λ_n> 0.

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

h~η_i, ~η_ji₂ = (λ_iλ_j)^−1/2hA~ξ_i, A~ξ_ji₂ = (λ_iλ_j)^−1/2h(A^TA)~ξ_i, ~ξ_ji₂

= (λ_iλj)^−1/2λih~ξi, ~ξji₂ =

( 0, i 6= j, 1, i = j.

Wektory ~ηi, 1 ¬ i ¬ n, można uzupełnić m − n dodatkowymi wektorami tak, aby cały układ (~ηi)^m_i=1 był 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 diagonalna 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 Householdera 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 u_j,i wektora ~u_j dla j + 1 ¬ i ¬ m na miejscach a[i, k], a współrzędna u_j,j i liczba γ_j odpowiednio na u[j] i gam[j], 1 ¬ j ¬ n.

for k := 1 to n do begin

{ obliczanie kolejnego odbicia Householdera } 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

5.3. ALGORYTM DLA ZADANIA WYGŁADZANIA 43 uu := aa + norm;

akk := −norm end else

begin

uu := aa − norm;

akk := norm 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 Householdera 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^ν (reprezentowana 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; mianowicie przez zastosowanie ortogonalizacji Grama-Schmidta do wektorów kolumn ~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 ~q₁, ~q₂, . . . , ~q_j
= span ~a₁, ~a₂, . . . , ~a_j
, 1 ¬ j ¬ n. (5.5) Zauważmy, że jeśli z wektorów ~qj stworzymy macierz Q = (~q1, ~q2, . . . , ~qn) ∈ 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 Householdera, macierz Q nie jest tu kwadratowa, ale za to kwadratowa jest macierz R.)

Wektory ortonormalne ~qj oraz współczynniki r_i,j macierzy R można wyznaczyć z układu równań

~a_j =

j

X

s=1

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

Mnożąc j-te równanie skalarnie przez ~q_i, 1 ¬ i ¬ j, otrzymujemy r_i,j = h~a_j, ~q_ii₂ = ~q_i^T~a_j, a stąd wzory rekurencyjne

44 ROZDZIAŁ 5. ZADANIE WYGŁADZANIA LINIOWEGO

for j := 1 to n do begin

for i := 1 to j − 1 do r_i,j := h~q_i, ~a_ji₂;

~

p_j := ~a_j − ^Pj−1_s=1r_s,j~q_s; rj,j := k~pjk₂;

~

qj := ~pj/rj,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 (~q_j)ⁿ_j=1 można formalnie 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, . . . , ~qm) 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 Householdera. Postępując dalej tak, jak w analizie algorytmu Householdera, otrzymujemy żądany wynik.

U. 5.6 Okazuje się, że algorytm ortogonalizacyjny Grama-Schmidta z U. 5.5 ma niedobre własności nume-ryczne, 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ą ortogona-lizację, najpierw do wektorów ~a_j, a potem do ~q_j⁽¹⁾, 1 ¬ j ¬ n. (Zob. Ćw. 5.8.)

Ćwiczenia

Ćw. 5.1 Pokazać, że dla macierzy pseudoodwrotnej A⁺do danej macierzy A ∈ R^m×n, rank(A) = n, macierz A⁺A = Injest identycznością w R^n×n, natomiast AA⁺jest rzutem prostopadłym w R^m×m na podprzestrzeń rozpiętą na wektorach-kolumnach macierzy A.

Ć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 wartości własne, a podprzestrzenie własne im odpowiadające mają ten sam wymiar.

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

Ćw. 5.5 Niech ~u będzie niezerowym wektorem w R^m oraz γ = k~uk²₂/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ła-dzania liniowego.

5.3. ALGORYTM DLA ZADANIA WYGŁADZANIA 45 Ćw. 5.7 Rozpatrzmy dwa rozkłady ortogonalno-trójkątne macierzy A ∈ R^m×n, rank(A) = n. Pierwszy to A = QR, gdzie Q ∈ R^m×m i R ∈ R^m×n (np. pochodzący z algorytmu Householdera), a drugi to A = Q₁R1, gdzie Q₁ ∈ R^m×n i R₁ ∈ R^n×n (np. pochodzący z algorytmu Grama-Schmidta z U. 5.5). Pokazać, że n pierwszych kolumn macierzy Q różnią się od n pierwszych kolumn macierzy Q₁ co najwyżej znakami; tzn.

istnieje macierz diagonalna D ∈ R^m×n taka, że d_i,i= ±1 dla 1 ¬ i ¬ n oraz Q₁ = D Q.

Ćw. 5.8 Zastosujmy ortogonalizację Grama-Schmidta do dwóch wektorów liniowo niezależnych ~a i ~b o nor-mach k~ak₂ = 1 = k~bk₂. Załóżmy dla uproszczenia, że w obliczeniach jedynie iloczyn skalarny tych wektorów s = h~a,~bi₂ 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, ~c i₂ = −ε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.

Ćw. 5.9 Zaprogramować algorytm Grama-Schmidta z U. 5.5 (z podwójną ortogonalizacją) rozwiązujący zadanie wygładzania.

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