Zagadnienia

4.1 Obroty Givensa
- 4.1.1 Definicja obrotu
- 4.1.2 Algorytm Gentlemana
4.2 Metoda Jacobiego
4.3 QR dla macierzy trójdiagonalnej
4.4 Rozkład według wartości szczególnych (SVD)

4. Zagadnienie własne IV

W tym rozdziale zajmiemy się jeszcze innymi nich dotychczas poznanymi metodami znajdowania wartości i wektorów własnych. Ponieważ duże znaczenie mają w nich obroty płaszczyzn, zaczniemy właśnie od definicji tych przekształceń, ich własności i pierwszych zastosowań.

4.1. Obroty Givensa

4.1.1. Definicja obrotu

Macierz obrotu o kąt φ w płaszczyźnie rozpiętej na wersorach e→i,e→j, gdze i<j, jest postaci

Oi,j=1⋱cos⁡φsin⁡φ⋱-sin⁡φcos⁡φ⋱1.

Macierz ta różni się od jednostkowej jedynie wyrazami oi,i=c=oj,j oraz oi,j=s=-oj,i. łatwo sprawdzić, że Oi,j jest przekształceniem ortogonalnym, Oi,j⁢⁢Oi,jT=I.

W obliczeniach numerycznych macierze obrotu służą najczęściej do wyzerowania określonej współrzędnej danego wektora. Rzeczywiście, obracając wektor a→=a1,…,anT w płaszczyźnie rozpiętej na e→i i e→j zmieniamy jedynie współrzędne i-tą i j-tą. Dokładniej, Oi,j⁢⁢a→=b→, gdzie

bk=⁡ai⁢cos⁡φ+aj⁢sin⁡φ,k=i,-ai⁢sin⁡φ+aj⁢cos⁡φ,k=j,ak,k≠i,j.⁢

Stąd, przyjmując

cos⁡φ=aiai2+aj2,⁢sin⁡φ=ajai2+aj2,

albo liczby do nich przeciwne, dostajemy bi=a→22 albo bi=-a→22, oraz bj=0. Z powyższych wzorów wynika, że kąt obrotu φ można zawsze tak dobrać, aby φ≤π/2.

Przekształcenia tak określone nazywamy obrotami Givensa.

Zauważmy, że stosując obroty Givensa kolejno w płaszczyznach rozpiętych na wektorach e→1 i e→j dla j=2,3,…,n-1 możemy przekształcić dany wektor a→ na kierunek pierwszego wersora, tzn.

O1,n⁢⁢⋯⁢⁢O1,3⁢⁢O1,2⁢⁢a→=±a→22⋅e→1.

Takie operacje wymagają wykonania 6⁢n mnożeń/dodawań i n-1 pierwiastkowań. Przypomnijmy, że ten sam efekt można uzyskać kosztem około 2⁢n mnożeń/dodawań i jednego pierwiastka stosując odbicia Householdera. Wydaje się więc, że obroty Givensa opłaca się stosować jedynie wtedy, gdy wektor a→ ma niewiele (proporcjonalnie do n) współrzędnych niezerowych. Nie jest to jednak do końca prawdą. Istnieje bowiem pewna modyfikacja podanych wzorów autorstwa Gentlemana pozwalająca istotnie zmniejszyć koszt, którą teraz zaprezentujemy.

4.1.2. Algorytm Gentlemana

Przypuśćmy, że chcemy zastosować serię obrotów Givensa Oik,jk dla k=1,2,…,m, gdzie k-ty obrót jest wybrany tak, że jk-ta współrzędna wektora

a→k:=Oik,jk⁢⁢⋯⁢⁢Oi1,j1⁢⁢a→

jest zerowa. Pomysł modyfikacji polega na przedstawieniu każdego z obrotów w postaci

Oik,jk⁢=⁢Dk⁢⁢Sik,jk⁢⁢Dk-1-1,

gdzie D0=I, a Dk są pewnymi nieosobliwymi macierzami diagonalnymi wybranymi tak, by mnożenie przez Sik,jk było prostsze niż mnożenie przez Oik,jk. Wtedy mamy

Oim,jm⁢⁢⋯⁢⁢Oi1,j1
	=	Dm⁢⁢Sim,jm⁢⁢Dm-1-1⁢⁢⋯⁢⁢D2⁢⁢Si2,j2⁢⁢D1-1⁢⁢D1⁢⁢Si1,j1⁢⁢D0-1
	=	Dm⁢⁢Sim,jm⁢⁢⋯⁢⁢Si2,j2⁢⁢Si1,j1.

Macierze diagonalne Dk, a tym samym i przekształcenia Sik,jk, zdefiniujemy indukcyjnie dla k=1,2,…,m. Załóżmy, że zdefiniowaliśmy już Dk-1=diag⁢d1,…,dn. Niech Oik,jk będzie obrotem o kąt φ oraz c=cos⁡φ, s=sin⁡φ.

Dla uproszczenia zapisu podstawimy D:=Dk-1, D⌃:=Dk=diag⁢d⌃,⁢…,d⌃n oraz p:=ik, q:=jk. Dalej mamy Op,q=Oik,jk=oi,ji,j=1n, Sp,q=Sik,jk=si,ji,j=1n. Ponieważ Sp,q=D⌃-1⁢⁢Op,q⁢⁢D to

