M5960 Vybrané partie z aplikované matematiky – seminář
QR-
-ROZKLAD
Petr Zemánek petr.zemanek@mail.muni.cz
Obsah
1. Trocha historie . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . 2
2. Motivace . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . 3
3. QR-rozklad . . . . . .. . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . 4
4. Konstrukce QR-rozkladu . . .. . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . 5
4.1 QR-rozklad pomocí Gram-Schmidtova algoritmu . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
4.2 Householderova matice zrcadlení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
4.3 QR-rozklad pomocí Householderovy matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.4 Givensova matice rovinné rotace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.5 QR-rozklad pomocí Givensovy matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.6 Srovnání algoritmů . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5. Aplikace QR-rozkladu . . . .. . . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . 22
5.1 Řešení systémů rovnic pomocí QR-rozkladu. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2 QR-rozklad a Moore-Penroseova zobecněná inverse. . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.3 QR-rozklad a vlastní čísla matice A – QR-algoritmus . . . . . . . . . . . . . . . . . . . . . . . 24
Seznam použité literatury . . . . . .. . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . 26
1
1. Trocha historie
RUTISHAUSER, H.: Solution of Eigenvalue Problems with the LR-transformation. Nat. Bur.
Standards Appl. Math. Ser., No. 49, str. 47-81, 1958.
FRANCIS, J.G.F.: The QR Transformation: a Unitary Analogue to the LR Transformation I.
The Computer Journal, 4, str. 265-271, 1961/1962.
KUBLANOVSKAYA, V.N.: On Some Algorithms for the Solution of the Complete Eigenvalue
Problem. U.S.S.R. Comput. Math. and Math. Phys., 3, str. 637-657, 1961.
FRANCIS, J.G.F.: The QR Transformation: a Unitary Analogue to the LR Transformation II.
The Computer Journal, 4, str. 332-345, 1961/1962.
PARLETT, B.: The Development and Use of Methods of LR Type. SIAM Review, Vol. 6(3),
str. 275-295, 1964.
PARLETT, B.: Convergence of the QR Algorithm. Numerische Mathematik, 7, str. 187-193, 1965.
TOP 10
DONGARRA, J. and SULLIVAN, F.: The top ten algorithms. Computing in Science and Engineering,
2000:
• Algoritmus Metropolis pro Monte Carlo.
• Simplexová metoda pro lineární programování.
• Iterace Krylovových podprostorů.
• Dekompenzační přístup k maticovým výpočtům – Choleského metoda, QR-rozklad, LU-rozklad
s pivotem, spektrální rozklad, Schurův rozklad, singulární rozklad.
• Optimalizační kompilátor Fortranu.
• QR-algoritmus pro výpočet vlastních hodnot.
• Quick sort.
• FFT – rychlá Fourierova transformace.
• Integer relation detection.
• FMM – rychlá multipólová metoda.
2
2. Motivace
Řešení rovnic pomocí Gaussovy eliminace se dá zapsat jako rozklad na součin dvou matic LU,
kde L je dolní trojúhleníková matice a U je ve schodovitém tvaru, tzn. nulové řádky jsou až na
konci a každý nenulový řádek začíná větším počtem nul než ten předchozí, např.
2 2 3 4
0 1 0 3
0 0 0 1
0 0 0 0
– spec.
pro regulární matice je U horní trojúhelníková, přičemž L lze chápat jako součin (řádkových)
elementárních matic I(i, j)
– což popisuje právě kroky v Gaussově eliminaci.
Podmíněnost soustavy (matice):
c(A) = A · A−1
.
Řekneme, že matice A je dobře podmíněná, pokud c(A) ≈ 1, a že je špatně podmíněná, pokud
c(A) ≫ 1.
Ovšem při Gaussově eliminaci platí pro podmíněnost c(A) ≤ c(L)c(U), takže se celková podmíněnost
ještě zhorší (výjimku tvoří pozitivně definitní matice, kde A 2 = L 2
2 a c2(A) = c2
2(L)).
Tomuto se můžeme vyhnout, pokud matici A rozložíme na součin ortogonální a horní trojúhelníkové
matice, tedy na tzv. QR-rozklad.
Potom při řešení soustavy rovnic Ax = QRx = b ⇒ Rx = QT
b ⇒ x = R−1
QT
b. Neboť ortogonální
transformace zachovává velikost euklidovské normy, bude Q 2 = 1 a c2(Q) = 1, takže
A 2 = R 2 a c2(A) = c2(R). Tedy nezhoršuje se podmíněnost dvou přechodných soustav
Rx = y a Qy = b ani nedochází k růstu prvků.
3
3. QR-rozklad
Definice 1. Dvojici matic Q a R nazveme QR-rozkladem matice A, pokud platí, že
A = QR,
přičemž Q je ortogonální matice a R je horní trojúhelníková matice.
Věta 1. O existenci QR-rozkladu. K libovolné reálné matici A ∈ Rm×n
, kde m ≥ n, existuje
ortogonální matice Q ∈ Rm×m
a horní trojúhelníková matice R ∈ Rm×n
tak, že platí A = QR.
Věta 2. O jednoznačnosti QR-rozkladu. Jsou-li sloupce matice A ∈ Rm×n
, m ≥ n, lineárně
nezávislé, potom v QR-rozkladu jsou matice R a prvních n sloupců matice Q určeny až na
znaménko jednoznačně.
Důkaz: QR-rozklad existuje vždy, a jak později uvidíme, lze jej získat různými algoritmy. To ovšem
nic nemění na tom, že pokud má matice A lineárně nezávislé sloupce, pak až na znaménka existuje
právě jeden QR-rozklad.
Mějme dva takové QR-rozklady matice A, tj. nechť
A = Q1R1 = Q2R2.
Neboť platí
Q1R1 = Q2R2,
dostaneme
R1R−1
2 = Q−1
1 Q2 = QT
1 Q2.
Dále také platí, že
(R1R−1
2 )−1
= (QT
1 Q2)−1
= Q−1
2 (QT
1 )−1
= QT
2 Q1 = (QT
1 Q2)T
= (R1R−1
2 )T
,
tedy matice R1R−1
2 je ortogonální a diagonální. Taková je však pouze matice |E|. Takže
|R1R−1
2 | = |E|,
což znamená, že |R1| = |R2|, tedy i |Q1| = |Q2|.
4
4. Konstrukce QR-rozkladu
4.1 QR-rozklad pomocí Gram-Schmidtova algoritmu
Věta 3. Gram-Schmidtův QR-rozklad. K libovolné reálné matici A ∈ Rm×n
, kde m ≥ n,
existuje ortogonální matice Q ∈ Rm×m
a horní trojúhelníková matice R ∈ Rm×n
s nezápornými
prvky na diagonále tak, že platí A = QR. V případě lineárně nezávislých sloupců matice A jsou
prvky na diagonále kladné.
Základní myšlenka důkazu: Máme-li matici A ∈ Rm×n
, pak aplikací zobecněného Gram-Schmidtova
ortogonalizačního procesu na sloupce matice A (ty mohou být lineárně závislé i nezávislé) a doplněním
těchto vektorů na bázi v Rm
získáme sloupce matice Q.
Uvažujme matici A = (a1| . . .|an) složenou ze sloupcových vektorů. Pak
u1 = a1, e1 =
u1
u1
,
u2 = a2 − proj e1
a2, e2 =
u2
u2
,
u3 = a3 − proj e1
a3 − proj e2
a3, e3 =
u3
u3
,
...
uk = ak −
k−1
j=1
proj ej
ak, ek =
uk
uk
,
kde proj u v = <v, u>
<v, v>
u. Po úpravě obdržíme vzorce pro vektory ai
a1 = e1 u1 ,
a2 = proj e1
a2 + e2 u2 ,
a3 = proj e1
a3 + proj e2
a3 + e3 u3 ,
...
ak =
k−1
j=1
proj ej
ak + ek uk .
Označme Q = (e1 | . . .| en). Nyní máme
R = QT
A =







