Factorizarea Gram-Schmidt#

O modalitate rapidă de a descompune o matrice în format \(A=QR\) este algoritmul Gram-Schmidt. Pentru a-l înțelege însă, este nevoie de o abordare mai matematică.

Explicația vectorială#

Să considerăm o mulțime de vectori din subspațiul \(\mathbb{R}^n\), \(n\in\mathbb{N}^*\), ce formează o bază \(B=\{\bm{a_1},\dots,\bm{a_n}\}\).

Bază

Amintim că o mulțime de vectori \(B=\{\bm{a_1},\dots,\bm{a_n}\}\) formează o bază pentru subspațiul alcătuit din vectorii \(\bm{u_1},\dots,\bm{u_m}\) dacă și numai dacă oricare dintre aceștia se poate scrie drept o combinație liniară folosind numai vectori din \(B\).

Riguros, \(\exists\alpha_1,\dots,\alpha_n\in\mathbb{R}\) astfel încât:

\[\bm{u_i}=\sum_{k=1}^n\alpha_k\bm{a_k},\,\forall i\in\overline{1,m}\]

Factorizarea Gram-Schmidt își propune să genereze o bază ortonormată \(Q=\{\bm{q_1},\dots,\bm{q_n}\}\) în funcție de baza \(B\), adică o bază cu proprietatea că fiecare vector component este ortonormat relativ față de ceilalți.

Fiecare vector \(\bm{a_i}\) din baza \(B\), \(\forall i\in\overline{1,n}\), va trece printr-un proces de transformare, descris foarte simplist prin intermediul acestor doi pași:

\[\boxed{\bm{a_i}\to\bm{v_i}\to\bm{q_i}}\,,\,\,\forall i\in\overline{1,n}\]

Am notat "\(\rightarrow\)" operația "trece în" - cu alte cuvinte, \(\bm{a_i}\) trece în \(\bm{v_i}\), care trece la rândul său în \(\bm{q_i}\). Vom vedea imediat cine sunt \(\bm{v_i}\) și \(\bm{q_i}\), important este să înțelegeți că întreaga procedură este algoritmică, deci va putea fi adaptată ulterior pentru un calculator.

Pentru ca această transformare să se poată realiza cu succes, trebuie să se respecte următoarele reguli:

  • Fiecare vector \(\bm{a_i}\), \(\forall i\in\overline{1,n}\), se va transforma pe rând, în ordinea crescătoare a indicilor; deci întâi \(\bm{a_i}\) și abia apoi \(\bm{a_j}\) dacă \(i\gt j\);

  • Vectorul \(\bm{v_i}\) va reprezenta un vector creat folosind \(\bm{a_i}\) și care este ortogonal cu fiecare vector \(\bm{v_j}\) aflat înaintea sa, adică \(\bm{v_i}\perp\bm{v_j}\), \(\forall i\in\overline{2,n}\) și \(j\in\overline{1,i-1}\).

    Bineînțeles, excepție de la această regulă o va face \(\bm{v_1}\), căci acesta nu are față de cine să fie ortogonal. Din acest motiv, vom considera chiar \(\bm{v_1} = \bm{a_1}\);

  • Vectorul \(\bm{q_i}\) va reprezenta normalizarea lui \(\bm{v_i}\), așadar \(\bm{q_i}=\frac{\bm{v_i}}{||\bm{v_i}||}\).

Corolar: \(\bm{v_i}\perp\bm{q_j}\)#

Deoarece vectorul \(\bm{v_i}\) va fi prin construcție ortogonal cu \(\bm{v_j}\), rezultă că acesta va fi ortogonal și cu \(\bm{q_j}\), adică:

\[ \boxed{\bm{v_i}\perp\bm{v_j}\Leftrightarrow\bm{v_i}\perp\bm{q_j}}\,,\,\, \forall i\in\overline{2,n},\,j\in\overline{1,i-1} \]
Demonstrație

Doi vectori \(\bm{v_i}\) și \(\bm{v_j}\), \(\forall i\in\overline{1,n},\,j\in\overline{1,i-1}\), sunt ortogonali dacă și numai dacă:

\[ \left\langle \bm{v}_i, \bm{v}_j\right\rangle = 0 \]

Cunoaștem faptul că \(\bm{q_j}\) reprezintă normalizarea lui \(\bm{v_j}\), adică:

\[ \bm{q_j}=\frac{\bm{v_j}}{||\bm{v_j}||} \Rightarrow \bm{v_j} = ||\bm{v_j}||\cdot \bm{q_j} \]

