MN12

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 $x$, które minimalizuje resztę,

$||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ę

Okazuje się bowiem, że jeśli np. potraktować $b$ jako dane eksperymentalne (obarczone pewnym losowym błędem pomiaru o rozkładzie normalnym), a $x$ --- parametrami zależności liniowej dla punktów pomiaru zadanych w macierzy $A$, to $x$ minimalizujący $||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 Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:[a,b]\toR} obserwujemy jej wartości $f_i$ (dokładne lub zaburzone) w punktach $t_i$, $1\le i\le m$. Funkcję tę chcielibyśmy przybliżyć inną funkcją $w$ należącą do pewnej $n$ wymiarowej przestrzeni liniowej $W$, np. przestrzeni wielomianów stopnia mniejszego niż $n$. Jakość przybliżenia mierzymy, sprawdzając, jak dokładnie spełniona jest przybliżona równość $f_i\approx w(t_i)$, dokładniej, badając tzw. błąd średniokwadratowy,

$\frac{1}{m}\sum _{i=1}^m(f_i-w(t_i){)}^2.$

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

$\sum _{i=1}^m{(f_i-\sum _{j=1}^n{c}_j{w}_j(t_i))}^2$

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

Rzeczywiście, kładąc Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle A=(a_{i,j})\inR^{m\times n}} z $a_{i,j}=w_j(t_i)$, $b=(f_i{)}_{i=1}^m$ i $x=(c_j{)}_{j=1}^n$, reszta jest równa $\Vert b-Ax{\Vert }_2^2$, a minimalizacja reszty jest oczywiście równoważna minimalizacji błędu średniokwadratowego.

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

Powyższe zadanie aproksymacji średniokwadratowej w zadanych węzłach $(x_i,y_i)$, $i=1,\dots ,m$. wielomianem stopnia co najwyżej $N$, realizuje w Octave funkcja polyfit(x,y,N). (Co dostaniemy, gdy $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 $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 $x$ modelu przy zachowaniu przybliżonego spełnienia warunków modelu, tzn. $Ax\approx b$.

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

$||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 $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 $A$ będzie daną macierzą o $m$ wierszach i $n$ kolumnach, Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle A\inR^{m\times n}} , taką, że

$m\ge n=\text{rank}(A),$

albo równoważnie, taką że jej wektory kolumny są liniowo niezależne. Niech także dany będzie wektor Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle b\inR^m} . Jasne jest, że wtedy układ równań $Ax=b$ nie zawsze ma rozwiązanie - mówimy, że układ jest nadokreślony.

Zadanie wygładzania liniowego polega na znalezieniu wektora Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle x^*\inR^n} , który minimalizuje wektor residualny (wektor reszty) $r=b-Ax$ w normie drugiej, tzn.

Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle \| b\,-\,A x^*\|_2\,=\,\min_{ x\inR^n} \| b\,-\,A x\|_2. }

Lemat

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

$A^{T}Ax=A^{T}b.$

Zauważmy, że jeśli macierz $A$ jest kwadratowa, $m=n$, to rozwiązaniem jest $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 $A^T$ przez $A$ i rozwiązać układ normalny. Zauważmy ponadto, że macierz $A^{T}A$ jest symetryczna i dodatnio określona, bo $(A^{T}A{)}^T=A^{T}A$ i dla $x\ne 0$ mamy $x^T(A^{T}A)x=(Ax{)}^T(Ax)=\Vert Ax{\Vert }_2>0$, przy czym ostatnia nierówność wynika z faktu, że kolumny macierzy $A$ są liniowo niezależne i dlatego $Ax\ne 0$. Przy mnożeniu $A^T$ przez $A$ wystarczy więc obliczyć tylko elementy na głównej przekątnej i pod nią, a do rozwiązania równania z macierzą $A^{T}A$ można zastosować algorytm Cholesky'ego-Banachiewicza. Jak łatwo się przekonać, koszt takiego algorytmu wynosi $n^2(m+n/3)$, przy czym dominuje koszt mnożenia obliczenia macierzy $A^{T}A$.

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

$A=\left(\begin{array}{cccc}1& 1& 1& 1\\ ϵ\\ & ϵ\\ & & ϵ\\ & & & ϵ\end{array}\right)$