< e1, a1 > < e1, a2 > < e1, a3 > · · · < e1, an >
0 < e2, a2 > < e2, a3 > · · · < e2, an >
0 0 < e3, a3 > · · · < e3, an >
...
...
...
...
...
0 0 0 . . . < en, an >







,
neboť QQT
= E a < ej, aj >= uj , < ej, ak >= 0 pro j > k.
5
Příklad 1. Proveďme QR-rozklad matice A =


12 −51 4
6 167 −68
−4 24 −41

. Gram-Schmidtovým procesem
dostaneme
U = (u1 | u2 | u3) =


12 −69 −58
6 158 6
−4 30 −165

 .
Matici Q potom získáme jako
Q =
u1
u1
u2
u2
u3
u3
=


6/7 −69/175 −58/175
3/7 158/175 6/175
−2/7 6/35 −33/35

 .
Pak A = QQT
A = QR, takže
R = QT
A =


14 21 −14
0 175 −70
0 0 35

 .
Algoritmus?
Mějme matici A. Položme
r11 = a1 , q1 =
a1
r11
,
pro k = 2, . . . , n spočítejme:
rjk = < qj, ak > pro j = 1, . . . , k − 1,
zk = ak −
k−1
j=1
rjkqj,
r2
kn = < zk, zn >
qk =
zk
rkk
.
Metodu lze také upravit tak, že zaměníme pořadí oprací. Tedy položme
A0
≡ A.
Pak pro k = 2, . . . , n spočtěme
rkk = a
(k−1)
k 2
, qk =
a
(k−1)
k
rkk
,
rki = qT
k a
(k−1)
i pro i = k + 1, . . . , n,
A(k)
= A(k−1)
− qkrT
k .
Z formálního hlediska jde o změnu pořadí operací, ovšem z numerického hlediska obdržíme kvalitativně
různé výsledky.
6
4.2 Householderova matice zrcadlení
Definice 2. Matice tvaru
H(u) : = E −
2uuT
uT u
= E −
2uuT
u 2
se nazývá Householderova matice (někdy též elementární zrcadlení (odraz), Householderova trans-
formace).
Vlastnosti:
• označení matice zrcadlení se používá proto, že aplikujeme-li matici H(u) pro nějaké u na
vektor x ∈ Rn
, pak je vektor H(u)x souměrný s vektorem x podle nadroviny ortogonální
k vektoru v;
u(uT
x)
uT x
u
x
P
−2u(uT
x)
uT u
Hx = (I−2uuT
)x
uT u
= x−2u(uT
x)
uT u
Obrázek 1: Householderova transformace – geometrický význam.
• matice E bývá považována za speciální případ Householderovy transformace – pro u = 0 je
H(0) = E;
• Hx 2 = x 2 pro každé x ∈ Rn
, tj. zrcadlení tedy nemění délku vektoru;
• Hy = y pro každé y ∈ P = {v ∈ Rn
| vT
u = 0}.
• H má jednoduchou vlastní hodnotu -1 a (n − 1)-násobnou vlastní hodnotu 1;
Důkaz: Neboť y ∈ P = {v ∈ Rn
| vT
u = 0} má n − 1 lineárně nezávislých vektorů
y1, . . . , yn−1 a Hyi = yi pro i = 1, 2, . . . , n − 1. Takže 1 je (n − 1)-násobná vlastní
hodnota. H také zrcadlí u na −u, tj. Hu = −u. Takže -1 je vlastní hodnota matice H, která
musí být jednoduchá, neboť H má pouze n vlastních hodnot.
• z věty o spektrálním rozkladu plyne
det(H) = (−1)1 · · ·1 = −1;
7
• Matice H je ortogonální a symetrická;
Důkaz: Symetrie plyne z
HT
(u) = ET
− 2
uuT
uT u
T
= E −
2uuT
u
= H(u).
Neboť platí
H2
(u) = E −
2uuT
uT u
E −
2uuT
uT u
= E2
− 4
uuT
u 2
+ 4
uuT
uuT
u 4
= E,
je matice H(u) ortogonální.
Historická poznámka 1. Poslední tvrzení, které má velmi široký význam především v aplikacích, dokázal v článku
Unitary Triangularization of a Nonsymmetric Matrix z roku 1958 Alston S. Householder (1904-1993), americký matematik,
který se zabýval numerickou analýzou.
Věta 4. Pro každé dva vektory y, z ∈ Rn
takové, že y = z a y 2 = z 2, platí
y = H(y − z)z.
Jinými slovy, každé dva různé vektory o stejné normě lze převést jeden na druhý Householderovou
transformací.
Důkaz: Platí
H(y − z)z = E −
2(y − z)(y − z)T
y − z 2
2
z = z − 2
yT
z − z 2
2
y − z 2
2
(y − z) =
= z +
y 2
2 + z 2
2 − 2yT
z
y − z 2
2
(y − z) = z +
y − z 2
2
y − z 2
2
(y − z) = y.
Důsledek 1. Jsou-li y, z dva vektory o stejné normě, potom existuje ortogonální matice Q taková,
že y = Qz.
Důkaz: Pro y = z stačí vzít Q = H(y − z), jinak Q = E.
Věta 5. Pro každé x ∈ Rn
je
H =
H(x + sgn(x1) x 2 e1) pro x1 = x 2,
E pro x1 = x 2
ortogonální matice s vlastností
Hx = x 2 e1 .
Neboli, aplikujeme-li vhodnou matici H na vektor x, dostaneme vektor, který má všechny složky
až na první nulové.
8
Důkaz: Je-li x1 = x 2, potom z x2
1 = x2
1 + · · · + x2
n plyne, že x2 = · · · = xn = 0. Tedy x = x1 e1 =
= x 2
e1 = Ex = Hx.
Je-li x1 = x 2, potom x + sgn(x1) x 2 e1 = 0, takže vektory y = sgn(x1) x 2 e1 a z = x jsou
různé a platí pro ně y 2 = x 2 = z 2, a odsud podle Věty 4 je
y = sgn(x1) x 2 e1 = H(y − z)z = H(−x − sgn(x1) x 2 e1)x.
Poznámka 1. Pro vektor určující Householderovu matici lze volit buď + x 2 e1 nebo − x 2 e1.
Z důvodu minimalizace numerických chyb volíme stejné znaménko jako u první složky vektoru x.
Věta 6. Pro každé x takové, že x 2 = 1, je
H =
H(x + sgn(x1) e1) pro x = e1,
E pro x = e1
ortogonální matice, jejímž prvním sloupcem je vektor x.
Důkaz: Pro x = e1 je zřejmý.
Nechť tedy x = e1. Protože x 2 = 1 = e1 2, je podle Věty 5
x = H(x + sgn(x1) e1) = He1 = H•1,
což je tvrzením věty.
Díky těmto větám tedy umíme najít vektor u tak, že daný nenulový vektor x se transformuje na
vektor, který má nenulovou pouze první složku.
Příklad 2. Lze
x = (−1, −2, 7)T H(u)
−−−→ (α, 0, 0)T
?
Protože x 2 = 3
√
6, položíme u = x − x 2 e1 = (−1 − 3
√
6, −2, 7)T
a u 2 = 6(18 +
√
6 ). Dále
uuT
=


