Metoda puterii directe#

Metoda puterii directe reprezintă o modalitate rapidă de a calcula valoarea proprie predominantă a unei matrice \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\). Vom studia descompunerea spectrală, apoi vom dezvolta de la zero un algoritm ce calculează valoarea predominantă, pornind de la un model pur teoretic și "reparându-l" pentru a putea fi utilizat în practică.

Descompunerea spectrală#

Teorema descompunerii spectrale leagă cumva noțiunile explicate în subcapitolul precedent de matricele ortogonale studiate anterior.

Orice matrice pătratică \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\) este simetrică dacă și numai dacă există o matrice ortogonală \(Q\in\mathbb{R}^{n\times n}\) și o matrice diagonală \(\Lambda\in\mathbb{R}^{n\times n}\) astfel încât să se formeze transformarea de asemănare:

\[ A=Q\Lambda Q^T \]

unde \(Q=\begin{bmatrix}\bm{q_1}&\dots&\bm{q_n}\end{bmatrix}\) este formată din vectorii proprii asociați matricei \(A\) ce formează o bază ortonormată, iar matricea \(\Lambda\) se calculează utilizând spectrul \(\lambda(A)\).

Proprietăți#

Vom enumera câteva proprietăți importante ale acestui tip de descompunere:

  • O matrice simetrică are exclusiv valori proprii reale:

    \[\begin{cases}A=A^T\\\lambda\in\lambda(A)\end{cases}\Rightarrow\lambda\in\mathbb{R}\]

  • Rangul lui \(A\) este dat de numărul de valori proprii nenule.

Ridicarea la o putere#

Considerăm descompunerea spectrală \(A=Q\Lambda Q^T\). Se poate demonstra prin inducție faptul că \(A^n=Q\Lambda^n Q^T\), \(\forall n\in\mathbb{N}\).

Demonstrație

Fie propoziția matematică: \(P(n): A^n=Q\Lambda^n Q^T,\forall n\in\mathbb{N}\).

Așadar, demonstrăm prin inducție:

  • Cazul de bază. Demonstrăm \(P(0)\) adevărat:

    \[\begin{split} P(0): A^0&=Q\Lambda^0Q^T= QQ^T\text{, dar } Q^T=Q^{-1}\\ A^0&= I_n \end{split}\]

    Așadar, \(P(0)\) adevărat.

  • Pasul inductiv. Demonstrăm că \(P(k)\rightarrow P(k+1)\), \(\forall k\in\mathbb{N}\).

    Considerăm \(P(k): A^k=Q\Lambda^k Q^T\) adevărat și demonstrăm \(P(k+1)\):

    \[\begin{split} P(k+1): A^{k+1}&=A^kA=(Q\Lambda^kQ^T)(Q\Lambda Q)\text{, dar } Q^T=Q^{-1}\\ A^{k+1}&= (Q\Lambda^{k})(\Lambda Q^T)=Q\Lambda^{k+1}Q^T \end{split}\]

    Așadar, \(P(k)\rightarrow P(k+1)\).

Coeficientul Rayleigh#

Fie o matrice pătratică \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\). Dacă \(\bm{x}\in\mathbb{R}^n\) este un vector oarecare, atunci coeficientul Rayleigh se definește drept scalarul:

\[ \boxed{r(\bm{x})=\frac{\bm{x}^TA\bm{x}}{||\bm{x}||^2}=\left(\frac{\bm{x}}{||\bm{x}||}\right)^TA\left(\frac{\bm{x}}{||\bm{x}||}\right)} \]

Coeficientul pentru vectori proprii#

Interesant este că, atunci când particularizăm \(\bm{x}\in\mathbb{R}^n\) drept un vector propriu al matricei \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\), coeficientul Rayleigh \(r(\bm{x})\) va fi valoarea proprie asociată acestuia, anume \(r(\bm{x})=\lambda\).

Demonstrație

Prin definiție, \(A\bm{x}=\lambda\bm{x}\). Dacă se înmulțește cu \(\bm{x}^T\) în partea stângă a egalulului, se obține următoarea formulă pentru valoarea proprie \(\lambda\):

