Metody numeryczne fizyki/Wyznaczanie wektorów własnych i wartości własnych dla dowolnej macierzy

Wyznaczanie wektorów własnych i wartości własnych dla dowolnej macierzy

---

Licencja

Autor: Mirosław Makowiecki Absolwent UMCS Fizyki Komputerowej Uniwersytetu Marii Curie-Skłodowskiej w Lublinie Email: miroslaw(kropka)makowiecki(małpa)gmail(kropka)pl Dotyczy: książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami. Użytkownika książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami nie zwalnia z odpowiedzialności prawnoautorskiej nieprzeczytanie warunków licencjonowania. Umowa prawna: Creative Commons: uznanie autorstwa, na tych samych warunkach, z możliwością obowiązywania dodatkowych ograniczeń. Autor tej książki dołożył wszelką staranność, aby informacje zawarte w książce były poprawne i najwyższej jakości, jednakże nie udzielana jest żadna gwarancja, czy też rękojma. Autor nie jest odpowiedzialny za wykorzystanie informacji zawarte w książce, nawet jeśli wywołaby jakąś szkodę, straty w zyskach, zastoju w prowadzeniu firmy, przedsiębiorstwa lub spółki bądź utraty informacji, niezależnie czy autor (a nawet Wikibooks) został powiadomiony o możliwości wystąpienie szkód. Informacje zawarte w książce mogą być wykorzystane tylko na własną odpowiedzialność.

---

---

Poznamy tutaj metody obliczania wektorów i wartości własnych dla macierzy o n wierszach i n kolumnach, czyli kwadratowej, a także poznamy metody obliczania wartości własnych dla tej macierzy znając definicję algebry, który można poznać ze standardowego kursu algebry.

Wstęp do obliczania wartości i wektorów własnych dla macierzy kwadratowej

• Wartości i wektory własne dla macierzy kwadratowej o stopniu n, czyli $\mathbf {A} _{n\times n}\;$, oznaczamy kolejno przez x, λ, wtedy te wielkości spełniają warunek Ax=λx.
• Wartości własne liczymy ze wzoru det(A-λI)=0, który przedstawia wielomian charakterystyczny, który przyrównamy tu do zera.
• Gdy rozpatrzymy macierz transponowaną A^T, to jego wektory i wartości własne nazywamy lewostronnymi wartościami i wektorami własnymi. Widmem macierzy $A_{n\times n}$, czyli Sp(A) nazywamy zbiór wartości własnych λ₁,λ₂,...,λ_n. Widma macierzy kwadratowej jego odpowiednika transponowanego są takie same.
• Jeśli przez λ_L oznaczymy wartości własne lewostronne, a przez λ_P wartości własne macierzy kwadratowej A, to wektory własne macierzy A i macierz transponowanej A^T są takie, że są do siebie prostopadłe według definicji iloczynu skalarnego, a dowód jego przebiega:<x_L,Ax_p>=λ_p(x_L,x_P), a także zachodzi <x_L,Ax_p>=<A^Tx_L,x_P=λ_L<x_L,x_P>, stąd odejmując te dwa ostatnie równania otrzymujemy (λ_L-λ_p)<x_L,x_P>=0.
• Jeśli wartości własne macierzy A są λ₁, λ₂,..,λ_n, to wartościami własnymi dla wielomianu p(t) przy t=A są liczby p(λ₁), p(λ₂),...,p(λ_n), to z definicji wspomnianego wielomianu wektory własne macierzy p(A) są takie same jak wektory własne macierzy A.
• Gdy λ jest wielomianem charakterystycznym macierzy A, który przyrównamy do zera, to A musi być macierzą zerową.
• Gdy zdefiniujemy macierz P, który jest macierzą podobieństwa i zachodzi P^-1AP=B, to macierze A i B są do siebie podobne.
• Macierz $\mathbf {Q} _{m\times n}\;$ i m≥n jest macierzą ortogonalną, gdy zachodzi warunek $\mathbf {Q} ^{T}\mathbf {Q} =I_{n\times n}\;$.
• Kolumny i wiersze macierzy ortogonalnej $\mathbf {Q} _{n\times n}\;$ są wektorami do siebie ortogonalnymi.
• Macierz diagonalna D jest macierzą podobną do macierzy A, gdy istnieje macierz ortogonalna Q, która jest macierzą ortogonalnie podobną, taką że zachodzi Q^TAQ=D.
• Macierze symetryczne A mają wektory własne o składowych rzeczywistych i rzeczywiste wartości własne.
• Macierz diagonalna D jest podobna do macierzy A na w sposób D=P^-1AP, które są do siebie podobne i mają takie same wartości własne zespolone.
• Każda macierz trójkątna, która jest w ogólności macierzą zespoloną, jest macierzą unitarnie podobną do dowolnej macierzy kwadratowej.

Zdefiniujmy macierz zwaną macierzą Jordana rzędu k, w której λ jest w ogólności liczbą zespoloną taką, że:

$\mathbf {J} ={\begin{bmatrix}\lambda &1&&0\\0&\ddots &\ddots &0\\0&0&\ddots &1\\0&0&&\lambda \end{bmatrix}}_{k\times k}\;$

(6.1)

• Jeśli macierze J₁,J₂,...,J_r będziemy nazywać klatkami Jordana (6.1), wtedy każdą macierz A przy definicji macierzy podobnej P, która w ogólności jest macierzą zespoloną, możemy przedstawić w postaci:

$\mathbf {P} ^{-1}\mathbf {A} \mathbf {P} ={\begin{bmatrix}\mathbf {J} _{1}&&&&&&\\&\mathbf {J} _{2}&&&&0&\\&0&&&&\ddots &\\&&&&&&\mathbf {J} _{r}\end{bmatrix}}\;$

(6.2)

• Kwadrat normy euklidesowej dla normy ||⋅||_E jest większy lub równy sumie kwadratów modułów wartości własnych macierzy A.

$\sum _{i=1}^{n}|\lambda _{i}|^{2}\leq ||\mathbf {A} ||_{E}^{2}\;$

(6.3)

• Widmo macierzy, które jest sumą macierzy A i skalaru pomnożonej przez macierz jednostkową jest równa sumie widma wspomnianej macierzy i skalaru c.

$\operatorname {Sp} (\mathbf {A} +c\mathbf {I} )=\operatorname {Sp} (\mathbf {A} )+c\;$

(6.4)

której powyższe wyrażenie jest zbiorem, którego elementy są sumą liczb λ_i i stałej "c", czyli λ₁+c,λ₂+c,...,λ_n+c.

• Koło domknięte o środku a_ii i promieniu równej sumie elementów na wierszu bez elementów na diagonali spełnia warunek:

$\mathbf {C} _{i}=\left\{z:|z-a_{ii}|\leq \sum _{\overset {j=1}{\underset {j\neq i}{}}}^{n}|a_{ij}|\right\}\;$

(6.5)

to wtedy widmo macierzy A zawiera się sumie zbiorów (6.5) dla wszystkich i od i=1 do n. Jeśli stworzymy zbiór dla k rozłącznych kół określonych według powyższego schematu, to w tak utworzonym zbiorze leży k wartości własnych należącej do macierzy A.

• Nierówność ρ(A)≥ρ(B) zachodzi, gdy moduły elementów macierzy B są mniejsze lub równe elementom macierzy A.
• Jeśli macierz A ma wartości własne λ, to macierz B=A-cx><x ma dokładnie takie same wartości i wektory własne jak macierz A dla wektorów własnych różnych od x, a dla wektora własnego równego x odpowiada jemu wartość własna równa λ-c.

Błędy zaokrągleń dowolnej macierzy wraz z jego wartościami i wektorami własnymi

Będziemy się tutaj zajmować błędami zaokrągleń macierzy A i jego wartości własnej, gdy błędy zaokrągleń naszej macierzy dążą do zera, czyli dla ${\overline {\lambda }}=\lambda +\delta \lambda \;$ błędy zaokrągleń wartości własnej dążą do zera, czyli dla macierze dokładnej ${\overline {A}}={\overline {A}}+\delta {A}\;$ przy zanikającym błędach zaokrągleń. Podajmy kilka definicji dotyczącej odległości pomiędzy zbiorami i odległości pomiędzy liczbą a zbiorem.

• Odległości pomiędzy widmem macierzy A a liczbą "z", czyli zbiorem wartości własnych Sp(A)={λ₁,λ₂,...,λ_n} nazywamy liczbę zdefiniowanej w postaci $\min _{i=1,..,n}|z-\lambda _{i}|\;$, które przedstawimy według oznaczenia dist(z-Sp(A)).
• Dla dwóch widm macierzy A i macierzy B, które są macierzami kwadratowymi, to odległość pomiędzy nimi przy zdefiniowanej wartości własnej dla tychże dwóch wspomnianych macierzy λ_i i γ_i^p jest ${\min _{p}\max _{i=1,...,n}|\lambda _{i}-\gamma _{i}^{(p)}|}\;$ dla kolejnych wskaźników i=1,2,..,n.
• Jeśli błędy zaokrągleń macierzy A dążą do zera, to widmo Sp(A+δA) dąży do widma Sp(A), wtedy powiemy ||Sp(A)-Sp(A)||_∞→0.
• Macierz A i macierz diagonalna D, jeśli są podobne do siebie, to zachodzi warunek P^-1AP=D, wtedy odległość wartości własnej macierzy przybliżonej A i macierzy dokładnej A+δA jest napisana poprzez oszacowanie dist(λ,Sp(A))≤α, przy której definicję stałej α jest opisywana poprzez
  α=||P^-1||₂||P||₂||δA||₂.
• Wskaźnik uwarunkowania macierzy A, która służy według zadania do wyliczania wartości własnej macierzy A podobnej do diagonalnej nazywamy liczbę K₂=||P^-1||₂||P||₂.
• Przy macierzy przekształcenia P i dowolnej macierzy symetrycznej A jest zadaniem dobrze uwarunkowanym, gdy zachodzi K₂(P)=1.
• Przekształcenie jednej macierzy w drugą nazywamy operację matematyczną w postaci formuły B=RAR^-1, która to macierz A jest macierzą nieosobliwa, co może to pogorszyć lub polepszyć uwarunkowanie wspomnianej macierzy, bo zachodzi K₂(P)=||P^-1||₂||P||₂, a dla macierzy drugiej podobnej do poprzedniej uwarunkowanie macierzy jest równe K₂=||(RP)^-1||₂||RP||₂.
• Gdy mamy przekształcenie B=QAQ^T przy macierzy ortogonalnej Q, to ono nie zmienia uwarunkowana wartości własnych macierzy, tzn. są spełnione warunki ||QP||₂=||P₂|| i ||(QP)^-1||₂=||P^-1Q^-1||₂=||P^-1||₂.
• Weźmy sobie dwie macierze A i D, która ta ostatnia macierz jest macierzą diagonalną, a wektory x_i i y_i dla i=1,...,n są odpowiednio wektorami własnymi i wektorami własnymi lewostronnymi wspomnianej macierzy, i oznaczając wartość własną macierzy A+δA przez γ_i, wtedy możemy napisać:

