MN12LAB

Liniowe zadanie najmniejszych kwadratów

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

Oglądaj wskazówki i rozwiązania __SHOWALL__
Ukryj wskazówki i rozwiązania __HIDEALL__

Ćwiczenie: Rozszerzony układ równań dla zadania najmniejszych kwadratów

Inną, oprócz sprowadzenia do układu równań normalnych, metodą transformacji zadania najmniejszych kwadratów do zadania rozwiązywania układu równań z macierzą kwadratową, jest zapisanie w formie układu dwóch układów równań (sic!). Dokładniej, możemy scharakteryzować zadanie wygładzania jako znalezienie dwóch wektorów $\displaystyle x\in R^{n}$ oraz $\displaystyle r\in R^{m}$ takich, że

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned r &= b - Ax,\\ A^T r &= 0. \endaligned}

Zapisz macierzowo ten układ równań. Wskaż, kiedy mogłoby być opłacalne stosowanie takiego podejścia. Porównaj koszt rozwiązania tego układu wprost metodą eliminacji Gaussa, z kosztem innych metod rozwiązywania zadania najmniejszych kwadratów.

Wskazówka

Układy równań można rozwiązywać nie tylko eliminacją Gaussa...

Rozwiązanie

Macierzowo, układ zapisuje się w postaci

$\displaystyle {\begin{pmatrix}I&A\\A^{T}&0\end{pmatrix}}{\begin{pmatrix}r\\x\end{pmatrix}}={\begin{pmatrix}b\\0\end{pmatrix}}.$

Koszt rozwiązywania takiego układu równań to oczywiście $\displaystyle O((m+n)^{3})$, dużo więcej niż innych poznanych metod. Ale zalety takiego podejścia mogą objawić się, gdy macierz $\displaystyle A$ jest rozrzedzona wielkiego wymiaru (i $\displaystyle m\approx n$), bo wtedy możemy zastosować znany nam arsenał metod iteracyjnych.

Ćwiczenie

W twierdzeniu o uwarunkowaniu zadania najmniejszych kwadratów mówi się, że

$\displaystyle {\frac {||b-Ax^{*}||_{2}}{||b||_{2}}}<1.$

Wyjaśnij, dlaczego rzeczywiście tak jest.

Rozwiązanie

Rozwiązanie zadania najmniejszych kwadratów minimalizuje $\displaystyle ||b-Ax||_{2}$, tzn. dla każdego $\displaystyle x$,

$\displaystyle ||b-Ax^{*}||_{2}\leq ||b-Ax||_{2}.$

W szczególności dla $\displaystyle x=0$ dostajemy $\displaystyle ||b-Ax^{*}||_{2}\leq ||b-A\cdot 0||_{2}=||b||_{2}$. Ostra nierówność wynika z jednoznaczności: $\displaystyle A^{T}Ax^{*}=A^{T}b\neq 0$, stąd $\displaystyle x^{*}\neq 0$.

Ćwiczenie: Dopasowanie liniowych parametrów funkcji do danych

Znajdź $\displaystyle a$ i $\displaystyle b$ takie, że funkcja $\displaystyle f(x)=a+b\,e^{-x}$ minimalizuje błąd średniokwadratowy dla danych:

x	f(x)
0.00	4.00000000000000e+00
1.25	3.28650479686019e+00
2.50	3.08208499862390e+00
3.75	3.02351774585601e+00
5.00	3.00673794699909e+00
6.25	3.00193045413623e+00
7.50	0.00055308437015e+00
8.75	3.00015846132512e+00
10.00	3.00004539992976e+00

Wskazówka

Musisz zadanie wyrazić w terminach liniowego zadania najmniejszych kwadratów. Zauważ, że $\displaystyle x_{i+1}=x_{i}+1.25$.

Rozwiązanie

Ma być

$\displaystyle \sum _{i=0}^{10}|a+b\,e^{-x_{i}}-f(x_{i})|^{2}\rightarrow \min !$

A więc, macierzowo,

$\displaystyle ||{\begin{pmatrix}1&e^{-x_{0}}\\\vdots \\1&e^{-x_{10}}\end{pmatrix}}\cdot {\begin{pmatrix}a\\b\end{pmatrix}}-{\begin{pmatrix}f(x_{0})\\\vdots \\f(x_{1}0)\end{pmatrix}}||_{2}^{2}\rightarrow \min !$

