Factorizarea Cholesky#

Schimbăm un pic macazul și discutăm despre o factorizare LU ce necesită o prezumție suplimentară, și anume că matricea de factorizat este simetrică și pozitivă.

Dacă nu ați parcurs deja secțiunea introductivă, acum este momentul!

Explicația algoritmului#

De data aceasta, vom începe cu o matrice pătratică, simetrică și pozitiv-definită \(A\in\mathbb{R}^{n\times n}\) cu \(n\) linii și \(n\) coloane, \(n\in\mathbb{N}^*\), pentru care vrem să calculăm o factorizare LU.

Evident, fiind o matrice simetrică, \(A=A^T\), așadar, o factorizare \(A=LU\) va avea proprietatea că:

Este însă evident că, prin transpunerea unei matrice superior triunghiulare, se va obține o matrice inferior triunghiulară (și invers). Ajungem deci la concluzia că:

Așadar, factorizarea Cholesky își propune să găsească o singură matrice inferior triunghiulară \(L\in\mathbb{R}^{n\times n}\) astfel încât \(A=LL^T\):

Exemplificarea algoritmului#

Din nefericire, nu putem exemplifica algoritmul cu aceeași matrice în încercarea de a compara rezultatul factorizărilor diferite, deoarece matricea respectivă nu era simetrică.

Vom utiliza deci o nouă matrice \(A = \begin{bmatrix} 4 & 12 & 16 \\ 12 & 37 & 43 \\ 16 & 43 & 98 \end{bmatrix}\) și \(A=LL^T\) descompunerea LU de tip Cholesky, unde \(L\) este o matrice inferior triunghiulară. Se formează așadar ecuația:

Asemănător cu ce am studiat deja, calculul pe hârtie presupune executarea înmuțirii:

Sărind peste calcule, ajungem la următorul rezultat:

Deducerea algoritmului#

Înarmați cu intuiția dobândită pentru celelalte factorizări, vom observa că Cholesky este mult mai ușor de realizat. Să mai aruncăm un ochi pe descompunerea precedentă:

Se observă cu ochiul liber că avem două tipuri de valori:

• Valori aflate pe diagonala principală;

• Valori afalte sub diagonala principală.

Tot ce se află deasupra diagonalei principale poate fi ignorat, întrucât factorizarea este simetrică. Dacă urmărim dependențele variabilelor, constatăm că acestea trebuiesc tratate în ordinea:

Mult mai simplu, nu?

Să încercăm acum să scriem un algoritm pentru factorizarea Cholesky.

Sintetizarea algoritmului#

Pentru fiecare linie \(i\in\overline{1,n}\) (în ordinea lor firească), ne plimbăm de la stânga la dreapta pe fiecare element \(l_{ij}\), \(j\in\overline{1,i}\):

• Dacă \(i=j\) (valoare diagonală), obținem:

  \[l_{ii}^2=a_{ii}-\sum_{j=1}^{i-1}l_{ij}^2\]

  Aici apare o discuție referitoare la faptul că există două valori acceptate pentru \(l_{ii}\), anume \(\sqrt{l_{ii}^2}\) și \(-\sqrt{l_{ii}^2}\). Pentru simplitate, cerând ca diagonala principală să conțină doar valori pozitive, garantăm că factorizarea Cholesky este unică pentru orice matrice simetrică și pozitiv-definită:

  \[\boxed{l_{ii}=\sqrt{a_{ii}-\sum_{j=1}^{i-1}l_{ij}^2}}\]

• Dacă \(i\neq j\) (valoare non-diagonală), ajungem la următoarea relație:

  \[\boxed{l_{ij}=\frac{1}{l_{jj}}\left(a_{ij}-\sum_{k=1}^{j-1}l_{ik} l_{jk}\right)}\]

Cod ilustrativ#

Avem tot ce ne trbuie pentru a implementa factorizarea Cholesky (sub forma funcției Cholesky). Presupunând că notăm cu A matricea pe care încercăm să o factorizăm, vom folosi codul de mai jos.

Funcția returnează matricea inferior triunghiulară L, ce are ca proprietate că factorizează LU matricea \(A\), adică \(A=LL^T\).

function L = Cholesky(A)
    [n, n] = size(A);
    L = zeros(n);

    for i = [1 : n]
        % Valorile non-diagonale
        for j = [1 : i-1]
            S = 0;
            for k = [1 : j-1]
                S += L(i,k) * L(j,k);
            end
            L(i,j) = (A(i,j) - S) / L(j,j);
        end
        % Valorile diagonale
        S = 0;
        for j = [1 : i-1]
            S += L(i,j) * L(i,j);
        end
        L(i,i) = sqrt(A(i,i) - S);
    end
end

Licență#

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