si,j⁢=⁢djd⌃i⁢⁢oi,j⁢=⁢⁡c⁢dp/d⌃p,i=j=p,c⁢dq/d⌃q,i=j=q,s⁢dq/d⌃p,i=p,j=q,-s⁢dp/d⌃q,i=q,j=p,dj/d⌃i,i=j≠p,q,0,wpp.⁢

Przypomnijmy, że obrót Oik,jk jest wybrany tak, że aqk=0. Zauważmy, że dla

z→:=D-1⁢⁢a→k-1=Sik-1,jk-1⁢⁢⋯⁢⁢Si1,j1⁢⁢a→

mamy

a→k⁢=⁢Op,q⁢⁢a→k-1⁢=⁢D⌃⁢⁢Sp,q⁢⁢D-1⁢⁢D⁢⁢z→⁢=⁢D⌃⁢⁢Sp,q⁢⁢z→.

Stąd warunek aqk=0 jest równoważny warunkowi Sp,q⁢⁢z→q=0, czyli

s⁢⁢zp⁢dp=c⁢⁢zq⁢dq.

Jeśli teraz zq=0 to przyjmujemy D⌃:=D i Sp,q:=I. W przeciwnym przypadku, uwzględniając równość s2+c2=1, dostajemy

c2=1/t,⁢t:=1+dq2⁢zq2dp2⁢zp2,⁢⁢s2=1/t′,⁢t′:=1+dp2⁢zp2dq2⁢zq2.

Mamy teraz dwa przypadki.

(i) Jeśli 0<zq2⁢dq2≤zp2⁢dp2 to kładziemy d⌃p=c⁢dp, d⌃q=c⁢dq i d⌃l=dl dla l≠p,q. Wtedy sl,l=1 dla l≠p,q oraz

sp,psp,qsq,psq,q=1α-β1,

gdzie

α=zq⁢dq2zp⁢dp2,⁢β=zqzp.

(ii) Jeśli 0≤zp2⁢dp2<zq2⁢dq2 to kładziemy d⌃p=s⁢dq, d⌃q=s⁢dp i d⌃l=dl dla l≠p,q. Wtedy sl,l=1 dla l≠p,q oraz

sp,psp,qsq,psq,q=α′1-1β′,

gdzie

α′=zp⁢dp2zq⁢dq2,⁢β′=zpzq.

Macierz Sp,q ma podobną strukturę jak Op,q. Obecność dwóch jedynek pozwalaja jednak zaoszczędzić dwa mnożenia przy operacji Sp,q⁢⁢x→ w porównaniu do Op,q⁢⁢x→. Ponadto, co jest ważniejsze, obliczanie wielkości α i β (ew. α′ i β′) definiujących Sp,q nie wymaga pierwiastkowania. Rzeczywiście, w kolejnych krokach procesu obliczeniowego wystarczy pamiętać współczynniki d⌃l2. Współczynniki te wyliczane są rekurencyjnie zgodnie ze wzorami d⌃l2=dl2 dla l≠p,q, a dla l=p,q

d⌃l2=⁡dl2/tjeśli (i),dl2/t′jeśli (ii).⁢

Zauważmy, jeszcze, że takie rozbicie na przypadki (i) i (ii) pozwala nie tylko uniknąć dzielenia przez zero, ale również uniknąć możliwego niedomiaru. Wielkość t jeśli zachodzi (i) oraz t′ jeśli zachodzi (ii) jest w przedziale 1/2,1.

Możnaby się spytać gdzie jest zysk obliczeniowy. Na końcu procesu dysponujemy bowiem jedynie macierzą Dm2 i aby pomnożyć przez Dm należy najpierw wykonać n pierwiastkowań. Nie jest to jednak zawsze prawdą.

Na przykład, jeśli chcemy przekształcić dany wektor a→ na kierunek wersora e→1 przy pomocy obrotów Givensa to po wykonaniu serii m=n-1 przekształceń

S1,n⁢⁢⋯⁢⁢S1,3⁢⁢S1,2⁢⁢a→

jedynie pierwsza współrzędna otrzymanego wektora nie jest zerowa. Mnożenie przez Dm sprowadza się więc do operacji na jego pierwszej współrzędnej, to zaś wymaga obliczenia tylko jednego pierwiastka kwadratowego.

Podobnie, używając obrotów Givensa można kosztem porównywalnym do algorytmu bazującego na odbiciach Householdera przekształcić daną nieosobliwą macierz kwadratową A do postaci trójkątnej górnej (albo, równoważnie, znaleźć jej rozkład ortogonalno-trójkątny). Rzeczywiście, dokonuje się tego wykonując najpierw na macierzy serię m=n-1⁢n/2 przekształceń

Sn-1,n⁢⁢Sn-2,n⁢⁢Sn-2,n-1⁢⁢⋯⁢⁢S2,n⁢⁢⋯⁢⁢S2,3⁢⁢S1,n⁢⁢⋯⁢⁢S1,2⁢⁢A,

