Gramova–Schmidtova ortogonalizace

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$. Ú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}^{T}q_{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_{1}r_{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_{1}r_{1,2}+q_{2}r_{2,2},$

kde $q_{2}$ je nějaký nový vektor takový, že $q_{1}^{T}q_{2}=0,\;\|q_{2}\|_{2}=1$. Po pronásobení předchozího vztahu $q_{1}^{T}$ zleva,

$q_{1}^{T}a_{2}=q_{1}^{T}q_{1}r_{1,2}+q_{1}^{T}q_{2}r_{2,2}=r_{1,2}$

(připomeňme, že $q_{1}^{T}q_{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_{1}r_{1,2}=q_{2}r_{2,2},\qquad {\text{kde}}\quad r_{2,2}=\|a_{2}-q_{1}r_{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_{1}q_{1}^{T}a_{2}=(I_{m}-q_{1}q_{1}^{T})a_{2}=q_{2}r_{2,2},\qquad {\text{kde matice}}\quad (I_{m}-q_{1}q_{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}^{T}p=\langle p,q_{j}\rangle$

07:     end

08:     for $j:=1,\ldots ,k-1$

09:        $p:=p-q_{j}r_{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. 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}^{T}Q_{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_{k}q_{k}^{T})\ldots (I_{m}-q_{2}q_{2}^{T})(I_{m}-q_{1}q_{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}$. 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\times 10^{7},$

kde $\kappa _{2}(A)$ je podmíněnost matice $A$. Uvažujeme-li standardní aritmetiku se strojovou přesností $\epsilon _{M}\approx 2.22\times 10^{-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\times 10^{-2}$	$\kappa _{2}^{2}(A)\epsilon _{M}$
MGS	$2.2\times 10^{-9}$	$\kappa _{2}(A)\epsilon _{M}$
ICGS	$2.4\times 10^{-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^{T}Q=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.

Literatura

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