Metode QR#

Algoritmul QR (a nu se confunda cu factorizarea QR) este utilizat pentru a calcula simultan valorile proprii ale unei matrice pătratice \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\).

Vom încerca mai întâi să explicăm proveniența algoritmului, iar mai apoi vom descoperi împreună algoritmul, alături de îmbunătățirile sale.

Evoluția din MPD#

Rugăm cititorii să își amintească cum funcționa metoda puterii directe înainte de a continua.

Fie \(A\in\mathbb{R}^{n\times n}\) o matrice pătratică, \(n\in\mathbb{N}^*\). Vom nota cu \(\bm{y_i'}^{(k)}\) vectorul propriu (la a \(k\)-a iterație) ce este asociat valorii proprii \(\lambda_i\). Suplimentar, vom considera momentan cunoscut vectorul propriu \(\bm{q_1}\) asociat cu \(\lambda_1\).

Vrem acum să studiem vectorul \(\bm{y_2'}^{(k)}\).

Întrucât am stabilit deja că \(\bm{y_1}^{(k)}\) va fi un vector cu norma \(1\), dorim să îl scriem pe \(\bm{y_2'}^{(k)}\) în aceeași manieră. Pentru a simplifica și mai mult calculele, vom forța ca \(\bm{y_1'}^{(k)}\perp\bm{y_2'}^{(k)}\).

Asemănător cu Gram-Schmidt, proiecția lui \(\bm{y_2'}^{(k)}\) pe \(\bm{q_1}\) poate fi scrisă sub forma \(\left\langle\bm{q_1},\bm{\bm{y_2'}^{(k)}}\right\rangle\bm{q_1}\), putem asigura \(\bm{y_1'}^{(k)}\perp\bm{y_2'}^{(k)}\) prin eliminarea acestei proiecții din \(\bm{y_2'}^{(k)}\), adică urmând trecerea:

\[ \boxed{\bm{y_2'}^{(k)}\rightarrow\bm{y_2'}^{(k)}-\left\langle\bm{q_1},\bm{\bm{y_2'}^{(k)}}\right\rangle\bm{q_1}} \]

Așadar, în urma analizei anterioare, am găsit o formulă pentru \(\bm{y_2'}^{(k)}\), însă problema este că am presupun deja cunoscută valoarea \(\bm{q_1}\), ceea ce face algoritmul pur teoretic.

Metoda practică#

Pentru a elimina dependența de \(\bm{q_1}\), putem rula în paralel MPD pentru a afla \(\bm{y_1'}^{(k)}\) (așa cum fusese el conceput inițial). În acest caz, întrucât \(\bm{y_1'}^{(k)}\approx\bm{q_1}\), poate fi înlocuită cu totul dependența de \(\bm{q_1}\) urmând acești pași simpli:

\(\bm{y_1'}^{(0)}\) și \(\bm{y_2'}^{(0)}\) primesc o aproximare de început;
Se ortogonalizează \(\bm{y_2'}^{(0)}\) în raport cu \(\bm{y_1'}^{(0)}\), adică:
\[ \bm{y_2'}^{(0)}\rightarrow\bm{y_2'}^{(0)}-\left\langle\bm{q_1},\bm{\bm{y_2'}^{(0)}}\right\rangle\bm{q_1} \]
Se normalizează \(\bm{y_1'}^{(0)}\) și \(\bm{y_2'}^{(0)}\);
Pentru fiecare pas \(k\geq 1\), se va obține o aproximare tot mai bună astfel:
- Se calculează \(\bm{y_1'}^{(k)}\):
  
  Se va folosi un vector auxiliar, \(\bm{t_1}=A\bm{y_1'}^{(k-1)}\), ce va fi apoi normalizat și reținut, adică \(\bm{y_1'}^{(k)}=\frac{\bm{t_1}}{||\bm{t_1}||}\);
- Se calculează \(\bm{y_2'}^{(k)}\):
  
  Se va folosi un vector auxiliar, \(\bm{t_{21}}=A\bm{y_2}^{(k-1)}\), a cărui componentă de pe \(\bm{y_1'}^{(k)}\) va fi eliminată, adică \(\bm{t_{22}}=\bm{t_{21}}-\left\langle\bm{q_1},\bm{\bm{y_2'}^{(k)}}\right\rangle\bm{q_1}\).
  
  Vectorul \(\bm{y_2'}^{(k)}\) este \(\bm{t_{22}}\) normalizat, adică \(\bm{y_2'}^{(k)}=\frac{\bm{t_{22}}}{||\bm{t_{22}}||}\);
- Se calculează noua aproximare a lui \(\lambda_1\) și a lui \(\lambda_2\) folosind coeficientul Rayleigh, adică:
  \[ \begin{cases} \lambda_1^{(k)}=\left(\bm{y_1'}^{(k)}\right)^TA\left(\bm{y_1'}^{(k)}\right)\\ \lambda_2^{(k)}=\left(\bm{y_2'}^{(k)}\right)^TA\left(\bm{y_2'}^{(k)}\right) \end{cases} \]

Trecerea la factorizare QR#

În baza algoritmului prezentat anterior, putem face câteva observații.

În primul rând, pașii de indice 2 și 3 sunt echivalenți calculării matricei ortogonale \(\tilde{Y}^{(0)}=\begin{bmatrix}\bm{y_1'}^{(0)} & \bm{y_2'}^{(0)}\end{bmatrix}\in\mathbb{R}^{n\times 2}\).

Ne amintim că există o descompunere specifică, descompunerea QR, care ne permite să ajungem la un produs între o matrice ortogonală și o matrice superior triunghiulară. Așadar, obținem:

\[ \left(\tilde{Y}^{(0)},R\right)=QR\left(\begin{bmatrix}\bm{y_1'}^{(0)} & \bm{y_2'}^{(0)}\end{bmatrix}\right) = QR\left(Y^{(0)}\right) \]

Suplimentar, ne îndreptăm atenția către pasul 4, și anume către calcularea valorilor \(\bm{y_1'}^{(k)}\) și \(\bm{y_2'}^{(k)}\) - aceste calcule pot fi rescrise în formă matriceală în felul următor (devine cumva evidentă această scriere odată ce se atrage atenția asupra vectorilor auxiliari obișnuiți, \(\bm{t_1}\) și \(\bm{t_{21}}\)):

\[ \left(\tilde{Y}^{(k)},R\right)=QR\left(A\begin{bmatrix}\bm{y_1'}^{(k-1)} & \bm{y_2'}^{(k-1)}\end{bmatrix}\right) = QR\left(A\tilde{Y}^{(k-1)}\right) \]

Vom folosi matricea \(\Lambda^{(k)}\in\mathbb{R}^2\) pentru a dispune valorile proprii calculate la a \(k\)-a iterație pe diagonala principală. Astfel, putem restrânge calcularea valorilor \(\lambda_1\), respectiv \(\lambda_2\), la următorul produs:

\[ \begin{bmatrix} \lambda_1^{(k)} & 0\\ 0 & \lambda_2^{(k)}\\ \end{bmatrix} = \begin{bmatrix} \bm{y_1'}^{(k)}\\ \bm{y_2'}^{(k)} \end{bmatrix}A\begin{bmatrix} \bm{y_1'}^{(k)}&\bm{y_2'}^{(k)} \end{bmatrix} \Rightarrow \boxed{\Lambda^{(k)}=\left(\tilde{Y}^{(k)}\right)^TA\left(\tilde{Y}^{(k)}\right)} \]

Acceptăm (fără demonstrație) că această metodă poate fi generalizată pentru a afla oricâte dintre primele \(p\in\overline{1,n}\) valori proprii ale matricei \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\).

Sinteza algoritmului#

Vom ilustra în câteva rânduri adevărata esență a algoritmului general:

Matricea \(Y^{(0)}=\begin{bmatrix}\bm{y_1}^{(0)} & \bm{y_2}^{(0)} & \dots & \bm{y_p}^{(0)}\end{bmatrix}\) primește o aproximare de început - aceasta poate fi complet aleatoare;
Se calculează \(\left(\tilde{Y}^{(0)},R\right)=QR\left(Y^{(0)}\right)\) (valoarea lui \(R\) nu este utilă);
Pentru fiecare pas \(k\geq 1\), se va obține o aproximare tot mai bună astfel:
- Se calculează: \(\left(\tilde{Y}^{(k)},R\right)=QR\left(A\tilde{Y}^{(k-1)}\right)\) (se ignoră valoarea \(R\), aceasta nu este utilă);
- Se aproximează valorile proprii, adică \(\Lambda^{(k)}=\left(\tilde{Y}^{(k)}\right)^TA\left(\tilde{Y}^{(k)}\right)\).

Este necesară o modalitate de a verifica pe parcurs că eroarea relativă este suficient de mică astfel încât să se poată încheia la un moment dat algoritmul. Acest lucru se poate implementa în mai multe modalități și nu vom propune noi o soluție.

Metoda QR (fără deplasare)#

Algoritmul anterior calculează pe parcurs matrice \(R\) pe care nu le folosește la nimic. Există însă o metodă QR ce se folosește de aceste valori. Vom explica mai multe detalii în paragrafele de mai jos.

Sinteza algoritmului#

Vom ilustra mai întâi algoritmul:

Pornim de la inițializările \(\Lambda^{(0)}=A\), respectiv \(Y^{(0)}=Q^{(0)}=R^{(0)}=I_n\);
Pentru fiecare pas \(k\geq 1\), se va obține o aproximare tot mai bună astfel:

Se calculează: \(\left(Q^{(k)},R^{(k)}\right)=QR\left(\Lambda^{(k-1)}\right)\);
Se aproximează valorile proprii, adică \(\Lambda^{(k)}=R^{(k)}Q^{(k)}\);
Se calculează vectorii proprii, deci \(Y^{(k)}=Y^{(k-1)}Q^{(k)}\).

Evident, se va folosi o tehnică de calcul a erorii relative pentru a încheia algoritmul.

Disecarea algoritmului#

Nu vom intra prea mult în detaliile algoritmului, însă vom prezenta câteva proprietăți întâmpinate pe parcurs. Pentru detalii suplimentare, consultați subcapitolul Lectură Suplimentară (sau direct cartea lui Jun-Lu[27]). Cu toate acestea...

Deoarece \(\left(Q^{(k)},R^{(k)}\right)=QR\left(\Lambda^{(k-1)}\right)\), știm că \(Q^{(k)}R^{(k)}=\Lambda^{(k-1)}\), deci \(R^{(k)}=\left(Q^{(k)}\right)^{-1}\Lambda^{(k-1)}\). Înlocuind în calculul valorilor proprii, se obține că:

\[ \Lambda^{(k)}=\left(Q^{(k)}\right)^{-1}\Lambda^{(k-1)}Q^{(k)} \]

Așadar, folosind definiția transformărilor de asemănare, observăm că \(\Lambda^{(k)}\) este o transformare de asemănare a lui \(\Lambda^{(k-1)}\). Prin continuarea acestui raționament, se va ajunge la faptul că \(\Lambda^{(k)}\) este o transformare de asemănare a lui \(\Lambda^{(0)}\), adică a lui \(A\), ceea ce înseamnă că \(\Lambda^{(k)}\) și \(A\) vor avea aceleași valori proprii.

Suplimentar, utilizând ecuația de mai devreme, putem rescrie \(\Lambda^{(k)}\) în felul următor:

\[\begin{split} \Lambda^{(k)}&=\left(Q^{(k)}\right)^{T}\Lambda^{(k-1)}Q^{(k)}\\ &=\left(Q^{(k)}\right)^{T}\left(Q^{(k-1)}\right)^{T}\Lambda^{(k-2)}Q^{(k-1)}Q^{(k)}\\ &=\dots\\ &=\left[\left(Q^{(k)}\right)^{T}\dots\left(Q^{(0)}\right)^{T}\right]\cdot A\cdot\left[Q^{(0)}\dots Q^{(k)}\right]\\ &=\left(Y^{(k)}\right)^TAY^{(k)} \end{split}\]

Am ajuns așadar la aceeași formulă pe care am prezentat-o inițial, ceea ce ar trebui să valideze corectitudinea algoritmului.

Metoda QR (cu deplasare)#

Privind lucrurile din perspectiva metodei puterii inverse, singura modificare față de algoritmul nemodificat intervine prin utlizarea coeficientului \(\mu\in\mathbb{R}\) acolo unde apare calcul ce îl conține pe \(\Lambda^{(k)}\).

Sinteza algoritmului#

Vom ilustra mai întâi algoritmul:

Pornim de la inițializările \(\Lambda^{(0)}=A\), respectiv \(Y^{(0)}=Q^{(0)}=R^{(0)}=I_n\);
Pentru fiecare pas \(k\geq 1\), se va obține o aproximare tot mai bună astfel:
- Se calculează deplasarea la acest pas, de obicei \(\mu^{(k)}=A^{(k-1)}_{ii}\);
- Se calculează: \(\left(Q^{(k)},R^{(k)}\right)=QR\left(\Lambda^{(k-1)}-\mu^{(k)}I_n\right)\);
- Se aproximează valorile proprii, adică \(\Lambda^{(k)}=R^{(k)}Q^{(k)}+\mu^{(k)}I_n\);
- Se calculează vectorii proprii, deci \(Y^{(k)}=Y^{(k-1)}Q^{(k)}\).

Evident, se va folosi o tehnică de calcul a erorii relative pentru a încheia algoritmul.

Licență#

The book "Metode Numerice", written by Valentin-Ioan Vintilă, is licensed under CC BY-NC-SA 4.0