gdzie Si,j jest tak dobrana, aby wyzerować współrzędną j,i aktualnej macierzy, a następnie mnożąc kolejne wiersze otrzymanej macierzy trójkątnej górnej przez pierwiastki z kolejnych elementów diagonalnych macierzy Dm2. Całość wymaga więc obliczenia n pierwiastków kwadratowych, podobnie jak w algorytmie Householdera. (Należy przy tym pamiętać, aby nie wykonywać niepotrzebnie operacji na elementach zerowych!)

4.2. Metoda Jacobiego

Metoda Jacobiego jest najstarszą z istniejących metod numerycznych znajdowania wartości własnych macierzy symetrycznych. Ma jednak nie tylko znaczenie historyczne. Chociaż jest stosunkowo wolno zbieżna to można ją stosować gdy chcemy obliczyć wartości własne z bardzo dobrą dokładnością w obecności błędów zaokrągleń.

Zakładamy, że A∈Rn,n jest macierzą symetryczną. Podobnie jak w metodzie QR, metoda Jacobiego startuje z macierzy A0=A i produkuje ciąg macierzy podobnych Ak dla k=1,2,3,…, który zbiega do macierzy diagonalnej

Λ=diag⁢λ1,λ2,…,λn

składającej się z wartości własnych A. Różnica polega na tym, że pojedyncze kroki metody wykonuje się używając obrotów Givensa,

	Ak	=	Oik,jkT⁢⁢Ak-1⁢⁢Oik,jk⁢=⁢Oik,jkT⁢⁢Oik-1,jk-1T⁢⁢Ak-2⁢⁢Oik,jk⁢⁢Oik-1,jk-1
		=	…⁢=⁢Oik,jkT⁢⁢⋯⁢⁢Oi1,j1T⁢⁢A0⁢⁢Oi1,j1⁢⁢⋯⁢⁢Oik,jk.

Oznaczając Ok=Oi1,j1⁢⁢⋯⁢⁢Oik,jk mamy więc Ak=OkT⁢⁢A⁢⁢Ok.

Jako miarę odchylenia od macierzy diagonalnej Λ w naturalny sposób przyjmiemy połowę sumy kwadratów elementów pozadiagonalnych, czyli

ωk:=ω⁢Ak=∑i<jai,jk2,

gdzie ai,jki,j=1n. Wielkość ωk dobrze charakteryzuje również błąd przybliżenia wartości własnych macierzy A przez elementy diagonalne macierzy Ak. Rzeczywiście, niech

Dk=diag⁢a1,1k,…,an,nk

i Rk=Ak-Dk. Wtedy

Dk=OkT⁢⁢A⁢⁢Ok-OkT⁢⁢Ok⁢⁢Rk⁢⁢OkT⁢⁢Ok=OkT⁢⁢A+Ek⁢⁢Ok,

gdzie Ek=-Ok⁢⁢Rk⁢⁢OkT oraz

EkF=RkF=2⁢ωk.

Wyrazy diagonalne macierzy Dk są więc wartościami własnymi macierzy A+Ek, którą możemy potraktować jako zaburzoną macierz A. Korzystając z twierdzenia Weyla 1.3 dostajemy, że istnieje takie uporządkowanie λp⁢1,…,λp⁢n wartości własnych macierzy A, że dla wszystkich i mamy

ai,ik-λp⁢i≤Ek2≤EkF=2⁢ωk.

Stąd wniosek, że poszczególne obroty Oik,jk powinniśmy wybierać tak, aby maksymalizować różnicę ωk-1-ωk. Aby odpowiedzieć na pytanie który obrót osiąga ten cel, przeanalizujmy pojedynczy krok. Dla prostoty zapisu oznaczymy p:=ik, q:=jk, A:=Ak-1, A⌃:=Ak=Op,qT⁢⁢A⁢⁢Op,q, gdzie Op,q jest obrotem o kąt φ, c=cos⁡φ, s=sin⁡φ, c2+s2=1.

Zauważmy, że mnożenie z lewej strony przez obrót Op,qT zmienia jedynie wiersze p-ty i q-ty macierzy, a mnożenie z prawej strony przez Op,q zmienia kolumny p-tą i q-tą. Dlatego dla i≠p,q i j≠p,q wyrazy a⌃i,j i ai,j są sobie równe. Jeśli zaś i≠p,q i j=p,q to

a⌃i,p=a⌃p,i=c⁢ai,p-s⁢ai,q,⁢⁢a⌃i,q=a⌃q,i=s⁢ai,p+c⁢ai,q,

a stąd

a⌃i,p2+a⌃i,q2=ai,p2+ai,q2

i podobnie a⌃p,i2+a⌃q,i2=ap,i2+aq,i2. Otrzymujemy

ωk-1-ωk=ω⁢A-ω⁢A⌃=ap,q2-a⌃p,q2.

Dla danych p<q wielkości c i s definiujące obrót należy więc wybrać tak, aby

a⌃p,q=0.

Jeśli ap,q=0 to powyższą równość realizuje obrót zerowy. Załóżmy więc, że ap,q≠0. Bezpośredni rachunek pokazuje, że wtedy