Înlocuind în ecuația \(\left\langle \bm{v}_i, \bm{v}_j\right\rangle = 0\) și folosind proprietățile produsului scalar, se obține că:

\[0=\left\langle \bm{v_i}, \bm{v_j}\right\rangle=\left\langle \bm{v_i},\,||\bm{v_j}||\cdot \bm{q_j}\right\rangle = ||\bm{v_j}||\cdot\left\langle \bm{v_i}, \bm{q_j}\right\rangle\]

Dar \(||\bm{v_j}||\neq0\), deci \(\langle \bm{v_i}, \bm{q_j}\rangle=0\), așadar \(\bm{v_i}\perp \bm{q_j}\).

Exemplificarea algoritmului#

Cea mai simplă modalitate de a înțelege algoritmul este să pornim de la un exemplu. Să considerăm baza \(B=\left\{\begin{bmatrix} 3 \\ 4 \\ 0 \end{bmatrix}, \begin{bmatrix} 6 \\ 0 \\ 8 \end{bmatrix}, \begin{bmatrix} 0 \\ 7 \\ 0 \end{bmatrix}\right\}\) pentru subspațiul vectorial \(\mathbb{R}^3\).

Ne propunem să ajungem la baza ortonormată \(Q=\{\bm{q_1},\bm{q_2},\bm{q_3}\}\), ținând cont de ordinea \(\bm{a_i}\to\bm{v_i}\to\bm{q_i}\) pe care am descris-o anterior.

  • Pasul 1. Calcularea lui \(\bm{q_1}\):

    Fiind excepție de la regulă, se acceptă \(\bm{v_1}=\bm{a_1}=\begin{bmatrix} 3 \\ 4 \\ 0 \end{bmatrix}\).

    Pentru a normaliza vectorul, se calculează norma sa euclidiană, \(||\bm{v_1}||=\sqrt{3^2+4^2+0^2}=5\), și se împarte \(\bm{v_1}\) la aceasta, adică:

    \[ \bm{q_1}=\frac{\bm{v_1}}{||\bm{v_1}||}=\begin{bmatrix} 3/5 \\ 4/5 \\ 0 \end{bmatrix}=\begin{bmatrix} 0.6 \\ 0.8 \\ 0 \end{bmatrix} \]

  • Pasul 2. Calcularea lui \(\bm{q_2}\):

    Valoarea lui \(\bm{v_2}\) depinde atât de \(\bm{a_2}\) (de la care vom porni), cât și de \(\bm{v_1}\), prin faptul că \(\bm{v_2}\perp\bm{v_1}\).

    Vom folosi corolarul descris anterior pentru a simplifica partea matematică, adică vom cere \(\bm{v_2}\perp \bm{q_1}\) în locul \(\bm{v_2}\perp\bm{v_1}\); cu alte cuvinte:

    \[ \bm{v_2}\perp \bm{q}_1\Rightarrow \langle \bm{v_2}, \bm{q}_1\rangle=0 \]

    Știm însă că \(\bm{a_2}\) nu este (neapărat) ortogonal cu \(\bm{q_1}\), adică există posibilitatea scrierii lui \(\bm{a_2}\) ca o combinație liniară în funcție de \(\bm{v_2}\) și \(\bm{q_1}\), sub forma:

    \[ \bm{a_2}=\bm{v_2}+\alpha_{21}\bm{q_1}\Rightarrow \bm{v_2}=\bm{a_2}-\alpha_{21}\bm{q_1} \]

    unde \(\alpha_{21}\in\mathbb{R}\) este o constantă.

    Folosind aceste două ecuații, se obține:

    \[ \langle \bm{v_2}, \bm{q_1}\rangle=0\Rightarrow \langle \bm{a_2}-\alpha_{21}\bm{q_1}, \bm{q_1}\rangle=0\Rightarrow \langle \bm{a_2}, \bm{q_1}\rangle - \alpha_{21}\cdot\langle \bm{q_1}, \bm{q_1}\rangle=0 \]

    Dar \(||\bm{q_1}||=1\), deci \(\langle \bm{q_1},\bm{q_1}\rangle=1\). Așadar, putem afla \(\alpha_{21}\):

    \[ \alpha_{21}=\langle \bm{a_2},\bm{q_1}\rangle=\bm{a_2}^T\bm{q_1}=\begin{bmatrix} 6 \\ 0 \\ 8\end{bmatrix}^T\begin{bmatrix} 0.6 \\ 0.8 \\ 0\end{bmatrix}=3.6 \]

    Se obține deci formula finală pentru \(\bm{v_2}\):

    \[ \boxed{\bm{v_2}=\bm{a_2}-\langle \bm{a_2},\bm{q_1}\rangle\cdot \bm{q_1}} \]

    Făcând înlocuirile numerice, se obține că:

    \[ \bm{v_2}=\begin{bmatrix} 6 \\ 0 \\ 8 \end{bmatrix}-3.6\cdot\begin{bmatrix} 0.6 \\ 0.8 \\ 0 \end{bmatrix}=\begin{bmatrix} 3.84 \\ -2.88 \\ 8 \end{bmatrix} \]

    Nu în ultimul rând, se calculează \(\bm{q_2}\) împărțindu-l pe \(\bm{v_2}\) la norma sa, \(||\bm{v_2}||\approx9.3295\):

    \[ \bm{q_2}=\frac{\bm{v_2}}{||\bm{v_2}||}\approx\begin{bmatrix} 0.4116\\ -0.3078\\0.8575 \end{bmatrix} \]

  • Pasul 3. Calcularea lui \(\bm{q_3}\):

    Asemănător, valoarea lui \(\bm{v_3}\) depinde atât de \(\bm{a_3}\), cât și de faptul că \(\bm{v_3}\perp \bm{v_1}\) și că \(\bm{v_3}\perp \bm{v_2}\), sau, prin corolar, \(\bm{v_3}\perp \bm{q_1}\) și \(\bm{v_3}\perp \bm{q_2}\).

    Deoarece \(\bm{a_3}\) nu este (neapărat) ortogonal cu \(\bm{q_1}\) și cu \(\bm{q_2}\), înseamnă că există posibilitatea scrierii lui \(\bm{q_3}\) sub forma unei combinații liniare în funcție de \(\bm{v_3}\), \(\bm{q_1}\) și \(\bm{q_2}\), așadar:

    \[ \bm{a_3}=\bm{v_3}+\alpha_{31}\bm{q_1}+\alpha_{32}\bm{q_2}\Rightarrow \bm{v_3}=\bm{a_3}- \alpha_{31}\bm{q_1}-\alpha_{32}\bm{q_2} \]

    unde \(\alpha_{31},\alpha_{32}\in\mathbb{R}\) sunt constante. Se obține deci sistemul:

    \[ \begin{cases} \langle \bm{v_3},\bm{q_1}\rangle = 0\\ \langle \bm{v_3},\bm{q_2}\rangle = 0 \end{cases}\Rightarrow \begin{cases} \langle \bm{a_3},\bm{q_1}\rangle - \alpha_{31}\cdot\langle \bm{q_1},\bm{q_1} \rangle - \alpha_{32}\cdot\langle \bm{q_2},\bm{q_1} \rangle = 0\\ \langle \bm{a_3},\bm{q_2}\rangle - \alpha_{31}\cdot\langle \bm{q_1},\bm{q_2} \rangle - \alpha_{32}\cdot\langle \bm{q_2},\bm{q_2} \rangle = 0\\ \end{cases} \]

    Este însă important să nu pierdem din vedere faptul că \(\bm{q_1}\perp \bm{q_2}\), deci \(\langle \bm{q_1}, \bm{q_2}\rangle = \langle \bm{q_2}, \bm{q_1}\rangle = 0\).

    Sistemul de ecuații se simplifică astfel:

    \[ \begin{cases} \alpha_{31}=\langle \bm{a_3}, \bm{q_1}\rangle=\langle\bm{a_3},\bm{q_1}\rangle=\bm{a_3}^T\bm{q_1}=\begin{bmatrix} 0 \\ 7 \\ 0\end{bmatrix}^T\begin{bmatrix} 0.6 \\ 0.8 \\ 0\end{bmatrix}=5.6 \\ \alpha_{32}=\langle \bm{a_3}, \bm{q_2}\rangle=\bm{a_3}^T\bm{q_2}\approx\begin{bmatrix} 0 \\ 7 \\ 0\end{bmatrix}^T\begin{bmatrix} 0.4116 \\ -0.3078 \\ 0.8575\end{bmatrix}\approx2.1609 \end{cases} \]

    Se obține așadar formula finală pentru \(\bm{v_2}\):

    \[ \boxed{ \bm{v_3}=\bm{a_3}-\langle \bm{a_3}, \bm{q_1}\rangle\cdot \bm{q_1}-\langle \bm{a_3}, \bm{q_2}\rangle\cdot \bm{q_2}} \]

    Făcând înlocuirile numerice, se obține că:

    \[ \bm{v_3}\approx\begin{bmatrix} 0\\7\\0 \end{bmatrix}-5.6\cdot\begin{bmatrix} 0.6\\0.8\\0 \end{bmatrix}-2.1609\cdot\begin{bmatrix} 0.4116\\-0.3078\\0.8757 \end{bmatrix}\approx\begin{bmatrix} -2.4732\\1.8568\\1.8476 \end{bmatrix} \]

    Într-un final, se calculează \(\bm{q_3}\), împărțindu-l pe \(\bm{v_3}\) la norma sa, \(||\bm{v_3}||\approx3.6025\):

    \[ \bm{q_3}=\frac{\bm{v_3}}{||\bm{v_3}||}\approx\begin{bmatrix} -0.6865\\ 0.5154\\0.5129 \end{bmatrix} \]
