Zagadnienia

3.1 Iteracje podprzestrzeni
- 3.1.1 Algorytm ogólny
- 3.1.2 Iteracje ortogonalne
3.2 Metoda QR
3.3 QR z przesunięciami i deflacja

3. Zagadnienie własne III

Metodę potęgową można traktować jako iteracje podprzestrzeni jednowymiarowej startujące ze span⁢x→0. Naturalnym uogólnieniem tego procesu są iteracje podprzestrzeni o wyższych wymiarach. Właśnie iteracje podprzestrzeni są punktem wyjścia konstrukcji obecnie najpopularniejszego algorytmu QR, którego zadaniem jest obliczenie jednocześnie wszystkich wartości własnych.

Chociaż metody prezentowane w tym rozdziale mogą być stosowane w większej ogólności, dla uproszczenia będziemy zakładać, że macierz A jest nieosobliwa, rzeczywista i symetryczna, tzn. A=AT∈Rn,n i det⁢A≠0. Przypomnijmy, że wtedy istnieje baza ortonormalna ξ→ii=1n wektorów własnych macierzy,

A⁢⁢ξ→i=λi⁢⁢ξ→i,⁢1≤i≤n.

Będziemy również zakładać, że wartości własne są uporządkowane tak, że

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

3.1. Iteracje podprzestrzeni

3.1.1. Algorytm ogólny

Niech 1≤p≤n-1 oraz Y0⊆Rn będzie podprzestrzenią o wymiarze p. Dla k=1,2,3,… rozpatrzmy następujące iteracje:

Yk⁢:=⁢A⁢Yk-1⁢=⁢⁢A⁢⁢x→:⁢x→∈Yk-1⁢.

Oczywiście, wobec nieosobliwości A wszystkie podprzestrzenie Yk mają też wymiar p. Pokażemy, że przy pewnych niekrępujących założeniach ciąg podprzestrzeni Yk zbiega do podprzestrzeni własnej

Pp:=span⁢ξ→1,ξ→2,…,ξ→p,

przy czym przez ,,zbieżność” rozumiemy tu zbieżność do zera kąta pomiędzy tymi podprzestrzeniami.

Przypomnijmy, że kąt pomiędzy dwiema podprzestrzeniami Y,Z⊆Rn definiujemy jako

∠⁢Yk,Pp:=maxy→∈Yk⁡minz→∈Pp⁡∠⁢y→,z→.

(Jako ćwiczenie można wykazać, że jeśli dim⁢Y=dim⁢Z to ∠⁢Y,Z=∠⁢Z,Y.)

Twierdzenie 3.1

Niech

ρp:=λp+1λp<1⁢⁢ oraz ⁢⁢γk:=∠⁢Yk,Pp,⁢k=0,1,2,….

Jeśli γ0<π/2 to

γk⁢≤⁢ρp⋅γk-1⁢≤⁢ρpk⋅γ0.

Dowód

Niech y→ będzie wektorem spełniającym ∠⁢y→,Pp<π/2. Wtedy

y→=∑i=1nv→i,⁢v→i∈span⁢ξ→i,⁢1≤i≤n,

gdzie nie wszystkie v→i dla 1≤i≤p, są są zerowe, oraz A⁢⁢y→=∑i=1nλi⁢v→i. Stąd

∠⁢A⁢⁢y→,Pp=∑i=p+1nλi2⁢v→i22∑i=1pλi2⁢v→i22≤λp+1λp⋅∑i=p+1nv→i22∑i=1pv→i22=ρp⋅∠⁢y→,Pp⁢p.

(3.1)

Jeśli teraz γk+1=∠⁢y→k+1,Pp dla y→k+1∈Yk+1 oraz y→k+1=A⁢⁢y→k to z (3.1) wynika, że

γk+1≤ρp⋅∠⁢y→k,Pp≤ρp⋅γk,

co kończy dowód.

∎

3.1.2. Iteracje ortogonalne

W praktyce, iteracje podprzestrzeni realizujemy reprezentując kolejne Yk przez ich bazy ortonormalne. Wtedy algorytm, zwany w tym przypadku również iteracjami ortogonalnymi, przybiera następującą postać. Wybieramy macierz kolumnowo-ortogonalną Z0 formatu n×p jako macierz startową, a następnie dla k=1,2,3,… iterujemy:

(IO)⁢⁢⁡Yk⁢:=⁢A⁢⁢Zk-1;Yk⁢:=⁢Zk⁢⁢Rk;⁢(ortonormalizacja)⁢

przy czym w drugiej linii następuje ortonormalizacja macierzy Yk realizowana poprzez rozkład tej macierzy na iloczyn macierzy kolumnowo-ortonormalnej Zk formatu n×p i macierzy trójkątnej górnej Rk formatu p×p. Dokonuje się tego znanym nam już algorytmem wykorzystującym odbicia Householdera.

łatwo zauważyć, że jeśli przyjmiemy, iż wektory-kolumny macierzy Z0 tworzą bazę ortonormalną podprzestrzeni Y0 to w kolejnych krokach wektory-kolumny Zk tworzą bazę ortonormalną podprzestrzeni Yk. Rzeczywiście, ponieważ Yk=A⁢⁢Zk-1 to kolumny Yk stanowią bazę Yk. Kolumny te są z kolei liniową kombinacją kolumn Zk, czyli rozpinają tą samą podprzestrzeń co kolumny Zk.

Poczynimy teraz kluczową obserwację. Przypuśćmy, że wszystkie wartości własne λi są parami różne. Gdybyśmy rozpoczynali iteracje od macierzy Y¯0 składającej się jedynie z p¯ początkowych kolumn macierzy Y0, przy czym 1≤p¯<p, to otrzymane w procesie iteracyjnym macierze Y¯k i Z¯k (formatu n×p¯ i p¯×p¯) byłyby odpowiednio ,,okrojonymi” macierzami Yk i Zk. Stąd, jeśli wszystkie wartości własne macierzy A mają różne moduły, to dla każdego takiego p¯ podprzestrzeń Ykp¯ rozpięta na początkowych p¯ wektorach macierzy Zk ,,zbiega” do podprzestrzeni własnej Pp¯, przy czym

∠⁢Ykp¯,Pp¯≤ρp¯k⋅∠⁢Y0p¯,Pp¯.

Z powyższej obserwacji wynika, że gdybyśmy realizowali iteracje podprzestrzeni przyjmując p=n i startując z macierzy identycznościowej Z0:=I formatu n×n to mielibyśmy zbieżność Ykp do podprzestrzeni własnych Pp dla wszystkich p=1,2,…,n.

W przypadku gdy pracujemy na bazach ortogonalnych i startujemy z Z0=I, zbieżność tą można wyrazić również w następujący sposób. Oznaczmy

V=ξ→1,ξ→2,…,ξ→n,⁢Zk=z→1k,z→2k,…,z→nk.

Lemat 3.1

Niech Bk będzie macierzą przejścia od bazy ξ→1,…,ξ→n do z→1k,…,z→nk,

Zk=V⁢⁢Bk.

Dla 1≤p≤n-1, niech

Bk=B1,1kB1,2kB2,1kB2,2k,

gdzie B1,1k jest podmacierzą formatu p×p, a pozostałe podmacierze odpowiednio formatów p×n-p, n-p×p, n-p×n-p. Jeśli ρp=λp+1/λp<1 i γ0<π/2 to

B1,2k2,B2,1k2⁢≤⁢ρpk⋅γ0.

Dowód

Niech Zk=Z1k,Z2k, V=V1,V2, gdzie Z1k i V1 są formatu n×p. Wtedy

Z1k=V1⁢⁢B1,2k+V2⁢⁢B2,1k.

Mnożąc to równanie z lewej strony przez V2T dostajemy B2,1k=V2T⁢⁢Z1k. Stąd

B2,1k2⁢=⁢V2T⁢⁢Z1k2⁢=⁢supx→2=1⁡V2T⁢⁢Z1k⁢⁢x→2.

