MN12

From Studia Informatyczne


Spis treści

Nadokreślone układy równań liniowych

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

Zajmiemy się zadaniem wygładzania liniowego, nazywanym też liniowym zadaniem najmniejszych kwadratów. Jest ono uogólnieniem zadania rozwiązywania kwadratowych układów równań liniowych do przypadku, gdy układ jest nadokreślony --- to znaczy, jest więcej równań niż niewiadomych. W takim przypadku nie należy liczyć na to, że uda się nam wskazać rozwiązanie spełniające wszystkie równania (jest ich za dużo!), dlatego będziemy szukać rozwiązania \displaystyle x, które minimalizuje resztę,

\displaystyle ||b-Ax||_2.

Jest to praktycznie bardzo często pojawiające się zadanie, a autorem pierwszego rozwiązania był nie kto inny jak sam wielki Gauss.

Carl Friedrich Gauss  Zobacz biografię
Enlarge
Carl Friedrich Gauss
Zobacz biografię

Okazuje się bowiem, że jeśli np. potraktować \displaystyle b jako dane eksperymentalne (obarczone pewnym losowym błędem pomiaru o rozkładzie normalnym), a \displaystyle x --- parametrami zależności liniowej dla punktów pomiaru zadanych w macierzy \displaystyle A, to \displaystyle x minimalizujący \displaystyle ||b-Ax||_2 (właśnie w tej normie!) jest jednocześnie najbardziej prawdopodobnym zestawem współczynników tej zależności. W języku statystyki takie zadanie nazywa się zadaniem regresji liniowej i jest w tym kontekście bardzo często znajdowane w najrozmaitszych gałęziach nauki --- wszędzie tam, gdzie zachodzi potrzeba dopasowania parametrów liniowego modelu do wyników uzyskanych na drodze eksperymentu.

Stąd zresztą nazwa zadania: wygładzanie liniowe, bo chodzi nam o to, by dopasowując parametry krzywej do wyników eksperymentu, wygładzić ewentualne błędy pomiarowe.

Dopasowanie krzywej minimalizującej błąd średniokwadratowy

Przykład

Przypuśćmy, że dla pewnej funkcji \displaystyle f:[a,b]\toR obserwujemy jej wartości \displaystyle f_i (dokładne lub zaburzone) w punktach \displaystyle t_i, \displaystyle 1\le i\le m. Funkcję tę chcielibyśmy przybliżyć inną funkcją \displaystyle w należącą do pewnej \displaystyle n wymiarowej przestrzeni liniowej \displaystyle W, np. przestrzeni wielomianów stopnia mniejszego niż \displaystyle n. Jakość przybliżenia mierzymy, sprawdzając, jak dokładnie spełniona jest przybliżona równość \displaystyle f_i \approx w(t_i), dokładniej, badając tzw. błąd średniokwadratowy,

\displaystyle    \frac{1}{m}\sum_{i=1}^m (f_i-w(t_i))^2.

Wybierając pewną bazę \displaystyle (w_j)_{j=1}^n w \displaystyle W i rozwijając \displaystyle w w tej bazie, \displaystyle w(t)=\sum_{j=1}^n c_jw_j(t), sprowadzamy problem do minimalizacji

\displaystyle \sum_{i=1}^m\left(f_i-\sum_{j=1}^n c_jw_j(t_i)\right)^2

względem \displaystyle c_j, a więc do zadania wygładzania liniowego.

Rzeczywiście, kładąc \displaystyle A=(a_{i,j})\inR^{m\times n} z \displaystyle a_{i,j}=w_j(t_i), \displaystyle  b=(f_i)_{i=1}^m i \displaystyle  x=(c_j)_{j=1}^n, reszta jest równa \displaystyle \| b-A x\|_2^2, a minimalizacja reszty jest oczywiście równoważna minimalizacji błędu średniokwadratowego.

Wielomian  (czerwony) stopnia 3, aproksymujący 7 zadanych wartości (zaznaczone na zielono) danej funkcji  w sensie minimalizacji błędu średniokwadratowego
Enlarge
Wielomian \displaystyle w (czerwony) stopnia 3, aproksymujący 7 zadanych wartości (zaznaczone na zielono) danej funkcji \displaystyle f w sensie minimalizacji błędu średniokwadratowego

