Aktualna wersja na dzień 11:17, 12 wrz 2023

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 - A x | |_{2}

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

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 - A x | |_{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 $f : [a, b] \to R$ obserwujemy jej wartości $f_{i}$ (dokładne lub zaburzone) w punktach $t_{i}$ , $1 \leq i \leq 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 $A = (a_{i, j}) \in R^{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 $‖ b - A x ‖_{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. $A x \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, $A \in R^{m \times n}$ , taką, że

m \geq n = rank (A),

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

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

‖ b - A x^{*} ‖_{2} = \min_{x \in R^{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} A x = 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 \neq 0$ mamy $x^{T} (A^{T} A) x = (A x)^{T} (A x) = ‖ A x ‖_{2} > 0$ , przy czym ostatnia nierówność wynika z faktu, że kolumny macierzy $A$ są liniowo niezależne i dlatego $A x \neq 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_{ν}$ powstanie po drodze dodatkowych błędów, które mogą nawet zmienić rząd macierzy. Na przykład, dla macierzy

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

mamy

A^{T} A = (\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})

Jeśli $ϵ^{2} < ν$ to $f l_{ν} (1 + ϵ^{2}) = 1$ , co implikuje $rank (f l_{ν} (A^{T} A)) = 1$ , podczas, gdy $rank (f l_{ν} (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 $w \in R^{m}$ o normie $‖ w ‖_{2} = \sqrt{w^{T} w} = 1$ , odbicie (macierz) Householdera zdefiniowane jest jako

H = I - 2 w w^{T}

Zauważmy, że

H x = 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 $H x$ 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 - 2 w w^{T})^{2} = I - 4 w w^{T} + 4 w (w^{T} w) w^{T} = I

,

oraz

H^{T} = (I - 2 w 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,

‖ H x ‖_{2} = \sqrt{(H x)^{T} (H x)} = \sqrt{x^{T} (H^{T} H) x} = \sqrt{x^{T} x} = ‖ x ‖_{2}

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

H x = (I - 2 w w^{T}) x = α e

Załóżmy dla uproszczenia, że $‖ e ‖_{2} = 1$ . Aby wyznaczyć $H$ zauważmy, że

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

a ponieważ $α = \pm ‖ x ‖_{2}$ i $‖ w ‖_{2} = 1$ to

w = \frac{x \mp ‖ x ‖_{2} e}{‖ x \mp ‖ x ‖_{2} e ‖_{2}}

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

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

gdzie

u_{i} = {\begin{cases} x_{1} \mp ‖ x ‖_{2} & i = 1, \\ x_{i} & 2 \leq i \leq m, \end{cases}

oraz

\begin{aligned} γ & = \frac{1}{2} ‖ u ‖_{2}^{2} = \frac{1}{2} ((x_{1} \mp ‖ x ‖_{2})^{2} + \sum_{i = 2}^{m} x_{i}^{2}) \\ = \frac{1}{2} (\sum_{i = 1}^{m} x_{i}^{2} + ‖ x ‖_{2}^{2} \mp 2 x_{1} ‖ x ‖_{2}) = ‖ x ‖_{2}^{2} \mp x_{1} ‖ x ‖_{2} . \end{aligned}

Otrzymaliśmy dwa odbicia Householdera przekształcające dany wektor $x$ na kierunek pierwszego wersora, w zależności od wybranego znaku przy $‖ x ‖_{2}$ . Ustalimy ten znak na plus gdy $x_{1} \geq 0$ oraz na minus gdy $x_{1} < 0$ , co pozwoli na obliczenie $u_{1}$ i $γ$ z małym błędem względem w $f l_{ν}$ . Wtedy bowiem mamy

u_{1} = {\begin{cases} x_{1} + ‖ x ‖_{2} & x_{1} \geq 0, \\ x_{1} - ‖ x ‖_{2} & x_{1} < 0, \end{cases}

oraz $γ = ‖ x ‖_{2}^{2} + | x_{1} | ‖ x ‖_{2}$ , czyli zawsze dodajemy liczby tych samych znaków. Ponadto pierwsza współrzędna wektora $H x$ jest równa $- ‖ x ‖_{2}$ , gdy $x_{1} \geq 0$ , a $+ ‖ x ‖_{2}$ jeśli $x_{1} < 0$ .

Rozkład QR

Odbić Householdera można użyć do rozkładu macierzy $A \in R^{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} / γ_{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 ${\bar{H}}_{2} = I_{m - 1} - v_{2} v_{2}^{T} / γ_{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 $v_{2} \in R^{m - 1}$ do wektora $u_{2} \in R^{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} / γ_{2}$ w $R^{m}$ postaci

H_{2} = (\begin{array}{cccc} 1 & 0^{T} \\ 0 & {\bar{H}}_{2} \end{array})

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} \dots H_{2} H_{1} A = R

,

gdzie $R \in R^{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} \dots H_{n}$ , dostajemy rozkład macierzy na iloczyn ortogonalno-trójkątny

A = Q \cdot R .

Rzeczywiście, macierz $Q \in R^{m \times m}$ jest ortogonalna, bo

\begin{aligned} Q^{- 1} & = (H_{1} H_{2} \dots H_{n})^{- 1} = H_{n}^{- 1} \dots H_{2}^{- 1} H_{1}^{- 1} \\ = H_{n}^{T} \dots H_{2}^{T} H_{1}^{T} = (H_{1} H_{2} \dots H_{n})^{T} = Q^{T} . \end{aligned}

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

\begin{aligned} ‖ r ‖_{2} & = ‖ b - A x ‖_{2} = ‖ b - Q R x ‖_{2} \\ = ‖ Q (Q^{T} b - R x) ‖_{2} = ‖ c - R x ‖_{2}, \end{aligned}

gdzie $c = Q^{T} b = H_{n} \dots H_{2} H_{1} b$ . Rozbijając wektor $c$ na $c = (c_{I}, c_{I I})^{T}$ , gdzie $c_{I} \in R^{n}$ i $c_{I I} \in R^{m - n}$ , oraz macierz $R$ na

R = (\begin{array}{c} R_{I} \\ 0 \end{array})

,

gdzie $R_{I} \in R^{n \times n}$ jest macierzą trójkątną górną, a $0$ jest macierzą zerową wymiaru $(m - n) \times n$ , otrzymujemy

‖ r ‖_{2}^{2} = ‖ c_{I} - R_{I} x ‖_{2}^{2} + ‖ c_{I I} ‖_{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 $‖ r^{*} ‖_{2} = ‖ b - A x^{*} ‖_{2} = ‖ c_{I I} ‖_{2}$ .

Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde z kolejnych przekształceń Householdera $H_{k}$ wyznaczamy przez obliczenie $γ_{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 \leq i \leq 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 / γ_{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 = Q R$ 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} + 2 n^{2} (m - n) = 2 n^{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 $γ$ --- 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}{γ} 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>\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:

H x = (I - \frac{1}{γ} u u^{T}) x = x - \frac{1}{γ} u u^{T} x = x - \frac{1}{γ} u (u^{T} x)

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

<math>\omega</math> = u'*x;
y = x - <math>\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 $Σ (A)$ będzie zbiorem wartości własnych macierzy $A^{T} A$ . Definiujemy

{cond}_{2} (A) = \sqrt{\frac{\max {λ : λ \in Σ (A)}}{\min {λ : λ \in Σ (A)}}}

(Jeśli w mianowniku pojawiłoby się zero, kładziemy ${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 $x$ będzie rozwiązaniem zadania najmniejszych kwadratów dla niezerowej prawej strony $b$ ,

| | b - A x | |_{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}} \leq ϵ

,

gdzie $ϵ$ jest dostatecznie małe.

Oznaczmy

\sin (θ) = \frac{| | b - A x | |_{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}} ≲ (\frac{2 {cond}_{2} (A)}{\cos (θ)} + \tan (θ) {cond}_{2}^{2} (A)) \cdot ϵ

Generalnie więc, jeśli reszta $| | b - A x | |_{2}$ jest mała, wrażliwość na zaburzenia jest na poziomie $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 A x$ w sensie średniokwadratowym), wtedy wrażliwość może być daleko większa.

Wniosek

W przypadku, gdy $m ≫ 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 (θ) ≪ 1$ , rozwiązanie obliczone (kosztowniejszym) rozkładem QR będzie miało błąd na poziomie ${cond}_{2} (A)$ , a tymczasem rozwiązanie wyznaczone z układu równań normalnych będzie obarczone błędem na poziomie ${cond}_{2}^{2} (A) > {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 - A x | |_{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.

MN12: Różnice pomiędzy wersjami

Aktualna wersja na dzień 11:17, 12 wrz 2023

Spis treści

Nadokreślone układy równań liniowych

Dopasowanie krzywej minimalizującej błąd średniokwadratowy

Układ równań normalnych

Odbicia Householdera

Rozkład QR

Implementacja

Uwarunkowanie

Biblioteki

Literatura

Menu nawigacyjne

Działania na stronie

Opcje strony

Narzędzia osobiste

Nawigacja

Szukaj

Narzędzia

@@ Linia 1: / Linia 1: @@
- <span id="sec:eigenvalue" \>
-Niech będzie dana rzeczywista kwadratowa macierz <math>A</math> wymiaru <math>N</math>. Wektorem własnym <math>x\in C^N</math> oraz odpowiadającą mu wartością własną <math>\lambda \in C</math> nazwiemy taką parę, dla której <span id=""/> <math>
+<!--
+Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
+Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki
+wprowadzi� przykry@mimuw.edu.pl
+-->
+=Nadokreślone układy równań liniowych=
+{{powrot |Metody numeryczne | do strony głównej
+przedmiotu <strong>Metody numeryczne</strong>}}
+Zajmiemy się zadaniem wygładzania liniowego,
+nazywanym też <strong>liniowym zadaniem najmniejszych kwadratów</strong>.
+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 <math>x</math>, które minimalizuje resztę,
-A x = \lambda x, </math> przy czym <math>x\neq 0</math>.
+<center><math>||b-Ax||_2</math></center>
-Zadanie wyznaczania  dupa '''hopla'''  wartości własnych  '''hopla'''  i wektorów własnych macierzy ma bardzo szerokie zastosowania w tak odległych do siebie dziedzinach jak np. analiza odporności konstrukcji mechanicznych (wieżowce, mosty, wagony kolejowe) na wibracje, czy też rankingowanie stron internetowych w wyszukiwarce Google.
+Jest to praktycznie bardzo często pojawiające się zadanie, a autorem pierwszego
+rozwiązania był nie kto inny jak sam wielki Gauss.
-{{przyklad|Odporność budynku na trzęsienie ziemi||
+[[grafika:Gauss.jpg|thumb|right||Carl Friedrich Gauss<br>  [[Biografia Gauss|Zobacz biografię]]]]
+Okazuje się bowiem, że jeśli np. potraktować <math>b</math> jako dane eksperymentalne (obarczone pewnym losowym błędem pomiaru o rozkładzie normalnym), a <math>x</math> --- parametrami zależności liniowej dla punktów pomiaru zadanych w macierzy <math>A</math>, to <math>x</math> minimalizujący <math>||b-Ax||_2</math> (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.
-Rozważmy prosty układ mechaniczny opisujący, naturalnie w pewnym jedynie przybliżeniu, zachowanie się układu <math>N</math> ciężkich płyt połączonych ze sobą relatywnie elatycznymi dźwigarami --- co może np. modelować konstrukcję wieżowca.
+==Dopasowanie krzywej minimalizującej błąd średniokwadratowy==
-\rysunek{}{Model wieżowca poddanego drganiom poprzecznym}
+<div style="margin-top:1em; padding-top,padding-bottom:1em;">
+<span  style="font-variant:small-caps;">Przykład</span>
+<div class="solution" style="margin-left,margin-right:3em;">
-Wiadomo, że jeśli częstotliwości drgań własnych tego wieżowca będą bliskie częstotliwości siły wymuszającej (o niewielkiej amplitudzie), to konstrukcja wpadnie w rezonans i w końcu rozpadnie się wskutek zbyt wielkich przemieszczeń. Wychylenia naszych płyt z położenia równowagi są opisywane układem pewnych równań różniczkowych. Teoria matematyczna takich równań różniczkowych pokazuje, że częstotliwości drgań własnych to nic innego jak <em>wartości własne</em> pewnej \CHECK{niesymetrycznej} macierzy wymiaru <math>2N</math>, która powstaje ze współczynników  równania różniczkowego opisującego dynamikę tego układu.  }}
+Przypuśćmy, że dla pewnej funkcji
+<math>f:[a,b]\to R</math> obserwujemy jej wartości <math>f_i</math> (dokładne lub
+zaburzone) w punktach <math>t_i</math>, <math>1\le i\le m</math>. Funkcję tę
+chcielibyśmy przybliżyć inną funkcją <math>w</math> należącą do
+pewnej <math>n</math> wymiarowej przestrzeni liniowej <math>W</math>, np. przestrzeni
+wielomianów stopnia mniejszego niż <math>n</math>. Jakość przybliżenia
+mierzymy, sprawdzając, ''jak dokładnie spełniona jest przybliżona równość <math>f_i \approx w(t_i)</math>'', dokładniej, badając tzw. <strong>błąd średniokwadratowy</strong>,
-{{przyklad|Macierz Google'a||
+<center><math>
+  \frac{1}{m}\sum_{i=1}^m (f_i-w(t_i))^2</math></center>
+Wybierając pewną bazę <math>(w_j)_{j=1}^n</math> w <math>W</math> i rozwijając <math>w</math>
+w tej bazie, <math>w(t)=\sum_{j=1}^n c_jw_j(t)</math>, sprowadzamy problem
+do minimalizacji
+<center><math>\sum_{i=1}^m\left(f_i-\sum_{j=1}^n c_jw_j(t_i)\right)^2
+</math></center>
-Podstawowy algorytm rankingowania stron WWW w  [http://www.wikipedia.org/pagerank[wyszukiwarce Google]]  sprowadza się do znalezienia rzeczywistego <em>wektora własnego</em> <math>\pi</math> pewnej silnie rozrzedzonej macierzy <math>A</math> (gigantycznego rozmiaru, równego liczbie indeksowanych stron, czyli w chwili pisania tego tekstu około XXXX stron), odpowiadającego wartości własnej równej 1:
+względem <math>c_j</math>, a więc do zadania wygładzania liniowego.
-\[ A \pi = \pi. \]
+Rzeczywiście, kładąc
+<math>A=(a_{i,j})\in R^{m\times n}</math> z <math>a_{i,j}=w_j(t_i)</math>,
+<math>b=(f_i)_{i=1}^m</math> i <math>x=(c_j)_{j=1}^n</math>, reszta jest równa <math>\| b-A x\|_2^2</math>, a minimalizacja reszty jest oczywiście równoważna minimalizacji błędu średniokwadratowego.
-Współrzędne wektora <math>\pi</math> interpretuje się jako wartość rankingową kolejnych stron WWW. Aby wszystko miało sens, współrzędne wektora muszą być  z przedziału [0,1].  Pewne twierdzenia matematyczne i subtelny dobór macierzy <math>A</math> gwarantują, że taki wektor <math>\pi</math> zawsze istnieje i jest jedyny! Co więcej, wartość 1 jest dominującą wartością własną <math>A</math>, a to z kolei ma ważne znaczenie dla tzw.  [[#sec:metoda-potegowa|metody potęgowej]]  numerycznego wyznaczania takiego wektora. }}
+[[Image:MNaproksymacjal2.png|thumb|550px|center|Wielomian <math>w</math> (czerwony) stopnia 3, aproksymujący 7 zadanych wartości (zaznaczone na zielono) danej funkcji <math>f</math> w sensie minimalizacji błędu średniokwadratowego]]
-{{przyklad|Wyznaczanie miejsc zerowych wielomianu||
+Powyższe zadanie aproksymacji średniokwadratowej w zadanych węzłach <math>(x_i,y_i)</math>, <math>i=1,\ldots,m</math>. wielomianem stopnia ''co najwyżej'' <math>N</math>, realizuje w Octave funkcja <code style="color: #006">polyfit(x,y,N)</code>. (Co dostaniemy, gdy <math>N=m-1</math>?)
+Można pokazać, że rozwiązanie minimalizujące błąd średniokwadratowy jest najbardziej prawdopodobnym zestawem parametrów naszego (liniowego) modelu, gdy zmierzone wartości <math>f_i</math> mogą być zaburzone losowym błędem pomiarowym.
+</div></div>
-Jak wiadomo, wartości własne to miejsca zerowe wielomianu charakterystycznego macierzy <math>P(\lambda) = \det(A - \lambda I)</math>. Zachodzi także fakt odwrotny, to znaczy miejsca zerowe wielomianu są wartościami pewnej macierzy, np. miejsca zerowe wielomianu  \[ p(\lambda) = p_1 \lambda^N + \ldots + p_N \lambda + p_{N+1} \] są wartościami własnymi m.in. macierzy stowarzyszonej, \[ A = \begin{pmatrix}  -p_2/p_1 	& -p_3/p_1 	& \cdots 	& -p_{N+1}/p_1\\ 1 		& 		& 		& 	\\ & 1		& 		&	\\ & 		& \ddots	&	\\ & 		& 		& 1 \end{pmatrix} \] Funkcja Octave'a  '''compan(p)'''  wyznacza macierz stowarzyszoną dla zadanego wielomianu o współczynnikach w wektorze <math>p = [p_1,\ldots,p_N, p_{N+1}]^T</math>. Z tej macierzy korzysta następnie funkcja Octave'a  '''roots'''  właśnie w taki sposób wyznacza pierwiastki wielomianów: jako wartości własne macierzy stowarzyszonej. }}
+W kontekście nie-statystycznym, możemy myśleć o zadaniu wygładzania liniowego jako sposobie skrócenia listy parametrów <math>x</math> modelu przy zachowaniu przybliżonego spełnienia warunków modelu, tzn. <math>Ax\approx b</math>.
-{{przyklad|||
+Dodajmy, że spotyka się uogólnienie tego zadania w formie następującej: dla danych wartości <math>b\in R^m</math>, i danej funkcji <math>F:R^n\rightarrow R^m</math>, znaleźć <math>x\in R^n</math> minimalizujący resztę:
-Praktyczne zadanie z macierzą symetryczną }}
+<center><math>||b-F(x)||_2</math></center>
-W praktyce obliczeniowej spotyka się zazwyczaj kilka typów zagadnień:
+Właśnie tego typu <strong>nieliniowe zadanie najmniejszych kwadratów</strong> rozwiązują np. nasze przenośne [  odbiorniki GPS]... Na marginesie zauważmy, że gdy <math>F</math> 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 [[Metody optymalizacji|wykładzie z metod optymalizacji]].
+==Układ równań normalnych==
-*Wyznaczenie dominującej wartości własnej (to znaczy: największej co do modułu) i odpowiadającego jej wektora własnego (a może kilku wektorów?)
+Niech <math>A</math> będzie daną macierzą o <math>m</math> wierszach i <math>n</math> kolumnach,
-*Wyznaczenie najmniejszej co do modułu wartości własnej i wektorów jej odpowiadających (zauważmy, że to jest np. zadanie wyznaczenia {\em jądra macierzy osobliwej} --- wtedy wiemy a priori, że szukana najmniejsza co do modułu wartość własna to zero)
+<math>A\in R^{m\times n}</math>, taką, że
-*Wyznaczenie wartości własnej najbliższej zadanej liczbie (to jest właśnie odpowiedź na pytanie jak blisko częstości wymuszającej są częstości drgań własnych budynku)
-*Wyznaczenie wszystkich wartości własnych
-*Wyznaczenie wszystkich wartości i wektorów własnych (tzw. pełne zagadnienie własne)
+<center><math>m\,\ge\,n\,=\, \mbox{rank} (A),
+</math></center>
-Jak domyślamy się, dla macierzy rozrzedzonych dużego wymiaru pełne zagadnienie własne jest zbyt kosztowne, gdyż najczęściej macierz wektorów własnych --- nawet dla macierzy rzadkiej --- jest gęsta.
+albo równoważnie, taką że jej wektory kolumny są liniowo
+niezależne. Niech także dany będzie wektor <math>b\in R^m</math>.
+Jasne jest, że wtedy układ równań <math>A x= b</math> nie zawsze
+ma rozwiązanie - mówimy, że układ jest <strong>nadokreślony</strong>.
-Ponieważ w zastosowaniach bardzo często pojawiają się macierze rzeczywiste symetryczne (powyższe przykłady pokazują, że nie tylko!) szczegółową analizę metod numerycznych ograniczymy do tego przypadku, gdyż wtedy zachodzi
+<strong>Zadanie wygładzania liniowego</strong> polega na znalezieniu wektora
+<math>x^*\in R^n</math>, który minimalizuje <strong>wektor residualny</strong> (wektor reszty)
+<math>r= b-A x</math> w normie drugiej, tzn.
-{{twierdzenie|o symetrycznym zadaniu włanym|thm:symetric-eig|
+<center><math>\| b\,-\,A x^*\|_2\,=\,\min_{ x\in R^n}
+    \| b\,-\,A x\|_2</math></center>
+{{lemat|||
+Zadanie wygładzania liniowego ma jednoznaczne
+rozwiązanie <math>x^*</math>, które można scharakteryzować jako rozwiązanie układu równań
+<center><math>
+A^TA x\,=\,A^T\, b</math></center>
+}}
+Zauważmy, że jeśli macierz <math>A</math> jest kwadratowa, <math>m=n</math>, to
+rozwiązaniem jest <math>x^*=A^{-1} b</math> i residuum jest zerem.
+Zadanie wygładzania liniowego jest więc uogólnieniem
+rozwiązywania kwadratowych układów równań liniowych.
-Każda macierz rzeczywista symetryczna <math>A</math> wymiaru <math>N</math> ma rozkład \[ A = Q\Lambda Q^T, \] gdzie <math>Q\in R^{N\times N}</math> jest ortogonalna (tzn. <math>Q^TQ = I</math>), a jej kolumnami są wektory własne <math>A</math>, natomiast <math>\Lambda\in R^N</math> jest diagonalna z  wartościami własnymi <math>A</math> na diagonali: \[ \Lambda = \begin{pmatrix}\lambda_1 & & \\ & \ddots & \\ & & \lambda_N\end{pmatrix}. \]  }}
+Równanie powyższe nazywa się <strong>układem równań normalnych</strong>.
+Może ono nam sugerować sposób rozwiązania zadania wygładzania
+liniowego. Wystarczy bowiem pomnożyć macierz <math>A^T</math> przez <math>A</math> i
+rozwiązać układ normalny. Zauważmy ponadto, że macierz <math>A^TA</math>
+jest symetryczna i dodatnio określona, bo <math>(A^TA)^T=A^TA</math> i dla
+<math>x\ne 0</math> mamy
+<math>x^T(A^TA) x=(A x)^T(A x)=\|A x\|_2>0</math>, przy
+czym ostatnia nierówność wynika z faktu, że kolumny macierzy <math>A</math>
+są liniowo niezależne i dlatego <math>A x\ne 0</math>. Przy mnożeniu
+<math>A^T</math> przez <math>A</math> wystarczy więc obliczyć tylko elementy na głównej
+przekątnej i pod nią, a do rozwiązania równania z macierzą
+<math>A^TA</math> można zastosować [[MN05LAB|algorytm Cholesky'ego-Banachiewicza]]. Jak łatwo się przekonać, koszt takiego
+algorytmu wynosi <math>n^2(m+n/3)</math>, przy czym dominuje koszt mnożenia
+obliczenia macierzy <math>A^TA</math>.
+Ma on jednak pewne wady. Mnożenie macierzy powoduje w <math>fl_\nu</math>
+powstanie po drodze dodatkowych błędów, które mogą nawet
+zmienić rząd macierzy. Na przykład, dla macierzy
+<center><math>A\,=\,\left(\begin{array} {cccc}
+  &  1  &  1  &  1  \\
+  \epsilon  \\
+        &\epsilon \\
+        &     &\epsilon \\
+        &     &      &\epsilon \end{array} \right)
+</math></center>
-{{twierdzenie|Bauer-Fike|thm:Bauer-Fike|
+mamy
+<center><math>A^TA\,=\,\left(\begin{array} {cccc}
++\epsilon^2 & 1 & 1 & 1 \\
+& 1+\epsilon^2 & 1 & 1 \\
+& 1 & 1+\epsilon^2 & 1 \\
+& 1 & 1 & 1+\epsilon^2 \end{array} \right)</math></center>
+Jeśli <math>\epsilon^2<\nu</math> to <math>fl_\nu(1+\epsilon^2)=1</math>, co implikuje
+<math>\mbox{rank} (fl_\nu(A^TA))=1</math>, podczas, gdy <math>\mbox{rank} (fl_\nu(A))=4</math>.
+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==
-Niech <math>A\in R^{N\times N}</math> będzie diagonalizowalna, to znaczy dla pewnej macierzy <math>X</math> zachodzi  \[ X^{-1}  A X = \begin{pmatrix} \lambda_1 & & \\ & \ddots & \\ & & \lambda_N\end{pmatrix}, \]
+Dla danego wektora <math>w\in R^m</math> o normie
+<math>\| w\|_2=\sqrt{ w^T w}=1</math>,
+<strong>odbicie</strong> (macierz) <strong>Householdera</strong> zdefiniowane jest jako
-a więc (gdyż macierz po prawej stronie jest podobna do <math>A</math>) <math>\lambda_i\in C</math>, <math>i=1,\ldots,N</math> są wartościami własnymi <math>A</math>. Rozważmy macierz zaburzoną <math>\tilde{A}</math> i jakąś jej wartość własną <math>\tilde{\lambda}</math>. Wtedy istnieje wartość własna <math>\lambda_j</math> macierzy <math>A</math> taka, że
+<center><math>H\,=\,I\,-\,2 w w^T</math></center>
-\[ |\lambda_j - \tilde{\lambda}| \leq \mbox{cond}_2(X) ||A - \tilde{A}||_2. \]
+Zauważmy, że
-}}
+<center><math>H x\,=\, x\,-\,2( w^T x) w</math>,</center>
-Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia <math>X</math> jest ortogonalna, <math>X^{-1} = X^T</math>, to mamy <math>\mbox{cond}_2(X) = 1</math> i w konsekwencji zachodzi
+a ponieważ <math>( w^T x) w=( x, w)_2 w</math>
+jest rzutem prostopadłym <math>x</math> na kierunek wektora <math>w</math>
+(<math>(\cdot,\cdot)_2</math> oznacza iloczyn skalarny), to <math>H x</math> jest
+odbiciem lustrzanym wektora <math>x</math> względem hiperpłaszczyzny
+(wymiaru <math>m-1</math>) prostopadłej do <math>w</math>.
-{{wniosek|Wartości własne macierzy symetrycznej są doskonale uwarunkowane||
+Odbicia Householdera są przekształceniami nieosobliwymi
+spełniającymi
+<center><math>H^{-1}\,=\,H\,=\,H^T</math></center>
+Rzeczywiście, ponieważ <math>w</math> ma normę jednostkową, mamy
-Przy oznaczeniach jak  [[#thm:Bauer-Fike|twierdzeniu Bauera-Fike'a]] , jeśli dodatkowo założymy, że macierz <math>A</math> jest rzeczywista i symetryczna, to
+<center><math>H^2 \,=\, (I-2 w w^T)^2\,=\,
+   I-4 w w^T+4 w( w^T w) w^T \,=\, I</math>,</center>
-\[ \min_{j=1,\ldots,N}|\lambda_j - \tilde{\lambda}| \leq ||A - \tilde{A}||_2. \] }}
+oraz
-Z drugiej strony, dla macierzy niediagonalizowalnych, uwarunkowanie wartości własnych może być dowolnie duże, co ilustruje poniższy
+<center><math>H^T\,=\,(I-2 w w^T)^T\,=\,I-2( w^T)^T w^T\,=\,I.
+</math></center>
-{{przyklad|||
+W szczególności <math>H</math> jest więc przekształceniem <strong>ortogonalnym</strong>,
+<math>H^{-1}=H^T</math>, czyli nie zmienia długości wektora,
-\[ A_\epsilon = \begin{pmatrix} a & 1 \\ \epsilon & a \end{pmatrix} \] Weźmy dla uproszczenia <math>a=0</math>. Wartości własne <math>A_\epsilon</math> to zera wielomianu <math>p_\epsilon(\lambda) = \lambda^2 - \epsilon</math>, zatem <math>\lambda_\epsilon = \pm \sqrt{\epsilon}</math> i w konsekwencji
+<center><math>\|H x\|_2\,=\,\sqrt{(H x)^T(H x)}\,=\,
+    \sqrt{ x^T(H^TH) x}\,=\,\sqrt{ x^T x}\,=\,
+    \| x\|_2</math></center>
-\[ |\lambda_\epsilon - \lambda_0| / ||A_\epsilon - A_0|| = \sqrt{\epsilon}/\epsilon \rightarrow \infty, \]
+Odbicia Householdera zastosujemy do przeprowadzenia danego wektora
+<math>x\ne 0</math> na kierunek innego niezerowego wektora, powiedzmy
+<math>e</math>, tzn.
-gdy <math>\epsilon \rightarrow 0^+</math>, a więc uwarunkowanie takiego zadania jest nieskończone: dowolnie mała zmiana macierzy powoduje zaburzenie wartości własnych niewspółmiernie wielkie wobec zaburzenia danych. Dodatkowo, wartości własne i wektory własne macierzy <math>A</math> dla ujemnego parametru <math>\epsilon</math> są zespolone!
+<center><math>H x\,=\,(I-2 w w^T) x\,=\,\alpha\, e</math></center>
-\rysunek{eigencond.png}{Zachowanie się wartości własnych macierzy <math>A</math> (z parametrem <math>a=1</math>) w otoczeniu <math>\delta = 0</math>}
+<!--
+[[Image:MNhouseholderidea.png|thumb|300px|center|Odbicie Householdera]]
+-->
+Załóżmy dla uproszczenia, że <math>\| e\|_2=1</math>.
+Aby wyznaczyć <math>H</math> zauważmy, że
-}}
+<center><math>w\,=\,\frac{ x-\alpha e}{2( w^T x)},
+</math></center>
-Bardziej spektakularny przykład pochodzi od Wilkinsona:
+a ponieważ <math>\alpha=\pm\| x\|_2</math> i <math>\| w\|_2=1</math> to
-{{przyklad|Perfidny wielomian Wilkinsona||
+<center><math>w\,=\,\frac{ x\mp\| x\|_2 e}
+                {\| x\mp\| x\|_2 e\|_2}</math></center>
+W szczególności, jeśli <math>e= e_1</math> jest pierwszym
+wersorem, powyższe wzory dają
+<center><math>H\,=\,I\,-\,\frac{ u u^T}{\gamma},
+</math></center>
-Niech \[ p(\lambda) = (\lambda -1)(\lambda - 2) \cdots (\lambda - 20). \]
+gdzie
-Zmiana współczynnika przy <math>\lambda^{19}</math> o <math>10^{-7}</math> skutkuje presunięciem niektórych miejsc zerowych nawet o kilka jednostek na płaszczyźnie zespolonej! Poniżej pokazujemy to na numerycznym przykładzie, gdzie prócz w/w zaburzenia mamy dodatkowo z zaburzeniami powstałymi wskutek wyznaczenia współczynników wielomianu w arytmetyce zmiennoprzecinkowej.
+<center><math>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.</math></center>
-\rysunek{wilkinson.png}{Zera oryginalnego i lekko zaburzonego perfidnego wielomianu Wilkinsona.}
+oraz
-Jak widzimy, zera bardzo mało zaburzonego wielomianu mogą stać się wyraźnie nie-rzeczywiste!
+<center><math>\begin{align} \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\,
+x_1\|x\|_2\Big) \,=\,\|x\|_2^2\,\mp\,x_1 \|x\|_2.
+\end{align}</math></center>
-}}
+Otrzymaliśmy dwa odbicia Householdera przekształcające dany wektor
+<math>x</math> na kierunek pierwszego wersora, w zależności od wybranego
+znaku przy <math>\| x\|_2</math>. Ustalimy ten znak na plus gdy <math>x_1\ge 0</math>
+oraz na minus gdy <math>x_1<0</math>, co pozwoli na obliczenie <math>u_1</math> i <math>\gamma</math>
+z małym błędem względem w <math>fl_\nu</math>. Wtedy bowiem mamy
-Jeśli chodzi o wektory własne, ich wrażliwość na zaburzenia macierzy jest bardziej skomplikowana i zależy m.in. od  uwarunkowania wartości własnych (czego łatwo się domyślić) oraz od tego, jak blisko siebie leżą wartości własne.
+<center><math>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.</math></center>
+oraz <math>\gamma=\| x\|_2^2+|x_1|\,\| x\|_2</math>, czyli zawsze
+dodajemy liczby tych samych znaków. Ponadto pierwsza współrzędna
+wektora <math>H x</math> jest równa <math>-\| x\|_2</math>, gdy <math>x_1\ge 0</math>, a
+<math>+\| x\|_2</math> jeśli <math>x_1<0</math>.
+==Rozkład QR==
-Jak okaże się za chwilę, czasem warto mieć ogólne rozeznanie o tym, gdzie {\em z grubsza} leżą wartości własne danej macierzy <math>A</math>. W tym celu mogą być nam pomocne dwa fakty:
+Odbić Householdera można użyć do rozkładu macierzy
+<math>A\in R^{m\times n}</math> na iloczyn ortogonalno-trójkątny.
-{{fakt|||
+Niech <math>A=( a_1, a_2,\ldots, a_n)</math>, gdzie <math>a_j</math> są
+wektorami-kolumnami macierzy <math>A</math>. Wybierzmy pierwsze odbicie
+Householdera <math>H_1=I_m- u_1 u_1^T/\gamma_1</math> tak, aby
+przekształcało pierwszy wektor-kolumnę macierzy <math>A</math> na kierunek
+<math>e_1</math>. Efektem pomnożenia macierzy <math>A</math> z lewej strony przez
+<math>H_1</math> będzie wtedy macierz
-Dowolna wartość własna <math>\lambda\in C</math> macierzy <math>A</math> spełnia \[ |\lambda| \leq ||A||, \] gdzie <math>||A||</math> jest dowolną normą macierzową indukowaną przez normę wektorową. }}
+<center><math>A^{(1)}\,=\,( a^{(1)}_1,\ldots, a^{(1)}_n)
+   \,=\,(H_1 a_1,\ldots, H_1 a_n)</math>,</center>
-Rzeczywiście, skoro istnieje wektor <math>x\neq 0</math> taki, że <math>Ax = \lambda x</math>, to stąd <math>||Ax||/||x|| = |\lambda|</math>, więc fakt powyższy wynika już z definicji normy macierzy: <span id=""/> <math>
+w której pierwsza kolumna <math>a^{(1)}_1</math> ma niezerową tylko
+pierwszą współrzędną. W następnym kroku wybieramy drugie
+przekształcenie Householdera
+<math>\bar H_2=I_{m-1}- v_2 v_2^T/\gamma_2</math> wymiaru <math>m-1</math> tak,
+aby przeprowadzało wektor <math>(a^{(1)}_{i,2})_{i=2}^m</math> na kierunek
+pierwszego wersora w <math>R^{m-1}</math>. Rozszerzając <math>v_2\in R^{m-1}</math>
+do wektora <math>u_2\in R^m</math> przez dodanie zera jako pierwszej
+współrzędnej, <math>u_2=(0, v_2)^T</math>, otrzymujemy
+przekształcenie (macierz) Householdera
+<math>H_2=I_m- u_2 u_2^T/\gamma_2</math> w <math>R^m</math> postaci
+<center><math>H_2\,=\,\left(\begin{array} {cccc}
+&  0^T \\
+& \bar H_2  \end{array} \right)</math></center>
+Pomnożenie macierzy <math>A^{(1)}</math> z lewej strony przez <math>H_2</math> spowoduje
+teraz wyzerowanie drugiej kolumny macierzy pod elementem
+<math>a^{(1)}_{2,2}</math>, przy czym pierwszy wiersz i pierwsza kolumna
+pozostaną niezmienione. Postępując tak dalej <math>n</math> razy
+(albo <math>n-1</math> razy gdy <math>m=n</math>) otrzymujemy
-||A|| = \max_{y\neq 0}\frac{||Ay||}{||y||} \geq ||Ax||/||x||. </math>
+<center><math>H_nH_{n-1}\cdots H_2H_1A\,=\,R</math>,</center>
-Drugie twierdzenie jest równie proste w dowodzie, ale daje trochę więcej informacji o lokalizacji widma. {{twierdzenie|Gerszgorina||
+gdzie <math>R\in R^{m\times n}</math> jest uogólnioną macierzą trójkątną
+górną, tzn. <math>r_{i,j}=0</math> dla <math>i>j</math>. Stąd, podstawiając
+<math>Q=H_1H_2\cdots H_n</math>, dostajemy rozkład macierzy na iloczyn
+ortogonalno-trójkątny
+<center><math>
+  A\,=\,Q\cdot R.
+</math></center>
+Rzeczywiście, macierz <math>Q\in R^{m\times m}</math> jest ortogonalna, bo
-Wartości własne macierzy <math>A</math> leżą w sumie (teoriomnogościowej) dysków <math>K_i</math> na płaszczyźnie zespolonej, \[ K_i = \{z \in C: |z - a_{ii}| \leq \sum_{j\neq i} |a_{ij}| \}, \qquad i = 1,\ldots N. \] }}
+<center><math>\begin{align} 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.
+\end{align}</math></center>
-{{przyklad|Koła Gerszgorina||
+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
+<center><math>\begin{align} \| r\|_2 &= \| b-A x\|_2\;=\;\| b-QR x\|_2 \\
+    &= \|Q(Q^T b-R x)\|_2 \;=\;\| c-R x\|_2,
+\end{align}</math></center>
+gdzie <math>c=Q^T b=H_n\cdots H_2H_1 b</math>.
+Rozbijając wektor <math>c</math> na <math>c=( c_I, c_{II})^T</math>,
+gdzie <math>c_I\in R^n</math> i <math>c_{II}\in R^{m-n}</math>, oraz macierz
+<math>R</math> na
-Niech   \[ A = \begin{pmatrix} 1.08930   & 1.38209  & -1.00037  &  0.69355  &  2.32178 \\ 0.14211  &  1.74696 &   1.68440 &   0.30664 &   1.26718 \\ -0.74620  &  2.02686 &  -0.68293 &   0.19684 &   0.35854 \\ 0.83517  &  0.74987 &   1.71331 &   1.09765 &  -0.44321 \\ 1.02132  & -2.62155 &   0.79247 &   1.11408 &   0.48076 \\ \end{pmatrix} \]
+<center><math>R\,=\,\left(\begin{array} {c} R_I \\ 0\end{array} \right)</math>,</center>
-\rysunek{gershgorindisks.png}{Lokalizacja wartości własnych macierzy <math>A</math> kołami Gerszgorina oraz zgrubna lokalizacja wewnątrz okręgu o promieniu równym <math>||A||_1</math>. Dokładne wartości własne zaznaczone trójkącikami.}
+gdzie <math>R_I\in R^{n\times n}</math> jest macierzą trójkątną górną, a
+<math>0</math> jest macierzą zerową wymiaru <math>(m-n)\times n</math>, otrzymujemy
-}}
+<center><math>\| r\|_2^2\;=\;\| c_I-R_I x\|_2^2\,+\,
+     \| c_{II}\|_2^2.
+</math></center>
-{{przyklad|Widmo macierzy jednowymiarowego Laplasjanu||
+Rozwiązanie <math>x^*</math> zadania wygładzania jest więc
+rozwiązaniem układu liniowego trójkątnego,
+<center><math>x^*\,=\,R_I^{-1} c_I</math>,</center>
+oraz <math>\| r^*\|_2=\| b-A x^*\|_2=\| c_{II}\|_2</math>.
-Norma daje:
+Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde
+z kolejnych przekształceń Householdera <math>H_k</math> wyznaczamy przez
+obliczenie <math>\gamma_k</math> oraz współrzędnych wektora <math>u_k</math>.
+Wektor ten ma tylko <math>m-k+1</math> współrzędnych niezerowych, a ponadto
+<math>u_{k,i}=a^{(k-1)}_{i,k}</math> dla <math>k+1\le i\le m</math>. Dzięki takiej
+reprezentacji <math>H_k</math>, mnożenia <math>H_k x</math> możemy dla dowolnego
+<math>x</math> realizować według wzoru
-Tw. Gerszgorina daje:
+<center><math>(H_k x)_i\,=\,x_i\,-\,s\,u_{k,i},
+</math></center>
-W rzeczywistości,
+gdzie <math>s= u_k^T x/\gamma_k</math>.
+Uwzględnizjąc obecność zerowych elementów w <math>u_k</math>,
+przejście od macierzy <math>A^{(k-1)}</math> do <math>A^{(k)}</math> kosztuje rzędu
+<math>4(m-k+1)(n-k)</math> operacji arytmetycznych i obliczenie jednego
+pierwiastka kwadratowego. Cały rozkład <math>A=QR</math> kosztuje więc
+rzędu (dla dużych <math>m</math> i <math>n</math>)
+<center><math>\sum_{k=1}^n 4(m-k+1)(n-k)\,\approx\,\frac 43n^3+2n^2(m-n)
+   \,=\,2n^2(m-n/3)
+</math></center>
-}}
+operacji arytmetycznych i <math>n</math> pierwiastków kwadratowych. Zauważmy,
+że w przypadku <math>m=n</math>, a więc dla kwadratowego układu równań,
+koszt ten wynosi <math>(4/3)n^3</math> i jest dwa razy większy od kosztu
+eliminacji Gaussa.
+===Implementacja===
+Cała informacja o przekształceniu Householdera znajduje się w wektorze <math>u</math> oraz czynniku skalującym <math>\gamma</math> --- 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 <math>x</math> przez macierz Householdera <math>H = I - \frac{1}{\gamma}uu^T</math>.
+Nie popełnijmy jednak częstego błędu, prostodusznie implementując to mnożenie (przykładowo, w Octave) jako
+ <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>H = eye(length(u)) - (u*u') / <math>\gamma</math>;
+y = H*x;
+</pre></div>
+Gdybyśmy użyli takiej implementacji, potrzebowalibyśmy aż <math>O(N^2)</math> miejsc w pamięci (chociaż, przypomnijmy raz jeszcze, cała informacja o <math>H</math> to tylko <math>O(N)</math> liczb). Ponadto, mnożenie przez macierz to aż <math>O(N^2)</math> działań arytmetycznych.
- <span id="sec:metoda-potegowa" \>
+Aby znacznie lepiej skorzystać z bardzo specyficznej postaci macierzy <math>H</math>, która jest po prostu zaburzeniem macierzy identyczności macierzą rzędu co najwyżej 1, wystarczy w odpowiednim miejscu wstawić nawiasy:
-Przypuśćmy, że wartości własne macierzy <math>A\in R^{N\times N}</math> spełniają <span id=""/> <math>
+<center><math>
+Hx = \left(I - \frac{1}{\gamma}uu^T\right) \, x = x - \frac{1}{\gamma}uu^Tx =
+x - \frac{1}{\gamma}u(u^Tx)</math></center>
+Stąd <strong>prawidłowa</strong> implementacja mnożenia przez macierz Householdera:
+ <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre><math>\omega</math> = u'*x;
+y = x - <math>\frac{\omega}{\gamma}</math>*u;
+</pre></div>
+Tym razem wcale nie potrzeba dodatkowej pamięci, a koszt algorytmu jest liniowy(!) względem <math>N</math>, a więc uzyskaliśmu <math>N</math>-krotne przyspieszenie w porównaniu z poprzednim!
-|\lambda_1| > |\lambda_2| \geq \ldots \geq |\lambda_N|, </math> (to znaczy, istnieje dokładnie jedna <em>dominująca</em> wartość własna macierzy <math>A</math>.
+Jest to całkiem typowe w numeryce:
-Załóżmy także, że istnieje baza złożona z wektorów własnych <math>q_1,\ldots,q_N</math> tej macierzy (tak jest np. dla macierzy symetrycznej na mocy \link{thm:symetric-eig}{twierdzenia o własnościach symetrycznego zadania własnego}).
+<blockquote  style="background-color: #fefeee; padding:1em;  margin-left,margin-right:2em;  margin-top,margin-bottom: 1em;">
+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).
+</blockquote>
-Kierunek własny <math>q_k</math> jakiejś macierzy <math>A</math> ma taką własność, że poddany działaniu przekształcenia <math>A</math> wydłuża się <math>\lambda_k</math> razy, wobec tego, dowolny wektor <math>x\in R^N</math> poddany działaniu <math>A</math> najbardziej wydłuży się w kierunku <math>q_1</math>. Iterując tę procedurę, powinniśmy dostawać w wyniku wektory, w których coraz bardziej dominuje kierunek <math>q_1</math>. Formalnie, niech
+==Uwarunkowanie==
-<span id=""/> <math>
+Ł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|Uwarunkowanie macierzy prostokątnej w normie euklidesowej|
-x = \alpha_1q_1 + \ldots + \alpha_Nq_N, </math>
+Niech <math>\Sigma(A)</math> będzie zbiorem wartości własnych macierzy <math>A^TA</math>. Definiujemy
-wtedy
+<center><math>\mbox{cond} _2(A) = \sqrt{\frac{\max\{\lambda: \lambda \in \Sigma(A)\}}{\min\{\lambda: \lambda \in \Sigma(A)\}}}</math></center>
-<span id=""/> <math>
+(Jeśli w mianowniku pojawiłoby się zero, kładziemy <math>\mbox{cond} _2(A) = +\infty</math>).
+}}
+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|O uwarunkowaniu zadania wygładzania liniowego|
-Ax = A \left( \sum_i \alpha_iq_i \right) = \sum_i \alpha_i A q_i  =  \sum_i \alpha_i \lambda_i q_i </math>
+Niech <math>x</math> będzie rozwiązaniem zadania najmniejszych kwadratów dla niezerowej prawej strony <math>b</math>,
+<center><math>
+||b-Ax||_2\rightarrow \min{} !
+</math></center>
-i w konsekwencji  <span id=""/> <math>
+i niech <math>\widetilde{x}</math> będzie rozwiązaniem zadania zaburzonego <center><math>
+||\widetilde{b}-\widetilde{A}\widetilde{x}||_2\rightarrow \min{} !</math>,</center>
+przy czym zakładamy, że
+<center><math>
+\frac{||\widetilde{b}-b||_2}{||b||_2}, \quad \frac{||\widetilde{A}-A||_2}{||A||_2} \leq \epsilon</math>,</center>
+gdzie <math>\epsilon</math> jest dostatecznie małe.
+Oznaczmy
-A^kx = \sum_i \alpha_i \lambda_i^k q_i = \lambda_1^k\left(\alpha_1q_1 + \alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^kq_2 + \ldots  + \alpha_N\left(\frac{\lambda_N}{\lambda_1}\right)^kq_N \right). </math>
+<center><math>\sin(\theta) = \frac{||b-Ax||_2}{||b||_2} < 1
+</math></center>
-Ponieważ z założenia, że istnieje dokładnie jedna dominująca wartość własna, <math>\left|\frac{\lambda_N}{\lambda_1}\right| < 1</math>, to wyrażenie w nawiasie dąży do <math>\alpha_1q_1</math> i w konsekwencji wektory <math>x_k = A^kx</math> dążą, gdy <math>k\rightarrow\infty</math>, do kierunku wektora własnego <math>q_1</math>, to znaczy wektora odpowiadającego dominującej wartości własnej <math>A</math>  (o ile tylko <math>\alpha_1 \neq 0</math>).
+--- będzie to miara, jak bardzo jesteśmy w stanie zminimalizować resztę w oryginalnym zadaniu.
-Szybkość zbieżności metody potęgowej jest liniowa, o współczynniku zależnym od stosunku <math>\lambda_2/\lambda_1|</math>. W patologicznym przypadku, gdy <math>|\lambda_1| \approx |\lambda_2|</math>, może więc okazać się, że metoda praktycznie nie jest zbieżna.
+Wtedy
+<center><math>\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</math></center>
+}}
-W praktyce nie wyznaczamy wzorem <math>x_k = (A^k)\cdot x</math>, lecz raczej korzystamy z metody iteracyjnej
+Generalnie więc, jeśli reszta <math>||b-Ax||_2</math> jest mała, wrażliwość na zaburzenia jest na poziomie <math>\mbox{cond} (A)</math>. Ale jeśli reszta jest duża (tzn. prawa strona jest taka, że nie można dobrze spełnić równania <math>b\approx Ax</math> w sensie średniokwadratowym), wtedy wrażliwość może być daleko większa.
+{{wniosek|||
+W przypadku, gdy <math>m \gg n</math>, 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 <math>\sin(\theta) \ll 1</math>, rozwiązanie obliczone (kosztowniejszym) rozkładem QR będzie miało błąd na poziomie <math>\mbox{cond} _2(A)</math>, a tymczasem rozwiązanie wyznaczone z układu równań normalnych będzie obarczone błędem na poziomie <math>\mbox{cond} _2^2(A) >  \mbox{cond} _2(A)</math>.
+}}
- $x_0$ = dowolny wektor startowy; k = 0;
+==Biblioteki==
- '''while''' ( !stop )
- {
-   $y_k$ = $Ax_{k-1}$;
-   $x_k$ = $y_k/||y_k||_\infty$;
-   k++;
- }
+W Octave, zadanie najmniejszych kwadratów rozwiązujemy praktycznie tak samo, jak równanie liniowe:
-Warunek normowania ma m.in. na celu zapobieżenie powstawania nadmiaru i niedomiaru (gdy <math>|\lambda_1| < 1</math>, to <math>||A^kx|| \rightarrow 0</math>, a gdy  <math>|\lambda_1| > 1</math>, to <math>||A^kx|| \rightarrow \infty</math>). Przy okazji, <math>||y_k||_\infty \rightarrow |\lambda_1|</math>, a więc mamy także sposób na wyznaczenie przybliżenia dominującej wartości własnej.
+ <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>x = A \ b;
+</pre></div>
+Dla zadania najmniejszych kwadratów mamy dwie podstawowe funkcje LAPACKa: <code style="color: #903">DGELS</code>, która rozwiązuje dokładnie zadanie takie, jak postawiliśmy w wykładzie, to znaczy w przypadku, gdy macierz <math>A</math> jest pełnego rzędu --- wykorzystując rozkład QR, który omówiliśmy.
-Zazwyczaj jako warunek stopu wybiera się kryterium małej poprawki, <math>||x_k - x_{k-1}|| \leq \epsilon</math>, lub warunek małego residuum, <math>||Ax_k - \lambda_{1,k} x_k||\leq \epsilon</math>, gdzie <math>\lambda_{1,k}</math> jest przybliżeniem <math>\lambda_1</math> dostępnym na <math>k</math>-tej iteracji.
+Natomiast dla przypadku, <strong>gdy macierz nie jest pełnego rzędu</strong>, działa funkcja <code style="color: #903">DGELSS</code>. 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 <math>x</math> miało <strong>minimalną normę</strong> euklidesową spośród wszystkich spełniających warunek <math>||b-Ax||_2 \rightarrow \min !</math>, 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==
-\rysunek{}{Zasada działania metody potęgowej}
+W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 5.3</b> w
+* D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.
-Metoda potęgowa doskonale sprawdza się, gdy macierz <math>A</math> jest macierzą rozrzedzoną --- np. w przypadku macierzy Google'a.
+Bardzo dużo na temat rozwiązywania liniowego zadania najmniejszych kwadratów można dowiedzieć się z książki
+* <span style="font-variant:small-caps">A. Kiełbasiński, H. Schwetlick</span>, <cite>Numeryczna algebra liniowa</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa, 1992.
-Zauważmy, że dla dowolnej macierzy kwadratowej <math>A</math> o wartościach własnych <math>\lambda_k</math> i odpowiadających im wektorach własnych <math>q_k</math>, mamy:
-*Macierz <math>A-\sigma I</math> ma wartości własne <math>\lambda_k - \sigma</math> oraz wektory własne <math>q_k</math>,
-*Jeśli dodatkowo <math>A</math> jest nieosobliwa, to macierz <math>A^{-1}</math> ma wartości własne <math>1/\lambda_k</math> oraz wektory własne <math>q_k</math>
-Łącząc te dwie własności mamy, że
-{{stwierdzenie|Transformacja widma macierzy||
-Macierz <math>(A-\sigma I)^{-1}</math> (o ile istnieje), to ma wartości własne równe <math>\frac{1}{\lambda_k - \sigma}</math> i wektory własne identyczne z <math>A</math>. }}
-Skoro tak, to jeśli najbliższą <math>\sigma</math>  wartością własną <math>A</math> jest <math>\lambda_j</math>, wówczas metoda potęgowa zastosowana do macierzy <math>(A-\sigma I)^{-1}</math> zbiegnie do <math>q_j</math>. To prowadzi do następującego algorytmu, odwrotnej metody potęgowej:
- $x_0$ = dowolny wektor startowy; k = 0;
- '''while''' ( !stop )
- {
-   $y_k$ = $Ax_{k-1}$;
-   $x_k$ = $y_k/||y_k||_\infty$;
-   k++;
- }
- $x_0$ = dowolny wektor startowy; k = 0;
- while( !stop )
- {
-   Rozwiąż układ równań $(A-\sigma I)y_k = x_{k-1}$;
-   $x_k$ = $y_k/||y_k||_\infty$;
-   k++;
- }
-Z własności metody potęgowej, metoda odwrotna potęgowa jest zbieżna tym szybciej, im bliżej <math>\lambda_j</math> jest przesunięcie <math>\sigma</math> (w stosunku do pozostałych wartości własnych). Dlatego dobrze byłoby --- dla zwiększenia szybkości zbieżności iteracji --- poprawiać wartość przesunięcia <math>\sigma</math>, korzystając z dotychczas wyznaczonego wektora <math>x_k \approx q_j</math> i ilorazu Rayleigh:
-<span id=""/> <math>
-\lambda_j = \frac{q_j^TAq_j}{q_j^Tq_j} \approx \frac{x_k^TAx_k}{x_k^Tx_k} </math>
- $x_0$ = dowolny wektor startowy; k = 0;
- '''while''' ( !stop )
- {
-   $y_k$ = $Ax_{k-1}$;
-   $x_k$ = $y_k/||y_k||_\infty$;
-   k++;
- }
- $x_0$ = dowolny wektor startowy; k = 0;
- while( !stop )
- {
-   Rozwiąż układ równań $(A-\sigma I)y_k = x_{k-1}$;
-   $x_k$ = $y_k/||y_k||_\infty$;
-   k++;
- }
- $x_0$ = dowolny wektor startowy; $\sigma_0$ = przybliżenie $\lambda_j$; k = 0;
- while( !stop )
- {
-   Rozwiąż układ równań $(A-\sigma_k I)y_k = x_{k-1}$;
-   $x_k$ = $y_k/||y_k||_2$;
-   $\sigma_{k+1}$ = $x_k^TAx_k$;
-   k++;
- }
-(wybierając normowanie wektora <math>x</math> w normie euklidesowej upraszczamy co nieco algorytm).
-Wielką zaletą metody RQI jest jej szybkość zbiezności: kwadratowa gdy wartość własna jest pojedyncza, a nawet sześcienna w przypadku macierzy symetrycznej.
-Wadą metody RQI jest to, że na każdym jej kroku należy rozwiązywać układ równań z <em>inną</em> macierzą.
-{{uwaga|Gdy złe uwarunkowanie pomaga...||
-Przez pewien czas numerycy odnosili się do tej metody z rezerwą, twierdząc, i słusznie, że im lepszym przybliżeniem <math>q_j</math> będzie <math>\sigma_k</math>, tym bardziej rośnie uwarunkowanie <math>A-\sigma_k I</math>, a tym samym --- błąd numerycznego rozwiązywania układu z tą macierzą będzie coraz większy i metoda będzie tracić stabilność. Tymczasem okazuje się, że --- choć rzeczywiście tak jest --- wektor błędu ma kierunek praktycznie zgodny z kierunkiem poszukiwanego wektora <math>q_j</math>, a tym samym tylko <em>pomaga</em> w zbieżności metody! }}
-\rysunek{}{Secular equation}
-==%s==