Descompunerea QR

De la Wikipedia, enciclopedia liberă

Salt la: Navigare, căutare

În algebra liniară, descompunerea QR (numită şi factorizarea QR) a unei matrice este o descompunere a acelei matrice într-un produs dintre o matrice ortogonală şi una triunghiulară. Descompunerea QR este adesea folosită pentru a rezolva problema celor mai mici pătrate. Descomunerea QR stă şi la baza unui algoritm de aflare a valorilor proprii, algoritmul QR.

Cuprins

[modifică] Definiţie

O descomunere QR a unei matrici pătrate reale A este o descompunere a lui A de forma

 A = QR, \,

unde Q este o matrice ortonormală (cu proprietatea că QTQ = I ) şi R este o matrice superior triunghiulară. Analog, se pot defini descompunerile QL, RQ şi LQ ale lui A.

Mai general, se poate factoriza o matrice complexă m×n (cu mn) sub forma unui produs dintre o matrice unitară m×n (în sensul că QQ = I ) şi o matrice n×n superior triunghiulară.

Dacă A este nesingulară, atunci această factorizare este unică dacă se pune condiţia ca elementele diagonale ale lui R să fie pozitive.

[modifică] Calculul descompunerii QR

Există câteva metode pentru calculul efectiv al descompunerii QR, cum ar fi cele cu ajutorul procedeului Gram–Schmidt, transformărilor Householder, sau al rotaţiilor Givens. Fiecare metodă are avantaje şi dezavantaje.

[modifică] Descompunerea QR prin procedeul Gram-Schmidt

Se consideră procedeul Gram–Schmidt, unde vectorii consideraţi în procedeu sunt coloanele matricei A=(\mathbf{a}_1| \cdots|\mathbf{a}_n). Se defineşte \mathrm{proj}_{\mathbf{e}}\mathbf{a} 
= \frac{\left\langle\mathbf{e},\mathbf{a}\right\rangle}{\left\langle\mathbf{e},\mathbf{e}\right\rangle}\mathbf{e}
unde \left\langle\mathbf{v},\mathbf{w}\right\rangle
=\mathbf{v}^T\mathbf{w}.

Atunci


\mathbf{u}_1 = \mathbf{a}_1, \qquad\mathbf{e}_1 = {\mathbf{u}_1 \over \|\mathbf{u}_1\|}

\mathbf{u}_2 = \mathbf{a}_2-\mathrm{proj}_{\mathbf{e}_1}\,\mathbf{a}_2, \qquad\mathbf{e}_2 = {\mathbf{u}_2 \over \|\mathbf{u}_2\|}

\mathbf{u}_3 = \mathbf{a}_3-\mathrm{proj}_{\mathbf{e}_1}\,\mathbf{a}_3-\mathrm{proj}_{\mathbf{e}_2}\,\mathbf{a}_3, \qquad\mathbf{e}_3 = {\mathbf{u}_3 \over \|\mathbf{u}_3\|}
\vdots

\mathbf{u}_k = \mathbf{a}_k-\sum_{j=1}^{k-1}\mathrm{proj}_{\mathbf{e}_j}\,\mathbf{a}_k,\qquad\mathbf{e}_k = {\mathbf{u}_k\over\|\mathbf{u}_k\|}

Atunci se rearanjează ecuaţiile de mai sus astfel încât \mathbf{a}_is să fie în stânga şi rezultă următoarele ecuaţii.

\mathbf{a}_1 = \mathbf{e}_1\|\mathbf{u}_1\|
\mathbf{a}_2 = \mathrm{proj}_{\mathbf{e}_1}\,\mathbf{a}_2+\mathbf{e}_2\|\mathbf{u}_2\|
\mathbf{a}_3 = \mathrm{proj}_{\mathbf{e}_1}\,\mathbf{a}_3+\mathrm{proj}_{\mathbf{e}_2}\,\mathbf{a}_3+\mathbf{e}_3\|\mathbf{u}_3\|
\vdots
\mathbf{a}_k = \sum_{j=1}^{k-1}\mathrm{proj}_{\mathbf{e}_j}\,\mathbf{a}_k+\mathbf{e}_k\|\mathbf{u}_k\|