\[ \bm{x}^TA\bm{x}=\bm{x}^T\lambda\bm{x}\Rightarrow\lambda=\frac{\bm{x}^TA\bm{x}}{\bm{x}^T\bm{x}} \]

Numitorul acestei ecuații este chiar \(\bm{x}^T\bm{x}=||\bm{x}||^2\), așadar \(r(\bm{x})=\lambda\).

Descoperirea algoritmului#

Vrem să găsim o metodă (pur teoretică, imposibilă în practică) de a calcula valoarea proprie predominantă a unei matrice pătratice.

Fie matricea simetrică \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\), aleasă astfel încât să formeze spectrul \(\lambda(A)= \{\lambda_1,\dots,\lambda_n\}\), unde considerăm valorile sortate descrescător din punct de vedere al modulului, adică \(|\lambda_1|\gt |\lambda_2|\geq\dots\geq|\lambda_n|\).

Atenție!

Mare grijă - se acceptă o singură valoare maximă, adică \(\lambda_1\) este strict mai mare decât \(\lambda_2\), care însă poate fi mai mare sau egală față de \(\lambda_3\) etc.

Notație

În continuare, asemănător cu ce am discutat în capitolele anterioare, voi folosi notația \(\psi^{(k)}\) pentru a reprezenta a \(k\)-a iterație a lui \(\psi\), NU puterea sa.

Am stabilit anterior faptul că vectorii proprii formează baza proprie, ceea ce ne permite scrierea oricărui vector \(\bm{y}^{(0)}\) drept o combinație liniară de vectori proprii:

\[\begin{split} \bm{y}^{(0)}&=x_1\bm{q_1}+x_2\bm{q_2}+\dots+x_n\bm{q_n}\\ &=\sum_{i=1}^n x_i\bm{q_i} \end{split}\]

În această situație, \(\bm{x}=\begin{bmatrix}x_1\\\vdots\\x_n\end{bmatrix}\in\mathbb{R}^n\). Alternativ, dacă notăm \(Q=\begin{bmatrix}\bm{q_1}&\dots&\bm{q_n}\end{bmatrix}\), putem rescrie ecuația astfel:

\[ \bm{y}^{(0)}=Q\bm{x}\Leftrightarrow \bm{x}=Q^{-1}\bm{y}^{(0)} \]

Dat fiind faptul că \(A\) este o matrice reală simetrică, aceasta poate fi descompusă spectral, ceea ce presupune:

\[ A^k=Q\Lambda^k Q^T\Rightarrow A^kQ=Q\Lambda^k\Rightarrow \boxed{A^k\bm{q_i}=\lambda_i^k\bm{q_i}}\,,\,i\in\overline{1,n} \]

Ultima ecuație poate fi deci rescrisă sub forma:

\[ \boxed{\bm{q_i}=\frac{A^k}{\lambda_i^k}\bm{q_i}} \]

Considerând \(\bm{y}^{(0)}\) aproximația inițială a vectorului propriu asociat valorii proprii predominante pe care o căutăm, respectiv \(\bm{y}^{(k)}\) îmbunătățirile sale treptate, \(k\geq 1\), remarcăm următoarea ecuație:

\[ \boxed{\bm{y}^{(k)}=\frac{A}{\lambda_1}\bm{y}^{(k-1)}} \]

Folosind această relație, putem obține oricând valoarea proprie asociată utilizând coeficientul Rayleigh - valoarea \(\bm{y}^{(k)}\) încearcă să aproximeze un vector propriu, așadar, aplicând \(r\left(\bm{y}^{(k)}\right)\), se va obține o îmbunătățire a lui \(\lambda_1\).

Această recurență este însă impractică deoarece ar presupune cunoscută valoare lui \(\lambda_1\) (pe care încercăm să îl calculăm!) încă de la pasul inițial.

Modificarea pe care o vom aduce pentru a repara algoritmul va fi normalizarea valorii. Cu alte cuvinte, definim:

