MN13

Wektory i wartości własne

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

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

$Ax=\lambda x$,

przy czym $x\ne 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 $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 $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 $\pi$ pewnej silnie rozrzedzonej macierzy $A$ (gigantycznego rozmiaru, równego liczbie indeksowanych stron, czyli w chwili pisania tego tekstu około $2.5\cdot 10^{10}$ stron), odpowiadającego wartości własnej równej 1:

$A\pi =\pi$

Współrzędne wektora $\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 $A$ gwarantują, że taki wektor $\pi$ zawsze istnieje i jest jedyny! Co więcej, wartość 1 jest dominującą wartością własną $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 $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

$p(\lambda )=p_1{\lambda }^N+\dots +p_N\lambda +p_{N+1}$

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

$A=\left(\begin{array}{cccc}-p_2/p_1& -p_3/p_1& \cdots & -p_{N+1}/p_1\\ 1& & & \\ & 1& & \\ & & \ddots & \\ & & & 1\end{array}\right)$

Funkcja Octave'a compan(p) wyznacza macierz stowarzyszoną dla zadanego wielomianu o współczynnikach w wektorze $p=[p_1,\dots ,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 $A$ wymiaru $N$ ma rozkład

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

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

$\mathrm{\Lambda }=\left(\begin{array}{ccc}{\lambda }_1& & \\ & \ddots & \\ & & {\lambda }_N\end{array}\right)$

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 $A$. Przy tym mogą być nam pomocne dwa fakty:

Fakt

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

$|\lambda |\le ||A||$,

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

Dowód

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

$||A||=\underset{y\ne 0}{\max }\frac{||Ay||}{||y||}\ge ||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 $A$ leżą w sumie (teoriomnogościowej) dysków $K_i$ na płaszczyźnie zespolonej,

$K_i=\{z\in C:|z-a_{ii}|\le \sum _{j\ne i}|a_{ij}|\},i=1,\dots N$

Przykład: Koła Gerszgorina

Niech

$A=\left(\begin{array}{ccccc}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{array}\right)$

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

Wyznaczanie pojedynczej pary własnej

Jak wiemy z algebry, nawet gdy $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 $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 $A\in R^{N\times N}$ spełniają

$|{\lambda }_1|>|{\lambda }_2|\ge \dots \ge |{\lambda }_N|$,

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

Załóżmy także, że istnieje baza złożona z wektorów własnych $q_1,\dots ,q_N$ tej macierzy (tak jest np. dla macierzy symetrycznej na mocy twierdzenia o własnościach symetrycznego zadania własnego).

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

$x={\alpha }_1{q}_1+\dots +{\alpha }_N{q}_N$,

wtedy

$Ax=A\left(\sum _i{\alpha }_i{q}_i\right)=\sum _i{\alpha }_{i}A{q}_i=\sum _i{\alpha }_i{\lambda }_i{q}_i$

i w konsekwencji

$A^{k}x=\sum _i{\alpha }_i{\lambda }_i^k{q}_i={\lambda }_1^k({\alpha }_1{q}_1+{\alpha }_2{\left(\frac{{\lambda }_2}{{\lambda }_1}\right)}^k{q}_2+\dots +{\alpha }_N{\left(\frac{{\lambda }_N}{{\lambda }_1}\right)}^k{q}_N)$

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

