MN12

Nadokreślone układy równań liniowych

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 przy\-pa\-dku, gdy układ jest nadokreślony.

Jest to praktycznie bardzo często pojawiające się zadanie (pewien jego wariant rozwiązują np. nasze przenośne odbiorniki GPS), a autorem pierwszego rozwiązania był nie kto inny jak sam wielki Gauss.

Carl Friedrich Gauss
Zobacz biografię

Układ normalny

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

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

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

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

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

Przykład

Przypuśćmy, że dla pewnej funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:[a,b]\toR} obserwujemy jej wartości $f_i$ (dokładne lub zaburzone) w punktach $t_i$, $1\le i\le m$. Funkcję tą chcielibyśmy przybliżyć inną funkcją $w$ należącą do pewnej $n$ wymiarowej przestrzeni liniowej $W$, np. przestrzeni wielomianów stopnia mniejszego niż $n$. Jakość przybliżenia mierzymy wielkością

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

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

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

względem $c_j$, a więc do zadania wygładzania liniowego. Rzeczywiście, kładąc Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle A=(a_{i,j})\inR^{m\times n}} z $a_{i,j}=w_j(t_i)$, $b=(f_i{)}_{i=1}^m$ i $x=(c_j{)}_{j=1}^n$, wielkość (Uzupelnic: unorm ) jest równa $\Vert b-Ax{\Vert }_2^2$.

Lemat

Zadanie wygładzania liniowego ma jednoznaczne rozwiązanie $x^{*}$, które spełnia układ równań

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

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

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

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

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

mamy

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

Jeśli ${ϵ}^2<\nu$ to $f{l}_{\nu }(1+{ϵ}^2)=1$, co implikuje $\text{rank}(f{l}_{\nu }(A^{T}A))=1$, podczs gdy $\text{rank}(f{l}_{\nu }(A))=4$.

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

Odbicia Householdera

Dla danego wektora Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle w\inR^m} o normie $\Vert w{\Vert }_2=\sqrt{w^{T}w}=1$, odbicie (macierz) Householdera zdefiniowane jest jako

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

Zauważmy, że

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

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

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

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

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

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

oraz

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

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

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

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

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

<flash>file=householderidea.swf</flash><div.thumbcaption>Odbicie Househodera

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

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

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

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

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

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

gdzie

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

oraz

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \gamma &= \frac 12\| u\|_2^2\,=\, \frac 1 2\Big((x_1\mp\| x\|_2)^2+\sum_{i=2}^m x_i^2\Big) \\ &= \frac 1 2 \Big(\sum_{i=1}^m x_i^2\,+\,\| x\|_2^2\,\mp\, 2 x_1\|x\|_2\Big) \,=\,\|x\|_2^2\,\mp\,x_1 \|x\|_2. \endaligned}

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

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

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

Rozkład QR

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

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

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

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

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

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

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

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

$A=Q\cdot R.$

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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned Q^{-1} &= (H_1H_2\cdots H_n)^{-1}\,=\, H_n^{-1}\cdots H_2^{-1}H_1^{-1} \\ &= H_n^T\cdots H_2^TH_1^T \,=\, (H_1H_2\cdots H_n)^T\,=\,Q^T. \endaligned}

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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \| r\|_2 &= \| b-A x\|_2\;=\;\| b-QR x\|_2 \\ &= \|Q(Q^T b-R x)\|_2 \;=\;\| c-R x\|_2, \endaligned}

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

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

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

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

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

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

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

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

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

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

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

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

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