−1 − 3
√
6
−2
7

 (−1 − 3
√
6, −2, 7) =


55 + 6
√
6 2 + 6
√
6 −7 − 21
√
6
2 + 6
√
6 4 −14
−7 − 21
√
6 −14 49

 ,
takže
H(u) =
1
3(18 +
√
6 )


−1 − 3
√
6 −2 − 6
√
6 7 + 21
√
6
−2 − 6
√
6 50 + 3
√
6 14
7 + 21
√
6 14 5 + 3
√
6

 .
Snadno lze ověřit, že
H(u)x = (3
√
6, 0, 0)T
.
9
4.3 QR-rozklad pomocí Householderovy matice
Věta 7. Householderův QR-rozklad. Každou matici A ∈ Rm×n
lze pomocí s = min{n, m−1}
Householderových matic rozložit na součin QR, a to tak, že platí
Hs · · · H2H1A = QT
A =



R1
0 m > n,
(R1, 0) m < n,
R m = n.
Mějme reálnou matici A
A =





a11 · · · a1n
a21 · · · a2n
...
...
...
am1 · · · amn





.
Krok 1.: Zkonstruujme Householderovu matici H1 tak, aby H1A měla v prvním sloupci pouze
samé 0 s výjimkou pozice (1, 1), tj. aby
H1A =





⊞ · · ·
0 · · ·
...
...
0 · · ·





.
K tomu stačí získat vektor un (dle předchozího) tak, že pro
H1 = E − 2
unuT
n
uT
n un
platí
H1





a11
a21
...
am1





=





⊞
0
...
0





.
Označme A(1)
: = H1A. Ta je tvaru
A(1)
=