Se observă că deoarece \mathbf{e}_i sunt vectori unitate, avem următoarele.

\mathbf{a}_1 = \mathbf{e}_1\|\mathbf{u}_1\|
\mathbf{a}_2 = \left\langle\mathbf{e}_1,\mathbf{a}_2\right\rangle\mathbf{e}_1
+\mathbf{e}_2\|\mathbf{u}_2\|
\mathbf{a}_3 = \left\langle\mathbf{e}_1,\mathbf{a}_3\right\rangle\mathbf{e}_1
+\left\langle\mathbf{e_2},\mathbf{a}_3\right\rangle\mathbf{e}_2
+\mathbf{e}_3\|\mathbf{u}_3\|
\vdots
\mathbf{a}_k = \sum_{j=1}^{k-1}\left\langle\mathbf{e}_j,\mathbf{a}_k\right\rangle\mathbf{e}_j
+\mathbf{e}_k\|\mathbf{u}_k\|

Aceste ecuaţii pot fi scrise sub formă matriceală după cum urmează.

\left(\mathbf{e}_1\left|\ldots\right|\mathbf{e}_n\right)
\begin{pmatrix} 
\|\mathbf{u}_1\| & \langle\mathbf{e}_1,\mathbf{a}_2\rangle &  \langle\mathbf{e}_1,\mathbf{a}_3\rangle  & \ldots \\
0                & \|\mathbf{u}_2\|                        &  \langle\mathbf{e}_2,\mathbf{a}_3\rangle  & \ldots \\
0                & 0                                       & \|\mathbf{u}_3\|                          & \ldots \\
\vdots           & \vdots                                  & \vdots                                    & \ddots \end{pmatrix}

Dar produsul fiecărui rând şi coloană al matricelor de mai sus ne dau o coloană corespunzătoare a matricei A iniţiale, şi împreună, ne dau matricea A, deci am factorizat pe A într-o matrice ortogonală Q (matricea formată din ek), via Gram Schmidt, şi evident, matricea superior triunghiulară este restul R.

Altfel, \begin{matrix} R \end{matrix} poate fi calculată după cum urmează:

Dat fiind că 
\begin{matrix}Q\end{matrix} = \left(\mathbf{e}_1\left|\ldots\right|\mathbf{e}_n\right).
avem


\begin{matrix} R = Q^{T}A = \end{matrix} 
\begin{pmatrix} 
\langle\mathbf{e}_1,\mathbf{a}_1\rangle & \langle\mathbf{e}_1,\mathbf{a}_2\rangle &  \langle\mathbf{e}_1,\mathbf{a}_3\rangle  & \ldots 
\\
0                & \langle\mathbf{e}_2,\mathbf{a}_2\rangle                        &  \langle\mathbf{e}_2,\mathbf{a}_3\rangle  & \ldots 
\\
0                & 0                                       & \langle\mathbf{e}_3,\mathbf{a}_3\rangle                          & \ldots 
\\
\vdots           & \vdots                                  & \vdots                                    & \ddots \end{pmatrix}.

Se observă că \langle\mathbf{e}_j,\mathbf{a}_j\rangle = \|\mathbf{u}_j\|, \langle\mathbf{e}_j,\mathbf{a}_k\rangle = 0 \mathrm{~~pentru~~} j > k, şi QQT = I, deci QT = Q − 1.

[modifică] Exemplu

Se cere descompunerea lui

A = 
\begin{pmatrix}
12 & -51 & 4 \\
6 & 167 & -68 \\
-4 & 24 & -41 
\end{pmatrix}
.

Fie matricea ortogonală Q astfel încât


\begin{matrix}
 Q\,Q^{T} = I.
\end{matrix}

Atunci putem calcula Q prin Gram-Schmidt astfel:


U = 
\begin{pmatrix}
\mathbf u_1 & \mathbf u_2 & \mathbf u_3
\end{pmatrix}
=
\begin{pmatrix}
12 & -69 & -58/5 \\
6  & 158 & 6/5 \\
-4 &  30 & -33 
\end{pmatrix};

