MN13

Wektory i wartości własne

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Niech będzie dana rzeczywista kwadratowa macierz $\displaystyle A$ wymiaru $\displaystyle N$. Wektorem własnym $\displaystyle x\in C^{N}$ oraz odpowiadającą mu wartością własną $\displaystyle \lambda \in C$ nazwiemy taką parę, dla której

$\displaystyle Ax=\lambda x,$

przy czym $\displaystyle x\neq 0$.

Zadanie wyznaczania wartości własnych i wektorów własnych macierzy ma bardzo szerokie zastosowania w tak odległych od siebie dziedzinach, jak np. analiza odporności konstrukcji mechanicznych (wieżowce, mosty, wagony kolejowe) na wibracje, czy też rankingowanie stron internetowych w wyszukiwarce Google.

Przykład: Odporność budynku na trzęsienie ziemi

Rozważmy prosty układ mechaniczny opisujący, naturalnie w pewnym jedynie przybliżeniu, zachowanie się układu $\displaystyle N$ ciężkich płyt połączonych ze sobą relatywnie elatycznymi dźwigarami --- co może np. modelować konstrukcję wieżowca.

Wiadomo, że jeśli częstotliwości drgań własnych tego wieżowca będą bliskie częstotliwości siły wymuszającej (o niewielkiej amplitudzie), konstrukcja wpadnie w rezonans i w końcu rozpadnie się wskutek zbyt wielkich przemieszczeń. Wychylenia naszych płyt z położenia równowagi są opisywane układem pewnych równań różniczkowych. Teoria matematyczna takich równań różniczkowych pokazuje, że częstotliwości drgań własnych to nic innego jak wartości własne pewnej symetrycznej macierzy wymiaru $\displaystyle 2N$, która powstaje ze współczynników równania różniczkowego opisującego dynamikę tego układu.

Przykład: Macierz Google'a

Podstawowy algorytm rankingowania stron WWW w wyszukiwarce Google sprowadza się do znalezienia rzeczywistego wektora własnego $\displaystyle \pi$ pewnej silnie rozrzedzonej macierzy $\displaystyle A$ (gigantycznego rozmiaru, równego liczbie indeksowanych stron, czyli w chwili pisania tego tekstu około $\displaystyle 2.5\cdot 10^{10}$ stron), odpowiadającego wartości własnej równej 1:

$\displaystyle A\pi =\pi .$

Współrzędne wektora $\displaystyle \pi$ interpretuje się jako wartość rankingową kolejnych stron WWW. Aby wszystko miało sens, współrzędne wektora muszą być z przedziału [0,1]. Pewne twierdzenia matematyczne i subtelny dobór macierzy $\displaystyle A$ gwarantują, że taki wektor $\displaystyle \pi$ zawsze istnieje i jest jedyny! Co więcej, wartość 1 jest dominującą wartością własną $\displaystyle A$, a to z kolei ma ważne znaczenie dla tzw. metody potęgowej numerycznego wyznaczania takiego wektora.

Przykład: Wyznaczanie miejsc zerowych wielomianu

Jak wiadomo, wartości własne to miejsca zerowe wielomianu charakterystycznego macierzy $\displaystyle P(\lambda )=\det(A-\lambda I)$. Zachodzi także fakt odwrotny, to znaczy miejsca zerowe wielomianu są wartościami pewnej macierzy, np. miejsca zerowe wielomianu

$\displaystyle p(\lambda )=p_{1}\lambda ^{N}+\ldots +p_{N}\lambda +p_{N+1}$

są wartościami własnymi m.in. macierzy stowarzyszonej,

$\displaystyle A={\begin{pmatrix}-p_{2}/p_{1}&-p_{3}/p_{1}&\cdots &-p_{N+1}/p_{1}\\1&&&\\&1&&\\&&\ddots &\\&&&1\end{pmatrix}}$

Funkcja Octave'a compan(p) wyznacza macierz stowarzyszoną dla zadanego wielomianu o współczynnikach w wektorze $\displaystyle p=[p_{1},\ldots ,p_{N},p_{N+1}]^{T}$. Z tej macierzy korzysta następnie funkcja Octave'a roots, która właśnie w taki sposób wyznacza pierwiastki wielomianów: jako wartości własne macierzy stowarzyszonej.

W praktyce obliczeniowej spotyka się zazwyczaj kilka typów zagadnień:

• Wyznaczenie dominującej wartości własnej (to znaczy: największej co do modułu) i odpowiadającego jej wektora własnego (a może kilku wektorów, gdy wartość własna jest wielokrotna?)
• Wyznaczenie najmniejszej co do modułu wartości własnej i wektorów jej odpowiadających (zauważmy, że to jest np. zadanie wyznaczenia jądra macierzy osobliwej --- wtedy wiemy a priori, że szukana najmniejsza co do modułu wartość własna to zero)
• Wyznaczenie wartości własnej najbliższej zadanej liczbie (to jest właśnie odpowiedź na pytanie jak blisko częstości wymuszającej są częstości drgań własnych budynku)
• Wyznaczenie wszystkich wartości własnych (na przykład, w celu znalezienia wszystkich pierwiastków zadanego wielomianu)
• Wyznaczenie wszystkich wartości i wektorów własnych (tzw. pełne zagadnienie własne)

Jak domyślamy się, dla macierzy rozrzedzonych dużego wymiaru pełne zagadnienie własne jest zbyt kosztowne, gdyż najczęściej macierz wektorów własnych --- nawet dla macierzy rzadkiej --- jest gęsta.

Ponieważ w zastosowaniach bardzo często pojawiają się macierze rzeczywiste symetryczne (powyższe przykłady pokazują, że nie tylko!) szczegółową analizę metod numerycznych ograniczymy do tego przypadku, gdyż wtedy zachodzi