a11 · · ·
0 · · ·
...
...
0 · · ·





.
Krok 2.: Zkonstruujme Householderovu matici H2 tak, že H2A(1)
má ve druhém sloupci 0 „pod
pozicí (2, 2) při zachování požadavku prvního kroku, tj.
A(2)
: = H2A(1)
=







⊞ · · ·
0 ⊞ · · ·
0 0
...
...
...
0 0 · · ·







.
10
Matici H2 získáme takto: Nejdříve zkonstruujeme Householderovu matici o rozměru (m−1)×(n−1)
H2 : = En−1 − 2
un−1uT
n−1
uT
n−1un−1
takovou, že
H2





a22
a32
...
am2





=





⊞
0
...
0





,
a definujme
H2 : =





1 0 · · · 0
0
... H2
0





.
Tím získáme matici A(2)
= H2A(1)
.
Analogicky pokračujeme dále.
Pro k ≤ s.
Krok k-tý: Obecně vytváříme Householderovu matici
Hk : = En−k+1 − 2
un−k+1uT
n−k+1
uT
n−k+1un−k+1
o rozměru (m − k + 1) × (n − k + 1) takovou, že
Hk



akk
...
amk


 =





⊞
0
...
0





.
Definujeme
Hk : =
Ek−1 0
0 Hk
,
čili můžeme spočítat A(k)
= HkA(k−1)
.
Tímto způsobem po s krocích obdržíme matici A(s)
, která bude v horním trojúhelníkovém tvaru
a bude právě maticí R.
Protože
A(k)
= HkA(k−1)
k = 2, . . . , s,
máme
R = A(s)
= HsA(s−1)
= HsHs−1A(s−2)
= · · · = HsHs−1 · · ·H2H1A.
Položme
QT
= HsHs−1 · · · H2H1.
11
Máme naši hledanou ortogonální matici (neboť každá z Hi je ortogonální). Celkem
R = QT
A,
tj.
A = QR.
(Uvědomme si, že Q = HT
1 HT
2 · · · HT
s = H1H2 · · · Hs.)
Příklad 3. Uvažme matici
A =


0 1 1
1 2 3
1 1 1

 .
Krok 1.: Konstrukce H1.
H1


0
1
1

 =


⊞
0
0

 .
Potom tedy dle Příkladu 2 spočteme
u3 =


0
1
1

 +
√
2


1
0
0

 =


√
2
1
1

 ,
takže
H1 = E3 − 2
u3uT
3
uT
3 u3
=


1 0 0
0 1 0
0 0 1

 −



1 1√
2
1√
2
1√
2
1
2
1
2
1√
2
1
2
1
2


 =



0 − 1√
2
− 1√
2
− 1√
2
1
2
−1
2
− 1√
2
−1
2
1
2


 .
Určeme
A(1)
= H1A =



−
√
2 −3
√
2
2
2
√
2
0 −1−
√
2
2
−2−
√
2
2
0 −1+
√
2
2
−2+
√
2
2


 .
Krok 2.: Zkonstruujeme
H2 =
−0, 2071
−1, 2071
=
⊞
0
,
u2 =
−0, 2071
−1, 2071
− 1, 2247
1
0
=
−1, 4318
−1, 2071
,
H2 =
−0, 1691 −0, 9856
−0, 9856 0, 1691
,
tzn.
H2 =


1 0 0
0 −0, 1691 −0, 9856
0 −0, 9856 0, 1691

 ,
a spočítáme
A(2)
= H2A(1)
= H2H1A =


−1, 4142 −2, 1213 −2, 8284
0 1, 2247 1, 6330
0 0 −0, 5774

 = R
12
Pro Q nyní platí
Q = H2H1 =


0 0, 8165 0, 5774
−0, 7071 0, 4082 −0, 5774
−0, 7071 −0, 4082 0, 5774

 .
Celkem tedy
A =


0 1 1
1 2 3
1 1 1

 =


0 0, 8165 0, 5774
−0, 7071 0, 4082 −0, 5774
−0, 7071 −0, 4082 0, 5774




−1, 4142 −2, 1213 −2, 8284
0 1, 2247 1, 6330
0 0 −0, 5774

 = QR.
Příklad 4. Uvažme nyní obdélníkovou matici
A =


1 1
0, 0001 0
0 0, 0001

 .
Platí s = min{2, 2}.
Krok 1.: Konstrukce H1.
u2 =


1
0, 0001
0

 + 1 + 0, 00012


1
0
0

 =


2
0, 0001
0

 ,
H1 = E3 − 2
u2uT
2
uT
2 u2
=


−1 −0, 0001 0
−0, 0001 1 0
0 0 1


A(1)
= H1A =


−1 −1
0 −0, 0001
0 0, 0001

 .
Krok 2.: Zkonstruujeme H2.
u1 =
−0, 0001
0, 0001
− (−0, 0001)2 + 0, 00012 1
0
= 10−4 −2, 4141
0, 1000
,
H2 =
−0, 7071 0, 7071
0, 7071 0, 7071
,
H2 =


1 0 0
0 −0, 7071 0, 7071
0 0, 7071 0, 7071

 ,
A(2)
= H2A(1)
=


−1 −1
0 0, 0001
0 0

 .
Q = H2H1 =


−1 0, 0001 −0, 0001
−0, 0001 −0, 7071 0, 7071
0 0, 7071 0, 7071


13
a
R1 =
−1 −1
0 0, 0001
.
Celkem tedy
A =


1 1
0, 0001 0
0 0, 0001

 =


−1 0, 0001 −0, 0001
−0, 0001 −0, 7071 0, 7071
0 0, 7071 0, 7071




−1 −1
0 0, 0001
0 0

 = QR.
