Gramova–Schmidtova ortogonalizace

Šablona:Upravit 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 ${ℝ}^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,\dots ,a_n\in {ℝ}^m$ jsou, opět pro jednoduchost, lineárně nezávislé, tedy $n\le m$. Úkolem je nalézt ortonormální bázi $q_1,\dots ,q_n$ $n$-rozměrného podprostoru ${ℝ}^m$, který vektory $a_i$ generují; má tedy platit

$\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{a_1,\dots ,a_n\}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{q_1,\dots ,q_n\},⟨q_k,q_j⟩=q_j^T{q}_k={\delta }_{k,j}$

kde $\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}$ 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 $\Vert q_1\Vert =1$, musí platit

$a_1=q_1{r}_{1,1},\text{kde}{r}_{1,1}=\Vert a_1{\Vert }_2,$

a dostáváme vztah pro výpočet prvního vektoru ortonormální báze $q_1=a_1/\Vert a_1{\Vert }_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,\Vert q_2{\Vert }_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=\Vert q_1{\Vert }_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},\text{kde}{r}_{2,2}=\Vert a_2-q_1{r}_{1,2}{\Vert }_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},\text{kde matice}(I_m-q_1{q}_1^T)$

není nic jiného, než ortogonální projektor do ortogonálního doplňku $\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{q_1{\}}^{\perp }$ lineárního obalu vektoru $q_1$ v ${ℝ}^m$.

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

Algoritmicky zapsáno:

00: vstup $a_1,\dots ,a_n$

01: $r_{1,1}:=\Vert a_1{\Vert }_2$

02: $q_1:=a_1/r_{1,1}$

03: for $k:=2,\dots ,n$

04:     $p:=a_k$

05:     for $j:=1,\dots ,k-1$

06:        $r_{j,k}:=q_j^{T}p=⟨p,q_j⟩$

07:     end

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

09:        $p:=p-q_j{r}_{j,k}$

10:     end

11:     $r_{k,k}:=\Vert p{\Vert }_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,\dots ,q_n$.

Označíme-li $Q_k=[q_1,\dots ,q_k]\in {ℝ}^{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)\dots (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),\dots ,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

$\Vert {\hat{Q}}_n^T{\hat{Q}}_n-I_n{\Vert }_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& \dots & 1\\ \rho & & 0\\ & \ddots & \\ 0& & \rho \end{array}\right]\in {ℝ}^{(n+1)\times n},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í ${ϵ}_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){ϵ}_M$
MGS	$2.2\times 10^{-9}$	${\kappa }_2(A){ϵ}_M$
ICGS	$2.4\times 10^{-16}$	${ϵ}_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,\dots ,a_n],Q=[q_1,\dots ,q_n]\in {ℝ}^{m\times n},R=\left[\begin{array}{cccc}{r}_{1,1}& r_{1,2}& \dots & r_{1,n}\\ 0& r_{2,2}& \dots & r_{2,n}\\ ⋮& \ddots & \ddots & ⋮\\ 0& \dots & 0& r_{n,n}\end{array}\right]\in {ℝ}^{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ů $𝒫(a,b)$ se skalárním součinem

$⟨p(x),q(x){⟩}_w=\underset{a}{\overset{b}{\int }}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,\dots ,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. Šablona:ISBN. (Kapitola 3, Ortogonální transformace a QR rozklady, str. 53-88.)