Twierdzenie O symetrycznym zadaniu własnym

Każda macierz rzeczywista symetryczna $\displaystyle A$ wymiaru $\displaystyle N$ ma rozkład

$\displaystyle A=Q\Lambda Q^{T},$

gdzie $\displaystyle Q\in R^{N\times N}$ jest ortogonalna (tzn. $\displaystyle Q^{T}Q=I$), a jej kolumnami są wektory własne $\displaystyle A$, natomiast $\displaystyle \Lambda \in R^{N}$ jest diagonalna z wartościami własnymi $\displaystyle A$ na diagonali:

$\displaystyle \Lambda ={\begin{pmatrix}\lambda _{1}&&\\&\ddots &\\&&\lambda _{N}\end{pmatrix}}.$

Lokalizacja wartości własnych

Jak okaże się za chwilę, czasem warto mieć ogólne rozeznanie o tym, gdzie z grubsza leżą wartości własne danej macierzy $\displaystyle A$. Przy tym mogą być nam pomocne dwa fakty:

Fakt

Dowolna wartość własna $\displaystyle \lambda \in C$ macierzy $\displaystyle A$ spełnia

$\displaystyle |\lambda |\leq ||A||,$

gdzie $\displaystyle ||A||$ jest dowolną normą macierzową indukowaną przez normę wektorową.

Dowód

Rzeczywiście, skoro istnieje wektor $\displaystyle x\neq 0$ taki, że $\displaystyle Ax=\lambda x$, to stąd $\displaystyle ||Ax||/||x||=|\lambda |$, więc fakt powyższy wynika już z definicji normy macierzy:

$\displaystyle ||A||=\max _{y\neq 0}{\frac {||Ay||}{||y||}}\geq ||Ax||/||x||.$

Drugie twierdzenie jest równie proste w dowodzie, ale daje trochę więcej informacji o lokalizacji widma.

Twierdzenie Gerszgorina, o lokalizacji widma macierzy

Wartości własne macierzy $\displaystyle A$ leżą w sumie (teoriomnogościowej) dysków $\displaystyle K_{i}$ na płaszczyźnie zespolonej,

$\displaystyle K_{i}=\{z\in C:|z-a_{ii}|\leq \sum _{j\neq i}|a_{ij}|\},\qquad i=1,\ldots N.$

Przykład: Koła Gerszgorina

Niech

$\displaystyle A={\begin{pmatrix}1.08930&1.38209&-1.00037&0.69355&2.32178\\0.14211&1.74696&1.68440&0.30664&1.26718\\-0.74620&2.02686&-0.68293&0.19684&0.35854\\0.83517&0.74987&1.71331&1.09765&-0.44321\\1.02132&-2.62155&0.79247&1.11408&0.48076\\\end{pmatrix}}$

Lokalizacja wartości własnych macierzy $\displaystyle A$ kołami Gerszgorina oraz zgrubna lokalizacja wewnątrz okręgu o promieniu równym $\displaystyle ||A||_{1}$. Dokładne wartości własne zaznaczone trójkącikami.

Wyznaczanie pojedynczej pary własnej

Jak wiemy z algebry, nawet gdy $\displaystyle A$ jest macierzą rzeczywistą, jej widmo może być zespolone! Analizując poniższe metody, będziemy zakładać, że poszukiwane wartości i wektory własne $\displaystyle A$ są rzeczywiste. Iterując na liczbach rzeczywistych nie mamy wszak szansy, by dotrzeć do liczb zespolonych!...

Metoda potęgowa

Przypuśćmy, że wartości własne macierzy $\displaystyle A\in R^{N\times N}$ spełniają

$\displaystyle |\lambda _{1}|>|\lambda _{2}|\geq \ldots \geq |\lambda _{N}|,$

