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

Zagadnienia

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.

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.