Concluzia exemplului

Astfel, urmând tot acest algoritm, s-a obținut din baza \(A\) baza ortonormată:

\[ Q=\left\{\begin{bmatrix} 0.6\\0.8\\0 \end{bmatrix},\begin{bmatrix} 0.4116\\-0.3078\\0.8757 \end{bmatrix},\begin{bmatrix} -0.6865\\ 0.5154\\0.5129 \end{bmatrix}\right\} \]

Un exemplu cam lung, știu, dar care ar fi trebuit să vă ajute să dobândiți o intuiție asupra modului în care va arăta algoritmul final.

Generalizarea algoritmului#

Exemplificarea anterioară nu explică (cu totul) algoritmul, așadar acesta va fi generalizat în următorii \(n\) pași (unde baza \(B\) va conține \(n\) elemente, adică \(|B|=n\)), iar fiecare pas \(i\in\overline{1,n}\) este definit astfel:

  • Pasul 1. Calcularea lui \(\bm{v_1}\) și a lui \(\bm{q_1}\):

    Deoarece \(\bm{v_1}\) este de fapt o excepție, se acceptă din start:

    \[ \bm{v_1} = \bm{a_1}\Rightarrow \bm{q_1}=\frac{\bm{v_1}}{||\bm{v_1}||} \]

  • Pasul \(i\) (\(\forall i\in\overline{2,n}\)). Calcularea lui \(\bm{v_i}\) și a lui \(\bm{q_i}\):

    Vectorul \(\bm{a_i}\) nu va fi (neapărat) ortogonal în raport cu \(\bm{q_1},\dots,\bm{q_{i-1}}\), deci poate fi scris sub forma unei combinații liniare:

    \[ \boxed{\bm{a_i}=\bm{v_i}+\sum_{k=1}^{i-1}\alpha_{ik}\bm{q_k} \Rightarrow \bm{v_i}=\bm{a_i}-\sum_{k=1}^{i-1}\alpha_{ik}\bm{q_k}} \]

    unde \(\alpha_{i1}, \dots, \alpha_{i,i-1}\in\mathbb{R}\) sunt constante. Plecând de la proprietatea impusă \(\langle \bm{v_i},\bm{v_j}\rangle = \langle \bm{v_i},\bm{q_j}\rangle = 0\), se obține sistemul:

    \[ \boxed{ \langle \bm{v_i}, \bm{q_j}\rangle = 0 \Rightarrow \langle \bm{a_i}, \bm{q_j}\rangle - \sum_{k=1}^{i-1}\alpha_{ik}\cdot\langle \bm{q_k}, \bm{q_j}\rangle = 0 }\,,\,\,\forall j\in\overline{1,i-1} \]

    Dar, întrucât \(\bm{q_i}\perp \bm{q_j}\), \(\forall i\neq j\), se obține că:

    \[ \bm{a_i},\bm{q_j}-\alpha_{ij}\langle \bm{q_j},\bm{q_j}\rangle=0\Rightarrow\boxed{\alpha_{ij}=\langle \bm{a_i},\bm{q_j}\rangle} \]

    Cunoscând \(\bm{v_i}\), este ușor de calculat \(\bm{q_i}=\frac{\bm{v_i}}{||\bm{v_i}||}\).