a⌃p,pa⌃p,qa⌃q,pa⌃q,q=c2⁢ap,p-2⁢s⁢c⁢ap,q+s2⁢aq,qc2-s2⁢ap,q+s⁢c⁢ap,p-aq,qc2-s2⁢ap,q+s⁢c⁢ap,p-aq,qs2⁢ap,p+2⁢s⁢c⁢ap,q+c2⁢aq,q.

Stąd równość a⌃p,q=0 zachodzi wtedy gdy

c2-s2⁢ap,q+s⁢c⁢ap,p-aq,q=0.

Jeśli podstawimy t=s/c=φ\tg to ostatnia równość jest równoważna równaniu kwadratowemu

t2+2⁢⁢β⁢⁢t-1⁢= 0,⁢⁢β=aq,q-ap,p2⁢ap,q.

Rozwiązując otrzymujemy następujące wzory:

c=11+t2,⁢s=t⋅c,⁢gdzie⁢⁢t=sign⁢ββ+1+β2.

Zauważmy jeszcze, że taka postać rozwiązania pozwala obliczyć s i c z małym błędem względnym w arytmetyce maszyny cyfrowej.

Pozostaje do ustalenia jak wybierać kolejne p,q. W oryginalnej wersji algorytmu p,q wybierane są tak, że

ap,q=maxi<j⁡ai,j.

Ponieważ wtedy ap,q2≥ω⁢A/N, gdzie N=n⁢n-1/2 to

ω⁢A⌃=ω⁢A-ap,q2≤ω⁢A⁢1-1/N.

(4.1)

Mamy więc zapewnioną zbieżność liniową, chociaż iloraz zbieżności jest bliski jedynce dla dużych wymiarów n. Taki wybór jest jednak zbyt kosztowny, bo wymaga znalezienia w każdym kroku elementu największego co do modułu wśród N elementów. Dlatego stosowane są inne strategie. Na przykład,

1,2,1,3,2,3,1,4,2,4,3,4,…,1,i,…,i-1,i,…,1,n,…,n-1,n,

i dalej cyklicznie, aż do osiągnięcia żądanej dokładności. W tym przypadku kłopot polega na tym, że obroty wykonuje się również wtedy, gdy ap,q2 jest znacznie mniejsze od ω⁢A, co nie jest opłacalne. Aby temu zapobiec, można wprowadzić pewien próg δ i nie wykonywać obrotów gdy ap,q2<δ. Na przykład, możemy wziąć δ=ω⁢A/N. Wtedy, po wykonaniu każdego kroku spełniona jest nierówność (4.1).

4.3. QR dla macierzy trójdiagonalnej

Metodę Jacobiego stosuje się do pełnej, rzeczywistej macierzy symetrycznej. W tym podrozdziale zatrzymamy się na przypadku, gdzie A=AT∈Rn,n jest dodatkowo trójdiagonalna,

A=T=a1b1b1⋱⋱⋱⋱bn-1bn-1an.

Przypomnijmy, że każdą macierz symetryczną można sprowadzić do postaci trójdiagonalnej przez podobieństwa ortogonalne. W rozdziale 2.1.1 pokazaliśmy jak to zrobić przy pomocy odbić Householdera. Ale równie dobrze możemy użyć obrotów Givensa. Rzeczywiście, najpierw wybieramy obroty O2,i dla i=3,4,…,n tak, aby mnożenie

O2,nT⋯O2,4TO2,3TA=:A⌃

spowodowało wyzerowanie kolejno elementów 3,1,4,1,…,n,1. (Pamiętamy o algorytmie Gentlemana!) Ponieważ mnożenie z prawej strony przez O2,i zmienia jedynie kolumny 2-gą oraz i-tą to po operacji

A⌃⁢⁢O2,3⁢⁢O2,4⁢⁢⋯⁢⁢O2,n

wyzerowane elementy w pierwszej kolumnie się nie zmienią. Dalej postępujemy indukcyjnie dla kolejnych kolumn, od 2-giej do n-2-giej, i ostatecznie dostajemy macierz trójdiagonalną T=OT⁢⁢A⁢⁢O, gdzie

O⁢=⁢O2,3⁢⁢⋯⁢⁢O2,n⁢⁢O3,4⁢⁢⋯⁢⁢O3,n⁢⁢⋯⁢⁢On-2,n-1⁢⁢On-2,n⁢⁢On-1,n,

przy czym mnożenie z lewej przez Oi,jT zeruje element i-1,j.

Oczywiście, wszystkie metody dla macierzy pełnej można stosować również dla macierzy trójdiagonalnej. Jednak ze względu na prostszą strukturę macierzy niektóre metody znacznie się upraszczają. Pokażemy to na przykładzie metody QR z rozdziału 3.2.

W oryginalnej metodzie QR (z przesunięciami lub bez) rozkład ortogonalno-trójkątny macierzy A=Q⁢⁢R dokonuje się przy pomocy odbić Householdera. Dla macierzy trójdiagonalnej T wygodnie jest użyć do tego celu obrotów Givensa Oi,i+1T zerujących kolejno elementy i+1,i dla i=1,2,…,n-1. Wtedy

