Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Matematyka obliczeniowa II – 4. Zagadnienie własne IV – MIM UW

Zagadnienia

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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.