Gramova–Schmidtova ortogonalizace

Technology
12 hours ago
8
4
2
Avatar
Author
Albert Flores

Gramův-Schmidtův proces neboli Gramova-Schmidtova ortogonalizace (nesprávně Gram-Schmidtova ortogonalizace) je metoda, která v daném unitárním prostoru (neboli vektorovém prostoru se skalárním součinem) umožňuje pro zadanou konečnou množinu vektorů nalézt ortonormální bázi podprostoru jimi generovaného.

Algoritmus

Uvažujme pro jednoduchost \mathbb{R}^m reálný lineární vektorový prostor sloupcových vektorů o m složkách (se standardním skalárním součinem). Nechť a_1,\ldots,a_n\in\mathbb{R}^m jsou, opět pro jednoduchost, lineárně nezávislé, tedy n\leq m. +more Úkolem je nalézt ortonormální bázi q_1,\ldots,q_n n-rozměrného podprostoru \mathbb{R}^m, který vektory a_i generují; má tedy platit.

:\mathrm{span}\{a_1,\ldots,a_n\}=\mathrm{span}\{q_1,\ldots,q_n\},\qquad \langle q_k,q_j\rangle=q_j^Tq_k = \delta_{k,j}

kde \mathrm{span} značí lineární obal množiny v závorce.

Algoritmus danou sadu vektorů prochází postupně přičemž v každém kroku vygeneruje nový vektor hledané báze. Omezíme-li se pouze na první vektor, a protože požadujeme aby \|q_1\|=1, musí platit

:a_1 = q_1r_{1,1},\qquad \text{kde}\quad r_{1,1}=\|a_1\|_2,

a dostáváme vztah pro výpočet prvního vektoru ortonormální báze q_1=a_1/\|a_1\|_2. Protože a_2 je lineárně nezávislý na a_1 a tedy i na q_1, můžeme ho vyjádřit jako

:a_2 = q_1r_{1,2} + q_2r_{2,2},

kde q_2 je nějaký nový vektor takový, že q_1^Tq_2=0,\;\|q_2\|_2=1. Po pronásobení předchozího vztahu q_1^T zleva,

:q_1^Ta_2 = q_1^Tq_1r_{1,2} + q_1^Tq_2r_{2,2}=r_{1,2}

(připomeňme, že q_1^Tq_1=\|q_1\|_2^2=1), dostaneme vztah pro výpočet r_{1,2} (ortogonalizační koeficient; tj. velikost projekce a_2 do směru q_1). +more Protože známe a_2,\,q_1,\,r_{1,2} dostáváme.

:a_2-q_1r_{1,2} = q_2r_{2,2}, \qquad\text{kde}\quad r_{2,2}=\|a_2-q_1r_{1,2}\|_2

je norma zbytku vektoru a_2 po ortogonalizaci proti q_1. Všimněme si, že po dosazení za r_{1,2} dostáváme

:a_2-q_1q_1^Ta_2 = (I_m-q_1q_1^T)a_2 = q_2r_{2,2},\qquad\text{kde matice}\quad (I_m-q_1q_1^T)

není nic jiného, než ortogonální projektor do ortogonálního doplňku \mathrm{span}\{q_1\}^\perp lineárního obalu vektoru q_1 v \mathbb{R}^m.

Tento postup lze zřejmě opakovat do vyčerpání všech vektorů a_k.

Algoritmicky zapsáno:

:00: vstup a_1,\ldots,a_n :01: r_{1,1} := \|a_1\|_2 :02: q_1 := a_1/r_{1,1} :03: for k := 2,\ldots,n :04: p := a_k :05: for j := 1,\ldots,k-1 :06: r_{j,k} := q_j^Tp = \langle p,q_j\rangle :07: end :08: for j := 1,\ldots,k-1 :09: p := p - q_jr_{j,k} :10: end :11: r_{k,k} := \|p\|_2 :12: q_k := p/r_{k,k} :13: end

Tato varianta algoritmu se nazývá klasický Gramův-Schmidtův algoritmus (CGS) a je novější než varianta původní, dnes zvaná modifikovaný Gramův-Schmidtův algoritmus (MGS). MGS získáme z výše popsaného CGS prostým vypuštěním řádků 07 a 08, tedy, spojením obou vnitřních cyklů.

Varianty algoritmu a jejich chování

Oba algoritmy CGS i MGS jsou matematicky ekvivalentní, jejich reálné implementace mají výrazně odlišné chování.

CGS algoritmus je výrazně paralelní. Výpočet první vnitřní smyčky (tj. +more výpočet koeficientů r_{j,k}) lze provádět nezávisle pro jednotlivá j; tedy, jednotlivá r_{j,k} mohu počítat na různých procesorech, jejich výpočet se neovlivňuje, nezávisí na sobe a může probíhat paralelně. Zatímco MGS je z tohoto pohledu sekvenční.

Na druhou stranu výpočet pomocí MGS je numericky výrazně stabilnější než výpočet pomocí CGS, kde může, vlivem zaokrouhlovacích chyb, dojít k úplné ztrátě ortogonality mezi vektory q_1,\ldots,q_n.

Označíme-li Q_{k}=[q_1,\ldots,q_{k}]\in\mathbb{R}^{m\times k},\;Q_{k}^TQ_{k}=I_{k}, lze vztah pro ortogonalizaci vektoru a_{k+1} psát pomocí projektorů dvěma matematicky ekvivalentními způsoby