R⁢=⁢On-1,nT⁢⁢On-2,n-1T⁢⁢⋯⁢⁢O2,3T⁢⁢O1,2T⁢⁢T

jest macierzą trójkątną górną. Teraz istotne jest, że po wykonaniu drugiego kroku metody QR, czyli mnożenia

O1,2⁢⁢O2,3⁢⁢⋯⁢⁢On-1,n⁢⁢R

otrzymana macierz T⌃ pozostaje trójdiagonalna. Rzeczywiście, skoro R jest trójkątna górna, a mnożenie z lewej strony przez Oi,i+1 kombinuje jedynie wiersze i-ty i i+1-szy to T⌃ jest macierzą górną Hessenberga, a ponieważ jest również symetryczna, to jest trójdiagonalna.

Koszt jednego kroku metody QR dla symetrycznej macierzy trójdiagonalnej T jest więc proporcjonalny do n.

Mamy własność ogólniejszą. W procesie iteracyjnym metody QR macierz trójdiagonalna pozostaje w każdym kroku trójdiagonalna niezależnie od zastosowanego algorytmu rozkładu na iloczyn ortogonalno-trójkątny. Wynika to w szczególności stąd, że, jak pamiętamy, rozkład taki jest wyznaczony jednoznacznie z dokładnością do macierzy diagonalnej z elementami ±1. (Jako ćwiczenie, można to pokazać bezpośrednio w przypadku gdy rozkład dokonywany jest przy użyciu odbić Householdera.)

Oczywiście wszystkie pozostałe fakty na temat metody QR pokazane wcześniej dla macierzy pełnych, w szczególności te dotyczące wyboru przesunięcia i szybkości zbieżności, pozostają w mocy gdy macierz jest trójdiagonalna.

4.4. Rozkład według wartości szczególnych (SVD)

Na końcu serii wykładów na temat zadania własnego zajmiemy się znajdowaniem wartości szczególnych danej macierzy, a to dlatego, że zadanie to i metody jego rozwiązania wiążą się ściśle z zadaniem i metodami znajdowania wartości własnych macierzy symetrycznych.

4.4.1. Twierdzenie o rozkładzie

Najpierw przypomnimy twierdzenie o rozkładzie według wartości szczególnych, czyli SVD (ang. Singular Value Decomposition). Będziemy zakładać, bez zmniejszenia ogólności, że liczba kolumn macierzy A nie przekracza liczby jej wierszy, co najczęściej występuje w praktyce. W przeciwnym przypadku wystarczy zastosować poniższe twierdzenie dla macierzy transponowanej AT.

Twierdzenie 4.1

Dla dowolnej macierzy A∈Rm,n, gdzie m≥n, istnieją macierze ortogonalne U∈Rm,m i V∈Rn,n takie, że

A⁢=⁢U⁢⁢Σ⁢⁢VT,

gdzie Σ∈Rm,n jest macierzą diagonalną,

Σ=Σ⌃0,⁢Σ⌃=diag⁢σ1,σ2,…,σn∈Rn,n,

σ1≥σ2≥⋯≥σk>σk+1=⋯=σn=0, a k jest rzędem macierzy A.

Dowód

Ponieważ macierz AT⁢⁢A jest symetryczna i nieujemnie określona to istnieje baza ortonormalna v→ii=1n jej wektorów własnych, a odpowiadające jej wartości własne są nieujemne (patrz twierdzenia 1.1). Oznaczmy te wartości własne odpowiednio przez σi2,

σ1≥⋯≥σk>σk+1=…=σn=0,

przy czym k=rank⁢AT⁢⁢A=rank⁢A. Zauważmy, że

A⁢⁢v→iT⁢⁢A⁢⁢v→j=v→iT⁢⁢AT⁢⁢A⁢⁢v→j=σj2⋅v→iT⁢⁢v→j=⁡0i≠j,σj2i=j.⁢

To oznacza, że wektory u→i:=A⁢⁢v→i/σi, 1≤i≤k, tworzą układ ortonormalny i można go uzupełnić wektorami u→k+1,…,u→m do bazy ortonormalnej w Rm.

Zdefiniujmy teraz macierze ortogonalne

V=V1,V2∈Rn,n,⁢V1=v→1,…,v→k∈Rn,k,⁢V2=v→k+1,…,v→n∈Rn,n-k,

U=U1,U2∈Rm,m,⁢U1=u→1,…,u→k∈Rm,k,⁢U2=u→k+1,…,u→n∈Rm,m-k.

Wtedy, oznaczając Σ′=diag⁢σ1,…,σk∈Rk,k, mamy

UT⁢⁢A⁢⁢V=U1TU2T⁢⁢A⁢⁢V1,A⁢⁢V2=U1TU2T⁢⁢U1⁢⁢Σ′,0=Σ′000=Σ,

albo, równoważnie, A=U⁢⁢Σ⁢⁢VT.

∎

