iloczyn macierzy \(\displaystyle{ A^{T}\cdot A.}\)
Dla macierzy małych rozmiarów - dobrze uwarunkowanych ta metoda jest efektywna - nie wprowadza dużych błędów zaokrągleń.
Dla problemów, w których występują macierze dużych rozmiarów, nawet macierzy rozrzedzonych nie jest to metoda efektywna.
Ponadto to liniowe rozwiązanie jest bardzo podatne na błędy zaokrągleń.
Aby się o tym przekonać rozważmy dla uproszczenia kwadratową, symetryczną, nieosobliwą macierz \(\displaystyle{ M }\) i obliczmy wskaźnik
uwarunkowania iloczynu \(\displaystyle{ M^{T}\cdot M.}\)
\(\displaystyle{ cond_{2}(M^{T}\cdot M) = \parallel M\cdot M^{T}\parallel_{2} \cdot \parallel M\cdot M^{-1}\parallel_{2} = \parallel M\parallel_{2}\cdot \parallel M^{-1}\parallel{2} = (cond_{2}(M))^{2}.}\)
Stąd, jeśli macierz \(\displaystyle{ M }\) jest nawet lekko źle uwarunkowana, to macierz \(\displaystyle{ M^{T}\cdot M }\) jest bardzo źle uwarunkowana.
Pokażemy sposób rozwiązania problemu \(\displaystyle{ MNK }\) bez wykonywania liniowego rozwiązania z potencjalnie źle uwarunkowaną macierzą metodą ortogonalizacji Gramma-Schmidta polegającą na przekształcaniu zbioru \(\displaystyle{ m }\) liniowo niezależnych wektorów w ortonormalny zbiór \(\displaystyle{ m }\) wektorów, takiego samego wymiaru.
Przypomnijmy z analizy wektorowej pojęcie rzutu wektora \(\displaystyle{ \vec{v} }\) w kierunku wektora \(\displaystyle{ \vec{u} }\)
\(\displaystyle{ rzut_{\vec{u}} (\vec{v}) = \frac{\vec{u}\cdot \vec{v}}{\parallel u \parallel^{2}}\vec{u}.}\)
Przekształcamy zbiór wektorów liniowo - niezależnych \(\displaystyle{ \{ \vec{a_{1}}, \vec{a_{2}}, \vec{a_{3}}, ... \} }\) na zbiór wektorów
ortogonalnych:
\(\displaystyle{ \vec{a_{1}} = \vec{a_{1}}, }\)
\(\displaystyle{ \vec{a_{2}} = \vec{a_{2}} - rzut_{\vec{a_{1}}}(\vec{a_{2}}) }\)
\(\displaystyle{ \vec{a_{3}} = \vec{a_{3}} - rzut_{\vec{a_{1}}}(\vec{a_{3}}) - rzut_{\vec{a_{2}}}(\vec{a_{3}}),}\)
..............................................
Każdy z tych wektorów normalizujemy:
\(\displaystyle{ \vec{q_{1}} = \vec{a_{1}}, }\)
\(\displaystyle{ \vec{q_{1}} = \frac{\vec{q_{1}}}{\parallel \vec{q_{1}}\parallel},}\)
\(\displaystyle{ \vec{q_{2}} = \vec{a_{2}}, }\)
\(\displaystyle{ \vec{q_{2}} = \vec{q_{2}}- (\vec{q_{1}}\cdot \vec{q_{2}})\cdot \vec{q_{1}}, }\)
\(\displaystyle{ \vec{q_{3}} = \vec{a_{3}}, }\)
\(\displaystyle{ \vec{q_{3}} = \vec{q_{3}} - (\vec{q_{2}}\cdot \vec{q_{3}})\vec{q_{2}} }\)
\(\displaystyle{ \vec{q_{3}} = \frac{\vec{q_{3}}}{\parallel \vec{q_{3}} \parallel },}\)
itd.
Otrzmaliśmy zbiór wektorów ortonormalnych
\(\displaystyle{ \{\vec{q_{1}}, \ \ \vec{q_{2}}, \ \ \vec{q_{3}}, \ \ ... \} }\)
Traktując ten zbiór wektorów jako kolumny macierzy \(\displaystyle{ A }\) można powiedzieć, że proces Grama-Schmidta przekształca ją w macierz
\(\displaystyle{ Q,}\) w której kolumny są ortonormalne mają taki sam wymiar jak kolumny macierzy \(\displaystyle{ A.}\)
Zauważmy, że jeśli \(\displaystyle{ A }\) jest macierzą o kolumnach \(\displaystyle{ \vec{a_{i}}, \ \ i =1,2,..., n, }\)
to każdy krok algorytmu Grama-Schmidta jest operacją wykorzystującą tylko kolumny macierzy \(\displaystyle{ A. }\)
Jest to dość prosty algorytm który podajemy poniżej w MATLAB.
Najpierw normalizujemy
\(\displaystyle{ \vec{q_{i}}, \ \ i =1,2,...,}\) po procesie ortogonalizacji, przed przejściem do następnej kolumny.
Następnie obliczamy mnożniki (rzuty skalarne)
\(\displaystyle{ r_{ij} = \vec{q_{i}}\cdot \vec{a_{j}}, }\) które pozwalają na uzyskanie iloczynu \(\displaystyle{ Q R. }\)
Proces rozkładu macierzy \(\displaystyle{ A }\) możemy przedstawić w postaci macierzowej:
\(\displaystyle{ A = Q\cdot R \ \ (1)}\)
Ponieważ rozwiązanie \(\displaystyle{ MNK }\) spełnia równanie:
\(\displaystyle{ A^{T}A \vec{x} = A^{T}\vec{b}, }\) uwzględniając równanie \(\displaystyle{ (1) }\) otrzymujemy:
\(\displaystyle{ R^{T}Q^{T}Q R\vec{x} = R^{T}Q^{T}\vec{b}, }\)
\(\displaystyle{ Q^{T}Q = I, }\)
\(\displaystyle{ R \vec{x} = Q^{T} \vec{b}. }\)
Rozwiązanie z macierzą \(\displaystyle{ R }\) jest szybkie i łatwe, ponieważ jest to macierz trójkątna - górna.
Ponadto otrzymujemy znacznie lepszy wskaźnik uwarunkowania zadania, bo \(\displaystyle{ cond_{2}(R) = cond_{2}(A), }\)
z równości \(\displaystyle{ Q^{-1}= Q^{T}, }\) i \(\displaystyle{ cond_{2}(Q) = 1. }\)
Rozwiązanie problemu \(\displaystyle{ MNK }\) jest już zadaniem dobrze uwarunkowanym, chyba że macierz \(\displaystyle{ A }\) jest macierzą źle
uwarunkowaną, ale to się zdarza rzadko.
Program w MATLAB
Kod: Zaznacz cały
function [Q,R] = gramschmidt_QR (A)
Q = A ;
n = size (Q ,2);
R = zeros (n , n );
for j = 1 : n
q = A (:,j);
for k = 1 : j -1
qk = Q (:,k );
R (k , j) = (qk'*q);
q =q -qk*R (k,j );
end
R (j,j)= norm (q,2);
Q (:,j)=q /R (j,j);
end
Kod: Zaznacz cały
>> A = rand (5 ,3)
A =
0.8147 0.0975 0.1576
0.9058 0.2785 0.9706
0.1270 0.5469 0.9572
0.9134 0.9575 0.4854
0.6324 0.9649 0.8003
>> Q = gramschmidt_QR ( A)
Q =
0.4927 -0.4806 -0.1780
0.5478 -0.3583 0.5777
0.0768 0.4754 0.6343
0.5523 0.3391 -0.4808
0.3824 0.5473 -0.0311
>> Q'*Q
ans =
1.0000 0.0000 -0.0000
0.0000 1.0000 -0.0000
-0.0000 -0.0000 1.0000