Zagadnienia

2.1 Uwagi wstępne
- 2.1.1 Sprowadzanie macierzy do postaci Hessenberga
2.2 Metoda Hymana
2.3 Metoda potęgowa
- 2.3.1 Definicja metody
- 2.3.2 Analiza zbieżności
2.4 Odwrotna metoda potęgowa i deflacja

2. Zagadnienie własne II

Kolejne dwa wykłady poświęcimy konkretnym metodom numerycznym rozwiązywania zagadnienia własnego.

2.1. Uwagi wstępne

Ponieważ wartości własne macierzy $A$ są zerami jej wielomianu charakterystycznego, narzucająca się metoda poszukiwania wartości własnych polegałaby na wyliczeniu współczynników tego wielomianu, na przykład w bazie potęgowej, a następnie na zastosowaniu jednej ze znanych metod znajdowania zer wielomianu. Dla wielomianów o współczynnikach rzeczywistych mogłaby to być np. metoda bisekcji, metoda Newtona (siecznych), albo pewna kombinacja obu method. Należy jednak przestrzec przed dokładnie takim postępowaniem. Rzecz w tym, że zera wielomianu mogą być bardzo wrażliwe na zaburzenia jego współczynnikówi. Dobrze ilustruje to następujący przykład.

Przykład 2.1

Załóżmy, że $A$ jest macierzą symetryczną formatu $20\times 20$ o wartościach własnych $1,2,\ldots,20$ . Zapisując wielomian charakterystyczny w postaci potęgowej mamy

${\rm det}(A-\lambda\, I)\,=\,\prod _{{k=1}}^{{20}}(\lambda _{k}-k)\,=\,\sum _{{k=0}}^{{20}}w_{k}\lambda^{k}.$

Okazuje się, że zaburzenie współczynnika $w_{{19}}$ mogą powodować $10^{{10}}$ razy większe zaburzenie wartości własnej $\lambda _{{16}}=16$ .

Nie znaczy to oczywiście, że metody bazujące na obliczaniu wielomianu charakterystycznego, czy jego pochodnych trzeba z gruntu odrzucić. Trzeba tylko pamiętać, aby przy obliczeniach korzystać bezpośrednio ze współczynników $a_{{i,j}}$ macierzy $A$ , ponieważ, jak wiemy z poprzedniego rozdziału, zadanie jest ze względu na te współczynniki dobrze uwarunkowane.

Z drugiej strony zauważmy, że policzenie wyznacznika macierzy pełnej wprost z definicji jest raczej kosztowne. Dlatego w praktyce metody wyznacznikowe są zwykle poprzedzone prekomputingiem polegającym na sprowadzeniu macierzy przez podobieństwa ortogonalne do prostszej postaci, z której wyznacznik można już obliczyć dużo tańszym kosztem niż dla macierzy pełnej. Tego typu precomputing ma również zastosowanie w innych metodach, w którym np. trzeba wykonywać mnożenie macierzy przez wektor.

Jedną z takich wygodnych postaci macierzy jest postać Hessenberga, a dla macierzy hermitowskich postać trójdiagonalna.

2.1.1. Sprowadzanie macierzy do postaci Hessenberga

Macierz $A=(a_{{i,j}})\in\mathbb{C}^{{n,n}}$ jest postaci Hessenberga (,,prawie” trójkątną górną) jeśli wszystkie wyrazy $a_{{i,j}}$ dla $i\ge j+2$ są zerami, tzn.

$A\,=\,\left[\begin{array}[]{ccccc}a_{{1,1}}&a_{{1,2}}&\cdots&a_{{1,n-1}}&a_{{1,n}}\\ a_{{2,1}}&a_{{2,2}}&\cdots&a_{{2,n-1}}&a_{{2,n}}\\ 0&a_{{3,2}}&\cdots&a_{{3,n-1}}&a_{{3,n}}\\ \vdots&\vdots&&\vdots&\vdots\\ 0&0&\cdots&a_{{n,n-1}}&a_{{n,n}}\end{array}\right].$

Oczywiście, jeśli macierz jest hermitowska, $A=A^{H}$ , to postać Hessenberga jest równoważna postaci trójdiagonalnej

$A\,=\,\left[\begin{array}[]{ccccc}c_{1}&b_{2}&0&\cdots&0\\ \overline{b}_{2}&c_{2}&b_{3}&\cdots&0\\ 0&\overline{b}_{3}&c_{3}&\cdots&0\\ \vdots&\vdots&\vdots&&\vdots\\ 0&0&\cdots&\cdots&c_{n}\end{array}\right].$

(2.1)