Poznámka 2. Pokud matice A nemá zaručenu plnou hodnost, je vhodné upravit výpočet tak,
abychom případné nulové (resp. „malé ) diagonální prvky rkk dostali až na závěr výpočtu. Můžeme
modifikovat QR-rozklad výběrem největšího sloupce v normě – tzv. QR-rozklad s pivotem.
Znamená to, že QR-rozklad matice AP = QR, kde P = P1P2 · · · Ps a Pk je permutační matice
výměny k-tého a j0-tého sloupce v k-tém kroku metody (v případě, že řešíme soustavu rovnic,
je třeba učinit příslušné permutace i složek vektoru řešení). V této modifikaci bude pro matici R
platit
r2
kk ≥
j
i=k
|rij|2
, j = k + 1, . . . , n,
tedy také |r11| > · · · > |rnn|. Nulové, resp. „téměř lineárně závislé sloupce se takto přesunou
na „konec matice a snadno se budou moci v případě potřeby oddělit. Později uvidíme výhody
této modifikace. Stejné výsledky obdržíme i v modifikovaném Gram-Schmidtově postupu s výběrem
pivota, pokud ještě před k-tým krokem najdeme index i0 tak, že ai0 2 = max a
(k−1)
i 2 a vyměníme
i0-tý a k-tý sloupec matice A(k−1)
.
4.4 Givensova matice rovinné rotace
Definice 3. Matice tvaru
G(i, j, c, s) : =












1 · · · 0 · · · 0 · · · 0
...
...
...
...
...
0 · · · c · · · s · · · 0
...
...
...
...
...
0 · · · −s · · · c · · · 0
...
...
...
...
...
0 · · · 0 · · · 0 · · · 1












= E + (c − 1)(ei eT
i + ej eT
j ) + s(ei eT
i − ej eT
j ),
kde c2
+ s2
= 1, se nazývá Givensova matice.
Historická poznámka 2. Matice nese jméno amerického matematika Wallace Givense (1910-1993), který se zabýval
především numerickou analýzou.
Můžeme volit c = cos Θ a s = sin Θ pro nějaké Θ. Pak značíme Givensovu matici jako G(i, j, Θ).
Geometricky můžeme matici G(i, j, Θ) chápat jako rotaci (otočení) Rn
kolem počátku o úhel
Θ v rovině (ei, ej). To je také důvod, proč se matice nazývá i Givensova matice rotace nebo
Givensova matice rovinné rotace v (ei, ej) rovině.
14
v =
cos(α − Θ)
sin(α − Θ)
=
cos Θ sin Θ
− sin Θ cos Θ
u
u =
cos α
sin α
α
Θ
ej
ei
Obrázek 2: Givensova transformace – geometrický význam.
Givensova matice je ortogonální, neboť z c2
+ s2
= 1 plyne
GT
(i, j, Θ)G(i, j, Θ) = G(i, j, Θ)GT
(i, j, Θ) = E.
Vzhledem k tvaru matice G platí, že pokud jí vynásobíme nějaký vektor
x =





x1
x2
...
xn





,
tak se změní pouze i-tá a j-tá komponenta vektoru x, ostatní se nezmění.
Pokud máme vektor x = ( x1
x2 ), pak lze snadno ověřit, že při volbě
c =
x1
x2
1 + x2
2
, s =
x2
x2
1 + x2
2
(1)
(pak totiž platí −x1 sin Θ + x2 cos Θ = 0) bude pro
G(1, 2, Θ) =
c s
−s c
platit
G(1, 2, Θ)
x1
x2
=
⊞
0
,
tedy „vynulujeme druhou složku vektoru x.
Chceme-li „vynulovat první složku vektoru x, položíme
c = −
x2
x2
1 + x2
2
, s =
x1
x2
1 + x2
2
.
Příklad 5. Nechť
x =
3
4
.
15
Položme (dle vztahu (1)) c = 3
5
a s = 4
5
. Transformací Givensovou maticí dostaneme
G(1, 2, Θ)x =
3
5
4
5
−4
5
3
5
3
4
=
1
0
.
Při násobení vektoru xT
zprava obdržíme
xT
GT
(1, 2, Θ) = 1 0 .
Při volbě c = −4
5
a s = 3
5
dostaneme
G(1, 2, Θ)x =
0
−1
.
a
xT
GT
(1, 2, Θ) = 0 −1 .
„Nulování na specifických pozicích
Givensova rotace je zvláště užitečná, pokud chceme „vytvořit nulu na určité pozici ve vektoru.
Nechť
x =














x1
x2
...
xi
...
xk
...
xn














a požadujme pouze nulu na pozici xk. Můžeme sestrojit rotaci G(i, k, Θ) (i < k) tak, že
G(i, k, Θ) · x bude mít nulu na k-té pozici.
Konstrukce: Nejdříve sestavíme 2 × 2 Givensovu matici G(i, k, Θ) = ( c s
−s c ), tak aby
c s
−s c
xi
xk
=
⊞
0
.
Tvar matice G(i, k, Θ) získáme tak, že na pozice (i, i) a (k, k) vložíme „c , na pozici (i, k)
vložíme „s a na pozici (k, i) vložíme „−s , zbytek doplníme na identickou matici.
Příklad 6. Buď
x =


1
−1
3


Sestavme Givensovu matici tak, aby „vytvořila nulu na třetí pozici, tj. k = 3.
Zvolme i = 3.
Krok 1.:
c s
−s c
−1
3
=
⊞
0
,
16
kde c = − 1√
10
, s = 3√
10
.
Krok 2.:
G(2, 3, Θ) =


1 0 0
0 − 1√
10
3√
10
0 − 3√
10
− 1√
10

 .