Wielkości σi, 1≤i≤n, to właśnie wartości szczególne macierzy A. Są one wyznaczone jednoznacznie. Rzeczywiście, jeśli mamy dwa rozkłady, U1⁢⁢Σ1⁢⁢V1T=A=U2⁢⁢Σ2⁢⁢V2T, to

V⁢⁢Σ1T⁢⁢Σ1=Σ2T⁢⁢Σ2⁢⁢V,⁢V=V2T⁢⁢V1.

Porównując wyrazy diagonalne po obu stronach powyższej równości dostajemy σi1=σi2 dla wszsytkich i.

Zanotujmy jeszcze, że mając rozkład SVD można np. łatwo rozwiązać liniowe zadanie najmniejszych kwadratów, tzn. znaleźć wektor x→* minimalizujący b→-A⁢⁢x→2 po wszystkich x→∈Rn. (Jeśli istnieje wiele takich wektorów to bierzemy ten o najmniejszej normie.) Mamy bowiem

b→-A⁢⁢x→2=U⁢⁢UT⁢⁢b→-U⁢⁢Σ⁢⁢VT⁢⁢x→2=c→-Σ⁢⁢y→2,

c→=UT⁢⁢b→, y→=VT⁢⁢x→. Stąd x→*=VT⁢⁢y→*, gdzie

yi*=⁡ci/σi,1≤i≤k,0,1≤i≤n.⁢

4.4.2. Dlaczego nie pomnożyć AT⁢⁢A?

Z dowodu twierdzenia 4.1 wynika, że σ12,…,σn2 są wartościami własnymi macierzy symetrycznej i nieujemnie określonej AT⁢⁢A, a kolumny macierzy ortogonalnej V tworzą odpowiednio bazę ortonormalną w Rn wektorów własnych tej macierzy. Podobnie, σ12,…,σn2,0,…,0︸m-n są wartościami własnymi A⁢⁢AT, a kolumny U tworzą bazę ortonormalną w Rm wektorów własnych. Wydaje się więc, że wartości szczególne (i w razie potrzeby cały rozkład SVD) można łatwo wyznaczyć wykonując najpierw mnożenie A⌃:=AT⁢⁢A lub A⌃:=A⁢⁢AT, a następnie aplikując do macierzy A⌃ jedną z rozpatrzonych wcześniej metod dla zadania własnego. Należy jednak przestrzec przed takim mechanicznym działaniem.

Przykład 4.1

Niech A=110ε. Dla ,,małych” ε wartości szczególne tej macierzy są bliskie 2 i ε/2. Jeśli ε jest na tyle małe, że 1+ε2 jest w arytmetyce fl reprezentowane przez 1 to macierz AT⁢⁢A=1111+ε2 jest reprezentowana przez macierz 1111, która jest już osobliwa - jej druga wartość szczególna jest zerowa.

Z uwagi na możliwą zmianę rzędu macierzy, jak w przytoczonym przykładzie, stosuje się nieco zmodyfikowane algorytmy znajdowania wartości własnych, które unikają jawnych mnożeń AT⁢⁢A lub A⁢⁢AT. Pokażemy to najpierw, nie wdając się w szczegóły, na przykładzie metody Jacobiego z rozdziału 4.2.

Przypomnijmy, że metoda Jacobiego zastosowana bezpośrednio do macierzy A⌃=AT⁢⁢A konstruuje ciąg macierzy podobnych A⌃=A⌃0,A⌃1,…,A⌃k,…, gdzie A⌃k=Oik,jkT⁢⁢A⌃k-1⁢⁢Oik,jk, a Oik,jk jest tak dobranym obrotem Givensa, aby wyzerować element ik,jk i jednocześnie, wobec symetrii, także jk,ik. Oznaczmy A0=A i Ak=Ak-1⁢⁢Oik,jk dla k≥1. Wtedy A⌃k=AkT⁢⁢Ak. Pomysł polega na tym, aby zamiast obliczać jawnie A⌃k=Oik,jkT⁢⁢Ak-1T⁢⁢Ak-1⁢⁢Oik,jk, w kolejnych krokach obliczać i pamiętać jedynie Ak=Ak-1⁢⁢Oik,jk. Jest to możliwe, bowiem wyznaczenie rotacji Oik,jk wymaga jedynie znajomości elementów

a⌃ik,ikk-1=∑l=1mal,ikk-12,⁢a⌃jk,jkk-1=∑l=1mal,jkk-12,⁢a⌃ik,jkk-1=∑l=1mal,ikk-1⁢al,jkk-1,

które można obliczyć korzystając z powyższych wzorów.

4.4.3. SVD dla macierzy dwudiagonalnych

W tym podrozdziale pokażemy algorytm, który można traktować jako wariant metody QR zastosowanej do macierzy AT⁢⁢A. Zanim jednak przejdziemy do samego algorytmu, poczynimy kilka pomocniczych uwag teoretycznych.

Niech T0 będzie macierzą symetryczną i dodanio okteśloną, T0=T0T>0. Rozpatrzmy proces iteracyjny, nazwijmy go LR, który startuje z macierzy T0 i w kolejnych krokach k=1,2,…

(i) wybiera przesunięcie τk, mniejsze od najmniejszej wartości własnej Tk-1 (np. τk=0),
(ii) dokonuje rozkładu Banachiewicza-Cholesky'ego

