PKow

Wyznaczanie wektorów i wartości własnych

Przykład Moj przykład

Tekst

Algorytm Nie robiący nic

---

 Leż

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.

{}{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 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 {ilu} 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 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,

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \displaystyle A = \beginpmatrix -p_2/p_1 & -p_3/p_1 & \cdots & -p_{N+1}/p_1\\ 1 & & & \\ & 1 & & \\ & & \ddots & \\ & & & 1 \endpmatrix }

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:

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

Uwarunkowanie zadania

Twierdzenie 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 \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$ cond ${}_2(X)||A-\tilde{A}|{|}_2.$

Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia $X$ jest ortogonalna, $X^{-1}=X^T$, to mamy cond ${}_2(X)=1$ i w konsekwencji zachodzi

Wniosek 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

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \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 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

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

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \displaystyle A = \beginpmatrix 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 \\ \endpmatrix }

{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 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 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 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}

Biblioteki

LAPACK zawiera w sobie kolekcję doskonałych narzędzi do rozwiązywania zadania własnego. {WYmienic}

ARPACK rozwiązuje zadanie własne dla macierzy rozrzedzonych, zwłaszcza symetrycznych.

Funkcja !eig! w Octave wyznacza wszystkie wartości własne (lub pary własne) zadaniej gęstej macierzy. Jak dotąd, tylko MATLAB potrafi skorzystać z ARPACKa dla wyznaczenia fragmentów widma macierzy rzadkiej, za pomocą funkcji !eigs!.