MN13

From Studia Informatyczne

Spis treści

1 Wektory i wartości własne

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 A x = \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^TQ = 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} .$

[Edytuj]

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 kołami Gerszgorina oraz zgrubna lokalizacja wewnątrz okręgu o promieniu równym . Dokładne wartości własne zaznaczone trójkącikami.

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.

[Edytuj]

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!...

[Edytuj]

Metoda potęgowa

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

(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_1q_1 + \ldots + \alpha_Nq_N,$

wtedy

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

i w konsekwencji

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

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_1q_1$ i w konsekwencji wektory $\displaystyle x_k = A^kx$ 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$ ).

Zbieżność metody potęgowej w zależności od stosunku . Na osi pionowej zaznaczono błąd przybliżenia dominującego wektora własnego.

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^kx|| \rightarrow 0$ , a gdy $\displaystyle |\lambda_1| > 1$ , to $\displaystyle ||A^kx|| \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 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 . Na osi pionowej zaznaczono błąd przybliżenia dominującego wektora własnego. Oczywiście, gdy na zbieżność praktycznie nie ma szans.

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.

[Edytuj]

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++;	
}

[Edytuj]

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^TAq_j}{q_j^Tq_j} \approx \frac{x_k^TAx_k}{x_k^Tx_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!

[Edytuj]

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 A Q_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}_1a = 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.

[Edytuj]

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 dla macierzy jednowymiarowego laplasjanu rozmiaru 10. Zwróć uwagę na asymptoty pionowe tej funkcji oraz jej przedziałową monotoniczność.

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łą.

[Edytuj]

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}^TA_kQ_{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}^TAV_{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.

[Edytuj]

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 A J_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.

[Edytuj]

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} A X = \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 (z parametrem ) w otoczeniu

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.

[Edytuj]

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.

[Edytuj]

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.

Źródło: "http://wazniak.mimuw.edu.pl/index.php?title=MN13"

MN13

From Studia Informatyczne

Spis treści

Wektory i wartości własne

Lokalizacja wartości własnych

Wyznaczanie pojedynczej pary własnej

Metoda potęgowa

Odwrotna metoda potęgowa

Metoda Rayleigh (RQI)

Metody rozwiązywania pełnego zadania własnego

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

Metoda QR

Metoda Jacobiego dla macierzy symetrycznej

Uwarunkowanie

Biblioteki

Literatura

Nawigacja

Szukaj

Napisz do nas