MN12LAB

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 {hopla} wartości własnych !hopla! 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 [Uzupelnij]

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

{}{Model wieżowca poddanego drganiom poprzecznym}

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 {niesymetrycznej} macierzy wymiaru $2N$, która powstaje ze współczynników równania różniczkowego opisującego dynamikę tego układu.

Przykład [Uzupelnij]

[Macierz Google'a] Podstawowy algorytm rankingowania stron WWW w {http://www.wikipedia.org/pagerank}{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 XXXX 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. {sec:metoda-potegowa}{metody potęgowej} numerycznego wyznaczania takiego wektora.

Przykład [Uzupelnij]

[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! właśnie w taki sposób wyznacza pierwiastki wielomianów: jako wartości własne macierzy stowarzyszonej.

Przykład [Uzupelnij]

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

• 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 [Uzupelnij]

[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:

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \Lambda = \beginpmatrix \lambda_1 & & \\ & \ddots & \\ & & \lambda_N\endpmatrix . }

Uwarunkowanie zadania

Twierdzenie [Uzupelnij]

[Bauer-Fike]

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

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle X^{-1} A X = \beginpmatrix \lambda_1 & & \\ & \ddots & \\ & & \lambda_N\endpmatrix , }

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 [Uzupelnij]

[Wartości własne macierzy symetrycznej są doskonale uwarunkowane] Przy oznaczeniach jak {thm:Bauer-Fike}{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 [Uzupelnij]

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle A_\epsilon = \beginpmatrix a & 1 \\ \epsilon & a \endpmatrix }

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!

{eigencond.png}{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 [Uzupelnij]

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

{wilkinson.png}{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 [Uzupelnij]

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 [Uzupelnij]

[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 [Uzupelnij]

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

{gershgorindisks.png}{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.}

Przykład [Uzupelnij]

[Widmo macierzy jednowymiarowego Laplasjanu] Norma daje:

Tw. Gerszgorina daje:

W rzeczywistości,

Metoda potęgowa, odwrotna potęgowa, RQI

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 {thm:symetric-eig}{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

[title=Metoda potęgowa] $x_0$ = dowolny wektor startowy; k = 0; while ( !stop ) { $y_k$ = $A{x}_{k-1}$; $x_k$ = $y_k/||y_k|{|}_{\mathrm{\infty }}$; 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.

{}{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 [Uzupelnij]

[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:

[title=Odwrotna metoda potęgowa] $x_0$ = dowolny wektor startowy; k = 0; while( !stop ) { Rozwiąż układ równań $(A-\sigma I)y_k=x_{k-1}$; $x_k$ = $y_k/||y_k|{|}_{\mathrm{\infty }}$; 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}$

[title=Metoda RQI (Rayleigh Quotient Iteration)] $x_0$ = dowolny wektor startowy; ${\sigma }_0$ = przybliżenie ${\lambda }_j$; k = 0; while( !stop ) { Rozwiąż układ równań $(A-{\sigma }_{k}I)y_k=x_{k-1}$; $x_k$ = $y_k/||y_k|{|}_2$; ${\sigma }_{k+1}$ = $x_k^{T}A{x}_k$; 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 [Uzupelnij]

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

Metoda dziel i rządź, metoda QR

{}{Secular equation}