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 $\varphi$ w płaszczyźnie rozpiętej na wersorach $\vec{e}_{i},\vec{e}_{j}$ , gdze $i<j$ , jest postaci

$O_{{i,j}}=\left[\begin{array}[]{ccccccccccc}1\\ &\ddots\\ &&\cos\varphi&&\sin\varphi\\ &&&\ddots\\ &&-\sin\varphi&&\cos\varphi\\ &&&&&\ddots\\ &&&&&&1\end{array}\right].$

Macierz ta różni się od jednostkowej jedynie wyrazami $o_{{i,i}}=c=o_{{j,j}}$ oraz $o_{{i,j}}=s=-o_{{j,i}}$ . łatwo sprawdzić, że $O_{{i,j}}$ jest przekształceniem ortogonalnym, $O_{{i,j}}\, O_{{i,j}}^{T}=I.$

W obliczeniach numerycznych macierze obrotu służą najczęściej do wyzerowania określonej współrzędnej danego wektora. Rzeczywiście, obracając wektor $\vec{a}=[a_{1},\ldots,a_{n}]^{T}$ w płaszczyźnie rozpiętej na $\vec{e}_{i}$ i $\vec{e}_{j}$ zmieniamy jedynie współrzędne $i$ -tą i $j$ -tą. Dokładniej, $O_{{i,j}}\,\vec{a}=\vec{b}$ , gdzie

