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: @@
-''Uwaga: przekonwertowane latex2mediawiki;
+<!--
-prawdopodobnie trzeba wprowadzi� poprawki!''
+Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
+Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki
-==Wyznaczanie wektorów i wartości własnych==
+wprowadzi� przykry@mimuw.edu.pl
+-->
-Niech będzie dana rzeczywista kwadratowa macierz <math>\displaystyle A</math> wymiaru <math>\displaystyle N</math>. Wektorem własnym <math>\displaystyle x\in C^N</math> oraz
-odpowiadającą mu wartością własną <math>\displaystyle \lambda \in C</math> nazwiemy taką parę, dla której
+=Nadokreślone układy równań liniowych=
-<center><math>\displaystyle A x = \lambda x,
-</math></center>
-przy czym <math>\displaystyle x\neq 0</math>.
+{{powrot |Metody numeryczne | do strony głównej
+przedmiotu <strong>Metody numeryczne</strong>}}
-Zadanie wyznaczania wartości własnych i wektorów własnych macierzy ma bardzo
+Zajmiemy się zadaniem wygładzania liniowego,
-szerokie zastosowania w tak odległych do siebie dziedzinach jak np. analiza
+nazywanym też <strong>liniowym zadaniem najmniejszych kwadratów</strong>.
-odporności konstrukcji mechanicznych (wieżowce, mosty, wagony kolejowe) na
+Jest ono uogólnieniem zadania rozwiązywania kwadratowych układów
-wibracje, czy też rankingowanie stron internetowych w wyszukiwarce Google.
+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ę,
-{{przyklad|Odporność budynku na trzęsienie ziemi||
+<center><math>||b-Ax||_2</math></center>
-Rozważmy prosty układ mechaniczny opisujący, naturalnie w pewnym
+Jest to praktycznie bardzo często pojawiające się zadanie, a autorem pierwszego
-jedynie przybliżeniu, zachowanie się układu <math>\displaystyle N</math> ciężkich płyt połączonych ze
+rozwiązania był nie kto inny jak sam wielki Gauss.
-sobą relatywnie elatycznymi dźwigarami --- co może np. modelować konstrukcję
-wieżowca.
-Wiadomo, że jeśli częstotliwości drgań własnych tego wieżowca będą bliskie
+[[grafika:Gauss.jpg|thumb|right||Carl Friedrich Gauss<br>  [[Biografia Gauss|Zobacz biografię]]]]
-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 ''wartości własne'' pewnej
- macierzy
-wymiaru <math>\displaystyle 2N</math>,
-która powstaje ze współczynników  równania różniczkowego opisującego dynamikę
-tego układu.
-}}
-{{przyklad|Macierz Google'a||
+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.
-Podstawowy algorytm rankingowania stron WWW w [http://www.wikipedia.org/pagerank  wyszukiwarce Google]
+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.
-sprowadza się do znalezienia rzeczywistego ''wektora własnego'' <math>\displaystyle \pi</math> pewnej silnie
-rozrzedzonej macierzy <math>\displaystyle A</math> (gigantycznego rozmiaru, równego liczbie indeksowanych
-stron, czyli w chwili pisania tego tekstu około <math>\displaystyle 2.5\cdot 10^{10}</math> stron), odpowiadającego wartości własnej równej 1:
-<center><math>\displaystyle A \pi = \pi.
+==Dopasowanie krzywej minimalizującej błąd średniokwadratowy==
-</math></center>
-Współrzędne wektora <math>\displaystyle \pi</math>
+<div style="margin-top:1em; padding-top,padding-bottom:1em;">
-interpretuje się jako wartość rankingową kolejnych stron WWW. Aby wszystko miało
+<span  style="font-variant:small-caps;">Przykład</span>
-sens, współrzędne wektora muszą być  z przedziału [0,1].  Pewne
+<div class="solution" style="margin-left,margin-right:3em;">
-twierdzenia matematyczne i subtelny dobór macierzy <math>\displaystyle A</math> gwarantują, że taki
-wektor <math>\displaystyle \pi</math> zawsze istnieje i jest jedyny! Co więcej, wartość 1 jest
-dominującą wartością własną <math>\displaystyle A</math>, a to z kolei ma ważne znaczenie dla tzw.
-[[sec:metoda-potegowa|Uzupe�nij: metody potęgowej]] numerycznego wyznaczania takiego wektora.
-}}
-{{przyklad|Wyznaczanie miejsc zerowych wielomianu||
+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>,
-Jak wiadomo, wartości własne to miejsca zerowe wielomianu charakterystycznego
+<center><math>
-macierzy <math>\displaystyle P(\lambda) = \det(A - \lambda I)</math>. Zachodzi także fakt odwrotny, to
+  \frac{1}{m}\sum_{i=1}^m (f_i-w(t_i))^2</math></center>
-znaczy miejsca zerowe wielomianu są wartościami pewnej macierzy, np. miejsca
-zerowe wielomianu
-<center><math>\displaystyle p(\lambda) = p_1 \lambda^N + \ldots + p_N \lambda + p_{N+1}
-</math></center>
-są wartościami własnymi m.in. macierzy stowarzyszonej,
+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>\displaystyle A = \begin{pmatrix}
+<center><math>\sum_{i=1}^m\left(f_i-\sum_{j=1}^n c_jw_j(t_i)\right)^2
--p_2/p_1 	& -p_3/p_1 	& \cdots 	& -p_{N+1}/p_1\\
-		& 		& 		& 	\\
-		& 1		& 		&	\\
-		& 		& \ddots	&	\\
-		& 		& 		& 1
-\end{pmatrix}
 </math></center>
-Funkcja Octave'a <code>compan(p)</code> wyznacza macierz stowarzyszoną dla zadanego
+względem <math>c_j</math>, a więc do zadania wygładzania liniowego.
-wielomianu o współczynnikach w wektorze <math>\displaystyle p = [p_1,\ldots,p_N, p_{N+1}]^T</math>. Z tej
-macierzy korzysta następnie funkcja Octave'a <code>roots</code>, która  właśnie w taki
-sposób wyznacza pierwiastki wielomianów: jako wartości własne macierzy
-stowarzyszonej.
-}}
-{{przyklad|||
+Rzeczywiście, kładąc
-Praktyczne zadanie z macierzą symetryczną
+<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.
-W praktyce obliczeniowej spotyka się zazwyczaj kilka typów zagadnień:
+[[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]]
-* 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?)
-* Wyznaczenie najmniejszej co do modułu wartości własnej i wektorów jej
-odpowiadających (zauważmy, że to jest np. zadanie wyznaczenia ''jądra
-macierzy osobliwej'' --- wtedy wiemy a priori, że szukana najmniejsza co do modułu
-wartość własna to zero)
-* 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 (na przykład, w celu znalezienia
-wszystkich pierwiastków zadanego wielomianu)
-* Wyznaczenie wszystkich wartości i wektorów własnych (tzw. pełne
-zagadnienie własne)
-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.
-Ponieważ w zastosowaniach bardzo często pojawiają się macierze rzeczywiste
+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>?)
-symetryczne (powyższe przykłady pokazują, że nie tylko!) szczegółową analizę
-metod numerycznych ograniczymy do tego przypadku, gdyż wtedy zachodzi
-{{twierdzenie|o symetrycznym zadaniu włanym||
+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.
-Każda macierz rzeczywista symetryczna <math>\displaystyle A</math> wymiaru <math>\displaystyle N</math> ma rozkład
+</div></div>
-<center><math>\displaystyle A = Q\Lambda Q^T,
-</math></center>
-gdzie <math>\displaystyle Q\in R^{N\times N}</math> jest ortogonalna (tzn. <math>\displaystyle Q^TQ = I</math>), a jej kolumnami są
+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>.
-wektory własne <math>\displaystyle A</math>, natomiast <math>\displaystyle \Lambda\in
-R^N</math> jest diagonalna z
-wartościami własnymi <math>\displaystyle A</math> na diagonali:
-<center><math>\displaystyle \Lambda = \begin{pmatrix} \lambda_1 & & \\ & \ddots & \\ & &
+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ę:
-\lambda_N\end{pmatrix} .
-</math></center>
-}}
+<center><math>||b-F(x)||_2</math></center>
-===Uwarunkowanie zadania===
+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]].
-{{twierdzenie|Bauer-Fike||
+==Układ równań normalnych==
-Niech <math>\displaystyle A\in R^{N\times N}</math> będzie diagonalizowalna, to
+Niech <math>A</math> będzie daną macierzą o <math>m</math> wierszach i <math>n</math> kolumnach,
-znaczy dla pewnej macierzy <math>\displaystyle X</math> zachodzi
+<math>A\in R^{m\times n}</math>, taką, że
-<center><math>\displaystyle X^{-1}  A X = \begin{pmatrix}  \lambda_1 & & \\ & \ddots & \\ & &
+<center><math>m\,\ge\,n\,=\, \mbox{rank} (A),
-\lambda_N\end{pmatrix} ,
 </math></center>
-a więc (gdyż macierz po prawej stronie jest podobna do <math>\displaystyle A</math>) <math>\displaystyle \lambda_i\in C</math>,
+albo równoważnie, taką że jej wektory kolumny są liniowo
-<math>\displaystyle i=1,\ldots,N</math> są
+niezależne. Niech także dany będzie wektor <math>b\in R^m</math>.
-wartościami własnymi <math>\displaystyle A</math>. Rozważmy macierz zaburzoną <math>\displaystyle \widetilde{A}</math> i jakąś jej
+Jasne jest, że wtedy układ równań <math>A x= b</math> nie zawsze
-wartość własną <math>\displaystyle \widetilde{\lambda}</math>. Wtedy istnieje wartość własna <math>\displaystyle \lambda_j</math>
+ma rozwiązanie - mówimy, że układ jest <strong>nadokreślony</strong>.
-macierzy <math>\displaystyle A</math> taka, że
-<center><math>\displaystyle |\lambda_j - \widetilde{\lambda}| \leq  \mbox{cond} _2(X) ||A - \widetilde{A}||_2.
+<strong>Zadanie wygładzania liniowego</strong> polega na znalezieniu wektora
-</math></center>
+<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.
-}}
+<center><math>\| b\,-\,A x^*\|_2\,=\,\min_{ x\in R^n}
+    \| b\,-\,A x\|_2</math></center>
-Ponieważ dla rzeczywistej macierzy symetrycznej macierz przejścia <math>\displaystyle X</math> jest
+{{lemat|||
-ortogonalna,
+Zadanie wygładzania liniowego ma jednoznaczne
-<math>\displaystyle X^{-1} = X^T</math>, to mamy <math>\displaystyle  \mbox{cond} _2(X) = 1</math> i w konsekwencji zachodzi
+rozwiązanie <math>x^*</math>, które można scharakteryzować jako rozwiązanie układu równań
-{{wniosek|Wartości własne macierzy symetrycznej są doskonale uwarunkowane||
+<center><math>
+A^TA x\,=\,A^T\, b</math></center>
-Przy oznaczeniach jak [[thm:Bauer-Fike|Uzupe�nij: twierdzeniu Bauera-Fike'a]], jeśli
-dodatkowo założymy, że macierz <math>\displaystyle A</math> jest rzeczywista i symetryczna, to
-<center><math>\displaystyle \min_{j=1,\ldots,N}|\lambda_j - \widetilde{\lambda}| \leq ||A - \widetilde{A}||_2.
-</math></center>
 }}
-Z drugiej strony, dla macierzy niediagonalizowalnych, uwarunkowanie wartości
+Zauważmy, że jeśli macierz <math>A</math> jest kwadratowa, <math>m=n</math>, to
-własnych może być
+rozwiązaniem jest <math>x^*=A^{-1} b</math> i residuum jest zerem.
-dowolnie duże, co ilustruje poniższy
+Zadanie wygładzania liniowego jest więc uogólnieniem
+rozwiązywania kwadratowych układów równań liniowych.
-{{przyklad|||
+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>.
-<center><math>\displaystyle A_\epsilon = \begin{pmatrix}  a & 1 \\ \epsilon & a \end{pmatrix}
+Ma on jednak pewne wady. Mnożenie macierzy powoduje w <math>fl_\nu</math>
-</math></center>
+powstanie po drodze dodatkowych błędów, które mogą nawet
+zmienić rząd macierzy. Na przykład, dla macierzy
-Weźmy dla uproszczenia <math>\displaystyle a=0</math>.
+<center><math>A\,=\,\left(\begin{array} {cccc}
-Wartości własne <math>\displaystyle A_\epsilon</math> to zera wielomianu <math>\displaystyle p_\epsilon(\lambda) = \lambda^2 - \epsilon</math>,
+  &  1  &  1  &  1  \\
-zatem <math>\displaystyle \lambda_\epsilon = \pm \sqrt{\epsilon}</math> i w konsekwencji
+  \epsilon  \\
+        &\epsilon \\
-<center><math>\displaystyle |\lambda_\epsilon - \lambda_0| / ||A_\epsilon - A_0|| = \sqrt{\epsilon}/\epsilon
+        &     &\epsilon \\
-\rightarrow \infty,
+        &     &      &\epsilon \end{array} \right)
 </math></center>
-gdy <math>\displaystyle \epsilon \rightarrow 0^+</math>, a więc uwarunkowanie takiego zadania jest
+mamy
-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>\displaystyle A</math> dla
-ujemnego parametru <math>\displaystyle \epsilon</math> są zespolone!
-[[Image:MNeigencond.png|frame|200px|Zachowanie się wartości własnych macierzy <math>\displaystyle A</math> (z
+<center><math>A^TA\,=\,\left(\begin{array} {cccc}
-parametrem <math>\displaystyle a=1</math>) w otoczeniu <math>\displaystyle \delta = 0</math>]]
++\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.
-Bardziej spektakularny przykład pochodzi od Wilkinsona:
+Poniżej przedstawimy inną metodę rozwiązywania zadania
+wygładzania liniowego, która oparta jest na specjalnych
+przekształceniach zwanych odbiciami Householdera.
-{{przyklad|Perfidny wielomian Wilkinsona||
+==Odbicia Householdera==
-Niech
+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
-<center><math>\displaystyle p(\lambda) = (\lambda -1)(\lambda - 2) \cdots (\lambda - 20).
+<center><math>H\,=\,I\,-\,2 w w^T</math></center>
-</math></center>
-Zmiana współczynnika przy <math>\displaystyle \lambda^{19}</math> o <math>\displaystyle 10^{-7}</math> skutkuje presunięciem niektórych
+Zauważmy, że
-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.
-[[Image:MNwilkinson.png|frame|200px|Zera oryginalnego i lekko zaburzonego perfidnego wielomianu
+<center><math>H x\,=\, x\,-\,2( w^T x) w</math>,</center>
-Wilkinsona.]]
-Jak widzimy, zera bardzo mało zaburzonego wielomianu mogą stać się wyraźnie nie-rzeczywiste!
+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>.
-}}
+Odbicia Householdera są przekształceniami nieosobliwymi
+spełniającymi
-Jeśli chodzi o wektory własne, ich wrażliwość na zaburzenia macierzy jest
+<center><math>H^{-1}\,=\,H\,=\,H^T</math></center>
-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.
-===Lokalizacja wartości własnych===
+Rzeczywiście, ponieważ <math>w</math> ma normę jednostkową, mamy
-Jak okaże się za chwilę, czasem warto mieć ogólne rozeznanie o tym, gdzie ''z
+<center><math>H^2 \,=\, (I-2 w w^T)^2\,=\,
-grubsza'' leżą wartości własne danej macierzy <math>\displaystyle A</math>. W tym celu mogą być nam
+   I-4 w w^T+4 w( w^T w) w^T \,=\, I</math>,</center>
-pomocne dwa fakty:
-{{fakt|||
+oraz
-Dowolna wartość własna <math>\displaystyle \lambda\in C</math> macierzy <math>\displaystyle A</math> spełnia
-<center><math>\displaystyle |\lambda| \leq ||A||,
+<center><math>H^T\,=\,(I-2 w w^T)^T\,=\,I-2( w^T)^T w^T\,=\,I.
 </math></center>
-gdzie <math>\displaystyle ||A||</math> jest dowolną normą macierzową indukowaną przez normę wektorową.
+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,
-Rzeczywiście, skoro istnieje wektor <math>\displaystyle x\neq 0</math> taki, że <math>\displaystyle Ax = \lambda x</math>, to stąd
+<center><math>\|H x\|_2\,=\,\sqrt{(H x)^T(H x)}\,=\,
-<math>\displaystyle ||Ax||/||x|| = |\lambda|</math>, więc fakt powyższy wynika już z definicji normy
+    \sqrt{ x^T(H^TH) x}\,=\,\sqrt{ x^T x}\,=\,
-macierzy:
+    \| x\|_2</math></center>
-<center><math>\displaystyle ||A|| = \max_{y\neq 0}\frac{||Ay||}{||y||} \geq ||Ax||/||x||.
+Odbicia Householdera zastosujemy do przeprowadzenia danego wektora
-</math></center>
+<math>x\ne 0</math> na kierunek innego niezerowego wektora, powiedzmy
+<math>e</math>, tzn.
-Drugie twierdzenie jest równie proste w dowodzie, ale daje trochę więcej
+<center><math>H x\,=\,(I-2 w w^T) x\,=\,\alpha\, e</math></center>
-informacji o lokalizacji widma.
-{{twierdzenie|Gerszgorina||
+<!--
+[[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
-Wartości własne macierzy <math>\displaystyle A</math> leżą w sumie (teoriomnogościowej) dysków <math>\displaystyle K_i</math> na
+<center><math>w\,=\,\frac{ x-\alpha e}{2( w^T x)},
-płaszczyźnie zespolonej,
-<center><math>\displaystyle K_i = \{z \in C: |z - a_{ii}| \leq \sum_{j\neq i} |a_{ij}| \}, \qquad i =
-,\ldots N.
 </math></center>
-}}
+a ponieważ <math>\alpha=\pm\| x\|_2</math> i <math>\| w\|_2=1</math> to
-{{przyklad|Koła Gerszgorina||
+<center><math>w\,=\,\frac{ x\mp\| x\|_2 e}
+                {\| x\mp\| x\|_2 e\|_2}</math></center>
-Niech
+W szczególności, jeśli <math>e= e_1</math> jest pierwszym
+wersorem, powyższe wzory dają
-<center><math>\displaystyle A = \begin{pmatrix}
+<center><math>H\,=\,I\,-\,\frac{ u u^T}{\gamma},
-.08930   & 1.38209  & -1.00037  &  0.69355  &  2.32178 \\
-.14211  &  1.74696 &   1.68440 &   0.30664 &   1.26718 \\
-  -0.74620  &  2.02686 &  -0.68293 &   0.19684 &   0.35854 \\
-.83517  &  0.74987 &   1.71331 &   1.09765 &  -0.44321 \\
-.02132  & -2.62155 &   0.79247 &   1.11408 &   0.48076 \\
-\end{pmatrix}
 </math></center>
-[[Image:MNgershgorindisks.png|frame|200px|Lokalizacja wartości własnych macierzy <math>\displaystyle A</math> kołami Gerszgorina oraz zgrubna
+gdzie
-lokalizacja wewnątrz okręgu
-o promieniu równym <math>\displaystyle ||A||_1</math>. Dokładne wartości własne zaznaczone trójkącikami.]]
-}}
+<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>
-===Metoda potęgowa, odwrotna potęgowa, RQI===
+oraz
-Jak wiemy z algebry, nawet gdy <math>\displaystyle A</math> jest macierzą rzeczywistą, jej
+<center><math>\begin{align} \gamma &= \frac 12\| u\|_2^2\,=\,
-widmo może być zespolone! Analizując poniższe metody, będziemy zakładać, że poszykiwane wartości i wektory
+      \frac 1 2\Big((x_1\mp\| x\|_2)^2+\sum_{i=2}^m x_i^2\Big) \\
- własne <math>\displaystyle A</math> są
+  &= \frac 1 2 \Big(\sum_{i=1}^m x_i^2\,+\,\| x\|_2^2\,\mp\,
-rzeczywiste. Iterując na liczbach rzeczywistych nie mamy wszak szansy, by
+x_1\|x\|_2\Big) \,=\,\|x\|_2^2\,\mp\,x_1 \|x\|_2.
-dotrzeć do liczb zespolonych!...
+\end{align}</math></center>
-====Metoda potęgowa====
+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
-Przypuśćmy, że wartości własne macierzy <math>\displaystyle A\in R^{N\times N}</math> spełniają
+<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>
-<center><math>\displaystyle |\lambda_1| > |\lambda_2| \geq \ldots \geq |\lambda_N|,
+oraz <math>\gamma=\| x\|_2^2+|x_1|\,\| x\|_2</math>, czyli zawsze
-</math></center>
+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>.
- (to znaczy, istnieje dokładnie jedna ''dominująca'' wartość własna macierzy
+==Rozkład QR==
- <math>\displaystyle A</math>.
-Załóżmy także, że istnieje baza złożona z wektorów własnych <math>\displaystyle q_1,\ldots,q_N</math> tej
-macierzy (tak jest np. dla macierzy symetrycznej na mocy
-[[thm:symetric-eig|Uzupe�nij: twierdzenia o własnościach symetrycznego zadania
-własnego]]).
-Kierunek własny <math>\displaystyle q_k</math> jakiejś macierzy <math>\displaystyle A</math> ma taką własność, że poddany działaniu przekształcenia
+Odbić Householdera można użyć do rozkładu macierzy
-<math>\displaystyle A</math> wydłuża się <math>\displaystyle \lambda_k</math> razy, wobec tego, dowolny wektor <math>\displaystyle x\in R^N</math> poddany
+<math>A\in R^{m\times n}</math> na iloczyn ortogonalno-trójkątny.
-działaniu <math>\displaystyle A</math> najbardziej wydłuży się w kierunku <math>\displaystyle q_1</math>. Iterując tę procedurę,
-powinniśmy dostawać w wyniku wektory, w których coraz bardziej dominuje kierunek
-<math>\displaystyle q_1</math>. Formalnie, niech
-<center><math>\displaystyle x = \alpha_1q_1 + \ldots + \alpha_Nq_N,
+Niech <math>A=( a_1, a_2,\ldots, a_n)</math>, gdzie <math>a_j</math> są
-</math></center>
+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
-wtedy
+<center><math>A^{(1)}\,=\,( a^{(1)}_1,\ldots, a^{(1)}_n)
+   \,=\,(H_1 a_1,\ldots, H_1 a_n)</math>,</center>
-<center><math>\displaystyle Ax = A \left( \sum_i \alpha_iq_i \right) = \sum_i \alpha_i A q_i
+w której pierwsza kolumna <math>a^{(1)}_1</math> ma niezerową tylko
-=  \sum_i \alpha_i \lambda_i q_i
+pierwszą współrzędną. W następnym kroku wybieramy drugie
-</math></center>
+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
-i w konsekwencji
+<center><math>H_2\,=\,\left(\begin{array} {cccc}
+&  0^T \\
+& \bar H_2  \end{array} \right)</math></center>
-<center><math>\displaystyle A^kx = \sum_i \alpha_i \lambda_i^k q_i = \lambda_1^k\left(\alpha_1q_1 +
+Pomnożenie macierzy <math>A^{(1)}</math> z lewej strony przez <math>H_2</math> spowoduje
-\alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^kq_2 + \ldots  +
+teraz wyzerowanie drugiej kolumny macierzy pod elementem
-\alpha_N\left(\frac{\lambda_N}{\lambda_1}\right)^kq_N \right).
+<math>a^{(1)}_{2,2}</math>, przy czym pierwszy wiersz i pierwsza kolumna
-</math></center>
+pozostaną niezmienione. Postępując tak dalej <math>n</math> razy
+(albo <math>n-1</math> razy gdy <math>m=n</math>) otrzymujemy
-Ponieważ z założenia, że istnieje dokładnie jedna dominująca wartość własna,
+<center><math>H_nH_{n-1}\cdots H_2H_1A\,=\,R</math>,</center>
-<math>\displaystyle \left|\frac{\lambda_N}{\lambda_1}\right| < 1</math>, to wyrażenie w nawiasie dąży do
-<math>\displaystyle \alpha_1q_1</math> i w konsekwencji wektory <math>\displaystyle x_k = A^kx</math> dążą, gdy
-<math>\displaystyle k\rightarrow\infty</math>, do kierunku wektora własnego <math>\displaystyle q_1</math>, to znaczy wektora
-odpowiadającego dominującej wartości własnej <math>\displaystyle A</math>  (o ile tylko <math>\displaystyle \alpha_1
-\neq 0</math>).
-Szybkość zbieżności metody potęgowej jest liniowa, o współczynniku zależnym od
+gdzie <math>R\in R^{m\times n}</math> jest uogólnioną macierzą trójkątną
-stosunku <math>\displaystyle \lambda_2/\lambda_1|</math>. W patologicznym przypadku, gdy <math>\displaystyle |\lambda_1|
+górną, tzn. <math>r_{i,j}=0</math> dla <math>i>j</math>. Stąd, podstawiając
-\approx |\lambda_2|</math>, może więc okazać się, że metoda praktycznie nie jest
+<math>Q=H_1H_2\cdots H_n</math>, dostajemy rozkład macierzy na iloczyn
-zbieżna.
+ortogonalno-trójkątny
-W praktyce nie wyznaczamy wzorem <math>\displaystyle x_k = (A^k)\cdot x</math>, lecz raczej korzystamy z
+<center><math>
-metody iteracyjnej
+  A\,=\,Q\cdot R.
+</math></center>
-{{algorytm|Metoda potęgowa||
+Rzeczywiście, macierz <math>Q\in R^{m\times m}</math> jest ortogonalna, bo
-<pre>
-<math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; k <nowiki>=</nowiki> 0;
+<center><math>\begin{align} Q^{-1} &= (H_1H_2\cdots H_n)^{-1}\,=\,
-while( !stop )
+     H_n^{-1}\cdots H_2^{-1}H_1^{-1} \\
-{
+   &= H_n^T\cdots H_2^TH_1^T \,=\,
-	<math>\displaystyle y_k</math> <nowiki>=</nowiki> <math>\displaystyle Ax_{k-1}</math>;
+     (H_1H_2\cdots H_n)^T\,=\,Q^T.
-	<math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_\infty</math>;
+\end{align}</math></center>
-	k++;
-}
-</pre>}}
-Warunek normowania ma m.in. na celu zapobieżenie powstawania nadmiaru i
+Dyspunując rozkładem QR, zadanie wygładzania liniowego
-niedomiaru (gdy <math>\displaystyle |\lambda_1| < 1</math>, to <math>\displaystyle ||A^kx|| \rightarrow 0</math>, a gdy
+można rozwiązać następująco. Ponieważ mnożenie przez macierz
-<math>\displaystyle |\lambda_1| > 1</math>, to <math>\displaystyle ||A^kx|| \rightarrow \infty</math>). Przy okazji,
+ortogonalną nie zmienia normy drugiej wektora, mamy
-<math>\displaystyle ||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.
-Zazwyczaj jako warunek stopu wybiera się kryterium małej poprawki, <math>\displaystyle ||x_k -
+<center><math>\begin{align} \| r\|_2 &= \| b-A x\|_2\;=\;\| b-QR x\|_2 \\
-x_{k-1}|| \leq \epsilon</math>, lub warunek małego residuum, <math>\displaystyle ||Ax_k - \lambda_{1,k}
+    &= \|Q(Q^T b-R x)\|_2 \;=\;\| c-R x\|_2,
-x_k||\leq \epsilon</math>, gdzie <math>\displaystyle \lambda_{1,k}</math> jest przybliżeniem <math>\displaystyle \lambda_1</math>
+\end{align}</math></center>
-dostępnym na <math>\displaystyle k</math>-tej iteracji.
-[[Image:MNXXX.png|frame|200px|Zasada działania metody potęgowej]]
+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
-Metoda potęgowa doskonale sprawdza się, gdy macierz <math>\displaystyle A</math> jest macierzą
+<center><math>R\,=\,\left(\begin{array} {c} R_I \\ 0\end{array} \right)</math>,</center>
-rozrzedzoną --- np. w przypadku macierzy Google'a.
-====Odwrotna metoda potęgowa====
+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
-Zauważmy, że dla dowolnej macierzy kwadratowej <math>\displaystyle A</math> o wartościach własnych
+<center><math>\| r\|_2^2\;=\;\| c_I-R_I x\|_2^2\,+\,
-<math>\displaystyle \lambda_k</math> i odpowiadających im wektorach własnych <math>\displaystyle q_k</math>, mamy:
+     \| c_{II}\|_2^2.
-* Macierz <math>\displaystyle A-\sigma I</math> ma wartości własne <math>\displaystyle \lambda_k - \sigma</math> oraz wektory
+</math></center>
-własne <math>\displaystyle q_k</math>,
-* Jeśli dodatkowo <math>\displaystyle A</math> jest nieosobliwa, to macierz <math>\displaystyle A^{-1}</math> ma wartości
-własne <math>\displaystyle 1/\lambda_k</math> oraz wektory własne <math>\displaystyle q_k</math>
-Łącząc te dwie własności mamy, że
+Rozwiązanie <math>x^*</math> zadania wygładzania jest więc
+rozwiązaniem układu liniowego trójkątnego,
-{{stwierdzenie|Transformacja widma macierzy||
+<center><math>x^*\,=\,R_I^{-1} c_I</math>,</center>
-Macierz <math>\displaystyle (A-\sigma I)^{-1}</math> (o ile istnieje),
+oraz <math>\| r^*\|_2=\| b-A x^*\|_2=\| c_{II}\|_2</math>.
-to ma wartości własne równe <math>\displaystyle \frac{1}{\lambda_k - \sigma}</math> i wektory własne
-identyczne z <math>\displaystyle A</math>.
-}}
-Skoro tak, to jeśli najbliższą <math>\displaystyle \sigma</math>  wartością własną <math>\displaystyle A</math> jest <math>\displaystyle \lambda_j</math>,
+Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde
-wówczas metoda potęgowa zastosowana do macierzy <math>\displaystyle (A-\sigma I)^{-1}</math> zbiegnie do
+z kolejnych przekształceń Householdera <math>H_k</math> wyznaczamy przez
-<math>\displaystyle q_j</math>. To prowadzi do następującego algorytmu, odwrotnej metody potęgowej:
+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
-{{algorytm|Odwrotna metoda potęgowa||
+<center><math>(H_k x)_i\,=\,x_i\,-\,s\,u_{k,i},
-<pre>
+</math></center>
-<math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; k <nowiki>=</nowiki> 0;
+gdzie <math>s= u_k^T x/\gamma_k</math>.
-while( !stop )
-{
-	Rozwiąż układ równań <math>\displaystyle (A-\sigma I)y_k = x_{k-1}</math>;
-	<math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_\infty</math>;
-	k++;
-}
-</pre>}}
-====Metoda Rayleigh====
+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>)
-Z własności metody potęgowej, metoda odwrotna potęgowa jest zbieżna tym
+<center><math>\sum_{k=1}^n 4(m-k+1)(n-k)\,\approx\,\frac 43n^3+2n^2(m-n)
-szybciej, im bliżej <math>\displaystyle \lambda_j</math> jest przesunięcie <math>\displaystyle \sigma</math> (w stosunku do
+   \,=\,2n^2(m-n/3)
-pozostałych wartości własnych). Dlatego dobrze byłoby --- dla zwiększenia
-szybkości zbieżności iteracji --- poprawiać wartość przesunięcia <math>\displaystyle \sigma</math>,
-korzystając z dotychczas wyznaczonego wektora <math>\displaystyle x_k \approx q_j</math> i ilorazu
-Rayleigh:
-<center><math>\displaystyle \lambda_j = \frac{q_j^TAq_j}{q_j^Tq_j} \approx \frac{x_k^TAx_k}{x_k^Tx_k}
 </math></center>
-{{algorytm|Metoda RQI (Rayleigh Quotient Iteration)||
+operacji arytmetycznych i <math>n</math> pierwiastków kwadratowych. Zauważmy,
-<pre>
+ż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.
-<math>\displaystyle x_0</math> <nowiki>=</nowiki> dowolny wektor startowy; <math>\displaystyle \sigma_0</math> <nowiki>=</nowiki> przybliżenie <math>\displaystyle \lambda_j</math>; k <nowiki>=</nowiki> 0;
+===Implementacja===
-while( !stop )
-{
-	Rozwiąż układ równań <math>\displaystyle (A-\sigma_k I)y_k = x_{k-1}</math>;
-	<math>\displaystyle x_k</math> <nowiki>=</nowiki> <math>\displaystyle y_k/||y_k||_2</math>;
-	<math>\displaystyle \sigma_{k+1}</math> <nowiki>=</nowiki> <math>\displaystyle x_k^TAx_k</math>;
-	k++;
-}
-</pre>}}
- (wybierając normowanie wektora <math>\displaystyle x</math> w normie euklidesowej upraszczamy co nieco
+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>.
- algorytm).
-Wielką zaletą metody RQI jest jej szybkość zbiezności: kwadratowa gdy wartość
+Nie popełnijmy jednak częstego błędu, prostodusznie implementując to mnożenie (przykładowo, w Octave) jako
-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ń
+ <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>H = eye(length(u)) - (u*u') / <math>\gamma</math>;
-z ''inną'' macierzą.
+y = H*x;
+</pre></div>
-{{uwaga|Gdy złe uwarunkowanie pomaga...||
+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.
-Przez pewien czas numerycy odnosili się do tej metody z rezerwą,
+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:
-twierdząc, i słusznie, że im lepszym przybliżeniem <math>\displaystyle q_j</math> będzie <math>\displaystyle \sigma_k</math>, tym
-bardziej rośnie uwarunkowanie <math>\displaystyle 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>\displaystyle q_j</math>, a tym samym tylko ''pomaga'' w zbieżności metody!
-}}
-===Metody rozwiązywania pełnego zadania własnego===
+<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>
-Najszybszą obecnie znaną metodą rozwiązywania pełnego zadania własnego (to
+Stąd <strong>prawidłowa</strong> implementacja mnożenia przez macierz Householdera:
-znaczy znajdowania wszystkich wartości i wektorów własnych) macierzy ''symetrycznej'' jest metoda '''"dziel i rządź"'''.
-Dla macierzy niesymetrycznych, najbardziej dopracowanym i przetestowanym, a więc
+ <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre><math>\omega</math> = u'*x;
-godnym zaufania algorytmem jest '''metoda QR z przesunięciami'''
+y = x - <math>\frac{\omega}{\gamma}</math>*u;
-(wykorzystująca, jak łatwo się domyślić, rozkład QR macierzy). Metoda QR
+</pre></div>
-przewyższa także metodę dziel i rządź w przypadku symetrycznym, gdy wymiar
-macierzy jest mały (mniej więcej <math>\displaystyle N \leq 25</math>).
-Obie metody są oczywiście metodami iteracyjnymi, jednak przyjęło się nazywać je
-metodami bezpośrednimi, gdyż praktycznie zawsze potrzebują z góry ograniczonej
-liczby iteracji do tego, by zbiec do wyniku o (niemal) maksymalnej rozsądnej
-dokładności.
-Dla efektywności obu metod kluczowy jest ''preprocessing'' macierzy,
-pozwalający niezbyt wygórowanym kosztem <math>\displaystyle O(N^3)</math> operacji sprowadzić przez
-ortogonalne podobieństwo zadanie z
-macierzą gęstą <math>\displaystyle A</math> do zadania z macierzą Hessenberga (w przypadku
-niesymetrycznym)
-<center><math>\displaystyle \begin{pmatrix}
-* & * & * & *       & \cdots & * \\
-* & * & * & *       & \cdots & * \\
-  & * & * & *       & \cdots & * \\
-  &   & \ddots  & \ddots  & \ldots &  \vdots \\
-  &   &   & \ddots  & \ddots & * \\
-  &   &   &         &   *    & *
-\end{pmatrix}
-</math></center>
-bądź wręcz trójdiagonalną, gdy <math>\displaystyle A</math> była symetryczna.
-Każdą macierz kwadratową <math>\displaystyle A</math> da się sprowadzić do postaci Hessenberga sekwencją
-przekształceń postaci
-<center><math>\displaystyle
-A := Q_k A Q_k^T,
-</math></center>
-gdzie <math>\displaystyle Q_k</math> jest pewnym przekształceniem Householdera. Niech <center><math>\displaystyle
-A = \begin{pmatrix}
-d_1 & * & * & *       & \cdots & * \\
-a_1  & * & * & *       & \cdots & * \\
-a_2  & * & * & *       & \cdots & * \\
-\vdots  & * & \ddots  & \ddots  & \ldots &  \vdots \\
-\vdots  & * & * & \ddots  & \ddots & * \\
-a_{N-1}  & * & * &   *     &   *    & *
-\end{pmatrix}
-</math></center>
-i oznaczmy <math>\displaystyle a = [a_1,\ldots,a_{N-1}]^T</math>. Możemy
+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!
-wziąć na początek przekształcenie Householdera <math>\displaystyle \widetilde{Q}_1</math> takie, że
-<math>\displaystyle \widetilde{Q}_1a = c\cdot
-e_1</math>, gdzie <math>\displaystyle e_1 = [1,0,\ldots, 0]^T</math>. Wtedy
-<center><math>\displaystyle
-\begin{pmatrix}
-& \\
-  & \widetilde{Q}_1
-\end{pmatrix}
-\cdot A = \begin{pmatrix}
-d_1 & * & * & *       & \cdots & * \\
-c  & * & * & *       & \cdots & * \\
-  & * & * & *       & \cdots & * \\
-  & * & \ddots  & \ddots  & \ldots &  \vdots \\
-  & * & * & \ddots  & \ddots & * \\
-  & * & * &   *     &   *    & *
-\end{pmatrix}
-</math></center>
-To samo przekształcenie przyłożone z prawej strony zachowa pierwszą kolumnę i
+Jest to całkiem typowe w numeryce:
-w efekcie  nie zmieni struktury macierzy:
-<center><math>\displaystyle
-\begin{pmatrix}
-& \\
-  & \widetilde{Q}_1
-\end{pmatrix}
-\cdot A
-\cdot
-\begin{pmatrix}
-& \\
-  & \widetilde{Q}_1
-\end{pmatrix}
- = \begin{pmatrix}
-d_1 & * & * & *       & \cdots & * \\
-c  & * & * & *       & \cdots & * \\
-  & * & * & *       & \cdots & * \\
-  & * & \ddots  & \ddots  & \ldots &  \vdots \\
-  & * & * & \ddots  & \ddots & * \\
-  & * & * &   *     &   *    & *
-\end{pmatrix} .
-</math></center>
-Dalej stosujemy tę samą metodę do podmacierzy wymiaru <math>\displaystyle N-1</math>, itd. aż dochodzimy
+<blockquote  style="background-color: #fefeee; padding:1em;  margin-left,margin-right:2em;  margin-top,margin-bottom: 1em;">
-do macierzy Hessenberga.
+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>
-Gdy wyjściowa macierz <math>\displaystyle A</math> jest symetryczna, to z definicji, macierz wynikowa
+==Uwarunkowanie==
-<math>\displaystyle \begin{pmatrix}
-I & \\
-  & \widetilde{Q}_{N-2}\end{pmatrix}
-\cdots \begin{pmatrix}
-& \\
-  & \widetilde{Q}_1
-\end{pmatrix}
- A \begin{pmatrix}
-& \\
-  & \widetilde{Q}_1\end{pmatrix}
-\cdots \begin{pmatrix}
-I & \\
-  & \widetilde{Q}_{N-2}
-\end{pmatrix} </math> też jest symetryczna i jednocześnie
-Hessenberga --- a więc musi być trójdiagonalna! Ponadto, macierz wynikowa będzie
-miała te same wartości własne, co <math>\displaystyle A</math>; wektory własne macierzy <math>\displaystyle A</math> także można
-łatwo (jak?) odzyskać z wektorów własnych macierzy wynikowej.
-====Metoda dziel i rządź====
+Ł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.
-Jest to obecnie najefektywniejsza metoda rozwiązywania zagadnienia własnego
+Najpierw jednak musimy rozszerzyć pojęcie uwarunkowania macierzy na macierze prostokątne.
-macierzy symetrycznej wymiaru powyżej kilkudziesięciu. Omówimy w zarysie jej
-najprostszy wariant (obarczony pewnymi wadami, usuniętymi w wersji
-bibliotecznej --- <code>DSYEVD</code> w LAPACKu).
-Startując z symetrycznej macierzy <math>\displaystyle A</math> już w postaci trójdiagonalnej, łatwo
+{{definicja|Uwarunkowanie macierzy prostokątnej w normie euklidesowej|Uwarunkowanie macierzy prostokątnej w normie euklidesowej|
-widzieć, że "prawie" rozpada się ona na dwie mniejsze macierze trójdiagonalne:
-dokładniej,
-<center><math>\displaystyle
+Niech <math>\Sigma(A)</math> będzie zbiorem wartości własnych macierzy <math>A^TA</math>. Definiujemy
-\begin{pmatrix}
-a_1 & b_1 &       &  \\
-b_1 & a_2 & \ddots & \\
-    & \ddots & \ddots & b_{N-1} \\
-    &      &   b_{N-1} & a_N
-\end{pmatrix}
-=
-\begin{pmatrix}
-T_1 & \\
-    & T_2
-\end{pmatrix}
-+ b_{m} uu^T,
-</math></center>
-gdzie <math>\displaystyle T_1 = \begin{pmatrix}
+<center><math>\mbox{cond} _2(A) = \sqrt{\frac{\max\{\lambda: \lambda \in \Sigma(A)\}}{\min\{\lambda: \lambda \in \Sigma(A)\}}}</math></center>
-a_1 & b_1 &       &  \\
-b_1 & a_2 & \ddots & \\
-    & \ddots & \ddots & b_{m-1} \\
-    &      &   b_{m-1} & a_m - b_m
-\end{pmatrix}
-</math>, <math>\displaystyle T_2 = \begin{pmatrix}
-a_{m+1} - b_m & b_{m+1} &       &  \\
-b_{m+1} & a_{m+2} & \ddots & \\
-    & \ddots & \ddots & b_{N-1} \\
-    &      &   b_{N-1} & a_N
-\end{pmatrix}
-</math> są --- tak jak <math>\displaystyle A</math> --- macierzami
-trójdiagonalnymi i symetrycznymi (jako podmacierze <math>\displaystyle A</math>), tylko o połowę mniejszego
-wymiaru, gdy <math>\displaystyle m \approx N/2</math>. Natomiast <math>\displaystyle u = e_{m} + e_{m+1}</math>, tak więc macierz
-<math>\displaystyle b_{m} uu^T</math> ma tylko cztery niezerowe elementy, każdy równy <math>\displaystyle b_m</math>.
-Zgodnie ze swoją nazwą, metoda dziel i rządź sprowadza zadanie znajdowania par
+(Jeśli w mianowniku pojawiłoby się zero, kładziemy <math>\mbox{cond} _2(A) = +\infty</math>).
-własnych macierzy wymiaru <math>\displaystyle N</math> do dwóch takich zadań dla macierzy dwa razy
+}}
-mniejszych. Te z kolei można potraktować w taki sam sposób i iteracyjnie
-zmniejszyć wymiar macierzy do tak małego (około 25), by opłacało się zastosować
-metodę QR (teoretycznie, można byłoby oczywiście doprowadzić podział do momentu,
-gdy macierze trójdiagonalne są rozmiaru <math>\displaystyle 1\times 1</math> --- dla których rozwiązanie
-zadania włanego jest trywialne --- ale taki algorytm byłby bardziej kosztowny od
-wariantu z udziałem QR).
-Rzeczywiście, przypuśćmy, że dla obu macierzy trójdiagonalnych <math>\displaystyle T_1,T_2</math> umiemy
+Zauważmy, że jest to  rozszerzenie definicji zgodne z tym, co wcześniej definiowaliśmy dla macierzy kwadratowych.
-rozwiązać zadanie własne tak, że znamy macierze: <math>\displaystyle Q_i</math> --- ortogonalną oraz
-<math>\displaystyle D_i</math> --- diagonalną, takie, że
-<center><math>\displaystyle
-Q_i^T T_i Q_i = D_i \qquad i=1,2.
-</math></center>
-Wtedy łatwo widzieć, że dla łatwo wyznaczalnego wektora <math>\displaystyle v</math>,
+{{twierdzenie|O uwarunkowaniu zadania wygładzania liniowego|O uwarunkowaniu zadania wygładzania liniowego|
-<center><math>\displaystyle
+Niech <math>x</math> będzie rozwiązaniem zadania najmniejszych kwadratów dla niezerowej prawej strony <math>b</math>,
-\begin{pmatrix}
+<center><math>
-Q_1^T & \\
+||b-Ax||_2\rightarrow \min{} !
-      & Q_2^T
-\end{pmatrix}
-\left(
-\begin{pmatrix}
-T_1 & \\
-    & T_2
-\end{pmatrix}
-+ b_{m} uu^T
-\right)
-\begin{pmatrix}
-Q_1 & \\
-      & Q_2
-\end{pmatrix}
-=
-\begin{pmatrix}
-D_1 & \\
-      & D_2
-\end{pmatrix}
-+ b_{m} vv^T.
 </math></center>
-W ten sposób zadanie własne dla oryginalnej macierzy <math>\displaystyle T</math> wymiaru <math>\displaystyle N</math> jest
+i niech <math>\widetilde{x}</math> będzie rozwiązaniem zadania zaburzonego <center><math>
-równoważne zadaniu własnemu macierzy diagonalnej zaburzonej o macierz rzędu 1.
+||\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>
-Na szczęście łatwo pokazać, że jeśli <math>\displaystyle \lambda</math> nie jest wartością własną
+gdzie <math>\epsilon</math> jest dostatecznie małe.
-macierzy diagonalnej <math>\displaystyle D = \begin{pmatrix}
-D_1 & \\
-      & D_2
-\end{pmatrix} </math>, to  wartości własne <math>\displaystyle \lambda</math> macierzy <math>\displaystyle D+ b_{m} vv^T</math>
-spełniają równanie
+Oznaczmy
-<center><math>\displaystyle
+<center><math>\sin(\theta) = \frac{||b-Ax||_2}{||b||_2} < 1
-f(\lambda) \equiv 1 + b_{m} \sum_{j=1}^N\frac{v_j^2}{d_j - \lambda} = 0,
 </math></center>
-gdzie <math>\displaystyle d_j</math> są elementami na diagonali macierzy <math>\displaystyle D</math>.
+--- będzie to miara, jak bardzo jesteśmy w stanie zminimalizować resztę w oryginalnym zadaniu.
-[[Image:MNsecular.png|frame|200px|Wykres <math>\displaystyle f(\lambda)</math> dla macierzy jednowymiarowego
+Wtedy
-laplasjanu rozmiaru 10. Zwróć uwagę na asymptoty pionowe tej funkcji oraz jej
-monotoniczność.]]
-W typowym przypadku <math>\displaystyle f</math> będzie miała dokładnie <math>\displaystyle N</math> pojedynczych miejsc zerowych
+<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>
-i wykres zachęcający do stosowania do niej metody Newtona. Okazuje się, że
-ogólny przypadek nie jest istotnie trudniejszy, choć wymaga ważnych modyfikacji,
-zarówno w celu szybszego rozwiązywania powyższego równania nieliniowego, jak i w
-celu zapewnienia lepszej stabilności algorytmu.
-Ostateczny koszt wyznaczenia wszystkich wektorów i wartości własnych jest rzędu <math>\displaystyle O(N^3)</math> z
+}}
-małą stałą.
-====Metoda QR====
+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.
-Dla zadania własnego z macierzą niesymetryczną najczęściej stosuje się metodę
+{{wniosek|||
-QR.
+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>.
-Jakkolwiek ostateczna wersja metody <math>\displaystyle QR</math> działa dla macierzy niesymetrycznych,
-wygodnie będzie nam założyć dla przejrzystości ekspozycji, że macierz jest
-symetryczna i w konsekwencji ma rzeczywiste widmo.
-W najprostszym wariancie (bez przesunięć), algorytm QR ma postać:
-{{algorytm|Metoda QR||
-<pre>
-<math>\displaystyle A_1 = A</math>;
-for k <nowiki>=</nowiki> 1, 2, ...
-{
-	wykonaj rozkład <math>\displaystyle A_k = Q_{k}R_{k}</math>;
-	<math>\displaystyle A_{k+1} = R_{k}\cdot Q_{k}</math>;
-}
-</pre>}}
-Można sprawdzić, że <math>\displaystyle A, A_1, A_2,\ldots</math> mają te same wartości własne, bo <math>\displaystyle A_{k+1} =
-Q_{k+1}^TA_kQ_{k+1}</math>. Co więcej, powyższy algorytm (gdy <math>\displaystyle A</math> jest nieosobliwa) w
-zasadzie jest równoważny teoretycznemu algorytmowi iteracji prostej
-zastosowanemu nie do pojedynczego wektora, ale do <math>\displaystyle N</math> wektorów naraz:
-{{algorytm|Iteracja prosta na przestrzeni||
-<pre>
-<math>\displaystyle V_1 = I</math>;
-for k <nowiki>=</nowiki> 1, 2, ...
-{
-	<math>\displaystyle W_{k+1} = A\cdot V_k</math>;
-	wyznacz rozkład QR <math>\displaystyle W_{k+1} = V_{k+1} R_{k+1}</math>, gdzie <math>\displaystyle V_{k+1}</math> jest ortogonalna;
-}
-</pre>}}
-Drugi krok w istocie ortogonalizuje kolumny <math>\displaystyle W_{k+1}</math>. Gdyby nie ortogonalizować
-zestawu wektorów <math>\displaystyle W_{k+1}</math>, oczywiście dostalibyśmy w efekcie zbieżność
-wszystkich kolumn macierzy do tego samego wektora --- odpowiadającego
-dominującej wartości własnej <math>\displaystyle A</math>. Zapewniając sobie ortogonalność <math>\displaystyle V_{k+1}</math>,
-możemy liczyć na to, że kolejne kolumny macierzy <math>\displaystyle V_k</math> będą dążyć do wektorów
-własnych odpowiadających kolejnym wartościom własnym <math>\displaystyle A</math> (przy stosownych
-założeniach o <math>\displaystyle A</math>, m.in. że
-wszystkie wartości własne <math>\displaystyle A</math> spełniają <math>\displaystyle |\lambda_i| \neq |\lambda_j|</math> dla <math>\displaystyle i
-\neq j</math>). Jeśli założyć dla uproszczenia, że oba używane rozkłady QR mają
-jednoznacznie określone czynniki rozkładu (na przykład, wymuszając, by
-diagonalne elementy macierzy <math>\displaystyle R</math> były dodatnie) mamy zależności <math>\displaystyle V_{k+1} =
-Q_1\cdots Q_k</math> oraz <math>\displaystyle A_{k+1} = V_{k+1}^TAV_{k+1}</math>.
-Tak więc, w sprzyjających warunkach, metoda QR, jako równoważna
-iteracji prostej na podprzestrzeni, będzie zbieżna: <math>\displaystyle A_k \rightarrow A_\infty</math>,
-gdzie <math>\displaystyle A_\infty</math> jest macierzą trójkątną (bo wektory własne odpowiadające różnym
-wartościom własnym są ortogonalne), a tym samym wartościami własnymi <math>\displaystyle A_\infty</math>
-(a więc także <math>\displaystyle A</math>)
-będą liczby na diagonali <math>\displaystyle A_\infty</math>.
-{{twierdzenie|Zbieżność metody QR w szczególnym przypadku||
-Niech wartości własne <math>\displaystyle A\in R^{N\times N}</math> spełniają <math>\displaystyle |\lambda_1|,\ldots,
-|\lambda_N| > 0</math> oraz macierz <math>\displaystyle T = [x_1,\ldots,x_N]</math> o kolumnach <math>\displaystyle x_i</math> złożonych z
-kolejnych wektorów własnych <math>\displaystyle A</math> ma taką własność, że <math>\displaystyle T^{-1}</math> ma rozkład LU,
-<math>\displaystyle T^{-1} = LU</math>.
-Wtedy w metodzie QR, ciąg macierzy <math>\displaystyle Q_k</math> jest zbieżny do macierzy diagonalnej, a
-ciąg <math>\displaystyle A_k</math> ma podciąg zbieżny do macierzy trójkątnej, której elementy diagonalne
-<math>\displaystyle u_{ii}</math> są równe <math>\displaystyle \lambda_i</math> dla <math>\displaystyle i = 1,\ldots, N</math>.
 }}
-Powyższa wersja algorytmu QR jest mało praktyczna, m.in. jest zbieżna wolno i
+==Biblioteki==
-przy poważnych ograniczaniach na <math>\displaystyle A</math>. Sprytna modyfikacja algorytmu wyjściowego
-daje  w wyniku tzw. metodę QR z przesunięciami, która jest praktycznie
-niezawodna dla dowolnej macierzy.
-{{algorytm|Metoda QR z przesunięciami||
+W Octave, zadanie najmniejszych kwadratów rozwiązujemy praktycznie tak samo, jak równanie liniowe:
-<pre>
-<math>\displaystyle A_1 = A</math>;
+ <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>x = A \ b;
-for k <nowiki>=</nowiki> 1, 2, ...
+</pre></div>
-{
-	wybierz sprytnie przesunięcie <math>\displaystyle \sigma_k</math>;
+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.
-	wykonaj rozkład <math>\displaystyle A_k - \sigma_kI = Q_{k}R_{k}</math>;
-	<math>\displaystyle A_{k+1} = R_{k}\cdot Q_{k} + \sigma_kI</math>;
-}
-</pre>}}
-Koszt wyznaczenia wszystkich wektorów i wartości własnych jest rzędu <math>\displaystyle O(N^3)</math> ze
+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.
-stałą równą około 30.
-===Biblioteki===
+Funkcje biblioteczne rozwiązujące zadanie wygładzania liniowego są oczywistym składnikiem wszystkich szanujących się pakietów statystycznych.
-LAPACK zawiera w sobie kolekcję doskonałych narzędzi do rozwiązywania różnych
+==Literatura==
-wariantów zadania własnego, m.in. <code>DGEEV</code> dla macierzy niesymetrycznych
-oraz <code>DSYEV</code> dla macierzy symetrycznych rozwiązują pełne zagadnienie własne
-(wyznaczając wszystkie wartości własne i, opcjonalnie, wektory własne). Dla
-macierzy symetrycznych mamy jeszcze m.in. funkcje <code>DSYEVX</code> (dla wybranych
-wartości własnych) i <code>DSYEVD</code> (z algorytmem dzieli i rządź)
-Fortranowska biblioteka ARPACK rozwiązuje zadanie własne dla macierzy
+W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 5.3</b> w
-rozrzedzonych, znajdując kilka wybranych (np. największych co do modułu)
+* D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.
-wartości i wektorów własnych.
-Funkcja <code>eig</code> w Octave i MATLABie wyznacza wszystkie wartości własne (i
+Bardzo dużo na temat rozwiązywania liniowego zadania najmniejszych kwadratów można dowiedzieć się z książki
-opcjonalnie wektory własne) zadaniej gęstej macierzy --- oczywiście korzystając
+* <span style="font-variant:small-caps">A. Kiełbasiński, H. Schwetlick</span>, <cite>Numeryczna algebra liniowa</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa, 1992.
-z LAPACKa. Jak dotąd, tylko MATLAB potrafi skorzystać z ARPACKa dla wyznaczenia
-fragmentów widma macierzy rzadkiej, za pomocą funkcji <code>eigs</code>.