Metody pro řídké matice ⬩ 01MRMMI ⬩ Adamátorovy zápisky

Na konci se bude psát zápočtový test prověřující praktické znalosti. Není potřeba umět důkazy.

Reprezentace řídkých matic v počítači

Chceme řídkou matici uložit tak, aby zabrala co nejméně paměti a efektivně se s ní počítalo.

Souřadnicový formát

Jedno pole reálných řísel obsahující hodnoty všech $k$ nenulových prvků;
dvě celočíselná pole délky $k$ obsahující jejich indexy.

(\begin{matrix} 1 & 0 & 0 & 2 & 0 \\ 3 & 4 & 0 & 5 & 0 \\ 6 & 0 & 7 & 8 & 9 \\ 0 & 0 & 10 & 11 & 0 \\ 0 & 0 & 0 & 0 & 12 \end{matrix})

Nevýhoda je, že v takové reprezentaci je těžké například najít prvek na dané pozici.

Komprimované řídké řádky (CSR)

Jedno pole reálných řísel obsahující hodnoty všech $k$ nenulových prvků seřazené po řádcích;
jedno celočíselné pole délky $k$ obsahující jejich sloupcové indexy;
jedno celočíselné pole délky $n + 1$ , kde $i$ -tý prvek obsahuje index začátku $i$ -tého řádku v předchozích dvou polích. Ten jeden navíc ukazuje za konec pole (tím se zjednoduší algoritmy iterující přes daný řádek).

(\begin{matrix} 1 & 0 & 0 & 2 & 0 \\ 3 & 4 & 0 & 5 & 0 \\ 6 & 0 & 7 & 8 & 9 \\ 0 & 0 & 10 & 11 & 0 \\ 0 & 0 & 0 & 0 & 12 \end{matrix})

Oproti souřadnicovému formátu to šetří paměť. Nevýhodou ovšem je, že je obtížné přidávat nové nenulové prvky, což může být důležité třeba u metod založených na Gaussově eliminaci.

Analogicky můžeme mít komprimované řídké sloupce (CSC).

Modifikované komprimované řídké řádky (MSR)

Motivací je, že často potřebujeme efektivnější přístup k prvkům na hlavní diagonále.

Jedno pole reálných řísel obsahující hodnoty všech $k$ nenulových prvků počínaje diagonálními, poté zarážku a zbylé prvky po řádcích;
jedno celočíselné pole délky $k + 1$ obsahující nejprve indexy začátků řádků a poté sloupcové indexy nediagonálních prvků.

(\begin{matrix} 1 & 0 & 0 & 2 & 0 \\ 3 & 4 & 0 & 5 & 0 \\ 6 & 0 & 7 & 8 & 9 \\ 0 & 0 & 10 & 11 & 0 \\ 0 & 0 & 0 & 0 & 12 \end{matrix})

Kromě efektivního přístupu k diagonále má stejné výhody a nevýhody jako CSR.

Ellpack

Tento formát je populární na vektorových počítačích.

Využívá se pro matice, kde počet nenulových prvků v jednom řádku je omezený nějakou konstantou $N_{d}$ .

Jedno pole reálných čísel velikosti $n \times N_{d}$ obsahující hodnoty nenulových prvků v jednotlivých řádcích;
jedno celočíselné pole velikosti $n \times N_{d}$ obsahující jejich indexy.

(\begin{matrix} 1 & 0 & 2 & 0 & 0 \\ 3 & 4 & 0 & 5 & 0 \\ 0 & 6 & 7 & 0 & 8 \\ 0 & 0 & 9 & 10 & 0 \\ 0 & 0 & 0 & 11 & 12 \end{matrix})

Formát je jednoduchý a většinou se s ním dobře pracuje, akorát nejde použít, když má matice nějaké husté řádky.

Diagonální uložení

Vhodné pro matice, které mají nenulové prvky pouze podél několika málo diagonál (například ty vzniklé z metody konečných diferencí).

Prostě uložíme hodnoty těchto diagonál a k tomu si poznamenáme, kolikáté jsou.