Aby daną macierz $A$ sprowadzić do postaci Hessenberga przez podobieństwa ortogonalne (a więc nie zmieniając wartości własnych), możemy zastosować znane nam odbicia Householdera. Przypomnijmy, że dla dowolnego niezerowego wektora $\vec{a}\in\mathbb{C}^{n}$ istnieje macierz ortogonalna (odbicie Householdera) $P=P^{H}=P^{{-1}}\in\mathbb{C}^{{n,n}}$ postaci $P=I-2\,\vec{u}\,\vec{u}^{{\, H}}$ , gdzie $\|\vec{u}\| _{2}=1$ , która przekształca $\vec{a}$ na kierunek pierwszego wersora, tzn. $P\,\vec{a}=\ell\,\vec{e}_{1}$ , $|\ell|=\|\vec{a}\| _{2}$ . (Odbicia Householdera zostały dokładnie przedstawione na wykładzie Analiza Numeryczna I.)

Niech $A=(a_{{i,j}})\in\mathbb{C}^{{n,n}}$ będzie dowolną macierzą. Algorytm konstruuje najpierw odbicie $P_{1}\in\mathbb{C}^{{n-1,n-1}}$ dla wektora $[a_{{2,1}},a_{{3,1}},\ldots,a_{{n,1}}]^{T}$ , a następnie biorąc macierz

$\widehat{P}_{1}=\left[\begin{array}[]{cc}1&\vec{0}^{{\, T}}\\ \vec{0}&P_{1}\end{array}\right]\in\mathbb{C}^{{n,n}}$

oblicza $A^{{(1)}}=\widehat{P}_{1}\, A\,\widehat{P}_{1}^{H}$ . łatwo zobaczyć, że wtedy elementy $(i,1)$ dla $i=3,4,\ldots,n$ w macierzy $A^{{(1)}}$ są równe zeru.

Postępując indukcyjnie załóżmy, że dostaliśmy już macierz $A^{{(k)}}=(a^{{(k)}}_{{i,j}})$ , $1\le k\le n-3$ , w której elementy $a^{{(k)}}_{{i,j}}$ dla $i\ge j+2$ , $1\le j\le k$ , są wyzerowane. W kroku $(k+1)$ -szym algorytm konstruuje odbicie $P_{{k+1}}\in\mathbb{C}^{{n-k,n-k}}$ dla wektora $[a^{{(k)}}_{{k+2,k+1}},\ldots,a^{{(k)}}_{{n,k+1}}]^{T}$ , a następnie biorąc macierz

$\widehat{P}_{{k+1}}=\left[\begin{array}[]{cc}I_{k}&0^{{\, T}}\\ 0&P_{{k+1}}\end{array}\right]\in\mathbb{C}^{{n,n}}$

oblicza $A^{{(k+1)}}=\widehat{P}_{{k+1}}\, A^{{(k)}}\,\widehat{P}_{{k+1}}^{H}$ . Wtedy wszystkie elementy $a^{{(k+1)}}_{{i,j}}$ dla $i\ge j+2$ , $1\le j\le k+1$ , są zerami. Po $(n-2)$ krokach otrzymana macierz $\widetilde{A}=A^{{(n-1)}}$ jest postaci Hessenberga.

Jasne jest, że jeśli $A=A^{H}$ to opisany algorytm prowadzi do macierzy trójdiagonalnej, $\widetilde{A}=\widetilde{A}^{H}$ . (Obliczenia można wtedy w każdym kroku wykonywać jedynie na głównej diagonali i pod nią.)

Dodajmy jeszcze, że koszt sprowadzenia macierzy przez podobieństwa ortogonalne do postaci Hessenberga/trójdiagonalnej jest proporcjonalny do $n^{3}$ .

2.2. Metoda Hymana

Przy pomocy metody Hymana możemy w zręczny sposób obliczyć wartości i pochodne wielomianu charakterystycznego $\det(A-\lambda\, I)$ dla macierzy Hessenberga $A=(a_{{i,j}})\in\mathbb{C}^{{n,n}}.$ Bez zmniejszenia ogólności będziemy zakładać, że wszystkie elementy $a_{{i+1,i}}\ne 0$ . W przeciwnym przypadku wyznacznik można obliczyć jako iloczyn wyznaczników macierzy Hessenberga niższych rzędów.

Metoda Hymana polega na dodaniu do pierwszego wiersza macierzy $A-\lambda\, I$ wiersza $i$ -tego pomnożonego przez pewien współczynnik $q_{i}=q_{i}(\lambda)$ , dla $i=2,3,\ldots,n$ tak, aby wyzerować elementy $(1,i)$ dla $i=1,2,\ldots,n-1$ . Ponieważ

