MN13

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania


Wektory i wartości własne

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Niech będzie dana rzeczywista kwadratowa macierz wymiaru . Wektorem własnym oraz odpowiadającą mu wartością własną nazwiemy taką parę, dla której

przy czym .

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 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 , 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 pewnej silnie rozrzedzonej macierzy (gigantycznego rozmiaru, równego liczbie indeksowanych stron, czyli w chwili pisania tego tekstu około stron), odpowiadającego wartości własnej równej 1:

Współrzędne wektora 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 gwarantują, że taki wektor zawsze istnieje i jest jedyny! Co więcej, wartość 1 jest dominującą wartością własną , 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 . Zachodzi także fakt odwrotny, to znaczy miejsca zerowe wielomianu są wartościami pewnej macierzy, np. miejsca zerowe wielomianu

są wartościami własnymi m.in. macierzy stowarzyszonej,

Funkcja Octave'a compan(p) wyznacza macierz stowarzyszoną dla zadanego wielomianu o współczynnikach w wektorze . 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 wymiaru ma rozkład

gdzie jest ortogonalna (tzn. ), a jej kolumnami są wektory własne , natomiast jest diagonalna z wartościami własnymi na diagonali:

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 . Przy tym mogą być nam pomocne dwa fakty:

Fakt

Dowolna wartość własna macierzy spełnia

gdzie jest dowolną normą macierzową indukowaną przez normę wektorową.

Dowód

Rzeczywiście, skoro istnieje wektor taki, że , to stąd , więc fakt powyższy wynika już z definicji normy macierzy:

End of proof.gif

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 leżą w sumie (teoriomnogościowej) dysków na płaszczyźnie zespolonej,

Przykład: Koła Gerszgorina

Niech

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.

Wyznaczanie pojedynczej pary własnej

Jak wiemy z algebry, nawet gdy 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 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 spełniają

(to znaczy, istnieje dokładnie jedna dominująca wartość własna macierzy .

Załóżmy także, że istnieje baza złożona z wektorów własnych tej macierzy (tak jest np. dla macierzy symetrycznej na mocy twierdzenia o własnościach symetrycznego zadania własnego).

Kierunek własny jakiejś macierzy ma taką własność, że poddany działaniu przekształcenia wydłuża się razy, wobec tego dowolny wektor poddany działaniu najbardziej wydłuży się w kierunku . Iterując tę procedurę, powinniśmy dostawać w wyniku wektory, w których coraz bardziej dominuje kierunek . Formalnie, niech

wtedy

i w konsekwencji

Założenia, że istnieje dokładnie jedna dominująca wartość własna, , wynika, że wyrażenie w nawiasie dąży do i w konsekwencji wektory dążą, gdy , do kierunku wektora własnego , to znaczy wektora odpowiadającego dominującej wartości własnej (o ile tylko ).

Szybkość zbieżności metody potęgowej jest liniowa, o współczynniku zależnym od stosunku . W patologicznym przypadku, gdy , może więc okazać się, że metoda praktycznie nie jest zbieżna.

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

W praktyce nie wyznaczamy wzorem , 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 , to , a gdy , to ). Przy okazji, zauważ, że , 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, lub warunek małego residuum, , gdzie jest przybliżeniem dostępnym na -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.

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

Odwrotna metoda potęgowa

Zauważmy, że dla dowolnej macierzy kwadratowej o wartościach własnych i odpowiadających im wektorach własnych , mamy:

  • Macierz ma wartości własne oraz wektory własne ,
  • Jeśli dodatkowo jest nieosobliwa, to macierz ma wartości własne oraz wektory własne

Z połączenia tych dwóch własności wynika, że

Stwierdzenie O transformacji widma macierzy

Macierz (o ile istnieje), to ma wartości własne równe i wektory własne identyczne z .

Skoro tak, to jeśli najbliższą wartością własną jest , wówczas metoda potęgowa zastosowana do macierzy zbiegnie do . 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++;	
} 

Metoda Rayleigh (RQI)

Z własności metody potęgowej wynika, że metoda odwrotna potęgowa jest zbieżna tym szybciej, im bliżej jest przesunięcie (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 , korzystając z dotychczas wyznaczonego wektora i ilorazu Rayleigh:

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 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 będzie , tym bardziej rośnie uwarunkowanie , 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 , a tym samym złe uwarunkowanie macierzy i skończona precyzja arytmetyki pomagają 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 ).

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 operacji sprowadzić przez ortogonalne podobieństwo zadanie z macierzą gęstą (w przypadku niesymetrycznym) do zadania z macierzą Hessenberga, czyli macierzą, której element jest zerowy gdy tylko :

bądź wręcz trójdiagonalną, gdy była symetryczna.

Każdą macierz kwadratową da się sprowadzić do postaci Hessenberga sekwencją przekształceń postaci

gdzie jest pewnym przekształceniem Householdera. Rzeczywiście, niech

i oznaczmy . Możemy wziąć na początek przekształcenie Householdera takie, że , gdzie . Wtedy