Oznaczając y→=Z1k⁢⁢x→ zauważamy, że y→∈Yk:=span⁢z→1k,…,z→pk i y→2=1. Stąd wektor V2T⁢⁢y→ zawiera współczynniki rzutu ortogonalnego y→ na Pp⟂:=span⁢ξ→p+1,…,ξ→n w bazie ξ→p+1,…,ξ→n. Z definicji normy drugiej i twierdzenia 3.1 dostajemy więc

B2,1k2=supx→2=1⁡sin⁡∠⁢Z1k⁢⁢x→,Pp⟂≤ρpk⋅γ0.

Rozumując analogicznie i wykorzystując fakt, że

sin⁡∠⁢Yk,Pp=sin⁡∠⁢Yk⟂,Pp⟂

dostajemy taką samą nierówność dla B1,2k2.

∎

Z lematu 3.1 wynika natychmiast, że jeśli moduły wartości własnych macierzy A są parami różne to wyrazy pozadiagonalne bi,jk macierzy Bk zbiegają do zera, a ponieważ Bk jest ortogonalna, to bi,ik zbiegają do jedynki. Algorytm (IO) można więc uzupełnić o obliczanie w każdym kroku macierzy

Ak⁢:=⁢ZkT⁢⁢A⁢⁢Zk,

która powinna zbiegać do Λ=diag⁢λ1,…,λn. Fakt ten, łącznie z analizą szybkości zbieżności, pokażemy formalnie w następnym podrozdziale.

3.2. Metoda QR

3.2.1. Wyprowadzenie metody

Załóżmy, że realizujemy algorytm iteracji ortogonalnych (IO) z macierzą początkową Z0=I. Wtedy A⁢⁢Zk-1=Zk⁢⁢Rk, skąd

Ak-1=Zk-1T⁢⁢A⁢⁢Zk-1=Zk-1T⁢⁢Zk⁢⁢Rk.

Oznaczając Qk:=Zk-1T⁢⁢Zk widzimy, że ostatnia równość to nic innego jak rozkład ortogonalno-trójkątny macierzy Ak-1,

Ak-1=Qk⁢⁢Rk.

Prawdziwe są więc związki rekurencyjne

Ak	=	ZkT⁢⁢A⁢⁢Zk=ZkT⁢⁢Zk-1⁢⁢Zk-1T⁢⁢A⁢⁢Zk-1⁢⁢Zk-1T⁢⁢Zk
	=	QkT⁢⁢Ak-1⁢⁢Qk,
Zk	=	Zk-1⁢⁢Zk-1T⁢⁢Zk=Zk-1⁢⁢Qk.

Zależności te prowadzą do ALGORYTMU QR, który oblicza kolejne macierze Ak i Zk w następujący sposób.

			A0:=A;⁢Z0:=I;
			dla⁢⁢k=1,2,3,…
	(QR)		⁡Ak-1:=Qk⁢⁢Rk;⁢(ortonormalizacja),Ak:=Rk⁢⁢Qk;⁢Zk:=Zk-1⁢⁢Qk;⁢

3.2.2. QR a iteracje ortogonalne

Przypomnijmy, że w rozkładzie ortogonalno-trójkątnym macierz ortogonalna jest wyznaczona jednoznacznie z dokładnością do znaków jej wektorów-kolumn. Dlatego macierze Zk (a tym samym także Ak) powstałe w wyniku realizacji algorytmów (IO) i (QR) nie muszą być takie same. Algorytmy te są jednak równoważne w następującym sensie.

Lemat 3.2

Niech Z~k, A~k i Zk, Ak będą macierzami powstałymi w wyniku realizacji odpowiednio algorytmów (IO) i (QR). Wtedy

Z~k⁢=⁢Zk⁢⁢Dk⁢⁢oraz⁢⁢A~k=Dk⁢⁢Ak⁢⁢Dk,

gdzie Dk są pewnymi macierzami diagonalnymi z elementami ±1 na głównej przekątnej.

Dowód

Zastosujemy indukcję względem k. Dla k=1 mamy Z~0=I=Z0 oraz A~0=A=A0. Załóżmy, że lemat jest prawdziwy dla k-1, tzn. Z~k-1=Zk⁢⁢Dk-1. Wtedy, z jednej strony

A~k-1⁢=⁢Dk-1⁢⁢Ak-1⁢⁢Dk-1⁢=⁢Dk-1⁢⁢Qk⁢⁢Rk⁢⁢Dk-1,

