MN12: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Przykry (dyskusja | edycje)
Przykry (dyskusja | edycje)
mNie podano opisu zmian
Linia 1: Linia 1:
''Uwaga: przekonwertowane latex2mediawiki;
prawdopodobnie trzeba wprowadzi� poprawki!''


==Wyznaczanie wektorów i wartości własnych==
Niech będzie dana rzeczywista kwadratowa macierz <math>\displaystyle A</math> wymiaru <math>\displaystyle N</math>. Wektorem własnym <math>\displaystyle x\in C^N</math> oraz
odpowiadającą mu wartością własną <math>\displaystyle \lambda \in C</math> nazwiemy taką parę, dla której
<center><math>\displaystyle A x = \lambda x,
</math></center>
przy czym <math>\displaystyle x\neq 0</math>.
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.
{{przyklad|Odporność budynku na trzęsienie ziemi||
Rozważmy prosty układ mechaniczny opisujący, naturalnie w pewnym
jedynie przybliżeniu, zachowanie się układu <math>\displaystyle N</math> 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 <math>\displaystyle 2N</math>,
która powstaje ze współczynników  równania różniczkowego opisującego dynamikę
tego układu.
}}
{{przyklad|Macierz Google'a||
Podstawowy algorytm rankingowania stron WWW w [http://www.wikipedia.org/pagerank  wyszukiwarce Google]
sprowadza się do znalezienia rzeczywistego ''wektora własnego'' <math>\displaystyle \pi</math> pewnej silnie
rozrzedzonej macierzy <math>\displaystyle A</math> (gigantycznego rozmiaru, równego liczbie indeksowanych
stron, czyli w chwili pisania tego tekstu około <math>\displaystyle 2.5\cdot 10^{10}</math> stron), odpowiadającego wartości własnej równej 1:
<center><math>\displaystyle A \pi = \pi.
</math></center>
Współrzędne wektora <math>\displaystyle \pi</math>
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 <math>\displaystyle A</math> gwarantują, że taki
wektor <math>\displaystyle \pi</math> zawsze istnieje i jest jedyny! Co więcej, wartość 1 jest
dominującą wartością własną <math>\displaystyle A</math>, a to z kolei ma ważne znaczenie dla tzw.
[[sec:metoda-potegowa|Uzupe�nij: metody potęgowej]] numerycznego wyznaczania takiego wektora.
}}
{{przyklad|Wyznaczanie miejsc zerowych wielomianu||
Jak wiadomo, wartości własne to miejsca zerowe wielomianu charakterystycznego
macierzy <math>\displaystyle P(\lambda) = \det(A - \lambda I)</math>. Zachodzi także fakt odwrotny, to
znaczy miejsca zerowe wielomianu są wartościami pewnej macierzy, np. miejsca
zerowe wielomianu
<center><math>\displaystyle p(\lambda) = p_1 \lambda^N + \ldots + p_N \lambda + p_{N+1}
</math></center>
są wartościami własnymi m.in. macierzy stowarzyszonej,
<center><math>\displaystyle A = \begin{pmatrix} 
-p_2/p_1 & -p_3/p_1 & \cdots & -p_{N+1}/p_1\\
1 & & & \\
& 1 & & \\
& & \ddots & \\
& & & 1
\end{pmatrix}
</math></center>
Funkcja Octave'a <code>compan(p)</code> wyznacza macierz stowarzyszoną dla zadanego
wielomianu o współczynnikach w wektorze <math>\displaystyle p = [p_1,\ldots,p_N, p_{N+1}]^T</math>. Z tej
macierzy korzysta następnie funkcja Octave'a <code>roots</code>, która  właśnie w taki
sposób wyznacza pierwiastki wielomianów: jako wartości własne macierzy
stowarzyszonej.
}}
{{przyklad|||
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 <math>\displaystyle A</math> wymiaru <math>\displaystyle N</math> ma rozkład
<center><math>\displaystyle A = Q\Lambda Q^T,
</math></center>
gdzie <math>\displaystyle Q\in R^{N\times N}</math> jest ortogonalna (tzn. <math>\displaystyle Q^TQ = I</math>), a jej kolumnami są
wektory własne <math>\displaystyle A</math>, natomiast <math>\displaystyle \Lambda\in
R^N</math> jest diagonalna z
wartościami własnymi <math>\displaystyle A</math> na diagonali:
<center><math>\displaystyle \Lambda = \begin{pmatrix} \lambda_1 & & \\ & \ddots & \\ & &
\lambda_N\end{pmatrix} .
</math></center>
}}
===Uwarunkowanie zadania===
{{twierdzenie|Bauer-Fike||
Niech <math>\displaystyle A\in R^{N\times N}</math> będzie diagonalizowalna, to
znaczy dla pewnej macierzy <math>\displaystyle X</math> zachodzi
<center><math>\displaystyle X^{-1}  A X = \begin{pmatrix}  \lambda_1 & & \\ & \ddots & \\ & &
\lambda_N\end{pmatrix} ,
</math></center>
a więc (gdyż macierz po prawej stronie jest podobna do <math>\displaystyle A</math>) <math>\displaystyle \lambda_i\in C</math>,
<math>\displaystyle i=1,\ldots,N</math> są
wartościami własnymi <math>\displaystyle A</math>. Rozważmy macierz zaburzoną <math>\displaystyle \widetilde{A}</math> i jakąś jej
wartość własną <math>\displaystyle \widetilde{\lambda}</math>. Wtedy istnieje wartość własna <math>\displaystyle \lambda_j</math>
macierzy <math>\displaystyle A</math> taka, że
<center><math>\displaystyle |\lambda_j - \widetilde{\lambda}| \leq  \mbox{cond} _2(X) ||A - \widetilde{A}||_2.
</math></center>
}}
Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia <math>\displaystyle X</math> jest
ortogonalna,
<math>\displaystyle X^{-1} = X^T</math>, to mamy <math>\displaystyle  \mbox{cond} _2(X) = 1</math> i w konsekwencji zachodzi
{{wniosek|Wartości własne macierzy symetrycznej są doskonale uwarunkowane||
Przy oznaczeniach jak [[thm:Bauer-Fike|Uzupe�nij: twierdzeniu Bauera-Fike'a]], jeśli
dodatkowo założymy, że macierz <math>\displaystyle A</math> jest rzeczywista i symetryczna, to
<center><math>\displaystyle \min_{j=1,\ldots,N}|\lambda_j - \widetilde{\lambda}| \leq ||A - \widetilde{A}||_2.
</math></center>
}}
Z drugiej strony, dla macierzy niediagonalizowalnych, uwarunkowanie wartości
własnych może być
dowolnie duże, co ilustruje poniższy
{{przyklad|||
<center><math>\displaystyle A_\epsilon = \begin{pmatrix}  a & 1 \\ \epsilon & a \end{pmatrix}
</math></center>
Weźmy dla uproszczenia <math>\displaystyle a=0</math>.
Wartości własne <math>\displaystyle A_\epsilon</math> to zera wielomianu <math>\displaystyle p_\epsilon(\lambda) = \lambda^2 - \epsilon</math>,
zatem <math>\displaystyle \lambda_\epsilon = \pm \sqrt{\epsilon}</math> i w konsekwencji
<center><math>\displaystyle |\lambda_\epsilon - \lambda_0| / ||A_\epsilon - A_0|| = \sqrt{\epsilon}/\epsilon
\rightarrow \infty,
</math></center>
gdy <math>\displaystyle \epsilon \rightarrow 0^+</math>, 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 <math>\displaystyle A</math> dla
ujemnego parametru <math>\displaystyle \epsilon</math> są zespolone!
[[Image:MNeigencond.png|frame|200px|Zachowanie się wartości własnych macierzy <math>\displaystyle A</math> (z
parametrem <math>\displaystyle a=1</math>) w otoczeniu <math>\displaystyle \delta = 0</math>]]
}}
Bardziej spektakularny przykład pochodzi od Wilkinsona:
{{przyklad|Perfidny wielomian Wilkinsona||
Niech
<center><math>\displaystyle p(\lambda) = (\lambda -1)(\lambda - 2) \cdots (\lambda - 20).
</math></center>
Zmiana współczynnika przy <math>\displaystyle \lambda^{19}</math> o <math>\displaystyle 10^{-7}</math> 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.
[[Image:MNwilkinson.png|frame|200px|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 <math>\displaystyle A</math>. W tym celu mogą być nam
pomocne dwa fakty:
{{fakt|||
Dowolna wartość własna <math>\displaystyle \lambda\in C</math> macierzy <math>\displaystyle A</math> spełnia
<center><math>\displaystyle |\lambda| \leq ||A||,
</math></center>
gdzie <math>\displaystyle ||A||</math> jest dowolną normą macierzową indukowaną przez normę wektorową.
}}
Rzeczywiście, skoro istnieje wektor <math>\displaystyle x\neq 0</math> taki, że <math>\displaystyle Ax = \lambda x</math>, to stąd
<math>\displaystyle ||Ax||/||x|| = |\lambda|</math>, więc fakt powyższy wynika już z definicji normy
macierzy:
<center><math>\displaystyle ||A|| = \max_{y\neq 0}\frac{||Ay||}{||y||} \geq ||Ax||/||x||.
</math></center>
Drugie twierdzenie jest równie proste w dowodzie, ale daje trochę więcej
informacji o lokalizacji widma.
{{twierdzenie|Gerszgorina||
Wartości własne macierzy <math>\displaystyle A</math> leżą w sumie (teoriomnogościowej) dysków <math>\displaystyle K_i</math> na
płaszczyźnie zespolonej,
<center><math>\displaystyle K_i = \{z \in C: |z - a_{ii}| \leq \sum_{j\neq i} |a_{ij}| \}, \qquad i =
1,\ldots N.
</math></center>
}}
{{przyklad|Koła Gerszgorina||
Niech 
<center><math>\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}
</math></center>
[[Image:MNgershgorindisks.png|frame|200px|Lokalizacja wartości własnych macierzy <math>\displaystyle A</math> kołami Gerszgorina oraz zgrubna
lokalizacja wewnątrz okręgu
o promieniu równym <math>\displaystyle ||A||_1</math>. Dokładne wartości własne zaznaczone trójkącikami.]]
}}
===Metoda potęgowa, odwrotna potęgowa, RQI===
Jak wiemy z algebry, nawet gdy <math>\displaystyle A</math> 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 <math>\displaystyle A</math> 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 <math>\displaystyle A\in R^{N\times N}</math> spełniają
<center><math>\displaystyle |\lambda_1| > |\lambda_2| \geq \ldots \geq |\lambda_N|,
</math></center>
(to znaczy, istnieje dokładnie jedna ''dominująca'' wartość własna macierzy
<math>\displaystyle A</math>.
Załóżmy także, że istnieje baza złożona z wektorów własnych <math>\displaystyle q_1,\ldots,q_N</math> tej
macierzy (tak jest np. dla macierzy symetrycznej na mocy
[[thm:symetric-eig|Uzupe�nij: twierdzenia o własnościach symetrycznego zadania
własnego]]).
Kierunek własny <math>\displaystyle q_k</math> jakiejś macierzy <math>\displaystyle A</math> ma taką własność, że poddany działaniu przekształcenia
<math>\displaystyle A</math> wydłuża się <math>\displaystyle \lambda_k</math> razy, wobec tego, dowolny wektor <math>\displaystyle x\in R^N</math> poddany
działaniu <math>\displaystyle A</math> najbardziej wydłuży się w kierunku <math>\displaystyle q_1</math>. Iterując tę procedurę,
powinniśmy dostawać w wyniku wektory, w których coraz bardziej dominuje kierunek
<math>\displaystyle q_1</math>. Formalnie, niech
<center><math>\displaystyle x = \alpha_1q_1 + \ldots + \alpha_Nq_N,
</math></center>
wtedy
<center><math>\displaystyle Ax = A \left( \sum_i \alpha_iq_i \right) = \sum_i \alpha_i A q_i
=  \sum_i \alpha_i \lambda_i q_i
</math></center>
i w konsekwencji
<center><math>\displaystyle A^kx = \sum_i \alpha_i \lambda_i^k q_i = \lambda_1^k\left(\alpha_1q_1 +
\alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^kq_2 + \ldots  +
\alpha_N\left(\frac{\lambda_N}{\lambda_1}\right)^kq_N \right).
</math></center>
Ponieważ z założenia, że istnieje dokładnie jedna dominująca wartość własna,
<math>\displaystyle \left|\frac{\lambda_N}{\lambda_1}\right| < 1</math>, to wyrażenie w nawiasie dąży do
<math>\displaystyle \alpha_1q_1</math> i w konsekwencji wektory <math>\displaystyle x_k = A^kx</math> dążą, gdy
<math>\displaystyle k\rightarrow\infty</math>, do kierunku wektora własnego <math>\displaystyle q_1</math>, to znaczy wektora
odpowiadającego dominującej wartości własnej <math>\displaystyle A</math>  (o ile tylko <math>\displaystyle \alpha_1
\neq 0</math>).
Szybkość zbieżności metody potęgowej jest liniowa, o współczynniku zależnym od
stosunku <math>\displaystyle \lambda_2/\lambda_1|</math>. W patologicznym przypadku, gdy <math>\displaystyle |\lambda_1|
\approx |\lambda_2|</math>, może więc okazać się, że metoda praktycznie nie jest
zbieżna.
W praktyce nie wyznaczamy wzorem <math>\displaystyle x_k = (A^k)\cdot x</math>, lecz raczej korzystamy z
metody iteracyjnej
{{algorytm|Metoda potęgowa||
<pre>
<math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; k <nowiki>=</nowiki> 0;
while( !stop )
{
<math>\displaystyle y_k</math> <nowiki>=</nowiki> <math>\displaystyle Ax_{k-1}</math>;
<math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_\infty</math>;
k++;
}
</pre>}}
Warunek normowania ma m.in. na celu zapobieżenie powstawania nadmiaru i
niedomiaru (gdy <math>\displaystyle |\lambda_1| < 1</math>, to <math>\displaystyle ||A^kx|| \rightarrow 0</math>, a gdy
<math>\displaystyle |\lambda_1| > 1</math>, to <math>\displaystyle ||A^kx|| \rightarrow \infty</math>). Przy okazji,
<math>\displaystyle ||y_k||_\infty \rightarrow |\lambda_1|</math>, 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, <math>\displaystyle ||x_k -
x_{k-1}|| \leq \epsilon</math>, lub warunek małego residuum, <math>\displaystyle ||Ax_k - \lambda_{1,k}
x_k||\leq \epsilon</math>, gdzie <math>\displaystyle \lambda_{1,k}</math> jest przybliżeniem <math>\displaystyle \lambda_1</math>
dostępnym na <math>\displaystyle k</math>-tej iteracji.
[[Image:MNXXX.png|frame|200px|Zasada działania metody potęgowej]]
Metoda potęgowa doskonale sprawdza się, gdy macierz <math>\displaystyle A</math> jest macierzą
rozrzedzoną --- np. w przypadku macierzy Google'a.
====Odwrotna metoda potęgowa====
Zauważmy, że dla dowolnej macierzy kwadratowej <math>\displaystyle A</math> o wartościach własnych
<math>\displaystyle \lambda_k</math> i odpowiadających im wektorach własnych <math>\displaystyle q_k</math>, mamy:
* Macierz <math>\displaystyle A-\sigma I</math> ma wartości własne <math>\displaystyle \lambda_k - \sigma</math> oraz wektory
własne <math>\displaystyle q_k</math>,
* Jeśli dodatkowo <math>\displaystyle A</math> jest nieosobliwa, to macierz <math>\displaystyle A^{-1}</math> ma wartości
własne <math>\displaystyle 1/\lambda_k</math> oraz wektory własne <math>\displaystyle q_k</math>
Łącząc te dwie własności mamy, że
{{stwierdzenie|Transformacja widma macierzy||
Macierz <math>\displaystyle (A-\sigma I)^{-1}</math> (o ile istnieje),
to ma wartości własne równe <math>\displaystyle \frac{1}{\lambda_k - \sigma}</math> i wektory własne
identyczne z <math>\displaystyle A</math>.
}}
Skoro tak, to jeśli najbliższą <math>\displaystyle \sigma</math>  wartością własną <math>\displaystyle A</math> jest <math>\displaystyle \lambda_j</math>,
wówczas metoda potęgowa zastosowana do macierzy <math>\displaystyle (A-\sigma I)^{-1}</math> zbiegnie do
<math>\displaystyle q_j</math>. To prowadzi do następującego algorytmu, odwrotnej metody potęgowej:
{{algorytm|Odwrotna metoda potęgowa||
<pre>
<math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; k <nowiki>=</nowiki> 0;
while( !stop )
{
Rozwiąż układ równań <math>\displaystyle (A-\sigma I)y_k = x_{k-1}</math>;
<math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_\infty</math>;
k++;
}
</pre>}}
====Metoda Rayleigh====
Z własności metody potęgowej, metoda odwrotna potęgowa jest zbieżna tym
szybciej, im bliżej <math>\displaystyle \lambda_j</math> jest przesunięcie <math>\displaystyle \sigma</math> (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 <math>\displaystyle \sigma</math>,
korzystając z dotychczas wyznaczonego wektora <math>\displaystyle x_k \approx q_j</math> i ilorazu
Rayleigh:
<center><math>\displaystyle \lambda_j = \frac{q_j^TAq_j}{q_j^Tq_j} \approx \frac{x_k^TAx_k}{x_k^Tx_k}
</math></center>
{{algorytm|Metoda RQI (Rayleigh Quotient Iteration)||
<pre>
<math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; <math>\displaystyle \sigma_0</math> <nowiki>=</nowiki> przybliżenie <math>\displaystyle \lambda_j</math>; k <nowiki>=</nowiki> 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> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_2</math>;
<math>\displaystyle \sigma_{k+1}</math> <nowiki>=</nowiki> <math>\displaystyle x_k^TAx_k</math>;
k++;
}
</pre>}}
(wybierając normowanie wektora <math>\displaystyle x</math> 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 <math>\displaystyle q_j</math> będzie <math>\displaystyle \sigma_k</math>, tym
bardziej rośnie uwarunkowanie <math>\displaystyle A-\sigma_k I</math>, 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
<math>\displaystyle q_j</math>, 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 <math>\displaystyle N \leq 25</math>).
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 <math>\displaystyle O(N^3)</math> operacji sprowadzić przez
ortogonalne podobieństwo zadanie z
macierzą gęstą <math>\displaystyle A</math> do zadania z macierzą Hessenberga (w przypadku
niesymetrycznym)
<center><math>\displaystyle \begin{pmatrix}
* & * & * & *      & \cdots & * \\
* & * & * & *      & \cdots & * \\
  & * & * & *      & \cdots & * \\
  &  & \ddots  & \ddots  & \ldots &  \vdots \\
  &  &  & \ddots  & \ddots & * \\
  &  &  &        &  *    & *
\end{pmatrix}
</math></center>
bądź wręcz trójdiagonalną, gdy <math>\displaystyle A</math> była symetryczna.
Każdą macierz kwadratową <math>\displaystyle A</math> da się sprowadzić do postaci Hessenberga sekwencją
przekształceń postaci
<center><math>\displaystyle
A := Q_k A Q_k^T,
</math></center>
gdzie <math>\displaystyle Q_k</math> jest pewnym przekształceniem Householdera. Niech <center><math>\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}
</math></center>
i oznaczmy <math>\displaystyle a = [a_1,\ldots,a_{N-1}]^T</math>. Możemy
wziąć na początek przekształcenie Householdera <math>\displaystyle \widetilde{Q}_1</math> takie, że
<math>\displaystyle \widetilde{Q}_1a = c\cdot
e_1</math>, gdzie <math>\displaystyle e_1 = [1,0,\ldots, 0]^T</math>. Wtedy
<center><math>\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}
</math></center>
To samo przekształcenie przyłożone z prawej strony zachowa pierwszą kolumnę i
w efekcie  nie zmieni struktury macierzy:
<center><math>\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} .
</math></center>
Dalej stosujemy tę samą metodę do podmacierzy wymiaru <math>\displaystyle N-1</math>, itd. aż dochodzimy
do macierzy Hessenberga.
Gdy wyjściowa macierz <math>\displaystyle A</math> jest symetryczna, to z definicji, macierz wynikowa
<math>\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} </math> 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 <math>\displaystyle A</math>; wektory własne macierzy <math>\displaystyle A</math> 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 --- <code>DSYEVD</code> w LAPACKu).
Startując z symetrycznej macierzy <math>\displaystyle A</math> już w postaci trójdiagonalnej, łatwo
widzieć, że "prawie" rozpada się ona na dwie mniejsze macierze trójdiagonalne:
dokładniej,
<center><math>\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,
</math></center>
gdzie <math>\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}
</math>, <math>\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}
</math> są --- tak jak <math>\displaystyle A</math> --- macierzami
trójdiagonalnymi i symetrycznymi (jako podmacierze <math>\displaystyle A</math>), tylko o połowę mniejszego
wymiaru, gdy <math>\displaystyle m \approx N/2</math>. Natomiast <math>\displaystyle u = e_{m} + e_{m+1}</math>, tak więc macierz
<math>\displaystyle b_{m} uu^T</math> ma tylko cztery niezerowe elementy, każdy równy <math>\displaystyle b_m</math>.
Zgodnie ze swoją nazwą, metoda dziel i rządź sprowadza zadanie znajdowania par
własnych macierzy wymiaru <math>\displaystyle N</math> 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 <math>\displaystyle 1\times 1</math> --- 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 <math>\displaystyle T_1,T_2</math> umiemy
rozwiązać zadanie własne tak, że znamy macierze: <math>\displaystyle Q_i</math> --- ortogonalną oraz
<math>\displaystyle D_i</math> --- diagonalną, takie, że
<center><math>\displaystyle
Q_i^T T_i Q_i = D_i \qquad i=1,2.
</math></center>
Wtedy łatwo widzieć, że dla łatwo wyznaczalnego wektora <math>\displaystyle v</math>,
<center><math>\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.
</math></center>
W ten sposób zadanie własne dla oryginalnej macierzy <math>\displaystyle T</math> wymiaru <math>\displaystyle N</math> jest
równoważne zadaniu własnemu macierzy diagonalnej zaburzonej o macierz rzędu 1.
Na szczęście łatwo pokazać, że jeśli <math>\displaystyle \lambda</math> nie jest wartością własną
macierzy diagonalnej <math>\displaystyle D = \begin{pmatrix}
D_1 & \\
      & D_2