To samo przekształcenie przyłożone z prawej strony zachowa pierwszą kolumnę i w efekcie nie zmieni struktury macierzy:

Dalej stosujemy tę samą metodę do podmacierzy wymiaru , itd. aż dochodzimy do macierzy Hessenberga.

Gdy wyjściowa macierz jest symetryczna, to z definicji, macierz wynikowa 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 ; wektory własne macierzy także można łatwo (jak?) odzyskać z wektorów własnych macierzy wynikowej.

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 już w postaci trójdiagonalnej, łatwo widzieć, że "prawie" rozpada się ona na dwie mniejsze macierze trójdiagonalne: dokładniej,

gdzie , są --- tak jak --- macierzami trójdiagonalnymi i symetrycznymi (jako podmacierze ), tylko o połowę mniejszego wymiaru, gdy Natomiast , więc macierz ma tylko cztery niezerowe elementy, każdy równy .

Zgodnie ze swoją nazwą, metoda "dziel i rządź" sprowadza zadanie znajdowania par własnych macierzy wymiaru 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 --- 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 umiemy rozwiązać zadanie własne tak, że znamy macierze: --- ortogonalną oraz --- diagonalną, takie, że

Wtedy łatwo widzieć, że dla łatwo wyznaczalnego wektora ,

W ten sposób zadanie własne dla oryginalnej macierzy wymiaru jest równoważne zadaniu własnemu macierzy diagonalnej zaburzonej o macierz rzędu 1.

Na szczęście łatwo pokazać, że jeśli nie jest wartością własną macierzy diagonalnej , to wartości własne macierzy

spełniają równanie

gdzie są elementami na diagonali macierzy .

Wykres dla macierzy jednowymiarowego laplasjanu rozmiaru 10. Zwróć uwagę na asymptoty pionowe tej funkcji oraz jej przedziałową monotoniczność.

W typowym przypadku będzie miała dokładnie 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 z małą stałą.

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 mają te same wartości własne, bo . Co więcej, powyższy algorytm (gdy jest nieosobliwa) w zasadzie jest równoważny hipotetycznemu algorytmowi iteracji prostej zastosowanemu nie do pojedynczego wektora, ale do 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 . Gdyby nie ortogonalizować zestawu wektorów , oczywiście dostalibyśmy w efekcie zbieżność wszystkich kolumn macierzy do tego samego wektora --- odpowiadającego dominującej wartości własnej . Zapewniając sobie ortogonalność , możemy liczyć na to, że kolejne kolumny macierzy będą dążyć do wektorów własnych odpowiadających kolejnym wartościom własnym (przy stosownych założeniach o , m.in. że wszystkie wartości własne spełniają dla ). 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 były dodatnie) mamy zależności oraz .

Tak więc, w sprzyjających warunkach, metoda QR, jako równoważna metodzie potęgowej na przestrzeni, będzie zbieżna: , gdzie 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 (a więc także ) będą liczby na diagonali .

Twierdzenie O zbieżności metody QR

Niech wartości własne spełniają oraz macierz o kolumnach złożonych z kolejnych wektorów własnych ma taką własność, że ma rozkład LU, .

Wtedy w metodzie QR ciąg macierzy jest zbieżny do macierzy diagonalnej, a ciąg ma podciąg zbieżny do macierzy trójkątnej, której elementy diagonalne są równe dla .

Powyższa wersja algorytmu QR jest mało praktyczna, m.in. jest zbieżna wolno i przy poważnych ograniczaniach na . 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 ze stałą równą około 30. Omówienie sposobów wyboru "sprytnego przesunięcia" wykracza niestety poza ramy wykładu.

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 sprowadzić wyjściową macierz symetryczną do postaci (prawie) diagonalnej:

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

Twierdzenie O zbieżności metody Jacobiego

Klasyczna metoda Jacobiego jest zbieżna co najmniej liniowo, tzn.

gdzie oznacza miarę niediagonalności na -tym kroku iteracji.

Można pokazać, że asymptotycznie (tzn. dostatecznie blisko granicy) zbieżność metody Jacobiego jest nawet kwadratowa.

Uwarunkowanie

Twierdzenie Bauera-Fike'a, o uwarunkowaniu wartości własnych

Niech będzie diagonalizowalna, to znaczy dla pewnej macierzy zachodzi

a więc (gdyż macierz po prawej stronie jest podobna do ) , są wartościami własnymi . Rozważmy macierz zaburzoną i jakąś jej wartość własną . Wtedy istnieje wartość własna macierzy taka, że

Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia jest ortogonalna, , to mamy 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 jest rzeczywista i symetryczna, to

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

Przykład

Weźmy dla uproszczenia . Wartości własne to zera wielomianu , zatem i w konsekwencji

gdy , 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 dla ujemnego parametru są zespolone!

Zachowanie się wartości własnych macierzy (z parametrem ) w otoczeniu

Bardziej spektakularny przykład pochodzi od Wilkinsona:

Przykład: Perfidny wielomian Wilkinsona

Niech

Zmiana współczynnika przy o 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.

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.

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.