a z drugiej strony

A~k-1⁢=⁢Z~k-1T⁢⁢Z~k⁢⁢R~k⁢⁢=⁢Dk-1⁢⁢Zk-1T⁢⁢Z~k⁢⁢R~k.

Mamy więc

Dk-1⁢⁢A~k-1⁢=⁢Zk-1T⁢⁢Z~k⁢⁢R~k⁢=⁢Qk⁢⁢Rk⁢⁢Dk-1,

przy czym równości te przedstawiają dwa rozkłady ortogonalno-trójkątne tej samej macierzy Dk-1⁢⁢A~k-1. A jeśli tak, to czynniki ortogonalne tych rozkładów różnią się jedynie znakami kolumn. Równoważnie, istnieje macierz diagonalna Dk z elementami ±1 na głównej przekątnej taka, że

Zk-1T⁢⁢Z~k⁢=⁢Q⁢⁢Dk⁢oraz⁢⁢R~k=Dk⁢⁢Rk⁢⁢Dk-1.

Stąd

Z~k⁢=⁢Zk-1⁢⁢Qk⁢⁢Dk⁢=⁢Zk⁢⁢Dk.

Dowód kończy następujący rachunek:

A~k⁢=⁢Z~kT⁢⁢A⁢⁢Z~k⁢=⁢DkT⁢⁢ZkT⁢⁢A⁢⁢Zk⁢⁢Dk⁢=⁢Dk⁢⁢Ak⁢⁢Dk.

∎

Zauważmy jeszcze, że jeśli przez a~i,jk i ai,jk oznaczymy odpowiednio elementy macierzy A~k i Ak oraz Dk=diag⁢d1k,…,dnk z dik=±1 to

a~i,jk=dik⁢djk⁢ai,jk.

To oznacza, że dla i≠j elementy te różnią się jedynie znakiem, a dla i=j są sobie równe.

3.2.3. Analiza zbieżności

Jesteśmy już gotowi do przedstawienia twierdzeń o zbieżności metody QR.

Twierdzenie 3.2

Dla 1≤p≤n, niech ρp:=λp+1/λp<1 i γk:=∠⁢Yk,Pp. Niech dalej

Ak=A1,1kA1,2kA2,1kA2,2k,⁢Qk=Q1,1kQ1,2kQ2,1kQ2,2k,⁢Rk=R1,1kR1,2k0R2,2k,

gdzie podmacierze A1,1k,Q1,1k,R1,1k są formatu p×p, a pozostałe podmacierze formatów odpowiednio p×n-p, n-p×p i n-p×n-p. Wtedy, przyjmując

εpk=ρpk⋅γ0

mamy γk≤εpk\tg oraz

A1,2k2=A2,1k2	≤	2⁢⁢εpk⁢A2,	(3.2)
Q1,2k2,⁢Q2,1k2	≤	2⁢⁢εpk,	(3.3)
R1,2k2	≤	4⁢⁢εpk⁢A2.	(3.4)

Dowód

Wobec wykazanej w lemacie \tg (teoretycznej) równoważności metody QR i iteracji (IO), nierówność γk≤εpk została już udowodniona w twierdzeniu 3.1. Aby pokazać (3.3), zapiszemy Zk w postaci Zk=Z1k,Z2k, gdzie Z1k i Z2k składają się odpowiednio z pierwszych p i z ostatnich n-p kolumn macierzy Zk. Wtedy

Q2,1=Z2k-1⁢⁢Z1k=VT⁢⁢Z2k-1T⁢⁢VT⁢⁢Z1k,

gdzie V=ξ→1,ξ→2,…,ξ→n. Pisząc V=V1,V2, podobnie jak dla Zk, mamy teraz

VTZ1k=[V1T⁢⁢Z1kV2T⁢⁢Z1k]=:[FG].

Ponieważ V12=Z1k2=1 to F2≤1, natomiast G2≤εpk, jak pokazaliśmy w dowodzie lematu 3.1.

Pisząc z kolei

VTZ2k-1=:[F′G′],

