Hierarkisk matris

I numerisk matematik används hierarkiska matriser (H-matriser) som dataglesa approximationer av icke-glesa matriser . Medan en gles matris med dimensionen $n$ effektivt kan representeras i $O(n)$ lagringsenheter genom att endast lagra dess poster som inte är noll, skulle en icke-gles matris kräva $O(n^{2})$ lagringsenheter, och att använda denna typ av matriser för stora problem skulle därför bli oöverkomligt dyrt i termer av lagring och beräkningstid. Hierarkiska matriser ger en approximation som endast kräver $O(nk\,\log(n))$ lagringsenheter, där $k$ är en parameter som styr noggrannheten av approximationen. I typiska tillämpningar, t.ex. när man diskretiserar integralekvationer, förkonditionerar de resulterande systemen av linjära ekvationer eller löser elliptiska partiella differentialekvationer, en rangordning proportionell mot $\log(1/\epsilon )^ {\gamma }$ med en liten konstant $\gamma$ är tillräckligt för att säkerställa en noggrannhet på $\epsilon$ . Jämfört med många andra data-glesa representationer av icke-glesa matriser erbjuder hierarkiska matriser en stor fördel: resultaten av matrisaritmetiska operationer som matrismultiplikation, faktorisering eller inversion kan approximeras i $O(nk^{\alpha }\,\log(n)^{\beta })$ operationer, där $\alpha ,\beta \in \{1,2,3\}.$

Grundläggande idé

