LZNK. Rozkład QR. Metoda Householdera

W tym rozdziale zajmiemy się liniowym zadaniem najmniejszych kwadratów (LZNK).

Dla danej macierzy $ A $ wymiaru $ M\times N $ i wektora $ b $ wymiaru $ M $ chcemy znaleźć wektor $ x $ wymiaru $ N $ taki, że

$$<br />
  \|Ax-b\|_2=\min_y\|Ay-b\|_2.<br />
$$

Jeśli $ A $ jest macierzą kolumnami regularną (rząd $ A $ jest maksymalny równy $ N $), to zadanie to ma jednoznaczne rozwiązanie i nazywamy je regularnym LZNK (RLZNK).

Podstawowym algorytmem służącym jego rozwiązaniu jest metoda Householdera, czyli znalezienia rozkładu $ QR $ macierzy $ A $, gdzie $ Q $ to macierz ortogonalna - iloczyn macierzy Householdera, a $ R $ to macierz górnotrójkątna.

W tym rozdziale przetestujemy podstawowy operator octave'a służacy do rozwiązywania dowolnego LZNK, tzn. operator \. Zauważmy, że jeśli $ M=N $ i $ A $ jest macierzą kolumnami regularną, to $ A $ jest nieosobliwa i to regularne LZNK jest równoważne rozwiązaniu układu równań liniowych $ Ax=b $.

Przetestujemy w zadaniach funkcję octave'a qr() służącą znajdowaniu rozkładu $ QR $ macierzy, kilka własności macierzy (przekształceń) Householdera i
rozwiążemy kilka konkretnych LZNK.

  1. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych i LZNK w octavie.
    • Przetestuj operator octave'a \ rozwiązując RLZNK dla macierzy $ A $ ze znanym konkretnym rozwiązaniem $ x $:
      $$A=[1,1;1,-1;1,3],\qquad x=[1;2]$$

      przyjmując, że $ f=Ax $.

      Czy y=A\f jest rozwiązaniem tego układu?

      Policz normy residualną $ \|Ay-f\|_2 $ i normę błędu $ \|x-y\|_2 $.

    • Przetestuj ten operator dla nieregularnego LZNK dla macierzy $ A^T $ z $ x=[1;1;1] $ i $ f=Ax $.

      Czy y=A'\f jest rozwiązaniem tego układu?

      Policz normy residualną $ \|A^Ty-f\|_2 $ i normę błędu $ \|x-y\|_2 $.

  2. Rozkład QR w octave. Funkcje qr().
    1. Zapoznaj się z pomocą do funkcji: qr().
    2. Dla macierzy $ A=[1,1;1,-1;1,3] $ znajdź jej rozkład $ A=QR $ z pomocą funkcji qr().
    3. Sprawdź, czy uzyskana $ Q $ jest ortogonalna - policz normy Frobeniusa $ QQ^T-I $ i $ Q^TQ-I $.
    4. Sprawdź ten rozkład licząc normy macierzowe: normę drugą i normę Frobeniusa błędu $ A-QR $.
    5. Zastosuj ten rozkład do znalezienia rozwiązania
      LZNK $ Ax=f $ ze znanym rozwiązaniem, np. $ x=[1;0] $ i $ f=[1;1;1] $.
    6. Policz normę drugą wektorową pomiędzy $ x $, a wynikiem algorytmu $ w $, polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy residualne, tzn. normy drugie $ Aw-f $ oraz $ Rw-Q^Tf $.
  3. \textrm{Układ równań normalnych, a rozkład QR}.

    Rozpatrzmy macierz $ A_{2n,k} $ - pod-macierz wymiaru $ 2n\times k $ macierzy Vandermonde'a $ A_{2n,2n} $ dla $ 2n $ węzłów równoodległych na $ [0,1] $.

    LZNK z $ A_{2n,k} $ z wektorem prawej strony $ f $ równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) rozwiąż trzema sposobami:

    1. używając operator \,
    2. używając rozkład $ QR $ uzyskanym funkcją qr(),
    3. poprzez rozwiązanie układu równań normalnych:
      $$Bx=g$$

      dla

      $$B=A_{2n,k}^TA_{2n,k}, \qquad g=A_{2n,k}^Tf,$$

      tzn.
      tworzymy macierz układu równań normalnych $ B $, wektor prawej strony $ g $ układu równań normalnych, a następnie
      rozwiązujemy układ
      równań normalnych operatorem \ .

    Macierz Vandermonde'a można w octave'ie utworzyć za pomocą funkcji vander().

    Przeprowadź testy dla $ N=10,20,40,80 $ i $ k=2,4,n $.
    Porównaj

    • czas obliczeń
    • błąd - $ \|x-y\|_2 $
    • błąd residualny $ \|Ax-f\|_2 $

    dla
    $ x $ rozwiązania dokładnego LZNK, $ f $ wektora prawej strony LZNK,
    $ y $ przybliżenia rozwiązania uzyskanego daną metodą.

  4. (#)
    Rozkład QR a operator \ przy rozwiązywaniu układów równań liniowych,

    Rozpatrzmy macierz $ A_{n,n} $ Vandermonde'a dla $ n $ węzłów równoodległych na $ [0,1] $.

    Układ równań liniowych z wektorem prawej strony równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) rozwiąż dwoma sposobami:

    1. operatorem \,
    2. rozkładem QR uzyskanym funkcją qr(),

    Przetestuj dla $ N=10,20,40,80 $.
    Porównaj

    • czas obliczeń,
    • błąd $ \|x-y\|_2 $,
    • błąd rezydualny: $ \|Ax-f\|_2 $,

    dla
    $ x $ rozwiązania dokładnego tego układu równań, $ f $ wektora prawej strony układu i
    $ y $ przybliżenia rozwiązania uzyskanego daną metodą.

  5. Krzywa najlepiej pasująca do danych punktów.

    Zastosuj octave'a do znalezienia współczynników $ a,b $ krzywej najlepiej pasującej do zadanych punktów: $ (x_k,y_k) $, tzn. znajdź takie $ a,b $, że

    $$<br />
  \sum_k|ax_k^2+by_k^2-1|^2=\min_{c,d} \sum_k|cx_k^2+dy_k^2-1|^2.<br />