(\begin{matrix} 1 & 0 & 2 & 0 & 0 \\ 3 & 4 & 0 & 5 & 0 \\ 0 & 6 & 7 & 0 & 8 \\ 0 & 0 & 9 & 10 & 0 \\ 0 & 0 & 0 & 11 & 12 \end{matrix})

Vznik zaplnění

Mějme matici $A$ , na kterou aplikujeme první krok podmaticového algoritmu. Tím se změní její zaplnění:

A = (\begin{matrix} \times & \times & \times & \times & 0 \\ \times & \times & 0 & 0 & 0 \\ \times & 0 & \times & 0 & 0 \\ \times & 0 & 0 & \times & \times \\ 0 & 0 & 0 & \times & \times \end{matrix}) \mapsto A^{(1)} = (\begin{matrix} \times & \times & \times & \times & 0 \\ \times & \times & \times & \times & 0 \\ \times & \times & \times & \times & 0 \\ \times & \times & \times & \times & \times \\ 0 & 0 & 0 & \times & \times \end{matrix}) .

Obecně to funguje tak, že v $k$ -tém kroku se všechny sousední vrcholy $k$ propojí do kliky. To je formálně vyjádřeno v následujícím lemmatu:

A \in ℝ^{n \times n}

A \in ℝ^{n \times n}

Eliminační strom

A \in ℝ^{n \times n}

A = (\begin{matrix} \times & 0 & 0 & \times & 0 \\ 0 & \times & 0 & 0 & \times \\ 0 & 0 & \times & \times & \times \\ \times & 0 & \times & \times & 0 \\ 0 & \times & \times & 0 & \times \end{matrix}) \mapsto L = (\begin{matrix} \times & 0 & 0 & 0 & 0 \\ 0 & \times & 0 & 0 & 0 \\ 0 & 0 & \times & 0 & 0 \\ \times & 0 & \times & \times & 0 \\ 0 & \times & \times & \times & \times \end{matrix})

n

L_{i, j} \neq 0

T_{i}

A \in ℝ^{n \times n}

j

L

Přeuspořádání matic

A ≔ (\begin{matrix} \times & \times & \times & \times & \times & \times & \times \\ \times & \times & 0 & 0 & 0 & 0 & 0 \\ \times & 0 & \times & 0 & 0 & 0 & 0 \\ \times & 0 & 0 & \times & 0 & 0 & 0 \\ \times & 0 & 0 & 0 & \times & 0 & 0 \\ \times & 0 & 0 & 0 & 0 & \times & 0 \\ \times & 0 & 0 & 0 & 0 & 0 & \times \end{matrix}) .

A

Je otázkou, jakou permutaci zvolit, abychom si usnadnili práci co nejvíce. K tomu je několik různých přístupů.

Algoritmy pro vytvoření pásu/profilu

A \in ℝ^{n \times n}

x

Uspořádání minimalizující zaplněnost

Chceme najít permutační matici $P$ takovou, že $| struct (L_{P A P^{𝖳}}) |$ je minimální.

Není známý efektivní algoritmus, který by vždy nalezl optimální řešení.

G_{0} ≔ G (A), A_{0} ≔ A

Uspořádání přizpůsobující matici počítačové architektuře

Máme paralelní počítač a chceme vrcholy eliminovat v takovém pořadí, aby se to dalo dobře paralelizovat.

Mějme graf, který se po odstranění nějaké malé podmnožiny vrcholů (separátoru $S$ ) rozdělí na dvě podobně velké komponenty $C_{1}$ a $C_{2}$ . Potom coloki, co budeme dělat s jednou z komponent, ovlivní jen tuto komponentu a separátor. Pro paralelizaci se tedy vyplatí očíslovat vrcholy v pořadí $C_{1} \to C_{2} \to S$ , takže $C_{1}$ a $C_{2}$ potom můžeme eliminovat paralelně.

Rozdělování můžeme rekurzivně opakovat, dokud jsou komponenty dostatečně velké, aby se to vyplatilo. Tomu se říká metoda vnořených řezů (nested dissection method).