$\gamma _{i}=\lambda _{i}+\delta ^{1}\lambda _{i}+\delta ^{2}\lambda _{i}+o(||\delta \mathbf {A} ||^{2})\;$

(6.6)

wtedy poszczególne składniki w sumie (6.6) przepisujemy według dwóch schematów:

$\delta ^{1}\lambda _{i}={{\langle \mathbf {y} _{i},\delta \mathbf {A} \mathbf {x} _{i}\rangle } \over {\langle \mathbf {y} _{i},\mathbf {x} _{i}\rangle }}\;$

(6.7)

$\delta ^{2}\lambda _{i}={{1} \over {\langle \mathbf {y} _{i},\mathbf {x} _{i}\rangle }}\sum _{\overset {j=1}{\underset {i\neq i}{}}}^{n}{{\langle \mathbf {y} _{j},\delta \mathbf {A} \mathbf {x} _{i}\rangle } \over {\langle \mathbf {y} _{i},\mathbf {x} _{i}\rangle (\lambda _{i}-\lambda _{j})}}\;$

(6.8)

Dla dowolnego zaburzenia dla dowolnej macierzy widmo macierzy A γ∈Sp(A+δA) przy przekształceniu P^-1AP=J jest tak napisane, by było spełnione oszacowanie przy definicji α=K₂(P)||δA||₂, a także napiszemy wzór dla dostatecznie małych zaburzeń macierzy A poprzez oszacowanie:

$\operatorname {dist} (\gamma ,Sp(\mathbf {A} ))\leq k\max(\alpha ,{\sqrt[{k}]{\alpha }})\;$

(6.9)

$\operatorname {dist} (\gamma ,Sp(\mathbf {A} ))\leq (1+(k-1){\sqrt[{k}]{\alpha }}){\sqrt[{k}]{\alpha }}\;$

(6.10)

Przejdźmy teraz do dowodu wzoru (6.10), które przedstawimy przy pomocy macierzy Jordana. Weźmy sobie teraz dowolną macierz A, która jest podobna do macierzy Jordana J według przekształcenia P^-1AP=J. Weźmy sobie wartość własną macierzy A+δA, czyli γ, którego równanie własne przepisujemy wedle schematu (A+δA)x=γx, wtedy równanie własne macierzy Jordana jest (J+P^-1(δAP)y=γy, jeśli przyjmować będziemy E=P^-1(δA)P, wtedy możemy napisać (J-γI)y=-Ey, wtedy lewą stronę ostatnio wspomnianego równania przepisujemy w postaci rozwiniętej:

$(\mathbf {J} -\gamma \mathbf {I} )={\begin{bmatrix}\mathbf {J} _{1}-\gamma \mathbf {I} _{1}&&&\\&\ddots &\\&&\mathbf {J} _{1}-\gamma \mathbf {I} _{2}&&\\&&&\ddots &\\&&&&\mathbf {J} _{n}-\gamma \mathbf {I} _{3}&\\\end{bmatrix}}\;$

(6.11)

Poszczególne macierze występujące w macierzy (6.11), czyli J_m-γI_m, możemy przepisać według schematu:

${\begin{bmatrix}\lambda _{m}-\gamma &1\\0&\lambda _{m}-\gamma &1&&0\\&&\ddots &\ddots &\\&0&&&\lambda _{m}-\gamma \\\end{bmatrix}}\;$

(6.12)

Macierz odwrotną do macierzy prostej J-γI piszemy wychodząc od wzoru (6.12), tzn.:

$(\mathbf {J} _{m}-\gamma \mathbf {I} _{m})={\begin{bmatrix}{{1} \over {\lambda _{m}-\gamma }}&{{1} \over {(\lambda _{m}-\gamma )^{2}}}&\cdots &{{(-1)^{k_{m}-1}} \over {(\lambda _{m}-\gamma )^{k_{m}-1}}}\\0&{{1} \over {(\lambda _{m}-\gamma )}}&\cdots &{{(-1)^{k_{m}-2}} \over {(\lambda _{m}-\gamma )^{k_{m}-1}}}\\&&\ddots &\vdots \\&0&&{{1} \over {(\lambda _{m}-\gamma )}}\\\end{bmatrix}}\;$

(6.13)

Macierz odwrotną do macierzy (6.11) przy definicji poszczególnych jego kratek (6.13) piszemy według:

$(\mathbf {J} -\gamma \mathbf {I} )^{-1}={\begin{bmatrix}(\mathbf {J} _{1}-\gamma \mathbf {I} _{1})^{-1}&\\0&\ddots &\\0&(\mathbf {J} _{2}-\gamma \mathbf {I} _{2})^{-1}&\\&0&\ddots \\&&&(\mathbf {J} _{m}-\gamma \mathbf {I} )^{-1}\\\end{bmatrix}}\;$

(6.14)

Ze wzoru powyżej podanego możemy napisać równość y=-(J-γI)^-1Ey, wtedy z definicji normy możemy napisać nierówność, która wynika z ostatniej równości 1≤||(J-γI)^-1||₂||E₂, wtedy na podstawie tego wyrażenie (6.14) i $||\mathbf {M} ||_{2}\leq {\sqrt {||\mathbf {M} ||_{1}||\mathbf {M} _{\infty }||}}\;$, piszemy:

$||(\mathbf {J} -\gamma \mathbf {I} )^{-1}||_{2}\leq {{1} \over {|\lambda -\gamma |}}+{{1} \over {|\lambda -\gamma |^{2}}}+...+{{1} \over {|\lambda -\gamma |^{k}}}\;$

(6.15)

Jeśli oznaczymy ρ=dist(γ,Sp(A)), wtedy wyrażenie (6.14) przepisujemy wedle:

$\rho ^{k}\leq (1+\rho +\rho ^{2}+...+\rho ^{k-1})||\mathbf {E} ||_{2}\;$

(6.16)

Jeśli oznaczymy przez ρ^k=k||E||₂, wtedy na pewno dla ρ<1 powiemy $\rho \leq {\sqrt[{k}]{k}}{\sqrt[{k}]{||\mathbf {E} ||}}|_{2}\leq k{\sqrt {\alpha }}\;$, a jeśli ρ>1, wtedy możemy napisać 1+ρ+..+ρ^k≤kρ^k-1, co za tym idzie ρ<kα, wtedy na podstawie tego otrzymujemy ostateczny wzór (6.9). Gdy ρ jest dostatecznie małe, tzn. zachodzi warunek ρ<1, powiemy:

$1+\rho +\rho ^{2}+..+\rho ^{k-1}\leq 1+(k-1)\rho \leq \left(1+{{(k-1)\rho } \over {k}}\right)^{k}\;$

(6.17)

Możemy połączyć wzory (6.17) ze wzorem (6.15) i w ostateczności mamy (6.10). Weźmy sobie wartości własne λ₁,λ₂,...,λ_m dla macierzy A, a także mamy wartości własne γ_i dla macierzy dokładnej A+δA. Gdy zaburzenie dąży do zera, to wartości własne macierzy dokładnej i zaokrąglonej dążą do siebie.

• Niech wektor własny odpowiadający macierzy zaokrąglonej oznaczymy przez x_i, a wektor własny macierzy dokładnej oznaczmy przez y_i, wtedy minimalna odległość pomiędzy tymi wektorami oznaczymy przez:

$\rho (\mathbf {x} _{i},\mathbf {x} _{i})=\min _{x\in X_{j}}||\mathbf {y} _{i}-\mathbf {x} _{i}||\;$

(6.18)

Jeśli oznaczymy zbiór wartości własnych przez dim X_j=1, wtedy iloczyn skalarny pomiędzy nimi określamy poprzez kąt θ pomiędzy tymi wektorami, który możemy wyliczyć ze wzoru <y_i,x_i>/||x_i||⋅||y_i||=sinθ.

Twierdzenie
Weźmy sobie macierz symetryczną y, a także jej wartości własne λ_i, wtedy odległość pomiędzy wartościami własnymi określamy jako: $d=\min _{\overset {j=1,...,n}{\underset {\lambda _{j}\neq \lambda _{i}}{}}}\left|\lambda _{i}-\lambda _{j}\right|\;$, wtedy dla zaokrągleń wspomnianej macierzy δA, dla którego zachodzi ||δA||₂ dla której to macierzy dokładnej A+δA odpowiadają wartości własne γ_i, wtedy odległość (6.10) pomiędzy wektorami własnymi macierzy zaokrąglonej i dokładnej jest równa:

$\rho (\mathbf {y} _{i},\mathbf {x} _{i})\leq {{||\delta \mathbf {A} ||} \over {d-||\delta \mathbf {A} ||_{2}}}\;$

(6.19)

Dowód
Równanie własne macierzy dokładnej możemy przepisać wedle schematu (A+δA)y_i=γ_iy_i, co po przekształceniu mamy: (A-γ_ii)y_i=-δAy_i. Możemy napisać dolne oszacowanie pomiędzy wartościami własnymi macierzy dokładnej i zaokrąglonej |γ_i-λ_i|≥||δA||₂, wtedy na podstawie tego piszemy ||(A-γ_iI)y^⊥||₂≥(d-||δA||₂)||||y_ii^⊥||, co na podstawie tego powiemy: ||y_i||₂(d-||δA||₂)≥||(A-γ_iI)y^⊥||₂≥||y_i||₂ ≥||A||₂). Jeśli oznaczymy przez ρ(y_i,x_i), co na podstawie tego możemy napisać (6.19).

• Weźmy sobie macierz symetryczną A, która jest symetryczna, a także możemy napisać ||r||=ρ, to wtedy znajduje się co najmniej jedna wartość własna w przedziale $\langle {\overline {\lambda }}-\rho ,{\overline {\lambda }}+\rho \rangle \;$. Dla wartości własnej macierzy, które znajdują się poza tym przedziałem o odległośc co najmniej d>ρ, wtedy rysujemy ${\overline {\rho }}({\overline {\mathbf {x} }},\mathbf {x} )\geq {{\rho } \over {d}}\;$, co jest błędem wektora własnego ${\overline {\mathbf {x} }}\;$ obliczonego zgodnie z (6.18), przy którym X_j=X jest podprzestrzenią rozpiętą na wektorach własnych macierzy A, które odpowiadają wartością własnym z wyżej podanego przedziału.