$b_{k}=\left\{\begin{array}[]{rl}a_{i}\cos\varphi+a_{j}\sin\varphi,&k=i,\\ -a_{i}\sin\varphi+a_{j}\cos\varphi,&k=j,\\ a_{k},&k\ne i,j.\end{array}\right.$

Stąd, przyjmując

$\cos\varphi=\frac{a_{i}}{\sqrt{a_{i}^{2}+a_{j}^{2}}},\qquad\sin\varphi=\frac{a_{j}}{\sqrt{a_{i}^{2}+a_{j}^{2}}},$

albo liczby do nich przeciwne, dostajemy $b_{i}=\|\vec{a}\| _{2}^{2}$ albo $b_{i}=-\|\vec{a}\| _{2}^{2}$ , oraz $b_{j}=0$ . Z powyższych wzorów wynika, że kąt obrotu $\varphi$ można zawsze tak dobrać, aby $|\varphi|\le\pi/2$ .

Przekształcenia tak określone nazywamy obrotami Givensa.

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

$O_{{1,n}}\,\cdots\, O_{{1,3}}\, O_{{1,2}}\,\vec{a}=\pm\|\vec{a}\| _{2}^{2}\cdot\vec{e}_{1}.$

Takie operacje wymagają wykonania $6n$ mnożeń/dodawań i $n-1$ pierwiastkowań. Przypomnijmy, że ten sam efekt można uzyskać kosztem około $2n$ 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 $\vec{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 $O_{{i_{k},j_{k}}}$ dla $k=1,2,\ldots,m$ , gdzie $k$ -ty obrót jest wybrany tak, że $j_{k}$ -ta współrzędna wektora

$\vec{a}^{{(k)}}:=O_{{i_{k},j_{k}}}\,\cdots\, O_{{i_{1},j_{1}}}\,\vec{a}$

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

$O_{{i_{k},j_{k}}}\,=\, D_{k}\, S_{{i_{k},j_{k}}}\, D_{{k-1}}^{{-1}},$

gdzie $D_{0}=I$ , a $D_{k}$ są pewnymi nieosobliwymi macierzami diagonalnymi wybranymi tak, by mnożenie przez $S_{{i_{k},j_{k}}}$ było prostsze niż mnożenie przez $O_{{i_{k},j_{k}}}$ . Wtedy mamy

$\displaystyle O_{{i_{m},j_{m}}}\,\cdots\, O_{{i_{1},j_{1}}}$
	$\displaystyle=$	$\displaystyle(D_{m}\, S_{{i_{m},j_{m}}}\, D_{{m-1}}^{{-1}})\,\cdots\,(D_{2}\, S_{{i_{2},j_{2}}}\, D_{1}^{{-1}})\,(D_{1}\, S_{{i_{1},j_{1}}}\, D_{0}^{{-1}})$
	$\displaystyle=$	$\displaystyle D_{m}\, S_{{i_{m},j_{m}}}\,\cdots\, S_{{i_{2},j_{2}}}\, S_{{i_{1},j_{1}}}.$

Macierze diagonalne $D_{k}$ , a tym samym i przekształcenia $S_{{i_{k},j_{k}}}$ , zdefiniujemy indukcyjnie dla $k=1,2,\ldots,m$ . Załóżmy, że zdefiniowaliśmy już $D_{{k-1}}={\rm diag}(d_{1},\ldots,d_{n})$ . Niech $O_{{i_{k},j_{k}}}$ będzie obrotem o kąt $\varphi$ oraz $c=\cos\varphi$ , $s=\sin\varphi$ .

Dla uproszczenia zapisu podstawimy $D:=D_{{k-1}}$ , $\hat{D}:=D_{k}={\rm diag}(\hat{d}_{,}\ldots,\hat{d}_{n})$ oraz $p:=i_{k}$ , $q:=j_{k}$ . Dalej mamy $O_{{p,q}}=O_{{i_{k},j_{k}}}=(o_{{i,j}})_{{i,j=1}}^{n}$ , $S_{{p,q}}=S_{{i_{k},j_{k}}}=(s_{{i,j}})_{{i,j=1}}^{n}$ . Ponieważ $S_{{p,q}}=\hat{D}^{{-1}}\, O_{{p,q}}\, D$ to

$s_{{i,j}}\,=\,\frac{d_{j}}{\hat{d}_{i}}\, o_{{i,j}}\,=\,\left\{\begin{array}[]{rl}cd_{{p}}/\hat{d}_{{p}},&i=j=p,\\ cd_{{q}}/\hat{d}_{{q}},&i=j=q,\\ sd_{{q}}/\hat{d}_{{p}},&i=p,j=q,\\ -sd_{{p}}/\hat{d}_{{q}},&i=q,j=p,\\ d_{j}/\hat{d}_{i},&i=j\ne p,q,\\ 0,&\mbox{wpp}.\end{array}\right.$

Przypomnijmy, że obrót $O_{{i_{k},j_{k}}}$ jest wybrany tak, że $a^{{(k)}}_{{q}}=0$ . Zauważmy, że dla

$\vec{z}:=D^{{-1}}\,\vec{a}^{{(k-1)}}=S_{{i_{{k-1}},j_{{k-1}}}}\,\cdots\, S_{{i_{1},j_{1}}}\,\vec{a}$

mamy

$\vec{a}^{{(k)}}\,=\, O_{{p,q}}\,\vec{a}^{{(k-1)}}\,=\,(\hat{D}\, S_{{p,q}}\, D^{{-1}})\,(D\,\vec{z})\,=\,\hat{D}\, S_{{p,q}}\,\vec{z}.$

Stąd warunek $a^{{(k)}}_{{q}}=0$ jest równoważny warunkowi $\left(S_{{p,q}}\,\vec{z}\right)_{{q}}=0$ , czyli

$s\, z_{{p}}d_{{p}}=c\, z_{{q}}d_{{q}}.$

Jeśli teraz $z_{{q}}=0$ to przyjmujemy $\hat{D}:=D$ i $S_{{p,q}}:=I$ . W przeciwnym przypadku, uwzględniając równość $s^{2}+c^{2}=1$ , dostajemy

$c^{2}=1/t,\quad t:=1+\frac{d_{{q}}^{2}z_{{q}}^{2}}{d_{{p}}^{2}z_{{p}}^{2}},\qquad\qquad s^{2}=1/t^{{\prime}},\quad t^{{\prime}}:=1+\frac{d_{{p}}^{2}z_{{p}}^{2}}{d_{{q}}^{2}z_{{q}}^{2}}.$

Mamy teraz dwa przypadki.

(i) Jeśli $0<z_{{q}}^{2}d_{{q}}^{2}\le z_{{p}}^{2}d_{{p}}^{2}$ to kładziemy $\hat{d}_{{p}}=cd_{{p}}$ , $\hat{d}_{{q}}=cd_{{q}}$ i $\hat{d}_{l}=d_{l}$ dla $l\ne p,q$ . Wtedy $s_{{l,l}}=1$ dla $l\ne p,q$ oraz

$\left[\begin{array}[]{cc}s_{{p,p}}&s_{{p,q}}\\ s_{{q,p}}&s_{{q,q}}\end{array}\right]=\left[\begin{array}[]{cc}1&\alpha\\ -\beta&1\end{array}\right],$

gdzie

$\alpha=\frac{z_{{q}}d_{{q}}^{2}}{z_{{p}}d_{{p}}^{2}},\qquad\beta=\frac{z_{{q}}}{z_{{p}}}.$

(ii) Jeśli $0\le z_{{p}}^{2}d_{{p}}^{2}<z_{{q}}^{2}d_{{q}}^{2}$ to kładziemy $\hat{d}_{{p}}=sd_{{q}}$ , $\hat{d}_{{q}}=sd_{{p}}$ i $\hat{d}_{l}=d_{l}$ dla $l\ne p,q$ . Wtedy $s_{{l,l}}=1$ dla $l\ne p,q$ oraz

$\left[\begin{array}[]{cc}s_{{p,p}}&s_{{p,q}}\\ s_{{q,p}}&s_{{q,q}}\end{array}\right]=\left[\begin{array}[]{cc}\alpha^{{\prime}}&1\\ -1&\beta^{{\prime}}\end{array}\right],$

gdzie

$\alpha^{{\prime}}=\frac{z_{{p}}d_{{p}}^{2}}{z_{{q}}d_{{q}}^{2}},\qquad\beta^{{\prime}}=\frac{z_{{p}}}{z_{{q}}}.$

Macierz $S_{{p,q}}$ ma podobną strukturę jak $O_{{p,q}}$ . Obecność dwóch jedynek pozwalaja jednak zaoszczędzić dwa mnożenia przy operacji $S_{{p,q}}\,\vec{x}$ w porównaniu do $O_{{p,q}}\,\vec{x}$ . Ponadto, co jest ważniejsze, obliczanie wielkości $\alpha$ i $\beta$ (ew. $\alpha^{{\prime}}$ i $\beta^{{\prime}}$ ) definiujących $S_{{p,q}}$ nie wymaga pierwiastkowania. Rzeczywiście, w kolejnych krokach procesu obliczeniowego wystarczy pamiętać współczynniki $\hat{d}_{l}^{2}$ . Współczynniki te wyliczane są rekurencyjnie zgodnie ze wzorami $\hat{d}_{l}^{2}=d_{l}^{2}$ dla $l\ne p,q$ , a dla $l=p,q$

$\hat{d}_{l}^{2}=\left\{\begin{array}[]{rl}d_{l}^{2}/t&\mbox{jeśli (i)},\\ d_{l}^{2}/t^{{\prime}}&\mbox{jeśli (ii)}.\end{array}\right.$

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^{{\prime}}$ 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ą $D_{m}^{2}$ i aby pomnożyć przez $D_{m}$ należy najpierw wykonać $n$ pierwiastkowań. Nie jest to jednak zawsze prawdą.

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

$S_{{1,n}}\,\cdots\, S_{{1,3}}\, S_{{1,2}}\,\vec{a}$

jedynie pierwsza współrzędna otrzymanego wektora nie jest zerowa. Mnożenie przez $D_{m}$ 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ń

$S_{{n-1,n}}\,(S_{{n-2,n}}\, S_{{n-2,n-1}})\,\cdots\,(S_{{2,n}}\,\cdots\, S_{{2,3}})\,(S_{{1,n}}\,\cdots\, S_{{1,2}})\, A,$

gdzie $S_{{i,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 $D_{m}^{2}$ . 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\in\mathbb{R}^{{n,n}}$ jest macierzą symetryczną. Podobnie jak w metodzie QR, metoda Jacobiego startuje z macierzy $A_{0}=A$ i produkuje ciąg macierzy podobnych $A_{k}$ dla $k=1,2,3,\ldots$ , który zbiega do macierzy diagonalnej

$\Lambda={\rm diag}(\lambda _{1},\lambda _{2},\ldots,\lambda _{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,

	$\displaystyle A_{k}$	$\displaystyle=$	$\displaystyle O_{{i_{k},j_{k}}}^{T}\, A_{{k-1}}\, O_{{i_{k},j_{k}}}\,=\, O_{{i_{k},j_{k}}}^{T}\, O_{{i_{{k-1}},j_{{k-1}}}}^{T}\, A_{{k-2}}\, O_{{i_{k},j_{k}}}\, O_{{i_{{k-1}},j_{{k-1}}}}$
		$\displaystyle=$	$\displaystyle\ldots\,=\, O_{{i_{k},j_{k}}}^{T}\,\cdots\, O_{{i_{1},j_{1}}}^{T}\, A_{0}\, O_{{i_{1},j_{1}}}\,\cdots\, O_{{i_{k},j_{k}}}.$

Oznaczając $O_{k}=O_{{i_{1},j_{1}}}\,\cdots\, O_{{i_{k},j_{k}}}$ mamy więc $A_{k}=O_{k}^{T}\, A\, O_{k}$ .

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

$\omega _{k}:=\omega(A_{k})=\sum _{{i<j}}(a_{{i,j}}^{{(k)}})^{2},$

gdzie $(a_{{i,j}}^{{(k)}})_{{i,j=1}}^{n}$ . Wielkość $\omega _{k}$ dobrze charakteryzuje również błąd przybliżenia wartości własnych macierzy $A$ przez elementy diagonalne macierzy $A_{k}$ . Rzeczywiście, niech

$D_{k}={\rm diag}(a_{{1,1}}^{{(k)}},\ldots,a_{{n,n}}^{{(k)}})$

i $R_{k}=A_{k}-D_{k}$ . Wtedy

$D_{k}=O_{k}^{T}\, A\, O_{k}-O_{k}^{T}\, O_{k}\, R_{k}\, O_{k}^{T}\, O_{k}=O_{k}^{T}\,(A+E_{k})\, O_{k},$

gdzie $E_{k}=-O_{k}\, R_{k}\, O_{k}^{T}$ oraz

$\| E_{k}\| _{F}=\| R_{k}\| _{F}=\sqrt{2\omega _{k}}.$

Wyrazy diagonalne macierzy $D_{k}$ są więc wartościami własnymi macierzy $A+E_{k}$ , którą możemy potraktować jako zaburzoną macierz $A$ . Korzystając z twierdzenia Weyla 1.3 dostajemy, że istnieje takie uporządkowanie $\lambda _{{p(1)}},\ldots,\lambda _{{p(n)}}$ wartości własnych macierzy $A$ , że dla wszystkich $i$ mamy

$|a_{{i,i}}^{{(k)}}-\lambda _{{p(i)}}|\le\| E_{k}\| _{2}\le\| E_{k}\| _{F}=\sqrt{2\omega _{k}}.$

Stąd wniosek, że poszczególne obroty $O_{{i_{k},j_{k}}}$ powinniśmy wybierać tak, aby maksymalizować różnicę $\omega _{{k-1}}-\omega _{k}$ . Aby odpowiedzieć na pytanie który obrót osiąga ten cel, przeanalizujmy pojedynczy krok. Dla prostoty zapisu oznaczymy $p:=i_{k}$ , $q:=j_{k}$ , $A:=A_{{k-1}}$ , $\hat{A}:=A_{k}=O_{{p,q}}^{T}\, A\, O_{{p,q}}$ , gdzie $O_{{p,q}}$ jest obrotem o kąt $\varphi$ , $c=\cos\varphi$ , $s=\sin\varphi$ , $c^{2}+s^{2}=1$ .

Zauważmy, że mnożenie z lewej strony przez obrót $O_{{p,q}}^{T}$ zmienia jedynie wiersze $p$ -ty i $q$ -ty macierzy, a mnożenie z prawej strony przez $O_{{p,q}}$ zmienia kolumny $p$ -tą i $q$ -tą. Dlatego dla $i\ne p,q$ i $j\ne p,q$ wyrazy $\hat{a}_{{i,j}}$ i $a_{{i,j}}$ są sobie równe. Jeśli zaś $i\ne p,q$ i $j=p,q$ to

$\hat{a}_{{i,p}}=\hat{a}_{{p,i}}=ca_{{i,p}}-sa_{{i,q}},\qquad\quad\hat{a}_{{i,q}}=\hat{a}_{{q,i}}=sa_{{i,p}}+ca_{{i,q}},$

a stąd

$(\hat{a}_{{i,p}})^{2}+(\hat{a}_{{i,q}})^{2}=(a_{{i,p}})^{2}+(a_{{i,q}})^{2}$

i podobnie $(\hat{a}_{{p,i}})^{2}+(\hat{a}_{{q,i}})^{2}=(a_{{p,i}})^{2}+(a_{{q,i}})^{2}$ . Otrzymujemy

$\omega _{{k-1}}-\omega _{k}=\omega(A)-\omega(\hat{A})=(a_{{p,q}})^{2}-(\hat{a}_{{p,q}})^{2}.$

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

$\hat{a}_{{p,q}}=0.$

Jeśli $a_{{p,q}}=0$ to powyższą równość realizuje obrót zerowy. Załóżmy więc, że $a_{{p,q}}\ne 0$ . Bezpośredni rachunek pokazuje, że wtedy

$\left[\begin{array}[]{cc}\hat{a}_{{p,p}}&\hat{a}_{{p,q}}\\ \hat{a}_{{q,p}}&\hat{a}_{{q,q}}\end{array}\right]=\left[\begin{array}[]{ll}c^{2}a_{{p,p}}-2sca_{{p,q}}+s^{2}a_{{q,q}}&(c^{2}-s^{2})a_{{p,q}}+sc(a_{{p,p}}-a_{{q,q}})\\ (c^{2}-s^{2})a_{{p,q}}+sc(a_{{p,p}}-a_{{q,q}})&s^{2}a_{{p,p}}+2sca_{{p,q}}+c^{2}a_{{q,q}}\end{array}\right].$

Stąd równość $\hat{a}_{{p,q}}=0$ zachodzi wtedy gdy

$(c^{2}-s^{2})a_{{p,q}}+sc(a_{{p,p}}-a_{{q,q}})=0.$

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

$t^{2}+2\,\beta\, t-1\,=\, 0,\qquad\quad\beta=\frac{a_{{q,q}}-a_{{p,p}}}{2a_{{p,q}}}.$

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

$c=\frac{1}{\sqrt{1+t^{2}}},\qquad s=t\cdot c,\qquad\mbox{gdzie}\quad t=\frac{{\rm sign}\beta}{|\beta|+\sqrt{1+\beta^{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

$|a_{{p,q}}|=\max _{{i<j}}|a_{{i,j}}|.$

Ponieważ wtedy $a_{{p,q}}^{2}\ge\omega(A)/N$ , gdzie $N=n(n-1)/2$ to

$\omega(\hat{A})=\omega(A)-a_{{p,q}}^{2}\le\omega(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),\ldots,(1,i),\ldots,(i-1,i),\ldots,(1,n),\ldots,(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 $a_{{p,q}}^{2}$ jest znacznie mniejsze od $\omega(A)$ , co nie jest opłacalne. Aby temu zapobiec, można wprowadzić pewien próg $\delta$ i nie wykonywać obrotów gdy $a_{{p,q}}^{2}<\delta$ . Na przykład, możemy wziąć $\delta=\omega(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=A^{T}\in\mathbb{R}^{{n,n}}$ jest dodatkowo trójdiagonalna,

$A=T=\left[\begin{array}[]{cccc}a_{1}&b_{1}&&\\ b_{1}&\ddots&\ddots&\\ &\ddots&\ddots&b_{{n-1}}\\ &&b_{{n-1}}&a_{n}\end{array}\right].$

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 $O_{{2,i}}$ dla $i=3,4,\ldots,n$ tak, aby mnożenie

$O_{{2,n}}^{T}\,\cdots\, O_{{2,4}}^{T}\, O_{{2,3}}^{T}\, A=:\hat{A}$

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

$\hat{A}\, O_{{2,3}}\, O_{{2,4}}\,\cdots\, O_{{2,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=O^{T}\, A\, O$ , gdzie

$O\,=\, O_{{2,3}}\,\cdots\, O_{{2,n}}\, O_{{3,4}}\,\cdots\, O_{{3,n}}\,\cdots\, O_{{n-2,n-1}}\, O_{{n-2,n}}\, O_{{n-1,n}},$

przy czym mnożenie z lewej przez $O_{{i,j}}^{T}$ 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 $O_{{i,i+1}}^{T}$ zerujących kolejno elementy $(i+1,i)$ dla $i=1,2,\ldots,n-1$ . Wtedy

$R\,=\,(O_{{n-1,n}}^{T}\, O_{{n-2,n-1}}^{T}\,\cdots\, O_{{2,3}}^{T}\, O_{{1,2}}^{T})\, T$

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

$(O_{{1,2}}\, O_{{2,3}}\,\cdots\, O_{{n-1,n}})\, R$

otrzymana macierz $\hat{T}$ pozostaje trójdiagonalna. Rzeczywiście, skoro $R$ jest trójkątna górna, a mnożenie z lewej strony przez $O_{{i,i+1}}$ kombinuje jedynie wiersze $i$ -ty i $(i+1)$ -szy to $\hat{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 $\pm 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 $A^{T}$ .

Twierdzenie 4.1

Dla dowolnej macierzy $A\in\mathbb{R}^{{m,n}}$ , gdzie $m\ge n$ , istnieją macierze ortogonalne $U\in\mathbb{R}^{{m,m}}$ i $V\in\mathbb{R}^{{n,n}}$ takie, że

$A\,=\, U\,\Sigma\, V^{T},$

gdzie $\Sigma\in\mathbb{R}^{{m,n}}$ jest macierzą diagonalną,

$\Sigma=\left[\begin{array}[]{c}\hat{\Sigma}\\ 0\end{array}\right],\qquad\hat{\Sigma}={\rm diag}(\sigma _{1},\sigma _{2},\ldots,\sigma _{n})\in\mathbb{R}^{{n,n}},$

$\sigma _{1}\ge\sigma _{2}\ge\cdots\ge\sigma _{k}>\sigma _{{k+1}}=\cdots=\sigma _{n}=0$ , a $k$ jest rzędem macierzy $A$ .

Dowód

Ponieważ macierz $A^{T}\, A$ jest symetryczna i nieujemnie określona to istnieje baza ortonormalna $\{\vec{v}_{i}\} _{{i=1}}^{n}$ 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 $\sigma _{i}^{2}$ ,

$\sigma _{1}\ge\cdots\ge\sigma _{k}>\sigma _{{k+1}}=\ldots=\sigma _{n}=0,$

przy czym $k={\rm rank}(A^{T}\, A)={\rm rank}(A)$ . Zauważmy, że

$(A\,\vec{v}_{i})^{T}\,(A\,\vec{v}_{j})=\vec{v}_{i}^{T}\,(A^{T}\, A)\,\vec{v}_{j}=\sigma^{2}_{j}\cdot(\vec{v}_{i}^{T}\,\vec{v}_{j})=\left\{\begin{array}[]{rl}0&i\ne j,\\ \sigma _{j}^{2}&i=j.\end{array}\right.$

To oznacza, że wektory $\vec{u}_{i}:=A\,\vec{v}_{i}/\sigma _{i}$ , $1\le i\le k$ , tworzą układ ortonormalny i można go uzupełnić wektorami $\vec{u}_{{k+1}},\ldots,\vec{u}_{m}$ do bazy ortonormalnej w $\mathbb{R}^{m}$ .

Zdefiniujmy teraz macierze ortogonalne

$V=[V_{1},V_{2}]\in\mathbb{R}^{{n,n}},\qquad V_{1}=[\vec{v}_{1},\ldots,\vec{v}_{k}]\in\mathbb{R}^{{n,k}},\quad V_{2}=[\vec{v}_{{k+1}},\ldots,\vec{v}_{n}]\in\mathbb{R}^{{n,n-k}},$

$U=[U_{1},U_{2}]\in\mathbb{R}^{{m,m}},\qquad U_{1}=[\vec{u}_{1},\ldots,\vec{u}_{k}]\in\mathbb{R}^{{m,k}},\quad U_{2}=[\vec{u}_{{k+1}},\ldots,\vec{u}_{n}]\in\mathbb{R}^{{m,m-k}}.$

Wtedy, oznaczając $\Sigma^{{\prime}}={\rm diag}(\sigma _{1},\ldots,\sigma _{k})\in\mathbb{R}^{{k,k}}$ , mamy

$U^{T}\, A\, V=\left[\begin{array}[]{c}U_{1}^{T}\\ U_{2}^{T}\end{array}\right]\,[A\, V_{1},A\, V_{2}]=\left[\begin{array}[]{c}U_{1}^{T}\\ U_{2}^{T}\end{array}\right]\,[U_{1}\,\Sigma^{{\prime}},0]=\left[\begin{array}[]{cc}\Sigma^{{\prime}}&0\\ 0&0\end{array}\right]=\Sigma,$

albo, równoważnie, $A=U\,\Sigma\, V^{T}$ .

∎

Wielkości $\sigma _{i}$ , $1\le i\le n$ , to właśnie wartości szczególne macierzy $A$ . Są one wyznaczone jednoznacznie. Rzeczywiście, jeśli mamy dwa rozkłady, $U_{1}\,\Sigma _{1}\, V_{1}^{T}=A=U_{2}\,\Sigma _{2}\, V_{2}^{T}$ , to

$V\,(\Sigma _{1}^{T}\,\Sigma _{1})=(\Sigma _{2}^{T}\,\Sigma _{2})\, V,\qquad V=V_{2}^{T}\, V_{1}.$

Porównując wyrazy diagonalne po obu stronach powyższej równości dostajemy $\sigma _{i}^{{(1)}}=\sigma _{i}^{{(2)}}$ 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 $\vec{x}^{*}$ minimalizujący $\|\vec{b}-A\,\vec{x}\| _{2}$ po wszystkich $\vec{x}\in\mathbb{R}^{n}$ . (Jeśli istnieje wiele takich wektorów to bierzemy ten o najmniejszej normie.) Mamy bowiem

$\|\vec{b}-A\,\vec{x}\| _{2}=\| U\, U^{T}\,\vec{b}-U\,\Sigma\, V^{T}\,\vec{x}\| _{2}=\|\vec{c}-\Sigma\,\vec{y}\| _{2},$

$\vec{c}=U^{T}\,\vec{b}$ , $\vec{y}=V^{T}\,\vec{x}$ . Stąd $\vec{x}^{*}=V^{T}\,\vec{y}^{*}$ , gdzie

$y_{i}^{*}=\left\{\begin{array}[]{rl}c_{i}/\sigma _{i},&1\le i\le k,\\ 0,&1\le i\le n.\end{array}\right.$

4.4.2. Dlaczego nie pomnożyć $A^{T}\, A$ ?

Z dowodu twierdzenia 4.1 wynika, że $\sigma _{1}^{2},\ldots,\sigma _{n}^{2}$ są wartościami własnymi macierzy symetrycznej i nieujemnie określonej $A^{T}\, A$ , a kolumny macierzy ortogonalnej $V$ tworzą odpowiednio bazę ortonormalną w $\mathbb{R}^{n}$ wektorów własnych tej macierzy. Podobnie, $\sigma _{1}^{2},\ldots,\sigma _{n}^{2},\underbrace{0,\ldots,0}_{{m-n}}$ są wartościami własnymi $A\, A^{T}$ , a kolumny $U$ tworzą bazę ortonormalną w $\mathbb{R}^{m}$ 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 $\hat{A}:=A^{T}\, A$ lub $\hat{A}:=A\, A^{T}$ , a następnie aplikując do macierzy $\hat{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=\left[\begin{array}[]{cc}1&1\\ 0&\varepsilon\end{array}\right].$ Dla ,,małych” $\varepsilon$ wartości szczególne tej macierzy są bliskie $\sqrt{2}$ i $|\varepsilon|/\sqrt{2}$ . Jeśli $\varepsilon$ jest na tyle małe, że $1+\varepsilon^{2}$ jest w arytmetyce fl reprezentowane przez $1$ to macierz $A^{T}\, A=\left[\begin{array}[]{cc}1&1\\ 1&1+\varepsilon^{2}\end{array}\right]$ jest reprezentowana przez macierz $\left[\begin{array}[]{cc}1&1\\ 1&1\end{array}\right]$ , 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ń $A^{T}\, A$ lub $A\, A^{T}$ . 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 $\hat{A}=A^{T}\, A$ konstruuje ciąg macierzy podobnych $\hat{A}=\hat{A}_{0},\hat{A}_{1},\ldots,\hat{A}_{k},\ldots$ , gdzie $\hat{A}_{k}=O^{T}_{{i_{k},j_{k}}}\,\hat{A}_{{k-1}}\, O_{{i_{k},j_{k}}}$ , a $O_{{i_{k},j_{k}}}$ jest tak dobranym obrotem Givensa, aby wyzerować element $(i_{k},j_{k})$ i jednocześnie, wobec symetrii, także $(j_{k},i_{k})$ . Oznaczmy $A_{0}=A$ i $A_{k}=A_{{k-1}}\, O_{{i_{k},j_{k}}}$ dla $k\ge 1$ . Wtedy $\hat{A}_{k}=A_{k}^{T}\, A_{k}$ . Pomysł polega na tym, aby zamiast obliczać jawnie $\hat{A}_{k}=O_{{i_{k},j_{k}}}^{T}\, A_{{k-1}}^{T}\, A_{{k-1}}\, O_{{i_{k},j_{k}}}$ , w kolejnych krokach obliczać i pamiętać jedynie $A_{k}=A_{{k-1}}\, O_{{i_{k},j_{k}}}$ . Jest to możliwe, bowiem wyznaczenie rotacji $O_{{i_{k},j_{k}}}$ wymaga jedynie znajomości elementów

$\hat{a}_{{i_{k},i_{k}}}^{{(k-1)}}=\sum _{{l=1}}^{m}(a_{{l,i_{k}}}^{{(k-1)}})^{2},\qquad\hat{a}_{{j_{k},j_{k}}}^{{(k-1)}}=\sum _{{l=1}}^{m}(a_{{l,j_{k}}}^{{(k-1)}})^{2},\qquad\hat{a}_{{i_{k},j_{k}}}^{{(k-1)}}=\sum _{{l=1}}^{m}a_{{l,i_{k}}}^{{(k-1)}}a_{{l,j_{k}}}^{{(k-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 $A^{T}\, A$ . Zanim jednak przejdziemy do samego algorytmu, poczynimy kilka pomocniczych uwag teoretycznych.

Niech $T_{0}$ będzie macierzą symetryczną i dodanio okteśloną, $T_{0}=T^{T}_{0}>0$ . Rozpatrzmy proces iteracyjny, nazwijmy go LR, który startuje z macierzy $T_{0}$ i w kolejnych krokach $k=1,2,\ldots$

(i) wybiera przesunięcie $\tau _{k}$ , mniejsze od najmniejszej wartości własnej $T_{{k-1}}$ (np. $\tau _{k}=0$ ),
(ii) dokonuje rozkładu Banachiewicza-Cholesky'ego

$T_{{k-1}}-\tau _{k}^{2}\cdot I=B_{k}^{T}\, B_{k},$

gdzie $B_{k}$ jest macierzą trójkątną górną z dodatnimi elementami na głównej przekątnej,
(iii) produkuje $T_{k}=B_{k}\, B_{k}^{T}+\tau _{k}^{2}\cdot I$ .

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

$T_{k}=B_{k}\, B_{k}^{T}+\tau _{k}^{2}\cdot I=B_{k}\,(T_{{k-1}}-\tau _{k}^{2}\cdot I)\, B_{k}^{{-1}}+\tau _{k}^{2}\cdot I=B_{k}\, T_{{k-1}}\, B_{k}^{{-1}}.$

Zachodzi ciekawsza własność.

Lemat 4.1

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

Dowód

Bez zmniejszenia ogólności przyjmijmy, że $k=0$ i $\tau _{1}=\tau _{2}=\tau$ . Z jednej strony, z rozkładu QR mamy $(T_{0}-\tau^{2}\cdot I)^{2}=(Q\, R)^{T}\,(Q\, R)=R^{T}\, R$ , gdzie elementy na głównej przekątnej macierzy $R$ są dodatnie, a z drugiej z rozkładu LR

	$\displaystyle(T_{0}-\tau^{2}\cdot I)^{2}$	$\displaystyle=$	$\displaystyle(B_{1}^{T}\, B_{1})\,(B_{1}^{T}\, B_{1})\,=\, B_{1}^{T}\,(T_{1}-\tau^{2}\cdot I)\, B_{1}$
		$\displaystyle=$	$\displaystyle B_{1}^{T}\, B_{2}^{T}\, B_{2}\, B_{1}\,=\,(B_{2}\, B_{1})^{T}\,(B_{2}\, B_{1}).$

Wobec jednoznaczności rozkładu Banachiewicza-Cholesky'ego mamy więc $R=B_{2}\, B_{1}$ . Stąd macierz powstała w wyniku jednego kroku iteracji QR wynosi

$\displaystyle\hat{T}$	$\displaystyle=$	$\displaystyle R\, Q+\tau^{2}\cdot I\,=\, R\,(Q\, R)\, R^{{-1}}+\tau^{2}\cdot I$
	$\displaystyle=$	$\displaystyle R\,(T_{0}-\tau^{2}\cdot I)\, R^{{-1}}+\tau^{2}\cdot I\,=\, R\, T_{0}\, R^{{-1}}$
	$\displaystyle=$	$\displaystyle(B_{2}\, B_{1})\,(B_{1}^{T}\, B_{1}+\tau^{2}\cdot I)\,(B_{2}\, B_{1})^{{-1}}$
	$\displaystyle=$	$\displaystyle B_{2}\, B_{1}\, B_{1}^{T}\, B_{1}\, B_{1}^{{-1}}\, B_{2}^{{-1}}+\tau^{2}(B_{2}\, B_{1})\,(B_{2}\, B_{1})^{{-1}}$
	$\displaystyle=$	$\displaystyle B_{2}\,(B_{1}\, B_{1}^{T})\, B_{2}^{{-1}}+\tau^{2}\cdot I\,=\, B_{2}\,(T_{1}-\tau^{2}\cdot I)\, B_{2}^{{-1}}+\tau^{2}\cdot I$
	$\displaystyle=$	$\displaystyle B_{2}\,(B_{2}^{T}\, B_{2})\, B_{2}^{{-1}}+\tau^{2}\cdot I\,=\, B_{2}\, B_{2}^{T}+\tau^{2}\cdot I$
	$\displaystyle=$	$\displaystyle(T_{2}-\tau^{2}\cdot I)+\tau^{2}\cdot I\,=\, T_{2},$

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\in\mathbb{R}^{{m,n}}$ jest kolumnowo regularna, czyli ${\rm rank}(A)=n$ , oraz dwudiagonalna. Dokładniej, $a_{{i,j}}=0$ dla $i\ge j+1$ i dla $j\ge i+2$ ,

$A=\left[\begin{array}[]{cccc}a_{1}&b_{1}&&\\ &a_{2}&\ddots&\\ &&\ddots&b_{{n-1}}\\ &&&a_{n}\\ &&&\\ &&&\end{array}\right].$

Oczywiście, wobec kolumnowej regularności macierzy mamy $a_{j}\ne 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 $T_{k}$ , podobne do $A^{T}\, A$ , pracuje bezpośrednio na macierzach $B_{k}$ . Przyjmujemy $B_{0}=D\, A$ , gdzie $D\in\mathbb{R}^{{n,m}}$ jest macierzą diagonalną z elementami na głównej przekątnej $d_{j}={\rm sign}(a_{j})$ , $1\le j\le n$ . Wtedy $B_{0}\in\mathbb{R}^{{n,n}}$ jest macierzą dwudiagonalną z dodatnimi elementami na głównej przekątnej oraz

$T_{0}=B_{0}^{T}\, B_{0}=A^{T}\,(D^{T}\, D)\, A=A^{T}\, A$

jest macierzą trójdiagonalną. Zobaczmy teraz jak mając macierz dwudiagonalną $B_{k}$ możemy skonstruować $B_{{k+1}}$ . Ponieważ $T_{k}$ jest trójdiagonalna to $B_{{k+1}}$ musi być dwudiagonalna. (W szczególności, wyniknie to również z dalszych rachunków.) Dla uproszczenia zapisu, oznaczmy elementy przekątniowe macierzy $B_{k}$ i $B_{{k+1}}$ odpowiednio przez $a_{j}$ i $\hat{a}_{j}$ , $1\le j\le n$ , a elementy nad główną przekątną przez $b_{j}$ i $\hat{b}_{j}$ , $1\le j\le n-1$ . Przyjmijmy dodatkowo $b_{0}=\hat{b}_{0}=b_{n}=\hat{b}_{n}=0$ . Macierze $B_{k}$ i $B_{{k+1}}$ związane są równaniem

$B_{{k+1}}^{T}\, B_{{k+1}}+\tau _{{k+1}}^{2}\cdot I\,=\, T_{k}\,=\, B_{k}\, B_{k}^{T}+\tau _{k}^{2}\cdot I.$

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

$\hat{a}_{j}^{2}+\hat{b}_{{j-1}}^{2}+\tau _{{k+1}}^{2}=a_{j}^{2}+b_{j}^{2}+\tau _{k}^{2},\qquad 1\le j\le n,$

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

$\hat{a}_{j}\hat{b}_{j}=a_{{j+1}}b_{j}^{2},\qquad 1\le j\le n-1.$

Stąd, podstawiając $\delta=\tau _{{k+1}}^{2}-\tau _{k}^{2}$ , mamy następujące wzory na $\hat{a}_{j}$ i $\hat{b}_{j}$ . Dla $j=1,2,\ldots,n-1$ obliczamy

$\hat{a}_{j}:=\sqrt{a_{j}^{2}+b_{j}^{2}-\hat{b}_{{j-1}}^{2}-\delta},\qquad\hat{b}_{j}:=a_{{j+1}}b_{j}/\hat{a}_{j}$

i na końcu $\hat{a}_{n}:=\sqrt{a_{n}^{2}-\hat{b}_{{n-1}}^{2}-\delta}$ .

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 $a_{j}$ i $b_{j}$ . Rzeczywiście, wprowadzając zmienne $p_{j}=a_{j}^{2}$ i $q_{j}=b_{j}^{2}$ , otrzymujemy następującą procedurę dla jednego kroku iteracyjnego:

.	begin
.	.	for $j=1$ to $n-1$ do
.	.	begin
.	.	.	$\hat{p}_{j}:=p_{j}+q_{j}-\hat{q}_{j}-\delta;$
.	.	.	$\hat{q}_{j}:=q_{j}\cdot(q_{{j+1}}/\hat{q}_{j})$
.	.	end;
.	.	$\hat{p}_{n}:=p_{n}-\hat{q}_{{n-1}}-\delta$
.	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 $H_{{L,1}}\in\mathbb{R}^{{m,m}}$ przekształcające pierwszą kolumnę $A$ na kierunek $\vec{e}_{1}$ . Oznaczmy $\tilde{A}_{1}=(\tilde{a}_{{i,j}})=H_{{L,1}}\, A$ . Następnie wybieramy odbicie $\tilde{H}_{{R,1}}\in\mathbb{R}^{{n-1,n-1}}$ tak, aby przekształcało wektor $[\tilde{a}_{{1,2}},\tilde{a}_{{1,3}}\ldots,\tilde{a}_{{1,n}}]^{T}\in\mathbb{R}^{{n-1}}$ na kierunek $\vec{e}_{1}\in\mathbb{R}^{{n-1}}$ i mnożymy $\tilde{A}_{1}$ z prawej przez $H^{T}_{{R,1}}=\left[\begin{array}[]{cc}1&\vec{0}^{T}\\ \vec{0}&\tilde{H}^{T}_{{R,1}}\end{array}\right]\in\mathbb{R}^{{n,n}}$ . Ponieważ to mnożenie nie zmienia pierwszej kolumny $\tilde{A}_{1}$ , w powstałej macierzy $A_{1}=H_{{L,1}}\, A\, H^{T}_{{R,1}}$ 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

$\hat{A}=U_{1}\, A\, V_{1}^{T},\quad\mbox{gdzie}\quad U_{1}=H_{{L,n-1}}\,\cdots\, H_{{L,1}},\quad V_{1}=H_{{R,n-2}}\,\cdots\, H_{{R,1}},$

która jest już dwudiagonalna. Jeśli teraz $\hat{A}=U_{2}\,\Sigma\, V_{2}^{T}$ to $A=(U_{1}^{T}\, U_{2})\,\Sigma\,(V_{1}^{T}\, V_{2})^{T}$ , a więc $A$ i $\hat{A}$ mają te same wartości szczególne. Ponadto,

$A=U\,\Sigma\, V^{T},\quad\mbox{gdzie}\quad U=U_{1}^{T}\, U_{2},\quad V=V_{1}^{T}\, V_{2}.$

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\cdot n^{2}$ .

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ć ?

Przykład 4.1

4.4.3. SVD dla macierzy dwudiagonalnych

Lemat 4.1

Dowód

4.4.2. Dlaczego nie pomnożyć $A^{T}\, A$ ?