QR algoritmus

Technology

12 hours ago

Author

QR algoritmus je numerická metoda pro výpočet vlastních čísel obecné čtvercové \boldsymbol A založená na principu QR rozkladu. Výhodou algoritmu je numerická stabilita.

Využití rozkladu

Nechť \boldsymbol A = \boldsymbol{QR}, kde \boldsymbol Q je ortogonální matice a \boldsymbol R je horní trojúhelníková matice. Pak matice \boldsymbol B = \boldsymbol{RQ} je podobná matici \boldsymbol A, protože: :\boldsymbol B = \boldsymbol{RQ} \implies \boldsymbol{BQ}^{-1} =\boldsymbol R \implies\boldsymbol A = \boldsymbol{QBQ}^{-1} = \boldsymbol{QBQ}^\mathrm{T} Stejným způsobem je možné vyjádřit matici \boldsymbol B = \boldsymbol Q^\mathrm{T}\boldsymbol{AQ}. +more Z podobnosti plyne, že obě matice mají shodná vlastní čísla.

Algoritmus

Finitní transformace na horní Hessenberův tvar (preprocessing)

Iteračnímu QR algoritmu zpravidla předchází transformace matice typicky do horního Hessenbergova tvaru :\boldsymbol A = \boldsymbol{QHQ}^\mathrm{T}, \qquad \text{kde} \qquad \boldsymbol H = \begin{pmatrix} h_{11}&h_{12}&h_{13}&\cdots&h_{1n}\\ h_{21}&h_{22}&h_{23}&\cdots&h_{2n}\\ &h_{32}&h_{33}&\cdots&h_{3n}\\ &&\ddots&\ddots&\vdots\\ &&&h_{n,n-1}&h_{n,n}\\ \end{pmatrix}, tedy do „skoro trojúhelníkového“ tvaru, až na obecně nenulovou první poddiagonálu. Tuto transformaci lze provést v konečném počtu kroků (tj. +more pomocí konečného počtu elementárních aritmetických operací sčítání, odečítání, násobení, dělení a výpočtu druhé odmocniny), např. pomocí tzv. Arnoldiho algoritmu.

Vlastní iterační algoritmus

Samotný iterační QR algoritmus konverguje limitně (pokud konverguje) # k := 0, \boldsymbol A_0:= \boldsymbol A je zadaná matice, případně matice \boldsymbol H v horním Hessenbergově tvaru, # vypočteme QR rozklad \boldsymbol A_{k} = \boldsymbol Q_k\boldsymbol R_k, # vypočteme novou matici \boldsymbol A_{k+1} := \boldsymbol R_k\boldsymbol Q_k (z předchozích úvah je zřejmé, že \boldsymbol A_{k+1} = \boldsymbol Q_k^\mathrm{T}\boldsymbol A_k\boldsymbol Q_k), # pokud je \boldsymbol A_{k+1} trojúhelníková matice, má vlastní čísla na diagonále. Jinak položíme k := k+1 a pokračujeme druhým krokem algoritmu.

V právě popsaném algoritmu je matice \boldsymbol Q_k přímo maticí z QR rozkladu, v jistém smyslu tedy eliminuje všechny poddiagonální prvky matice \boldsymbol A_k. Jednotlivé iterace je ale možné dělat „jemněji“, resp. +more postupně. Matice \boldsymbol Q_k může být volena tak, aby byl eliminovala jeden nenulový prvek matice \boldsymbol A_k pod diagonálou. (Možným postupem je volit \boldsymbol Q_k jako Givensovu rotaci \boldsymbol G(i,j,\theta), kde i a j jsou indexy řádku a sloupce eliminovaného prvku; k tomu je potřeba nalézt úhel \theta, pod kterým lze prvek eliminovat. ).

Horního Hessenbergova tvar má hned několik výhod: Tento tvar je invariantní vůči transformaci \boldsymbol A_k\;\longmapsto\;\boldsymbol A_{k+1} popsané v bodech 2 a 3, tedy je-li \boldsymbol A_0 v horním Hessenbergově tvaru, jsou také všechny matice \boldsymbol A_k v horním Hessenbergově tvaru. V matici navíc potřebujeme eliminovat pouze n-1 nenulových poddiagonálních prvků, k výpočtu QR rozkladu tedy potřebujeme pouze n-1 Givensových rotací.

Pořadí eliminovaných prvků může být např. podle indexů nebo podle velikosti prvku v absolutní hodnotě (vždy ten největší) atp. +more Další výhodou horního Hessenbergova tvaru je triviální pozorování, že se s každou úspěšnou eliminací poddiagonálního prvku:.

* buď zmenší efektivní rozměr matice o jedna přičemž zároveň dojde k identifikaci jednoho vlastního čísla (když dojde k eliminaci prvního, nebo posledního poddiagonálního prvku), * nebo se úloha rozpadne na dvě menší zcela nezávislé podúlohy, což je možné využít k paralelizaci (ve všech ostatních případech).

Poznámka ke konvergenci

Zde prezentovaný QR algoritmus nemusí vždy konvergovat ke trojúhelníkové matici. V praxi se (nejen) proto používá řada sofistikovaných „triků“, které chování algoritmu vylepšují.

Pro příklad, kdy shora uvedený algoritmus nebude konvergovat, stačí uvažovat reálnou čtvercovou matici \boldsymbol A, která má alespoň jedno vlastní číslo s nenulovou imaginární složkou (pro jednoduchost budeme říkat komplexní vlastní číslo). V důsledku jsou zřejmě všechny matice \boldsymbol A_k, \boldsymbol Q_k, \boldsymbol R_k reálné (nezávisle na tom, zda matice \boldsymbol Q_k a \boldsymbol R_k (viz krok 2) pocházejí přímo z QR rozkladu, nebo zda mají původ jen v jediné Givensově rotaci; součiny reálných matic \boldsymbol A_{k+1} (viz krok 3) jsou opět pouze reálné matice). +more Trojúhelníková matice, ke které bychom rádi dokonvergovali (viz krok 4), má ale na diagonále vlastní čísla matice \boldsymbol A, je tedy nutně komplexní.

Jacobiova diagonalizace

Speciálním případem QR algoritmu je Jacobiova diagonalizace pojmenovaná po matematikovi Carlu Gustavu Jacobu Jacobiovi. Tento případ nastává, pokud je matice \boldsymbol A symetrická. +more Při volbě \boldsymbol Q_k = \boldsymbol G(i,j,\theta) dochází k eliminaci nejen prvku \boldsymbol A_{ij}, ale také prvku \boldsymbol A_{ji}. Výsledkem je pak diagonální matice podobná \boldsymbol A.