A więc mamy sformułowane zadanie w języku liniowego zadania najmniejszych kwadratów. Reszta jest liczeniem...

Wyznaczone najlepsze dopasowanie naszego modelu do danych.

Jak widzisz, dane punkty pasują --- za wyjątkiem "dziwnego" $\displaystyle x=7.5$ --- do modelu $\displaystyle f^{*}(x)=3+e^{-x}$. Duże i niespodziewane zaburzenie w $\displaystyle x=7.5$ spowodowało, że dopasowanie w sensie najmniejszych kwadratów ma istotnie inne parametry od $\displaystyle f^{*}$. Sposobem zmniejszenia wpływu takiego zaburzenia na ogólny wynik może być wprowadzenie do zadania wag (relatywnie małej dla $\displaystyle x=7.5$), i minimalizacja $\displaystyle \sum _{i=0}^{10}\omega _{i}\cdot |a+b\,e^{-x_{i}}-f(x_{i})|^{2}\rightarrow \min !$. Zobacz też następne zadanie.

Ćwiczenie: Ważone zadanie najmniejszych kwadratów

Niech $\displaystyle A$ będzie macierzą $\displaystyle m\times n$ pełnego rzędu, przy czym $\displaystyle m\geq n$. Podaj algorytm rozwiązywania ważonego zadania najmniejszych kwadratów:

$\displaystyle \sum _{i=1}^{n}\omega _{i}|b_{i}-(Ax)_{i}|^{2}\rightarrow \min !,$

gdzie zakładamy, że $\displaystyle 0<\omega _{i}\leq 1$ są danymi wagami. (Gdy wszystkie $\displaystyle \omega _{i}=1$, zadanie sprowadza się do zwykłego zadania najmniejszych kwadratów.)

Rozwiązanie

Innymi słowy, szukamy $\displaystyle x$, minimalizującego

$\displaystyle \sum _{i=1}^{n}|{\sqrt {\omega _{i}}}(b_{i}-(Ax)_{i})|^{2}=\sum _{i=1}^{n}|({\widetilde {b}}_{i}-({\widetilde {A}}x)_{i})|^{2},$

gdzie $\displaystyle {\widetilde {b}}=Db$, $\displaystyle {\widetilde {A}}=DA$ i

$\displaystyle D={\begin{pmatrix}\omega _{1}&&\\&\ddots &\\&&\omega _{n}\end{pmatrix}}.$

A więc, zadanie sprowadza się do zadania najmniejszych kwadratów bez wag dla zmodyfikowanej macierzy i wektora prawej strony, $\displaystyle ||{\widetilde {b}}-{\widetilde {A}}x||_{2}\rightarrow \min !$.

Ćwiczenie

Opisz szczegółowo sposób rozwiązywania układu $\displaystyle N$ równań z $\displaystyle N$ niewiadomymi

$\displaystyle Ax=b,$

korzystający z rozkładu QR metodą Householdera.

Rozwiązanie

Stosując rozkład QR metodą Householdera do macierzy $\displaystyle A$, dostajemy w rezultacie

$\displaystyle A=QR=H_{N-1}\cdot H_{N-2}\cdots H_{1}R.$

Ponieważ $\displaystyle H_{j}^{-1}=H_{j}^{T}=H_{j}$, to $\displaystyle Q^{-1}=H_{1}\cdots H_{N-2}\cdot H_{N-1}$. Dlatego

$\displaystyle x=A^{-1}b=R^{-1}\,Q^{-1}=R^{-1}\cdot H_{1}\cdots H_{N-2}\cdot H_{N-1}b.$

Powyższą równość implementujemy korzystając z operacji mnożenia $\displaystyle H_{j}$, opisanej w poprzednim zadaniu. W szczególności pamiętamy, by nie wyznaczać pełnej macierzy $\displaystyle H_{j}$. Pseudokod procedury byłby następujący:

Algorytm Metoda rozwiązywania układu równań przy użyciu przekształceń Householdera

---

wyznacz macierze Householdera <math>\displaystyle H_1,\ldots, H_{N-1}</math> oraz macierz trójkątną <math>\displaystyle R</math>,
	okreslające rozkład QR macierzy <math>\displaystyle A</math>;

y = b;
for i = 1:N-1
	 y = <math>\displaystyle H_{N-i}</math>*y;
end
	 
rozwiąż układ z macierzą trójkątną <math>\displaystyle Rx = y</math>;

