Liniowe zadania najmniejszych kwadratów

Liniowe zadania najmniejszych kwadratów polega na znalezieniu $ x^*\in \mathbb{R}^n $, który minimalizuje

$$<br />
 \|Ax-b\|_2<br />
$$

dla danej macierzy $ A\in \mathbb{R}^{m,n} $ i wektora $ b\in \mathbb{R}^m $.

Zauważmy, że jeśli $ m=n $ i $ A $ jest nieosobliwa, to LZNK jest równoważne rozwiązaniu układu $ Ax=b $.

Funkcją octave'a rozwiązującą dowolne LZNK jest operator \.

Wywołanie

 x=A\b

powinno zwrócić rozwiązanie LZNK.

W przypadku niejednoznaczności LZNK rozwiązanie jest wybrane według odpowiedniego kryterium, o którym napiszemy poniżej w rozdziale  [link].

Przetestujmy ten operator na kilku bardzo prostych przykładach macierzy
wymiaru $ 3\times 2 $.

Pierwszy przykład z $ b $ będącym pierwszą kolumną:

 A=[1 1 1; 1 1 2]'
 b=[1 1 1]'
 x=A\b

Tutaj rozwiązanie jest jednoznaczne.

Weźmy macierz o niezerowym jądrze (przestrzeni zerowej), wtedy oczywiście dodanie wektora z jądra macierzy do $ x $ nie zmieni normy $ Ax-b $.

 A=[1 1 1; -1 -1 -1]'
 b=[1 1 1]'
 x=A\b

Jądro macierzy:

 y=null(A)
 norm(A*(x+y)-b,2)

Regularne LZNK

Jeśli macierz $ A $ jest kolumnami regularna (tzn. kolumny stanowią układ niezależny liniowo), to LZNK nazywamy regularnym i tylko wtedy ma ono jednoznaczne rozwiązanie.

Rozwiążmy najbardziej klasyczne zadanie najmniejszych kwadratów. Dla danych punktów $ (x_k,y_k)_{k=1}^m $ dopasuj krzywą postaci $ \sum_{j=1}^na_jf_j(x) $
dla znanych funkcji $ \{f_j\} $ tak, aby

$$<br />
  \min_{a_j}\sum_k |y_k-\sum_ja_jf_j(x_k)|^2<br />
$$

była minimalna.

To jest oczywiście LZNK z macierzą:

$$<br />
  A=(f_j(x_k))_{k,j}<br />
$$

i wektorem prawej strony $ b=(y_k) $.

mo-NALrLZNKdop1.jpg
Rozwiążmy to zadanie dla dowolnych punktów zadanych poprzez dwa pionowe wektory $ x=(x_k) $ i $ y=(y_k) $ oraz krzywych postaci $ a_1\cos(x)+a_2\sin(x) $:

 function a=dopasuj(x,y,N=length(x))
   A=[cos(x),sin(x)];
   a=A\y;
 end

Przetestujmy na punktach leżących na wykresie funkcji $ \sin $.

 x=linspace(0,2*pi,10)';
 y=sin(x);
 a=dopasuj(x,y)

Czyli $ a_1=0 $ i $ a_2=1 $ (tak, jak się spodziewaliśmy).