Powyższe zadanie aproksymacji średniokwadratowej w zadanych węzłach \displaystyle (x_i,y_i), \displaystyle i=1,\ldots,m. wielomianem stopnia co najwyżej \displaystyle N, realizuje w Octave funkcja polyfit(x,y,N). (Co dostaniemy, gdy \displaystyle N=m-1?)

Można pokazać, że rozwiązanie minimalizujące błąd średniokwadratowy jest najbardziej prawdopodobnym zestawem parametrów naszego (liniowego) modelu, gdy zmierzone wartości \displaystyle f_i mogą być zaburzone losowym błędem pomiarowym.

W kontekście nie-statystycznym, możemy myśleć o zadaniu wygładzania liniowego jako sposobie skrócenia listy parametrów \displaystyle x modelu przy zachowaniu przybliżonego spełnienia warunków modelu, tzn. \displaystyle Ax\approx b.

Dodajmy, że spotyka się uogólnienie tego zadania w formie następującej: dla danych wartości \displaystyle b\in R^m, i danej funkcji \displaystyle F:R^n\rightarrow R^m, znaleźć \displaystyle x\in R^n minimalizujący resztę:

\displaystyle ||b-F(x)||_2.

Właśnie tego typu nieliniowe zadanie najmniejszych kwadratów rozwiązują np. nasze przenośne [ odbiorniki GPS]... Na marginesie zauważmy, że gdy \displaystyle F jest liniowa, zadanie sprowadza się do poprzedniego. W niniejszym wykładzie ograniczymy się wyłącznie do liniowego zadania najmniejszych kwadratów, nieliniowe jest omówiane na wykładzie z metod optymalizacji.

Układ równań normalnych

Niech \displaystyle A będzie daną macierzą o \displaystyle m wierszach i \displaystyle n kolumnach, \displaystyle A\inR^{m\times n}, taką, że

\displaystyle m\,\ge\,n\,=\, \mbox{rank} (A),

albo równoważnie, taką że jej wektory kolumny są liniowo niezależne. Niech także dany będzie wektor \displaystyle  b\inR^m. Jasne jest, że wtedy układ równań \displaystyle A x= b nie zawsze ma rozwiązanie - mówimy, że układ jest nadokreślony.

Zadanie wygładzania liniowego polega na znalezieniu wektora \displaystyle  x^*\inR^n, który minimalizuje wektor residualny (wektor reszty) \displaystyle  r= b-A x w normie drugiej, tzn.

\displaystyle \| b\,-\,A x^*\|_2\,=\,\min_{ x\inR^n}     \| b\,-\,A x\|_2.

Lemat

Zadanie wygładzania liniowego ma jednoznaczne rozwiązanie \displaystyle  x^*, które można scharakteryzować jako rozwiązanie układu równań

\displaystyle  A^TA x\,=\,A^T\, b.

Zauważmy, że jeśli macierz \displaystyle A jest kwadratowa, \displaystyle m=n, to rozwiązaniem jest \displaystyle  x^*=A^{-1} b i residuum jest zerem. Zadanie wygładzania liniowego jest więc uogólnieniem rozwiązywania kwadratowych układów równań liniowych.

Równanie powyższe nazywa się układem równań normalnych. Może ono nam sugerować sposób rozwiązania zadania wygładzania liniowego. Wystarczy bowiem pomnożyć macierz \displaystyle A^T przez \displaystyle A i rozwiązać układ normalny. Zauważmy ponadto, że macierz \displaystyle A^TA jest symetryczna i dodatnio określona, bo \displaystyle (A^TA)^T=A^TA i dla \displaystyle  x\ne 0 mamy \displaystyle  x^T(A^TA) x=(A x)^T(A x)=\|A x\|_2>0, przy czym ostatnia nierówność wynika z faktu, że kolumny macierzy \displaystyle A są liniowo niezależne i dlatego \displaystyle A x\ne 0. Przy mnożeniu \displaystyle A^T przez \displaystyle A wystarczy więc obliczyć tylko elementy na głównej przekątnej i pod nią, a do rozwiązania równania z macierzą \displaystyle A^TA można zastosować algorytm Cholesky'ego-Banachiewicza. Jak łatwo się przekonać, koszt takiego algorytmu wynosi \displaystyle n^2(m+n/3), przy czym dominuje koszt mnożenia obliczenia macierzy \displaystyle A^TA.