Szybkość zbieżności metody potęgowej jest liniowa, o współczynniku zależnym od stosunku $|{\lambda }_2/{\lambda }_1|$. W patologicznym przypadku, gdy $|{\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 $|{\lambda }_2/{\lambda }_1|$. Na osi pionowej zaznaczono błąd przybliżenia dominującego wektora własnego.

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

Algorytm Metoda potęgowa

---

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

Warunek normowania ma m.in. na celu zapobieżenie powstawania nadmiaru i niedomiaru (gdy $|{\lambda }_1|<1$, to $||A^{k}x||\to 0$, a gdy $|{\lambda }_1|>1$, to $||A^{k}x||\to \mathrm{\infty }$). Przy okazji, zauważ, że $||y_k|{|}_{\mathrm{\infty }}\to |{\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, $||x_k-x_{k-1}||\le ϵ$ lub warunek małego residuum, $||A{x}_k-{\lambda }_{1,k}{x}_k||\le ϵ$, gdzie ${\lambda }_{1,k}$ jest przybliżeniem ${\lambda }_1$ dostępnym na $k$-tej iteracji.

Gdy $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 $|{\lambda }_2/{\lambda }_1|$. Na osi pionowej zaznaczono błąd przybliżenia dominującego wektora własnego. Oczywiście, gdy ${\lambda }_2={\lambda }_1$ na zbieżność praktycznie nie ma szans.

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

Odwrotna metoda potęgowa

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

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

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

Stwierdzenie O transformacji widma macierzy

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

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

Algorytm Odwrotna metoda potęgowa

---

<math>x_0</math> = dowolny wektor startowy; k = 0;
while( !stop )
{
	Rozwiąż układ równań <math>(A-\sigma I)y_k = x_{k-1}</math>;
	<math>x_k</math> = <math>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 ${\lambda }_j$ jest przesunięcie $\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 $\sigma$, korzystając z dotychczas wyznaczonego wektora $x_k\approx q_j$ i ilorazu Rayleigh:

${\lambda }_j=\frac{q_j^{T}A{q}_j}{q_j^T{q}_j}\approx \frac{x_k^{T}A{x}_k}{x_k^T{x}_k}$

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

Algorytm Metoda RQI

---

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

(wybierając normowanie wektora $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 $q_j$ będzie ${\sigma }_k$, tym bardziej rośnie uwarunkowanie $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 $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 $N\le 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 $O(N^3)$ operacji sprowadzić przez ortogonalne podobieństwo zadanie z macierzą gęstą $A$ (w przypadku niesymetrycznym) do zadania z macierzą Hessenberga, czyli macierzą, której element $(i,j)$ jest zerowy gdy tylko $i-j>1$:

$\left(\begin{array}{cccccc}*& *& *& *& \cdots & *\\ *& *& *& *& \cdots & *\\ & *& *& *& \cdots & *\\ & & \ddots & \ddots & \dots & ⋮\\ & & & \ddots & \ddots & *\\ & & & & *& *\end{array}\right)$

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

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

$A:=Q_{k}A{Q}_k^T$,

gdzie

$Q_k$

jest pewnym przekształceniem Householdera. Rzeczywiście, niech

$A=\left(\begin{array}{cccccc}{d}_1& *& *& *& \cdots & *\\ a_1& *& *& *& \cdots & *\\ a_2& *& *& *& \cdots & *\\ ⋮& *& \ddots & \ddots & \dots & ⋮\\ ⋮& *& *& \ddots & \ddots & *\\ a_{N-1}& *& *& *& *& *\end{array}\right)$

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

$\left(\begin{array}{cc}1& \\ & {\tilde{Q}}_1\end{array}\right)\cdot A=\left(\begin{array}{cccccc}{d}_1& *& *& *& \cdots & *\\ c& *& *& *& \cdots & *\\ & *& *& *& \cdots & *\\ & *& \ddots & \ddots & \dots & ⋮\\ & *& *& \ddots & \ddots & *\\ & *& *& *& *& *\end{array}\right)$

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

$\left(\begin{array}{cc}1& \\ & {\tilde{Q}}_1\end{array}\right)\cdot A\cdot \left(\begin{array}{cc}1& \\ & {\tilde{Q}}_1\end{array}\right)=\left(\begin{array}{cccccc}{d}_1& *& *& *& \cdots & *\\ c& *& *& *& \cdots & *\\ & *& *& *& \cdots & *\\ & *& \ddots & \ddots & \dots & ⋮\\ & *& *& \ddots & \ddots & *\\ & *& *& *& *& *\end{array}\right)$

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

Gdy wyjściowa macierz $A$ jest symetryczna, to z definicji, macierz wynikowa $\left(\begin{array}{cc}I& \\ & {\tilde{Q}}_{N-2}\end{array}\right)\cdots \left(\begin{array}{cc}1& \\ & {\tilde{Q}}_1\end{array}\right)A\left(\begin{array}{cc}1& \\ & {\tilde{Q}}_1\end{array}\right)\cdots \left(\begin{array}{cc}I& \\ & {\tilde{Q}}_{N-2}\end{array}\right)$ 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 $A$; wektory własne macierzy $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 $A$ już w postaci trójdiagonalnej, łatwo widzieć, że "prawie" rozpada się ona na dwie mniejsze macierze trójdiagonalne: dokładniej,

$\left(\begin{array}{cccc}{a}_1& b_1& & \\ b_1& a_2& \ddots & \\ & \ddots & \ddots & b_{N-1}\\ & & b_{N-1}& a_N\end{array}\right)=\left(\begin{array}{cc}{T}_1& \\ & T_2\end{array}\right)+b_{m}u{u}^T$,

gdzie $T_1=\left(\begin{array}{cccc}{a}_1& b_1& & \\ b_1& a_2& \ddots & \\ & \ddots & \ddots & b_{m-1}\\ & & b_{m-1}& a_m-b_m\end{array}\right)$, $T_2=\left(\begin{array}{cccc}{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{array}\right)$ są --- tak jak $A$ --- macierzami trójdiagonalnymi i symetrycznymi (jako podmacierze $A$), tylko o połowę mniejszego wymiaru, gdy $m\approx N/2$. Natomiast $u=e_m+e_{m+1}$, więc macierz $b_{m}u{u}^T$ ma tylko cztery niezerowe elementy, każdy równy $b_m$.

Zgodnie ze swoją nazwą, metoda "dziel i rządź" sprowadza zadanie znajdowania par własnych macierzy wymiaru $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 $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 $T_1,T_2$ umiemy rozwiązać zadanie własne tak, że znamy macierze: $Q_i$ --- ortogonalną oraz $D_i$ --- diagonalną, takie, że

$Q_i^T{T}_i{Q}_i=D_{i}i=1,2$

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

$\left(\begin{array}{cc}{Q}_1^T& \\ & Q_2^T\end{array}\right)(\left(\begin{array}{cc}{T}_1& \\ & T_2\end{array}\right)+b_{m}u{u}^T)\left(\begin{array}{cc}{Q}_1& \\ & Q_2\end{array}\right)=\left(\begin{array}{cc}{D}_1& \\ & D_2\end{array}\right)+b_{m}v{v}^T$

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

Na szczęście łatwo pokazać, że jeśli $\lambda$ nie jest wartością własną macierzy diagonalnej $D=\left(\begin{array}{cc}{D}_1& \\ & D_2\end{array}\right)$, to wartości własne $\lambda$ macierzy $D+b_{m}v{v}^T$

spełniają równanie

$f(\lambda )\equiv 1+b_m\sum _{j=1}^N\frac{v_j^2}{d_j-\lambda }=0$,

gdzie $d_j$ są elementami na diagonali macierzy $D$.

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

W typowym przypadku $f$ będzie miała dokładnie $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 $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>A_1 = A</math>;
for k = 1, 2, ...
{
	wykonaj rozkład <math>A_k = Q_{k}R_{k}</math>;
	<math>A_{k+1} = R_{k}\cdot Q_{k}</math>;
}

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

Algorytm Metoda potęgowa na przestrzeni

---

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

Drugi krok w istocie ortogonalizuje kolumny $W_{k+1}$. Gdyby nie ortogonalizować zestawu wektorów $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 $A$. Zapewniając sobie ortogonalność $V_{k+1}$, możemy liczyć na to, że kolejne kolumny macierzy $V_k$ będą dążyć do wektorów własnych odpowiadających kolejnym wartościom własnym $A$ (przy stosownych założeniach o $A$, m.in. że wszystkie wartości własne $A$ spełniają $|{\lambda }_i|\ne |{\lambda }_j|$ dla $i\ne 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 $R$ były dodatnie) mamy zależności $V_{k+1}=Q_1\cdots Q_k$ oraz $A_{k+1}=V_{k+1}^{T}A{V}_{k+1}$.

Tak więc, w sprzyjających warunkach, metoda QR, jako równoważna metodzie potęgowej na przestrzeni, będzie zbieżna: $A_k\to A_{\mathrm{\infty }}$, gdzie $A_{\mathrm{\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 $A_{\mathrm{\infty }}$ (a więc także $A$) będą liczby na diagonali $A_{\mathrm{\infty }}$.

Twierdzenie O zbieżności metody QR

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

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

Powyższa wersja algorytmu QR jest mało praktyczna, m.in. jest zbieżna wolno i przy poważnych ograniczaniach na $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>A_1 = A</math>;
for k = 1, 2, ...
{
	wybierz sprytnie przesunięcie <math>\sigma_k</math>; 
	wykonaj rozkład <math>A_k - \sigma_kI = Q_{k}R_{k}</math>;
	<math>A_{k+1} = R_{k}\cdot Q_{k} + \sigma_kI</math>;
}

Koszt wyznaczenia wszystkich wektorów i wartości własnych jest rzędu $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 $J_0,J_1,\dots$ sprowadzić wyjściową macierz symetryczną $A$ do postaci (prawie) diagonalnej:

$J_k^T{J}_{k-1}^T\cdots J_0^{T}A{J}_0\cdots J_{k-1}{J}_k\approx \left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_N\end{array}\right)$

Tak więc iteracja Jacobiego ma postać:

Algorytm Metoda Jacobiego

---

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

Macierze $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:

$\omega =\sum _{j<i}{a}_{ij}^2$

Twierdzenie O zbieżności metody Jacobiego

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

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

gdzie ${\omega }_k$ oznacza miarę niediagonalności na $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 $A\in R^{N\times N}$ będzie diagonalizowalna, to znaczy dla pewnej macierzy $X$ zachodzi

$X^{-1}AX=\left(\begin{array}{ccc}{\lambda }_1& & \\ & \ddots & \\ & & {\lambda }_N\end{array}\right)$,

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

$|{\lambda }_j-\tilde{\lambda }|\le {\text{cond}}_2(X)||A-\tilde{A}|{|}_2$

Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia $X$ jest ortogonalna, $X^{-1}=X^T$, to mamy ${\text{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 $A$ jest rzeczywista i symetryczna, to

$\underset{j=1,\dots ,N}{\min }|{\lambda }_j-\tilde{\lambda }|\le ||A-\tilde{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

$A_{ϵ}=\left(\begin{array}{cc}a& 1\\ ϵ& a\end{array}\right)$

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

$|{\lambda }_{ϵ}-{\lambda }_0|/||A_{ϵ}-A_0||=\sqrt{ϵ}/ϵ\to \mathrm{\infty }$,

gdy $ϵ\to 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 $A$ dla ujemnego parametru $ϵ$ są zespolone!

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

Bardziej spektakularny przykład pochodzi od Wilkinsona:

Przykład: Perfidny wielomian Wilkinsona

Niech

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

Zmiana współczynnika przy ${\lambda }^{19}$ o $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.