Frontální metoda

Místo abychom uspořádání stanovili předem, budeme ho hledat dynamicky za běhu eliminace.

V každém kroku přehodíme řádky a sloupce tak, aby nenulové prvky šly co nejblíž k diagonále.

Často tím vzniknou malé husté podmatice, které můžeme řešit standardními metodami.

Nevýhoda je, že permutování za běhu může být náročné.

Výhoda je, že můžeme konstruovat Choleského rozklad matice, i když ji ještě nemáme celou spočtenou. To se může vyplatit například u metody konečných prvků, kde počítání matice je náročné.

Poznámky k obecnějším systémům

Zatím jsme se zabývali metodami pro pozitivně definitní matice. Jak se tyto poznatky dají zobecnit na jiné třídy matic?

Indefinitní symetrické matice

(\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix})

Symetrické indefinitní matice v praxi vznikají při použití smíšené formulace metody konečných prvků na eliptickou úlohu druhého řádu nebo také při hledání vázaných extrémů metodou Lagrangeových multiplikátorů. V obou těchto případech vznikne matice tvaru $A = (\begin{matrix} B & C \\ C^{𝖳} & 0 \end{matrix})$ , kde $B \in ℝ^{m \times m}$ je pozitivně definitní, $C \in ℝ^{m \times (n - m)}, m \geq n - m, h (C) = n - m$ . Soustava je tedy ve tvaru

\begin{align} B x + C y & = b_{1}, \\ C^{𝖳} x & = b_{2} . \end{align}

Jelikož $B$ je pozitivně definitní, můžeme najít její Choleského rozklad $B = L \cdot L^{𝖳}$ . Podle první rovnice je $x = B^{- 1} (b_{1} - C y)$ . Dosazením do druhé rovnice máme

C^{𝖳} B^{- 1} (b_{1} - C y) = b_{2};

\underset{S}{\underset{⏟}{C^{𝖳} B^{- 1} C}} y = - b_{2} + C^{𝖳} B^{- 1} b_{1},

kde $S \in ℝ^{(n - m) \times (n - m)}$ je Schurův doplněk. Dokážeme, že $S$ je pozitivně definitní. Pro $y \in ℝ^{m} s - {0}$ je

⟨ y | S y ⟩ = ⟨ C^{𝖳} B^{- 1} C y | y ⟩ = ⟨ B^{- 1} C y | C y ⟩ = ⟨ L^{- 𝖳} L^{- 1} C y | C y ⟩ = ⟨ L^{- 1} C y | L^{- 1} C y ⟩ = {‖ L^{- 1} C y ‖}^{2} > 0 .

Soustavu s maticí $S$ tedy můžeme řešit Choleského rozkladem. Akorát je problém v tom, že může být hustá, i když $A$ je řídká. Vyplatí se to tedy, jenom když je malá. V opačném případě je možné to nějak prokletě zkombinovat s iterační metodou (typicky konjugované gradienty). Potom nemusíme $S$ explicitně sestavovat, ale stačí s ní umět násobit vektory. K tomu nám stačí na vektor postupně aplikovat $C, L^{- 1}, L^{- 𝖳}, C^{𝖳}$ .

Obecné regulární matice

U matice, která není symetrická, nemáme Choleského rozklad, ale jen LU rozklad. Ten navíc nnexistuje, když matice není silně regulární; v takovém případě ještě navíc potřebujeme pivotaci. Pivotace by ideálně neměla způsobit velké zaplnění, což je ale často v rozporu s požadavky výše uvedených metod. Navíc výběr pivota obvyklým způsobem závisí na numerických hodnotách, takže není možné analyzovat strukturu předem.

Místo pivotace za běhu můžeme zkusit najít takovou permutaci, která přesune absolutně velké prvky na diagonálu.

(\begin{matrix} \times & 0 & 0 & 0 & 0 & 0 & \times \\ \times & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \times & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \times & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \times & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \times & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \times & \times \end{matrix}) .

Metody pro řídké matice ⬩ 01MRMMI

Úvod