$A-\lambda\, I\,=\,\left[\begin{array}[]{ccccc}a_{{1,1}}-\lambda&a_{{1,2}}&\cdots&a_{{1,n-1}}&a_{{1,n}}\\ a_{{2,1}}&a_{{2,2}}-\lambda&\cdots&a_{{2,n-1}}&a_{{2,n}}\\ 0&a_{{3,2}}&\cdots&a_{{3,n-1}}&a_{{3,n}}\\ \vdots&\vdots&&\vdots&\vdots\\ 0&0&\cdots&a_{{n,n-1}}&a_{{n,n}}-\lambda\end{array}\right],$

mamy następujące równania na $q_{i}$ :

	$\displaystyle(a_{{1,1}}-\lambda)\,+\, a_{{1,2}}q_{2}$	$\displaystyle=$	$\displaystyle 0,$
	$\displaystyle a_{{1,1}}\,+\, a_{{2,1}}q_{2}\,+\,\cdots\,+\, a_{{i-1,i}}q_{{i-1}}\,+\,(a_{{i,i}}-\lambda)q_{i}\,+\, a_{{i+1,i}}q_{{i+1}}$	$\displaystyle=$	$\displaystyle 0,$		(2.2)

dla $i=2,3,\ldots,n-1.$ Stąd, definiując dodatkowo $q_{1}=1$ , dostajemy następujące równania rekurencyjne:

	$\displaystyle q_{2}(\lambda)$	$\displaystyle=$	$\displaystyle\frac{-(a_{{1,1}}-\lambda)}{a_{{2,1}}},$
	$\displaystyle q_{{i+1}}(\lambda)$	$\displaystyle=$	$\displaystyle\frac{-\left(a_{{1,i}}q_{1}+\cdots+a_{{i-1,i}}q_{{i-1}}+(a_{{i,i}}-\lambda)q_{i}\right)}{a_{{i+1,i}}}.$

Po opisanym przekształceniu macierzy $A-\lambda\, I$ zmieni się jedynie jej pierwszy wiersz; wyrazy $a_{{1,i}}$ , $1\le i\le n-1$ , zostaną wyzerowane, a $a_{{1,n}}$ zostanie przekształcony do

$\widehat{a}_{{1,n}}\,=\, a_{{1,n}}q_{1}+\cdots+a_{{n-1,n}}q_{{n-1}}+(a_{{n,n}}-\lambda)q_{n}.$

Rozwijając wyznacznik względem (przekształconego) pierwszego wiersza otrzymujemy

$\det(A-\lambda\, I)\,=\,(-1)^{{n+1}}a_{{2,1}}a_{{3,2}}\cdots a_{{n,n-1}}\left(a_{{1,n}}q_{1}+\cdots+a_{{n-1,n}}q_{{n-1}}+(a_{{n,n}}-\lambda)q_{n}\right),$

a stąd zera wielomianu charakterystycznego są równe zerom wielomianu

$q_{{n+1}}(\lambda)\,=\,-\left(a_{{1,n}}q_{1}+\cdots+a_{{n-1,n}}q_{{n-1}}+(a_{{n,n}}-\lambda)q_{n}\right).$

Oczywiście, wartości $q_{{n+1}}(\lambda)$ można obliczać stosując wzory rekurencyjne (2.2) przedłużając je o $i=n$ , przy dodatkowym formalnym podstawieniu $a_{{n+1,n}}=1.$

Aby móc zastowoać netodę Newtona (stycznych) do znalezienia zer wielomianu $q_{{n+1}}$ potrzebujemy również wiedzieć jak obliczać jego pochodne. To też nie jest problem, bo rekurencyjne wzory na pochodne można uzyskać po prostu różniczkując wzory (2.2). Dodajmy, że koszt obliczenia wartości $q_{{n+1}}(\lambda)$ i $q^{\prime}_{{n+1}}(\lambda)$ jest proporjonalny do $n^{2}$ , co jest istotnym zyskiem w porównaniu z kosztem $n^{3}$ obliczania wyznacznika pełnej macierzy za pomocą faktoryzacji metodą eliminacji Gaussa.

Formuły na obliczanie wyznacznika $\det(A-\lambda\, I)$ i jego pochodnych uproszczają się jeszcze bardziej gdy macierz $A$ jest dodatkowo hermitowska, czyli jest macierzą trójdiagonalną postaci (2.1). Wtedy, stosując zwykłe rozumowanie rekurencyjne, łatwo się przekonać, że kolejne minory główne (czyli wyznaczniki macierzy kątowych) spełniają zależności

			$\displaystyle p_{0}(\lambda)\,=\, 1,\qquad\quad p_{1}(\lambda)\,=\, c_{1}-\lambda,$
			$\displaystyle p_{i}(\lambda)\,=\,(c_{i}-\lambda)p_{{i-1}}(\lambda)-\|b_{i}\|^{2}p_{{i-2}}(\lambda)\qquad\mbox{dla}\quad i=2,3,\ldots,n.$