w analogiczny sposób pokazujemy, że F′2≤1 i G′2≤εpk. Stąd mamy

Q2,1k2=GT⁢⁢F+FT⁢⁢G2≤G′2⁢F+F′2⁢G2≤εk-1+εpk≤2⋅εpk.

Dokładnie tak samo pokazujemy oszacowanie dla Q1,22.

Aby pokazać (3.2) zauważmy, że wobec Ak=Qk⁢⁢Rk mamy też A2,1k=Q2,1k⁢⁢R1,1k. Ale

R1,12≤Rk2=A22=A2,

co w połączeniu z (3.3) dowodzi (3.2).

W końcu, wobec Rk=Ak⁢⁢QkT mamy

R1,2k=A1,1k⁢⁢Q2,1kT+A1,2k⁢⁢Q2,2kT.

Nierówność (3.4) wynika teraz z nierówności A1,1k2≤Ak2=A2, (3.2), (3.3), oraz Q2,2k2≤1.

∎

Z twierdzenia 3.2 wynika, że elementy pozadiagonalne macierzy Ak zbiegają do zera liniowo, przy czym iloraz zbieżności zależy od wartości ρp=λp+1/λp. Dokładniej, dla danego elementu ai,jk macierzy Ak, gdzie 1≤i<j≤n, prawdziwe jest oszacowanie

ai,jk≤min⁡εik,εi+1k,…⁢εj-1k

(gdzie skorzystaliśmy także z faktu, że moduł pojedynczego elementu macierzy jest nie większy od normy drugiej tej macierzy).

A jaka jest zbieżność elementów diagonalnych Ak do wartości własnych macierzy A? Oznaczmy

ε¯pk=⁡ε1kp=1,max⁡εp-1k,εpk2≤p≤n-1,εn-1kp=n.⁢

Twierdzenie 3.3

Dla wszystkich 1≤p≤n mamy

ap,pk-λp⁢≤ 4⋅ε¯pk2⋅A2.

Dowód

Wobec Zk=V⁢⁢Bk mamy

ap,pk	=	z→pkT⁢⁢A⁢⁢z→pk⁢=⁢V⁢⁢b→pkT⁢⁢A⁢⁢V⁢⁢b→pk
	=	b→pk⁢⁢VT⁢⁢A⁢⁢V⁢⁢b→pk⁢=⁢b→pk⁢⁢Λ⁢⁢b→pk
	=	λp⁢bp,pk2+∑i≠pλi⁢bi,jk2.

Stąd i z lematu 3.1 otrzymujemy

	ap,pk-λp	=	-λp⁢1-bp,pk2+∑i≠pλi⁢bi,pk2⁢=⁢∑i≠pbi,pk2⁢λi-λp
		≤	2⋅A2⋅∑i<pbi,pk2+∑p<ibi,pk2⁢≤ 4⋅ε¯pk2⋅A2.

∎

Elementy diagonalne ap,pk macierzy Ak zbiegają więc do wartości własnych λp macierzy A co najmniej jak λ2/λ12⁢k dla p=1, λn/λn-12⁢k dla p=n, oraz min⁡λp/λp-12⁢k,λp+1/λp2⁢k dla 2≤p≤n-1.

3.3. QR z przesunięciami i deflacja

Tak jak w przypadku metody potęgowej, można spróbować przyspieszyć zbieżność metody QR poprzez przesunięcie widma macierzy A, czyli zastosowanie QR do macierzy A¯=A-σ⁢⁢I. Oznaczmy przez A¯k kolejne macierze powstające w tym procesie. Z wcześniejszej analizy metody QR wynika, że jeśli

λi*-σ<mini≠i*⁡λi-σ

to wyraz a¯n,nk w prawym dolnym rogu macierzy A¯k zbiega do λi*-σ i zbieżność jest tym lepsza im σ jest bliższe λi*. Dlatego ważne jest, aby σ wybrać tak, aby dobrze aproksymowała jedną z wartości własnych macierzy A. Oczywiście, przesunięcie możnaby w każdym kroku wybierać w zależności od aktualnie dostępnej informacji. Metoda wyglądałaby wtedy następująco:

			A0:=A;⁢Z0:=I;
			dla⁢⁢k=1,2,3,…
	(QR-shift)		⁡Ak-1-σk⁢⁢I:=Qk⁢⁢Rk;⁢(ortonormalizacja)Ak:=Rk⁢⁢Qk+σk⁢⁢I;⁢Zk:=Zk-1⁢⁢Qk;⁢