Korzystając z definicji równania własnego, w którym poszczególne obliczenia są dokonane dokładnie, czyli dla wyrażenia Ax-λx=0, wtedy wektor reszt obliczamy, gdy poszczególne obliczenia w ostatnio wspomnianym równaniu są dokonane w przybliżeniem (zaokrągleniem) sposobem:

$\operatorname {fl} (\langle \mathbf {a} ^{(i)},{\overline {\mathbf {x} }}\rangle -\lambda {\overline {\mathbf {x} }})=a_{1}^{(i)}x_{1}(1+\epsilon _{1})(1+\epsilon _{2})(1+\epsilon _{3})+a_{2}^{(i)}x_{2}(1+\epsilon _{4})(1+\epsilon _{5})(1+\epsilon _{6})+...+\;$
$+a_{n}^{(i)}x_{n}(1+\epsilon _{3n+1})(1+\epsilon _{3n+2})(1+\epsilon _{3n+3})-\lambda {\overline {\mathbf {x} }}(1+\epsilon _{3n+4})(1+\epsilon _{3n+5})\;$

(6.20)

W obliczeniach (6.20) dokonaliśmy przybliżenia ${\overline {x}}_{j}=x_{j}(1+\epsilon _{j})\;$. Będziemy przyjmować, że długość słowa binarnego jest bardzo duża, tzn. zachodzi (n+2)⋅2^-t<0,1, wtedy na podstawie tego napiszemy:

$||\mathbf {r} _{i}||={\overline {\epsilon }}[(n+2)||\mathbf {A} ||_{\infty }+3|\lambda |]\leq {\overline {\epsilon }}(n+5)||\mathbf {A} ||_{\infty }\;$

(6.21)

Jeśli będziemy przyjmować ${\overline {\epsilon }}=1,06\epsilon \;$, i jeśli będziemy przyjmować, że $||\mathbf {r} ||_{2}\leq {\sqrt {n}}||\mathbf {r} ||_{\infty }\;$, wtedy na podstawie rozważań przedstawionym w tym zdaniu i wzoru (6.21) możemy napisać:

$||\mathbf {r} ||_{2}\leq {\overline {\epsilon }}\left(n^{{3} \over {2}}+5n^{{1} \over {2}}\right)||\mathbf {A} ||_{\infty }\;$

(6.22)

Gdy posiadamy zaokraglone wartości ${\overline {\lambda }}\;$ i ${\overline {x}}\;$ i obliczenia dokonujemy z użyciem podwójnej precyzji, wtedy dowiadujemy się przy oszacowaniu błędu przy liczeniu wyrażenia Ax-λx, że:

$||\mathbf {r} ||_{2}\leq {\overline {\epsilon }}\left(n^{{1} \over {2}}||\mathbf {A} ||_{\infty }+2n^{{1} \over {2}}|\lambda |\right)\leq 3{\overline {\epsilon }}n^{{1} \over {2}}||\mathbf {A} ||_{\infty }\;$

(6.23)

Znajdowanie miejsc zerowych dowolnej macierzy

W przypadku dowolnej macierzy wartości własne leżą w kole C_i o środku w punkcie a_ii i mającej promień równy promieniowi spektralnemu ρ(A). Dla dowolnej macierzy A zachodzi ρ(A)≤||A||_p przy którym p=1,2,∞,E. Wartości własne zatem należą do zbioru {z:|z|≤||A||_p}. Gdy p=∞ to oszacowanie wartości własnych piszemy przez |z|≤||A||_∞, ale ponieważ zachodzi Sp(A)=Sp(A^T), to również możemy powiedzieć |λ_i|≤||A||₁.

Wykorzystanie metod potęgowych przy wyznaczaniu poszczególnych wartości własnych i wektorów własnych dla dowolnej macierzy

Metoda potęgowa polega na wyznaczaniu wartości i wektorów własnych dowolnej macierzy A. Ta metoda jest bardzo skuteczna. Algorytm ten polega na obraniu dla i=0 dowolnego wektora początkowego x_i, którego norma dla ∞ jest równa jeden, tzn. ||x₀||_∞ i dalej ustalamy dalszy krok iteracji. Następnie obliczamy v_i+1=Ax_i, a także m_i+1=||v_i+1||_∞. Jesli m_i+1=0 przerywamy krok iteracji, jeśli nie, to dalszym krokiem iteracji jest:

$\mathbf {x} _{i+1}={{\mathbf {v} _{i+1}} \over {m_{i+1}}}\;$

(6.24)

Jeszli dalej i+1≥ITM, to wtedy algorytm zatrzymujemy, jeśli nie, to przechodzimy do etapu pierwszego liczenia wektora v_i+1. Okazuje się, że ciąg wektorów x₀,x₂,...,x_2j jest ciągiem zbieżnym do wektora x, który jest wektorem własnym równania własnego macierzy A, również ciąg {m_i} jest zbieżny do pewnej wartości m. Jeśli zachodzi ||x_j+1-x_j||→0, to wtedy mamy Ax=mx, a jeśli w przeciwnym wypadku ||x_j+1-x_j||_∞→2||x||, to wtedy zachodzi Ax=-mx. Dla tych przypadków λ=m lub λ=-m znajdujemy wektor własny w wyniku powyższych iteracji dla wektora x, którego norma jest równa jeden dla p=∞. W ogólnym przypadku dla uzyskanego ciągu wektorów własnych x_2i dla macierzy A, to może być on rozbieżny. Gdy jednak istnieje przekształcenie P, dla której macierz A jest podobna do macierzy diagonalnej λ₁,λ₂,..,λ_n, to wtedy możemy napisać twierdzenie:

Twierdzenie
Ogólnie rzecz mówiąc, jeśli mamy wartości własne λ₁, λ₂,..,λ_n, które są wartościami własnymi macierzy A, w których poszczególne elementy są różne, ale mają równe za to moduły, wtedy przy dowolnym wyborze elementu startowego x₀ uzyskany ciąg wektorów x₀, x₂, x₄,... jest ciągiem zbieżnym do pewnego wektora będącego wektorem własnym macierzy A.

Dowód
Dowód przeprowadźmy dla przypadku macierzy symetrycznej A, której wektory własne są wektorami bazy ortogonalnej w przestrzeni Rⁿ, której wektory bazy oznaczymy przez e₁, e₂, e₃,...,e_n, której odpowiadają wartości własne |λ₁|≥|λ₂|≥...≥|λ_n|. Oznaczmy wartość początkową wektora iteracji przez x₀=α₁e₁+α₂e₂+...+α_ne_n. Może się zdarzyć, że wszystkie współczynniki liniowości są równe zero oprócz jednej lub wartość własna może być wielokrotna. Weźmy sobie λ=0, wtedy możemy obrać tak x₀, by było można napisać Ax₀=0 dla ||x₀||=1. Natomiast, gdy λ≠0, to wtedy przyjmując oznaczenia

a=α_ie_i+...+α_i+ke_i+k, δ₀=α_i+k+1e_i+k+1+...+α_ne_n z powodów oczywistych otrzymujemy:

$\mathbf {A} ^{L}\mathbf {x} _{0}=\mathbf {A} ^{L}(\alpha _{i}\mathbf {e} _{i}+....+\alpha _{n}\mathbf {e} _{n})=\mathbf {A} (\mathbf {a} +\alpha _{i+k+1}\mathbf {e} _{i+k+1}+...+\alpha _{n}\mathbf {e} _{i+k+1})=\;$
$=\lambda ^{L}\mathbf {a} +\lambda _{i+k+1}^{L}\alpha _{i+k+1}\mathbf {e} _{i+k+1}+...+\lambda _{n}^{L}\alpha _{n}\mathbf {e} _{n}=\;$
$=\lambda ^{L}\left(\mathbf {a} +\underbrace {\left({{\lambda _{i+k+1}} \over {\lambda }}\right)^{L}\alpha _{i+k+1}\mathbf {e} _{i+k+1}+...+\left({{\lambda _{n}} \over {\lambda }}\right)^{L}\alpha _{n}\mathbf {e} _{n}} _{\mathbf {\delta } _{L}}\right)\;$

(6.25)

Przyjmując oznaczenia, co do δ₀, wówczas możemy napisać z założenia, że wektor δ₀ dąży do zera, a także że wartości własne macierzy A nie rosną wraz z rosnącym wskaźnikiem:

$||\mathbf {\delta } _{L}||_{2}=\left[\left({{\lambda _{i+k+1}} \over {\lambda }}\right)^{2L}\alpha _{i+k+1}^{2}+...+\left({{\lambda _{n}} \over {\lambda }}\right)^{2L}\alpha _{n}^{2}\right]^{{1} \over {2}}\leq \;$
$\leq \left|{{\lambda _{i+k+1}} \over {\lambda }}\right|^{L}\left(\alpha _{i+k+1}^{2}+...+\alpha _{n}^{2}\right)=\left|{{\lambda _{i+k+1}} \over {\lambda }}\right|^{L}||\mathbf {\delta } _{0}||_{2}\rightarrow 0\;$

(6.26)

Według obliczeń (6.25) i (6.26) od razu dochodzimy do wniosku, że zachodzą dwie tożsamości, tzn. ${{\mathbf {A} ^{L}\mathbf {x} _{0}} \over {\lambda ^{L}}}\rightarrow \mathbf {a} \;$ oraz ${{||\mathbf {A} ^{L}\mathbf {x} _{0}||_{\infty }} \over {|\lambda |^{L}}}\rightarrow ||\mathbf {a} ||_{\infty }\;$, wtedy na podstawie tego powiemy $\mathbf {y} _{L}={{\mathbf {A} ^{L}\mathbf {x} _{0}} \over {||\mathbf {A} ^{L}\mathbf {x} _{0}||_{\infty }}}(\operatorname {sqn} \lambda )^{L}\;$, że jest zawsze zbieżny do ${{\mathbf {a} } \over {||\mathbf {a} ||_{\infty }}}\;$, co jest wektorem własnym macierzy A dla wartości własnej λ. Oznaczmy wzór na y_L w postaci y_L=x_L(sqn λ), i wtedy mamy dla λ>0, że y_j=x_j, i dla λ<0 mamy y_2j+1=-x_2j+1. Dla pierwszego przypadku mamy ||x_i+1-x_i||→0, a dla drugiego przypadku dostajemy ||x_i+1-x_i||→2||x||.

Wykorzystanie metod Hessenberga przy wyznaczaniu poszczególnych wartości własnych dla dowolnej macierzy