Pak
G(2, 3, Θ)x =


1 0 0
0 − 1√
10
3√
10
0 − 3√
10
− 1√
10




1
−1
3

 =


1√
10
0

 .
Nulování všech prvků ve vektoru pod pozicí (1, 1)
Mějme n-vektor x. Jestliže požadujeme nulu na všech pozicích v x s výjimkou první, musíme
sestrojit
G(1, 2, Θ),
x(1)
= G(1, 2, Θ)x,
x(2)
= G(1, 3, Θ)x(1)
,
x(3)
= G(1, 4, Θ)x(2)
,
...
x(n)
= G(1, n, Θ)x(n−1)
.
Pak pro matici
P : = G(1, n, Θ) · · · G(1, 3, Θ)G(1, 2, Θ)
platí, že Px je násobek vektoru e1. Protože každá rotace je orotgonální, je ortogonální i matice P.
Příklad 7. Buď
x =


1
−1
2

 .
Nyní
G(1, 2, Θ) =


1√
2
− 1√
2
0
1√
2
1√
2
0
0 0 1

 .
x(1)
= G(1, 2, Θ)x =


√
2
0
2

 ,
G(1, 3, Θ) =




2
6
0 2√
6
0 1 0
− 2√
6
0 2
6



 .
17
Dále
G(1, 3, Θ)x(1)
=


√
6
0
0

 .
Tedy
P = G(1, 3, Θ)G(1, 2, Θ) =



1√
6
− 1√
6
2√
6
1√
2
1√
2
0
− 2
6
2
6
2
6



a
Px =


√
6
0
0

 .
Vytváření nul v maticích
Myšlenku nulování určitého prvku ve vektoru lze triviálně převést na vytváření nul na specifických
pozicích matice.
Pokud chceme vytvořit nulu na pozici (j, i), j > i, můžeme to udělat tak, že zkonstruujeme
Givensovu matici G(i, j, Θ) působící pouze na i-tý a j-tý řádek tak, že G(i, j, Θ)A bude mít
nulu na pozici (j, i).
Příklad 8. V matici
A =


1 2 3
3 3 4
4 5 6


vytvořme nulu na pozici (2, 1) s pomocí matice G(1, 2, Θ).
Krok 1.: Najdeme c a s tak, aby platilo
c s
−s c
aii
aji
=
c s
−s c
1
3
=
⊞
0
.
Takže
c =
1
√
10
, s =
3
√
10
.
Krok 2.: Sestavení G(1, 2, Θ). Je
G(1, 2, Θ) =


1√
10
3√
10
0
− 3√
10
1√
10
0
0 0 1

 ,
pak
G(1, 2, Θ)A =


10√
10
11√
10
15√
10
0 − 3√
10
− 5√
10
4 5 6

 .
18
4.5 QR-rozklad pomocí Givensovy matice
Opět chceme setrojit matice Q1, Q2, . . ., Qs tentokrát však pomocí Givensových matic tak, aby
A(1)
= Q1A měla nuly pod prvkem (1, 1) v prvním sloupci, matice A(2)
= Q2A(1)
měla nuly pod
(2, 2) ve druhém sloupci, atd. Každou z matic Qi lze sestrojit jako součin Givensových matic –
ten je možné sestrojit takto:
Q1 : = G(1, m, Θ)G(1, m − 1, Θ) · · ·G(1, 3, Θ)G(1, 2, Θ)
Q2 : = G(2, m, Θ)G(2, m − 1, Θ) · · ·G(2, 3, Θ)
...
Buď s = min{m − 1, n}. Pak
R = A(s)
= QsA(s−1)
= · · · = QsQs−1 · · · Q2Q1A = QT
A.
Nyní máme A = QR, kde QT
= Qs · · · Q2Q1. To lze zformulovat do následující věty.
Věta 8. Givensův QR-rozklad. Buď A matice m × n a nechť s = min{m − 1, n}. Existuje
s ortogonálních matic Q1, . . ., Qs definovaných jako
Qi : = G(i, m, Θ)G(i, m − 1, Θ) · · · G(i, i + 1, Θ),
že pro
Q = QT
1 QT
2 · · · QT
s
platí
A = QR,
kde R je matice m × n s nulami pod hlavní diagonálou.
Znázorněme si schématicky Givensovu metodu redukce matice A ∈ R3×3
na horní trojúhelníkový
tvar (symbol • značí prvky, které se transformací nezměnily, a ± značí prvky, které se změnily):
A =


• • •
• • •
• • •

 G(1, 2, Θ)
−−−−−−→


± ± ±
0 ± ±
• • •

 G(1, 3, Θ)
−−−−−−→
G(1, 3, Θ)
−−−−−−→


± ± ±
0 • •
0 ± ±

 G(2, 3, Θ)
−−−−−−→


• • •
0 ± ±
0 0 ±

 = R.
Příklad 9. Nechť
A =


0 1 1
1 2 3
1 1 1

 .
Krok 1.: Najděme c a s tak, aby
c s
−s c
a11
a21
=
⊞
0
.
19
Neboť a11 = 0 a a21 = 1, musí být c = 0 a s = 1, tedy
G(1, 2, Θ) =


0 1 0
−1 0 0
0 0 1

 .
Pak dostneme
A = G(1, 2, Θ)A =


0 1 0
−1 0 0
0 0 1




0 1 1
1 2 3
1 1 1

 =


1 2 3
0 −1 −1
1 1 1

 .
Nyní najděme c a s tak, aby
c s
−s c
a11
a31
=
⊞
0
.
Neboť a11 = 1 a a31 = 1, bude c = 1√
2
a s = 1√
2
, tedy
G(1, 3, Θ) =


