109
Ortogonální transformace a QR rozklady Petr Tichý 9. října 2013 1
Ortogonální transformace a QR rozklady
Petr Tichý
9. října 2013
1
ÚvodUnitární (ortogonální) transformace, Gram-Schmidtova ortogonalizace
Příklad Schurovy věty → unitární transformace nezvětšujíchyby ve vstupních datech.
Tato kapitola: základními dva typy unitárních transformací →Givensovy rotace, Householderovy reflexe.
Využití: (numericky stabilní) transformace matice na maticis předem zvolenou strukturou.
Výpočet QR rozkladu → široké použití.
Gram-Schmidtova ortogonalizace, výpočetní náročnost,numerická stabilita.
2
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
3
Ortogonální projekce
Libovolný x ∈ Rn lze rozložit na součet ortogonálních vektorů
x = xq + y tak, že xq leží v předem daném směru q a 〈y, xq〉 = 0.
ei
ej
x
y
xq
q
4
Ortogonální projekce vektoru x do směru q
Dán q, ‖q‖ = 1. Hledáme rozklad vektoru x ∈ Rn,
x = xq + y
takový, že xq leží v předem daném směru q a y ⊥ q.
5
Ortogonální projekce vektoru x do směru q
Dán q, ‖q‖ = 1. Hledáme rozklad vektoru x ∈ Rn,
x = xq + y
takový, že xq leží v předem daném směru q a y ⊥ q.
Ukážeme, žexq = 〈x, q〉 q . �
xq nazveme ortogonální projekcí x do prostoru span{q}.
5
Ortogonální projekce vektoru x do směru q
Dán q, ‖q‖ = 1. Hledáme rozklad vektoru x ∈ Rn,
x = xq + y
takový, že xq leží v předem daném směru q a y ⊥ q.
Ukážeme, žexq = 〈x, q〉 q . �
xq nazveme ortogonální projekcí x do prostoru span{q}.
Ekvivalentní zápis,
xq = Pq x , kde Pq ≡ qqT . �
5
Ortogonální projekce vektoru x do směru q
Dán q, ‖q‖ = 1. Hledáme rozklad vektoru x ∈ Rn,
x = xq + y
takový, že xq leží v předem daném směru q a y ⊥ q.
Ukážeme, žexq = 〈x, q〉 q . �
xq nazveme ortogonální projekcí x do prostoru span{q}.
Ekvivalentní zápis,
xq = Pq x , kde Pq ≡ qqT . �
Pq je čtvercová symetrická matice řádu n,
Pq : x → Pqx ∈ span{q}, x − Pqx ⊥ span{q} .
Říkáme: Pq projektuje x do prostoru generovaného bází {q}.
5
Ortogonální projekce vektoru do podprostoru
Uvažujme ortonormální soubor vektorů {q1, . . . , qm} ⊂ Rn,
Q ≡ span{q1, . . . , qm}, Q ≡ [q1, . . . , qm] ∈ Rn×m.
6
Ortogonální projekce vektoru do podprostoru
Uvažujme ortonormální soubor vektorů {q1, . . . , qm} ⊂ Rn,
Q ≡ span{q1, . . . , qm}, Q ≡ [q1, . . . , qm] ∈ Rn×m.
Dán x ∈ Rn, hledáme rozklad x na složky
x = xQ + y, xQ ∈ Q, y ⊥ Q.
6
Ortogonální projekce vektoru do podprostoru
Uvažujme ortonormální soubor vektorů {q1, . . . , qm} ⊂ Rn,
Q ≡ span{q1, . . . , qm}, Q ≡ [q1, . . . , qm] ∈ Rn×m.
Dán x ∈ Rn, hledáme rozklad x na složky
x = xQ + y, xQ ∈ Q, y ⊥ Q.
xQ nazveme ortogonální projekcí x do Q. Platí
xQ = PQ x , kde PQ ≡ QQT . �
6
Ortogonální projekce vektoru do podprostoru
Uvažujme ortonormální soubor vektorů {q1, . . . , qm} ⊂ Rn,
Q ≡ span{q1, . . . , qm}, Q ≡ [q1, . . . , qm] ∈ Rn×m.
Dán x ∈ Rn, hledáme rozklad x na složky
x = xQ + y, xQ ∈ Q, y ⊥ Q.
xQ nazveme ortogonální projekcí x do Q. Platí
xQ = PQ x , kde PQ ≡ QQT . �
PQ je čtvercová symetrická matice hodnosti m, projektujelibovolný x do prostoru Q.
6
Ortogonální projektory
Definice: Ortogonální projektor je lineární operátor P , kterýje symetrický a idempotentní (P T = P a P 2 = P ).
7
Ortogonální projektory
Definice: Ortogonální projektor je lineární operátor P , kterýje symetrický a idempotentní (P T = P a P 2 = P ).
Matice PQ je ortogonálním projektorem. �
7
Ortogonální projektory
Definice: Ortogonální projektor je lineární operátor P , kterýje symetrický a idempotentní (P T = P a P 2 = P ).
Matice PQ je ortogonálním projektorem. �
Vektor y = x − xQ lze psát ve tvaru
x − xQ = ΠQ x , kde ΠQ ≡ I − QQT = I − PQ.
7
Ortogonální projektory
Definice: Ortogonální projektor je lineární operátor P , kterýje symetrický a idempotentní (P T = P a P 2 = P ).
Matice PQ je ortogonálním projektorem. �
Vektor y = x − xQ lze psát ve tvaru
x − xQ = ΠQ x , kde ΠQ ≡ I − QQT = I − PQ.
ΠQ je ortogonálním projektorem do prostoru Q⊥.ΠQ se nazývá projektor komplementární k projektoru PQ.
7
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
8
Givensovy rotace v Rn
Úloha v R2
Chceme sestrojit matici G(ϕ) ∈ R2×2, která realizuje pootočenílibovolného vektoru x o úhel ϕ proti směru hodinových ručiček.
Zapíšeme-li vektor x v bázi {e1, e2},
x =
[ξ1
ξ2
]= ξ1
[10
]+ ξ2
[01
]= ξ1 e1 + ξ2 e2 ,
je možné vektor G(ϕ)x vyjádřit ve tvaru
G(ϕ)x = ξ1 (G(ϕ) e1) + ξ2 (G(ϕ) e2).
9
Givensovy rotace v Rn
Úloha v R2
Chceme sestrojit matici G(ϕ) ∈ R2×2, která realizuje pootočenílibovolného vektoru x o úhel ϕ proti směru hodinových ručiček.
Zapíšeme-li vektor x v bázi {e1, e2},
x =
[ξ1
ξ2
]= ξ1
[10
]+ ξ2
[01
]= ξ1 e1 + ξ2 e2 ,
je možné vektor G(ϕ)x vyjádřit ve tvaru
G(ϕ)x = ξ1 (G(ϕ) e1) + ξ2 (G(ϕ) e2).
Otáčí-li G(ϕ) bázové vektory e1 a e2 o úhel ϕ, otáčí ilibovolný vektor x o úhel ϕ.
9
Rotace jednotkových vektorů o úhel ϕ
0
10
Matice Givensovy rotace v R2
G(ϕ)
[10
]=
[cos ϕsin ϕ
], G(ϕ)
[01
]=
[− sin ϕ
cos ϕ
],
tj.
G(ϕ) =
[cos ϕ − sin ϕsin ϕ cos ϕ
]
Matice G(ϕ) se nazývá matice Givensovy rotace.
11
Matice Givensovy rotace v R2
G(ϕ)
[10
]=
[cos ϕsin ϕ
], G(ϕ)
[01
]=
[− sin ϕ
cos ϕ
],
tj.
G(ϕ) =
[cos ϕ − sin ϕsin ϕ cos ϕ
]
Matice G(ϕ) se nazývá matice Givensovy rotace.
Použití: Jsou-li v R2 dány vektory x a y, ‖x‖ = ‖y‖ 6= 0.Potom y lze získat pootočením x,
y = G(ϕ)x.
11
Matice Givensovy rotace v Rn
Elementární Givensova rotace
Rotaci v rovině dané dvojicí jednotkových vektorů{ei , ej} , i < j, o úhel ϕ ve směru od ei k ej lze realizovat pomocí
Gi,j(ϕ) =
1. . .
1cos ϕ − sin ϕ
1. . .
1sin ϕ cos ϕ
1. . .
1
.
12
Matice Givensovy rotace v Rn
Vlastnosti
je ortonormální, det(Gi,j(ϕ)) = 1 (cvičení).
Gi,j(ϕ)x modifikuje pouze i-tý a j-tý prvek vektorux = [ξ1, . . . , ξn]T ,
Gi,j(ϕ) x =
ξ1...
ξi cos ϕ − ξj sin ϕ...
ξi sin ϕ + ξj cos ϕ...
ξn
· · · i-tý řádek
· · · j-tý řádek.
�
13
Nulování prvků vektoru pomocí Givensových rotacíNulování v R2
Požadujeme, aby platilo
y = G(ϕ) x =
[cos ϕ − sin ϕsin ϕ cos ϕ
] [ξ1
ξ2
]=
[±
√ξ2
1 + ξ22
0
].
14
Nulování prvků vektoru pomocí Givensových rotacíNulování v R2
Požadujeme, aby platilo
y = G(ϕ) x =
[cos ϕ − sin ϕsin ϕ cos ϕ
] [ξ1
ξ2
]=
[±
√ξ2
1 + ξ22
0
].
Ze vztahuξ1 sin ϕ + ξ2 cos ϕ = 0
dostaneme
sin ϕ = ∓ ξ2√ξ2
1 + ξ22
,
cos ϕ = ± ξ1√ξ2
1 + ξ22
.
Ověření �14
Nulování prvků vektoru pomocí Givensových rotacíNulování v Rn
Úloha: vynulovat n − 1 složek vektoru x ∈ Rn,
x −→ ±‖x‖e1 .
15
Nulování prvků vektoru pomocí Givensových rotacíNulování v Rn
Úloha: vynulovat n − 1 složek vektoru x ∈ Rn,
x −→ ±‖x‖e1 .
Opakovaně aplikujeme elementární Givensovy rotace:
x =
•∗...∗•
→
•∗...•0
→ . . . →
••0...0
→
±‖x‖0...
0
= y .
Nulujeme prvky na pozicích n, n − 1, . . . , 2 (volíme např.roviny rotace span{e1 , en}, . . . , span{e1 , e2})
15
Nulování prvků vektoru pomocí Givensových rotacíFormalizace
Označíme-li jednotlivé elementární Givensovy rotace jakoG1,2, . . . , G1,n−1, G1,n, potom
y = Γx , kde Γ ≡ G1,2 . . . G1,n−1G1,n .
Matici Γ budeme nazývat složenou Givensovou rotací.
16
Nulování prvků vektoru pomocí Givensových rotacíFormalizace
Označíme-li jednotlivé elementární Givensovy rotace jakoG1,2, . . . , G1,n−1, G1,n, potom
y = Γx , kde Γ ≡ G1,2 . . . G1,n−1G1,n .
Matici Γ budeme nazývat složenou Givensovou rotací.
Prvky nemusíme nutně nulovat v pořadí, které jsme naznačili.
16
Nulování prvků vektoru pomocí Givensových rotacíFormalizace
Označíme-li jednotlivé elementární Givensovy rotace jakoG1,2, . . . , G1,n−1, G1,n, potom
y = Γx , kde Γ ≡ G1,2 . . . G1,n−1G1,n .
Matici Γ budeme nazývat složenou Givensovou rotací.
Prvky nemusíme nutně nulovat v pořadí, které jsme naznačili.
Násobení elementárních (tedy i složených) Givensových rotacínení obecně komutativní (cvičení).
16
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
17
Nulování prvků vektoru pomocí Householderových reflexíNulování v Rn
Druhá základní unitární transformace je Householderovareflexe (zrcadlení, odraz).
18
Nulování prvků vektoru pomocí Householderových reflexíNulování v Rn
Druhá základní unitární transformace je Householderovareflexe (zrcadlení, odraz).
Úloha: Nechť je v Rn dána nadrovina dimenze n − 1, kterou
popíšeme jejím normálovým vektorem q, ‖q‖ = 1,
H(q) ≡ {z ∈ Rn : z ⊥ q} ,
a nechť je dán vektor x ∈ Rn.
18
Nulování prvků vektoru pomocí Householderových reflexíNulování v Rn
Druhá základní unitární transformace je Householderovareflexe (zrcadlení, odraz).
Úloha: Nechť je v Rn dána nadrovina dimenze n − 1, kterou
popíšeme jejím normálovým vektorem q, ‖q‖ = 1,
H(q) ≡ {z ∈ Rn : z ⊥ q} ,
a nechť je dán vektor x ∈ Rn.
Cíl: nalézt zrcadlový obraz vektoru x podle nadroviny H(q)(nadrovinu H(q) nazveme nadrovinou zrcadlení).
18
Zrcadlení
x
x xq−
y = x 2xq−
xq
q
0
xq = (qqT ) x,
y = (x − xq) − xq = x − 2xq = (I − 2qqT ) x .
19
Matice Householderovy reflexe
Nechť q ∈ Rn a ‖q‖ = 1. Pak matici
H(q) = I − 2qqT ∈ Rn×n
nazýváme maticí Householderovy reflexe vzhledemk nadrovině H(q) definované normálovým vektorem q.
20
Matice Householderovy reflexe
Nechť q ∈ Rn a ‖q‖ = 1. Pak matici
H(q) = I − 2qqT ∈ Rn×n
nazýváme maticí Householderovy reflexe vzhledemk nadrovině H(q) definované normálovým vektorem q.
Vlastnosti: H(q) je ortonormální a symetrická, platíH2(q) = I, det(H(q)) = −1 (cvičení).
20
Zrcadlení x na y
pomocí Householderovy reflexe
Úloha: Dány x a y stejné délky, nalézt H tak, aby y = Hx.
21
Zrcadlení x na y
pomocí Householderovy reflexe
Úloha: Dány x a y stejné délky, nalézt H tak, aby y = Hx.
Zrcadlení x na ±y v Rn
Nechť jsou dány dva různé vektory x ∈ Rn a y ∈ Rn, ‖x‖ = ‖y‖,a nechť
q1 ≡ x − y
‖x − y‖ , q2 ≡ x + y
‖x + y‖ .
PotomH(q1)x = y, H(q2)x = −y.
21
Zrcadlení x na y
pomocí Householderovy reflexe
Úloha: Dány x a y stejné délky, nalézt H tak, aby y = Hx.
Zrcadlení x na ±y v Rn
Nechť jsou dány dva různé vektory x ∈ Rn a y ∈ Rn, ‖x‖ = ‖y‖,a nechť
q1 ≡ x − y
‖x − y‖ , q2 ≡ x + y
‖x + y‖ .
PotomH(q1)x = y, H(q2)x = −y.
Důkaz: Vektor x − y je kolmý k nadrovině zrcadlení vektoru x navektor y. Podobně, vektor x + y je kolmý k nadrovině zrcadlenívektoru x na vektor −y.
21
Nulování prvků vektoru pomocí Householderových reflexí
Dán x ∈ Rn a
y = ±‖x‖e1.
Hledáme H.
22
Nulování prvků vektoru pomocí Householderových reflexí
Dán x ∈ Rn a
y = ±‖x‖e1.
Hledáme H.
Dle předchozího, pro
q1 =x − ‖x‖e1
‖x − ‖x‖e1‖ , q2 =x + ‖x‖e1
‖x + ‖x‖e1‖
jeH(q1)x = ‖x‖e1, H(q2)x = −‖x‖e1.
22
Nulování prvků vektoru pomocí Householderových reflexí
Dán x ∈ Rn a
y = ±‖x‖e1.
Hledáme H.
Dle předchozího, pro
q1 =x − ‖x‖e1
‖x − ‖x‖e1‖ , q2 =x + ‖x‖e1
‖x + ‖x‖e1‖
jeH(q1)x = ‖x‖e1, H(q2)x = −‖x‖e1.
Je lepší zvolit q1 nebo q2?
22
Nulování prvků vektoru pomocí Householderových reflexíVolba vektoru q
x ± ‖x‖e1 =
ξ1 ± ‖x‖ξ2...
ξn
Pokud ξ1 ≥ 0, ξ1 ≈ ‖x‖, potom při odečítání ξ1 − ‖x‖ můžedojít k vyrušení platných číslic. Navíc, ‖ξ1 − ‖x‖e1‖ ≈ 0,dělíme malou normou → relativní zvětšení chybyvypočteného vektoru ve srovnání s chybou obsaženou v x.
23
Nulování prvků vektoru pomocí Householderových reflexíVolba vektoru q
x ± ‖x‖e1 =
ξ1 ± ‖x‖ξ2...
ξn
Pokud ξ1 ≥ 0, ξ1 ≈ ‖x‖, potom při odečítání ξ1 − ‖x‖ můžedojít k vyrušení platných číslic. Navíc, ‖ξ1 − ‖x‖e1‖ ≈ 0,dělíme malou normou → relativní zvětšení chybyvypočteného vektoru ve srovnání s chybou obsaženou v x.
Kvůli numerické stabilitě, pokud ξ1 ≥ 0 pak
q ≡ x + ‖x‖e1
‖x + ‖x‖e1‖ .
23
Nulování prvků vektoru pomocí Householderových reflexíVolba vektoru q
x ± ‖x‖e1 =
ξ1 ± ‖x‖ξ2...
ξn
Pokud ξ1 ≥ 0, ξ1 ≈ ‖x‖, potom při odečítání ξ1 − ‖x‖ můžedojít k vyrušení platných číslic. Navíc, ‖ξ1 − ‖x‖e1‖ ≈ 0,dělíme malou normou → relativní zvětšení chybyvypočteného vektoru ve srovnání s chybou obsaženou v x.
Kvůli numerické stabilitě, pokud ξ1 ≥ 0 pak
q ≡ x + ‖x‖e1
‖x + ‖x‖e1‖ .
Analogicky pro ξ1 < 0 volíme q = x−‖x‖e1
‖x−‖x‖e1‖ .
23
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
24
QR rozklad
Nechť A ∈ Cn×m je obecně obdélníková matice. Rozklad
A = QR,
kde Q je matice s ortonormálními sloupci (Q∗Q = I) a R mávšechny prvky pod hlavní diagonálou nulové, nazvemeQR rozkladem matice A.
Pro n > m může QR rozklad vypadat schematicky jakoA
=
Q R
❅❅
❅
0,
kde Q ∈ Cn×n a R ∈ C
n×m.
25
QR rozklad
Vynecháním posledních n − m nulových řádků matice R apříslušných sloupců Q dostaneme
A
=
Q R
❅❅
❅
0,
kde Q ∈ Cn×m a R ∈ C
m×m. Tomuto QR rozkladu, jehožuložení vyžaduje méně paměťového místa než v prvnímpřípadě, se říká ekonomický QR rozklad.
26
QR rozklad
Vynecháním posledních n − m nulových řádků matice R apříslušných sloupců Q dostaneme
A
=
Q R
❅❅
❅
0,
kde Q ∈ Cn×m a R ∈ C
m×m. Tomuto QR rozkladu, jehožuložení vyžaduje méně paměťového místa než v prvnímpřípadě, se říká ekonomický QR rozklad.
Pro n ≤ m má QR rozklad, schematicky, tvarA
=
Q R
❅❅
❅
0,
kde Q ∈ Cn×n a R ∈ C
n×m.26
Obecně o QR rozkladu
Lze použít k numerickému řešení mnoha problémů, napříkladproblémů nejmenších čtverců Ax ≈ b nebo Ax = b.
Ax = b, kde A ∈ Cn×n je regulární matice,
Ax = b ⇐⇒ QRx = b,
⇐⇒ QRx = QQ∗b,
⇐⇒ Rx = Q∗b.
Díky ortogonalitě bývá výpočet numericky stabilnější než připoužití jiných druhů rozkladů (LU rozklad).
Problém pro velké řídké matice
požadavek ortogonality × ztráta řídkosti .
Algoritmů pro výpočet QR rozkladu je několik.
27
QR rozklad užitím Givensových rotací I
A ∈ Cn×m, A = [a1, . . . , am].
Cíl: vynulování všech prvků pod hlavní diagonálou za použitíGivensových rotací, získáme Givensův QR rozklad.
28
QR rozklad užitím Givensových rotací I
A ∈ Cn×m, A = [a1, . . . , am].
Cíl: vynulování všech prvků pod hlavní diagonálou za použitíGivensových rotací, získáme Givensův QR rozklad.
Symbolem Γ1 označíme složenou Givensovu rotaci, kterárealizuje transformaci
a1 =
a1,1
a2,1...
an,1
−→
r1,1
0...0
= Γ1a1 ≡ r1 .
28
QR rozklad užitím Givensových rotací II
Aplikujeme-li Γ1 na A, dostaneme (r1,1 = 0 ⇔ a1 = 0)
A =
• ∗ ∗ ∗• ∗ ∗ ∗• ∗ ∗ ∗• ∗ ∗ ∗• ∗ ∗ ∗
−→
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗
= Γ1A ≡ A(1) ,
29
QR rozklad užitím Givensových rotací II
Aplikujeme-li Γ1 na A, dostaneme (r1,1 = 0 ⇔ a1 = 0)
A =
• ∗ ∗ ∗• ∗ ∗ ∗• ∗ ∗ ∗• ∗ ∗ ∗• ∗ ∗ ∗
−→
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗
= Γ1A ≡ A(1) ,
Nyní nulujeme poddiagonální prvky druhého sloupce maticeA(1) = [r1, a
(1)2 , . . . , a
(1)m ], aplikujeme n − 2 rotací,
a(1)2 =
a(1)1,2
a(1)2,2
a(1)3,2...
a(1)n,2
−→
a(1)1,2
r2,2
0......0
= Γ2a(1)2 ≡ r2 .
29
QR rozklad užitím Givensových rotací III
Aplikací složené rotace Γ2 na matici A(1) = Γ1A dostaneme
A(1) =
∗ ∗ ∗ ∗0 • ∗ ∗0 • ∗ ∗0 • ∗ ∗0 • ∗ ∗
−→
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 ∗ ∗0 0 ∗ ∗
= Γ2Γ1A ≡ A(2) .
30
QR rozklad užitím Givensových rotací III
Aplikací složené rotace Γ2 na matici A(1) = Γ1A dostaneme
A(1) =
∗ ∗ ∗ ∗0 • ∗ ∗0 • ∗ ∗0 • ∗ ∗0 • ∗ ∗
−→
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 ∗ ∗0 0 ∗ ∗
= Γ2Γ1A ≡ A(2) .
V k-tém kroku, k < min{m, n}, konstruujeme složenou rotaci
tak, aby nulovala n − k poddiagonálních prvků sloupce a(k−1)k .
Příslušná složená rotace Γk je blokově diagonální se dvěmabloky, kde první blok je jednotková matice řádu k − 1.
30
QR rozklad užitím Givensových rotací III
Aplikací složené rotace Γ2 na matici A(1) = Γ1A dostaneme
A(1) =
∗ ∗ ∗ ∗0 • ∗ ∗0 • ∗ ∗0 • ∗ ∗0 • ∗ ∗
−→
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 ∗ ∗0 0 ∗ ∗
= Γ2Γ1A ≡ A(2) .
V k-tém kroku, k < min{m, n}, konstruujeme složenou rotaci
tak, aby nulovala n − k poddiagonálních prvků sloupce a(k−1)k .
Příslušná složená rotace Γk je blokově diagonální se dvěmabloky, kde první blok je jednotková matice řádu k − 1.
rk,k = 0 ⇔ je-li k-tý sloupec matice A lineární kombinacípředchozích k − 1 sloupců matice A (cvičení).
30
QR rozklad užitím Givensových rotací IVPoslední transformační maticí Γm−1 pro n > m resp. Γn−1 pro n ≤ m
Případ n > m, platí
A(m−2) =
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 0 •0 0 0 •
−→
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 0 ∗0 0 0 0
= A(m−1) ≡ R .
A(m−1) = Γm−1 . . . Γ1A.
31
QR rozklad užitím Givensových rotací IVPoslední transformační maticí Γm−1 pro n > m resp. Γn−1 pro n ≤ m
Případ n > m, platí
A(m−2) =
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 0 •0 0 0 •
−→
∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 0 ∗0 0 0 0
= A(m−1) ≡ R .
A(m−1) = Γm−1 . . . Γ1A.
Γi jsou unitární ⇒ součin je také unitární matice. Označíme-li
Q ≡ Γ∗1Γ∗
2 . . . Γ∗m−1,
platíA = QR .
31
Existence a jednoznačnost QR rozkladuExistence
Konstrukcí jsme získali rozklad typu
A
=
Q R
❅❅
❅
0,
resp.
A
=
Q R
❅❅
❅
0,
kde Q ∈ Cn×na R ∈ C
n×m.
Existence QR rozkladu (i ekonomického) plyne z konstrunkce.32
Existence a jednoznačnost QR rozkladuJednoznačnost
Není obecně jednoznačný: ∀ diagonální D, D∗D = I, platí
A = QR = (QD∗)(DR) = Q̃R̃.
Povolíme-li na diagonále R jen kladné prvky → jednoznačnost.
33
Existence a jednoznačnost QR rozkladuJednoznačnost
Není obecně jednoznačný: ∀ diagonální D, D∗D = I, platí
A = QR = (QD∗)(DR) = Q̃R̃.
Povolíme-li na diagonále R jen kladné prvky → jednoznačnost.
Jednoznačnost QR rozkladu
Nechť A ∈ Cn×m, n ≥ m, je matice s lineárně nezávislými sloupci.Pak existuje jediná dvojice matic Q ∈ Cn×m a R ∈ Cm×m taková,že Q má ortonormální sloupce a R je horní trojúhelníková matices kladnými diagonálními prvky a přitom
A = QR.
�
33
QR rozklad pomocí Householderových reflexí
Householderův QR rozklad – konstrukčně stejný jakoGivensův, pro n > m transformujeme A ∈ Cn×m
A → H1A → . . . → Hm−1 . . . H2H1A ≡ R ,
kde Hk je blokově diagonální matice, první blok je jednotkovámatice řádu k − 1, druhý blok je matice Householderovyreflexe nulující poddiagonální prvky aktuálního k-tého sloupce.
Označme
Q ≡ H∗1 H∗
2 . . . H∗m−1 ⇒ A = QR .
Analogicky QR rozklad matice A ∈ Cn×m pro n ≤ m.
34
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
35
QR rozklad a Gram-Schmidtův ortogonalizační proces
Pro jednoduchost případ n > m a podprostor Cn generovanýlineárně nezávislými vektory {a1, . . . , am}.
Gram-Schmidtův ortogonalizační proces spočteortonormální bázi {q1, . . . , qm} takovou, že
span{a1, . . . , ak} = span{q1, . . . , qk}, k = 1, . . . , m.
36
QR rozklad a Gram-Schmidtův ortogonalizační proces
Pro jednoduchost případ n > m a podprostor Cn generovanýlineárně nezávislými vektory {a1, . . . , am}.
Gram-Schmidtův ortogonalizační proces spočteortonormální bázi {q1, . . . , qm} takovou, že
span{a1, . . . , ak} = span{q1, . . . , qk}, k = 1, . . . , m.
k = 1: q1? span{a1} = span{q1} a ‖q1‖ = 1 ⇒
q1 =a1
‖a1‖ . Označme r1,1 = ‖a1‖.
36
QR rozklad a Gram-Schmidtův ortogonalizační proces
Pro jednoduchost případ n > m a podprostor Cn generovanýlineárně nezávislými vektory {a1, . . . , am}.
Gram-Schmidtův ortogonalizační proces spočteortonormální bázi {q1, . . . , qm} takovou, že
span{a1, . . . , ak} = span{q1, . . . , qk}, k = 1, . . . , m.
k = 1: q1? span{a1} = span{q1} a ‖q1‖ = 1 ⇒
q1 =a1
‖a1‖ . Označme r1,1 = ‖a1‖.
k = 2: q2? span{a1, a2} = span{q1, q2}, q2 ⊥ q1, ‖q2‖ = 1.Odečteme projekce a2 ve směru q1 a normalizujeme,
z = (I − q1q∗1) a2 = a2 − (q∗
1 a2)q1, r2,2 = ‖z‖ ,
q2 = z/r2,2 .
36
Gram-Schmidtův ortogonalizační procesk-tý krok
obecné k: odečteme od ak projekci na prostor generovanýsloupci [q1, . . . , qk−1] = Qk−1 a normalizujeme
z = (I − Qk−1Q∗k−1) ak = ak −
k−1∑
i=1
(q∗i ak)qi
rk,k = ‖z‖qk = z/rk,k .
37
Gram-Schmidtův ortogonalizační procesk-tý krok
obecné k: odečteme od ak projekci na prostor generovanýsloupci [q1, . . . , qk−1] = Qk−1 a normalizujeme
z = (I − Qk−1Q∗k−1) ak = ak −
k−1∑
i=1
(q∗i ak)qi
rk,k = ‖z‖qk = z/rk,k .
Proveditelnost: Jelikož jsou vektory a1, . . . , am lineárněnezávislé, je vektor z vždy nenulový, a tudíž rk,k 6= 0,k = 1, . . . , m.
37
Gram-Schmidtův ortogonalizační procesk-tý krok
obecné k: odečteme od ak projekci na prostor generovanýsloupci [q1, . . . , qk−1] = Qk−1 a normalizujeme
z = (I − Qk−1Q∗k−1) ak = ak −
k−1∑
i=1
(q∗i ak)qi
rk,k = ‖z‖qk = z/rk,k .
Proveditelnost: Jelikož jsou vektory a1, . . . , am lineárněnezávislé, je vektor z vždy nenulový, a tudíž rk,k 6= 0,k = 1, . . . , m.
Dosadíme-li za z vektor rk,k qk a označíme-li ri,k = q∗i ak,
ak =k−1∑
i=1
ri,k qi + rk,k qk, k = 1, . . . , m.
37
Gram-Schmidtův ortogonalizační procesQR rozklad
Po rozepsání
a1 = r1,1q1
a2 = r1,2q1 + r2,2q2
...
am = r1,mq1 + r2,mq2 + · · · + rm,mqm.
38
Gram-Schmidtův ortogonalizační procesQR rozklad
Po rozepsání
a1 = r1,1q1
a2 = r1,2q1 + r2,2q2
...
am = r1,mq1 + r2,mq2 + · · · + rm,mqm.
S označením A = [a1, . . . , am], Q = [q1, . . . , qm] a
R ≡
r1,1 r1,2 · · · r1,m
0 r2,2 r2,m
.... . . . . .
...0 · · · 0 rm,m
∈ Cm×m
dostávámeA = QR
38
Gram-Schmidtův ortogonalizační procesQR rozklad – komentáře
Narozdíl od QR rozkladu spočteného pomocí rotací či reflexí→ rozměr Q je n × m, R má rozměr m × m,→ pro n > m máme přímo ekonomický QR rozklad.
39
Gram-Schmidtův ortogonalizační procesQR rozklad – komentáře
Narozdíl od QR rozkladu spočteného pomocí rotací či reflexí→ rozměr Q je n × m, R má rozměr m × m,→ pro n > m máme přímo ekonomický QR rozklad.
Gram-Schmidtův ortogonalizační proces vzhledemk obecnému skalárnímu součinu 〈·, ·〉:
z = ak −k−1∑
i=1
〈ak, qi〉qi, rk,k =√
〈z, z〉,
qk = z/rk,k.
{q1, . . . , qk} je ortonormální báze prostoru span{a1, . . . , ak}.
39
Implementace Gram-Schmidtova procesuKlasický Gram-Schmidtův algoritmus
Gram-Schmidtův proces lze zapsat několika matematickyekvivalentními způsoby. Dvě základní varianty →klasický a modifikovaný algoritmus.
40
Implementace Gram-Schmidtova procesuKlasický Gram-Schmidtův algoritmus
Gram-Schmidtův proces lze zapsat několika matematickyekvivalentními způsoby. Dvě základní varianty →klasický a modifikovaný algoritmus.
Výpočet vektoru z v k-tém kroku lze zapsat ve tvaru
z = (I −k−1∑
i=1
qiq∗i ) ak ≡ CQ ak ,
kde ortogonální projektor
CQ = I −k−1∑
i=1
qiq∗i = I − Qk−1Q∗
k−1
je komplementárním projektorem k projektoru Qk−1Q∗k−1.
40
Implementace Gram-Schmidtova procesuKlasický Gram-Schmidtův algoritmus
Gram-Schmidtův proces lze zapsat několika matematickyekvivalentními způsoby. Dvě základní varianty →klasický a modifikovaný algoritmus.
Výpočet vektoru z v k-tém kroku lze zapsat ve tvaru
z = (I −k−1∑
i=1
qiq∗i ) ak ≡ CQ ak ,
kde ortogonální projektor
CQ = I −k−1∑
i=1
qiq∗i = I − Qk−1Q∗
k−1
je komplementárním projektorem k projektoru Qk−1Q∗k−1.
Tento zápis umožňuje paralelní implementaci.
40
Klasický Gram-Schmidtův algoritmus (CGS)Operátor CQ
input A = [a1, . . . , am]r11 := ‖a1‖q1 := a1/r11
Q1 := [q1]for k = 2 : m do
z := ak
[r1,k, . . . , rk−1,k]T := Q∗k−1z
z := z − Qk−1[r1,k, . . . , rk−1,k]T
rkk := ‖z‖qk := z/rkk
Qk := [Qk−1, qk]end for
41
Implementace Gram-Schmidtova procesuModifikovaný Gram-Schmidtův algoritmus
Z ortonormality vektorů q1, . . . , qk−1 plyne (cvičení)
I −k−1∑
i=1
qiq∗i
︸ ︷︷ ︸CQ
= (I − qk−1q∗k−1) . . . (I − q2q∗
2)(I − q1q∗1)
︸ ︷︷ ︸MQ
.
Vektor z = CQak lze tedy ekvivalentně spočít pomocí
z = MQ ak = (I − qk−1q∗k−1) . . . (I − q2q∗
2)(I − q1q∗1) ak .
Implementace využívající MQ je sekvenční
z1 = ak
z2 = z1 − (q∗1z1) q1 ,
...
zk = zk−1 − (q∗k−1zk−1) qk−1 .
42
Modifikovaný Gram-Schmidtův algoritmus (MGS)Operátor MQ
input A = [a1, . . . , am]r11 := ‖a1‖q1 := a1/r11
Q1 := [q1]for k = 2 : m do
z := ak
for i = 1 : k − 1 dorik := q∗
i zz := z − rikqi
end forrkk := ‖z‖qk := z/rkk
Qk := [Qk−1, qk]end for
43
Iterační zpřesnění
Vztahy používající CQ a MQ jsou matematicky ekvivalentní.Vedou však na algoritmy s různými numerickými vlastnostmi.
44
Iterační zpřesnění
Vztahy používající CQ a MQ jsou matematicky ekvivalentní.Vedou však na algoritmy s různými numerickými vlastnostmi.
CGS počítá méně přesněji, u MGS jsou chyby vznikající přivýpočtu z částečně eliminovány, výpočet ri,k := q∗
i z následujeaž po výpočtu aktualizace.
44
Iterační zpřesnění
Vztahy používající CQ a MQ jsou matematicky ekvivalentní.Vedou však na algoritmy s různými numerickými vlastnostmi.
CGS počítá méně přesněji, u MGS jsou chyby vznikající přivýpočtu z částečně eliminovány, výpočet ri,k := q∗
i z následujeaž po výpočtu aktualizace.
Ortogonalizační krok můžeme opakovat → CGS a MGSs iteračním zpřesněním → ICGS, IMGS.
44
Iterační zpřesnění
Vztahy používající CQ a MQ jsou matematicky ekvivalentní.Vedou však na algoritmy s různými numerickými vlastnostmi.
CGS počítá méně přesněji, u MGS jsou chyby vznikající přivýpočtu z částečně eliminovány, výpočet ri,k := q∗
i z následujeaž po výpočtu aktualizace.
Ortogonalizační krok můžeme opakovat → CGS a MGSs iteračním zpřesněním → ICGS, IMGS.
ICGS si zachovává hlavní výhodu CGS – velmi dobře separalelizuje a navíc redukuje ztrátu ortogonality díkyopakované ortogonalizaci.
44
Iterační zpřesnění
Vztahy používající CQ a MQ jsou matematicky ekvivalentní.Vedou však na algoritmy s různými numerickými vlastnostmi.
CGS počítá méně přesněji, u MGS jsou chyby vznikající přivýpočtu z částečně eliminovány, výpočet ri,k := q∗
i z následujeaž po výpočtu aktualizace.
Ortogonalizační krok můžeme opakovat → CGS a MGSs iteračním zpřesněním → ICGS, IMGS.
ICGS si zachovává hlavní výhodu CGS – velmi dobře separalelizuje a navíc redukuje ztrátu ortogonality díkyopakované ortogonalizaci.
Pro dosažení maximální přesnosti (ztráty ortogonality naúrovni úměrné strojové přesnosti) stačí jediné opakováníortogonalizace.
44
Iterační zpřesnění
Vztahy používající CQ a MQ jsou matematicky ekvivalentní.Vedou však na algoritmy s různými numerickými vlastnostmi.
CGS počítá méně přesněji, u MGS jsou chyby vznikající přivýpočtu z částečně eliminovány, výpočet ri,k := q∗
i z následujeaž po výpočtu aktualizace.
Ortogonalizační krok můžeme opakovat → CGS a MGSs iteračním zpřesněním → ICGS, IMGS.
ICGS si zachovává hlavní výhodu CGS – velmi dobře separalelizuje a navíc redukuje ztrátu ortogonality díkyopakované ortogonalizaci.
Pro dosažení maximální přesnosti (ztráty ortogonality naúrovni úměrné strojové přesnosti) stačí jediné opakováníortogonalizace.
Jak vypočítat zpřesněné koeficienty ri,k?
44
CGS algoritmus s opakovanou ortogonalizací
První ortogonalizace: projektujeme ak
z = ak − [(q∗k−1ak) qk−1 + . . . + (q∗
1ak) q1] .
45
CGS algoritmus s opakovanou ortogonalizací
První ortogonalizace: projektujeme ak
z = ak − [(q∗k−1ak) qk−1 + . . . + (q∗
1ak) q1] .
Opakovaná ortogonalizace: projektujeme z
w = z − [(q∗k−1z) qk−1 + . . . + (q∗
1z) q1] .
45
CGS algoritmus s opakovanou ortogonalizací
První ortogonalizace: projektujeme ak
z = ak − [(q∗k−1ak) qk−1 + . . . + (q∗
1ak) q1] .
Opakovaná ortogonalizace: projektujeme z
w = z − [(q∗k−1z) qk−1 + . . . + (q∗
1z) q1] .
Dosadíme-li za z do druhého vztahu,
w = ak − [(q∗k−1ak + q∗
k−1z) qk−1 + . . . + (q∗1ak + q∗
1z) q1] .
Zpřesněné koeficienty (prvky horní trojúhelníkové R)vzniknou součtem projekcí ak a z do jednotlivých směrů.
45
Iterovaný klasický Gram-Schmidtův algoritmus
input A = [a1, . . . , am]r11 := ‖a1‖q1 := a1/r11
Q1 := [q1]for k = 2 : m do
z := ak
r := 0for ℓ = 1 : 2 do
r̃ := Q∗k−1z
z := z − Qk−1r̃r := r + r̃
end for[r1,k, . . . , rk−1,k]T := rrkk := ‖z‖qk := z/rkk
Qk := [Qk−1, qk]end for
46
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
47
Cena výpočtu Gram-Schmidtova procesu a QR rozkladuČtvercová matice
Měříme počtem aritmetických operací.
Není třeba rozlišovat CGS a MGS (počet operací je stejný).ICGS vyžaduje přibližně dvakrát více operací než CGS a MGS.
Cena k-tého kroku (cvičení):
z = ak −k−1∑
i=1
(q∗i ak)qi , qk =
z
‖z‖ =z√z∗z
,
2n(k − 1) + 2n násobení, (2n − 1)(k − 1) + n − 1 sčítání neboodčítání, jednu odmocninu a jedno dělení.
Sečtením všech operací pro kroky k = 1, . . . , n zjistíme, žeCGS (MGS) stojí približně
2n3 aritmetických operací.
48
Cena výpočtu QR rozkladupomocí Gram-Schmidtova procesu
Operace k-tý krok celkem
× 2n(k − 1) + 2n n3 + n2
+, − (2n − 1)(k − 1) + n − 1 n3 − 12n2 − 1
2n: 1 n√
1 n
49
Výpočetní ceny různých implementací QR rozkladuObdélníková matice
Algoritmus celkový počet operací
n > m n = m
Householderův QR rozklad 2nm2 − 2m3/3 4/3n3
Givensův QR rozklad 3nm2 − m3 2n3
CGS, MGS 2nm2 2n3
Počet operací u Householderových reflexí je počet operacípotřebných ke spočtení rozkladu v součinovém tvaru
A = H∗1 H∗
2 . . . H∗m−1R,
kde Hi jsou jednoduché Householderovy reflexe.
50
Výpočetní ceny různých implementací QR rozkladuObdélníková matice
Algoritmus celkový počet operací
n > m n = m
Householderův QR rozklad 2nm2 − 2m3/3 4/3n3
Givensův QR rozklad 3nm2 − m3 2n3
CGS, MGS 2nm2 2n3
Počet operací u Householderových reflexí je počet operacípotřebných ke spočtení rozkladu v součinovém tvaru
A = H∗1 H∗
2 . . . H∗m−1R,
kde Hi jsou jednoduché Householderovy reflexe.
Explicitní znalost matice Q vyžaduje dalších
2nm2 − 2m3/3 operací.
50
Výpočetní ceny různých implementací QR rozkladuObdélníková matice
Algoritmus celkový počet operací
n > m n = m
Householderův QR rozklad 2nm2 − 2m3/3 4/3n3
Givensův QR rozklad 3nm2 − m3 2n3
CGS, MGS 2nm2 2n3
Počet operací u Householderových reflexí je počet operacípotřebných ke spočtení rozkladu v součinovém tvaru
A = H∗1 H∗
2 . . . H∗m−1R,
kde Hi jsou jednoduché Householderovy reflexe.
Explicitní znalost matice Q vyžaduje dalších
2nm2 − 2m3/3 operací.
Totéž u Givensových rotací.
50
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
51
Numerická stabilita QR rozkladu
Říkáme, že algoritmus je zpětně stabilní, pokud spočtenéřešení dané úlohy je přesným řešením modifikované úlohy sevstupními daty, která jsou blízká původním datům dané úlohy.
Nechť Q̂ a R̂ jsou spočtené faktory QR rozkladu maticeA ∈ Cn×m v aritmetice s konečnou přesností.
Nechť fl(UA) je spočtený výsledek aplikace elementárníGivensovy rotace nebo Householderovy reflexe U na matici Av aritmetice s konečnou přesností.
52
Numerická stabilita QR rozkladu
Říkáme, že algoritmus je zpětně stabilní, pokud spočtenéřešení dané úlohy je přesným řešením modifikované úlohy sevstupními daty, která jsou blízká původním datům dané úlohy.
Nechť Q̂ a R̂ jsou spočtené faktory QR rozkladu maticeA ∈ Cn×m v aritmetice s konečnou přesností.
Nechť fl(UA) je spočtený výsledek aplikace elementárníGivensovy rotace nebo Householderovy reflexe U na matici Av aritmetice s konečnou přesností.
Otázky:Jaká je přesnost vypočtené matice fl(UA)?Co lze říci o ortogonalitě spočtené matice Q̂?Jak přesně spočtený je faktor R̂?Do jaké míry odpovídá součin Q̂R̂ matici A?
52
Aplikace Givensových rotací a Householderových reflexív konečné aritmetice
Wilkinson, Turing: Je-li U ∈ Cn×n matice elementárníGivensovy rotace nebo Householderovy reflexe, potom
fl(UA) = U(A + E) , kde‖E‖‖A‖ ≤ γn2
u + O(u2) ,
γ malá konstanta.
Faktor n2 je velmi pesimistický a není důležitý z hlediskapochopení fungování algoritmu v konečné aritmetice počítače.
Výše zmíněný odhad normy E budeme zjednodušeně zapisovat
‖E‖ ≈ u‖A‖.
53
Numerická stabilita QR rozkladuOrtogonalita spočtené matice Q̂
Ztrátu ortogonality matic vypočtených v aritmetices konečnou přesností budeme měřit pomocí 2-normy matice
EQ ≡ Q̂∗Q̂ − I.
Velikosti normy EQ pro jednotlivé algoritmy:
Algoritmus ‖EQ‖Householderův QR rozklad u
Givensův QR rozklad u
CGS κ2(A) u
MGS κ(A) u
ICGS u
Lze matematicky dokázat.
54
Numerická stabilita QR rozkladuStabilita výpočtu faktoru R̂
Pomocí analýzy zaokrouhlovacích chyb lze ukázat, žeHouseholderův QR, Givensův QR a MGS jsou zpětněstabilní algoritmy v následujícím smyslu:
Nechť R̂ je spočtený trojúhelníkový faktor. Potom existujeunitární matice Q a perturbace E tak, že platí
Q∗(A + E) = R̂, ‖E‖F ≈ u‖A‖F .
Překvapivé hlavně pro MGS, kde není zaručená ztrátaortogonality na úrovni strojové přesnosti.
55
Numerická stabilita QR rozkladuNorma rezidua A − Q̂R̂
Matice Q̂ a R̂ vypočtené všemi uvažovanými variantami QRrozkladů splňují
‖A − Q̂R̂‖ ≈ u‖A‖.
Platí i pro CGS, kde je ztráta ortogonality největší.
56
Outline
1 Ortogonální projekce a ortogonální projektory
2 Givensovy rotace
3 Householderovy reflexe
4 QR rozklad
5 QR rozklad a Gram-Schmidtův ortogonalizační proces
6 Cena výpočtu QR rozkladu
7 Numerická stabilita QR rozkladu
8 Cvičení
57
Cvičení
3.3 Určete vlastní čísla a determinant elementární Givensovyrotace.
3.5 Určete vlastní čísla a determinant Householderovy reflexev reálném oboru.
3.8 Ukažte, že Householderovy reflexe odpovídající vzájemněortogonálním normálovým vektorům jsou komutativní.
58
Cvičení
3.11 Mějme dánu matici Givensovy rotace v R2
G(ϕ) =
[cos(ϕ) − sin(ϕ)sin(ϕ) cos(ϕ)
].
Definujme dva normalizované vektory
q1 =√
0.5
[ √1 − cos(ϕ)√1 + cos(ϕ)
], q2 =
[01
],
‖q1‖ =1
2(1 − cos(ϕ) + 1 + cos(ϕ)) = 1 = ‖q2‖ .
Ukažte, že
H(q2)H(q1) =
[cos(ϕ) − sin(ϕ)sin(ϕ) cos(ϕ)
]= G(ϕ) ,
tedy každou Givensovu rotaci lze popsaným způsobem složitze dvou Householderových reflexí.
59
Cvičení
3.13 Dokažte, že matici Householderovy reflexe nelze složitz Givensových rotací. Dále dokažte, že pomocí Givensovýchrotací nelze vyjádřit žádná matice, která je součinem lichéhopočtu Householderových reflexí.
3.15 Ukažte, že jsou-li vektory q1, . . . , qk navzájem ortogonální,pak platí
I − q1q∗1 − q2q∗
2 . . . − qkq∗k = (I − q1q∗
1)(I − q2q∗2) . . . (I − qkq∗
k),
tedy součin na pravé straně je komutativní.
3.16 Ukažte, že cena k-tého kroku algoritmu CGS je dána cenou2n(k − 1) + 2n součinů, (2n − 1)(k − 1) + n − 1 součtů neborozdílů, jednoho podílu a jednoho výpočtu odmocniny.Využijte sčítací vzorce
n∑
k=1
k =1
2n(n + 1) a
n∑
k=1
k2 =1
6n (2n + 1) (n + 1) .
60