Każdą macierz Hassenberga nazywamy macierz zapisaną w postaci sumy dwóch macierzy, tzn. macierzy trójdiagonalnej (5.93) T i macierzy trójkątnej górnej U, zdefiniowanej według H=T+U. Wiadomo jednak h_ij=0 dla j<i-1. Jeśli jednak potrafimy wyznaczyć wartości własne macierzy H, to również potrafimy wyznaczyć wartości własne dowolnej macierzy. Według prostego algorytmy QR każdą macierz A możemy rozłożyć na iloczyn macierzy ortogonalnej Q i macierzy trójkątnej R, którą zapisujemy według schematu A=QR. Gdy macierz A jest nieosobliwa, to poszczególne kolumnach macierzy Q możemy otrzymać dokonując ortogonalizacji macierzy A metodą Grama-Schmidta, wtedy kolumny macierzy R są zbudowane z współczynników rozwiniecie z ortogonalizowanej macierzy A. Zbudujmy teraz iterację, która jest algorytmem QR, która opisuje pewien ciąg macierzy H=H⁽ⁱ⁾,H⁽ⁱ⁾,..., która jest zbudowana w postaci dwóch wzorów:

$\mathbf {H} ^{(i)}=\mathbf {Q} ^{(i)}\mathbf {R} ^{(i)}\;$

(6.27)

$\mathbf {H} ^{(i+1)}=\mathbf {R} ^{(i)}\mathbf {Q} ^{(i)}\;$

(6.28)

We wzorach (6.27) i (6.28) wszystkie macierze Q⁽ⁱ⁾, a także H dla i=1,2,... są na pewno macierzami Hessenerga. Ogólnie rzecz biorąc powyższe dwa wzory możemy połączyć w jedną całość jako:

$\mathbf {H} ^{(i+1)}={\mathbf {Q} ^{(i)}}^{T}\mathbf {H} ^{(i)}\mathbf {Q} ^{(i)}\;$

(6.29)

Jeśli mamy wartości własne λ₁,λ₂,...,λ_n, które są parami sprzężone, tzn. mającej równe moduły, wtedy nie wszystkie elementy pod diagonalą dążą do zera. Macierz Hessenberga można otrzymać z ciągu macierzy H dla wskaźnika górnego górnego nieskończonego, wtedy otrzymujemy np. macierz:

$\mathbf {H} ^{(\infty )}={\begin{bmatrix}*&*&\cdot &\cdot &\cdot \\*&*&\cdot &\cdot &\cdot \\0&0&*&\cdot &\cdot \\0&0&0&*&*\\0&0&0&*&*\\\end{bmatrix}}\;$

(6.30)

We wzorze (6.30) gwiazdki (*) oznaczają elementy zbieżne, a kropki (⋅) oznaczają wartości dążące do ustalonych wartości. Mówiąc ogólnie macierz Hⁱ nie musi mieć ogólnie postaci (6.30), ale może mieć za to postać:

$\mathbf {H} ={\begin{bmatrix}0&0&0&1\\1&0&0&0\\0&1&0&0\\0&0&1&0\\\end{bmatrix}}\;$

(6.31)

Mając macierz (6.31) wykorzystując przy tym wzór (6.29), to w wyniku końcowych obliczeń otrzymujemy tożsamości H⁽¹⁾=H²=...=H⁽ⁱ⁾ przy macierzy H ortogonalnej, to w rozkładzie H=QR macierz R występująca w tymże ostatnim wzorze jest macierzą jednostkową. Zbieżność elementów występująca na diagonali w macierzy Hessenberga jest często powolna, aby temu zaradzić możemy dokonać przesunięcia wartości własnych, w ten sposób tworzymy ciąg macierzy H⁽¹⁾,H⁽²⁾,H⁽³⁾,..., które definiujemy przy pomocy dwóch wzorów:

$\mathbf {H} ^{(i)}-k_{i}\mathbf {I} =\mathbf {Q} ^{(i)}\mathbf {R} ^{(i)}\;$

(6.32)

$\mathbf {H} ^{(i+1)}=\mathbf {R} ^{(i)}\mathbf {Q} ^{(i)}+k_{i}\mathbf {I} \;$

(6.33)

Łącząc wzory (6.32) i (6.33) w jedną całość dostajemy końcowy iteracyjny wzór, który generuje elementy zbieżne do macierzy Hessenberga, którego schemat iteracyjny przepisujemy w postaci:

$\mathbf {H} ^{(i+1)}={\mathbf {Q} ^{(i)}}^{T}(\mathbf {H} ^{(i)}-k_{i}\mathbf {I} )\mathbf {Q} ^{(i)}+k_{i}\mathbf {I} ={\mathbf {Q} ^{(i)}}^{T}(\mathbf {H} ^{(i)}-k_{i}\mathbf {I} )\mathbf {Q} ^{(i)}+{\mathbf {Q} ^{(i)}}^{T}k_{i}\mathbf {I} \mathbf {Q} ^{(i)}={\mathbf {Q} ^{(i)}}^{T}\mathbf {H} ^{(i)}\mathbf {Q} ^{(i)}\;$

(6.34)

Gdy obierzemy odpowiednio k_i to otrzymujemy szybszą zbieżność do macierzy Hessenberga H⁽ⁱ⁾, współczynniki k_i dobiera je się w takiej postaci, by one były równe wartością własnym kolejnych elementu ciągu macierzy Hessenberga, tzn. $k_{i}=a_{nn}^{(i)}\;$, gdzie elementy k_i i k_i+1 są równe kolejnym elementom macierzy kwadratowej dwuwymiarowej usytuowanej w H⁽ⁱ⁾ w prawym dolnym rogu tejże wspomnianej macierzy. Ale te elementy mogą być również zespolone, wieć rachunki prowadzą do rachunku na macierzach zespolonych, wtedy należy wykonać dwie dalsze iteracje, ale ta metoda zawodzi, gdy mamy doczynienia z macierzami np. (6.31).

Znajdowanie wartości własnych dowolnej macierzy poprzez doprowadzenie jej do postaci Hessenberga

Istnieje wiele metod sprowadzania dowolnej macierzy A do macierzy Hessenberga poprzez wykorzystanie wiadomości o podobieństwie dwóch macierzy przy przekształceniu P. Istnieją trzy takie metody, tzn. metoda Householdera, metoda Givensa i na samym końcu metoda eliminacji Gaussa. Pierwsze dwie metody potrzebują ${{10} \over {3}}n^{3}+O(n^{2})\;$ i ${{5} \over {3}}n^{3}+O(n^{2})\;$ mnożeń w arytmetyce zmiennoprzecinkowej, dzięki której otrzymujemy macierz Hessenberga H, której błąd zaokrąglenia jest ||E||_E≤k₂εn²||A||_E dla pierwszej metody oraz ||E||_E≤k₁εn^1/2||A||_E dla drugiej metody. Stałe k₁ i k₂ występujące w dwóch ostatnich wzorach są oczywiście rzędu jedynki, które zależą od sposobu zaokrąglania w komputerze, czyli na maszynie cyfrowej. Natomiast liczba mnożeń w metodzie w eliminacji Gaussa jest ${{2} \over {3}}n^{3}+O(n^{2})\;$, którego ilość operacji jest istotnie mniejsza niż w poprzednich metodach, i ta metoda jest powszechnie stosowana w obliczeniach numerycznych przy wyznaczaniu wartości własnych dowolnej macierzy A. W tej ostatniej metodzie wykonujemy n-2 przekształceń uzyskując macierze A=A⁽¹⁾,A⁽²⁾,...,A^(n-1)=H, gdzie ta ostatnia macierz jest macierzą Hessenberga. Macierz A⁽ⁱ⁾ w której i-1 początkowych kolumn są również i-1 kolumnami macierzy Hessenberga H. Przekształcenie macierzy A ze wskaźnika górnego i do wskaźnika górnego i+1 postępujemy w taki sposób, by wybierać element o największym module z elementów $a_{i+1,i}^{(i)},a_{i+2,i}^{(i)},...,a_{n,i}^{(i)}\;$, a gdy tychże elementów jest więcej, to wybieramy element pierwszy z nich, a jeśli mamy doczynienia z zerowymi elementami, to powinno być A⁽ⁱ⁾=A⁽ⁱ⁺¹⁾, dalej wyznaczamy element A⁽ⁱ⁺²⁾. Jeśli element o największym module znajduje się na k-tym wierszu dla k≥i, to należy przedstawić k-ty wiersz i i+1-tą kolumną. Dalej liczymy wielkość $m_{j}^{(i)}={{a_{j,i}^{(i)}} \over {a_{i+1,i}^{(i)}}}\;$ dla j=i+2,i+3,..,n, które zgodnie z wyborem elementu maksymalnego zachodzi $|m_{j}^{(i)}|\leq 1\;$, co dzięki temu nasza metoda jest numerycznie stabilna. Następnie i+1-ty wiersz odejmujemy od j-tego wiersza pomnożonej przez $m_{j}^{(i)}\;$ i j-tą kolumnę dodajemy od i+1-tej kolumny pomnożonej przez $m_{j}^{(i)}\;$ dla j=1+2,i+3,..,n, wtedy przy takim postępowaniu uzyskana macierz A⁽ⁱ⁺¹⁾ jest podobna do macierzy A⁽ⁱ⁾ przy elementach $a_{i+2,i}^{(i+1)},a_{i+3,i}^{(i+1)},..,a_{n,i}^{(i+1)}\;$, które przy naszym postępowaniu są równe zero, co oznacza, że i-ta kolumna A⁽ⁱ⁺¹⁾ jest taka sama jak i-ta kolumna macierzy Hessenberga H.

Przegląd metod wyznaczania wektorów własnych metodą QR,LR i metodą iteracji odwrotnej

Niech macierz P jest macierzą przekształcenia, w którym macierz A jest podobna do macierzy B, jeśli potrafimy wyznaczyć wektory własne według równania Bx=λx, to wtedy wektor y=Bx jest wektorem własnym macierzy A. Jeśli założymy, że macierz B jest macierzą trójkątną górną, ale posiada różne wartości własne, wtedy wartość własna macierzy B odpowiada liczba λ_i=b_ii, a wektorowi własnemu wektor x⁽ⁱ⁾, które te wielkości wyliczamy z formuł:

$x_{j}^{(i)}=0{\mbox{ dla }}j=n,n-2,...,i+1{\mbox{ i }}x_{i}^{(i)}=1{\mbox{ i }}x_{j}^{(i)}=-{{\left(\sum _{k=j+1}^{i}b_{jk}x_{k}^{(i)}\right)} \over {b_{jj}-b_{ii}}}{\mbox{ j=i-1,i-2,...,1}}\;$

(6.35)

Wyznaczanie wektorów własnych metodą QR