Różniczkując otrzymane wzory otrzymujemy formuły na pochodne kolejnych minorów po $\lambda$ . Wartości wielomianu $\det(A-\lambda\, I)=p_{n}(\lambda)$ oraz jego pochodnych można więc obliczać kosztem liniowym w $n$ .

Sprawę konkretnej implementacji metod iteracyjnych bisekcji i/lub Newtona do wyznaczenia zer wielomianu $\det(A-\lambda\, I)$ tutaj pomijamy. Ograniczymy się jedynie do stwierdzenia, że nie jest to rzecz całkiem trywialna.

2.3. Metoda potęgowa

2.3.1. Definicja metody

Metoda potęgowa zdefiniowana jest w następujący prosty sposób. Rozpoczynając od dowolnego wektora $\vec{x}_{0}\in\mathbb{C}^{n}$ o normie $\|\vec{x}_{0}\| _{2}=1$ obliczamy kolejno dla $k=0,1,2,\ldots$

$\vec{y}_{k}\,:=\, A\,\vec{x}_{k},\qquad\vec{x}_{{k+1}}\,:=\,\frac{\vec{y}_{k}}{\|\vec{y}_{k}\| _{2}}.$

Równoważnie możemy napisać

$\vec{x}_{k}\,=\,\frac{A^{k}\,\vec{x}_{0}}{\| A^{k}\,\vec{x}_{0}\| _{2}}.$

Wektory $\vec{x}_{k}$ stanowią kolejne przybliżenia wektora własnego. Odpowiadającą mu wartość własną przybliżamy wzorem

$\eta _{k}\,:=\,\vec{x}_{k}^{{\, H}}\, A\,\vec{x}_{k}\,=\,\vec{x}_{k}^{{\, H}}\,\vec{y}_{k}.$

2.3.2. Analiza zbieżności

Analizę metody potęgowej przeprowadzimy przy założeniu, że macierz $A$ jest diagonalizowalna. Oznaczmy przez $\mu _{i}$ , $1\le i\le s$ , różne wartości własne macierzy $A$ , a przez ${\mathcal{V}}_{i}$ odpowiadające im podprzestrzenie własne,

${\mathcal{V}}_{1}\oplus{\mathcal{V}}_{2}\oplus\cdots\oplus{\mathcal{V}}_{s}\,=\,\mathbb{C}^{n}.$

(2.3)

Założymy, że $\mu _{1}$ jest dominującą wartością własną oraz

$|\mu _{1}|\,>\,|\mu _{2}|\,>\,\cdots\,>\,|\mu _{s}|.$

Przedstawmy $\vec{x}_{0}$ jednoznacznie w postaci

$\vec{x}_{0}\,=\,\sum _{{j=1}}^{s}\vec{v}_{j}$

gdzie $v_{j}\in{\mathcal{V}}_{j}$ . Dla uproszczenia, będziemy zakładać, że każda ze składowych

$\vec{v}_{j}\,\ne\,\vec{0}.$

Podkreślmy, że nie jest to założenie ograniczające, bo w przeciwnym przypadku składowe zerowe po prostu ignorujemy w poniższej analizie zbieżności. Poza tym, jeśli wektor początkowy jest wybrany losowo, to teoretyczne prawdopodobieństwo takiego zdarzenia jest zerowe. Co więcej, nawet jeśli jedna ze składowych znika to wskutek błędów zaokrągleń w procesie obliczeniowym składowa ta w wektorze $\vec{x}_{1}$ będzie z pewnością niezerowa. (Mamy tu do czynienia z ciekawym zjawiskiem, kiedy błędy zaokrągleń pomagają!)

Prosty rachunek pokazuje, że

$A^{k}\,\vec{x}_{0}\,=\,\sum _{{j=1}}^{s}A^{k}\,\vec{v}_{j}\,=\,\sum _{{j=1}}^{s}\mu _{j}^{k}\,\vec{v}_{j}\,=\,\mu _{1}^{k}\,\left(\vec{v}_{1}+\sum _{{j=2}}^{s}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\,\vec{v}_{j}\right).$

Ponieważ każdy z ilorazów $\mu _{j}/\mu _{1}<1$ to $\vec{x}_{k}$ zbiega do wektora własnego $\vec{v}_{1}/\|\vec{v}_{1}\| _{2}$ , a $\eta _{k}$ do dominującej wartości własnej $\mu _{1}$ . Przyjrzyjmy się od czego zależy szybkość zbieżności.