Teraz rozpatrzmy losowe punkty blisko wykresu $ \cos $, por. rysunek  [link]:

 y=cos(x)+0.5*(rand(10,1)-0.5);
 a=dopasuj(x,y)
 z=linspace(0,2*pi);
 plot(x,y,"+",z,a(1)*cos(z)+a(2)*sin(z),";...
 dopasowana krzywa;",z,cos(z),";cos;")

Teraz możemy np. dodać duże dodatkowe zaburzenie np. $ y_1 $ i zobaczyć jak zmieni się dopasowywana krzywa. Pozostawimy to jako zadanie.

Nieregularne LZNK i rozkład SVD

(#)

W programie Matematyki Obliczeniowej nie ma treści związanych z nieregularnym LZNK, niemniej my w tym skrypcie pokażemy jak można rozwiązywać takie LZNK w octave'ie.

W ogólności nie ma miejsca wtedy jednoznaczności rozwiązania - jak to już zauważyliśmy.

Do rozwiązywania nieregularnego LZNK służy rozkład macierzy SVD, czyli rozkład względem szczególnych wartości macierzy.

Rozkład SVD macierzy $ A $ (dowolnej) to

$$<br />
 A=U\Lambda V^T,<br />
$$

gdzie $ U,V $ to macierze ortogonalne, a $ \Lambda $ to macierz diagonalna
z nierosnącymi nieujemnymi wartościami na diagonali. Wszystkie te wartości z diagonali, które są nieujemne, oznaczmy przez $ \sigma_1,\ldots,\sigma_r $. Nazywamy je wartościami szczególnymi macierzy $ A $. Wyznaczone one są jednoznacznie, ale już sam rozkład SVD nie jest jednoznaczny.

Widać, że
\begin{eqnarray*}
\|Ax-b\|_2^2=\|U\Lambda V^Tx-b\|_2^2=\|\Lambda y- U^T b\|_2^2
\end{eqnarray*}
dla $ y=V^Tx $.
Jeśli więc znajdziemy rozwiązanie $ y $ LZNK dla macierzy diagonalnej $ \Lambda $ z prawą stroną
$ g=U^Tb $, to od razu $ x=Vy $.

Dalej otrzymujemy, że
\begin{eqnarray*}
\|\Lambda y- U^T b\|_2^2=\sum_{k=1}^r |\sigma_k y_k -g_k|^2 +
\sum_{k=r+1}^n |g_k|^2.
\end{eqnarray*}
Biorąc minimum po $ y_k $ widzimy, że rozwiązaniem tego LZNK jest dowolne $ y $ spełniające:

$$<br />
  y_k=\frac{g_k}{\sigma_k } \qquad k=1,\ldots,r.<br />
$$

$ y_{r+1},\ldots,y_n $ mogą przyjąć dowolne wartości.

Otrzymujemy więc, że rozwiązaniem wyjściowego LZNK jest podprzestrzeń:

$$<br />
    \left\{V\left(<br />
 \begin{array}{c}<br />
  w \\ z<br />
 \end{array}<br />
\right): z\in \mathbb{R}^{n-r}  \right\}=<br />
V\left(<br />
 \begin{array}{c}<br />
  w \\ 0<br />
 \end{array}<br />
\right)+\left\{V\left(<br />
 \begin{array}{c}<br />
  0\\ z<br />
 \end{array}<br />
\right): z\in \mathbb{R}^{n-r}  \right\}<br />
$$

dla $ w\in  \mathbb{R}^r $ takiego, że $ w_k=y_k=\frac{g_k}{\sigma_k } $.

Zbiorem rozwiązań jest podprzestrzeń afiniczna wymiaru $ n-r $.

Operator \ zwraca element zbioru rozwiązań o najmniejszej normie drugiej, czyli

$$<br />
V\left(<br />
 \begin{array}{c}<br />
  w \\ 0<br />
 \end{array}<br />
\right).<br />
$$

Do znalezienia rozkładu SVD macierzy służy funkcja SVD:

  [U, S, V] = svd(A)

gdzie $ U,V $ to macierze ortogonalne, $ S $ to macierz diagonalna z wartościami szczególnych macierzy $ A $ takie, że $ U*S*V^T=A $.
Można też wywołać tę funkcję jako:

  S = svd (A)

Wtedy funkcja zwróci tylko wartości szczególne macierzy w wektorze.

Przetestujmy tę funkcję:

 A=hilb(3);
 S=svd(A)
 [U, S, V] = svd(A)
 norm(U*S*V'-A,1)

Przetestujmy tę funkcję dla macierzy niekwadratowej:

 A=[1 1 1; 1 2 3];
 B=A';
 [U, S, V] = svd (A)
 norm(U*S*V'-A,1)
 [U, S, V] = svd (B)
 norm(U*S*V'-B,1)

Rozpatrzmy macierz kwadratową osobliwą:

 A=[1 1 ;2 2];
 [U, S, V] = svd (A)

Formalnie wynik jest fałszywy. Druga własność szczególna powinna być równa zero. Z uwagi na błędy zaokrągleń nie możemy się spodziewać tego, że otrzymamy zero, a wynik na poziomie $ 1e-16 $ możemy uznać za numeryczne zero.
Rząd macierzy nie jest numerycznie dobrze uwarunkowany, najmniejsze zaburzenie macierzy może zmienić rząd, a więc
badanie rzędu maciery w arytmetyce zmiennopozycyjnej jest właściwie niemożliwe.

Przetestujmy, jak działa operator \ oraz funkcja svd() zastosowane do rozwiązania nieregularnego LZNK. Na początek rozważymy LZNK z macierzą o maksymalnym rzędzie, ale o większej ilości kolumn niż wierszy:

 A=[1 1 1; 1 2 3]
 [U, S, V] = svd(A)
 s=[1;0;0];
 f=A*s;
 xx=A\f
 y=V(:,1:2)*(S(1:2,1:2)\(U'*f))
 norm(A*xx-f,2)
 norm(A*y-f,2)

Rozpatrzmy $ f $ nie będące w przestrzeni kolumn macierzy $ A $.

 f=rand(2,1);
 xx=A\f
 y=V(:,1:2)*(S(1:2,1:2)\(U'*f))
 norm(xx-y,2)

Zajmijmy się macierzą, która nie ma maksymalnego rzędu:

 A=[1 1 1 1; 2 2 2 2; 1 -1 2 3]'
 [U, S, V] = svd(A)
 s=[1;0;0];
 f=A*s;
 xx=A\f
 g=U'*f;
 x=V(:,1:3)*(S(1:3,1:3)\g(1:3))
 norm(A*xx-f,2)
 norm(A*x-f,2)

Znów rozpatrzmy $ f $ nie będące w przestrzeni kolumn macierzy $ A $:

 f=rand(4,1);
 xx=A\f
 g=U'*f;
 x=V(:,1:3)*(S(1:3,1:3)\g(1:3))
 norm(A*xx-f,2)
 norm(A*x-f,2)

Wynik jest inny niż tego oczekiwaliśmy. Przyjmijmy, że trzecia wartość szczególna wynosi zero:

 S(3,3)=0;
 x=V(:,1:2)*(S(1:2,1:2)\g(1:2))
 norm(A*xx-f,2)
 norm(A*y-f,2)
 norm(xx-x,2)

Otrzymaliśmy te same rozwiązania.

W praktyce pojawia się pytanie, jak ustalić epsilon takie, że
jeśli $ |\sigma_k|\leq \epsilon $, to przyjmiemy, że $ \sigma_k $ jest równa zero.
Jakiś próg musimy przyjąć, ale zawsze wrowadzi to pewien błąd do naszego przybliżonego rozwiązania LZNK.

Przy pomocy rozkładu SVD wyznacza się też współczynnik uwarunkowania macierzy. Porównajmy wyniki obliczania przybliżonego uwarunkowania macierzy otrzymane przez svd(A) i funkcję cond(A):

 H=hilb(15);
 S=svd(H);
 max(S)/min(S)
 cond(A)

Otrzymaliśmy ten sam wynik. Współczynnik uwarunkowania macierzy Hilberta dla $ N=15 $ jest równy $ 2.7878e+17 $.

W przypadku, gdy $ A $ jest symetryczna i dodatnio określona, to $ U=V $ w rozkładzie SVD tej macierzy, tzn. rozkład SVD jest równoważny
znalezieniu bazy złożonej z układu ortonormalnych wektorów własnych $ A $.

 H=hilb(4);
 [U, S, V] = svd(H);
 norm(U-V,1)

SAGE