Tk-1-τk2⋅I=BkT⁢⁢Bk,

gdzie Bk jest macierzą trójkątną górną z dodatnimi elementami na głównej przekątnej,
(iii) produkuje Tk=Bk⁢⁢BkT+τk2⋅I.

(Przypomnijmy, że rozkładu Banachiewicza-Cholesky'ego dokonujemy zmodyfikowanym algorytmem eliminacji Gaussa.) Oczywiście, macierze Tk są do siebie podobne, bo

Tk=Bk⁢⁢BkT+τk2⋅I=Bk⁢⁢Tk-1-τk2⋅I⁢⁢Bk-1+τk2⋅I=Bk⁢⁢Tk-1⁢⁢Bk-1.

Zachodzi ciekawsza własność.

Lemat 4.1

Dwa kroki iteracji LR z tym samym przesunięciem τ są równoważne jednemu krokowi iteracji QR z przesunięciem τ.

Dowód

Bez zmniejszenia ogólności przyjmijmy, że k=0 i τ1=τ2=τ. Z jednej strony, z rozkładu QR mamy T0-τ2⋅I2=Q⁢⁢RT⁢⁢Q⁢⁢R=RT⁢⁢R, gdzie elementy na głównej przekątnej macierzy R są dodatnie, a z drugiej z rozkładu LR

	T0-τ2⋅I2	=	B1T⁢⁢B1⁢⁢B1T⁢⁢B1⁢=⁢B1T⁢⁢T1-τ2⋅I⁢⁢B1
		=	B1T⁢⁢B2T⁢⁢B2⁢⁢B1⁢=⁢B2⁢⁢B1T⁢⁢B2⁢⁢B1.

Wobec jednoznaczności rozkładu Banachiewicza-Cholesky'ego mamy więc R=B2⁢⁢B1. Stąd macierz powstała w wyniku jednego kroku iteracji QR wynosi

T⌃	=	R⁢⁢Q+τ2⋅I⁢=⁢R⁢⁢Q⁢⁢R⁢⁢R-1+τ2⋅I
	=	R⁢⁢T0-τ2⋅I⁢⁢R-1+τ2⋅I⁢=⁢R⁢⁢T0⁢⁢R-1
	=	B2⁢⁢B1⁢⁢B1T⁢⁢B1+τ2⋅I⁢⁢B2⁢⁢B1-1
	=	B2⁢⁢B1⁢⁢B1T⁢⁢B1⁢⁢B1-1⁢⁢B2-1+τ2⁢B2⁢⁢B1⁢⁢B2⁢⁢B1-1
	=	B2⁢⁢B1⁢⁢B1T⁢⁢B2-1+τ2⋅I⁢=⁢B2⁢⁢T1-τ2⋅I⁢⁢B2-1+τ2⋅I
	=	B2⁢⁢B2T⁢⁢B2⁢⁢B2-1+τ2⋅I⁢=⁢B2⁢⁢B2T+τ2⋅I
	=	T2-τ2⋅I+τ2⋅I⁢=⁢T2,

jak w tezie lematu.

∎

Powyższy lemat pokazuje, że rozważania teoretyczne dotyczące iteracji QR można przenieść na iteracje LR. Dlatego dalej nie będziemy się już zajmować analizą teoretyczną LR, a przejdziemy do opisu algorytmu te iteracje wykorzystującego.

Zakładamy, że A∈Rm,n jest kolumnowo regularna, czyli rank⁢A=n, oraz dwudiagonalna. Dokładniej, ai,j=0 dla i≥j+1 i dla j≥i+2,

A=a1b1a2⋱⋱bn-1an.

Oczywiście, wobec kolumnowej regularności macierzy mamy aj≠0. Założenie o dwudiagonalności wydaje się mocno ograniczające. W istocie jednak tak nie jest, bo każdą macierz można sprowadzić do takiej postaci nie zmieniając jej wartości szczególnych i kontrolując odpowiednie wektory własne, co pokażemy na samym końcu.

Algorytm działa podobnie jak iteracje LR z przesunięciami, ale zamiast tworzyć jawnie macierze Tk, podobne do AT⁢⁢A, pracuje bezpośrednio na macierzach Bk. Przyjmujemy B0=D⁢⁢A, gdzie D∈Rn,m jest macierzą diagonalną z elementami na głównej przekątnej dj=sign⁢aj, 1≤j≤n. Wtedy B0∈Rn,n jest macierzą dwudiagonalną z dodatnimi elementami na głównej przekątnej oraz

T0=B0T⁢⁢B0=AT⁢⁢DT⁢⁢D⁢⁢A=AT⁢⁢A

jest macierzą trójdiagonalną. Zobaczmy teraz jak mając macierz dwudiagonalną Bk możemy skonstruować Bk+1. Ponieważ Tk jest trójdiagonalna to Bk+1 musi być dwudiagonalna. (W szczególności, wyniknie to również z dalszych rachunków.) Dla uproszczenia zapisu, oznaczmy elementy przekątniowe macierzy Bk i Bk+1 odpowiednio przez aj i a⌃j, 1≤j≤n, a elementy nad główną przekątną przez bj i b⌃j, 1≤j≤n-1. Przyjmijmy dodatkowo b0=b⌃0=bn=b⌃n=0. Macierze Bk i Bk+1 związane są równaniem

Bk+1T⁢⁢Bk+1+τk+12⋅I⁢=⁢Tk⁢=⁢Bk⁢⁢BkT+τk2⋅I.

Porównując wyrazy przekątniowe po obu stronach tej równości otrzymujemy

a⌃j2+b⌃j-12+τk+12=aj2+bj2+τk2,⁢1≤j≤n,

a porównując wyrazy nad przekątną

a⌃j⁢b⌃j=aj+1⁢bj2,⁢1≤j≤n-1.

Stąd, podstawiając δ=τk+12-τk2, mamy następujące wzory na a⌃j i b⌃j. Dla j=1,2,…,n-1 obliczamy

a⌃j:=aj2+bj2-b⌃j-12-δ,⁢b⌃j:=aj+1⁢bj/a⌃j

i na końcu a⌃n:=an2-b⌃n-12-δ.

Ponieważ obliczanie pierwiastków jest stosunkowo kosztowne, wzory te opłaca się zmodyfikować tak, aby nie obliczać pierwiastków w każdej iteracji, a jedynie na końcu całego procesu iteracyjnego. Można to osiągnąć pracując na kwadratach aj i bj. Rzeczywiście, wprowadzając zmienne pj=aj2 i qj=bj2, otrzymujemy następującą procedurę dla jednego kroku iteracyjnego:

.	begin
.	.	for j=1 to n-1 do
.	.	begin
.	.	.	p⌃j:=pj+qj-q⌃j-δ;
.	.	.	q⌃j:=qj⋅qj+1/q⌃j
.	.	end;
.	.	p⌃n:=pn-q⌃n-1-δ
.	end;

Koszt jednego kroku iteracyjnego jest stały. Stąd, jeśli macierz wyjściowa A jest dwudiagonalna to całkowity koszt algorytmu jest proporcjonalny do liczby iteracji.

A co jeśli A nie jest dwudiagonalna? Wtedy trzeba ją wstępnie przekształcić do postaci dwudiagonalnej nie zmieniając wartości szczególnych. W tym celu możemy użyć np. odbić Householdera. Najpierw zerujemy wyrazy w pierwszej kolumnie, oprócz wyrazu diagonalnego, poprzez pomnożenie macierzy A z lewej strony przez odbicie HL,1∈Rm,m przekształcające pierwszą kolumnę A na kierunek e→1. Oznaczmy A~1=a~i,j=HL,1⁢⁢A. Następnie wybieramy odbicie H~R,1∈Rn-1,n-1 tak, aby przekształcało wektor a~1,2,a~1,3⁢…,a~1,nT∈Rn-1 na kierunek e→1∈Rn-1 i mnożymy A~1 z prawej przez HR,1T=10→T0→H~R,1T∈Rn,n. Ponieważ to mnożenie nie zmienia pierwszej kolumny A~1, w powstałej macierzy A1=HL,1⁢⁢A⁢⁢HR,1T elementy w pierwszej kolumnie i w pierwszym wierszu, poza 1,1 i 1,2, są wyzerowane. Dalej postępujemy indukcyjnie zerując na przemian odpowiednie elementy w kolumnach i wierszach. Ostatecznie, po n-1 krokach otrzymujemy macierz

A⌃=U1⁢⁢A⁢⁢V1T,⁢gdzie⁢⁢U1=HL,n-1⁢⁢⋯⁢⁢HL,1,⁢V1=HR,n-2⁢⁢⋯⁢⁢HR,1,

która jest już dwudiagonalna. Jeśli teraz A⌃=U2⁢⁢Σ⁢⁢V2T to A=U1T⁢⁢U2⁢⁢Σ⁢⁢V1T⁢⁢V2T, a więc A i A⌃ mają te same wartości szczególne. Ponadto,

A=U⁢⁢Σ⁢⁢VT,⁢gdzie⁢⁢U=U1T⁢⁢U2,⁢V=V1T⁢⁢V2.

Zanotujmy jeszcze, że dowolną macierz można sprowadzić w podobny sposób do postaci dwudiagonalnej przy pomocy obrotów Givensa. Zarówno w przypadku zastosowania obrotów Givensa jak i odbić Householdera koszt jest proporcjonalny do m⋅n2.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Matematyka obliczeniowa II

Zagadnienia

4. Zagadnienie własne IV

4.1. Obroty Givensa

4.1.1. Definicja obrotu

4.1.2. Algorytm Gentlemana

4.2. Metoda Jacobiego

4.3. QR dla macierzy trójdiagonalnej

4.4. Rozkład według wartości szczególnych (SVD)

4.4.1. Twierdzenie o rozkładzie

Twierdzenie 4.1

Dowód

4.4.2. Dlaczego nie pomnożyć AT⁢⁢A?

Przykład 4.1

4.4.3. SVD dla macierzy dwudiagonalnych

Lemat 4.1

Dowód