mamy

$A^{T}A=\left(\begin{array}{cccc}1+{ϵ}^2& 1& 1& 1\\ 1& 1+{ϵ}^2& 1& 1\\ 1& 1& 1+{ϵ}^2& 1\\ 1& 1& 1& 1+{ϵ}^2\end{array}\right).$

Jeśli ${ϵ}^2<\nu$ to $f{l}_{\nu }(1+{ϵ}^2)=1$, co implikuje $\text{rank}(f{l}_{\nu }(A^{T}A))=1$, podczas, gdy $\text{rank}(f{l}_{\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 Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle w\inR^m} o normie $\Vert w{\Vert }_2=\sqrt{w^{T}w}=1$, odbicie (macierz) Householdera zdefiniowane jest jako

$H=I-2w{w}^T.$

Zauważmy, że

$Hx=x-2(w^{T}x)w,$

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

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

$H^{-1}=H=H^T.$

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

$H^2=(I-2w{w}^T{)}^2=I-4w{w}^T+4w(w^{T}w)w^T=I,$

oraz

$H^T=(I-2w{w}^T{)}^T=I-2(w^T{)}^T{w}^T=I.$

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

$\Vert Hx{\Vert }_2=\sqrt{(Hx{)}^T(Hx)}=\sqrt{x^T(H^{T}H)x}=\sqrt{x^{T}x}=\Vert x{\Vert }_2.$

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

$Hx=(I-2w{w}^T)x=\alpha e.$

Załóżmy dla uproszczenia, że $\Vert e{\Vert }_2=1$. Aby wyznaczyć $H$ zauważmy, że

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

a ponieważ $\alpha =\pm \Vert x{\Vert }_2$ i $\Vert w{\Vert }_2=1$ to

$w=\frac{x\mp \Vert x{\Vert }_{2}e}{\Vert x\mp \Vert x{\Vert }_{2}e{\Vert }_2}.$

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

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

gdzie

$u_i=\{\begin{array}{ll}{x}_1\mp \Vert x{\Vert }_2& i=1,\\ x_i& 2\le i\le m,\end{array}$

oraz

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $x$ na kierunek pierwszego wersora, w zależności od wybranego znaku przy $\Vert x{\Vert }_2$. Ustalimy ten znak na plus gdy $x_1\ge 0$ oraz na minus gdy $x_1<0$, co pozwoli na obliczenie $u_1$ i $\gamma$ z małym błędem względem w $f{l}_{\nu }$. Wtedy bowiem mamy

$u_1=\{\begin{array}{ll}{x}_1+\Vert x{\Vert }_2& x_1\ge 0,\\ x_1-\Vert x{\Vert }_2& x_1<0,\end{array}$

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

Rozkład QR

Odbić Householdera można użyć do rozkładu macierzy Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle A\inR^{m\times n}} na iloczyn ortogonalno-trójkątny.

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

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

w której pierwsza kolumna $a_1^{(1)}$ ma niezerową tylko pierwszą współrzędną. W następnym kroku wybieramy drugie przekształcenie Householdera ${\overline{H}}_2=I_{m-1}-v_2{v}_2^T/{\gamma }_2$ wymiaru $m-1$ tak, aby przeprowadzało wektor $(a_{i,2}^{(1)}{)}_{i=2}^m$ na kierunek pierwszego wersora w $R^{m-1}$. Rozszerzając Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle v_2\inR^{m-1}} do wektora Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle u_2\inR^m} przez dodanie zera jako pierwszej współrzędnej, $u_2=(0,v_2{)}^T$, otrzymujemy przekształcenie (macierz) Householdera $H_2=I_m-u_2{u}_2^T/{\gamma }_2$ w $R^m$ postaci

$H_2=\left(\begin{array}{cc}1& 0^T\\ 0& {\overline{H}}_2\end{array}\right).$

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

$H_n{H}_{n-1}\cdots H_2{H}_{1}A=R,$

gdzie Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle R\inR^{m\times n}} jest uogólnioną macierzą trójkątną górną, tzn. $r_{i,j}=0$ dla $i>j$. Stąd, podstawiając $Q=H_1{H}_2\cdots H_n$, dostajemy rozkład macierzy na iloczyn ortogonalno-trójkątny