Hierarkiska matriser förlitar sig på lokala lågrankade approximationer: låt $I,J$ vara indexuppsättningar och låt $G\in {\mathbb {R} }^{I\ gånger J}$ betecknar matrisen vi måste approximera. I många applikationer (se ovan) kan vi hitta delmängder $t\subseteq I,s\subseteq J$ så att $G|_{t\times s}$ kan approximeras av en rang- $k$ -matris. Denna approximation kan representeras i faktoriserad form $G|_{t\times s}\approx AB^{*}$ med faktorerna $A\in {\ mathbb {R} }^{t\times k},B\in {\mathbb {R} }^{s\times k}$ . Medan standardrepresentationen av matrisen $G|_{t\times s}$ kräver $O((\#t)(\#s))$ lagringsenheter, Faktoriserad representation kräver endast $O(k(\#t+\#s))$ enheter. Om $k$ inte är för stor, reduceras lagringskraven avsevärt.

För att approximera hela matrisen $G$ är den uppdelad i en familj av submatriser. Stora submatriser lagras i faktoriserad representation, medan små submatriser lagras i standardrepresentation för att förbättra effektiviteten.

Lågrankade matriser är nära besläktade med degenererade expansioner som används i panelklustring och den snabba multipolmetoden för att approximera integraloperatorer. I denna mening kan hierarkiska matriser betraktas som de algebraiska motsvarigheterna till dessa tekniker.

Applikation på integrerade operatörer

Hierarkiska matriser används framgångsrikt för att behandla integralekvationer, t.ex. enkel- och dubbellagers potentiella operatorer som visas i gränselementmetoden . En typisk operatör har formen

${\mathcal {G}}[u](x)=\int _{\Omega }\kappa (x,y)u(y)\,dy.$

Galerkinmetoden leder till matrisinmatningar av formuläret

$g_{ij}=\int _{\Omega }\int _{\Omega }\kappa (x,y)\varphi _{i}(x)\psi _{j}(y)\,dy\,dx,$

där $(\varphi _{i})_{i\in I}$ och $(\psi _{j})_{j\ i J}$ är familjer av finita elementbasfunktioner. Om kärnfunktionen $\kappa$ är tillräckligt jämn, kan vi approximera den genom polynominterpolation för att erhålla

${\tilde {\kappa }}(x,y)=\summa _{\nu = 1}^{k}\kappa (x,\xi _{\nu})\ell _{\nu}(y),$

där $(\xi _{\nu })_{\nu =1}^{k}$ är familjen av interpolationspunkter och ${\ displaystyle (\ell _{\nu })_{\nu =1}^{k}}$ är motsvarande familj av Lagrangepolynom . Att ersätta $\kappa$ med ${\tilde {\kappa }}$ ger en approximation

${\tilde {g}}_{ij}=\int _{\Omega } \int _{\Omega }{\tilde {\kappa }}(x,y)\varphi _{i}(x)\psi _{j}(y)\,dy\,dx=\summa _{\ nu =1}^{k}\int _{\Omega }\kappa (x,\xi _{\nu})\varphi _{i}(x)\,dx\int _{\Omega }\ell _ {\nu }(y)\psi _{j}(y)\,dy=\sum _{\nu =1}^{k}a_{i\nu }b_{j\nu }$

med koefficienterna

$a_{i\nu }=\int _{\Omega }\kappa (x,\xi _{\nu } )\varphi _{i}(x)\,dx,$

$b_{j\nu }=\int _{\Omega }\ell _{\nu}(y)\psi _{j}(y)\,dy.$

Om vi ​​väljer $t\subseteq I,s\subseteq J$ och använder samma interpolationspunkter för alla $i\in t,j\in s$ , får vi $G|_{t\times s}\approx AB^{*}$ .

Uppenbarligen skulle varje annan approximation som separerar variablerna $x$ och ${\displaystyle y},$ t.ex. multipolexpansionen, också tillåta oss att dela dubbelintegralen i två enkla integraler och därmed komma fram till en liknande faktoriserad låg- rangmatris.

Av särskilt intresse är korsapproximationstekniker som endast använder ingångarna i den ursprungliga matrisen $G$ för att konstruera en approximation med låg rangordning .

Tillämpning på elliptiska partiella differentialekvationer

Eftersom lösningsoperatorn för en elliptisk partiell differentialekvation kan uttryckas som en integraloperator som involverar Greens funktion , är det inte förvånande att inversen av styvhetsmatrisen som härrör från finita elementmetoden och spektralmetoden kan approximeras av en hierarkisk matris.

Gröns funktion beror på formen på beräkningsdomänen, därför är den vanligtvis inte känd. Icke desto mindre kan approximativa aritmetiska operationer användas för att beräkna en approximativ invers utan att explicit känna till funktionen.

Överraskande nog är det möjligt att bevisa att inversen kan approximeras även om differentialoperatorn involverar icke-jämna koefficienter och Greens funktion därför inte är jämn.

Aritmetiska operationer

Den viktigaste innovationen med den hierarkiska matrismetoden är utvecklingen av effektiva algoritmer för att utföra (ungefärliga) matrisaritmetiska operationer på icke-glesa matriser, t.ex. för att beräkna approximativa inverser, LU-dekomponeringar och lösningar på matrisekvationer.

Den centrala algoritmen är den effektiva matris-matrismultiplikationen, dvs. beräkningen av $Z=Z+\alpha XY$ för hierarkiska matriser $X,Y,Z$ och en skalär faktor $\alpha$ . Algoritmen kräver att submatriserna för de hierarkiska matriserna organiseras i en blockträdstruktur och drar fördel av egenskaperna hos faktoriserade lågrankade matriser för att beräkna den uppdaterade $Z$ i $O(nk^{2}\,\log(n)^{2})$ operationer.

Genom att dra fördel av blockstrukturen kan inversen beräknas genom att använda rekursion för att beräkna inverser och Schur-komplement av diagonala block och kombinera båda med hjälp av matris-matrismultiplikationen. På ett liknande sätt LU-sönderdelningen konstrueras med endast rekursion och multiplikation. Båda operationerna kräver också $O(nk^{2}\,\log(n)^{2})$ operationer.

H2 ^- matriser

För att behandla mycket stora problem kan strukturen av hierarkiska matriser förbättras: H ² -matriser ersätter blockens allmänna lågrankade struktur med en hierarkisk representation nära relaterad till den snabba multipolmetoden för att reducera lagringskomplexiteten till $O(nk)$ .

I samband med gränsintegraloperatorer leder ersättning av den fasta rankningen $k$ med blockberoende rankningar till approximationer som bevarar konvergenshastigheten för den underliggande gränselementmetoden vid en komplexitet av $O(n).$

Aritmetiska operationer som multiplikation, inversion och Cholesky- eller LR-faktorisering av H ² -matriser kan implementeras baserat på två fundamentala operationer: matris-vektormultiplikationen med submatriser och lågrankningsuppdateringen av submatriser. Även om matris-vektormultiplikationen är enkel, är det en betydande utmaning att implementera effektiva lågrankade uppdateringar med adaptivt optimerade klusterbaser.

programvara

HLib är ett C-programbibliotek som implementerar de viktigaste algoritmerna för hierarkiska och ${\mathcal {H}}^{2}$ -matriser.

AHMED är ett C++-programbibliotek som kan laddas ner för utbildningsändamål.

HLIBpro är en implementering av de kärna hierarkiska matrisalgoritmerna för kommersiella applikationer.

H2Lib är en öppen källkodsimplementering av hierarkiska matrisalgoritmer avsedda för forskning och undervisning.

awesome-hierarchical-matrices är ett arkiv som innehåller en lista över andra H-Matrices-implementationer.

HierarchicalMatrices.jl är ett Julia-paket som implementerar hierarkiska matriser.