Odległość ${\rm dist}(\vec{x}_{k},{\mathcal{V}}_{1})$ wektora $\vec{x}_{k}$ od podprzestrzeni ${\mathcal{V}}_{1}$ można oszacować z góry przez odległość $\vec{x}_{k}$ od ${\rm span}(\vec{v}_{1})$ , która z kolei jest równa długości rzutu $P_{1}\vec{x}_{k}$ tego wektora na podprzestrzeń ortogonalną do $\vec{v}_{1}$ ,

$P_{1}\vec{x}_{k}\,=\,\vec{x}_{k}\,-\,\frac{\vec{v}_{1}^{{\, H}}\,\vec{x}_{k}}{\|\vec{v}_{1}\| _{2}^{2}}\,\vec{v}_{1}.$

Oznaczając

$\beta _{k}\,:=\,\left\|\vec{v}_{1}+\sum _{{j=2}}^{s}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\,\vec{v}_{j}\right\| _{2}$

mamy $\lim _{{k\to\infty}}\beta _{k}=\|\vec{v}_{1}\| _{2}$ oraz

$\displaystyle P_{1}\vec{x}_{k}$	$\displaystyle=$	$\displaystyle\frac{1}{\beta _{k}}\,\left(\left(\vec{v}_{1}+\sum _{{j=2}}^{s}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\,\vec{v}_{j}\right)-\left(\\|\vec{v}_{1}\\| _{2}^{2}+\sum _{{j=2}}^{s}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\,\vec{v}_{1}^{{\, H}}\,\vec{v}_{j}\right)\,\frac{\vec{v}_{1}}{\\|\vec{v}_{1}\\| _{2}^{2}}\right)$
	$\displaystyle=$	$\displaystyle\frac{1}{\beta _{k}}\,\left(\sum _{{j=2}}^{s}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\,\left(\vec{v}_{j}-\frac{\vec{v}_{1}^{{\, H}}\,\vec{v}_{j}}{\\|\vec{v}_{1}\\| _{2}^{2}}\,\vec{v}_{1}\right)\right)$
	$\displaystyle\approx$	$\displaystyle\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\frac{1}{\\|\vec{v}_{1}\\| _{2}}\left(\vec{v}_{2}-\frac{\vec{v}_{1}^{{\, H}}\,\vec{v}_{2}}{\\|\vec{v}_{1}\\| _{2}^{2}}\,\vec{v}_{1}\right)\qquad\qquad(k\to\infty).$

Biorąc normę i stosując nierówność trójkąta dostajemy

${\rm dist}(\vec{x}_{k},{\mathcal{V}}_{1})\,=\,\| P_{1}\vec{x}_{k}\| _{2}\;\lessapprox\; 2\cdot\rho _{1}^{k}\cdot\frac{\|\vec{v}_{2}\| _{2}}{\|\vec{v}_{1}\| _{2}}\qquad\qquad(k\to\infty),$

gdzie

$\rho _{1}\,:=\,\frac{|\mu _{2}|}{|\mu _{1}|}\,<\, 1.$

Zbieżność jest więc liniowa z ilorazem $\rho _{1}$ .

Zobaczmy teraz jak bardzo $\eta _{k}$ różni się od $\mu _{1}$ . Mamy

	$\displaystyle\eta _{k}$	$\displaystyle=$	$\displaystyle\vec{x}_{k}^{{\, H}}\, A\,\vec{x}_{k}\,=\,\frac{\mu _{1}}{\beta^{2}}\,\left(\vec{v}_{1}^{{\, H}}+\sum _{{j=2}}^{s}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\,\vec{v}_{j}^{{\, H}}\right)\,\left(\vec{v}_{1}+\sum _{{j=2}}^{s}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{{k+1}}\,\vec{v}_{1}^{{\, H}}\,\vec{v}_{j}\right)$
		$\displaystyle=$	$\displaystyle\frac{\mu _{1}}{\beta^{2}}\left(\\|\vec{v}_{1}\\| _{2}^{2}+\sum _{{j=2}}^{2}\left(\frac{\mu _{j}}{\mu _{1}}\right)^{k}\left(\vec{v}_{j}^{{\, H}}\,\vec{v}_{1}+\frac{\mu _{j}}{\mu _{1}}\,\vec{v}_{1}^{{\, H}}\,\vec{v}_{j}\right)+\sum _{{i,j=2}}^{s}\left(\frac{\mu _{i}^{k}\mu _{j}^{{k+1}}}{\mu _{1}^{{2k+1}}}\right)\,\vec{v}_{i}^{{\, H}}\,\vec{v}_{j}\right).$

