Descompunerea QR

Î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. Descompunerea QR stă și la baza unui algoritm de aflare a valorilor proprii, algoritmul QR.

Definiție

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

A = Q R , {\displaystyle 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 {\displaystyle m} × n {\displaystyle n} (cu mn) sub forma unui produs dintre o matrice unitară m {\displaystyle m} × n {\displaystyle n} (în sensul că QQ = I ) și o matrice n {\displaystyle n} × n {\displaystyle 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.

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.

Descompunerea QR prin procedeul Gram-Schmidt

Se consideră procedeul Gram–Schmidt, unde vectorii considerați în procedeu sunt coloanele matricei A = ( a 1 | | a n ) {\displaystyle A=(\mathbf {a} _{1}|\cdots |\mathbf {a} _{n})} . Se definește p r o j e a = e , a e , e e {\displaystyle \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 v , w = v T w {\displaystyle \left\langle \mathbf {v} ,\mathbf {w} \right\rangle =\mathbf {v} ^{T}\mathbf {w} } .

Atunci

u 1 = a 1 , e 1 = u 1 u 1 {\displaystyle \mathbf {u} _{1}=\mathbf {a} _{1},\qquad \mathbf {e} _{1}={\mathbf {u} _{1} \over \|\mathbf {u} _{1}\|}}
u 2 = a 2 p r o j e 1 a 2 , e 2 = u 2 u 2 {\displaystyle \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}\|}}
u 3 = a 3 p r o j e 1 a 3 p r o j e 2 a 3 , e 3 = u 3 u 3 {\displaystyle \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}\|}}
{\displaystyle \vdots }
u k = a k j = 1 k 1 p r o j e j a k , e k = u k u k {\displaystyle \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 a i {\displaystyle \mathbf {a} _{i}} s să fie în stânga și rezultă următoarele ecuații.

a 1 = e 1 u 1 {\displaystyle \mathbf {a} _{1}=\mathbf {e} _{1}\|\mathbf {u} _{1}\|}
a 2 = p r o j e 1 a 2 + e 2 u 2 {\displaystyle \mathbf {a} _{2}=\mathrm {proj} _{\mathbf {e} _{1}}\,\mathbf {a} _{2}+\mathbf {e} _{2}\|\mathbf {u} _{2}\|}
a 3 = p r o j e 1 a 3 + p r o j e 2 a 3 + e 3 u 3 {\displaystyle \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}\|}
{\displaystyle \vdots }
a k = j = 1 k 1 p r o j e j a k + e k u k {\displaystyle \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 e i {\displaystyle \mathbf {e} _{i}} sunt versori, avem următoarele.

a 1 = e 1 u 1 {\displaystyle \mathbf {a} _{1}=\mathbf {e} _{1}\|\mathbf {u} _{1}\|}
a 2 = e 1 , a 2 e 1 + e 2 u 2 {\displaystyle \mathbf {a} _{2}=\left\langle \mathbf {e} _{1},\mathbf {a} _{2}\right\rangle \mathbf {e} _{1}+\mathbf {e} _{2}\|\mathbf {u} _{2}\|}
a 3 = e 1 , a 3 e 1 + e 2 , a 3 e 2 + e 3 u 3 {\displaystyle \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}\|}
{\displaystyle \vdots }
a k = j = 1 k 1 e j , a k e j + e k u k {\displaystyle \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ă.

( e 1 | | e n ) ( u 1 e 1 , a 2 e 1 , a 3 0 u 2 e 2 , a 3 0 0 u 3 ) {\displaystyle \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, R {\displaystyle {\begin{matrix}R\end{matrix}}} poate fi calculată după cum urmează:

Dat fiind că Q = ( e 1 | | e n ) . {\displaystyle {\begin{matrix}Q\end{matrix}}=\left(\mathbf {e} _{1}\left|\ldots \right|\mathbf {e} _{n}\right).} avem

R = Q T A = ( e 1 , a 1 e 1 , a 2 e 1 , a 3 0 e 2 , a 2 e 2 , a 3 0 0 e 3 , a 3 ) . {\displaystyle {\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ă e j , a j = u j , {\displaystyle \langle \mathbf {e} _{j},\mathbf {a} _{j}\rangle =\|\mathbf {u} _{j}\|,} e j , a k = 0     p e n t r u     j > k , {\displaystyle \langle \mathbf {e} _{j},\mathbf {a} _{k}\rangle =0\mathrm {~~pentru~~} j>k,} și Q Q T = I {\displaystyle QQ^{T}=I} , deci Q T = Q 1 {\displaystyle Q^{T}=Q^{-1}} .

Exemplu

Se cere descompunerea lui

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

Fie matricea ortogonală Q {\displaystyle Q} astfel încât

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

Atunci putem calcula Q {\displaystyle Q} prin Gram-Schmidt astfel:

U = ( u 1 u 2 u 3 ) = ( 12 69 58 / 5 6 158 6 / 5 4 30 33 ) ; {\displaystyle 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 = ( u 1 u 1 u 2 u 2 u 3 u 3 ) = ( 6 / 7 69 / 175 58 / 175 3 / 7 158 / 175 6 / 175 2 / 7 6 / 35 33 / 35 ) ; {\displaystyle 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:

A = Q Q T A = Q R ; {\displaystyle {\begin{matrix}A=Q\,Q^{T}A=QR;\end{matrix}}}
R = Q T A = ( 14 21 14 0 175 70 0 0 35 ) . {\displaystyle {\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: Q = ( 0.857142857142857 0.394285714285714 0.331428571428571 0.428571428571429 0.902857142857143 0.034285714285714 0.285714285714286 0.171428571428571 0.942857142857143 ) ; {\displaystyle {\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}};}

R = ( 14 21 14 1.11022302462516 × 10 16 175 70 1.77635683940025 × 10 15 5.32907051820075 × 10 14 35 ) . {\displaystyle {\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}}.}

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

α = e i arg x 1 x {\displaystyle \alpha =-\mathrm {e} ^{\mathrm {i} \arg x_{1}}\|\mathbf {x} \|}

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

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

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

Q {\displaystyle Q} este o matrice Householder și

Q x = ( α , 0 , , 0 ) T . {\displaystyle 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 1 A = [ α 1 0 A 0 ] {\displaystyle Q_{1}A={\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 = ( I k 1 0 0 Q k ) . {\displaystyle Q_{k}={\begin{pmatrix}I_{k-1}&0\\0&Q_{k}'\end{pmatrix}}.}

După t {\displaystyle t} iterații ale acestui proces, t = min ( m 1 , n ) {\displaystyle t=\min(m-1,n)} ,

R = Q t Q 2 Q 1 A {\displaystyle R=Q_{t}\cdots Q_{2}Q_{1}A}

este o matrice superior triunghiulară. Deci, cu

Q = Q 1 Q 2 Q t {\displaystyle Q=Q_{1}Q_{2}\cdots Q_{t}}

A = Q R {\displaystyle A=QR} is a QR decomposition of A {\displaystyle 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 ( n k + 1 ) 2 {\displaystyle 2(n-k+1)^{2}}
adunări ( n k + 1 ) 2 + ( n k + 1 ) ( n k ) + 2 {\displaystyle (n-k+1)^{2}+(n-k+1)(n-k)+2}
împărțiri 1 {\displaystyle 1}
rădăcini pătrate 1 {\displaystyle 1}

Adunând aceste numere pe cei ( n 1 ) {\displaystyle (n-1)} pași (pentru o matrice pătrată de dimensiune n), complexitatea algoritmului este dată de

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

Exemplu

Se va calcula descompunerea matricei

A = ( 12 51 4 6 167 68 4 24 41 ) . {\displaystyle 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 a 1 = ( 12 , 6 , 4 ) T {\displaystyle \mathbf {a} _{1}=(12,6,-4)^{T}} , în a 1 e 1 = ( 14 , 0 , 0 ) T . {\displaystyle \|\mathbf {a} _{1}\|\;\mathrm {e} _{1}=(14,0,0)^{T}.}

Acum,

u = x α e 1 , {\displaystyle \mathbf {u} =\mathbf {x} -\alpha \mathbf {e} _{1},}

și

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

Aici,

α = 14 {\displaystyle \alpha =14} and x = a 1 = ( 12 , 6 , 4 ) T {\displaystyle \mathbf {x} =\mathbf {a} _{1}=(12,6,-4)^{T}}

Deci

u = ( 2 , 6 , 4 ) T {\displaystyle \mathbf {u} =(-2,6,-4)^{T}} and v = 1 14 ( 1 , 3 , 2 ) T {\displaystyle \mathbf {v} ={1 \over {\sqrt {14}}}(-1,3,-2)^{T}} , și apoi
Q 1 = I 2 14 14 ( 1 3 2 ) ( 1 3 2 ) {\displaystyle 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 7 ( 1 3 2 3 9 6 2 6 4 ) {\displaystyle =I-{1 \over 7}{\begin{pmatrix}1&-3&2\\-3&9&-6\\2&-6&4\end{pmatrix}}}
= ( 6 / 7 3 / 7 2 / 7 3 / 7 2 / 7 6 / 7 2 / 7 6 / 7 3 / 7 ) . {\displaystyle ={\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 1 A = ( 14 21 14 0 49 14 0 168 77 ) , {\displaystyle Q_{1}A={\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 = ( 49 14 168 77 ) . {\displaystyle 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 = ( 1 0 0 0 7 / 25 24 / 25 0 24 / 25 7 / 25 ) {\displaystyle 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 1 Q 2 = ( 6 / 7 69 / 175 58 / 175 3 / 7 158 / 175 6 / 175 2 / 7 6 / 35 33 / 35 ) {\displaystyle Q=Q_{1}Q_{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 2 Q 1 A = Q A = ( 14 21 14 0 175 70 0 0 35 ) . {\displaystyle R=Q_{2}Q_{1}A=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ă.

Bibliografie

  • Horn, Roger A. (). Matrix Analysis. Cambridge University Press. ISBN 0-521-38632-2. . Section 2.8.
  • Stoer, Josef (). Introduction to Numerical Analysis (ed. 3rd). Springer. ISBN 0-387-95452-X. .