Ma on jednak pewne wady. Mnożenie macierzy powoduje w \displaystyle fl_\nu powstanie po drodze dodatkowych błędów, które mogą nawet zmienić rząd macierzy. Na przykład, dla macierzy

\displaystyle A\,=\,\left(\begin{array} {cccc}      1  &  1  &  1  &  1  \\   \epsilon  \\         &\epsilon \\         &     &\epsilon \\         &     &      &\epsilon \end{array} \right)

mamy

\displaystyle A^TA\,=\,\left(\begin{array} {cccc}      1+\epsilon^2 & 1 & 1 & 1 \\      1 & 1+\epsilon^2 & 1 & 1 \\      1 & 1 & 1+\epsilon^2 & 1 \\      1 & 1 & 1 & 1+\epsilon^2 \end{array} \right).

Jeśli \displaystyle \epsilon^2<\nu to \displaystyle fl_\nu(1+\epsilon^2)=1, co implikuje \displaystyle  \mbox{rank} (fl_\nu(A^TA))=1, podczas, gdy \displaystyle  \mbox{rank} (fl_\nu(A))=4. Inne potencjalne wady układu równań normalnych wymieniamy w dalszej części wykładu.

Poniżej przedstawimy inną metodę rozwiązywania zadania wygładzania liniowego, która oparta jest na specjalnych przekształceniach zwanych odbiciami Householdera.

Odbicia Householdera

Dla danego wektora \displaystyle  w\inR^m o normie \displaystyle \| w\|_2=\sqrt{ w^T w}=1, odbicie (macierz) Householdera zdefiniowane jest jako

\displaystyle H\,=\,I\,-\,2 w w^T.

Zauważmy, że

\displaystyle H x\,=\, x\,-\,2( w^T x) w,

a ponieważ \displaystyle ( w^T x) w=( x, w)_2 w jest rzutem prostopadłym \displaystyle  x na kierunek wektora \displaystyle  w (\displaystyle (\cdot,\cdot)_2 oznacza iloczyn skalarny), to \displaystyle H x jest odbiciem lustrzanym wektora \displaystyle  x względem hiperpłaszczyzny (wymiaru \displaystyle m-1) prostopadłej do \displaystyle  w.

Odbicia Householdera są przekształceniami nieosobliwymi spełniającymi

\displaystyle H^{-1}\,=\,H\,=\,H^T.

Rzeczywiście, ponieważ \displaystyle  w ma normę jednostkową, mamy

\displaystyle H^2 \,=\, (I-2 w w^T)^2\,=\,    I-4 w w^T+4 w( w^T w) w^T \,=\, I,

oraz

\displaystyle H^T\,=\,(I-2 w w^T)^T\,=\,I-2( w^T)^T w^T\,=\,I.

W szczególności \displaystyle H jest więc przekształceniem ortogonalnym, \displaystyle H^{-1}=H^T, czyli nie zmienia długości wektora,

\displaystyle \|H x\|_2\,=\,\sqrt{(H x)^T(H x)}\,=\,     \sqrt{ x^T(H^TH) x}\,=\,\sqrt{ x^T x}\,=\,     \| x\|_2.

Odbicia Householdera zastosujemy do przeprowadzenia danego wektora \displaystyle  x\ne 0 na kierunek innego niezerowego wektora, powiedzmy \displaystyle  e, tzn.

\displaystyle H x\,=\,(I-2 w w^T) x\,=\,\alpha\, e.


Załóżmy dla uproszczenia, że \displaystyle \| e\|_2=1. Aby wyznaczyć \displaystyle H zauważmy, że

\displaystyle w\,=\,\frac{ x-\alpha e}{2( w^T x)},

a ponieważ \displaystyle \alpha=\pm\| x\|_2 i \displaystyle \| w\|_2=1 to