Stąd wynika, że błąd $|\mu _{1}-\eta _{k}|$ zależy od iloczynów skalarnych $\vec{v}_{j}^{{\, H}}\,\vec{v}_{1}$ oraz od stosunków

$\rho _{j}\,:=\,\frac{|\mu _{{j+1}}|}{|\mu _{1}|},\qquad 1\le j\le s-1.$

Dokładniej, niech

$\ell\,:=\,\min\,\left\{ 2\le j\le s:\;\vec{v}_{j}^{{\, H}}\,\vec{v}_{1}\ne 0\right\}$

albo $\ell=s+1$ jeśli powyższy zbiór jest pusty. Jeśli teraz $\ell=s+1$ albo mamy jednocześnie $2\le\ell\le s$ i $\rho _{{\ell}}<\rho _{1}^{2}$ to

$|\mu _{1}-\eta _{k}|\,\lessapprox\,|\mu _{1}|\cdot\rho _{1}^{{2k+1}}\cdot\frac{\|\vec{v}_{2}\| _{2}}{\|\vec{v}_{1}\| _{2}}$

i zbieżność jest liniowa z ilorazem $\rho _{1}^{2}$ . Jeśli zaś $2\le\ell\le s$ i $\rho _{{\ell}}\ge\rho _{1}^{2}$ to

$|\mu _{1}-\eta _{k}|\,\lessapprox\,|\mu _{1}|\cdot\rho _{{\ell}}^{k}\cdot\frac{\|\vec{v}_{2}\| _{2}}{\|\vec{v}_{1}\| _{2}}$

i zbieżność jest liniowa z ilorazem $\rho _{{\ell}}$ . W szczególności, jeśli $\vec{v}_{2}^{{\, H}}\,\vec{v}_{1}\ne 0$ to zbieżość jest z ilorazem $\rho _{1}$ .

Dla macierzy rzeczywistych i symetrycznych możemy stosunkowo łatwo pokazać dokładne nierówności na błąd metody potęgowej.

Twierdzenie 2.1

Załóżmy, że macierz $A=A^{T}\in\mathbb{R}^{{n,n}}$ . Niech $\gamma _{k}$ będzie kątem pomiędzy wektorem $k$ -tego przybliżenia $\vec{x}_{k}$ otrzymanego metodą potęgową i podprzestrzenią własną ${\mathcal{V}}_{1}$ odpowiadającą dominującej wartości własnej $\mu _{1}$ . Wtedy

$\tan\gamma _{k}\,\le\,\rho _{1}\cdot\tan\gamma _{{k-1}}$

oraz

$|\mu _{1}-\eta _{k}|\,\le\,\max _{{2\le j\le s}}|\mu _{1}-\mu _{j}|\cdot\sin^{2}\gamma _{k}.$

Dowód

Ponieważ macierz jest symetryczna, podprzestrzenie własne ${\mathcal{V}}_{j}$ zdefiniowane w (2.3) są parami ortogonalne. Dlatego

$\vec{x}_{k}\,=\,\frac{A^{k}\,\vec{x}_{0}}{\| A^{k}\,\vec{x}_{0}\| _{2}}\,=\,\frac{\vec{v}_{1}+\sum _{{j=2}}^{s}\left|\frac{\mu _{j}}{\mu _{1}}\right|^{k}\,\vec{v}_{j}}{\sqrt{\|\vec{v}_{1}\| _{2}^{2}+\sum _{{j=2}}^{s}\left|\frac{\mu _{j}}{\mu _{1}}\right|^{{2k}}\,\|\vec{v}_{j}\| _{2}^{2}}}.$

Tangens kąta $\gamma _{k}$ jest równy stosunkowi długości składowej wektora $\vec{x}_{k}$ w podprzestrzeni ${\mathcal{V}}_{2}\oplus\cdots\oplus{\mathcal{V}}_{s}$ do długości składowej tego wektora w podprzestrzeni ${\mathcal{V}}_{1}$ . Mamy więc

$\tan\gamma _{k}\,=\,\left(\sum _{{j=2}}^{s}\left|\frac{\mu _{j}}{\mu _{1}}\right|^{{2k}}\frac{\|\vec{v}_{j}\| _{2}^{2}}{\|\vec{v}_{1}\| _{2}^{2}}\right)^{{1/2}}\,\le\,\max _{{2\le j\le s}}\left|\frac{\mu _{j}}{\mu _{1}}\right|\cdot\tan\gamma _{{k-1}}.$

Aby pokazać pozostałą część twierdzenia, przedstawmy $\vec{x}_{k}$ w postaci $\vec{x}_{k}=\sum _{{j=1}}^{s}\vec{z}_{j}$ , gdzie $\vec{z}_{j}\in{\mathcal{V}_{j}}$ . Wtedy $\sum _{{j=1}}^{s}\|\vec{z}_{j}\| _{2}^{2}=\|\vec{x}_{k}\| _{2}^{2}=1$ oraz

