MN12

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania


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

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ć jako dane eksperymentalne (obarczone pewnym losowym błędem pomiaru o rozkładzie normalnym), a --- parametrami zależności liniowej dla punktów pomiaru zadanych w macierzy , to minimalizujący (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 obserwujemy jej wartości (dokładne lub zaburzone) w punktach , . Funkcję tę chcielibyśmy przybliżyć inną funkcją należącą do pewnej wymiarowej przestrzeni liniowej , np. przestrzeni wielomianów stopnia mniejszego niż . Jakość przybliżenia mierzymy, sprawdzając, jak dokładnie spełniona jest przybliżona równość , dokładniej, badając tzw. błąd średniokwadratowy,

Wybierając pewną bazę w i rozwijając w tej bazie, , sprowadzamy problem do minimalizacji

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

Rzeczywiście, kładąc z , i , reszta jest równa , 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

Powyższe zadanie aproksymacji średniokwadratowej w zadanych węzłach , . wielomianem stopnia co najwyżej , realizuje w Octave funkcja polyfit(x,y,N). (Co dostaniemy, gdy ?)

Można pokazać, że rozwiązanie minimalizujące błąd średniokwadratowy jest najbardziej prawdopodobnym zestawem parametrów naszego (liniowego) modelu, gdy zmierzone wartości 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 modelu przy zachowaniu przybliżonego spełnienia warunków modelu, tzn. .

Dodajmy, że spotyka się uogólnienie tego zadania w formie następującej: dla danych wartości , i danej funkcji , znaleźć minimalizujący resztę:

Właśnie tego typu nieliniowe zadanie najmniejszych kwadratów rozwiązują np. nasze przenośne [ odbiorniki GPS]... Na marginesie zauważmy, że gdy 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 będzie daną macierzą o wierszach i kolumnach, , taką, że

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

Zadanie wygładzania liniowego polega na znalezieniu wektora , który minimalizuje wektor residualny (wektor reszty) w normie drugiej, tzn.

Lemat

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

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

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

mamy

Jeśli to , co implikuje , podczas, gdy . 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 o normie , odbicie (macierz) Householdera zdefiniowane jest jako

Zauważmy, że

a ponieważ jest rzutem prostopadłym na kierunek wektora ( oznacza iloczyn skalarny), to jest odbiciem lustrzanym wektora względem hiperpłaszczyzny (wymiaru ) prostopadłej do .

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

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

oraz

W szczególności jest więc przekształceniem ortogonalnym, , czyli nie zmienia długości wektora,

Odbicia Householdera zastosujemy do przeprowadzenia danego wektora na kierunek innego niezerowego wektora, powiedzmy , tzn.


Załóżmy dla uproszczenia, że . Aby wyznaczyć zauważmy, że

a ponieważ i to

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

gdzie

oraz

Otrzymaliśmy dwa odbicia Householdera przekształcające dany wektor na kierunek pierwszego wersora, w zależności od wybranego znaku przy . Ustalimy ten znak na plus gdy oraz na minus gdy , co pozwoli na obliczenie i z małym błędem względem w . Wtedy bowiem mamy

oraz , czyli zawsze dodajemy liczby tych samych znaków. Ponadto pierwsza współrzędna wektora jest równa , gdy , a jeśli .

Rozkład QR

Odbić Householdera można użyć do rozkładu macierzy na iloczyn ortogonalno-trójkątny.

Niech , gdzie są wektorami-kolumnami macierzy . Wybierzmy pierwsze odbicie Householdera tak, aby przekształcało pierwszy wektor-kolumnę macierzy na kierunek . Efektem pomnożenia macierzy z lewej strony przez będzie wtedy macierz

w której pierwsza kolumna ma niezerową tylko pierwszą współrzędną. W następnym kroku wybieramy drugie przekształcenie Householdera wymiaru tak, aby przeprowadzało wektor na kierunek pierwszego wersora w . Rozszerzając do wektora przez dodanie zera jako pierwszej współrzędnej, , otrzymujemy przekształcenie (macierz) Householdera w postaci

Pomnożenie macierzy z lewej strony przez spowoduje teraz wyzerowanie drugiej kolumny macierzy pod elementem , przy czym pierwszy wiersz i pierwsza kolumna pozostaną niezmienione. Postępując tak dalej razy (albo razy gdy ) otrzymujemy

gdzie jest uogólnioną macierzą trójkątną górną, tzn. dla . Stąd, podstawiając , dostajemy rozkład macierzy na iloczyn ortogonalno-trójkątny

Rzeczywiście, macierz jest ortogonalna, bo

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

gdzie . Rozbijając wektor na , gdzie i , oraz macierz na

gdzie jest macierzą trójkątną górną, a jest macierzą zerową wymiaru , otrzymujemy

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

oraz .

Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde z kolejnych przekształceń Householdera wyznaczamy przez obliczenie oraz współrzędnych wektora . Wektor ten ma tylko współrzędnych niezerowych, a ponadto dla . Dzięki takiej reprezentacji , mnożenia możemy dla dowolnego realizować według wzoru

gdzie .

Uwzględnizjąc obecność zerowych elementów w , przejście od macierzy do kosztuje rzędu operacji arytmetycznych i obliczenie jednego pierwiastka kwadratowego. Cały rozkład kosztuje więc rzędu (dla dużych i )

operacji arytmetycznych i pierwiastków kwadratowych. Zauważmy, że w przypadku , a więc dla kwadratowego układu równań, koszt ten wynosi i jest dwa razy większy od kosztu eliminacji Gaussa.

Implementacja

Cała informacja o przekształceniu Householdera znajduje się w wektorze oraz czynniku skalującym --- 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 przez macierz Householdera .

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ż miejsc w pamięci (chociaż, przypomnijmy raz jeszcze, cała informacja o to tylko liczb). Ponadto, mnożenie przez macierz to aż działań arytmetycznych.

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

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 , a więc uzyskaliśmu -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 będzie zbiorem wartości własnych macierzy . Definiujemy

(Jeśli w mianowniku pojawiłoby się zero, kładziemy ).

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 będzie rozwiązaniem zadania najmniejszych kwadratów dla niezerowej prawej strony ,

i niech będzie rozwiązaniem zadania zaburzonego

przy czym zakładamy, że

gdzie jest dostatecznie małe.

Oznaczmy

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

Wtedy

Generalnie więc, jeśli reszta jest mała, wrażliwość na zaburzenia jest na poziomie . Ale jeśli reszta jest duża (tzn. prawa strona jest taka, że nie można dobrze spełnić równania w sensie średniokwadratowym), wtedy wrażliwość może być daleko większa.

Wniosek

W przypadku, gdy , 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 , rozwiązanie obliczone (kosztowniejszym) rozkładem QR będzie miało błąd na poziomie , a tymczasem rozwiązanie wyznaczone z układu równań normalnych będzie obarczone błędem na poziomie .

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 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 miało minimalną normę euklidesową spośród wszystkich spełniających warunek , 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.