Gramova-Schmidtova ortogonalizace

Z Wikipedie, otevřené encyklopedie
Skočit na: Navigace, Hledání

Gramův-Schmidtův proces neboli Gramova-Schmidtova ortogonalizace (nesprávně[1] 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[editovat | editovat zdroj]

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. Ú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}

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). 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 zbytu 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í[editovat | editovat zdroj]

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. 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[editovat | editovat zdroj]

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. 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[2]

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.22\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

Algoritmus Ztráta ortogonality (Lauchliho matice) Ztráta ortogonality (obecně)
CGS 2.2\times10^{-2} \kappa_2^2(A)\epsilon_M
MGS 2.2\times10^{-9} \kappa_2(A)\epsilon_M
ICGS 2.4\times10^{-16} \epsilon_M

Vztah Gramova-Schmidtova algoritmu a QR rozkladu[editovat | editovat zdroj]

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[editovat | editovat zdroj]

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[editovat | editovat zdroj]

  1. Internetová jazyková příručka [online]. Ústav pro jazyk český Akademie věd ČR, [cit. 2011-10-31]. Dostupné online.  
  2. Lauchli matrix [online]. Matrix Market, [cit. 2011-11-03]. Dostupné online.  

Literatura[editovat | editovat zdroj]

  • Gene Howard Golub, Charles F. Van Loan: Matrix Computations, John 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. ISBN 978-80-7378-201-6. (Kapitola 3, Ortogonální transformace a QR rozklady, str. 53-88.)