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 ${\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.

	$\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}.$

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