$\vec{x}_{k}^{{\, T}}\, A\,\vec{x}_{k}\,=\,\sum _{{j=1}}^{s}\mu _{j}\,\|\vec{z}_{j}\| _{2}^{2}\,=\,\mu _{1}+\sum _{{j=2}}^{s}(\mu _{j}-\mu _{1})\,\|\vec{z}_{j}\| _{2}^{2},$

a stąd

$|\mu _{1}-\eta _{k}|\,\le\,\max _{{2\le j\le s}}|\mu _{1}-\mu _{j}|\cdot\sum _{{j=2}}^{s}\|\vec{z}_{j}\| _{2}^{2}\,=\,\max _{{2\le j\le s}}|\mu _{1}-\mu _{j}|\cdot\sin^{2}\gamma _{k}.$

∎

Przypomnijmy, że $\sin\gamma _{k}\le\tan\gamma _{k}$ przy czym mamy asymptotyczną równość gdy $\gamma _{k}\to 0$ . Ponieważ macierz $A$ jest symetryczna to mamy również $|\mu _{i}-\mu _{j}|\le 2\| A\| _{2}$ . Z twierdzenia 2.1 wynika więc, że

$\tan\gamma _{k}\,\le\,\rho _{1}^{k}\cdot\tan\gamma _{0},\qquad\mbox{oraz}\qquad|\mu _{1}-\eta _{k}|\,\le\, 2\rho _{1}^{{2k}}\cdot\| A\| _{2}\cdot\tan^{2}\gamma _{0}.$

Metoda potęgowa nie opłaca się gdy wymiar $n$ nie jest duży, albo macierz $A$ jest pełna. Inaczej jest, gdy $n$ jest ,,wielkie”, np. rzędu co najmniej kilkuset, a macierz $A$ jest rozrzedzona, tzn. ma jedynie proporcjonalnie do $n$ elementów niezerowych. Z taką sytuacją mamy do czynienia, gdy np. macierz jest pięciodiagonalna. Wtedy istotne dla metody mnożenie $A\,\vec{x}$ można wykonać kosztem proporcjonalnym do $n$ i poświęcając tyle samo pamięci (wyrazy zerowe można zignorować).

2.4. Odwrotna metoda potęgowa i deflacja

Odwrotna metoda potęgowa albo metoda Wielandta, polega na zastosowaniu (prostej) metody potęgowej do macierzy $(A-\sigma\, I)^{{-1}}$ . Dokładniej, startując z przybliżenia początkowego $\vec{x}_{0}$ o jednostkowej normie drugiej obliczamy kolejno dla $k=0,1,2,\ldots$

$\vec{y}_{k}\,:=\,(A-\sigma\, I)^{{-1}}\,\vec{x}_{k},\qquad\vec{x}_{{k+1}}\,:=\,\frac{\vec{y}_{k}}{\|\vec{y}_{k}\| _{2}},$

oraz $\eta _{k}:=\vec{x}^{{\, H}}_{k}\, A\,\vec{x}_{k}$ .

Od razu nasuwają się dwie uwagi. Po pierwsze, iteracja odwrotna ma sens tylko wtedy gdy parametr $\sigma$ jest tak dobrany, że macierz $A-\sigma\, I$ jest nieosobliwa, co jest równoważne warunkowi $\sigma\ne\mu _{i}$ dla wszystkich $i$ . Po drugie, w realizacji numerycznej, dla znalezienia wektora $\vec{y}_{k}$ nie odwracamy macierzy $(A-\sigma\, I)$ , ale rozwiązujemy układ równań $(A-\sigma\, I)\,\vec{y}=\vec{x}_{k}$ , co czyni metodę tak samo kosztowną co iteracja prosta.

Analiza iteracji odwrotnych przebiega podobnie do analizy iteracji prostych dla macierzy $(A-\sigma\, I)^{{-1}}$ . Nietrudno zauważyć, że macierz ta ma te same wektory własne co $A$ , a jej wartości własne wynoszą

$\mu _{i}^{{(\sigma)}}=\frac{1}{\mu _{i}-\sigma},\qquad 1\le i\le s.$

(To wynika bezpośrednio z równości $(A+\sigma\, I)\,\vec{v}_{i}=(\mu _{i}+\sigma)\,\vec{v}_{i}$ .) Dlatego iteracje odwrotne zbiegają do wektora własnego odpowiadającego wartości własnej $\mu _{{i^{*}}}$ takiej, że