W tej metodzie należy wyznaczyć macierze Q⁽¹⁾,Q⁽ⁱ⁾,..., a następnie znając macierz przekształcenia P wyznaczyć macierz Hessenberga H z macierzy A według P^-1AP=H. Dalej ze wzorów (6.27) i (6.28) lub ze wzorów (6.32) i (6.33) możemy wyznaczyć macierze H⁽ⁱ⁺¹⁾=Q^(i)TH⁽ⁱ⁾Q⁽ⁱ⁾, a także wyznaczamy kolejne macierze pomocnicze Z⁽ⁱ⁺¹⁾=Z⁽ⁱ⁾Q⁽ⁱ⁾ z warunkiem początkowym Z⁽¹⁾=P. Dalej musimy stwierdzić, czy macierz H^(N) ma elementy pod diagonalą dostatecznie bliskie zero, jeśli tak, to możemy przystąpić do wyznaczania wektorów własnych macierzy Hessenberga ze wzoru:

$\mathbf {y} ^{(k)}=\mathbf {Z} ^{(N)}\mathbf {x} ^{(k)}\;$

(6.36)

Wyznaczanie wektorów własnych metodą LR

Ta metoda jest bardzo podobna do poprzednio opisywanej metody, tutaj należy wyznaczyć ciąg macierzy A=A⁽¹⁾,A⁽²⁾,.., które liczymy z formuł dla i=1,2,..., które to L i U są odpowiednio macierzami trójkątnymi dolnymi i górnymi:

$\mathbf {A} ^{(i)}=\mathbf {L} ^{(i)}\mathbf {U} ^{(i)}\;$

(6.37)

$\mathbf {A} ^{(i+1)}=\mathbf {U} ^{(i)}\mathbf {L} ^{(i)}\;$

(6.38)

Kolejne macierze A⁽ⁱ⁾ możemy wyznaczyć ze wzoru A⁽ⁱ⁺¹⁾=(L⁽ⁱ⁾)^-1A⁽ⁱ⁾L⁽ⁱ⁾m, które są podobne do macierzy , który ten wzór wynika ze wzorów (6.37) i (6.38), gdzie macierze A⁽ⁱ⁾ dążą do macierzy trójkątnej. W celu wyznaczenia wektorów własnych dla macierzy A trzeba wyznaczyć obiekty Z⁽ⁱ⁺¹⁾=Z⁽ⁱ⁾L⁽ⁱ⁾ i Z⁽¹⁾=I. Aby otrzymać rozkład A⁽ⁱ⁾=L⁽ⁱ⁾U⁽ⁱ⁾ należy wykonać ${{1} \over {3}}n^{3}+O(n^{2})\;$ mnożeń, czyli trzy razy mniej niż przy rozkładzie A⁽ⁱ⁾=QR⁽ⁱ⁾, co sugeruje przewagę tejże metody nad metoda poprzednią, jednak stosując poprzednią metodę wykonujemy mniejszy błąd zaokrągleń, ale natomiast gdy macierz A jest macierzą symetryczną i dodatnio określoną to ta metoda jest stabilna numeryczna i szybsza od poprzednio opisywanej metody. Podobnie jak w poprzedniej metodzie stosuje się tutaj przesunięcie A⁽ⁱ⁾-k_iI, mając początkowe macierze A sprowadzane się je do macierzy Hessenberga, co zmniejsza liczbę obliczeń przy rozkładzie macierzy Hessenberga H=LU przy nakładzie pracy $_{{{1} \over {2}}n^{3}+O(n)}\;$.

Wyznaczanie wektorów własnych metodą iteracji odwrotnej

Metody QR, LR są bardzo podobne do metody potęgowej, dlatego że w nim następuje znikanie niektórych kolumn macierzy A na podstawie jego różnych wartości własnych, ale zbieżność w metodzie potęgowej jest jednak bardzo powolna, ale za to w tychże metodach przy właściwym wyborze przesunięcia k_i, która jest bliska wartości własnej λ_n, która jej dobrym przybliżeniem jest element $a_{nn}^{i}\;$ następuje zwiększenie zbieżności przy wyznaczaniu wektorów własnych. Szybkość zbieżności jest zależna od ${{|\lambda _{n}-k_{i}|} \over {|\lambda _{n-1}-k_{i}}}\;$, następuje wtedy zwiększenie szybkości zbieżności tejże metody. Wartościami własnymi macierzy A są zaś liczby λ₁, λ_n, ...,λ_n takie, że spełniają nierówności słabe |λ₁|≥|λ₂|≥....≥|λ_n|, to stąd wynika, że wartościami własnymi macierzy odwrotnej do A-k_iI są liczby ${{1} \over {\lambda _{1}-k_{1}}},{{1} \over {\lambda _{2}-k_{1}}},..,{{1} \over {\lambda _{n}-k_{1}}}\;$, wtedy można tak dobrać k_i, by uzyskać zbieżność wektorów zależnych od wskaźnika górnego "i" do wektora własnego macierzy A, czyli x z omówionym powyżej współczynnikiem. Omówmy metodę iteracji potęgowej. Weźmy i=0 i obierzmy taki wektor x₀ by był spełniony warunek ||x₀||=1 i ustalamy maksymalną wartość iteracji ITM. Dalej wyznaczamy wektor v_i+1=(A-k_iI)^-1x_i i skalar m_i+1=||v_i+1||_∞, a także również wektor x_i+1=v_i+1/m_i+1. Jeśli i+1≥ITM zatrzymujemy algorytm, a jeśli nie, to wtedy poprzedni krok wykonujemy dla i:=i+1. Liczby m₁, m₂,... są stosowane do wyznaczania wartości własnej λ_n, tak jak w metodzie potęgowej. Wartość k_i jest przybliżeniem wartości λ_n, gdy jednak jest ona dostatecznie bliska do niej, wtedy macierz (A-k_iI) jest oczywiście macierzą źle uwarunkowaną, co prowadzi do dużych błędów numerycznych przy wyznaczaniu wektora v_i+1, a jeśli nasz wektor własny jest dobrze uwarunkowany, wtedy zaokrąglenia przy wyznaczaniu układu równań zapisanej w postaci macierzowej (A-k_ii)v_i+1=x_i nie będą zbyt duże, bo musimy się zabezpieczyć się by wektor v_i+1 był taki by jego norma była dostatecznie duża. Weźmy sobie k_i=λ_n+e_n, by potem można było wyznaczyć wektor v_i+1, w której macierz A jest zapisana z pewnym zaokrągleniem, który zapisujemy w postaci układu równań zapisanej w postaci macierzowej (A+E-k_iI)v_i+1=x_i przy założeniu ||E||_∞=e. Napiszmy:

$(\mathbf {A} -\lambda _{n}\mathbf {I} )\mathbf {v} _{i+1}=\left(\mathbf {A} -(k_{i}-e_{n})\mathbf {I} +\mathbf {E} -\mathbf {E} \right)=\left(\mathbf {A} -k_{i}\mathbf {I} +\mathbf {E} +e_{n}\mathbf {I} -\mathbf {E} \right)\mathbf {v} _{i+1}=\;$
$=\mathbf {x} _{i}+(e_{n}\mathbf {I} -\mathbf {E} )\Rightarrow (\mathbf {A} -\lambda _{n}\mathbf {I} ){{\mathbf {v} _{i+1}} \over {||\mathbf {v} _{i+1}||_{\infty }}}=(e_{n}\mathbf {I} -\mathbf {E} ){{\mathbf {v} _{i+1}} \over {||v_{i+1}||_{\infty }}}+{{\mathbf {x} _{i}} \over {||\mathbf {v} _{i+1}||_{\infty }}}\Rightarrow \;$
$\Rightarrow ||(\mathbf {A} -\lambda _{n}\mathbf {I} )\mathbf {x} _{i+1}||_{\infty }\leq (|e_{n}|+e)+{{1} \over {||\mathbf {v} _{i+1}||_{\infty }}}\;$

(6.39)

Przy opisie (6.39) należy stosować taki układ równań w postaci macierzowej (A-k_iI)v_i+1=x_i, by norma błędu ||E||_∞ była mała, ale duży współczynnik uwarunkowania (A-k_iI) nie jest wtedy bardzo groźny, stąd wynika, że wektor reszt (A-λ_n)x_i+1 jest bardzo mały dla dobrze uwarunkowanego wektora x, stąd wynika, że dobrym przybliżeniem jest do niego x_i+1.

Ogólne metody rozkładu dowolnej macierzy na iloczyn macierzy QR

Podamy tutaj twierdzenia dotyczące rozkładu macierzy A na iloczyn macierzy QR, a później omówimy metodę rozkładu dowolnej macierzy Grama-Schmidta, Householdera i Givensa.

Twierdzenie
Dowolną macierz $\mathbf {A} _{m\times n}\;$ przy założeniu m≥n na k≤n niezależnych w nim liniowo kolumn, wtedy mamy rozkład A=QR, przy czym macierz $Q_{m\times k}\;$ jest macierzą ortogonalną, a macierz $\mathbf {R} _{k\times n}\;$ jest macierzą trapezoidalną górną, tzn. posiadająca elementy zerowe pod diagonalą (r_ij=0 dla i≤j), a jeśli mamy doczynienia z macierzą kwadratową A, to macierz R jest macierzą trójkątną górną. Należy pamiętać, że dowolną macierz $A_{m\times n}\;$ możemy rozłożyć w sposób jednoznaczny na iloczyn A=QR, to wtedy elementy na diagonali są zawsze dodatnie.

Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Grama-Schmidta

Poznamy tu rozkład dowolnej macierzy A na rozkład macierzy QR metodą ortogonalizacji Gramma-Schmidta, które zastosujemy dla naszej macierzy dla kolumn a₁, a₂,...,_n. Obierzmy sobie macierze A⁽¹⁾, A⁽²⁾,...,A⁽ⁿ⁾, które posiadają 1,2,...,.n kolumn macierzy A. Jeśli dodatkowo założymy, że pierwsza kolumna naszej macierzy jest niezerowa, to wtedy napiszmy wzory na Q⁽¹⁾ i na R⁽¹⁾, które wykorzystamy do liczenia macierzy A⁽¹⁾:

$\mathbf {Q} ^{(1)}=\left[{{\mathbf {a} _{1}} \over {||\mathbf {a} _{1}||_{2}}}\right]\;$

(6.40)

$\mathbf {R} ^{(1)}=\left[||\mathbf {a} _{1}||_{2}\right]\;$

(6.41)

$\mathbf {A} ^{(1)}=\mathbf {Q} ^{(1)}\mathbf {R} ^{(1)}\;$

(6.42)