Explicația algoritmică#

Algoritmul prezentat anterior funcționează foarte bine pentru vectori, însă cum facem legătura cu factorizarea QR?

Fie \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\), o matrice pătratică și \(A=QR\) factorizarea sa QR.

Atunci, pentru a calcula matricea \(Q\), se procedează exact ca în algoritmul prezentat anterior, considerând baza \(B=\{\bm{a_1},\bm{a_2},\dots,\bm{a_n}\}\), unde prin \(\bm{a_i}\) se înțelege a \(i\)-a coloană a matricei \(A\), \(\forall i\in\overline{1,n}\).

Se va obține baza ortonormată \(Q=\{\bm{q_1},\bm{q_2},\dots,\bm{q_n}\}\), unde fiecare element \(\bm{q_i}\) alcătuiește coloana \(i\) a matricei \(Q\), \(\forall i\in\overline{1,n}\).

Se poate demonstra că matricea \(R=(r_{ij})_{1\leq i,j\leq n}\), pentru a lua forma unei matrice superior triunghiulare ce rezolvă \(A=QR\), va fi clădită prin elementele:

\[ r_{ij}=\begin{cases} 0,&i\gt j\\ ||\bm{v_i}||,&i=j\\ \langle \bm{a_j}, \bm{q_i}\rangle=\alpha_{ji},&i\lt j \end{cases} \]
Atenție!

