MN12

Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi� poprawki!

Wyznaczanie wektorów i wartości własnych

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 do 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), to 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

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. Uzupe�nij: 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.

Przykład

Praktyczne zadanie z macierzą symetryczną

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?)

• 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łanym

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).$

Uwarunkowanie zadania

Twierdzenie Bauer-Fike

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 Uzupe�nij: 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 $\delta =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 presunię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 w/w zaburzenia mamy dodatkowo 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 i zależy m.in. od uwarunkowania wartości własnych (czego łatwo się domyślić) oraz od tego, jak blisko siebie leżą wartości własne.

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$. W tym celu 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ą.

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

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.

Metoda potęgowa, odwrotna potęgowa, RQI

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 poszykiwane 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 Uzupe�nij: 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).$

Ponieważ z założenia, że istnieje dokładnie jedna dominująca wartość własna, $\left|\frac{{\lambda }_N}{{\lambda }_1}\right|<1$, to 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.

W praktyce nie wyznaczamy wzorem $x_k=(A^k)\cdot x$, lecz raczej 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 $|{\lambda }_1|<1$, to $||A^{k}x||\to 0$, a gdy $|{\lambda }_1|>1$, to $||A^{k}x||\to \mathrm{\infty }$). Przy okazji, $||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.

Plik:MNXXX.png

Zasada działania metody potęgowej

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$

Łącząc te dwie własności mamy, że

Stwierdzenie Transformacja 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>\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

Z własności metody potęgowej, 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}$

Algorytm Metoda RQI (Rayleigh Quotient Iteration)

---

<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 $x$ w normie euklidesowej upraszczamy co nieco
algorytm).

Wielką zaletą metody RQI jest jej szybkość zbieznoś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 pomaga...

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 tak jest ---

 wektor błędu ma kierunek praktycznie zgodny z kierunkiem poszukiwanego wektora

$q_j$, a tym samym tylko pomaga 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$ do zadania z macierzą Hessenberga (w przypadku niesymetrycznym)

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

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}$, tak 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 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>\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 $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 teoretycznemu algorytmowi iteracji prostej zastosowanemu nie do pojedynczego wektora, ale do $N$ wektorów naraz:

Algorytm Iteracja prosta 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 $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 iteracji prostej na podprzestrzeni, 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), a tym samym wartościami własnymi $A_{\mathrm{\infty }}$ (a więc także $A$) będą liczby na diagonali $A_{\mathrm{\infty }}$.

Twierdzenie Zbieżność metody QR w szczególnym przypadku

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>\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 $O(N^3)$ ze stałą równą około 30.

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, opcjonalnie, wektory własne). Dla macierzy symetrycznych mamy jeszcze m.in. funkcje DSYEVX (dla wybranych wartości własnych) i DSYEVD (z algorytmem dzieli 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. Jak dotąd, tylko MATLAB potrafi skorzystać z ARPACKa dla wyznaczenia fragmentów widma macierzy rzadkiej, za pomocą funkcji eigs.