Calculul descompunerii QR
modificare
Descompunerea QR prin procedeul Gram-Schmidt
modificare
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})}
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} }
⟨
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}}
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}}
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 e k superior triunghiulară este restul R .
Altfel,
R
{\displaystyle {\begin{matrix}R\end{matrix}}}
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).}
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,}
Q
Q
T
=
I
{\displaystyle QQ^{T}=I}
Q
T
=
Q
−
1
{\displaystyle Q^{T}=Q^{-1}}
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}
Q
Q
T
=
I
.
{\displaystyle {\begin{matrix}Q\,Q^{T}=I.\end{matrix}}}
Atunci putem calcula
Q
{\displaystyle Q}
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
modificare
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} }
m -dimensional cu proprietatea că ||
x
{\displaystyle \mathbf {x} }
aritmetica în virgulă mobilă , atunci α trebuie să aibă semnul opus primei coordonate a lui
x
{\displaystyle \mathbf {x} }
pierderea de semnificație . Dacă
x
{\displaystyle \mathbf {x} }
α
=
−
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}}
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}
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 Q 1 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 Q 1 A ștergând primul rând și prima coloană), având ca rezultat o matrice Householder Q ′2 . Se observă că Q ′2 este mai mică decât Q 1 . Deoarece este de dorit ca ea să opereze asupra lui Q 1 A î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}
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}
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)}
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})}
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}}
‖
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}
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}}
v
=
1
14
(
−
1
,
3
,
−
2
)
T
{\displaystyle \mathbf {v} ={1 \over {\sqrt {14}}}(-1,3,-2)^{T}}
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ă.
Horn, Roger A. (1985 ). Matrix Analysis . Cambridge University Press. ISBN 0-521-38632-2 . Stoer, Josef (2002 ). Introduction to Numerical Analysis (ed. 3rd). Springer. ISBN 0-387-95452-X . 