Jeśli mamy już macierze początkowe, tzn. (6.40), (6.41) i (6.42), to wtedy możemy policzyć następne elementy tychże macierzy, tzn. Q⁽ⁱ⁺¹⁾, Rⁱ⁺¹, a na jej podstawie policzyć macierz A⁽ⁱ⁺¹⁾:

$\mathbf {Q} ^{(i+1)}=\left[\mathbf {Q} ^{(i)}\mathbf {q} _{i+1}\right]\;$

(6.43)

$\mathbf {R} ^{(i+1)}={\begin{bmatrix}\mathbf {R} ^{(i)}&\mathbf {r} _{i+1}\\0&\mathbf {s} _{i+1}\end{bmatrix}}\;$

(6.44)

$\mathbf {A} ^{(i+1)}=\mathbf {Q} ^{(i+1)}\mathbf {R} ^{(i+1)}=\left[\mathbf {A} ^{(i)}\mathbf {a} _{i+1}\right]\;$

(6.45)

We wzorach (6.43), (6.44) i (6.45) występują pewne nieznane wielkości, które to wyznaczymy wedle poniższych formuł dla i=1,2,...,n-1:

$\mathbf {r} _{i+1}=\left(\mathbf {Q} ^{(i)}\right)^{T}\mathbf {a} _{i+1}\;$

(6.46)

$\mathbf {p} _{i+1}=\mathbf {a} _{i+1}-\mathbf {Q} ^{(i)}\mathbf {r} _{i+1}\;$

(6.47)

$s_{i+1}=||\mathbf {p} _{i+1}||_{2}\;$

(6.48)

Jeśli s_i+1=0, to zachodzi równość Q⁽ⁱ⁺¹⁾=Q⁽ⁱ⁾, a jeśli mamy s_i+1≠0, to wtedy możemy napisać q_i+1=p_i+1/s_i+1. Łatwo udowodnić, że macierz (6.44) jest macierzą trapezoidalną, a macierz Q⁽ⁱ⁾ jest macierzą ortogonalną, a macierz A⁽ⁱ⁾ spełnia zależność (6.45). Wyniku użycia metody Grama-Schmidta jest możliwa utrata ortogonalności macierzy Q⁽ⁱ⁾ dla i=1,2,..,n w wyniku obliczeń numerycznych, które powstają wyniku zaokrąglenia wyników obliczeń. W tym celu należy przeprowadzić reortogonalizację p_i+1, ale to zwiększa nakład obliczeń numerycznych. Przy rozkładzie macierzy A=QR, która z założenia jest macierzą nieososbliwą należy wykonać $_{D=n^{3}-{{1} \over {2}}n^{2}-{{1} \over {2}}n}\;$ dodawań i $_{M=n^{3}+n}\;$ i "n" pierwiastkowań, które wykonujemy bez ponownej reorogonalizacji.

Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Householdera

Weźmy sobie macierz kwadratową H o k wierszach i k kolumnach przy zdefiniowanym wektorze u=z-α||z||₂e₁ i definicji liczby τ=||z||₂²-α(z,e₁):

$\mathbf {H} =\mathbf {I} -{{1} \over {\tau }}\mathbf {u} \rangle \langle \mathbf {u} \;$

(6.49)

Macierz zdefiniowana w punkcie (6.49) nazywamy macierzą Householdera. Podziałajmy macierzą (6.49) na wektor z=[z₁,z₂,...,z_k]^T przy wyżej podanych definicjach u i τ:

$\mathbf {H} \mathbf {z} =\mathbf {z} -{{1} \over {\tau }}\mathbf {u} (\mathbf {u} ,\mathbf {z} )=\mathbf {z} -\mathbf {u} {{||\mathbf {z} ||_{2}^{2}-\alpha (\mathbf {z} ,\mathbf {e} _{1})} \over {||\mathbf {z} ||_{2}^{2}-\alpha (\mathbf {z} ,\mathbf {e} _{1})}}=\mathbf {z} -\mathbf {u} =\mathbf {z} -(\mathbf {z} -\alpha ||\mathbf {z} ||_{2}\mathbf {e} _{1})=\alpha ||\mathbf {z} ||_{2}\mathbf {e} _{1}\;$

(6.50)

Aby pominąc błędy zaokrągleń spowodowane dzieleniem, która podowują nadmierny wzrost błędów numerycznych należy tak obrać odpowiednią liczbę τ, wtedy wprowadza się definicję stałej α w postaci: α=-sqn z₁, co powoduje α=-1 dla z₁≥0 i α=1 dla z₁<0. Macierz przekształcenia P⁽¹⁾=H⁽¹⁾ ustala pierwszą kolumnę $\mathbf {a} _{1}^{(1)}\;$ macierzy A zapisaną w postaci schematu: $\alpha ^{(1)}||\mathbf {a} _{1}^{(1)}||_{2}[1,0,...,0]_{1\times n}^{T}\;$ do macierzy przekształcenia. Zdefiniujmy sobie macierz powstałej z iloczynu lewostronnego macierzy przekształceń o wskaźnikach i=1,...,n-1 i macierzy A wedle definicji P^(n-1)P^(n-2)....P⁽²⁾P⁽¹⁾A=R. Zdefiniujmy macierz przekształcenia o wskaźniku górnym dwa jako:$\mathbf {P} ^{(2)}={\begin{bmatrix}1&0\\0&\mathbf {H} ^{(2)}\end{bmatrix}}\;$, wtedy macierz H⁽²⁾ stworzona jest na podstawie wierszy i kolumn macierzy P⁽¹⁾A o numerach i=2,..,n, w której pierwszą kolumną jest $\mathbf {a} _{1}^{(2)}\;$, którą jeśli podziałamy macierzą H, to wtedy mamy $\mathbf {H} ^{(2)}\mathbf {a} _{1}=\alpha ^{(n-1)}||\mathbf {a} _{1}^{(2)}||_{2}[1,0,...,0]_{1\times n-1}^{T}\;$. Mając $\mathbf {z} =\mathbf {a} _{1}^{(1)}\;$, to możemy podstawić za z do równania (6.49) i w ten sposób otrzymać ${\hat {H}}^{(2)}\;$. Idąc dalej

$\mathbf {P} ^{(n-1)}={\begin{bmatrix}1&&&&\\0&\ddots \\&&1&\\&&&\mathbf {H} ^{(n-1)}\end{bmatrix}}\;$

(6.51)

wtedy macierz H^(n-1) stworzona jest na podstawie wierszy i kolumn macierzy P^(n-2)...P⁽²⁾A o numerach n-1 i n, w której pierwszą kolumną jest $\mathbf {z} =\mathbf {a} _{1}^{(n-1)}\;$, którą podziałamy macierzą H, otrzymujemy ${\hat {H}}^{n-1}\mathbf {z} =\alpha ^{(n-1)}||\mathbf {a} _{1}^{(n-1)}||_{2}[1,0]^{T}\;$, mając z podstawiamy go do (6.49), w ten sposób otrzymujemy H^(n-1). Na podstawie powyższych rozważań możemy rozłożyć macierz A na iloczyn macierzy Q i R, czyli A=QR, wtedy poszczególne czynniki w tym iloczynie piszemy według jego definicji Q=P⁽¹⁾P⁽²⁾⋅...⋅P^(n-1) oraz R=Q^TA.

Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Givensa

Metoda Givensa jest bardzo podobna do metody Householdera. Weźmy sobie wektor $\mathbf {z} =[z_{1},z_{2},..,z_{k}]^{T}\;$, a wynik działania macierzy Givensa na wektor z przedstawiamy α||z||₂⋅[1,0,...,0]^T, tzn. dla α=±1 mamy Gz=α||z||₂e₁. Macierz Givensa definiujemy sposobem:

$\mathbf {G} ^{(j)}={\begin{bmatrix}1&\\0&\ddots &&&&0&\\&&c^{(j)}&s^{(j)}&\\&&s^{(j)}&-c^{(j)}&\\&0&&&&\ddots &\\&&&&&&1\\\end{bmatrix}}\;$

(6.52)

W macierzy (6.52) są tak odpowiednio dobrane liczby c^(j), s^(j), by podczas działania macierzy G^(j) na wektor pionowy $\mathbf {z} ^{(j)}=[z_{1}^{(j)},...,z_{j}^{(j)},z_{j+1}^{(j)},...,z_{n}^{(j)}]^{T}\;$, by j+1 współrzędna wyniku działania była równa zero. Opisane tutaj liczby wyznaczamy według schematów:

$\rho _{1}={\sqrt {\left(z_{j}^{(j)}\right)^{2}+\left(z_{j+1}^{(j)}\right)^{2}}}\;$

(6.53)

$\rho _{2}=\alpha \rho _{1}{\mbox{ dla }}\alpha =\pm 1\;$

(6.54)

${\begin{cases}c^{(j)}=z_{j}^{(j)}/\rho _{2}\\s^{(j)}=z_{j+1}^{(j)}/\rho _{2}\end{cases}}{\mbox{ dla }}\rho _{1}\neq 0\;$

(6.55)

${\begin{cases}c^{(j)}=1\\s^{(j)}=0\end{cases}}{\mbox{ dla }}\rho _{1}=0\;$

(6.56)

Napiszmy macierz G jako iloczyn macierzy Givensa o wskaźnikach górnych i=2,...,n, którego definicja jest G=G⁽²⁾⋅...⋅G^(n-1)G⁽ⁿ⁾. Podczas działania macierzy G⁽ⁿ⁾ na wektor pionowy n wyzerowuje jego element z_n, i tak idąc dalej do końca macierz G⁽²⁾, który podczas działania na wektor z wyzerowuje element z₂ do zera. Utwórzmy sobie macierz $\mathbf {G} _{1}=\mathbf {G} _{1}^{(2)}\mathbf {G} _{1}^{(3)}\cdot ....\cdot \mathbf {G} _{1}^{(n)}\;$, który podczas działania na wektor $\mathbf {a} _{1}^{(1)}\;$ (jest to pierwsza kolumna macierzy A) otrzymujemy wynik $\alpha ^{(1)}||\mathbf {a} _{1}^{(2)}||_{2}[1,0,...,0]_{1\times n}^{t}\;$. Weźmy sobie macierz G₂, który to działamy na wektor pionowy $\mathbf {a} _{1}^{(2)}\;$ (jest to pierwsza kolumna utworzonej macierzy o n-1 wierszach i kolumnach położonej w prawym dolnym rogu macierzy G₁A) otrzymujemy wynik $\alpha ^{(2)}||\mathbf {a} _{1}^{(2)}||\cdot [1,0,....,0]_{1\times (n-1)}^{T}\;$. Na samym końcu zbudujmy sobie macierz G_n-1 podczas działania na wektor $\mathbf {a} _{1}^{(n-1)}\;$, który jest macierzą o wymiarze 2x2 utworzoną z wierszy i kolumn macierzy G_n-2G_n-1...GA otrzymujemy wektor $\alpha ^{(n-1)}||\mathbf {a} _{1}^{(n-1)}||_{2}\cdot [1,0]^{T}\;$. Podczas rozkładu macierzy A na iloczyn Q i R otrzymujemy przy dokonanych obliczeniach przy powyższym schemacie otrzymujemy następujące ich definicje $\mathbf {Q} =(\mathbf {G} _{1}^{(n)}\mathbf {G} _{1}^{(n-1)}\cdot ....\cdot \mathbf {G} _{1}^{2})(\mathbf {G} _{2}^{(n)}\mathbf {G} _{2}^{(n-1)}\cdot ....\cdot \mathbf {G} _{2}^{(2)})\cdot ...\cdot \mathbf {G} _{n-1}^{(2)}\;$, a $\mathbf {R} =\mathbf {Q} ^{T}\mathbf {A} \;$. Liczba operacji potrzebna do wyznaczenia macierzy Givensa i macierzy R jest w postaci $M={{4} \over {3}}n^{3}-{{4} \over {3}}n\;$ mnożeń i $D={{2} \over {3}}n^{3}-{{1} \over {2}}n^{2}-{{1} \over {6}}n\;$ dodawań i ${{1} \over {2}}n^{2}-{{1} \over {2}}n\;$ pierwiastkowań.