:p := (I_m-Q_{k}Q_{k}^T)a_{k+1}=(I_m-q_kq_k^T)\ldots(I_m-q_2q_2^T)(I_m-q_1q_1^T)a_{k+1}.

První projekce odpovídá výpočtu pomocí CGS, druhá postupná výpočtu pomocí MGS. Je zřejmé že CGS ortogonalizace (projekce) se počítá paralelně do všech směrů najednou, kdežto sekvenční ortogonalizace (projekce) MGS umožňuje v j-tém kroku částečně eliminovat chyby vzniklé zaokrouhlováním v předchozích krocích (j-1),\ldots,2,1.

Řešením v praxi používaným bývá tzv. klasický Gramův-Schmidtův algoritmus s iteračním zpřesněním (ICGS), který obsahuje dvě vnitřní smyčky jako CGS (je tedy paralelizovatelný), ale obě smyčky se provedou dvakrát (čímž se výrazně zlepší numerické vlastnosti, ztráta ortogonality mezi vektory q_i je pak dokonce menší než u MGS).

Ztráta ortogonality

Nechť \hat{Q}_n je matice vektorů spočtených pomocí některé varianty Gramova-Schmidtova algoritmu v počítači se standardní konečnou aritmetikou s plovoucí řádovou čárkou, tj. \hat{Q}_n=Q_n+E_n\approx Q_n a \hat{Q}_n^T\hat{Q}_n\approx I_n. +more Veličina.

:\|\hat{Q}_n^T\hat{Q}_n-I_n\|_2

se nazývá ztráta ortogonality a je jednou z klíčových veličin sloužících k posouzení kvality spočtené ortonormální báze.

Uvažujme tzv. Lauchliho matici

:A=\left[\begin{array}{ccc}1&\ldots&1\\\rho&&0\\&\ddots&\\0&&\rho\end{array}\right]\in\mathbb{R}^{(n+1)\times n}, \qquad n=20,\; \rho=10^{-7},\; \kappa_2(A)\approx 4.47\times10^7,

kde \kappa_2(A) je podmíněnost matice A. Uvažujeme-li standardní aritmetiku se strojovou přesností \epsilon_M\approx2. +more22\times10^{-16} (double), pak ztráta ortogonality odpovídající jednotlivým výše zmíněným algoritmům aplikovaným na danou Lauchliho matici, je ve druhém sloupci následující tabulky. Ve třetím sloupci je obecný vztah platný pro libovolnou matici A.

CGS2. 2\times10^{-2}\kappa_2^2(A)\epsilon_M
MGS2. +more2\times10^{-9}\kappa_2(A)\epsilon_M
ICGS2. 4\times10^{-16}\epsilon_M
.

Vztah Gramova-Schmidtova algoritmu a QR rozkladu

Srovnáním sloupcových vektorů a_k,\;q_j a koeficientů r_{j,k} do matic,

:A=[a_1,\ldots,a_n],\quad Q=[q_1,\ldots,q_n]\in\mathbb{R}^{m\times n},\quad R=\left[\begin{array}{cccc}r_{1,1}&r_{1,2}&\ldots&r_{1,n} \\ 0&r_{2,2}&\ldots&r_{2,n}\\\vdots&\ddots&\ddots&\vdots\\0&\ldots&0&r_{n,n}\end{array}\right]\in\mathbb{R}^{n\times n},

kde Q^TQ=I_n a R je čtvercová regulární matice dostáváme

:A=QR

tedy QR rozklad matice A (pro ověření stačí porovnat k-tý sloupec rovnosti, tedy a_k s k-tým sloupcem součinu QR).

Ortogonální polynomy

Gramův-Schmidtův algoritmus lze aplikovat na prvky libovolného prostoru se skalárním součinem. Uvažujeme například prostor polynomů \mathcal{P}(a,b) se skalárním součinem

:\langle p(x),q(x)\rangle_w = \int_a^b p(x)q(x)w(x)dx,

kde w(x) je nějaká váhová funkce. Aplikací Gramova-Schmidtova algoritmu na sadu vektorů (polynomů) 1,x,x^2,\ldots,x^{n-1} (v tomto pořadí) dostáváme, pro vhodně volené a,\,b a váhovou funkci w(x) libovolnou sadu ortogonálních polynomů (normalizovanou vzhledem k normě indukované daným skalárním součinem).

Pro a=-1,\,b=1,\,w(x)=1 Gramův-Schmidtův algoritmus generuje Legenderovy polynomy, pro a=-1,\,b=1,\,w(x)=(1-x^2)^{-1/2} dostaneme Čebyševovy polynomy prvního druhu, atd.

Reference

Literatura

Gene Howard Golub, Charles F. +more Van Loan: Matrix Computations, Johns Hopkins University Press, 1996 (3rd Ed. ). (Zejména kapitoly 5. 2. 7 CGS, 5. 2. 8 MGS a 5. 2. 9 Work and Accuracy. ) * J. Duintjer Tebbens, I. Hnětynková, M. Plešinger, Z. Strakoš, P. Tichý: Analýza metod pro maticové výpočty, základní metody. Matfyzpress 2012. . (Kapitola 3, Ortogonální transformace a QR rozklady, str. 53-88. ).

Kategorie:Lineární algebra

5 min read
Share this post:
Like it 8

Leave a Comment

Please, enter your name.
Please, provide a valid email address.
Please, enter your comment.
Enjoy this post? Join Cesko.wiki
Don’t forget to share it
Top