\[ \boxed{\bm{y'}^{(k)}=\frac{\bm{y}^{(k)}}{||\bm{y}^{(k)}||}=\frac{A\bm{y}^{(k-1)}/\lambda_1}{||A\bm{y}^{(k-1)}/\lambda_1||}=\frac{A\bm{y}^{(k-1)}}{||A\bm{y}^{(k-1)}||}\color{red}=\frac{A\bm{y'}^{(k-1)}}{||A\bm{y'}^{(k-1)}||}} \]

Deoarece \(\bm{y'}^{(k)}\) va avea norma euclidiană \(1\), calcularea valorii proprii \(\lambda_1\) se reduce la:

\[ \boxed{\lambda_1^{(k)}=r\left(\bm{y'}^{(k)}\right)=\frac{\left(\bm{y'}^{(k)}\right)^TA\left(\bm{y'}^{(k)}\right)}{\left|\left|\bm{y'}^{(k)}\right|\right|^2}=\left(\bm{y'}^{(k)}\right)^TA\left(\bm{y'}^{(k)}\right)} \]

Sinteza algoritmului#

Facem notația \(\bm{y^{(k)}}\rightarrow\bm{t}\).

Foarte pe scurt, algoritmul descris anterior se rezumă la acești pași:

  1. \(\bm{y}^{(0)}\) primește o aproximare de început;

  2. Pentru fiecare pas \(k\geq 1\), se va obține o aproximare tot mai bună astfel:

    • Se calculează un vector auxiliar: \(\bm{t}=A\bm{y}^{(k-1)}\);

    • Vectorul acesta normalizat va fi chiar \(\bm{y'}^{(k)}\), adică: \(\bm{y'}^{(k)}=\dfrac{\bm{t}}{||\bm{t}||}\);

    • Se calculează noua aproximare a lui \(\lambda_1\) folosind coeficientul Rayleigh, adică:

      \[\begin{split}\lambda_1^{(k)}&=r\left(\bm{y'}^{(k)}\right)=\frac{\left(\bm{y'}^{(k)}\right)^TA\left(\bm{y'}^{(k)}\right)}{\left|\left|\bm{y'}^{(k)}\right|\right|^2}\\&=\left(\bm{y'}^{(k)}\right)^TA\left(\bm{y'}^{(k)}\right)\end{split}\]

Pe parcursul algoritmului, se va verifica dacă eroarea relativă dintre \(\lambda_1^{(k)}\) și \(\lambda_1^{(k-1)}\) este suficient de mică, adică mai mică decât o valoare standard prestabilită, o toleranță aleasă.

Analiza convergenței#

Nu vom prezenta aici analiza convergenței metodei puterii directe. Cu toate acestea, două resurse de nădejde ce explica convergența sunt cărțile lui Jun Lu[26][27], evidențiate și în subcapitolul dedicat.

Cod ilustrativ#

Implementarea metodei puterii directe, sub forma funcției mpd, va utiliza notațiile:

  • A - matricea pentru care se caută valoarea proprie dominantă;

  • y0 - aproximația inițială pentru vectorul propriu dominant, echivalentul lui \(y_0\);

  • tol - toleranța acceptată, adică diferența minimă relativă dintre soluțiile a două iterații consecutive pentru ca algoritmul să continue;

  • max_iter - numărul maxim de iterații pe care îl poate executa algoritmul înainte să renunțe.

Funcția returnează tripletul [lambda, y, steps], alcătuit din:

  • lambda - valoarea proprie predominantă, echivalentul lui \(\lambda_1\);

  • y - vectorul asociat valorii proprii predominante;

  • steps - numărul de pași până se ajunge la convergență. 

function [lambda, y, steps] = MPD(A, y0, tol, max_iter)
    [m, n] = size(A);
    y = y0;
    lambda = inf;
    for steps = [1 : max_iter]
        t = A * y;
        y = t / norm(t);
        lambda_old = lambda;
        lambda = y' * A * y;
        err = abs((lambda - lambda_old) / lambda);
        if err < tol
            return
        end
    end
end

Licență#

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

Results will appear as you type