(to znaczy, istnieje dokładnie jedna dominująca wartość własna macierzy $\displaystyle A$.

Załóżmy także, że istnieje baza złożona z wektorów własnych $\displaystyle q_{1},\ldots ,q_{N}$ tej macierzy (tak jest np. dla macierzy symetrycznej na mocy twierdzenia o własnościach symetrycznego zadania własnego).

Kierunek własny $\displaystyle q_{k}$ jakiejś macierzy $\displaystyle A$ ma taką własność, że poddany działaniu przekształcenia $\displaystyle A$ wydłuża się $\displaystyle \lambda _{k}$ razy, wobec tego dowolny wektor $\displaystyle x\in R^{N}$ poddany działaniu $\displaystyle A$ najbardziej wydłuży się w kierunku $\displaystyle q_{1}$. Iterując tę procedurę, powinniśmy dostawać w wyniku wektory, w których coraz bardziej dominuje kierunek $\displaystyle q_{1}$. Formalnie, niech

$\displaystyle x=\alpha _{1}q_{1}+\ldots +\alpha _{N}q_{N},$

wtedy

$\displaystyle Ax=A\left(\sum _{i}\alpha _{i}q_{i}\right)=\sum _{i}\alpha _{i}Aq_{i}=\sum _{i}\alpha _{i}\lambda _{i}q_{i}$

i w konsekwencji

$\displaystyle A^{k}x=\sum _{i}\alpha _{i}\lambda _{i}^{k}q_{i}=\lambda _{1}^{k}\left(\alpha _{1}q_{1}+\alpha _{2}\left({\frac {\lambda _{2}}{\lambda _{1}}}\right)^{k}q_{2}+\ldots +\alpha _{N}\left({\frac {\lambda _{N}}{\lambda _{1}}}\right)^{k}q_{N}\right).$

Założenia, że istnieje dokładnie jedna dominująca wartość własna, $\displaystyle \left|{\frac {\lambda _{N}}{\lambda _{1}}}\right|<1$, wynika, że wyrażenie w nawiasie dąży do $\displaystyle \alpha _{1}q_{1}$ i w konsekwencji wektory $\displaystyle x_{k}=A^{k}x$ dążą, gdy $\displaystyle k\rightarrow \infty$, do kierunku wektora własnego $\displaystyle q_{1}$, to znaczy wektora odpowiadającego dominującej wartości własnej $\displaystyle A$ (o ile tylko $\displaystyle \alpha _{1}\neq 0$).

Szybkość zbieżności metody potęgowej jest liniowa, o współczynniku zależnym od stosunku $\displaystyle |\lambda _{2}/\lambda _{1}|$. W patologicznym przypadku, gdy $\displaystyle |\lambda _{1}|\approx |\lambda _{2}|$, może więc okazać się, że metoda praktycznie nie jest zbieżna.

Zbieżność metody potęgowej w zależności od stosunku $\displaystyle |\lambda _{2}/\lambda _{1}|$. Na osi pionowej zaznaczono błąd przybliżenia dominującego wektora własnego.

W praktyce nie wyznaczamy wzorem $\displaystyle x_{k}=(A^{k})\cdot x$, lecz korzystamy z metody iteracyjnej

Algorytm Metoda potęgowa

---

<math>\displaystyle x_0</math> = dowolny wektor startowy; k = 0;
while( !stop )
{
	<math>\displaystyle y_k</math> = <math>\displaystyle Ax_{k-1}</math>;
	<math>\displaystyle x_k</math> = <math>\displaystyle y_k/||y_k||_\infty</math>;
	k++;	
} 

Warunek normowania ma m.in. na celu zapobieżenie powstawania nadmiaru i niedomiaru (gdy $\displaystyle |\lambda _{1}|<1$, to $\displaystyle ||A^{k}x||\rightarrow 0$, a gdy $\displaystyle |\lambda _{1}|>1$, to $\displaystyle ||A^{k}x||\rightarrow \infty$). Przy okazji, zauważ, że $\displaystyle ||y_{k}||_{\infty }\rightarrow |\lambda _{1}|$, a więc mamy także sposób na wyznaczenie przybliżenia dominującej wartości własnej.

Zazwyczaj jako warunek stopu wybiera się kryterium małej poprawki, $\displaystyle ||x_{k}-x_{k-1}||\leq \epsilon$ lub warunek małego residuum, $\displaystyle ||Ax_{k}-\lambda _{1,k}x_{k}||\leq \epsilon$, gdzie $\displaystyle \lambda _{1,k}$ jest przybliżeniem $\displaystyle \lambda _{1}$ dostępnym na $\displaystyle k$-tej iteracji.

Gdy $\displaystyle x_{0}$ nie ma składowej w kierunku dominującego wektora własnego, w teorii nie powinniśmy obserwować zbieżności. Jednak w praktyce, dzięki błędom zaokrągleń, w ciągu iteracji taka składowa się pojawia i doprowadza do zbieżności. Na rysunku przedstawiono zbieżność metody potęgowej w zależności od stosunku $\displaystyle |\lambda _{2}/\lambda _{1}|$. Na osi pionowej zaznaczono błąd przybliżenia dominującego wektora własnego. Oczywiście, gdy $\displaystyle \lambda _{2}=\lambda _{1}$ na zbieżność praktycznie nie ma szans.

Metoda potęgowa doskonale sprawdza się, gdy macierz $\displaystyle A$ jest macierzą rozrzedzoną --- np. w przypadku macierzy Google'a.

Odwrotna metoda potęgowa

Zauważmy, że dla dowolnej macierzy kwadratowej $\displaystyle A$ o wartościach własnych $\displaystyle \lambda _{k}$ i odpowiadających im wektorach własnych $\displaystyle q_{k}$, mamy:

• Macierz $\displaystyle A-\sigma I$ ma wartości własne $\displaystyle \lambda _{k}-\sigma$ oraz wektory własne $\displaystyle q_{k}$,
• Jeśli dodatkowo $\displaystyle A$ jest nieosobliwa, to macierz $\displaystyle A^{-1}$ ma wartości własne $\displaystyle 1/\lambda _{k}$ oraz wektory własne $\displaystyle q_{k}$

Z połączenia tych dwóch własności wynika, że

Stwierdzenie O transformacji widma macierzy

Macierz $\displaystyle (A-\sigma I)^{-1}$ (o ile istnieje), to ma wartości własne równe $\displaystyle {\frac {1}{\lambda _{k}-\sigma }}$ i wektory własne identyczne z $\displaystyle A$.

Skoro tak, to jeśli najbliższą $\displaystyle \sigma$ wartością własną $\displaystyle A$ jest $\displaystyle \lambda _{j}$, wówczas metoda potęgowa zastosowana do macierzy $\displaystyle (A-\sigma I)^{-1}$ zbiegnie do $\displaystyle q_{j}$. To prowadzi do następującego algorytmu, odwrotnej metody potęgowej:

Algorytm Odwrotna metoda potęgowa

---

<math>\displaystyle x_0</math> = dowolny wektor startowy; k = 0;
while( !stop )
{
	Rozwiąż układ równań <math>\displaystyle (A-\sigma I)y_k = x_{k-1}</math>;
	<math>\displaystyle x_k</math> = <math>\displaystyle y_k/||y_k||_\infty</math>;
	k++;	
} 

Metoda Rayleigh (RQI)

Z własności metody potęgowej wynika, że metoda odwrotna potęgowa jest zbieżna tym szybciej, im bliżej $\displaystyle \lambda _{j}$ jest przesunięcie $\displaystyle \sigma$ (w stosunku do pozostałych wartości własnych). Dlatego dobrze byłoby --- dla zwiększenia szybkości zbieżności iteracji --- poprawiać wartość przesunięcia $\displaystyle \sigma$, korzystając z dotychczas wyznaczonego wektora $\displaystyle x_{k}\approx q_{j}$ i ilorazu Rayleigh:

$\displaystyle \lambda _{j}={\frac {q_{j}^{T}Aq_{j}}{q_{j}^{T}q_{j}}}\approx {\frac {x_{k}^{T}Ax_{k}}{x_{k}^{T}x_{k}}}$

Stąd nazwa metody, w skrócie RQI (Rayleigh Quotient Iteration).

Algorytm Metoda RQI

---

<math>\displaystyle x_0</math> = dowolny wektor startowy; <math>\displaystyle \sigma_0</math> = przybliżenie <math>\displaystyle \lambda_j</math>; k = 0;
while( !stop )
{
	Rozwiąż układ równań <math>\displaystyle (A-\sigma_k I)y_k = x_{k-1}</math>;
	<math>\displaystyle x_k</math> = <math>\displaystyle y_k/||y_k||_2</math>;
	<math>\displaystyle \sigma_{k+1}</math> = <math>\displaystyle x_k^TAx_k</math>;
	k++;	
} 

(wybierając normowanie wektora $\displaystyle x$ w normie euklidesowej upraszczamy co nieco algorytm).

Wielką zaletą metody RQI jest jej szybkość zbieżności: kwadratowa gdy wartość własna jest pojedyncza, a nawet sześcienna w przypadku macierzy symetrycznej.

Wadą metody RQI jest to, że na każdym jej kroku należy rozwiązywać układ równań z inną macierzą.

Uwaga: Gdy złe uwarunkowanie i skończona precyzja arytmetyki pomagają...

Przez pewien czas numerycy odnosili się do tej metody z rezerwą, twierdząc, i słusznie, że im lepszym przybliżeniem $\displaystyle q_{j}$ będzie $\displaystyle \sigma _{k}$, tym bardziej rośnie uwarunkowanie $\displaystyle A-\sigma _{k}I$, a tym samym błąd numerycznego rozwiązywania układu z tą macierzą będzie coraz większy i metoda będzie tracić stabilność. Tymczasem okazuje się, że --- choć rzeczywiście rozwiązanie układu jest obarczone wielkim błędem --- to wektor błędu ma kierunek praktycznie zgodny z kierunkiem poszukiwanego wektora $\displaystyle q_{j}$, a tym samym złe uwarunkowanie macierzy i skończona precyzja arytmetyki pomagają w zbieżności metody!

Metody rozwiązywania pełnego zadania własnego

Najszybszą obecnie znaną metodą rozwiązywania pełnego zadania własnego (to znaczy znajdowania wszystkich wartości i wektorów własnych) macierzy symetrycznej jest metoda dziel i rządź.

Dla macierzy niesymetrycznych najbardziej dopracowanym i przetestowanym, a więc godnym zaufania algorytmem, jest metoda QR z przesunięciami (wykorzystująca, jak łatwo się domyślić, rozkład QR macierzy). Metoda QR przewyższa także metodę "dziel i rządź" w przypadku symetrycznym, gdy wymiar macierzy jest mały (mniej więcej $\displaystyle N\leq 25$).

Obie metody są oczywiście metodami iteracyjnymi, jednak przyjęło się nazywać je metodami bezpośrednimi, gdyż praktycznie zawsze potrzebują z góry ograniczonej liczby iteracji do tego, by zbiec do wyniku o (niemal) maksymalnej rozsądnej dokładności.

Dla efektywności obu metod kluczowy jest preprocessing macierzy, pozwalający niezbyt wygórowanym kosztem $\displaystyle O(N^{3})$ operacji sprowadzić przez ortogonalne podobieństwo zadanie z macierzą gęstą $\displaystyle A$ (w przypadku niesymetrycznym) do zadania z macierzą Hessenberga, czyli macierzą, której element $\displaystyle (i,j)$ jest zerowy gdy tylko $\displaystyle i-j>1$:

$\displaystyle {\begin{pmatrix}*&*&*&*&\cdots &*\\*&*&*&*&\cdots &*\\&*&*&*&\cdots &*\\&&\ddots &\ddots &\ldots &\vdots \\&&&\ddots &\ddots &*\\&&&&*&*\end{pmatrix}}$

bądź wręcz trójdiagonalną, gdy $\displaystyle A$ była symetryczna.

Każdą macierz kwadratową $\displaystyle A$ da się sprowadzić do postaci Hessenberga sekwencją przekształceń postaci

$\displaystyle A:=Q_{k}AQ_{k}^{T},$

gdzie $\displaystyle Q_{k}$ jest pewnym przekształceniem Householdera. Rzeczywiście, niech

$\displaystyle A={\begin{pmatrix}d_{1}&*&*&*&\cdots &*\\a_{1}&*&*&*&\cdots &*\\a_{2}&*&*&*&\cdots &*\\\vdots &*&\ddots &\ddots &\ldots &\vdots \\\vdots &*&*&\ddots &\ddots &*\\a_{N-1}&*&*&*&*&*\end{pmatrix}}$

i oznaczmy $\displaystyle a=[a_{1},\ldots ,a_{N-1}]^{T}$. Możemy wziąć na początek przekształcenie Householdera $\displaystyle {\widetilde {Q}}_{1}$ takie, że $\displaystyle {\widetilde {Q}}_{1}a=c\cdot e_{1}$, gdzie $\displaystyle e_{1}=[1,0,\ldots ,0]^{T}$. Wtedy

$\displaystyle {\begin{pmatrix}1&\\&{\widetilde {Q}}_{1}\end{pmatrix}}\cdot A={\begin{pmatrix}d_{1}&*&*&*&\cdots &*\\c&*&*&*&\cdots &*\\&*&*&*&\cdots &*\\&*&\ddots &\ddots &\ldots &\vdots \\&*&*&\ddots &\ddots &*\\&*&*&*&*&*\end{pmatrix}}$

To samo przekształcenie przyłożone z prawej strony zachowa pierwszą kolumnę i w efekcie nie zmieni struktury macierzy:

$\displaystyle {\begin{pmatrix}1&\\&{\widetilde {Q}}_{1}\end{pmatrix}}\cdot A\cdot {\begin{pmatrix}1&\\&{\widetilde {Q}}_{1}\end{pmatrix}}={\begin{pmatrix}d_{1}&*&*&*&\cdots &*\\c&*&*&*&\cdots &*\\&*&*&*&\cdots &*\\&*&\ddots &\ddots &\ldots &\vdots \\&*&*&\ddots &\ddots &*\\&*&*&*&*&*\end{pmatrix}}.$

Dalej stosujemy tę samą metodę do podmacierzy wymiaru $\displaystyle N-1$, itd. aż dochodzimy do macierzy Hessenberga.

Gdy wyjściowa macierz $\displaystyle A$ jest symetryczna, to z definicji, macierz wynikowa $\displaystyle {\begin{pmatrix}I&\\&{\widetilde {Q}}_{N-2}\end{pmatrix}}\cdots {\begin{pmatrix}1&\\&{\widetilde {Q}}_{1}\end{pmatrix}}A{\begin{pmatrix}1&\\&{\widetilde {Q}}_{1}\end{pmatrix}}\cdots {\begin{pmatrix}I&\\&{\widetilde {Q}}_{N-2}\end{pmatrix}}$ też jest symetryczna i jednocześnie Hessenberga --- a więc musi być trójdiagonalna! Ponadto, macierz wynikowa będzie miała te same wartości własne co $\displaystyle A$; wektory własne macierzy $\displaystyle A$ także można łatwo (jak?) odzyskać z wektorów własnych macierzy wynikowej.

Metoda "dziel i rządź" dla macierzy symetrycznej

Jest to obecnie najefektywniejsza metoda rozwiązywania zagadnienia własnego macierzy symetrycznej wymiaru powyżej kilkudziesięciu. Omówimy w zarysie jej najprostszy wariant (obarczony pewnymi wadami, usuniętymi w wersji bibliotecznej --- DSYEVD w LAPACKu).

Startując z symetrycznej macierzy $\displaystyle A$ już w postaci trójdiagonalnej, łatwo widzieć, że "prawie" rozpada się ona na dwie mniejsze macierze trójdiagonalne: dokładniej,

$\displaystyle {\begin{pmatrix}a_{1}&b_{1}&&\\b_{1}&a_{2}&\ddots &\\&\ddots &\ddots &b_{N-1}\\&&b_{N-1}&a_{N}\end{pmatrix}}={\begin{pmatrix}T_{1}&\\&T_{2}\end{pmatrix}}+b_{m}uu^{T},$

gdzie $\displaystyle T_{1}={\begin{pmatrix}a_{1}&b_{1}&&\\b_{1}&a_{2}&\ddots &\\&\ddots &\ddots &b_{m-1}\\&&b_{m-1}&a_{m}-b_{m}\end{pmatrix}}$, $\displaystyle T_{2}={\begin{pmatrix}a_{m+1}-b_{m}&b_{m+1}&&\\b_{m+1}&a_{m+2}&\ddots &\\&\ddots &\ddots &b_{N-1}\\&&b_{N-1}&a_{N}\end{pmatrix}}$ są --- tak jak $\displaystyle A$ --- macierzami trójdiagonalnymi i symetrycznymi (jako podmacierze $\displaystyle A$), tylko o połowę mniejszego wymiaru, gdy $\displaystyle m\approx N/2.$ Natomiast $\displaystyle u=e_{m}+e_{m+1}$, więc macierz $\displaystyle b_{m}uu^{T}$ ma tylko cztery niezerowe elementy, każdy równy $\displaystyle b_{m}$.

Zgodnie ze swoją nazwą, metoda "dziel i rządź" sprowadza zadanie znajdowania par własnych macierzy wymiaru $\displaystyle N$ do dwóch takich zadań dla macierzy dwa razy mniejszych. Te z kolei można potraktować w taki sam sposób i iteracyjnie zmniejszyć wymiar macierzy do tak małego (około 25), by opłacało się zastosować metodę QR (teoretycznie, można byłoby oczywiście doprowadzić podział do momentu, gdy macierze trójdiagonalne są rozmiaru $\displaystyle 1\times 1$ --- dla których rozwiązanie zadania włanego jest trywialne --- ale taki algorytm byłby bardziej kosztowny od wariantu z udziałem QR).

Rzeczywiście, przypuśćmy, że dla obu macierzy trójdiagonalnych $\displaystyle T_{1},T_{2}$ umiemy rozwiązać zadanie własne tak, że znamy macierze: $\displaystyle Q_{i}$ --- ortogonalną oraz $\displaystyle D_{i}$ --- diagonalną, takie, że

$\displaystyle Q_{i}^{T}T_{i}Q_{i}=D_{i}\qquad i=1,2.$

Wtedy łatwo widzieć, że dla łatwo wyznaczalnego wektora $\displaystyle v$,

$\displaystyle {\begin{pmatrix}Q_{1}^{T}&\\&Q_{2}^{T}\end{pmatrix}}\left({\begin{pmatrix}T_{1}&\\&T_{2}\end{pmatrix}}+b_{m}uu^{T}\right){\begin{pmatrix}Q_{1}&\\&Q_{2}\end{pmatrix}}={\begin{pmatrix}D_{1}&\\&D_{2}\end{pmatrix}}+b_{m}vv^{T}.$

W ten sposób zadanie własne dla oryginalnej macierzy $\displaystyle T$ wymiaru $\displaystyle N$ jest równoważne zadaniu własnemu macierzy diagonalnej zaburzonej o macierz rzędu 1.

Na szczęście łatwo pokazać, że jeśli $\displaystyle \lambda$ nie jest wartością własną macierzy diagonalnej $\displaystyle D={\begin{pmatrix}D_{1}&\\&D_{2}\end{pmatrix}}$, to wartości własne $\displaystyle \lambda$ macierzy $\displaystyle D+b_{m}vv^{T}$

spełniają równanie

$\displaystyle f(\lambda )\equiv 1+b_{m}\sum _{j=1}^{N}{\frac {v_{j}^{2}}{d_{j}-\lambda }}=0,$

gdzie $\displaystyle d_{j}$ są elementami na diagonali macierzy $\displaystyle D$.

Wykres $\displaystyle f(\lambda )$ dla macierzy jednowymiarowego laplasjanu rozmiaru 10. Zwróć uwagę na asymptoty pionowe tej funkcji oraz jej przedziałową monotoniczność.

W typowym przypadku $\displaystyle f$ będzie miała dokładnie $\displaystyle N$ pojedynczych miejsc zerowych i wykres zachęcający do stosowania do niej metody Newtona. Okazuje się, że ogólny przypadek nie jest istotnie trudniejszy, choć wymaga ważnych modyfikacji, zarówno w celu szybszego rozwiązywania powyższego równania nieliniowego, jak i w celu zapewnienia lepszej stabilności algorytmu.

Ostateczny koszt wyznaczenia wszystkich wektorów i wartości własnych jest rzędu $\displaystyle O(N^{3})$ z małą stałą.

Metoda QR

Dla zadania własnego z macierzą niesymetryczną najczęściej stosuje się metodę QR.

Jakkolwiek ostateczna wersja metody QR działa dla macierzy niesymetrycznych, wygodnie będzie nam założyć dla przejrzystości ekspozycji, że macierz jest symetryczna i w konsekwencji ma rzeczywiste widmo.

W najprostszym wariancie (bez przesunięć), algorytm QR ma postać:

Algorytm Metoda QR

---

<math>\displaystyle A_1 = A</math>;
for k = 1, 2, ...
{
	wykonaj rozkład <math>\displaystyle A_k = Q_{k}R_{k}</math>;
	<math>\displaystyle A_{k+1} = R_{k}\cdot Q_{k}</math>;
}

Można sprawdzić, że $\displaystyle A,A_{1},A_{2},\ldots$ mają te same wartości własne, bo $\displaystyle A_{k+1}=Q_{k+1}^{T}A_{k}Q_{k+1}$. Co więcej, powyższy algorytm (gdy $\displaystyle A$ jest nieosobliwa) w zasadzie jest równoważny hipotetycznemu algorytmowi iteracji prostej zastosowanemu nie do pojedynczego wektora, ale do $\displaystyle N$ wektorów naraz:

Algorytm Metoda potęgowa na przestrzeni

---

<math>\displaystyle V_1 = I</math>;
for k = 1, 2, ...
{
	<math>\displaystyle W_{k+1} = A\cdot V_k</math>;
	wyznacz rozkład QR <math>\displaystyle W_{k+1} = V_{k+1} R_{k+1}</math>, gdzie <math>\displaystyle V_{k+1}</math> jest ortogonalna;
}

Drugi krok w istocie ortogonalizuje kolumny $\displaystyle W_{k+1}$. Gdyby nie ortogonalizować zestawu wektorów $\displaystyle W_{k+1}$, oczywiście dostalibyśmy w efekcie zbieżność wszystkich kolumn macierzy do tego samego wektora --- odpowiadającego dominującej wartości własnej $\displaystyle A$. Zapewniając sobie ortogonalność $\displaystyle V_{k+1}$, możemy liczyć na to, że kolejne kolumny macierzy $\displaystyle V_{k}$ będą dążyć do wektorów własnych odpowiadających kolejnym wartościom własnym $\displaystyle A$ (przy stosownych założeniach o $\displaystyle A$, m.in. że wszystkie wartości własne $\displaystyle A$ spełniają $\displaystyle |\lambda _{i}|\neq |\lambda _{j}|$ dla $\displaystyle i\neq j$). Jeśli założyć dla uproszczenia, że oba używane rozkłady QR mają jednoznacznie określone czynniki rozkładu (na przykład, wymuszając, by diagonalne elementy macierzy $\displaystyle R$ były dodatnie) mamy zależności $\displaystyle V_{k+1}=Q_{1}\cdots Q_{k}$ oraz $\displaystyle A_{k+1}=V_{k+1}^{T}AV_{k+1}$.

Tak więc, w sprzyjających warunkach, metoda QR, jako równoważna metodzie potęgowej na przestrzeni, będzie zbieżna: $\displaystyle A_{k}\rightarrow A_{\infty }$, gdzie $\displaystyle A_{\infty }$ jest macierzą trójkątną (bo wektory własne odpowiadające różnym wartościom własnym są ortogonalne. Tym samym, wartościami własnymi $\displaystyle A_{\infty }$ (a więc także $\displaystyle A$) będą liczby na diagonali $\displaystyle A_{\infty }$.

Twierdzenie O zbieżności metody QR

Niech wartości własne $\displaystyle A\in R^{N\times N}$ spełniają $\displaystyle |\lambda _{1}|,\ldots ,|\lambda _{N}|>0$ oraz macierz $\displaystyle T=[x_{1},\ldots ,x_{N}]$ o kolumnach $\displaystyle x_{i}$ złożonych z kolejnych wektorów własnych $\displaystyle A$ ma taką własność, że $\displaystyle T^{-1}$ ma rozkład LU, $\displaystyle T^{-1}=LU$.

Wtedy w metodzie QR ciąg macierzy $\displaystyle Q_{k}$ jest zbieżny do macierzy diagonalnej, a ciąg $\displaystyle A_{k}$ ma podciąg zbieżny do macierzy trójkątnej, której elementy diagonalne $\displaystyle u_{ii}$ są równe $\displaystyle \lambda _{i}$ dla $\displaystyle i=1,\ldots ,N$.

Powyższa wersja algorytmu QR jest mało praktyczna, m.in. jest zbieżna wolno i przy poważnych ograniczaniach na $\displaystyle A$. Sprytna modyfikacja algorytmu wyjściowego daje w wyniku tzw. metodę QR z przesunięciami, która jest praktycznie niezawodna dla dowolnej macierzy.

Algorytm Metoda QR z przesunięciami

---

<math>\displaystyle A_1 = A</math>;
for k = 1, 2, ...
{
	wybierz sprytnie przesunięcie <math>\displaystyle \sigma_k</math>; 
	wykonaj rozkład <math>\displaystyle A_k - \sigma_kI = Q_{k}R_{k}</math>;
	<math>\displaystyle A_{k+1} = R_{k}\cdot Q_{k} + \sigma_kI</math>;
}

Koszt wyznaczenia wszystkich wektorów i wartości własnych jest rzędu $\displaystyle O(N^{3})$ ze stałą równą około 30. Omówienie sposobów wyboru "sprytnego przesunięcia" wykracza niestety poza ramy wykładu.

Metoda Jacobiego dla macierzy symetrycznej

Na zakończenie wspomnijmy o (bardzo starej) metodzie Jacobiego, która działa na oryginalnej macierzy symetrycznej (bez konieczności uprzedniego sprowadzenia do postaci trójdiagonalnej).

Nie będziemy szczegółowo jej tu omawiać z braku miejsca, wymienimy tylko jej dwie najważniejsze cechy odróżniające ją od metod omawianych wcześniej:

• jest znacznie wolniejsza
• (niekiedy) wyznacza dokładniejsze wartości i wektory własne niż inne metody

Dodatkową zaletą metody Jacobiego jest to, że łatwo ją zrównoleglić (w czym jest podobna do metody dziel i rządź).

Pomysł metody Jacobiego jest stosunkowo prosty: należy sekwencją przekształceń ortogonalnych $\displaystyle J_{0},J_{1},\ldots$ sprowadzić wyjściową macierz symetryczną $\displaystyle A$ do postaci (prawie) diagonalnej:

$\displaystyle J_{k}^{T}J_{k-1}^{T}\cdots J_{0}^{T}AJ_{0}\cdots J_{k-1}J_{k}\approx {\begin{pmatrix}\lambda _{1}&&&\\&\lambda _{2}&&\\&&\ddots &\\&&&\lambda _{N}\end{pmatrix}}.$

Tak więc iteracja Jacobiego ma postać:

Algorytm Metoda Jacobiego

---

for k = 0, 1, ...
	wybierz <math>\displaystyle J_k</math>;
	<math>\displaystyle A</math> = <math>\displaystyle J_k^T A J_k</math>;

Macierze $\displaystyle J_{k}$ wybieramy jako tzw. obroty Jacobiego, tzn. obroty Givensa dobrane tak, by w danym kroku iteracji wyzerować kolejną parę pozadiagonalnych elementów macierzy. W klasycznej wersji, zerowaniu podlega pozadiagonalna para o największym module --- w ten sposób najbardziej zredukujemy miarę niediagonalności macierzy wyrażoną jako suma kwadratów elementów pozadiagonalnych:

$\displaystyle \omega =\sum _{j<i}a_{ij}^{2}.$

Twierdzenie O zbieżności metody Jacobiego

Klasyczna metoda Jacobiego jest zbieżna co najmniej liniowo, tzn.

$\displaystyle \omega _{k+1}\leq {\sqrt {1-{\frac {2}{N(N-1)}}}}\omega _{k},$

gdzie $\displaystyle \omega _{k}$ oznacza miarę niediagonalności na $\displaystyle k$-tym kroku iteracji.

Można pokazać, że asymptotycznie (tzn. dostatecznie blisko granicy) zbieżność metody Jacobiego jest nawet kwadratowa.

Uwarunkowanie

Twierdzenie Bauera-Fike'a, o uwarunkowaniu wartości własnych

Niech $\displaystyle A\in R^{N\times N}$ będzie diagonalizowalna, to znaczy dla pewnej macierzy $\displaystyle X$ zachodzi

$\displaystyle X^{-1}AX={\begin{pmatrix}\lambda _{1}&&\\&\ddots &\\&&\lambda _{N}\end{pmatrix}},$

a więc (gdyż macierz po prawej stronie jest podobna do $\displaystyle A$) $\displaystyle \lambda _{i}\in C$, $\displaystyle i=1,\ldots ,N$ są wartościami własnymi $\displaystyle A$. Rozważmy macierz zaburzoną $\displaystyle {\widetilde {A}}$ i jakąś jej wartość własną $\displaystyle {\widetilde {\lambda }}$. Wtedy istnieje wartość własna $\displaystyle \lambda _{j}$ macierzy $\displaystyle A$ taka, że

$\displaystyle |\lambda _{j}-{\widetilde {\lambda }}|\leq {\mbox{cond}}_{2}(X)||A-{\widetilde {A}}||_{2}.$

Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia $\displaystyle X$ jest ortogonalna, $\displaystyle X^{-1}=X^{T}$, to mamy $\displaystyle {\mbox{cond}}_{2}(X)=1$ i w konsekwencji zachodzi

Wniosek Wartości własne macierzy symetrycznej są doskonale uwarunkowane

Przy oznaczeniach jak twierdzeniu Bauera-Fike'a, jeśli dodatkowo założymy, że macierz $\displaystyle A$ jest rzeczywista i symetryczna, to

$\displaystyle \min _{j=1,\ldots ,N}|\lambda _{j}-{\widetilde {\lambda }}|\leq ||A-{\widetilde {A}}||_{2}.$

Z drugiej strony, dla macierzy niediagonalizowalnych, uwarunkowanie wartości własnych może być dowolnie duże, co ilustruje poniższy

Przykład

$\displaystyle A_{\epsilon }={\begin{pmatrix}a&1\\\epsilon &a\end{pmatrix}}$

Weźmy dla uproszczenia $\displaystyle a=0$. Wartości własne $\displaystyle A_{\epsilon }$ to zera wielomianu $\displaystyle p_{\epsilon }(\lambda )=\lambda ^{2}-\epsilon$, zatem $\displaystyle \lambda _{\epsilon }=\pm {\sqrt {\epsilon }}$ i w konsekwencji

$\displaystyle |\lambda _{\epsilon }-\lambda _{0}|/||A_{\epsilon }-A_{0}||={\sqrt {\epsilon }}/\epsilon \rightarrow \infty ,$

gdy $\displaystyle \epsilon \rightarrow 0^{+}$, a więc uwarunkowanie takiego zadania jest nieskończone: dowolnie mała zmiana macierzy powoduje zaburzenie wartości własnych niewspółmiernie wielkie wobec zaburzenia danych. Dodatkowo, wartości własne i wektory własne macierzy $\displaystyle A$ dla ujemnego parametru $\displaystyle \epsilon$ są zespolone!

Zachowanie się wartości własnych macierzy $\displaystyle A$ (z parametrem $\displaystyle a=1$) w otoczeniu $\displaystyle \epsilon =0$

Bardziej spektakularny przykład pochodzi od Wilkinsona:

Przykład: Perfidny wielomian Wilkinsona

Niech

$\displaystyle p(\lambda )=(\lambda -1)(\lambda -2)\cdots (\lambda -20).$

Zmiana współczynnika przy $\displaystyle \lambda ^{19}$ o $\displaystyle 10^{-7}$ skutkuje przesunięciem niektórych miejsc zerowych nawet o kilka jednostek na płaszczyźnie zespolonej! Poniżej pokazujemy to na numerycznym przykładzie, gdzie prócz wyżej wymienionego zaburzenia mamy dodatkowo do czynienia z zaburzeniami powstałymi wskutek wyznaczenia współczynników wielomianu w arytmetyce zmiennoprzecinkowej.

Zera oryginalnego i lekko zaburzonego perfidnego wielomianu Wilkinsona.

Jak widzimy, zera bardzo mało zaburzonego wielomianu mogą stać się wyraźnie nie-rzeczywiste!

Jeśli chodzi o wektory własne, ich wrażliwość na zaburzenia macierzy jest bardziej skomplikowana.

Biblioteki

LAPACK zawiera w sobie kolekcję doskonałych narzędzi do rozwiązywania różnych wariantów zadania własnego, m.in. DGEEV dla macierzy niesymetrycznych oraz DSYEV dla macierzy symetrycznych rozwiązują pełne zagadnienie własne, wyznaczając wszystkie wartości własne i wektory własne. Dla macierzy symetrycznych mamy jeszcze m.in. funkcje DSYEVX (dla wybranych wartości własnych) i DSYEVD (z algorytmem "dziel i rządź")

Fortranowska biblioteka ARPACK rozwiązuje zadanie własne dla macierzy rozrzedzonych, znajdując kilka wybranych (np. największych co do modułu) wartości i wektorów własnych.

Funkcja eig w Octave i MATLABie wyznacza wszystkie wartości własne (i opcjonalnie wektory własne) zadaniej gęstej macierzy --- oczywiście korzystając z LAPACKa.

octave:1> A = [0 1; 1e-5, 0] A = 0.00000 1.00000 0.00001 0.00000 octave:2> eig(A) ans = 0.0031623 -0.0031623 octave:3> [V, L] = eig(A) V = 0.9999950 -0.9999950 0.0031623 0.0031623 L = 0.0031623 0.0000000 0.0000000 -0.0031623

Jak dotąd, tylko MATLAB potrafi skorzystać z ARPACKa dla wyznaczenia fragmentów widma macierzy rzadkiej, za pomocą funkcji eigs.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 5.1, 5.2 i 5.5 w

• D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Część wykładu oparto na materiałach zawartych w bardzo ciekawym podręczniku

• J. Demmel, Numerical linear algebra, SIAM, 1997.

Od dziesięcioleci, wspaniałym przeżyciem jest lektura książki ojca nowoczesnej analizy numerycznej,

• J. H. Wilkinson, The algebraic eigenvalue problem, Clarendon Press, 1965,

a także

• B. Parlett, The symmetric eigenvalue problem, Prentice-Hall, 1980.