Gramova–Schmidtova ortogonalizace

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

${\displaystyle \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 ${\displaystyle \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 ${\displaystyle \|q_{1}\|=1}$, musí platit

${\displaystyle 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 ${\displaystyle q_{1}=a_{1}/\|a_{1}\|_{2}}$. Protože ${\displaystyle a_{2}}$ je lineárně nezávislý na ${\displaystyle a_{1}}$ a tedy i na ${\displaystyle q_{1}}$, můžeme ho vyjádřit jako

${\displaystyle a_{2}=q_{1}r_{1,2}+q_{2}r_{2,2},}$

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

${\displaystyle 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 ${\displaystyle q_{1}^{T}q_{1}=\|q_{1}\|_{2}^{2}=1}$), dostaneme vztah pro výpočet ${\displaystyle r_{1,2}}$ (ortogonalizační koeficient; tj. velikost projekce ${\displaystyle a_{2}}$ do směru ${\displaystyle q_{1}}$). Protože známe ${\displaystyle a_{2},\,q_{1},\,r_{1,2}}$ dostáváme

${\displaystyle 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 ${\displaystyle a_{2}}$ po ortogonalizaci proti ${\displaystyle q_{1}}$. Všimněme si, že po dosazení za ${\displaystyle r_{1,2}}$ dostáváme

${\displaystyle 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 ${\displaystyle \mathrm {span} \{q_{1}\}^{\perp }}$ lineárního obalu vektoru ${\displaystyle q_{1}}$ v ${\displaystyle \mathbb {R} ^{m}}$.

Tento postup lze zřejmě opakovat do vyčerpání všech vektorů ${\displaystyle a_{k}}$.

Algoritmicky zapsáno:

00: vstup ${\displaystyle a_{1},\ldots ,a_{n}}$

01: ${\displaystyle r_{1,1}:=\|a_{1}\|_{2}}$

02: ${\displaystyle q_{1}:=a_{1}/r_{1,1}}$

03: for ${\displaystyle k:=2,\ldots ,n}$

04:     ${\displaystyle p:=a_{k}}$

05:     for ${\displaystyle j:=1,\ldots ,k-1}$

06:        ${\displaystyle r_{j,k}:=q_{j}^{T}p=\langle p,q_{j}\rangle }$

07:     end

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

09:        ${\displaystyle p:=p-q_{j}r_{j,k}}$

10:     end

11:     ${\displaystyle r_{k,k}:=\|p\|_{2}}$

12:     ${\displaystyle 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ů.

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ů ${\displaystyle r_{j,k}}$) lze provádět nezávisle pro jednotlivá ${\displaystyle j}$; tedy, jednotlivá ${\displaystyle 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 ${\displaystyle q_{1},\ldots ,q_{n}}$.

Označíme-li ${\displaystyle 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 ${\displaystyle a_{k+1}}$ psát pomocí projektorů dvěma matematicky ekvivalentními způsoby

${\displaystyle 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 ${\displaystyle j}$-tém kroku částečně eliminovat chyby vzniklé zaokrouhlováním v předchozích krocích ${\displaystyle (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 ${\displaystyle q_{i}}$ je pak dokonce menší než u MGS).

Ztráta ortogonality

Nechť ${\displaystyle {\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. ${\displaystyle {\hat {Q}}_{n}=Q_{n}+E_{n}\approx Q_{n}}$ a ${\displaystyle {\hat {Q}}_{n}^{T}{\hat {Q}}_{n}\approx I_{n}}$. Veličina

${\displaystyle \|{\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]

${\displaystyle 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 ${\displaystyle \kappa _{2}(A)}$ je podmíněnost matice ${\displaystyle A}$. Uvažujeme-li standardní aritmetiku se strojovou přesností ${\displaystyle \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 ${\displaystyle A}$

Algoritmus	Ztráta ortogonality (Lauchliho matice)	Ztráta ortogonality (obecně)
CGS	${\displaystyle 2.2\times 10^{-2}}$	${\displaystyle \kappa _{2}^{2}(A)\epsilon _{M}}$
MGS	${\displaystyle 2.2\times 10^{-9}}$	${\displaystyle \kappa _{2}(A)\epsilon _{M}}$
ICGS	${\displaystyle 2.4\times 10^{-16}}$	${\displaystyle \epsilon _{M}}$

Srovnáním sloupcových vektorů ${\displaystyle a_{k},\;q_{j}}$ a koeficientů ${\displaystyle r_{j,k}}$ do matic,

${\displaystyle 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 ${\displaystyle Q^{T}Q=I_{n}}$ a ${\displaystyle R}$ je čtvercová regulární matice dostáváme

${\displaystyle A=QR}$

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

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

${\displaystyle \langle p(x),q(x)\rangle _{w}=\int _{a}^{b}p(x)q(x)w(x)dx}$,

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

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

• 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.)