MN12

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

\rysunek{}{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 \CHECK{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

{{{3}}}

Przykład Wyznaczanie miejsc zerowych wielomianu

{{{3}}}

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 {\em 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 o symetrycznym zadaniu włanym

{{{3}}}

Twierdzenie Bauer-Fike

\lambda_j - \tilde{\lambda}

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

Wniosek Wartości własne macierzy symetrycznej są doskonale uwarunkowane

\lambda_j - \tilde{\lambda}

Z drugiej strony, dla macierzy niediagonalizowalnych, uwarunkowanie wartości własnych może być dowolnie duże, co ilustruje poniższy

Przykład

\lambda_\epsilon - \lambda_0

Bardziej spektakularny przykład pochodzi od Wilkinsona:

Przykład Perfidny wielomian Wilkinsona

{{{3}}}

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.

Jak okaże się za chwilę, czasem warto mieć ogólne rozeznanie o tym, gdzie {\em 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 \[

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

z - a_{ii}

Przykład Koła Gerszgorina

{{{3}}}

Przykład Widmo macierzy jednowymiarowego Laplasjanu

Norma daje:

Tw. Gerszgorina daje:

W rzeczywistości,

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 \link{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

$x_0$ = dowolny wektor startowy; k = 0;
while ( !stop )
{
  $y_k$ = $Ax_{k-1}$;
  $x_k$ = $y_k/||y_k||_\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.

\rysunek{}{Zasada działania metody potęgowej}

Metoda potęgowa doskonale sprawdza się, gdy macierz $A$ jest macierzą rozrzedzoną --- np. w przypadku macierzy Google'a.

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:

$x_0$ = dowolny wektor startowy; k = 0;
while ( !stop )
{
  $y_k$ = $Ax_{k-1}$;
  $x_k$ = $y_k/||y_k||_\infty$;
  k++;  
} 

$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||_\infty$;
  k++;  
} 

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

$x_0$ = dowolny wektor startowy; k = 0;
while ( !stop )
{
  $y_k$ = $Ax_{k-1}$;
  $x_k$ = $y_k/||y_k||_\infty$;
  k++;  
} 

$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||_\infty$;
  k++;  
} 

$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^TAx_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!

\rysunek{}{Secular equation}