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 – 3. Zagadnienie własne III – MIM UW

Zagadnienia

3. Zagadnienie własne III

Metodę potęgową można traktować jako iteracje podprzestrzeni jednowymiarowej startujące ze {\rm span}(\vec{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=A^{T}\in\mathbb{R}^{{n,n}} i {\rm det}(A)\ne 0. Przypomnijmy, że wtedy istnieje baza ortonormalna \{\vec{\xi}_{i}\} _{{i=1}}^{n} wektorów własnych macierzy,

A\,\vec{\xi}_{i}=\lambda _{i}\,\vec{\xi}_{i},\quad 1\le i\le n.

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

|\lambda _{1}|\ge|\lambda _{2}|\ge\cdots\ge|\lambda _{n}|>0.

3.1. Iteracje podprzestrzeni

3.1.1. Algorytm ogólny

Niech 1\le p\le n-1 oraz {\mathcal{Y}_{0}}\subseteq\mathbb{R}^{n} będzie podprzestrzenią o wymiarze p. Dla k=1,2,3,\ldots rozpatrzmy następujące iteracje:

{\mathcal{Y}}_{k}\,:=\, A({\mathcal{Y}}_{{k-1}})\,=\,\{\, A\,\vec{x}:\,\vec{x}\in{\mathcal{Y}}_{{k-1}}\,\}.

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

{\mathcal{P}}^{{(p)}}:={\rm span}(\vec{\xi}_{1},\vec{\xi}_{2},\ldots,\vec{\xi}_{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 {\mathcal{Y}},{\mathcal{Z}}\subseteq\mathbb{R}^{n} definiujemy jako

\angle({\mathcal{Y}}_{k},{\mathcal{P}}^{{(p)}}):=\max _{{\vec{y}\in{\mathcal{Y}}_{k}}}\min _{{\vec{z}\in{\mathcal{P}}^{{(p)}}}}\angle(\vec{y},\vec{z}).

(Jako ćwiczenie można wykazać, że jeśli {\rm dim}({\mathcal{Y}})={\rm dim}({\mathcal{Z}}) to \angle({\mathcal{Y}},{\mathcal{Z}})=\angle({\mathcal{Z}},{\mathcal{Y}}).)

Twierdzenie 3.1

Niech

\rho _{p}:=\frac{|\lambda _{{p+1}}|}{|\lambda _{p}|}<1\quad\mbox{ oraz }\quad\gamma _{k}:=\angle({\mathcal{Y}}_{k},{\mathcal{P}}^{{(p)}}),\quad k=0,1,2,\ldots.

Jeśli \gamma _{0}<\pi/2 to

\tg\gamma _{k}\,\le\,\rho _{p}\cdot\tg\gamma _{{k-1}}\,\le\,\rho _{p}^{k}\cdot\tg\gamma _{0}.
Dowód

Niech \vec{y} będzie wektorem spełniającym \angle(\vec{y},{\mathcal{P}}^{{(p)}})<\pi/2. Wtedy

\vec{y}=\sum _{{i=1}}^{n}\vec{v}_{i},\qquad\vec{v}_{i}\in{\rm span}(\vec{\xi}_{i}),\quad 1\le i\le n,

gdzie nie wszystkie \vec{v}_{i} dla 1\le i\le p, są są zerowe, oraz A\,\vec{y}=\sum _{{i=1}}^{n}\lambda _{i}\vec{v}_{i}. Stąd

\tg\angle(A\,\vec{y},{\mathcal{P}}^{{(p)}})=\sqrt{\frac{\sum _{{i=p+1}}^{n}|\lambda _{i}|^{2}\|\vec{v}_{i}\| _{2}^{2}}{\sum _{{i=1}}^{p}|\lambda _{i}|^{2}\|\vec{v}_{i}\| _{2}^{2}}}\le\frac{|\lambda _{{p+1}}|}{|\lambda _{p}|}\cdot\sqrt{\frac{\sum _{{i=p+1}}^{n}\|\vec{v}_{i}\| _{2}^{2}}{\sum _{{i=1}}^{p}\|\vec{v}_{i}\| _{2}^{2}}}=\rho _{p}\cdot\tg\angle(\vec{y},{\mathcal{P}}^{{(p)}}p). (3.1)

Jeśli teraz \gamma _{{k+1}}=\angle(\vec{y}_{{k+1}},{\mathcal{P}}^{{(p)}}) dla \vec{y}_{{k+1}}\in{\mathcal{Y}}_{{k+1}} oraz \vec{y}_{{k+1}}=A\,\vec{y}_{k} to z (3.1) wynika, że

\tg\gamma _{{k+1}}\le\rho _{p}\cdot\tg\angle(\vec{y}_{k},{\mathcal{P}}^{{(p)}})\le\rho _{p}\cdot\tg\gamma _{k},

co kończy dowód.

3.1.2. Iteracje ortogonalne

W praktyce, iteracje podprzestrzeni realizujemy reprezentując kolejne {\mathcal{Y}}_{k} przez ich bazy ortonormalne. Wtedy algorytm, zwany w tym przypadku również iteracjami ortogonalnymi, przybiera następującą postać. Wybieramy macierz kolumnowo-ortogonalną Z_{0} formatu n\times p jako macierz startową, a następnie dla k=1,2,3,\ldots iterujemy:

\mbox{(IO)}\quad\left\{\begin{array}[]{l}Y_{k}\,:=\, A\, Z_{{k-1}};\\
Y_{k}\,:=\, Z_{k}\, R_{k};\quad\mbox{(ortonormalizacja)}\end{array}\right.

przy czym w drugiej linii następuje ortonormalizacja macierzy Y_{k} realizowana poprzez rozkład tej macierzy na iloczyn macierzy kolumnowo-ortonormalnej Z_{k} formatu n\times p i macierzy trójkątnej górnej R_{k} formatu p\times p. Dokonuje się tego znanym nam już algorytmem wykorzystującym odbicia Householdera.

łatwo zauważyć, że jeśli przyjmiemy, iż wektory-kolumny macierzy Z_{0} tworzą bazę ortonormalną podprzestrzeni {\mathcal{Y}}_{0} to w kolejnych krokach wektory-kolumny Z_{k} tworzą bazę ortonormalną podprzestrzeni {\mathcal{Y}}_{k}. Rzeczywiście, ponieważ Y_{k}=A\, Z_{{k-1}} to kolumny Y_{k} stanowią bazę {\mathcal{Y}}_{k}. Kolumny te są z kolei liniową kombinacją kolumn Z_{k}, czyli rozpinają tą samą podprzestrzeń co kolumny Z_{k}.

Poczynimy teraz kluczową obserwację. Przypuśćmy, że wszystkie wartości własne \lambda _{i} są parami różne. Gdybyśmy rozpoczynali iteracje od macierzy \overline{Y}_{0} składającej się jedynie z \overline{p} początkowych kolumn macierzy Y_{0}, przy czym 1\le\overline{p}<p, to otrzymane w procesie iteracyjnym macierze \overline{Y}_{k} i \overline{Z}_{k} (formatu n\times\overline{p} i \overline{p}\times\overline{p}) byłyby odpowiednio ,,okrojonymi” macierzami Y_{k} i Z_{k}. Stąd, jeśli wszystkie wartości własne macierzy A mają różne moduły, to dla każdego takiego \overline{p} podprzestrzeń {\mathcal{Y}}_{k}^{{(\overline{p})}} rozpięta na początkowych \overline{p} wektorach macierzy Z_{k} ,,zbiega” do podprzestrzeni własnej {\mathcal{P}}^{{(\overline{p})}}, przy czym

\tg\angle\left({\mathcal{Y}}_{k}^{{(\overline{p})}},{\mathcal{P}}^{{(\overline{p})}}\right)\le\rho _{{\overline{p}}}^{k}\cdot\tg\angle\left({\mathcal{Y}}_{0}^{{(\overline{p})}},{\mathcal{P}}^{{\overline{p}}}\right).

Z powyższej obserwacji wynika, że gdybyśmy realizowali iteracje podprzestrzeni przyjmując p=n i startując z macierzy identycznościowej Z_{0}:=I formatu n\times n to mielibyśmy zbieżność {\mathcal{Y}}_{k}^{{(p)}} do podprzestrzeni własnych {\mathcal{P}}^{{(p)}} dla wszystkich p=1,2,\ldots,n.

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

V=[\vec{\xi}_{1},\vec{\xi}_{2},\ldots,\vec{\xi}_{n}],\qquad Z_{k}=[\vec{z}_{1}^{{(k)}},\vec{z}_{2}^{{(k)}},\ldots,\vec{z}_{n}^{{(k)}}].
Lemat 3.1

Niech B_{k} będzie macierzą przejścia od bazy (\vec{\xi}_{1},\ldots,\vec{\xi}_{n}) do (\vec{z}_{1}^{{(k)}},\ldots,\vec{z}_{n}^{{(k)}}),

Z_{k}=V\, B_{k}.

Dla 1\le p\le n-1, niech

B_{k}=\left[\begin{array}[]{cc}B_{{1,1}}^{{(k)}}&B_{{1,2}}^{{(k)}}\\
B_{{2,1}}^{{(k)}}&B_{{2,2}}^{{(k)}}\end{array}\right],

gdzie B_{{1,1}}^{{(k)}} jest podmacierzą formatu p\times p, a pozostałe podmacierze odpowiednio formatów p\times(n-p), (n-p)\times p, (n-p)\times(n-p). Jeśli \rho _{p}=|\lambda _{{p+1}}|/|\lambda _{p}|<1 i \gamma _{0}<\pi/2 to

\| B_{{1,2}}^{{(k)}}\| _{2},\| B_{{2,1}}^{{(k)}}\| _{2}\,\le\,\rho _{p}^{k}\cdot\tg\gamma _{0}.
Dowód

Niech Z_{k}=[Z_{1}^{{(k)}},Z_{2}^{{(k)}}], V=[V_{1},V_{2}], gdzie Z_{1}^{{(k)}} i V_{1} są formatu n\times p. Wtedy

Z_{1}^{{(k)}}=V_{1}\, B_{{1,2}}^{{(k)}}+V_{2}\, B_{{2,1}}^{{(k)}}.

Mnożąc to równanie z lewej strony przez V_{2}^{T} dostajemy B_{{2,1}}^{{(k)}}=V_{2}^{T}\, Z_{1}^{{(k)}}. Stąd

\| B_{{2,1}}^{{(k)}}\| _{2}\,=\,\| V_{2}^{T}\, Z_{1}^{{(k)}}\| _{2}\,=\,\sup _{{\|\vec{x}\| _{2}=1}}\| V_{2}^{T}\,(Z_{1}^{{(k)}}\,\vec{x})\| _{2}.

Oznaczając \vec{y}=Z_{1}^{{(k)}}\,\vec{x} zauważamy, że \vec{y}\in{\mathcal{Y}}_{k}:={\rm span}(\vec{z}_{1}^{{(k)}},\ldots,\vec{z}_{p}^{{(k)}}) i \|\vec{y}\| _{2}=1. Stąd wektor V_{2}^{T}\,\vec{y} zawiera współczynniki rzutu ortogonalnego \vec{y} na ({\mathcal{P}}^{{(p)}})^{\perp}:={\rm span}(\vec{\xi}_{{p+1}},\ldots,\vec{\xi}_{n}) w bazie (\vec{\xi}_{{p+1}},\ldots,\vec{\xi}_{n}). Z definicji normy drugiej i twierdzenia 3.1 dostajemy więc

\| B_{{2,1}}^{{(k)}}\| _{2}=\sup _{{\|\vec{x}\| _{2}=1}}\sin\angle\left(Z_{1}^{{(k)}}\,\vec{x},({\mathcal{P}}^{{(p)}})^{\perp}\right)\le\rho _{p}^{k}\cdot\tg\gamma _{0}.

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

\sin\angle({\mathcal{Y}}_{k},{\mathcal{P}}^{{(p)}})=\sin\angle\left({\mathcal{Y}}_{k}^{\perp},({\mathcal{P}}^{{(p)}})^{\perp}\right)

dostajemy taką samą nierówność dla \| B_{{1,2}}^{{(k)}}\| _{2}.

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

A_{k}\,:=\, Z_{k}^{T}\, A\, Z_{k},

która powinna zbiegać do \Lambda={\rm diag}(\lambda _{1},\ldots,\lambda _{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ą Z_{0}=I. Wtedy A\, Z_{{k-1}}=Z_{k}\, R_{k}, skąd

A_{{k-1}}=Z_{{k-1}}^{T}\, A\, Z_{{k-1}}=(Z_{{k-1}}^{T}\, Z_{k})\, R_{k}.

Oznaczając Q_{k}:=Z_{{k-1}}^{T}\, Z_{k} widzimy, że ostatnia równość to nic innego jak rozkład ortogonalno-trójkątny macierzy A_{{k-1}},

A_{{k-1}}=Q_{k}\, R_{k}.

Prawdziwe są więc związki rekurencyjne

\displaystyle A_{k} \displaystyle= \displaystyle Z_{k}^{T}\, A\, Z_{k}=(Z_{k}^{T}\, Z_{{k-1}})\,(Z_{{k-1}}^{T}\, A\, Z_{{k-1}})\,(Z_{{k-1}}^{T}\, Z_{k})
\displaystyle= \displaystyle Q_{k}^{T}\, A_{{k-1}}\, Q_{k},
\displaystyle Z_{k} \displaystyle= \displaystyle Z_{{k-1}}\,(Z_{{k-1}}^{T}\, Z_{k})=Z_{{k-1}}\, Q_{k}.

Zależności te prowadzą do ALGORYTMU QR, który oblicza kolejne macierze A_{k} i Z_{k} w następujący sposób.

\displaystyle A_{0}:=A;\quad Z_{0}:=I;
\displaystyle\mbox{dla}\quad k=1,2,3,\ldots
(QR) \displaystyle\left\{\begin{array}[]{lll}A_{{k-1}}:=Q_{k}\, R_{k};\quad\mbox{(ortonormalizacja)},\\
A_{k}:=R_{k}\, Q_{k};\quad Z_{k}:=Z_{{k-1}}\, Q_{k};\end{array}\right.

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 Z_{k} (a tym samym także A_{k}) 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 \widetilde{Z}_{k}, \widetilde{A}_{k} i Z_{k}, A_{k} będą macierzami powstałymi w wyniku realizacji odpowiednio algorytmów (IO) i (QR). Wtedy

\widetilde{Z}_{k}\,=\, Z_{k}\, D_{k}\qquad\mbox{oraz}\qquad\widetilde{A}_{k}=D_{k}\, A_{k}\, D_{k},

gdzie D_{k} są pewnymi macierzami diagonalnymi z elementami \pm 1 na głównej przekątnej.

Dowód

Zastosujemy indukcję względem k. Dla k=1 mamy \widetilde{Z}_{0}=I=Z_{0} oraz \widetilde{A}_{0}=A=A_{0}. Załóżmy, że lemat jest prawdziwy dla k-1, tzn. \widetilde{Z}_{{k-1}}=Z_{k}\, D_{{k-1}}. Wtedy, z jednej strony

\widetilde{A}_{{k-1}}\,=\, D_{{k-1}}\, A_{{k-1}}\, D_{{k-1}}\,=\, D_{{k-1}}\, Q_{k}\, R_{k}\, D_{{k-1}},

a z drugiej strony

\widetilde{A}_{{k-1}}\,=\,\widetilde{Z}_{{k-1}}^{T}\,\widetilde{Z}_{k}\,\widetilde{R}_{k}\\
\,=\, D_{{k-1}}\, Z_{{k-1}}^{T}\,\tilde{Z}_{k}\,\tilde{R}_{k}.

Mamy więc

D_{{k-1}}\,\widetilde{A}_{{k-1}}\,=\,(Z_{{k-1}}^{T}\,\tilde{Z}_{k})\,\tilde{R}_{k}\,=\, Q_{k}\,(R_{k}\, D_{{k-1}}),

przy czym równości te przedstawiają dwa rozkłady ortogonalno-trójkątne tej samej macierzy D_{{k-1}}\,\widetilde{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 D_{k} z elementami \pm 1 na głównej przekątnej taka, że

Z_{{k-1}}^{T}\,\widetilde{Z}_{k}\,=\, Q\, D_{k}\qquad(\mbox{oraz}\quad\widetilde{R}_{k}=D_{k}\, R_{k}\, D_{{k-1}}).

Stąd

\widetilde{Z}_{k}\,=\, Z_{{k-1}}\, Q_{k}\, D_{k}\,=\, Z_{k}\, D_{k}.

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

\widetilde{A}_{k}\,=\,\widetilde{Z}_{k}^{T}\, A\,\widetilde{Z}_{k}\,=\, D_{k}^{T}\,(Z_{k}^{T}\, A\, Z_{k})\, D_{k}\,=\, D_{k}\, A_{k}\, D_{k}.

Zauważmy jeszcze, że jeśli przez \widetilde{a}_{{i,j}}^{{(k)}} i a_{{i,j}}^{{(k)}} oznaczymy odpowiednio elementy macierzy \widetilde{A}_{k} i A_{k} oraz D_{k}={\rm diag}(d_{1}^{{(k)}},\ldots,d_{n}^{{(k)}}) z d_{i}^{{(k)}}=\pm 1 to

\widetilde{a}_{{i,j}}^{{(k)}}=d_{i}^{{(k)}}d_{j}^{{(k)}}a_{{i,j}}^{{(k)}}.

To oznacza, że dla i\ne 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\le p\le n, niech \rho _{p}:=|\lambda _{{p+1}}/\lambda _{p}|<1 i \gamma _{k}:=\angle\left({\mathcal{Y}}_{k},{\mathcal{P}}^{{(p)}}\right). Niech dalej

A_{k}=\left[\begin{array}[]{cc}A_{{1,1}}^{{(k)}}&A_{{1,2}}^{{(k)}}\\
A_{{2,1}}^{{(k)}}&A_{{2,2}}^{{(k)}}\end{array}\right],\quad Q_{k}=\left[\begin{array}[]{cc}Q_{{1,1}}^{{(k)}}&Q_{{1,2}}^{{(k)}}\\
Q_{{2,1}}^{{(k)}}&Q_{{2,2}}^{{(k)}}\end{array}\right],\quad R_{k}=\left[\begin{array}[]{cc}R_{{1,1}}^{{(k)}}&R_{{1,2}}^{{(k)}}\\
0&R_{{2,2}}^{{(k)}}\end{array}\right],

gdzie podmacierze A_{{1,1}}^{{(k)}},Q_{{1,1}}^{{(k)}},R_{{1,1}}^{{(k)}} są formatu p\times p, a pozostałe podmacierze formatów odpowiednio p\times(n-p), (n-p)\times p i (n-p)\times(n-p). Wtedy, przyjmując

\varepsilon _{p}^{{(k)}}=\rho _{p}^{k}\cdot\tg\gamma _{0}

mamy \tg\gamma _{k}\le\varepsilon _{p}^{{(k)}}\tg oraz

\displaystyle\| A_{{1,2}}^{{(k)}}\| _{2}=\| A_{{2,1}}^{{(k)}}\| _{2} \displaystyle\le \displaystyle 2\,\varepsilon _{p}^{{(k)}}\| A\| _{2}, (3.2)
\displaystyle\| Q_{{1,2}}^{{(k)}}\| _{2},\,\| Q_{{2,1}}^{{(k)}}\| _{2} \displaystyle\le \displaystyle 2\,\varepsilon _{p}^{{(k)}}, (3.3)
\displaystyle\| R_{{1,2}}^{{(k)}}\| _{2} \displaystyle\le \displaystyle 4\,\varepsilon _{p}^{{(k)}}\| A\| _{2}. (3.4)
Dowód

Wobec wykazanej w lemacie \tg (teoretycznej) równoważności metody QR i iteracji (IO), nierówność \tg\gamma _{k}\le\varepsilon _{p}^{{(k)}} została już udowodniona w twierdzeniu 3.1. Aby pokazać (3.3), zapiszemy Z_{k} w postaci Z_{k}=[Z_{1}^{{(k)}},Z_{2}^{{(k)}}], gdzie Z_{1}^{{(k)}} i Z_{2}^{{(k)}} składają się odpowiednio z pierwszych p i z ostatnich (n-p) kolumn macierzy Z_{k}. Wtedy

Q_{{2,1}}=Z_{2}^{{(k-1)}}\, Z_{1}^{{(k)}}=\left(V^{T}\, Z_{2}^{{(k-1)}}\right)^{T}\,\left(V^{T}\, Z_{1}^{{(k)}}\right),

gdzie V=[\vec{\xi}_{1},\vec{\xi}_{2},\ldots,\vec{\xi}_{n}]. Pisząc V=[V_{1},V_{2}], podobnie jak dla Z_{k}, mamy teraz

V^{T}\, Z_{1}^{{(k)}}=\left[\begin{array}[]{cc}V_{1}^{T}\, Z_{1}^{{(k)}}\\
V_{2}^{T}\, Z_{1}^{{(k)}}\end{array}\right]=:\left[\begin{array}[]{cc}F\\
G\end{array}\right].

Ponieważ \| V_{1}\| _{2}=\| Z_{1}^{{(k)}}\| _{2}=1 to \| F\| _{2}\le 1, natomiast \| G\| _{2}\le\varepsilon _{p}^{{(k)}}, jak pokazaliśmy w dowodzie lematu 3.1.

Pisząc z kolei

V^{T}\, Z_{2}^{{(k-1)}}=:\left[\begin{array}[]{cc}F^{{\prime}}\\
G^{{\prime}}\end{array}\right],

w analogiczny sposób pokazujemy, że \| F^{{\prime}}\| _{2}\le 1 i \| G^{{\prime}}\| _{2}\le\varepsilon _{p}^{{(k)}}. Stąd mamy

\| Q_{{2,1}}^{{(k)}}\| _{2}=\| G^{T}\, F+F^{T}\, G\| _{2}\le\| G^{{\prime}}\| _{2}\| F\|+\| F^{{\prime}}\| _{2}\| G\| _{2}\le\varepsilon _{{k-1}}+\varepsilon _{p}^{{(k)}}\le 2\cdot\varepsilon _{p}^{{(k)}}.

Dokładnie tak samo pokazujemy oszacowanie dla \| Q_{{1,2}}\| _{2}.

Aby pokazać (3.2) zauważmy, że wobec A_{k}=Q_{k}\, R_{k} mamy też A_{{2,1}}^{{(k)}}=Q_{{2,1}}^{{(k)}}\, R_{{1,1}}^{{(k)}}. Ale

\| R_{{1,1}}\| _{2}\le\| R_{k}\| _{2}=\| A_{2}\| _{2}=\| A\| _{2},

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

W końcu, wobec R_{k}=A_{k}\, Q_{k}^{T} mamy

R_{{1,2}}^{{(k)}}=A_{{1,1}}^{{(k)}}\,{Q_{{2,1}}^{{(k)}}}^{T}+A_{{1,2}}^{{(k)}}\,{Q_{{2,2}}^{{(k)}}}^{T}.

Nierówność (3.4) wynika teraz z nierówności \| A_{{1,1}}^{{(k)}}\| _{2}\le\| A_{k}\| _{2}=\| A\| _{2}, (3.2), (3.3), oraz \| Q_{{2,2}}^{{(k)}}\| _{2}\le 1.

Z twierdzenia 3.2 wynika, że elementy pozadiagonalne macierzy A_{k} zbiegają do zera liniowo, przy czym iloraz zbieżności zależy od wartości \rho _{p}=|\lambda _{{p+1}}|/|\lambda _{p}|. Dokładniej, dla danego elementu a_{{i,j}}^{{(k)}} macierzy A_{k}, gdzie 1\le i<j\le n, prawdziwe jest oszacowanie

|a_{{i,j}}^{{(k)}}|\le\min\left(\varepsilon _{i}^{{(k)}},\varepsilon _{{i+1}}^{{(k)}},\ldots\varepsilon _{{j-1}}^{{(k)}}\right)

(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 A_{k} do wartości własnych macierzy A? Oznaczmy

\overline{\varepsilon}_{p}^{{(k)}}=\left\{\begin{array}[]{ll}\varepsilon _{1}^{{(k)}}&p=1,\\
\max\left(\varepsilon _{{p-1}}^{{(k)}},\varepsilon _{p}^{{(k)}}\right)&2\le p\le n-1,\\
\varepsilon _{{n-1}}^{{(k)}}&p=n.\end{array}\right.
Twierdzenie 3.3

Dla wszystkich 1\le p\le n mamy

|a_{{p,p}}^{{(k)}}-\lambda _{p}|\,\le\, 4\cdot\left(\overline{\varepsilon}_{p}^{{(k)}}\right)^{2}\cdot\| A\| _{2}.
Dowód

Wobec Z_{k}=V\, B_{k} mamy

\displaystyle a_{{p,p}}^{{(k)}} \displaystyle= \displaystyle(\vec{z}_{p}^{{(k)}})^{T}\, A\,\vec{z}_{p}^{{(k)}}\,=\,\left(V\,\vec{b}_{p}^{{(k)}}\right)^{T}\, A\,\left(V\,\vec{b}_{p}^{{(k)}}\right)
\displaystyle= \displaystyle\vec{b}_{p}^{{(k)}}\, V^{T}\, A\, V\,\vec{b}_{p}^{{(k)}}\,=\,\vec{b}_{p}^{{(k)}}\,\Lambda\,\vec{b}_{p}^{{(k)}}
\displaystyle= \displaystyle\lambda _{p}\left(b_{{p,p}}^{{(k)}}\right)^{2}+\sum _{{i\ne p}}\lambda _{i}\left(b_{{i,j}}^{{(k)}}\right)^{2}.

Stąd i z lematu 3.1 otrzymujemy

\displaystyle|a_{{p,p}}^{{(k)}}-\lambda _{p}| \displaystyle= \displaystyle\left|-\lambda _{p}\left(1-\left(b_{{p,p}}^{{(k)}}\right)^{2}\right)+\sum _{{i\ne p}}\lambda _{i}\left(b_{{i,p}}^{{(k)}}\right)^{2}\right|\,=\,\left|\sum _{{i\ne p}}\left(b_{{i,p}}^{{(k)}}\right)^{2}(\lambda _{i}-\lambda _{p})\right|
\displaystyle\le \displaystyle 2\cdot\| A\| _{2}\cdot\left(\sum _{{i<p}}\left(b_{{i,p}}^{{(k)}}\right)^{2}+\sum _{{p<i}}\left(b_{{i,p}}^{{(k)}}\right)^{2}\right)\,\le\, 4\cdot\left(\overline{\varepsilon}_{p}^{{(k)}}\right)^{2}\cdot\| A\| _{2}.

Elementy diagonalne a_{{p,p}}^{{(k)}} macierzy A_{k} zbiegają więc do wartości własnych \lambda _{p} macierzy A co najmniej jak |\lambda _{2}/\lambda _{1}|^{{2k}} dla p=1, |\lambda _{n}/\lambda _{{n-1}}|^{{2k}} dla p=n, oraz \min\left\{|\lambda _{p}/\lambda _{{p-1}}|^{{2k}},|\lambda _{{p+1}}/\lambda _{p}|^{{2k}}\right\} dla 2\le p\le 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 \overline{A}=A-\sigma\, I. Oznaczmy przez \overline{A}_{k} kolejne macierze powstające w tym procesie. Z wcześniejszej analizy metody QR wynika, że jeśli

|\lambda _{{i^{*}}}-\sigma|<\min _{{i\ne i^{*}}}|\lambda _{i}-\sigma|

to wyraz \overline{a}_{{n,n}}^{{(k)}} w prawym dolnym rogu macierzy \overline{A}_{k} zbiega do \lambda _{{i^{*}}}-\sigma i zbieżność jest tym lepsza im \sigma jest bliższe \lambda _{{i^{*}}}. Dlatego ważne jest, aby \sigma 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:

\displaystyle A_{0}:=A;\quad Z_{0}:=I;
\displaystyle\mbox{dla}\quad k=1,2,3,\ldots
(QR-shift) \displaystyle\left\{\begin{array}[]{lll}A_{{k-1}}-\sigma _{k}\, I:=Q_{k}\, R_{k};\quad\mbox{(ortonormalizacja)}\\
A_{k}:=R_{k}\, Q_{k}+\sigma _{k}\, I;\quad Z_{k}:=Z_{{k-1}}\, Q_{k};\end{array}\right.

Jasne, że

A_{k}=Q_{k}^{T}\,(A_{{k-1}}-\sigma _{k}\, I)\, Q_{k}+\sigma _{k}\, I\,=\, Q_{k}^{T}\, A_{{k-1}}\, Q_{k},

tzn. kolejne macierze A_{k} są ortogonalnie podobne. Okazuje się, że QR ma jeszcze jedną ważną własność. Oznaczmy przez \vec{z}_{n}^{{(k)}} ostatnią kolumnę macierzy Z_{k}, przez r_{{n,n}}^{{(k)}} wyraz w prawym dolnym rogu macierzy trójkątnej górnej R_{k}, oraz przez \vec{e}_{n} ostatni wersor. Wtedy dla każdego k mamy

\displaystyle(A-\sigma _{k}\, I)\,\vec{z}_{n}^{{(k)}} \displaystyle= \displaystyle Z_{{k-1}}\,(A_{{k-1}}-\sigma _{{k-1}}\, I)^{T}\, Z_{{k-1}}^{T}\,\vec{z}_{n}^{{(k)}}
\displaystyle= \displaystyle Z_{{k-1}}\,(Q_{k}\, R_{k})^{T}\, Z_{{k-1}}^{T}\,\vec{z}_{n}^{{(k)}}
\displaystyle= \displaystyle Z_{{k-1}}\, R_{k}^{T}\, Z_{k}^{T}\,\vec{z}_{n}^{{(k)}}\;=\; Z_{{k-1}}\, R_{k}^{T}\,\vec{e}_{n}
\displaystyle= \displaystyle r_{{n,n}}^{{(k)}}\, Z_{{k-1}}\,\vec{e}_{n}\;=\; r_{{n,n}}^{{(k)}}\,\vec{z}_{n}^{{(k-1)}},

czyli

\frac{z_{n}^{{(k)}}}{r_{{n,n}}^{{(k)}}}=(A-\sigma _{k}\, I)^{{-1}}\,\vec{z}_{n}^{{(k-1)}}.

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

\sigma _{k}=(\vec{z}_{n}^{{(k)}})^{T}\, A\,\vec{z}_{n}^{{(k)}}=\vec{e}_{n}^{T}\, Z_{k}^{T}\, A\, Z_{k}\,\vec{e}_{n}=\vec{e}_{n}\, A_{k}\,\vec{e}_{n}=a_{{n,n}}^{{(k)}}.

Przy tak dobranym \sigma _{k} dostajemy asymptotycznie zbieżność sześcienną elementu a_{{n,n}}^{{(k)}} do \lambda _{{i^{*}}}.

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

\left[\begin{array}[]{ll}a_{{n-1,n-1}}^{{(k)}}&a_{{n-1,n}}^{{(k)}}\\
a_{{n,n-1}}^{{(k)}}&a_{{n,n}}^{{(k)}}\end{array}\right],

która jest najbliższa a_{{n,n}}^{{(k)}}, a jeśli obie są jednakowo oddalone od a_{{n,n}}^{{(k)}} 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 \lambda _{{i^{*}}} (co zwykle osiąga się wykonując kilka kroków QR i sprawdza przez obliczenie residuum \| A\,\vec{z}_{n}^{{(k)}}-a_{{n,n}}^{{(k)}}\,\vec{z}_{n}^{{(k)}}\| _{2}) to uznajemy, że wyrazy w ostatniej kolumnie i ostatnim wierszu macierzy A_{k}, poza wyrazem a_{{n,n}}^{{(k)}}, są zerowe i redukujemy problem do macierzy formatu (n-1)\times(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.

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.