$$

    Za $ (x_k,y_k) $ przyjmiemy zaburzone punkty z danej elipsy

    $$y_k=\sqrt{1-4*x_k^2} + z_k,$$

    gdzie $ z_k $ to zaburzenie wylosowane z $ [0,10^{-2}] $ a $ x_k=1/k $ lub $ x_k=-1+h*k $ dla $ h=2/N $ $ k=1,...,N $.

    Czy obliczone $ a $ i $ b $ jest bliskie $ 4 $ i $ 1 $?

    W jednym oknie zaznacz punkty $ (x_k,y_k) $ plusami oraz narysuj fragmenty wykresów obu elips: pierwszej - dla $ a=4,b=1 $ i drugiej elipsy - dopasowanej do zaburzonych punktów.

    Powtórz obliczenia dla różnych zaburzeń $ z_k $.

  6. Zaprogramuj funkcję octave'a

     function y=H(x,w,nw)

    która dla danych wektorów $ \vec{x} $ i $ \vec{w} $ tego samego wymiaru $ N $ i skalaru $ nw=\|w\|_2 $
    zwróci wektor
    $ y=H_w \vec{x} $ dla

    $$H_w=I-2*\frac{1}{nw}\vec{w}\vec{w}^T$$

    czyli przekształcenia (macierzy) Householdera.

    Skalar może być parametrem opcjonalnym. Jeśli funkcja będzie wywołana z dwoma tylko parametrami, to normę $ w $ można obliczyć w tej funkcji.

    Przetestuj tę funkcję dla losowych wektorów $ \vec{x} $ i $ \vec{w} $ i sprawdź, czy

    $$<br />
  \|H_w \vec{x} \|_2=\|\vec{x} \|_2, \qquad H_w(H_wx)=x.<br />