Wyznaczenie wartości własnej dla symetrycznej trójdiagonalnej macierzy

Macierz trójdiagonalna jest definiowana wedle schematu (5.93), a dodatkowo jeśli jest symetryczna, to jego definicja jest:

$\mathbf {T} ={\begin{bmatrix}b_{1}&a_{1}&&0&\\a_{2}&b_{2}&a_{2}&&\\0&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &a_{n-1}\\&0&&a_{n}&b_{n}\end{bmatrix}}\;$

(6.57)

Omówimy tutaj dwie metody, tzn. metodę bisekcji i metodę QR jako metody wyznaczania odpowiednich wartości własnych.

Wyznaczanie wartości własnych dla trójdiagonalnej macierzy symetrycznej metodą bisekcji

Weźmy sobie dowolną liczbę λ, wtedy możemy napisać wielomian charakterystyczny ω(λ)=det(A-λI). Ponieważ co do założenia macierzy A wielomian ω(λ) możemy wyznaczyć ze związków rekurencyjnych zapisanej w postaci:

$\omega _{0}(\lambda )=1\;$
$\omega _{1}(\omega )=b_{1}-\lambda \;$
$\omega _{i}(\lambda )=(b_{i}-\lambda )\omega _{i-1}-a_{i}^{2}\omega _{i-2}(\lambda ){\mbox{ dla }}i=2,3,...,n\;$
$\omega (\lambda )=\omega _{n}(\lambda )\;$

(6.58)

• Jeśli wszystkie elementy a₂, a₃,...,a_n są niezerowymi wartościami, to macierz (6.57) ma wszystkie wartości własne różne od siebie, a także na podstawie (6.57) można wygenerować ciąg wartości zależnych od λ, tzn. ω₀(λ), ω₁(λ),...,ω_n(λ), które dla ω_i(λ), gdy i<n, spełniają warunki ω_i-1(λ)ω_i+1(λ)<0, lub gdy spełniony jest warunek ω_n(λ)≠0, to liczba zmian znaków opisywanego powyżej ciągu, dla i=1,..,n, jest równa liczbie wartości własnych macierzy (6.58). A ponadto gdy ω_n(λ)=0, to λ jest wartością własną macierzy (6.58), i tyle jest wartości własnych oczywiście mniejszych niż λ, ile nastąpiło zmian znaków ciągu dla naszego powyżej opisywanego ciągu dla i=0,..,n-1.

Powyższa metod chociaż jest metodą stabilną, ale poszczególne elementy ciągu ω₀(λ), ω₁(λ),...,ω_n(λ) mogą osiągać zbyt duże wartości, nawet dla małych liczb naturalnych "n", co następuje, gdy wartość λ odbiega zbytnio od wartości własnych macierzy (6.57), z tego względu ciąg generowany przez algorytm (6.58) przekształcamy w ciąg p_i: $p_{1}(\lambda )={{\omega _{1}(\lambda )} \over {\omega _{0}(\lambda )}},...p_{n}(\lambda )={{\omega _{n}(\lambda )} \over {\omega _{n-1}(\lambda )}}\;$, wtedy jego poszczególne elementy są generowane przez ciąg wynikającego z przekształconego algorytmu (6.58), jako:

$p_{1}(\lambda )=b_{1}-\lambda \;$

(6.59)

$p_{i}(\lambda )=(b_{i}-\lambda )-a_{i}^{2}/p_{i-1}(\lambda ){\mbox{ dla }}i=2,..,n\;$

(6.60)

Gdy ciąg generowany (6.60) ma element p_i-1 równy zerowy, to ten nasz algorytm musi ulec modyfikacji. Wyznaczmy dokładną wartość wyrażenia p_i(λ), który to przepisujemy znając zaokrągloną wartość b_i, a_i i p_i-1, na w sposób:

$\operatorname {fl} (p_{i}(\lambda ))=b_{i}(1+\epsilon _{1})(1+\epsilon _{3})-\lambda (1+\epsilon _{2})(1+\epsilon _{3})-{{a_{i}^{2}} \over {p_{i-1}(\lambda )}}(1+\epsilon _{4})(1+\epsilon _{5})(1+\epsilon _{6})\;$

(6.61)

Liczba zmian znaków elementu dokładnego fl(p_i(λ)) jest równa liczbie zmian znaków (6.61) podzielonej przez (1+ε₂)(1+ε₃), bo podzielenie przez to wyrażenie nie zmienia znaku powyższego wyrażenia, co wtedy tak otrzymane wyrażenie, a właściwie jej prawą stronę piszemy przez:

$b_{i}{{1+\epsilon _{1}} \over {1+\epsilon _{2}}}-\lambda -{{a_{i}^{2}} \over {p_{i-1}(\lambda )}}{{(1+\epsilon _{4})(1+\epsilon _{4})(1+\epsilon _{6})} \over {(1+\epsilon _{2})(1+\epsilon _{3})}}\;$

(6.62)

Oszacowania, co do błędów poszczególnych elementów (6.57), które powstają podczas zaokrąglania, są napisane jako:

$|\delta b_{i}|<2,04\epsilon |b_{i}|\;$

(6.63)

$|\delta a_{i}|<2,06\epsilon |a_{i}|\;$

(6.64)

Ze wzorów (6.63) i (6.64) wynika ||δT||_∞<2,06ε||T||_∞. Weźmy sobie początkowy przedział $\langle \alpha ^{(0)},\beta ^{(0)}\rangle \;$, który dzielimy pokolei na m części i obierzemy przez wartość t jako środek ostatniego przedziału i mając wartość własną λ_k możemy napisać jego oszacowanie:

$|t-\lambda _{k}|<||\delta \mathbf {T} ||_{\infty }+(\beta ^{(0)}-\alpha ^{(0)})\cdot 2^{-m-1}\;$

(6.65)

Jeśli dodatkowo przyjmować będziemy β⁽⁰⁾=||T||_∞ i α⁽⁰⁾=-||T||_∞, to wzór (6.65) na podstawie oszacowania δT mamy w postaci:

$|t-\lambda _{k}|<(2,06\epsilon +2^{-m})||\mathbf {T} ||_{\infty }\;$

(6.66)

Jeśli wyznaczymy wartości własne opisane metodą powyżej, to możemy przystąpić do liczenia poszczególnych składowych wektora własnego macierzy (6.57) według sposobu:

$x_{1}=1\;$

(6.67)

$x_{2}={{(\lambda -b_{1})} \over {a_{2}}}\;$

(6.68)

$x_{i+1}={{[(\lambda -b_{i})x_{i}-a_{i}x_{i-1}]} \over {a_{i+1}}}{\mbox{ dla }}i=2,3...,n-1\;$

(6.69)

Uzyskiwanie rozkładu macierzy trójdiagonalnej symetrycznej metodą QR

Weźmy sobie macierz trójdiagonalną T (6.57), którą będziemy rozkładać na iloczyn dwóch macierzy Q i R, wtedy możemy napisać iloczyn macierzy trójdiagonalnej P⁽¹⁾, P⁽²⁾,...,P^(n-1) i macierzy A, który jest macierzą R, zapisywanej jako P^(n-1)P^(n-2)⋅...⋅P⁽²⁾P⁽¹⁾T=R. Jeśli podziałamy macierzą P⁽¹⁾ na macierz trójdiagonalną T, to ona ma za zadanie wyeliminować w macierzy T element t₂₁ do zera. Macierzą P⁽¹⁾ może być zarówno macierzą Householdera H⁽¹⁾ (6.49), jak i macierzą Givensa G⁽²⁾. Wektor u jest tak zdefiniowany by miał dwie początkowe składowe. Obie te macierze, tzn. macierz Householdera lub Givensa, możemy zapisać znając elementy "a" i "b" w formie:

${\begin{bmatrix}a&b&&&\\b&-a&&0&\\&&1\\&0&&\ddots \\&&&&1\end{bmatrix}}\;$

(6.70)

którym składowe element a i b w obu metodach odpowiednio różnią się znakiem. A więc metody Householdera i Givensa są na pewno równoważne w tym przypadku. Dla metody Givensa macierze R i Q są zapisywane jako: R=G⁽ⁿ⁾G^(n-1)⋅...⋅G⁽³⁾T=R i Q=G⁽²⁾⋅...⋅G⁽ⁿ⁾, która ta ostatnia macierz jest macierzą Hessenberga. Jeśli wymnożymy wszystkie czynniki przez siebie w macierzy Q i w R, to w ostateczności otrzymamy macierze w jego pełnych formach:

$\mathbf {Q} ={\begin{bmatrix}p_{1}\beta _{1}&p_{2}\beta _{1}&\cdots &\cdots &p_{n}\beta _{1}\\s_{2}&p_{2}\beta _{2}&\cdot &\cdots &p_{n}\beta _{2}\\&s_{3}&\ddots &\\&0&\ddots &\ddots \\&&&s_{n}&p_{n}\beta _{n}\\\end{bmatrix}}\;$

(6.71)

$\mathbf {R} ={\begin{bmatrix}r_{1}&p_{1}&q_{1}&0&\\&\ddots &\ddots &\ddots &\\&&\ddots &\ddots &q_{n-2}\\&0&&\ddots &p_{n-1}\\&&&&r_{n}\end{bmatrix}}\;$