\end{pmatrix} </math>, to  wartości własne <math>\displaystyle \lambda</math> macierzy <math>\displaystyle D+ b_{m} vv^T</math>
spełniają równanie
<center><math>\displaystyle
f(\lambda) \equiv 1 + b_{m} \sum_{j=1}^N\frac{v_j^2}{d_j - \lambda} = 0,
</math></center>
gdzie <math>\displaystyle d_j</math> są elementami na diagonali macierzy <math>\displaystyle D</math>.
[[Image:MNsecular.png|frame|50|Wykres <math>\displaystyle f(\lambda)</math> dla macierzy jednowymiarowego
laplasjanu rozmiaru 10. Zwróć uwagę na asymptoty pionowe tej funkcji oraz jej
monotoniczność.]]
W typowym przypadku <math>\displaystyle f</math> będzie miała dokładnie <math>\displaystyle N</math> 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 <math>\displaystyle O(N^3)</math> z
małą stałą.
====Metoda QR====
Dla zadania własnego z macierzą niesymetryczną najczęściej stosuje się metodę
QR.
Jakkolwiek ostateczna wersja metody <math>\displaystyle QR</math> 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||
<pre>
<math>\displaystyle A_1 = A</math>;
for k <nowiki>=</nowiki> 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>;
}
</pre>}}
Można sprawdzić, że <math>\displaystyle A, A_1, A_2,\ldots</math> mają te same wartości własne, bo <math>\displaystyle A_{k+1} =
Q_{k+1}^TA_kQ_{k+1}</math>. Co więcej, powyższy algorytm (gdy <math>\displaystyle A</math> jest nieosobliwa) w
zasadzie jest równoważny teoretycznemu algorytmowi iteracji prostej
zastosowanemu nie do pojedynczego wektora, ale do <math>\displaystyle N</math> wektorów naraz:
{{algorytm|Iteracja prosta na przestrzeni||
<pre>
<math>\displaystyle V_1 = I</math>;
for k <nowiki>=</nowiki> 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;
}
</pre>}}
Drugi krok w istocie ortogonalizuje kolumny <math>\displaystyle W_{k+1}</math>. Gdyby nie ortogonalizować
zestawu wektorów <math>\displaystyle W_{k+1}</math>, oczywiście dostalibyśmy w efekcie zbieżność
wszystkich kolumn macierzy do tego samego wektora --- odpowiadającego
dominującej wartości własnej <math>\displaystyle A</math>. Zapewniając sobie ortogonalność <math>\displaystyle V_{k+1}</math>,
możemy liczyć na to, że kolejne kolumny macierzy <math>\displaystyle V_k</math> będą dążyć do wektorów
własnych odpowiadających kolejnym wartościom własnym <math>\displaystyle A</math> (przy stosownych
założeniach o <math>\displaystyle A</math>, m.in. że
wszystkie wartości własne <math>\displaystyle A</math> spełniają <math>\displaystyle |\lambda_i| \neq |\lambda_j|</math> dla <math>\displaystyle i
\neq j</math>). 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 <math>\displaystyle R</math> były dodatnie) mamy zależności <math>\displaystyle V_{k+1} =
Q_1\cdots Q_k</math> oraz <math>\displaystyle A_{k+1} = V_{k+1}^TAV_{k+1}</math>.
Tak więc, w sprzyjających warunkach, metoda QR, jako równoważna
iteracji prostej na podprzestrzeni, będzie zbieżna: <math>\displaystyle A_k \rightarrow A_\infty</math>,
gdzie <math>\displaystyle A_\infty</math> 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 <math>\displaystyle A_\infty</math>
(a więc także <math>\displaystyle A</math>)
będą liczby na diagonali <math>\displaystyle A_\infty</math>.
{{twierdzenie|Zbieżność metody QR w szczególnym przypadku||
Niech wartości własne <math>\displaystyle A\in R^{N\times N}</math> spełniają <math>\displaystyle |\lambda_1|,\ldots,
|\lambda_N| > 0</math> oraz macierz <math>\displaystyle T = [x_1,\ldots,x_N]</math> o kolumnach <math>\displaystyle x_i</math> złożonych z
kolejnych wektorów własnych <math>\displaystyle A</math> ma taką własność, że <math>\displaystyle T^{-1}</math> ma rozkład LU,
<math>\displaystyle T^{-1} = LU</math>.
Wtedy w metodzie QR, ciąg macierzy <math>\displaystyle Q_k</math> jest zbieżny do macierzy diagonalnej, a
ciąg <math>\displaystyle A_k</math> ma podciąg zbieżny do macierzy trójkątnej, której elementy diagonalne
<math>\displaystyle u_{ii}</math> są równe <math>\displaystyle \lambda_i</math> dla <math>\displaystyle i = 1,\ldots, N</math>.
}}
Powyższa wersja algorytmu QR jest mało praktyczna, m.in. jest zbieżna wolno i
przy poważnych ograniczaniach na <math>\displaystyle A</math>. 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||
<pre>
<math>\displaystyle A_1 = A</math>;
for k <nowiki>=</nowiki> 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>;
}
</pre>}}
Koszt wyznaczenia wszystkich wektorów i wartości własnych jest rzędu <math>\displaystyle O(N^3)</math> 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. <code>DGEEV</code> dla macierzy niesymetrycznych
oraz <code>DSYEV</code> 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 <code>DSYEVX</code> (dla wybranych
wartości własnych) i <code>DSYEVD</code> (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 <code>eig</code> 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 <code>eigs</code>.

Wersja z 18:50, 28 sie 2006