Metody pro řídké matice ⬩ 01MRMMI

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.

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).

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ů.

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.

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.

Eliminační strom

Algoritmy pro vytvoření pásu/profilu

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

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

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

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é.

Indefinitní symetrické matice

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=\left(\begin{array}{cc}B& C\\ C^{𝖳}& 0\end{array}\right)$, kde $B\in {ℝ}^{m\times m}$ je pozitivně definitní, $C\in {ℝ}^{m\times (n-m)},m\ge n-m,h\left(C\right)=n-m$. Soustava je tedy ve tvaru

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-Cy)$. Dosazením do druhé rovnice máme

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-\left\{0\right\}$ je

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.

K nalezení permutace řádků a sloupců, které nám umožní matici převést na blokově trojúhelníkový tvar, použijeme orientovaný graf matice. Najdeme silně souvislé komponenty, topologicky je seřadíme a z každé uděláme blok, přičemž je řadíme v opačném pořadí.

K rozkladu jednotlivých bloků se dá použít částečná pivotace (pivota hledáme jen v příslušném řádku/sloupci) nebo úplná pivotace (pivota hledáme v celé podmatici).

Předpodmiňování stacionárních metod

Nechť $A=M-N$, kde $M$ je regulární matice. Potom

Budeme hledat posloupnost $\left(x_k\right)$ takovou, že $M{x}_{k+1}=N{x}_k+b$. Z toho máme

kde $r_k$ je reziduum $k$-té aproximace. Stačí tedy vždy spočítat reziduum $r≔b-Ax$, spočítat opravu $Mz=r$ a opravit $x\leftarrow x+z$. Tato iterace konverguje pro $\rho \left(M^{-1}N\right)=\rho (𝐈-M^{-1}A)<1$.

Rychlost konvergence stacionárních metod

Nechť $A=M-N$ je regulární rozklad a $G≔M^{-1}N$. Máme aproximaci $x_{k+1}=G{x}_k+M^{-1}b$ a chceme určit, jak se liší od přesného řešení $x=Gx+M^{-1}b$. Označme $e_k≔x_k-x$. Podle předpisu iterace je $e_{k+1}=G{x}_k-Gx=G{e}_k$. Rozvinutím dostáváme

kde $e_0$ je počáteční chyba. Z toho plyne odhad $\Vert e_k\Vert \le \Vert G\Vert \Vert e_0\Vert$. Podle nějaké věty je $\underset{n\to \infty }{lim}\sqrt[n]{\Vert G^n\Vert }=\rho \left(G\right)$. Když to rozepíšeme z definice, máme

Jednoduchými úpravami dostaneme $\Vert G^n\Vert \le {(\rho \left(G\right)+\epsilon )}^n$. Pokud $\rho \left(G\right)<1$, můžeme zvolit $\epsilon$ takové, aby $\rho \left(G\right)+\epsilon <1$, takže od nějakého členu dál bude chyba exponenciálně klesat. Nevíme ovšem, jak se to bude chovat na začátku. Konvergence dokonce ani nemusí být monotónní, tedy může se stát, že na začátku bude chyba růst.

Speciálně je-li $G$ normální, máme $\rho \left(G\right)={\Vert G\Vert }_2$, takže v tomto případě máme monotónní konvergenci zaručenou. Počet iterací potřebný k dosažení chyby nanejvýš $\delta$ v poměru k počáteční chybě je potom ${log}_{{\Vert G\Vert }_2}\delta$.

Volba předpodmínění u stacionární metody

Vezmeme-li $M≔𝐈$, máme metodu prosté iterace, která konverguje, pokud $\rho (𝐈-A)<1$. Ideálně $\rho (𝐈-A)\ll 1$, tedy vlastní čísla $A$ jsou všechna blízká $1$. Pokud tomu tak není, musíme předpodmiňovat.

Richardsonovo předpodmínění

$M≔\frac{1}{\omega }𝐈$, kde $\omega >0$ je relaxační parametr. Potom o konvergenci rozhoduje $\rho (𝐈-\omega A)$. Často je možné vhodně zvolit $\omega$ tak, aby konvergence fungovala, i když může být pomalá.

Jacobiho předpodmínění

$M≔D$, kde $D$ je diagonální matice se stejnou diagonálou jako $A$. O konvergenci potom rozhoduje $\rho (𝐈-D^{-1}A)$. O tom obecně nemůžeme říct, ale občas se to hodí a je to hodně jednoduché na výpočet.

Gaussovo–Seidelovo předpodmínění

$M≔D-L$, kde $L$ je dolní trojúhelník (i s diagonálou) $A$. $M^{-1}$ potom získáme řešením soustavy s trojúhelníkovou maticí, což se taky počítá snadno, ale na rozdíl od Jacobiho nejde snadno paralelizovat.

SOR metoda

$M≔\frac{1}{\omega }D-L$, kde $\omega >0$ je relaxační parametr. Ten volíme tak, aby $\rho (𝐈-{(\frac{1}{\omega }D-L)}^{-1}A)$ bylo minimální. Ještě se dostaneme k tomu, jak to udělat.

Chceme-li předpodmínění použít u krylovských metod pracujících se symetrickou maticí (např. metoda sdružených gradientů), musíme ho zvolit tak, aby se symetrie zachovala. U Richardsona nebo Jacobiho to není problém, ale SOR musíme upravit, čímž vznikne SSOR (Symmetric SOR). Definujme $M_1≔\frac{1}{\omega }D-L,M_2≔\frac{1}{\omega }D-U$, kde $U$ je horní trojúhelníková část $A$. Potom budeme střídavě řešit

kde $N_{1,2}≔M_{1,2}-A$. Ve výsledku máme

Předpodmínění metody sdružených gradientů

Je-li $M$ symetrická, existuje Choleského rozklad $M=L\cdot L^{𝖳}$. Potom můžeme metodu sdružených gradientů aplikovat na rovnici

Mějme pozitivně definitní matici $A\in {ℝ}^{n\times n}$ a počáteční reziduum $r_0≔b-A{x}_0$. Vyřešíme si rovnici $M{z}_0=r_0$. Následně v metodě budeme na některých místech místo $r_i$ používat $z_i$ a po každém kroku si dopočteme další $z_i$. Rychlost konvergence bude

Všimněme si, že to vždy bude konvergovat, i když pro velké $\varkappa \left(A\right)$ hodně pomalu.

Neúplné LU rozklady

Při úplném LU rozkladu vzniká hodně zaplnění. Místo toho můžeme spočítat jen nějakou jeho aproximaci a tu následně použít pro předpodmínění.

Klasický LU rozklad (řádková varianta) vypadá takto:

1. Pro každé $i=2,\dots ,n$:
   1. Pro každé $k=1,\dots ,i-1$:
      1. $a_{i,k}\leftarrow \frac{a_{i,k}}{a_{k,k}}$
      2. Pro každé $j=k+1,\dots ,n$:
         1. $a_{i,l}\leftarrow a_{i,j}-a_{i,k}\cdot a_{k,j}$

Neúplný LU rozklad bude vypadat stejně, ale s tím, že přepisujeme jen některé prvky. Základní varianta je ILU(0), kde přepisujeme jen ty prvky, které už jsou nenulové. Obecněji si můžeme předepsat množinu pozic $P$, u kterých ignorujeme zaplnění. Jelikož diagonální pozice jsou důležité, do $P$ je nikdy nebudeme zahrnovat, tedy

Chceme-li přesnější neúplný LU rozklad, můžeme připustit i nějaké větší zaplnění.