Numeryczne zadanie własne

W tym rozdziale zajmiemy się symetrycznym zadaniem własnym, tzn. zadaniem znajdowania wartości i/lub wektorów własnych dla macierzy symetrycznej $ A=A^T $. W zadaniach zbadamy, jak działa podstawowa funkcja octave'a rozwiązująca zadanie własne, tzn. funkcja octave'a eig() oraz zbadamy zbieżność różnych wersji metody potęgowej.

Metoda potęgowa jest zdefiniowana następująco: dla dowolnego wektora
$ \hat{x}_0=x_0 $ o normie drugiej równej jeden (np. losowego), definiujemy
ciąg iteracji metody potęgowej następująco:
\begin{equation}(#)
\begin{array}{lclr}
\hat{x}_k&=&Ax_{k-1} &\qquad (\mathrm{iteracja})\\
x_k&=&\frac{ \hat{x}_k}{ \|\hat{x}_k\|_2} & \qquad (\mathrm{normowanie}) \\
r_k&=& x_k^TAx_k \qquad &(\mathrm{iloraz      Rayleigha})
\end{array} \qquad k>0
\end{equation}

Zgodnie z teorią, o ile $ x_0 $ nie jest ortogonalny do
wektora własnego (podprzestrzeni własnej) dla największej wartości własnej, to $ x_{2k} $ zbiegnie
do tego wektora własnego, a $ r_k $ - iloraz Rayleigha - do tej wartości własnej.

  1. Zapoznaj się z pomocą do funkcji octave'a eig().

    Za jej pomocą oblicz wartości własne macierzy

    $$<br />
 \left(<br />
 \begin{array}{rrr}<br />
   0& 1 &1\\<br />
   1& 0 & 1\\<br />
   1 & 1 & 0<br />
 \end{array}<br />
 \right).<br />
$$

    Sprawdź, czy dla $ \lambda $ wartości własnej $ A $ prawdą jest, że
    wyznacznik $ A-\lambda*I $ wynosi zero.

    Wyznacz macierz $ V $ taką, że jej kolumny to wektory własne $ A $.

    Przetestuj w normie macierzowej Frobeniusa, czy $ \|AV-\Lambda V\|_2 $ wynosi zero.
    Tutaj $ \Lambda $ - to macierz diagonalna z wartościami własnymi $ A $.

  2. Zastosuj metodę potęgową ( [link]) do zadania własnego z macierzą $ A $ z poprzedniego zadania.
    Za warunek stopu przyjmiemy, że iloraz Rayleigha $ r_k $ praktycznie przestanie się zmieniać, co może sugerować, że jest bliski granicy, czyli odpowiedniej wartości własnej, tzn. jeśli

    $$<br />
  |r_{k+1}-r_k|<tol<br />
$$

    to zatrzymujemy iterację. Przyjmijmy, że $ tol=10^{-12} $.

    Przetestuj, czy otrzymane $ r_k $ rzeczywiście jest dobrym przybliżeniem
    największej co do modułu wartości własnej $ A $ - czyli dwa, a $ x_k $ zbiegło do wektora własnego dla tej wartości własnej o normie drugiej równej jeden, tzn. do
    $ v_1=(\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})^T $, lub $ -v_1 $.

    Przetestuj algorytm biorąc za $ x_0 $:

    1. wektory losowe - otrzymane za pomocą funkcji losowych rand() lub randn().
    2. wektor własny dla drugiej wartości własnej, tzn. dla minus jeden, czyli unormowany wektor $ v_2=(-2,1,1)^T $
    3. wektor prawie ortogonalny do $ v_1=(\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})^T $, np. wektor otrzymany przez zaburzenie ostatniej składowej wektora $ v_2 $: $ (-2,1,1.01)^T $.

    Czy we wszystkich tych przypadkach metoda zbiegnie do $ v_1 $ lub $ -v_1 $?

    Czy wybór $ x_0 $ wpływa na ilość iteracji metody potęgowej?

  3. (#) Napisz funkcję octave'a z implementacją metody potęgowej ( [link]), tzn. napisz m-plik z funkcją:
     function [x,r,it]=power(A, tol, x0)

    która za parametry ma

    • $ A $ - macierz wymiary $ N\times N $, której wektora i wartości własnej szukamy
    • $ tol $ - warunek stopu, domyślnie przyjmujący wartość $ 10^{-12} $,
    • $ x_0 $ -wektor startowy niezerowy, domyślnie losowy,

    i która zwraca

    • przybliżenie - unormowanego w normie drugiej - wektora własnego $ x $ dla największej wartości własnej co do modułu macierzy $ A $,
    • $ r $ - przybliżenie tej wartości własnej,
    • $ it $ - ilość iteracji metody potęgowej.

    Warunek stopu przyjmiemy jak w poprzednim zadaniu, tzn. że iloraz Rayleigha $ r_k $ przestanie się zmieniać. Jeśli

    $$<br />
  |r_{k+1}-r_k|<tol<br />
$$

    to zatrzymujemy iterację.

  4. Utwórz macierz symetryczną (ale nie diagonalną) $ N\times N $ o wartościach własnych $ 1,\ldots,N $.
    Sprawdź dla $ N=10,100,1000 $ czy funkcja octave'a eig() znajdzie
    wartości własne tej macierzy.

  5. (#) Utwórz macierze symetryczne $ A_a $ (niediagonalne) wymiaru $ 11\times 11 $ o wartościach własnych $ 1,\ldots,10,a $ dla $ a=11,101,1001 $ i ze znanymi wektorami własnymi.

    Sprawdź, ile iteracji będzie potrzebowała metoda potęgowa dla losowego wektora startowego do znalezienia największej wartości własnej.

    Czy potwierdzają się wyniki teoretyczne, że szybkość zbieżności zależy od stosunku największej co do modułu wartości własnej do drugiej co do modułu wartości własnej, czyli w tym zadaniu $ a/10 $?

  6. Dla macierzy z zadania [link] dla $ a=11,101 $ ze znanym wektorem własnym $ v_1 $ dla wartości własnej $ a $, przetestuj zbieżność $ x_k $ obliczanego przy pomocy
    metody potęgowej do $ v_1 $. W szczególności sprawdź, czy

    $$\|v_1-x_k\|_2=O(|\lambda_2/\lambda_1|^k)$$

    dla $ \lambda_2=10 $ i $ \lambda_1=a $.

  7. Dla macierzy z zadania [link] dla $ a=11,101 $ dla wartości własnej $ a $, przetestuj zbieżność $ r_k $ ilorazu Rayleigha obliczanego przy pomocy
    metody potęgowej do $ \lambda_1 $. W szczególności sprawdź, czy

    $$|\lambda_1-r_k|=O(|\lambda_2/\lambda_1|^{2k})$$

    dla $ \lambda_2=10 $ i $ \lambda_1=a $.

  8. Dla macierzy z zadania [link] dla $ a=11,101 $ ze znanym wektorem własnym $ v_1 $ dla wartości własnej $ 11 $ przetestuj zbieżność
    metody potęgowej dla różnych wektorów startowych:

    • dla losowego wektora $ x_0 $ ortogonalnego do $ v_1 $ (jak go można utworzyć?)
    • dla wektora prawie ortogonalnego - za $ x_0 $ bierzemy $ w+10^{-6}*z $,
      gdzie $ w^Tv_1=0 $ i $ z $ to wektor losowy o normie drugiej równej jeden.
    • dla dowolnego losowego wektora $ x_0 $

    Proszę powtórzyć testy kilka razy dla każdego podpunktu - czy zdarzyło się, że
    wyniki były istotnie różne w zależności od wylosowanego wektora?

  9. Zaimplementuj funkcję analogiczną do tej z zadania  [link] z implementacją odwrotnej metody potęgowej, tzn. metody potęgowej zastosowanej do macierzy odwrotnej $ A^{-1} $ zamiast $ A $:
     function [x,r,it]=invpower(A, tol, x0)

    Wszystkie parametry
    i zwracane wartości będą takie same jak w funkcji power() z zadania  [link].

  10. Utwórz macierz z zadania [link] dla $ a=11,101,1001 $.
    Sprawdź, ile iteracji będzie potrzebowała odwrotna metoda potęgowa dla losowego wektora startowego
    do znalezienia przybliżenia najmniejszej wartości własnej przy ustalonym warunku stopu $ |r_k-r_{k+1}|\leq 10^{-12} $. Tutaj $ r_k $ - to odpowiedni iloraz Rayleigha.

  11. Zaimplementuj funkcję analogiczną do tej z zadania  [link] z implementacją odwrotnej metody potęgowej z przesunięciem (shift), tzn. metody potęgowej zastosowanej do macierzy $ (A- b*I)^{-1} $ dla parametru $ b $:
     function [x,r,it]=invpowwiel(A, b,tol, x0)

    Tutaj $ b $ będzie dodatkowym parametrem tej funkcji, a pozostałe parametry i zwracane wartości pozostają takie same jak w funkcji z zadania  [link].

  12. Utwórz macierz z zadania [link] dla $ a=11,101,1001 $.

    Sprawdź, ile iteracji będzie potrzebowała odwrotna metoda potęgowa z przesunięciem dla losowego wektora startowego
    i różnych wartości parametru $ b=1.3, 1.5,1.8,2.1, 2.01 $ do znalezienia różnych wartości własnych macierzy $ A_a $.

  13. Zaimplementuj funkcję analogiczną do tej z zadania  [link] z implementacją metody ilorazów Rayleigha:
     function [x,r,it]=raylegh(A, b,tol, x0)

    Metoda polega na tym, że dla danego wektora startowego $ x_0 $ o drugiej normie równej jeden, obliczamy iloraz Rayleigha $ r_0=x_0^TAx_0 $ i iterujemy:
    \begin{equation}(#)
    \begin{array}{lcl}
    x_k&=&\frac{(A-r_{k-1}I)^{-1}x_{k-1}}{\| (A-r_{k-1}I)^{-1}x_{k-1}\|_2} \\
    r_k&=& x_k^TAx_k
    \end{array} \qquad k>0
    \end{equation}
    Oczywiście w implementacji proszę nie obliczać macierzy odwrotnej
    do $ A-r_k*I $, tylko rozwiązywać układ równań liniowych z tą macierzą.

    Wszystkie parametry
    i zwracane wartości będą takie same jak w funkcji power() z zadania  [link].

  14. Przetestuj metodę ilorazów Rayleigha na przykładzie macierzy symetrycznej z zadania [link] dla $ a=11,101,1001 $ biorąc losowy wektor startowy. W kolejnych testach proszę wziąć wektory startowe bliskie wektorom własnym dla wartości własnych $ 2 $ i $ 7 $.
  15. Zaimplementuj funkcję octave'a sprowadzającą macierz symetryczną
    $ A $ wymiaru $ N\times N $ do macierzy podobnej trójdiagonalnej: $ Q^TBQ=A $ lub $ B=QAQ^T $ za pomocą odpowiednich przekształceń Householdera:

     function [B,Q]=sym2tridiag(A,typ)

    która za parametry ma

    • $ A $ - symetryczną macierz wymiary $ N\times N $
    • $ tol $ - typ zwracanej macierzy $ Q $, typ=0 oznacza, że
      $ Q $ jest zwracana jako macierz $ N\times N $ - wymnożone macierze Householdera (domyślna opcja),
      typ=1 oznacza, że w macierzy $ Q $ zwracamy kolumny odpowiednich macierzy Householdera

    i która zwraca

    • $ B $ - trójdiagonalną macierz symetryczną podobną do $ A $ w formacie rzadkim octave'a (sparse),
    • $ Q $ macierz będącą iloczynem przekształceń Householdera lub macierz, w której kolumnach są odpowiednie wektory Householdera w zależności od wartości parametru typ.
  16. Przetestuj funkcję z poprzedniego zadania dla jakiejś macierzy symetrycznej $ A $ (ani diagonalnej, ani trójdiagonalnej, ani trójkątnej) o znanych wartościach własnych, np. dla macierzy z zadania [link] dla $ a=11,101,1001 $.

    Policz normę Frobeniusa $ \|A-Q^TBQ \|_F $ i $ \|QAQ^T-B \|_F $.

    Sprawdź, czy funkcja octave'a eig() dla obu macierzy zwróci te same wartości własne, i czy odpowiednie wektory własne spełniają:

    $$<br />
  Qv_k=w_k, \qquad Q^Tw_k=v_k,<br />
$$

    gdzie $ w_k $ to wektor własny macierzy $ B $ dla wartości własnej $ \lambda_k $,
    a $ v_k $ - to wektor własny dla tej samej wartości własnej $ \lambda_k $ dla macierzy $ A $.