$|\mu _{{i^{*}}}-\sigma|=\min _{{1\le i\le s}}|\mu _{i}-\sigma|,$

przy czym szybkość zbieżności zależy teraz od wielkości

$\rho _{j}^{{(\sigma)}}:=\frac{|\mu _{{i^{*}}}-\sigma|}{\min _{{i\ne i^{*}}}|\mu _{i}-\sigma|}$

(a nie od $\rho _{j}=|\mu _{{j+1}}|/|\mu _{1}|$ , jak w iteracji prostej).

Niewątpliwą zalety odwrotnej metody potęgowej jest to, że zbiega do wartości własnej najbliższej $\sigma$ , przy czym im $\sigma$ bliższe $\mu _{{i^{*}}}$ tym lepiej. Metoda jest więc szczególnie efektywna w przypadku gdy znamy dobre przybliżenia wartości własnych macierzy $A$ . Niestety, taka informacja nie zawsze jest dana wprost.

Pamiętając, że kolejne przybliżenia wartości własnej w metodzie potęgowej są obliczane przy pomocy ilorazów Rayleigha, dobrym pomysłem na przesunięcie wydaje się być wzięcie w $k$ -tym kroku iteracji odwrotnej

$\sigma=\eta _{k}=\vec{x}^{H}_{k}\, A\,\vec{x}_{k}.$

Rzeczywiście. Niech $\gamma _{k}$ będzie kątem pomiędzy $\vec{x}_{k}$ a podprzestrzenią własną $V_{{i^{*}}}$ . Zakładając, że $\eta _{k}$ jest dostatecznie blisko $\mu _{{i^{*}}}$ oraz przyjmując $\delta=\min _{{i\ne i^{*}}}|\mu _{{i^{*}}}-\mu _{i}|$ , na podstawie twierdzenia 2.1 mamy

	$\displaystyle\tan\gamma _{{k+1}}$	$\displaystyle\le$	$\displaystyle\frac{\|\mu _{{i^{}}}-\eta _{k}\|}{\min _{{i\ne i^{}}}\|\mu _{i}-\eta _{k}\|}\cdot\tan\gamma _{k}\,\le\,\frac{\|\mu _{{i^{}}}-\eta _{k}\|}{\delta-\|\mu _{{i^{}}}-\eta _{k}\|}\cdot\tan\gamma _{k}$
		$\displaystyle\le$	$\displaystyle\frac{2\,\\| A\\| _{2}\tan^{3}\gamma _{k}}{\delta-2\\| A\\| _{2}\tan^{2}\gamma _{k}}\,\approx\,\frac{2}{\delta}\cdot\\| A\\| _{2}\tan^{3}\gamma _{k}.$

Zamiast zbieżności kwadratowej, jak w prostej metodzie potęgowej, dostaliśmy (asymptotycznie!) zbieżność sześcienną.

Metoda potęgowa, prosta lub odwrotna, pozwala wyznaczyć tylko jedną wartość własną, powiedzmy $\mu _{1}$ , oraz odpowiadający jej wektor własny. Naturalne jest teraz pytanie o to jak postępować, aby znaleźć inne pary własne. Jednym ze sposobów jest zastosowanie deflacji.

Jeśli macierz jest hermitowska, $A=A^{H}$ , oraz znaleźliśmy wektor własny $\vec{v}_{1}$ odpowiadający wartości własnej $\mu _{1}$ , to możemy ponowić proces metody potęgowej ograniczając go do podprzestrzeni prostopadłej do $\vec{v}_{1}$ wymiaru $n-1$ . Osiągamy to startując z wektora $\vec{x}_{0}$ prostopadłego do $\vec{v}_{1}$ , tzn.

$\vec{x}_{0}:=\frac{\vec{v}}{\|\vec{v}\| _{2}},\qquad\vec{v}:=\vec{w}-\vec{v}_{1}\,(\vec{v}_{1}^{{\, H}}\,\vec{w}),$

gdzie $\vec{w}\ne\vec{0}$ jest wybrany losowo. Przy idealnej realizacji procesu konstruowany ciąg $\{\vec{x}_{k}\}$ należałby do podprzestrzeni prostopadłej do $\vec{v}_{1}$ . W obecności błędów zaokrągleń należałoby to wymusić poprzez reortogonalizację,

$\vec{y}_{k}:=\vec{y}_{k}-\vec{v}_{1}\,(\vec{v}_{1}^{{\, H}}\,\vec{y}_{k}),$

wykonywaną np. co kilka kroków. Po znalezieniu drugiej pary własnej, proces deflacji możemy kontynuować, ograniczając się do odpowiedniej podprzestrzeni własnej wymiaru $n-2$ .

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

Matematyka obliczeniowa II