1√
2
0 1√
2
0 1 0
− 1√
2
0 1√
2

 .
Celkem
A(1)
= G(1, 3, Θ)A =


1√
2
0 1√
2
0 1 0
− 1√
2
0 1√
2




1 2 3
0 −1 −1
1 1 1

 =



√
2 3√
2
2
√
2
0 −1 −1
0 − 1√
2
−
√
2


 .
Krok 2.: Určeme c a s tak, aby
c s
−s c
a
(1)
22
a
(1)
32
=
⊞
0
.
Neboť a
(1)
22 = −1 a a
(1)
32 = − 1√
2
, bude c = − 2
3
a s = − 1√
3
, tedy
G(1, 3, Θ)A(1)
=




1 0 0
0 − 2
3
− 1√
3
0 − 1√
3
− 2
3







√
2 3√
2
2
√
2
0 −1 −1
0 − 1√
2
−
√
2


 =



√
2 3√
2
2
√
2
0 3
2
2 2
3
0 0 1√
3


 = R.
Což je tentýž výsledek jako v Příkladě 3.
4.6 Srovnání algoritmů
Při výpočtu QR-rozkladu pomocí Householderovy matice je počet provedených operací roven číslu
n2
(3 −
n
3
).
K explicitnímu vyjádření matice Q je navíc potřeba
2(m2
n − mn2
+
1
3
n3
)
20
operací, tedy celkem
2m2
n − mn2
+
1
3
n3
.
Zatímco pro QR-rozklad pomocí Givensovy matice je tento počet dvojnásobný, tj.
2n2
(3 −
n
3
).
Ovšem pokud v metodě s Givensovou maticí nahradíme matici rotace ( c s
−s c ) a matice odrazu
( c s
s −c ) maticemi 1 a
−a 1 a a 1
1 −a s ortogonálními sloupci, pak se nám podaří snížit počet operací
na úroveň metody využívající Householderovy matice – jedná se o tzv. matice rychlé Givensovy
transformace.
Householderova matice má však tu nevýhodu, že v matici, kterou ji násobíme, nám změní všechny
prvky (zatímco Givensova matice jen i-tý a k-tý řádek), takže nám například může z řídké matice
vytvořit matici plnou.
V modifikované metodě s Gram-Schmidtovým algoritmem je počet operací
mn2
.
21
5. Aplikace QR-rozkladu
5.1 Řešení systémů rovnic pomocí QR-rozkladu
Při řešení soustavy Ax = b pomocí Gaussovy eliminace (tj. pomocí LU-rozkladu) je nutné provést
mn2
2
− n3
6
operací. Tedy méně než při QR-rozkladu. Ovšem na rozdíl od LU-rozkladu existuje
QR-rozklad vždy, je jednoznačný, numericky stabilnější a nezvyšuje podmíněnost soustavy.
Mějme tedy soustavu
Ax = b
s regulární n × n maticí A a její rozklad, tj. nechť A = QR. Pak můžeme soustavu přepsat
QR = b,
tedy
Rx = QT
b = b
′
,
kde R je horní trojúhelníková matice, a její řešení lze nalézt pomocí zpětné substituce.
Příklad 10. Nechť
A =


0 1 1
1 2 3
1 1 1

 , b =


2
6
3

 .
Z Příkladu 3 víme, že
R =


−1, 4142 −2, 1213 −2, 8284
0 1, 2247 1, 6330
0 0 −0, 5770

 ,
H1 =



0 − 1√
2
− 1√
2
− 1√
2
1
2
−1
2
− 1√
2
−1
2
1
2


 =


0 −0, 7071 −0, 7071
−0, 7071 0, 5 −0, 5
−0, 7071 −0, 5 0, 5


a
H2 =


1 0 0
0 −0, 1691 −0, 9856
0 −0, 9856 0, 1691

 .
Vypočítejme b
′
. Dostáváme
y1 = b =


2
6
3

 ,
y2 = H1y1 =


−6, 3640
0, 0858
−2, 9145

 ,
a tedy
b
′
= y3 = H2y2 =


−6, 3640
2, 8577
−0, 5774

 .
22
Vektor b
′
lze spočítat bez explicitního vyjádření matice Q – proto je výhodnější použití QR-rozkladu
pomocí Householderových matic.
Vyřešme
Rx = b
′
.


−1, 4142 −2, 1213 −2, 8284
0 1, 2247 1, 6330
0 0 −0, 5770

 =


−6, 3640
2, 8577
−0, 5774

 ,
což dává řešení
x =


1
1
1

 .