Ćwiczenie: Obroty Givensa

Innym sposobem wyzerowania wybranych elementów zadanego wektora $\displaystyle x$ za pomocą przekształceń ortogonalnych jest zastosowanie tzw. obrotów Givensa,

$\displaystyle {\begin{pmatrix}c&s\\-s&c\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}=||x||_{2}{\begin{pmatrix}1\\0\end{pmatrix}}.$

Wskaż jak dobrać $\displaystyle c$ i $\displaystyle s$ tak, by macierz

$\displaystyle G={\begin{pmatrix}c&s\\-s&c\end{pmatrix}}$

była ortogonalna i przekształcała $\displaystyle x$ w zadany wyżej sposób. Jak zastosować sekwencję obrotów Givensa tak, by zadany wektor $\displaystyle N$-wymiarowy przeprowadzić na wektor o kierunku wektora jednostkowego? Porównaj koszt tej operacji z kosztem przekształcenia Householdera. Kiedy opłaca się stosować obroty Givensa w miejsce odbić Householdera?

Rozwiązanie

Prosty rachunek pokazuje, że

$\displaystyle c={\frac {x_{1}}{||x||_{2}}},\quad s={\frac {x_{2}}{||x||_{2}}}.$

(Zauważ, że $\displaystyle c^{2}+s^{2}=1$, więc $\displaystyle G$ faktycznie można traktować jako macierz obrotu o kąt $\displaystyle \theta$ taki, że $\displaystyle c=\cos(\theta )$ i $\displaystyle s=\sin(\theta )$.)

Jak widać, występuje tu zadanie obliczania normy euklidesowej i w związku z tym ryzyko niepotrzebnego nadmiaru bądź niedomiaru. Dlatego w praktyce obliczeniowej rozpatrujemy dwa przypadki:

Algorytm Wyznaczenie obrotu Givensa

---

if ( <math>\displaystyle |x_1|</math> > <math>\displaystyle |x_2|</math> )
{
	t = <math>\displaystyle x_2</math> / <math>\displaystyle x_1</math>;
	c = 1 / sqrt(1+t*t);
	s = t * c;
}
else
{
	t = <math>\displaystyle x_1</math> / <math>\displaystyle x_2</math>;
	s = 1 / sqrt(1+t*t);
	c = t * s;
}

Chcąc obrotami Givensa wyzerować wszystkie --- z wyjątkiem pierwszej --- współrzędne danego wektora $\displaystyle N$-wymiarowego, musimy zastosować sekwencję obrotów dotyczących kolejno: pierwszej i drugiej współrzędnej, pierwszej i trzeciej, itp. Po $\displaystyle N-1$ krokach dostaniemy wektor, o który nam chodziło.

Koszt jednego obrotu Givensa to 4 działania arytmetyczne i jedno pierwiastkowanie, zatem koszt wyzerowania wszystkich $\displaystyle N-1$ (tzn. oprócz pierwszej) współrzędnych wektora jest równy $\displaystyle 4N-4$ działań arytmetycznych oraz $\displaystyle N-1$ pierwiastkowań, a więc jest wyższy niż analogicznego przekształcenia Householdera (ech, te pierwiastki!...). Istnieje jednak sprytna modyfikacja, tzw. algorytm Gentlemana, praktycznie zrównujący koszty implementacji sekwencji obrotów Givensa i odbić Householdera.

Ponadto, jest ważna klasa macierzy, dla których stosowanie obrotów Givensa jest znacznie tańsze od odbić Householdera: gdy w wektorze $\displaystyle x$ już na starcie jest wiele współrzędnych zerowych, bo wtedy wystarczy obrotami Givensa wyzerować pozostałe niezerowe współrzędne.

Takim przypadkiem jest np. konstrukcja rozkładu QR dla macierzy Hessenberga, czyli macierzy górnej trójkątnej uzupełnionej o jedną niezerową poddiagonalę --- precyzyjniej, dla takiej macierzy $\displaystyle A$, której elementy spełniają $\displaystyle a_{ij}=0$ dla $\displaystyle i-j>1$. Rzeczywiście, wtedy w każdej kolumnie wystarczy wyzerować tylko jeden element! Zadanie znalezienia rozkładu QR niedużej i prawie-kwadratowej macierzy Hessenberga jest częścią składową metody GMRES iteracyjnego rozwiązywania wielkich układów równań liniowych z macierzą niesymetryczną.