Téma: Hauseholderova identita a možnost jejího využití Zpracoval Doc. RNDr. Zdeněk Hlaváč, CSc
1
Hauseholderova identita
Nechť K = [kij ] je matice typu (m, n), která je řídká v tom smyslu, že obsahuje samé nulové řádky s výjimkou r < m řádků pořadových čísel i1 , . . . , ir a samé nulové sloupce s výjimkou s < n sloupců pořadových čísel j1 , . . . , js . Všechny informace o této matici uchováme, budeme-li znát • její typ (m, n), • vektory i = [ia ] a j = [jb ] pořadových čísel nenulových řádků a sloupců včetně jejich dimenzí r a s, C • tzv. komprimovanou matici K C = [kab ] typu (r, s), pro kterou je C kab = kia jb ; a = 1, . . . , r ; b = 1, . . . , s .
Definujme tzv. Booleovskou matici B pl příslušející k vektoru pořadových čísel l dimenze r a k přirozenému číslu p ≥ r jako matici typu (r, p) složenou z nul a jedniček, ve které jedničky jsou pouze na pozicích (a, la ) pro a = 1, . . . , r. Pomocí této matice lze provádět některé maticové operace, jak jest obsahem následujícího tvrzení: Tvrzení 1: Nechť K C je matice typu (r, s) (komprimovaná, tj. ”malého” typu). Potom C C T • (B m i ) K je matice typu (m, s) s řádky matice K umístěnými v řádcích s pořadovými čísly i1 , . . . , ir . Ostatní řádky této matice jsou nulové.
• K C B nj je matice typu (r, n) se sloupci matice K C umístěnými ve sloupcích s pořadovými čísly j1 , . . . , js . Ostatní sloupce této matice jsou nulové. Nechť K je matice typu (m, n) (nekomprimovaná, to jest ”velkého” typu). Potom • Bm i K je matice typu (r, n), v níž na řádku pořadového čísla a ; a = 1, . . . , r je umístěn řádek pořadového čísla ia matice K a ostatní řádky matice K jsou vynechány. • K(B nj )T je matice typu (m, s), ve které ve sloupci pořadového čísla a ; a = 1, . . . , s je umístěn sloupec pořadového čísla ja matice K a ostatní sloupce matice K jsou vynechány. Ověření tohoto tvrzení plyne bezprostředně z definice Booleovské matice provedením ve větě naznačených maticových násobení. Tvrzení 2: Mezi maticemi K typu (m, n) a K C typu (r, s) definovanými výše platí následující vztahy 1
C n T K = (B m i ) K Bj ,
(1)
n T K C = Bm i K(B j ) .
(2)
Bez ověření formulujeme hlavní tvrzení tohoto odstavce. Tvrzení 3 (Hauseholderova identita:) Nechť A je regulární čtvercová matice řádu n, X a Y matice typu (p, n). Je-li matice A + X T Y regulární, platí (A + X T Y )−1 = A−1 − A−1 X T (E p + Y A−1 X T )−1 Y A−1 .
2
(3)
Využití Hauseholderovy identity
Použijme nyní identity (3) na inverzi matice D = E n − KG ”velkého” řádu n, přičemž matice G = [gij ] je obecně plná typu (p, n) a matice K typu (n, p) je řídká ve výše popsaném smyslu. Dosazením za K z (1) do definičního vztahu pro matici D, kdy m v (1) je nyní n a n v (1) je nyní p, dostáváme D = E n − (B ni )T K C B pj G . Uzávorkováním menšitele po dvojicích činitelů vidíme, že pro inverzi matice D prostřednictvím Hauseholderovy identity stačí ve vztahu (1) vzít A = E n ; X T = (B ni )T K C ; Y = B pj G .
(4)
Podle této identity, protože je A−1 = E −1 n = E n , platí D
−1
n T
= E n − (B i ) K
C
p
n T
E s + B j G(B i ) K
C
−1
B pj G .
(5)
Matice, kterou je nutno na pravé straně (5) invertovat, je ”malého” řádu s, protože podle tvrzení 1 matice B pj G je typu (s, n) a matice (B ni )T K C je typu (n, s). Proto i jednotková matice v invertovaném výrazu v (5) je pouze řádu s a nikoliv řádu p, jak by vyplývalo z přímé aplikace Hauseholderovy identity. Podle tvrzení 1 zřejmě je
(B ni )T K C
C k11 = ... C kr1
C , . . . , k1s i − tý řádek 1 ... ... C , . . . , krs ir − tý řádek .
(6)
Řádky mající jiné pořadí než i1 −té až ir −té této matice typu (n, s) jsou nulové. Podle téhož tvrzení dále platí
gj1 1 , . . . , gj1 n ... ... B pj G = ... . gjs 1 , . . . , gjs n
(7)
Odtud plyne, že matice řádu s, kterou je nutno na pravé straně vztahu (5) invertovat, má tvar 2
"
p
E s + B j G(B ni )T K c = δab +
r X
#
gja il klbC ; a, b = 1, . . . , s .
l=1
(8)
V tomto výrazu δab je Kroneckerův symbol δab = 0 pro a 6= b a δaa = 1. Jestliže tuto matici např. Gaussovou eliminací zinvertujeme a inverzní matici (”malého” řádu s) označíme C = [cij ], podle (7), (6) a (5) zřejmě pro matici D −1 vychází
D
−1
X s X s C k1l cll′ gjl′ 1 l=1 l′ =1 ... = En − s X s X C krl cll′ gjl′ 1 l=1 l′ =1
,..., ... ,...,
s X s X
l=1 l′ =1
s X s X
l=1 l′ =1
C k1l cll′ gjl′ n
... C krl cll′ gjl′ n
i1 − tý řádek ir − tý řádek
(9)
Řádky mající jiné pořadí než i1 −té až ir −té této matice řádu n v menšiteli jsou nulové. Matice D −1 tedy obsahuje čtyři typy prvků: 1. Prvky tvaru
s X s X
C kal cll′ gjl′ b na pozicích (ia , b) ; a = 1, . . . , r ; b = 1, . . . , n pro
l=1 l =1 ′
ia 6= b.
2. Prvky tvaru 1 −
s s X X
C kal cll′ gjl′ b na těchže pozicích, ovšem pro ia = b.
l=1 l′ =1
3. Nulové prvky na pozicích (j, b) pro j = 1, . . . , n ; b = 1, . . . , n ; j 6= ia ; j 6= b . 4. Jednotkové prvky na těchže pozicích jako v předchozím bodě, ovšem pro j = b. Pro dokonalý popis této matice tedy postačí znát (kromě jejího řádu n a souřadnic vektoru výběru i) prvky, které nejsou ani nulové, ani jednotkové. Tyto prvky však zaplňují pouze pole typu (r, n), tedy rn prvků místo pole typu (n, n) tedy n2 prvků. Pro výpočet matice D −1 je potřeba znát řídkou matici K (tedy její typ (n, p), vektor i, jehož souřadnice popisují nenulové řádky, vektor j, jehož souřadnice popisují nenulové sloupce této matice a příslušnou komprimovanou matici K C ) a plnou matici G. Pro výpočet matice D −1 potom postupujeme podle následujícího algoritmu: 1. Určíme matici E s + B pj G(B ni )T K C podle (8). 2. Tuto matici ”malého” řádu s invertujeme. Vznikne matice C = [cij ]. 3. Podle (9) určíme matici D −1 . Poznámka: Z (8) a (9) plyne, že k uskutečnění popsaného algoritmu si z plné matice G stačí pamatovat pouze řádky popsané vektorem výběru j, tedy pouze pole typu (s, n), což jest sn prvků. Místo matice K si stačí pamatovat pouze vektory i, j a komprimovanou matici K C , to jest rs + r + s prvků. Z matice D −1 si stačí pamatovat pouze řádky zaplněné nenulovými a nejednotkovými prvky, tedy pouze pole typu (r, n), čili rn prvků. Kromě toho je třeba si pamatovat prvky matice C (a eventuálně i prvky inverzní matice C −1 podle druhu algoritmu inverze matic). To odpovídá dvěma maticím řádu s, tedy 2s 3
prvkům. Pro uskutečnění celého algoritmu inverze matice D (včetně vstupů a výstupu) je třeba rezervovat místo pro rs + ns + 2s2 + nr + r + s = (r + s)(n + s + 1) + s2 prvků. Při nevyužití Hauseholderovy identity a řídkosti matice K by bylo ekvivalentně potřeba si pamatovat matici K (tedy np prvků), matici G (rovněž np prvků), matici D a matici D −1 (pro obě matice celkem n2 prvků). Pro uskutečnění celého algoritmu to znamená rezervovat paměť pro 2np + 2n2 = 2n(n + p) prvků. Např. pro n = p = 100 a r = s = 50 se jedná o 17600 prvků při využití popisovaného postupu, oproti 40000 prvků při nevyužití tohoto postupu. Je patrno, že jde o podstatnou paměťovou úsporu, která je tím markantnější, čím je matice K řidší (a tedy konstanty r a s menší). Zmíněná úspora umožňuje bez použití vnější paměti invertovat podstatně větší matice. Kromě úspory paměťové dochází rovněž k úspoře časové, takže při inverzi matice řádu s se vyplatí provést navíc násobení matic naznačená ve vztazích (8) a (9). Použijme v dalším identitu (3) na inverzi matice D = E n − KG, kde obecně plná matice G je nyní typu (p, n) a matice K je řídká výše popsaným způsobem a je typu (p, n). Oproti výše diskutovanému případu si pouze plná a řídká matice vyměnily pořadí. Dosazením za K ze (1) do definičního vztahu pro matici D, kdy m v (1) jest nyní p a n v (1) zůstává, dostáváme D = E n − G(B pi )T K C B nj . Uzávorkujeme-li menšitele tohoto výrazu po dvojicích činitelů, je patrno, že pro inverzi matice D použitím Hauseholderovy identity stačí v (3) brát A = E n ; X T = G(B pi )T ; Y = K C B nj .
(10)
Podle (3) pak platí D
−1
p T
= E n − G(B i )
C
p T
n
E r + K B j G(B i )
−1
K C B nj .
(11)
Matice v menšiteli výrazu na pravé straně (11) je čtvercová řádu r. Proto i jednotková matice v tomtéž menšiteli je řádu r. Podle tvrzení 1 zřejmě je
g1i1 , . . . , g1ir p T G(B i ) = . . . . . . . . . gni1 , . . . , gnir
K C B nj =
C k11 .. . C kr1
.. C . k1s .. .. . . .. C . krs
(12)
(13)
V matici na pravé straně (13) jsou zaplněny pouze sloupce pořadových čísel j1 až js tam uvedenými prvky. Ostatní sloupce této matice (která je typu (r, n)) jsou nulové. Odtud plyne, že matice řádu r, kterou je nutno na pravé straně (11) invertovat, je tvaru p T
"
E r + K C B nj G(B i ) = δab +
s X l=1
#
C gjl ib ; a, b = 1, . . . , r . kal
(14)
Inverzí této matice nechť vznikne matice C = [cij ] řádu r. Podle (11), (12) a (13) pak pro matici D −1 vychází 4
r r X r r X . .. X .. X C ′ g1il cll′ klC′ s .. g1il cll kl′ 1 . .
D −1
= En −
l=1 l′ =1
.. .. . . r X r X .. . gnil cll′ klC′ 1 l=1 l′ =1
l=1 l′ =1
.. .. . . r X r X .. . gnil cll′ klC′ s l=1 l′ =1
.. . .. .
.
(15)
V matici v menšiteli výrazu na pravé straně (15) jsou zaplněny pouze sloupce pořadových čísel j1 až js tam uvedenými prvky. Ostatní sloupce této matice (která je řádu n) jsou nulové. Poznámky: 1. Podle(14) a (15) postačí k inverzi matice D znát pouze j1 −tý až js −tý sloupec matice G, tedy pole o ns prvcích. 2. Pro uchování matice D −1 si stačí zapamatovat pouze prvky, které nemají ani nulovou ani jednotkovou hodnotu. Postačí proto pro výstup si pamatovat pouze pole typu (n, s) tedy ns prvků. 3. Paměťová úspora oproti stavu, kdy bychom nevyužili řídkosti matice K, je stejná jako výše odvozená pro případ opačného pořadí matic K a G v součinu. Jestliže ve výrazu E n − KG jsou obě matice řídké, nelze s větší výhodou využít řídkosti obou z nich. Využijeme vždy řídkosti toho z činitelů, který je z hlediska inverze matice nízkého řádu výhodnější. Druhý činitel pak považujeme za obecně plný. Jestliže tedy platí K = (B ni1 )T K C B pj , 1 G = (B pi )T GC B nj , 2 2 kde K C je komprimovaná matice typu (r1 , s1 ) a GC komprimovaná matice typu (r2 , s2 ), jeví se, pro případ inverze matice D = E n − KG kdy r2 < s1 , výhodnější uvažovat matici G jako řídkou a K jako plnou, zatímco při platnosti r2 > s1 je naopak výhodnější uvažovat G jako plnou a K jako řídkou. V prvním případě totiž použití metody vede na inverzi matice řádu r2 , zatímco v druhém případě na inverzi matice řádu s1 . Při inverzi matice D = E n − GK je rozhodující relace mezi čísly r1 a s2 . Pro r1 < s2 je výhodnější uvažovat matici K jako řídkou, protože použití metody vede na inverzi matice řádu r1 . Pro relaci r1 > s2 je naopak výhodnější uvažovat jako řídkou matici G, neboť použití metody vede na inverzi matice matice řádu s2 . V některých případech bývají matice K i G čtvercové ”velkého” řádu n a matice K bývá řídká ve výše uvedeném smyslu a navíc symetrická. Odtud plyne, že pro její popis stačí zadat • její řád n, • vektor výběru i = [ia ] pořadových čísel nenulových řádků (a vzhledem k symetrii i sloupců), včetně jeho dimenze r, C • symetrickou komprimovanou matici K C = [kab ] řádu r, pro kterou C kab = kia ib ; a, b = 1, . . . , r .
5
Při inverzi matice D = E n − KG potom postupujeme podle následujícího algoritmu: 1. Určíme matici řádu r o prvcích δab +
r X
gia il klbC ; a, b = 1, . . . , r .
l=1
2. Invertujeme tuto matici. Vznikne matice C = [cij ] řádu r. 3. Určíme matici typu (r, n) o prvcích r r X X
C kal cll′ gil′ b ; a = 1, . . . , r ; b = 1, . . . , n .
l=1 l′ =1
4. Jednotlivé řádky této matice umístíme na místa i1 −tého až ir −tého řádku matice řádu n, doplníme nulovými řádkami, změníme znaménka a přičteme jednotkovou matici řádu n. Výsledkem je matice D −1 . Při inverzi matice D = E n − GK má analogický algoritmus tvar: 1. Určíme matici řádu r o prvcích δab +
r X
C gil ib kal ; a, b = 1, . . . , r .
l=1
2. Invertujeme tuto matici. Vznikne matice C = [cij ] řádu r. 3. Určíme matici typu (n, r) o prvcích r X r X
klC′ b cll′ gail ; a = 1, . . . , n ; b = 1, . . . , r .
l=1 l′ =1
4. Jednotlivé sloupce této matice umístíme na místa i1 −tého až ir −tého sloupce matice řádu n, doplníme nulovými sloupci, změníme znaménka a přičteme jednotkovou matici řádu n. Výsledkem je matice D −1 . Poznámka: Tvrzení o paměťových úsporách zůstává v platnosti položíme-li v úvahách výše n = p ; r = s.
6