Cazul \(i\lt j\) utilizează corect \(\alpha_{ji}\), NU \(\alpha_{ij}\).

Intuiție

Scrisă desfășurat, ecuația factorizării QR, în contextul factorizării Gram-Schmidt, ia forma:

\[ \begin{bmatrix} \bm{a_1} & \bm{a_2} & \cdots & \bm{a_n} \end{bmatrix}=\begin{bmatrix} \bm{q_1} & \bm{q_2} & \cdots & \bm{q_n} \end{bmatrix}\begin{bmatrix} r_{11} & r_{12} & \cdots & r_{1n} \\ 0 & r_{22} & \cdots & r_{2n} \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & r_{nn} \end{bmatrix} \]

Cu alte cuvinte, orice element \(a_j\), \(\forall j\in\overline{1,n}\), poate fi scris sub forma:

\[ \boxed{\bm{a_j}=\sum_{i=1}^jr_{ij}\bm{q_i}=\sum_{i=1}^{j-1}r_{ij}\bm{q_i}+r_{jj}\bm{q_j}}\,, \,\,\forall j\in\overline{1,n} \]

Ecuația ar trebui să vă pară cunoscută, fiind cumva parte din algoritmul prezentat anterior pentru vectori. Este așadar ușor de intuit cum au fost calculate valorile \(r_{ij}\).

Cod ilustrativ#

Pentru a implementa algoritmul Gram-Schmidt sub forma funcției gs, vom considera matricea A, anume matricea pe care dorim să o factorizăm QR.

Funcția returnează perechea [Q, R], unde am notat:

  • Q - matricea ortogonală din factorizarea QR;

  • R - matricea superior triunghiulară a factorizării QR.

Evident, se respectă egalitatea \(A=QR\). 

function [Q, R] = gs(A)
    n = length(A);
    Q = R = zeros(n);
    for i = [1 : n]
        R(1:i-1, i) = Q(1:n, 1:i-1)' * A(1:n, i);
        y = A(1:n, i) - Q(1:n, 1:i-1) * R(1:i-1, i);
        R(i, i) = norm(y);
        Q(1:n, i) = y ./ R(i, i);
    end
end

Deoarece există o corespondență directă între \(\bm{a_i}\) și \(\bm{q_i}\), este imposibil de găsit un vector \(\bm{q_i}\) pentru \(i\gt n\) (adică \(A\in\mathbb{R}^{n\times m}\) și \(m\lt n\)), așadar această metodă poate fi utilizată exclusiv pentru matricele pătratice.

Analiza complexității#

Calcularea fiecărui vector \(\bm{q_i}\) se poate realiza în timp liniar odată ce se cunoaște \(\bm{v_i}\), \(\forall i\in\overline{1,m}\).

Cu toate acestea, aflarea lui \(\bm{v_i}\) presupune cunoașterea valorilor \(\alpha_{i1},\dots,\alpha_{i,i-1}\), iar fiecare dintre aceste operații presupune un produs scalar între vectori de lungime \(n\).

Dacă notăm cu \(GS_{i}(n)\) costul de a calcula \(\alpha_{i1},\dots,\alpha_{i,i-1}\), atunci:

\[\begin{split} GS(n)&=\sum_{i=2}^nGS_i(n)=\sum_{i=2}^nn(n-i+1)\\ &=n\left[n(n-1)-\left(\frac{n(n+1)}{2}-1\right)+(n-2+1)\right]\\ &=O(n^3)-O(n^2)+O(n^2) \end{split}\]

Evident, \(GS(n)=O(n^3)\).

Concluzii#

În ciuda ușurinței cu care se înțelege și implementează acest algoritm, el rămâne nefolosit în practică, întrucât este foarte instabil din punct de vedere numeric. Detalii suplimentare vor apărea în subcapitolele următoare.

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