(6.72)

Poszczególne elementy pełnej macierzy Q, tzn. p_i, β_i i s_i możemy zdefiniować wedle schematu w formie:

1)	$p_{1}=c_{2}\;$	$\beta _{1}=1\;$	$\alpha _{1}=s_{2}\;$	$\gamma _{1}=s_{2}\;$
2)	$p_{j}=c_{j+1}\alpha _{j-1}\;$	$\gamma _{j}=s_{j+1}\;$	$\beta _{j}=-c_{j}/\alpha _{j+1}\;$	$\alpha _{j}=s_{j+1}\alpha _{j-1}\;$ <-- -->|-	3)			$p_{n}=\alpha _{n-1}\;$	$\beta _{n}=-c_{n}p_{n}\;$

Odwróćmy teraz kolejność mnożenia w przedstawianiu macierzy T , wtedy otrzymujemy inną macierz ${\tilde {\mathbf {T} }}=\mathbf {R} \mathbf {Q} \;$, ale ona powinna być podobna ortogonalnie do macierzy T (6.57) według schematu RQ=Q^TTQ, wtedy na podstawie tego przedstawienie tej naszej macierzy:

${\tilde {\mathbf {T} }}={\begin{bmatrix}r_{1}c_{2}+p_{1}s_{2}&*&\\r_{2}s_{2}&\ddots &*&0&\\&\ddots &\ddots &*\\0&\ddots &\ddots &&*\\&&r_{n}s_{n}&-c_{n-1}r_{n-1}&c_{n}+p_{n-1}s_{n}\\\end{bmatrix}}\;$

(6.73)

Podczas rozkładu macierzy T (6.57) na iloczyn QR należy wykonać M=9n-9 mnożeń, D=3n-3 dodawań i n-1 pierwiastkowań. Dodatkowo, aby obliczyć macierz ${\tilde {\mathbf {T} }}=\mathbf {R} \mathbf {Q} \;$ wedle jego schematu (6.72) należy wykonać M=4n-1 mnożeń i D=n dodawań. Ilość pamięci jakie są zużywane do wyznaczenia macierzy R I g⁽²⁾,....,G⁽ⁿ⁾ jest 4n-3 komórek pamięci maszyny cyfrowej.

Sprowadzanie dowolnej macierzy symetrycznych do postaci trójdiagonalnej

Okazuje się, że można sprowadzić dowolną macierz A do postaci trójdiagonalnej symetrycznej, tzn. istnieje taka macierz przekształcenia P, dzięki któremu jest to możliwe.

Wyznaczanie macierzy trójdiagonalnej symetrycznej metodą Householdera

Obierzmy sobie macierze P⁽¹⁾, P⁽²⁾,...,P^(n-2). Pierwszą macierz tego ciągu macierzy przekształceń piszemy przez:

$\mathbf {P} ^{(1)}={\begin{bmatrix}1&0\\0&\mathbf {H} ^{(1)}\end{bmatrix}}\;$

(6.74)

Właściwości macierzy H⁽ⁱ⁾ zostały omówione w rozdziale Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Householdera. Zdefiniujmy macierz M⁽²⁾, który wychodzi z macierzy A, którą będziemy doprowadzać do postaci diagonalnej:

$\mathbf {M} ^{(2)}=\mathbf {P} ^{(1)}\mathbf {A} \mathbf {P} ^{(1)}={\begin{bmatrix}b_{1}&{\begin{matrix}a_{1}&0&\cdots &0\end{matrix}}\\{\begin{matrix}b_{1}\\a_{1}\\0\\\vdots \\0\end{matrix}}&\mathbf {A} _{(n-1)\times (n-1)}^{(2)}\\\end{bmatrix}}\;$

(6.75)

Obierzmy sobie drugą macierz przekształcenia zapisanej przy pomocy macierzy H⁽²⁾ obraną w formie:

$\mathbf {P} ^{(2)}={\begin{bmatrix}1&0\\0&1\\0&\mathbf {H} ^{(2)}\end{bmatrix}}\;$

(6.76)

wtedy możemy obliczyć macierz M⁽³⁾ otrzymaną z macierzy M⁽²⁾ poprzez macierz przekształcenia P⁽²⁾ otrzymywaną na w sposób:

$\mathbf {M} ^{(3)}=\mathbf {P} ^{(2)}\mathbf {M} ^{(2)}\mathbf {P} ^{(2)}={\begin{bmatrix}{\begin{matrix}b_{1}&a_{1}&|\\a_{1}&b_{2}&|\\0&a_{2}&|\\\vdots &0&|\\\vdots &\vdots &|\\0&0&|\\\end{matrix}}&{\begin{matrix}{\begin{matrix}0&\cdots &\cdots &0\\a_{2}&0&\cdots &0\end{matrix}}\\\\\mathbf {A} _{(n-2)\times (n-2)}^{(3)}\\\\\\\end{matrix}}\end{bmatrix}}\;$

(6.77)

Gdy wykonamy n-2 przekształceń macierzy A otrzymujemy w końcu M^(n-1)=T, czyli macierz trójdiagonalną symetryczną. Macierze H⁽ⁱ⁾ są macierzami Householdera zapisnych w postaci macierzy (6.49). Aby wyznaczyć macierz H⁽ⁱ⁾AH⁽ⁱ⁾ należy policzyć wektor $\mathbf {w} ^{(i)}={{1} \over {\tau ^{(i)}}}\mathbf {A} ^{(i)}\mathbf {u} ^{(i)}\;$, a także $\mathbf {v} ^{(i)}=\mathbf {w} ^{(i)}-{{1} \over {2\tau ^{(i)}}}\langle \mathbf {u} ^{(i)},\mathbf {w} \rangle \cdot \mathbf {u} ^{(i)}\;$. Wykorzystując definicję macierzy Householdera (6.49) możemy napisać:

$\mathbf {H} ^{(i)}\mathbf {A} ^{(i)}\mathbf {H} ^{(i)}=\left(\mathbf {I} -{{1} \over {\tau ^{(i)}}}\mathbf {u} ^{(i)}\rangle \langle \mathbf {u} ^{(i)}\right)\mathbf {A} ^{(i)}\left(\mathbf {I} -{{1} \over {\tau ^{(i)}}}\mathbf {u} ^{(i)}\rangle \langle \mathbf {u} ^{(i)}\right)=\;$
$=\mathbf {A} ^{(i)}-{{1} \over {\tau ^{(i)}}}\left(\mathbf {A} ^{(i)}\mathbf {u} ^{(i)}\right){\mathbf {u} ^{(i)}}^{T}-{{1} \over {\tau ^{(i)}}}({\mathbf {u} ^{(i)}})\left({\mathbf {u} ^{(i)}}^{T}\mathbf {A} ^{(i)}\right)+{{1} \over {(\tau ^{(i)})^{2}}}\mathbf {u} ^{(i)}({\mathbf {u} ^{i}}^{T}\mathbf {A} ^{(i)}\mathbf {u} ){\mathbf {u} ^{(i)}}^{T}=\;$
$=\mathbf {A} ^{(i)}-\mathbf {w} ^{(i)}{\mathbf {u} ^{(i)}}^{T}-\mathbf {u} ^{(i)}{\mathbf {w} ^{(i)}}^{T}+{{1} \over {\tau ^{(i)}}}\mathbf {u} ^{(i)}\left({\mathbf {u} ^{(i)}}^{T}\mathbf {w} ^{(i)}\right){\mathbf {u} ^{(i)}}^{T}=\;$
$=\mathbf {A} ^{(i)}-{\mathbf {u} ^{(i)}}\left(\mathbf {w} ^{(i)}-{{1} \over {2\tau ^{(i)}}}\langle \mathbf {u} ^{(i)},\mathbf {w} ^{(i)}\rangle \mathbf {u} ^{(i)}\right)^{T}-\left(\mathbf {w} ^{(i)}-{{1} \over {2\tau ^{(i)}}}\mathbf {u} ^{(i)}\langle \mathbf {u} ^{(i)},\mathbf {w} ^{(i)}\rangle \right){\mathbf {u} ^{(i)}}^{T}=\;$
$=\mathbf {A} ^{(i)}-\mathbf {u} ^{(i)}\rangle \langle \mathbf {v} ^{(i)}-\mathbf {v} ^{(i)}\rangle \langle \mathbf {u} ^{(i)}$

(6.78)

Obliczenia przeprowadzone w punkcie (6.78) bardzo są potrzebne do obliczania pokolei macierzy M⁽ⁱ⁾. Aby wyznaczyć trójdiagonalną symetryczną macierz należy dokonać $M={{2} \over {3}}n^{3}+O_{1}(n^{2})\;$ mnożeń i $D={{2} \over {3}}n^{3}n^{3}+O_{2}(n^{2})\;$ dodawań oraz n-2 pierwiastkowań. Błąd w tej metodzie wyznaczania macierzy dokładnej T w jego przybliżonej wersji ${\tilde {\mathbf {T} }}\;$ jest:

$||\mathbf {T} -{\tilde {\mathbf {T} }}||_{E}\leq 2k_{1}n^{2}\epsilon \left(1+k_{1}\epsilon n\right)^{2n}||\mathbf {A} ||_{E}\;$

(6.79)

Stała k₁ we wzorze (6.79) zależy od sposobu zaokrąglania na typowej maszynie cyfrowej i dla typowej maszyny cyfrowej jej wartość jest k₁=20.

Wyznaczanie macierzy trójdiagonalnej symetrycznej metodą Givensa

Przy wyznaczaniu macierzy P⁽ⁱ⁾, która jest macierzą przekształcenia, musimy skorzystać z macierzy Givensa opracowanej w punkcie Metoda rozkładu dowolnej macierzy A na iloczyn macierzy Q i R metodą Givensa, i na jej podstawie mamy macierze przekształcenia, co w rezultacie mamy ciąg M⁽²⁾,M⁽³⁾,...,M^(n-1)=T. W tej metodzie dokonaliśmy $M={{4} \over {3}}n^{3}+O_{1}(n^{2})\;$ mnożeń i $D={{2} \over {3}}n^{3}+O_{2}(n^{2})\;$ dodawań, a także ${{1} \over {2}}n(n-1)\;$ pierwiastkowań. Błąd w tej metodzie wyznaczania macierzy dokładnej T w jego przybliżonej wersji ${\tilde {\mathbf {T} }}\;$ jest:

$||\mathbf {T} -{\tilde {\mathbf {T} }}||_{E}\leq 12n^{{3} \over {2}}\epsilon \left(1+6\epsilon \right)^{4n-7}||\mathbf {A} ||_{E}\;$

(6.80)