Q = 
\begin{pmatrix}
\frac{\mathbf u_1}{\|\mathbf u_1\|} & 
\frac{\mathbf u_2}{\|\mathbf u_2\|} & 
\frac{\mathbf u_3}{\|\mathbf u_3\|}
\end{pmatrix}
=
\begin{pmatrix}
     6/7    &    -69/175   &   -58/175   \\
     3/7    &    158/175   &     6/175   \\
    -2/7    &      6/35    &   -33/35    
\end{pmatrix};

Deci avem:


\begin{matrix}
 A = Q\,Q^{T}A = Q R; 
\end{matrix}

\begin{matrix}
 R = Q^{T}A =
\end{matrix}
\begin{pmatrix}
    14  &  21          &            -14 \\
     0  & 175          &            -70 \\
     0  &   0          &             35
\end{pmatrix}.

Efectuând operaţia cu ajutorul MATLAB, admiţând erorile de rotunjire datorate operaţiilor cu precizie finită, se obţine: 
\begin{matrix}
 Q = 
\end{matrix}
\begin{pmatrix}
         0.857142857142857     &   -0.394285714285714  &      -0.331428571428571 \\
         0.428571428571429     &    0.902857142857143  &       0.034285714285714 \\
        -0.285714285714286     &    0.171428571428571  &      -0.942857142857143 
\end{pmatrix};


\begin{matrix}
 R = 
\end{matrix}
\begin{pmatrix}
                        14  &                      21           &            -14 \\
     1.11022302462516 \times 10^{-16}  &                     175           &            -70 \\
    -1.77635683940025 \times 10^{-15}  &  -5.32907051820075 \times 10^{-14}           &             35
\end{pmatrix}.

[modifică] Calculul QR cu ajutorul reflectorilor Householder

Un reflector Householder (sau transformare Householder) este o transformare operată asupra unui vector pe care îl reflectă faţă de un plan. Putem folosi această proprietate pentru a calcula factorizarea QR a unei matrice.

Q poate fi folosită pentru a reflecta un vector în aşa fel încât dispar toate coordonatele mai puţin una.

Fie \mathbf{x} un vector-coloană arbitrar m-dimensional cu proprietatea că ||\mathbf{x}|| = |α| pentru un scalar α. Dacă algoritmul este implementat folosind aritmetica în virgulă mobilă, atunci α trebuie să aibă semnul opus primei coordonate a lui \mathbf{x} pentru a evita pierderea de semnificaţie. Dacă \mathbf{x} e un vector complex, atunci definiţia

 \alpha = - \mathrm{e}^{\mathrm{i} \arg x_1} \|\mathbf{x}\|

ar trebui să fie utilizată (Stoer,Bulirsch,2002,p.225).

Atunci, unde \mathbf{e}_1 este vectorul (1,0,...,0)T, şi ||·|| norma euclidiană, fie

\mathbf{u} = \mathbf{x} - \alpha\mathbf{e}_1,
\mathbf{v} = {\mathbf{u}\over\|\mathbf{u}\|},
Q = I - 2 \mathbf{v}\mathbf{v}^T.

Q este o matrice Householder şi

Qx = (\alpha, 0, \cdots, 0)^T.\,

Aceasta se poate folosi treptat pentru a transforma o matrice m-pe-n A în forma superior triunghiulară. Întâi, se înmulţeşte A cu matricea Householder Q1 obţinută prin alegerea primei coloane pentru x. Aceasta are ca rezultat o matrice QA cu zerouri în coloana din stânga (cu excepţia primului rând).

Q_1A = \begin{bmatrix}
                   \alpha_1&\star&\dots&\star\\
                      0    &     &     &    \\
                   \vdots  &     &  A' &    \\
                      0    &     &     & \end{bmatrix}

Aceasta se poate repeta pentru A′ (obţinută din Q1A ştergând primul rând şi prima coloană), având ca rezultat o matrice Householder Q2. Se observă că Q2 este mai mică decât Q1. Deoarece este de dorit ca ea să opereze asupra lui Q1A în loc de A′ trebuie să fie extinsă spre sus şi stânga, completând-o cu un 1, sau în general:

Q_k = \begin{pmatrix}
                  I_{k-1} & 0\\
                   0  & Q_k'\end{pmatrix}.

După t iteraţii ale acestui proces, t = min(m − 1,n),

 R = Q_t \cdots Q_2Q_1A