\displaystyle w\,=\,\frac{ x\mp\| x\|_2 e}                 {\| x\mp\| x\|_2 e\|_2}.

W szczególności, jeśli \displaystyle  e= e_1 jest pierwszym wersorem, powyższe wzory dają

\displaystyle H\,=\,I\,-\,\frac{ u u^T}{\gamma},

gdzie

\displaystyle u_i\,=\,\left\{\begin{array} {ll}         x_1\mp\| x\|_2  &\quad i=1, \\                        x_i  &\quad 2\le i\le m,           \end{array} \right.

oraz

\displaystyle \aligned \gamma &= \frac 12\| u\|_2^2\,=\,       \frac 1 2\Big((x_1\mp\| x\|_2)^2+\sum_{i=2}^m x_i^2\Big) \\   &= \frac 1 2 \Big(\sum_{i=1}^m x_i^2\,+\,\| x\|_2^2\,\mp\,        2 x_1\|x\|_2\Big) \,=\,\|x\|_2^2\,\mp\,x_1 \|x\|_2.  \endaligned

Otrzymaliśmy dwa odbicia Householdera przekształcające dany wektor \displaystyle  x na kierunek pierwszego wersora, w zależności od wybranego znaku przy \displaystyle \| x\|_2. Ustalimy ten znak na plus gdy \displaystyle x_1\ge 0 oraz na minus gdy \displaystyle x_1<0, co pozwoli na obliczenie \displaystyle u_1 i \displaystyle \gamma z małym błędem względem w \displaystyle fl_\nu. Wtedy bowiem mamy

\displaystyle u_1\,=\,\left\{\begin{array} {ll}        x_1+\|x\|_2 & \quad x_1\ge 0, \\        x_1-\|x\|_2 & \quad x_1<0, \end{array} \right.

oraz \displaystyle \gamma=\| x\|_2^2+|x_1|\,\| x\|_2, czyli zawsze dodajemy liczby tych samych znaków. Ponadto pierwsza współrzędna wektora \displaystyle H x jest równa \displaystyle -\| x\|_2, gdy \displaystyle x_1\ge 0, a \displaystyle +\| x\|_2 jeśli \displaystyle x_1<0.

Rozkład QR

Odbić Householdera można użyć do rozkładu macierzy \displaystyle A\inR^{m\times n} na iloczyn ortogonalno-trójkątny.

Niech \displaystyle A=( a_1, a_2,\ldots, a_n), gdzie \displaystyle  a_j są wektorami-kolumnami macierzy \displaystyle A. Wybierzmy pierwsze odbicie Householdera \displaystyle H_1=I_m- u_1 u_1^T/\gamma_1 tak, aby przekształcało pierwszy wektor-kolumnę macierzy \displaystyle A na kierunek \displaystyle  e_1. Efektem pomnożenia macierzy \displaystyle A z lewej strony przez \displaystyle H_1 będzie wtedy macierz

\displaystyle A^{(1)}\,=\,( a^{(1)}_1,\ldots, a^{(1)}_n)    \,=\,(H_1 a_1,\ldots, H_1 a_n),

w której pierwsza kolumna \displaystyle  a^{(1)}_1 ma niezerową tylko pierwszą współrzędną. W następnym kroku wybieramy drugie przekształcenie Householdera \displaystyle \bar H_2=I_{m-1}- v_2 v_2^T/\gamma_2 wymiaru \displaystyle m-1 tak, aby przeprowadzało wektor \displaystyle (a^{(1)}_{i,2})_{i=2}^m na kierunek pierwszego wersora w \displaystyle R^{m-1}. Rozszerzając \displaystyle  v_2\inR^{m-1} do wektora \displaystyle  u_2\inR^m przez dodanie zera jako pierwszej współrzędnej, \displaystyle u_2=(0, v_2)^T, otrzymujemy przekształcenie (macierz) Householdera \displaystyle H_2=I_m- u_2 u_2^T/\gamma_2 w \displaystyle R^m postaci

\displaystyle H_2\,=\,\left(\begin{array} {cccc}      1 &  0^T \\       0 & \bar H_2  \end{array} \right).

Pomnożenie macierzy \displaystyle A^{(1)} z lewej strony przez \displaystyle H_2 spowoduje teraz wyzerowanie drugiej kolumny macierzy pod elementem \displaystyle a^{(1)}_{2,2}, przy czym pierwszy wiersz i pierwsza kolumna pozostaną niezmienione. Postępując tak dalej \displaystyle n razy (albo \displaystyle n-1 razy gdy \displaystyle m=n) otrzymujemy

\displaystyle H_nH_{n-1}\cdots H_2H_1A\,=\,R,

gdzie \displaystyle R\inR^{m\times n} jest uogólnioną macierzą trójkątną górną, tzn. \displaystyle r_{i,j}=0 dla \displaystyle i>j. Stąd, podstawiając \displaystyle Q=H_1H_2\cdots H_n, dostajemy rozkład macierzy na iloczyn ortogonalno-trójkątny

\displaystyle    A\,=\,Q\cdot R.

Rzeczywiście, macierz \displaystyle Q\inR^{m\times m} jest ortogonalna, bo

\displaystyle \aligned Q^{-1} &= (H_1H_2\cdots H_n)^{-1}\,=\,      H_n^{-1}\cdots H_2^{-1}H_1^{-1} \\    &= H_n^T\cdots H_2^TH_1^T \,=\,      (H_1H_2\cdots H_n)^T\,=\,Q^T. \endaligned

Dyspunując rozkładem QR, zadanie wygładzania liniowego można rozwiązać następująco. Ponieważ mnożenie przez macierz ortogonalną nie zmienia normy drugiej wektora, mamy

\displaystyle \aligned \| r\|_2 &= \| b-A x\|_2\;=\;\| b-QR x\|_2 \\     &= \|Q(Q^T b-R x)\|_2 \;=\;\| c-R x\|_2, \endaligned

gdzie \displaystyle  c=Q^T b=H_n\cdots H_2H_1 b. Rozbijając wektor \displaystyle  c na \displaystyle  c=( c_I, c_{II})^T, gdzie \displaystyle  c_I\inR^n i \displaystyle  c_{II}\inR^{m-n}, oraz macierz \displaystyle R na

\displaystyle R\,=\,\left(\begin{array} {c} R_I \\ 0\end{array} \right),

gdzie \displaystyle R_I\inR^{n\times n} jest macierzą trójkątną górną, a \displaystyle 0 jest macierzą zerową wymiaru \displaystyle (m-n)\times n, otrzymujemy

\displaystyle \| r\|_2^2\;=\;\| c_I-R_I x\|_2^2\,+\,      \| c_{II}\|_2^2.

Rozwiązanie \displaystyle  x^* zadania wygładzania jest więc rozwiązaniem układu liniowego trójkątnego,

\displaystyle x^*\,=\,R_I^{-1} c_I,

oraz \displaystyle \| r^*\|_2=\| b-A x^*\|_2=\| c_{II}\|_2.

Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde z kolejnych przekształceń Householdera \displaystyle H_k wyznaczamy przez obliczenie \displaystyle \gamma_k oraz współrzędnych wektora \displaystyle  u_k. Wektor ten ma tylko \displaystyle m-k+1 współrzędnych niezerowych, a ponadto \displaystyle u_{k,i}=a^{(k-1)}_{i,k} dla \displaystyle k+1\le i\le m. Dzięki takiej reprezentacji \displaystyle H_k, mnożenia \displaystyle H_k x możemy dla dowolnego \displaystyle  x realizować według wzoru

\displaystyle (H_k x)_i\,=\,x_i\,-\,s\,u_{k,i},

gdzie \displaystyle s= u_k^T x/\gamma_k.

Uwzględnizjąc obecność zerowych elementów w \displaystyle  u_k, przejście od macierzy \displaystyle A^{(k-1)} do \displaystyle A^{(k)} kosztuje rzędu \displaystyle 4(m-k+1)(n-k) operacji arytmetycznych i obliczenie jednego pierwiastka kwadratowego. Cały rozkład \displaystyle A=QR kosztuje więc rzędu (dla dużych \displaystyle m i \displaystyle n)

\displaystyle \sum_{k=1}^n 4(m-k+1)(n-k)\,\approx\,\frac 43n^3+2n^2(m-n)    \,=\,2n^2(m-n/3)

operacji arytmetycznych i \displaystyle n pierwiastków kwadratowych. Zauważmy, że w przypadku \displaystyle m=n, a więc dla kwadratowego układu równań, koszt ten wynosi \displaystyle (4/3)n^3 i jest dwa razy większy od kosztu eliminacji Gaussa.

Implementacja

Cała informacja o przekształceniu Householdera znajduje się w wektorze \displaystyle u oraz czynniku skalującym \displaystyle \gamma --- i w ten sposób najwygodniej przechowywać macierz Householdera. W żadnym miejscu algorytmu nie będzie nam potrzebne nic ponad umiejętność mnożenia zadanego wektora \displaystyle x przez macierz Householdera \displaystyle H = I - \frac{1}{\gamma}uu^T.

Nie popełnijmy jednak częstego błędu, prostodusznie implementując to mnożenie (przykładowo, w Octave) jako

H = eye(length(u)) - (u*u') / <math>\displaystyle \gamma</math>;
y = H*x;

Gdybyśmy użyli takiej implementacji, potrzebowalibyśmy aż \displaystyle O(N^2) miejsc w pamięci (chociaż, przypomnijmy raz jeszcze, cała informacja o \displaystyle H to tylko \displaystyle O(N) liczb). Ponadto, mnożenie przez macierz to aż \displaystyle O(N^2) działań arytmetycznych.

Aby znacznie lepiej skorzystać z bardzo specyficznej postaci macierzy \displaystyle H, która jest po prostu zaburzeniem macierzy identyczności macierzą rzędu co najwyżej 1, wystarczy w odpowiednim miejscu wstawić nawiasy:

\displaystyle  Hx = \left(I - \frac{1}{\gamma}uu^T\right) \, x = x - \frac{1}{\gamma}uu^Tx =  x - \frac{1}{\gamma}u(u^Tx).

Stąd prawidłowa implementacja mnożenia przez macierz Householdera:

<math>\displaystyle \omega</math> = u'*x;
y = x - <math>\displaystyle \frac{\omega}{\gamma}</math>*u;

Tym razem wcale nie potrzeba dodatkowej pamięci, a koszt algorytmu jest liniowy(!) względem \displaystyle N, a więc uzyskaliśmu \displaystyle N-krotne przyspieszenie w porównaniu z poprzednim!

Jest to całkiem typowe w numeryce:

Optymalizacja kodu źródłowego może być źródłem dużego przyspieszenia programu numerycznego. Ale największe przyspieszenie zazwyczaj jest efektem restrukturyzacji całego algorytmu (lub wręcz jego zmiany).

Uwarunkowanie

Łatwo domyślać się, że uwarunkowanie zadania wygładzania będzie miało jakieś cechy podobieństwa do uwarunkowania zadania rozwiązywania układu równań liniowych. Ale są także różnice, gdyż, w przeciwieństwie do układu równań liniowych, wrażliwość rozwiązania na zaburzenia będzie zależna nie tylko od samej macierzy układu, ale także od prawej strony.

Najpierw jednak musimy rozszerzyć pojęcie uwarunkowania macierzy na macierze prostokątne.

Definicja Uwarunkowanie macierzy prostokątnej w normie euklidesowej

Niech \displaystyle \Sigma(A) będzie zbiorem wartości własnych macierzy \displaystyle A^TA. Definiujemy

\displaystyle  \mbox{cond} _2(A) = \sqrt{\frac{\max\{\lambda: \lambda \in \Sigma(A)\}}{\min\{\lambda: \lambda \in \Sigma(A)\}}}.

(Jeśli w mianowniku pojawiłoby się zero, kładziemy \displaystyle  \mbox{cond} _2(A) = +\infty).

Zauważmy, że jest to rozszerzenie definicji zgodne z tym, co wcześniej definiowaliśmy dla macierzy kwadratowych.

Twierdzenie O uwarunkowaniu zadania wygładzania liniowego

Niech \displaystyle x będzie rozwiązaniem zadania najmniejszych kwadratów dla niezerowej prawej strony \displaystyle b,

\displaystyle  ||b-Ax||_2\rightarrow \min{} !
i niech \displaystyle \widetilde{x} będzie rozwiązaniem zadania zaburzonego
\displaystyle  ||\widetilde{b}-\widetilde{A}\widetilde{x}||_2\rightarrow \min{} !,

przy czym zakładamy, że

\displaystyle  \frac{||\widetilde{b}-b||_2}{||b||_2}, \quad \frac{||\widetilde{A}-A||_2}{||A||_2} \leq \epsilon,

gdzie \displaystyle \epsilon jest dostatecznie małe.

Oznaczmy

\displaystyle \sin(\theta) = \frac{||b-Ax||_2}{||b||_2} < 1

--- będzie to miara, jak bardzo jesteśmy w stanie zminimalizować resztę w oryginalnym zadaniu.

Wtedy

\displaystyle \frac{||\widetilde{x}-x||_2}{||x||_2} \lesssim \left( \frac{2 \mbox{cond} _2(A)}{\cos(\theta)} + \tan(\theta) \mbox{cond} _2^2(A)\right) \cdot \epsilon.

Generalnie więc, jeśli reszta \displaystyle ||b-Ax||_2 jest mała, wrażliwość na zaburzenia jest na poziomie \displaystyle  \mbox{cond} (A). Ale jeśli reszta jest duża (tzn. prawa strona jest taka, że nie można dobrze spełnić równania \displaystyle b\approx Ax w sensie średniokwadratowym), wtedy wrażliwość może być daleko większa.

Wniosek

W przypadku, gdy \displaystyle m \gg n, zdawać by się mogło --- zgodnie z popularnym, acz błędnym, jak za chwilę się okaże, poglądem --- że użycie układu równań normalnych jest najszybszym algorytmem, a skoro tak, to powinno dawać najmniejszą "akumulację błędu zaokrągleń". Tymczasem widzimy, że jest sens rozwiązywać nasze zadanie poprzez układ równań normalnych tylko wtedy, gdy reszta rozwiązania jest duża. W przeciwnym wypadku, gdy \displaystyle \sin(\theta) \ll 1, rozwiązanie obliczone (kosztowniejszym) rozkładem QR będzie miało błąd na poziomie \displaystyle  \mbox{cond} _2(A), a tymczasem rozwiązanie wyznaczone z układu równań normalnych będzie obarczone błędem na poziomie \displaystyle  \mbox{cond} _2^2(A) >  \mbox{cond} _2(A).

Biblioteki

W Octave, zadanie najmniejszych kwadratów rozwiązujemy praktycznie tak samo, jak równanie liniowe:

x = A \ b;

Dla zadania najmniejszych kwadratów mamy dwie podstawowe funkcje LAPACKa: DGELS, która rozwiązuje dokładnie zadanie takie, jak postawiliśmy w wykładzie, to znaczy w przypadku, gdy macierz \displaystyle A jest pełnego rzędu --- wykorzystując rozkład QR, który omówiliśmy.

Natomiast dla przypadku, gdy macierz nie jest pełnego rzędu, działa funkcja DGELSS. Wówczas, co łatwo sprawdzić, zadanie najmniejszych kwadratów tak, jak je postawiliśmy, nie musi mieć jednoznacznego rozwiązania. Jednak jeśli dołożyć wymaganie, by znalezione rozwiązanie \displaystyle x miało minimalną normę euklidesową spośród wszystkich spełniających warunek \displaystyle ||b-Ax||_2 \rightarrow \min !, to wtedy takie rozwiązanie jest już jedyne. Jednakże dla takiego zadania rozkład QR jest już niewystarczający i stosuje się inny rozkład, tzw. SVD, który wykracza poza ramy naszego wykładu.

Funkcje biblioteczne rozwiązujące zadanie wygładzania liniowego są oczywistym składnikiem wszystkich szanujących się pakietów statystycznych.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 5.3 w

  • D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Bardzo dużo na temat rozwiązywania liniowego zadania najmniejszych kwadratów można dowiedzieć się z książki

  • A. Kiełbasiński, H. Schwetlick, Numeryczna algebra liniowa, Wydawnictwa Naukowo-Techniczne, Warszawa, 1992.