$$
  7. Napisz ogólniejszą wersję funkcji z poprzedniego zadania, tzn.:
    funkcję:

     function Y=Hm(X,w,nw),

    gdzie
    $ X $ macierz $ N\times M $ i wtedy zwracany wynik to macierz $ Y=H_w*X $ . Pozostałe dwa parametry funkcji pozostaną bez zmian.

    Czy można zaimplementować taką funkcję bez użycia pętli?

    Sprawdź, wykorzystując tę funkcję, czy mnożenie przez macierz Householder nie zmienia norm macierzowych drugiej i Frobeniusa, tzn. czy:

    $$\|A\|_2=\|H_w*A\|_2=\|A*H_w\|_2$$

    i

    $$\|A\|_F=\|H_w*A\|_F=\|A*H_w\|_F$$

    dla losowej macierzy $ A $ i $ H_w $ macierzy Householdera dla losowego wektora $ w\not=0 $.

  8. Znajdź wektor Householdera $ \vec{w} $ taki, że odpowiednie przekształcenie Householdera
    przeprowadza dany wektor $ \vec{u}\not=0 $ na kierunek drugiego danego niezerowego wektora $ l*\vec{v}\not=0 $ dla
    $ l=\frac{\|v\|_2}{\|u\|_2} $.

    Przetestuj dla dowolnych dwóch różnych wektorów o tej samej długości, czy rzeczywiście $ H_w \vec{u}= \vec{v} $.

  9. Zastosuj metodę Householdera do rozwiązania zadania znalezienia prostej $ y=ax+b $ najlepiej przybliżającej
    $ N $ punktów $ (x_k,y_k)=(k,1+2*k+\epsilon_k) $ dla $ k=1,\ldots,N $ gdzie $ (\epsilon_1,\ldots,\epsilon_N) $ to
    losowy wektor za zakresu $ [-\epsilon,\epsilon] $, tzn.:

    $$<br />
  \sum_k|ax_k+b-y_k|^2=\min_{c,d}  \sum_k|cx_k+d-y_k|^2.<br />
$$

    Należy testować dla wartości $ \epsilon=[1,10^{-1},10^{-2},10^{-3}] $.

    Funkcja \lstinline+rand(n)+ generuje wektor losowy o rozkładzie jednostajnym na $ [0,1] $ w octavie.

    Porównaj z wynikami otrzymanymi za pomocą standardowej funkcji octave'a, tzn. \lstinline+ \ +, oraz przy wykorzystaniu funkcji octave'a \lstinline+qr(A)+.

  10. Zaprogramuj metodę Householdera rozwiązywania układu równań liniowych
    $ A\vec{x}=\vec{b} $ dla $ A $ macierzy trójdiagonalnej $ N\times N $, tzn.
    napisz funkcję octave'a:

    function [x]=hous3diag(a,b,c,f,N)

    Parametry funkcji:

    • $ a,b,c $ przekątna, pod-przekątna i nad-przekątna macierzy $ A $,
    • $ f $ - wektor prawej strony,
    • $ N $ - wymiar zadania - długość przekątnej $ a $.

    Funkcja zwraca $ x $ rozwiązanie $ Ax=f $.

    Przetestuj działanie funkcji analogicznie do zadania  [link] dla macierzy trójdiagonalnej o stałych diagonalach, np. takiej, że elementy na głównej diagonali są równe dwa, a elementy na pod- i nad-diagonalach są równe minus jeden dla $ N=10^p $ z $ p=1,2,\ldots,9 $.. Za
    wektor prawej strony $ f $ możemy przyjąć pierwszą kolumnę macierzy $ A $.