Jasne, że

Ak=QkT⁢⁢Ak-1-σk⁢⁢I⁢⁢Qk+σk⁢⁢I⁢=⁢QkT⁢⁢Ak-1⁢⁢Qk,

tzn. kolejne macierze Ak są ortogonalnie podobne. Okazuje się, że QR ma jeszcze jedną ważną własność. Oznaczmy przez z→nk ostatnią kolumnę macierzy Zk, przez rn,nk wyraz w prawym dolnym rogu macierzy trójkątnej górnej Rk, oraz przez e→n ostatni wersor. Wtedy dla każdego k mamy

A-σk⁢⁢I⁢⁢z→nk	=	Zk-1⁢⁢Ak-1-σk-1⁢⁢IT⁢⁢Zk-1T⁢⁢z→nk
	=	Zk-1⁢⁢Qk⁢⁢RkT⁢⁢Zk-1T⁢⁢z→nk
	=	Zk-1⁢⁢RkT⁢⁢ZkT⁢⁢z→nk⁢=⁢Zk-1⁢⁢RkT⁢⁢e→n
	=	rn,nk⁢⁢Zk-1⁢⁢e→n⁢=⁢rn,nk⁢⁢z→nk-1,

czyli

znkrn,nk=A-σk⁢⁢I-1⁢⁢z→nk-1.

To zaś oznacza, że w kolejnych krokach iteracji QR z przesunięciami σk ostatnia kolumna macierzy Zk jest taka sama jak wektor x→k w odwrotnej metodzie potęgowej z wektorem startowym e→n i kolejnymi przesunięciami σk. A jeśli tak, to możemy wykorzystać wyniki rozdziału 2.4 i stwierdzić, że dobrym wyborem przesunięcia jest

σk=z→nkT⁢⁢A⁢⁢z→nk=e→nT⁢⁢ZkT⁢⁢A⁢⁢Zk⁢⁢e→n=e→n⁢⁢Ak⁢⁢e→n=an,nk.

Przy tak dobranym σk dostajemy asymptotycznie zbieżność sześcienną elementu an,nk do λi*.

Dla informacji podamy jeszcze, że możliwy jest też wybór przesunięcia Wilkinsona, gdzie jako σk bierze się tą wartość własną macierzy 2×2,

an-1,n-1kan-1,nkan,n-1kan,nk,

która jest najbliższa an,nk, a jeśli obie są jednakowo oddalone od an,nk to mniejszą z nich. Wtedy co prawda nie zwiększamy szybkości zbieżności, ale za to metoda jest zawsze zbieżna.

Pozostałe wartości własne macierzy A możemy obliczyć stosując deflację. Mianowicie, gdy już jesteśmy dostatecznie blisko wartości własnej λi* (co zwykle osiąga się wykonując kilka kroków QR i sprawdza przez obliczenie residuum A⁢⁢z→nk-an,nk⁢⁢z→nk2) to uznajemy, że wyrazy w ostatniej kolumnie i ostatnim wierszu macierzy Ak, poza wyrazem an,nk, są zerowe i redukujemy problem do macierzy formatu n-1×n-1. Oczywiście, ten ogólny przepis wymaga właściwej implementacji.

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

Matematyka obliczeniowa II

Zagadnienia

3. Zagadnienie własne III

3.1. Iteracje podprzestrzeni

3.1.1. Algorytm ogólny

Twierdzenie 3.1

Dowód

3.1.2. Iteracje ortogonalne

Lemat 3.1

Dowód

3.2. Metoda QR

3.2.1. Wyprowadzenie metody

3.2.2. QR a iteracje ortogonalne

Lemat 3.2

Dowód

3.2.3. Analiza zbieżności

Twierdzenie 3.2

Dowód

Twierdzenie 3.3

Dowód

3.3. QR z przesunięciami i deflacja