$A=Q\cdot R.$

Rzeczywiście, macierz Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle Q\inR^{m\times m}} jest ortogonalna, bo

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $c=Q^{T}b=H_n\cdots H_2{H}_{1}b$. Rozbijając wektor $c$ na $c=(c_I,c_{II}{)}^T$, gdzie Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle c_I\inR^n} i Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle c_{II}\inR^{m-n}} , oraz macierz $R$ na

$R=\left(\begin{array}{c}{R}_I\\ 0\end{array}\right),$

gdzie Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle R_I\inR^{n\times n}} jest macierzą trójkątną górną, a $0$ jest macierzą zerową wymiaru $(m-n)\times n$, otrzymujemy

$\Vert r{\Vert }_2^2=\Vert c_I-R_{I}x{\Vert }_2^2+\Vert c_{II}{\Vert }_2^2.$

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

$x^{*}=R_I^{-1}{c}_I,$

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

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

$(H_{k}x{)}_i=x_i-s{u}_{k,i},$

gdzie $s=u_k^{T}x/{\gamma }_k$.

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

$\sum _{k=1}^{n}4(m-k+1)(n-k)\approx \frac{4}{3}{n}^3+2n^2(m-n)=2n^2(m-n/3)$

operacji arytmetycznych i $n$ pierwiastków kwadratowych. Zauważmy, że w przypadku $m=n$, a więc dla kwadratowego układu równań, koszt ten wynosi $(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 $u$ oraz czynniku skalującym $\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 $x$ przez macierz Householdera $H=I-\frac{1}{\gamma }u{u}^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ż $O(N^2)$ miejsc w pamięci (chociaż, przypomnijmy raz jeszcze, cała informacja o $H$ to tylko $O(N)$ liczb). Ponadto, mnożenie przez macierz to aż $O(N^2)$ działań arytmetycznych.

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

$Hx=(I-\frac{1}{\gamma }u{u}^T)x=x-\frac{1}{\gamma }u{u}^{T}x=x-\frac{1}{\gamma }u(u^{T}x).$

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 $N$, a więc uzyskaliśmu $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 $\mathrm{\Sigma }(A)$ będzie zbiorem wartości własnych macierzy $A^{T}A$. Definiujemy

${\text{cond}}_2(A)=\sqrt{\frac{\max \{\lambda :\lambda \in \mathrm{\Sigma }(A)\}}{\min \{\lambda :\lambda \in \mathrm{\Sigma }(A)\}}}.$

(Jeśli w mianowniku pojawiłoby się zero, kładziemy ${\text{cond}}_2(A)=+\mathrm{\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 $x$ będzie rozwiązaniem zadania najmniejszych kwadratów dla niezerowej prawej strony $b$,

$||b-Ax|{|}_2\to \min !$ i niech $\tilde{x}$ będzie rozwiązaniem zadania zaburzonego $||\tilde{b}-\tilde{A}\tilde{x}|{|}_2\to \min !,$

przy czym zakładamy, że

$\frac{||\tilde{b}-b|{|}_2}{||b|{|}_2},\frac{||\tilde{A}-A|{|}_2}{||A|{|}_2}\le ϵ,$

gdzie $ϵ$ jest dostatecznie małe.

Oznaczmy

$\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

$\frac{||\tilde{x}-x|{|}_2}{||x|{|}_2}\lesssim (\frac{2{\text{cond}}_2(A)}{\cos (\theta )}+\tan (\theta ){\text{cond}}_2^2(A))\cdot ϵ.$

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

Wniosek

W przypadku, gdy $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 $\sin (\theta )\ll 1$, rozwiązanie obliczone (kosztowniejszym) rozkładem QR będzie miało błąd na poziomie ${\text{cond}}_2(A)$, a tymczasem rozwiązanie wyznaczone z układu równań normalnych będzie obarczone błędem na poziomie ${\text{cond}}_2^2(A)>{\text{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 $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 $x$ miało minimalną normę euklidesową spośród wszystkich spełniających warunek $||b-Ax|{|}_2\to \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.