este o matrice superior triunghiulară. Deci, cu

 Q = Q_1Q_2 \cdots Q_t

A = QR is a QR decomposition of A.

Aceasta metodă are o stabilitate numerică superioară metodei Gram-Schmidt descrisă mai sus.

Tabelul următor dă numărul de operaţii în pasul k al descompunerii QR prin transformări Householder, presupunând o matrice pătrată de dimensiune n.

Operaţie Număr de operaţii în pasul k
înmulţiri 2(nk + 1)2
adunări (nk + 1)2 + (nk + 1)(nk) + 2
împărţiri 1
rădăcini pătrate 1

Adunând aceste numere pe cei (n − 1) paşi (pentru o matrice pătrată de dimensiune n), complexitatea algoritmului este dată de

\frac{4}{3}n^3+\frac{3}{2}n^2+\frac{19}{6}n-6=O(n^3)

[modifică] Exemplu

Se va calcula descompunerea matricei

A = \begin{pmatrix}
12 & -51 & 4 \\
6 & 167 & -68 \\
-4 & 24 & -41 \end{pmatrix}.

Întâi, trebuie să fie găsit un reflector care transformă prima coloană a lui A, vector \mathbf{a}_1 = (12, 6, -4)^T, în \|\mathbf{a}_1\| \;\mathrm{e}_1 = (14, 0, 0)^T.

Acum,

\mathbf{u} = \mathbf{x} - \alpha\mathbf{e}_1,

şi

\mathbf{v} = {\mathbf{u}\over\|\mathbf{u}\|},.

Aici,

α = 14 and \mathbf{x} = \mathbf{a}_1 = (12, 6, -4)^T

Deci

\mathbf{u} = (-2, 6, -4)^T and \mathbf{v} = {1 \over \sqrt{14}}(-1, 3, -2)^T, şi apoi
Q_1 = I - {2 \over \sqrt{14} \sqrt{14}} \begin{pmatrix} -1 \\ 3 \\ -2 \end{pmatrix}\begin{pmatrix} -1 & 3 & -2 \end{pmatrix}
 = I - {1 \over 7}\begin{pmatrix}
1 & -3  & 2 \\
-3 & 9 & -6 \\
2  & -6  & 4 
\end{pmatrix}
 = \begin{pmatrix}
6/7 & 3/7   &  -2/7 \\
3/7  &-2/7  &  6/7 \\
-2/7 & 6/7  &   3/7 \\
\end{pmatrix}.

Se observă că:

Q_1A = \begin{pmatrix}
14 & 21 & -14 \\
0 & -49 & -14 \\
0 & 168 & -77 \end{pmatrix},

deci avem deja o matrice aproape triunghiulară. Trebuie doar adusă la zero valoarea de pe poziţia (3, 2).

Se ia minorul (1, 1) minor, şi se aplică din nou procedeul pe

A' = M_{11} = \begin{pmatrix}
-49 & -14 \\
168 & -77 \end{pmatrix}.

Prin aceeaşi metodă ca mai sus, se obţine matricea de transformare Householder

Q_2 = \begin{pmatrix}
1 & 0 & 0 \\
0 & -7/25 & 24/25 \\
0 & 24/25 & 7/25 \end{pmatrix}

după efectuarea unei sume directe cu 1 pentru a ne asigura că următorul pas din procedeu funcţionează corect.

Se găseşte

Q=Q_1Q_2=\begin{pmatrix}
6/7 & -69/175 & 58/175\\
3/7 & 158/175 & -6/175 \\
-2/7 & 6/35 & 33/35
\end{pmatrix}
R=Q_2Q_1A=Q^\top A=\begin{pmatrix}
14 & 21 & -14 \\
0 & 175 & -70 \\
0 & 0 & -35
\end{pmatrix}.

Matricea Q este ortogonală iar R este superior triunghiulară, deci A = QR este descompunerea QR căutată.

[modifică] Bibliografie

  • Horn, Roger A.; Charles R. Johnson (1985). Matrix Analysis, Cambridge University Press. ISBN 0-521-38632-2. . Section 2.8.
  • Stoer, Josef; Roland Bulirsch (2002). Introduction to Numerical Analysis, 3rd, Springer. ISBN 0-387-95452-X. .
Unelte personale