5.2 QR-rozklad a Moore-Penroseova zobecněná inverse
Definice 4. Matici A†
∈ Rn×m
nazveme pseudoinversní maticí k matici A ∈ Rm×n
, jestliže pro
ni platí
• A†
AA†
= A†
,
• AA†
A = A,
• (A†
A)T
= A†
A,
• (AA†
)T
= AA†
.
Přičemž lze dokázat, že matice A†
je určena jednoznačně.
Věta 9. Má-li matice A lineárně nezávislé sloupce, je matice AT
A regulární a platí
A†
= (AT
A)−1
AT
.
Důkaz: I. Předpokládejme, že AT
A není regulární. Pak existuje nenulový vektor y tak, že AT
Ay =
= 0. Potom však
0 = yT
AT
Ay =< Ay, Ay >,
tedy Ay = 0, což je spor s tím, že sloupce matice A jsou lineárně nezávislé.
II. Stačí ověřit, že matice (AT
A)−1
AT
splňuje podmínky definice pseudoinversní matice.
Pro matici s lineárně nezávislými sloupci lze použít QR-rozklad i na výpočet pseudoinversní matice.
Věta 10. Nechť matice A má lineárně nezávislé sloupce a nechť A = QR je její QR-rozklad.
Potom pro pseudoinversní matici A†
platí
A†
= R−1
QT
.
Důkaz: Nechť A = QR. Podle předchozí věty platí, že A†
= (AT
A)−1
AT
, takže
A†
= (RT
QT
QR)−1
RT
QT
= R−1
(RT
)−1
RT
QT
= R−1
QT
.
23
5.3 QR-rozklad a vlastní čísla matice A – QR-algoritmus
Základní QR-algoritmus: Mějme matici A. Sestrojme její QR-rozklad, tj.
A = A0 = Q0R0,
určeme
A1 : = R0Q0.
Nyní sestrojme QR-rozklad matice A1, tj.
A1 = Q1R1,
a spočtěme
A2 = R1Q1.
Takto pokračujme analogicky dále.
Jistě platí
Ak+1 = RkQk = QT
k AkQk = QT
k Rk−1Qk−1Qk =
= QT
k QT
k−1Ak−1Qk−1Qk = · · · = (Q0Q1 · · · Qk)T
A(Q0Q1 · · · Qk).
Tedy matice A0, A1, . . . jsou kongruentní (tj. A ≡ B ⇐⇒ A = PT
BP). Navíc díky ortogonálnosti
matic Q0, Q1, . . . jsou matice A0, A1, . . . také podobné (tj. A ∼ B (též A ≈ B)
⇐⇒ A = P−1
BP). Tyto matice mají díky podobnosti stejná vlastní čísla jako matice A. Posloupnost
těchto matic konverguje za určitých předpokladů k horní trojúhelníkové (resp. horní
blokově trojúhelníkové) matici, která má vlastní čísla na diagonále (resp. diagonální bloky mají
vlastní čísla se stejnou absolutní hodnotou) seřazena podle velikosti počínaje největším vlastním
čísel. „Poddiagonální prvky (resp. „poddiagonální bloky ) konvergují k nule. Ovšem důkazy konvergence
existují jen pro některé speciální typy matic. Například má-li matice A kladná vlastní
čísla, pak Qk konverguje k jednotkové matici a posloupnost matic Ak k horní troúhelníkové matici,
přičemž diagonální prvky této matice jsou vlastní čísla matice A.
Příklad 11. Určeme vlastní čísla matice
A =


2 1
3
1
3 −5
3
1
0 11
9
5
3

 .
Lze snadno ověřit, že vlastní čísla matice A jsou λ1 = 1, λ2 = −2 a λ3 = 3. Výsledky získané
QR-algoritmem:
k a
(k)
11 a
(k)
22 a
(k)
33
1 2,0 -1,6666667 1,6666667
5 3,1781374 -2,2260322 1,0478949
10 2,9486278 -1,9471270 0,9984996
15 3,0003596 -2,0064061 1,0000468
20 2,9991547 -1,9991527 0,9999984
25 3,0001104 -2,0001098 0,9999999
24
Ke zrychlení konvergence lze využít tzv. posunutí a počítat nikoli rozklad matice Ak = QkRk,
nýbrž matice
Ak − σkE = QkRk.
Původní spektrum matice A se tímto posune o σk (je výhodné volit jej jako nějakou aproximaci
vlastního čísla; matice Ak − σkE má vlastní čísla λj − σk, jsou-li λj vlastní čísla matice Ak).
Zaznamenáváme-li velikosti posunití σk, snadno ze znalosti spektra matice Ak najdeme spektrum
matice A. Protože ještě (v metodě bez posunutí)
Ak = (Q0Q1 · · · Qk−1)T
A(Q0Q1 · · · Qk−1),
je
A = (Q0Q1 · · · Qk−1)T
Ak(Q0Q1 · · · Qk−1).
Vlastní vektory y matice A dostaneme z vlastních vektorů matice Ak podle vzorce
y = Q0Q1 · · ·Qk−1z.
Týž vzorec (s maticemi Qj) zůstává v platnosti i pro metodu s posunutími, protože při posunutí
se vlastní vektory nemění.
Volí-li se posunutí speciálně, dostáváme v některých důležitých případech i kubickou konvergenci
(tzn. zhruba řečeno, počet platných míst se v každém kroku přibližně ztrojnásobí).
Poznámka 3. Nejvýhodnější se jeví upravit nejdříve matici A do tzv. Hessenbergova tvaru (tj.
aij = 0 pro j < i−1, i, j = 1, . . . , n) pomocí Gaussovy eliminace a pak na tuto upravenou matici
použít QR-rozklad. Konvergence je potom rychlejší (obzvláště použijeme-li metodu posunu, kde
za σk volíme tzv. Rayleighův podíl
eT
k H ek
eT
k ek
, kde H je právě matice A v Hessenbergově tvaru). Navíc
platí, že je-li matice A v Hessenbergově tvaru, pak každá z matic Hk je také v Hessenbergově
tvaru, a to i při metodě posunutí.
25
Seznam použité literatury
[1] Miroslav Fiedler: Numerické metody lineární algebry. SNTL Praha, 1978. Str. 200-213,
242-249.
[2] Biswa Nath Datta: Numerical linear algebra and applications. Brooks and Cole Publishing
Company, Pacific Grove, California, 1995. Str. 133-183, 215-217, 236, 237.
[3] Anthony Rahlston: Základy numerické matematiky. Academia Praha, 1978. Str. 551-564,
576-577.
[4] Kubíček Milan, Dubcová Miroslava, Janovská Drahoslava: Numerické metody a algoritmy.
VŠCHT Praha, 2005. Str. 167-170.
26