Cholesky-decompositie

De Cholesky-decompositie van een positief-definiete Hermitische matrix, of in het reële geval, een positief-definiete symmetrische matrix, is een LU-decompositie van de vorm:

${\displaystyle A=LL^{T}}$

waarin ${\displaystyle L}$ een benedendriehoeksmatrix is. ${\displaystyle L^{T}}$ is de getransponeerde matrix van ${\displaystyle L}$. ${\displaystyle L}$ noemt men de Choleskyfactor van ${\displaystyle A}$.

De Cholesky-decompositie is genoemd naar de Franse militaire officier en wiskundige André-Louis Cholesky (1875-1918), die kort voor het einde van de Eerste Wereldoorlog sneuvelde. Het is niet exact bekend wanneer Cholesky zijn methode bedacht. Hij publiceerde ze zelf niet; ze werd wel indirect bekend dankzij een artikel van commandant Benoît in het Bulletin géodesique van 1924, waarin hij het "procédé du commandant Cholesky" beschreef. Later is een manuscript van Cholesky uit 1910 gevonden waarin hij zijn methode gedetailleerd beschrijft, onder de titel "Sur la résolution numérique des Systèmes d'équations linéaires."^[1]

Methode

De Cholesky-factor van een ${\displaystyle n\times n}$-matrix ${\displaystyle A}$ kan met het volgende algoritme recursief berekend worden in ${\displaystyle n}$ stappen. Elke stap levert een kolom van ${\displaystyle L}$ op.

De vergelijking ${\displaystyle A=LL^{T}}$ wordt uitgeschreven als:

${\displaystyle {\begin{pmatrix}a_{1,1}&A_{2,1}^{T}\\A_{2,1}&A_{2,2}\end{pmatrix}}={\begin{pmatrix}l_{1,1}&0\\L_{2,1}&L_{2,2}\end{pmatrix}}{\begin{pmatrix}l_{1,1}&L_{2,1}^{T}\\0&L_{2,2}^{T}\end{pmatrix}}={\begin{pmatrix}l_{1,1}^{2}&l_{1,1}L_{2,1}^{T}\\l_{1,1}L_{2,1}&L_{2,2}L_{2,2}^{T}+L_{2,1}L_{2,1}^{T}\end{pmatrix}}}$

Hieruit wordt de eerste kolom van ${\displaystyle L}$ bepaald:

${\displaystyle l_{1,1}={\sqrt {a_{1,1}}}}$ (het element op de diagonaal)

${\displaystyle L_{2,1}={\frac {A_{2,1}}{l_{1,1}}}}$ (dit is de rest van de kolom).

Rest nog de onbekende matrix ${\displaystyle L_{2,2}}$, die volgt uit de vergelijking:

${\displaystyle A_{2,2}-L_{2,1}L_{2,1}^{T}=L_{2,2}L_{2,2}^{T}}$

(${\displaystyle L_{2,1}L_{2,1}^{T}}$ is het kruisproduct van de kolomvector${\displaystyle L_{2,1}}$ en de rijvector ${\displaystyle L_{2,1}^{T}}$)

Dit is opnieuw een vergelijking van de vorm ${\displaystyle A=LL^{T}}$, maar nu met een matrix die een kolom en een rij minder telt. Het bovenstaande kan nu herhaald worden om de eerste kolom van ${\displaystyle L_{2,2}}$ te vinden, wat de tweede kolom van ${\displaystyle L}$ is, en zo verder tot er een 1x1-matrix overblijft.

De methode van Cholesky is numeriek stabiel. In de Cholesky-decompositie hoeft men geen rijen of kolommen te verwisselen om te voorkomen dat men een pivotgetal nul bekomt. Maar de volgorde waarin de kolommen en rijen worden geëlimineerd, kan wel een grote invloed hebben op de looptijd van het algoritme. Dit is zeker zo als de matrix schaarsbezet is (een ijle matrix met veel nullen). Om de symmetrie te behouden moet men in de Cholesky-methode steeds rijen samen met de corresponderende kolommen verwisselen. Hiervoor houdt men een permutatiematrix ${\displaystyle P}$ bij, zodat de decompositie wordt uitgedrukt als:

${\displaystyle PAP^{T}=LL^{T}}$

Door een zorgvuldige keuze van de eliminatievolgorde kan men een Cholesky-factor ${\displaystyle L}$ bekomen die zeer schaarsbezet is en snel berekend kan worden (zie verder bij Grafentheoretische interpretatie).

Voorbeeld

We zoeken de decompositie van de matrix

${\displaystyle A={\begin{pmatrix}4&12&-16\\12&37&-43\\-16&-43&98\end{pmatrix}}}$

Eerste stap

${\displaystyle {\begin{pmatrix}4&12&-16\\12&37&-43\\-16&-43&98\end{pmatrix}}={\begin{pmatrix}l_{1,1}&0&0\\l_{2,1}&l_{2,2}&0\\l_{3,1}&l_{3,2}&l_{3,3}\end{pmatrix}}{\begin{pmatrix}l_{1,1}&l_{2,1}&l_{3,1}\\0&l_{2,2}&l_{3,2}\\0&0&l_{3,3}\end{pmatrix}}}$

De eerste kolom van ${\displaystyle L}$ wordt:

${\displaystyle l_{1,1}={\sqrt {4}}=2}$

${\displaystyle l_{2,1}=12/2=6}$

${\displaystyle l_{3,1}=-16/2=-8}$

De 2x2-matrix voor de volgende stap is

${\displaystyle {\begin{pmatrix}37&-43\\-43&98\end{pmatrix}}-{\begin{pmatrix}6\\-8\end{pmatrix}}{\begin{pmatrix}6&-8\end{pmatrix}}={\begin{pmatrix}1&5\\5&34\end{pmatrix}}}$

Tweede stap

De tweede stap gaat op dezelfde manier als de eerste stap:

${\displaystyle l_{2,2}={\sqrt {1}}=1}$

${\displaystyle l_{3,2}=5}$

Blijft over de matrixvergelijking met 1×1-matrices:

${\displaystyle {\begin{pmatrix}34\end{pmatrix}}={\begin{pmatrix}5&l_{3,3}\end{pmatrix}}{\begin{pmatrix}5\\l_{3,3}\end{pmatrix}}={\begin{pmatrix}25+l_{3,3}^{2}\end{pmatrix}}}$

Derde en laatste stap

Uit de vorige vergelijking volgt eenvoudigweg dat

${\displaystyle l_{3,3}={\sqrt {9}}=3}$

De Cholesky-decompositie van ${\displaystyle A}$ is dus

${\displaystyle A={\begin{pmatrix}2&0&0\\6&1&0\\-8&5&3\\\end{pmatrix}}{\begin{pmatrix}2&6&-8\\0&1&5\\0&0&3\\\end{pmatrix}}}$

Symmetrische vorm

De symmetrische vorm van de Cholesky-decompositie is:

${\displaystyle A=L_{1}DL_{1}^{T}}$

Hierin is ${\displaystyle L_{1}}$ een benedendriehoeksmatrix met 1 op de hoofddiagonaal, en ${\displaystyle D}$ is een diagonaalmatrix.

Uit

${\displaystyle LL^{T}=L_{1}DL_{1}^{T}}$

volgt

${\displaystyle L=L_{1}D^{1/2}}$.

De elementen in ${\displaystyle D}$ zijn de kwadraten van de elementen op de hoofddiagonaal van ${\displaystyle L}$. De kolommen van ${\displaystyle L_{1}}$ zijn gelijk aan de kolommen van ${\displaystyle L}$ gedeeld door de elementen op de hoofddiagonaal van ${\displaystyle L}$.

In het bovenstaande voorbeeld wordt dit

${\displaystyle L={\begin{pmatrix}2&0&0\\6&1&0\\-8&5&3\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\3&1&0\\-4&5&1\\\end{pmatrix}}{\begin{pmatrix}2&0&0\\0&1&0\\0&0&3\\\end{pmatrix}}}$

Dus

${\displaystyle A={\begin{pmatrix}1&0&0\\3&1&0\\-4&5&1\\\end{pmatrix}}{\begin{pmatrix}4&0&0\\0&1&0\\0&0&9\\\end{pmatrix}}{\begin{pmatrix}1&3&-4\\0&1&5\\0&0&1\\\end{pmatrix}}}$

Oplossen van een stelsel vergelijkingen

Om een stelsel van lineaire vergelijkingen ${\displaystyle Ax=B}$ op te lossen via Cholesky-decompositie gaat men als volgt te werk:

1. Bereken de matrix ${\displaystyle L}$ en los ${\displaystyle LL^{T}x=B}$ op in twee stappen:
2. los eerst ${\displaystyle Lz=B}$ op met voorwaartse substitutie;
3. los dan ${\displaystyle L^{T}x=z}$ op met achterwaartse substitutie.

Inverse matrix

De diagonaalelementen van de driehoeksmatrix ${\displaystyle L}$ zijn verschillend van nul. Dat betekent dat ${\displaystyle L}$ inverteerbaar is. De determinant ervan is gelijk aan het product van de elementen op de hoofddiagonaal.

De inverse matrix van ${\displaystyle A=LL^{T}}$ wordt dan gegeven door:

${\displaystyle (LL^{T})^{-1}=(L^{T})^{-1}L^{-1}}$

De inverse matrix van een benedendriehoeksmatrix is ook een benedendriehoeksmatrix en kan snel, element per element, berekend worden met een recursieformule.

Interpretatie met grafen

Elke symmetrische vierkante ${\displaystyle n\times n}$-matrix ${\displaystyle A}$ kan beschouwd worden als de bogenmatrix van een graaf ${\displaystyle G}$ met ${\displaystyle n}$ knopen: er is een zijde tussen de knopen ${\displaystyle i}$ en ${\displaystyle j}$ als het matrixelement ${\displaystyle a_{ij}}$ verschilt van nul.

Merk op: de elementen op de diagonaal van ${\displaystyle A}$ komen overeen met een lus in elke knoop; deze worden verder niet meer beschouwd (ze zijn toch steeds verschillend van nul).

De buren van knoop ${\displaystyle i}$ zijn de knopen die door een zijde verbonden zijn met ${\displaystyle i.}$

De Cholesky-decompositie kan men nu beschrijven als een opeenvolging van eliminaties van een knoop uit ${\displaystyle G}$ zodat een reeks eliminatiegrafen ontstaat met telkens een knoop minder:

${\displaystyle G=G_{0},G_{1},G_{2},\ldots ,G_{n-1}}$

De ${\displaystyle i}$-de stap in de decompositie correspondeert met de afleiding van ${\displaystyle G_{i}}$ uit ${\displaystyle G_{i-1}}$ door verwijdering van de knoop ${\displaystyle i}$ en van de zijden die in die knoop samenkomen. Indien nodig moeten nieuwe zijden worden toegevoegd, zodat de buren van knoop ${\displaystyle i}$ paarsgewijs verbonden zijn. Deze nieuwe zijden noemt men "opvulling" (fill-in).

De nul/niet-nul-structuur van de Choleskyfactor ${\displaystyle L}$ kan men uit de opeenvolgende eliminatiegrafen afleiden: de elementen in de ${\displaystyle i}$-de kolom van ${\displaystyle L}$ die verschillen van nul, staan op de rijen met de nummers van de buren van knoop ${\displaystyle i}$ in de graaf ${\displaystyle G_{i-1}.}$ Het element op de diagonaal komt daar nog bij.

Men wenst de hoeveelheid opvulling zo laag mogelijk te houden en het aantal nulelementen in ${\displaystyle L}$ zo hoog mogelijk. Dat vermindert het latere rekenwerk en dit is vooral belangrijk voor zeer grote maar schaarsbezette ijle matrices, die bijvoorbeeld in de eindige-elementenmethode voorkomen. Men kan dit beïnvloeden door de rijen en kolommen in de matrix ${\displaystyle A}$ geschikt te rangschikken vooraleer de Cholesky-decompositie te maken, zoals het volgende voorbeeld duidelijk maakt:

Voorbeeld

Alternatief 1

Stel de matrix ${\displaystyle A}$ heeft de volgende nul/niet-nul-structuur

${\displaystyle A={\begin{pmatrix}\ *&0&*&*&0&0\\\ 0&*&0&0&*&*\\\ *&0&*&0&*&0\\\ *&0&0&*&0&0\\\ 0&*&*&0&*&0\\\ 0&*&0&0&0&*\end{pmatrix}}}$

De corresponderende graaf is:

In de eerste eliminatiestap verwijderen we knoop 1 en de zijden 1-3 en 1-4:

Knoop 1 heeft als buren 3 en 4 en omdat deze niet verbonden zijn, moeten we de graaf opvullen met een nieuwe zijde 3-4. Dit geeft de volgende eliminatiegraaf:

Knoop 2 heeft buren 5 en 6; dit betekent dat in de tweede kolom van de Choleskyfactor ${\displaystyle L}$ de elementen in rijen 5 en 6 niet nul zijn (naast die op de diagonaal). Knoop 2 en kde zijden 2-6 en 2-5 worden geëlimineerd en er moet een nieuwe zijde 5-6 worden toegevoegd, wat als resultaat oplevert:

Knoop 3 heeft buren 4 en 5. De eliminatie van knoop 3 vereist opnieuw de opvulling met een zijde 4-5:

Knoop 4 heeft als enige buur 5; er is geen opvulling nodig. De volgende eliminatiestappen zijn eenvoudig:

en ten slotte

De nul/niet-nul-structuur van de Choleskyfactor van ${\displaystyle A}$ is bijgevolg:

${\displaystyle L={\begin{pmatrix}\ *&0&0&0&0&0\\\ 0&*&0&0&0&0\\\ *&0&*&0&0&0\\\ *&0&*&*&0&0\\\ 0&*&*&*&*&0\\\ 0&*&0&0&*&*\end{pmatrix}}}$

Er zijn in totaal 3 "fill-in" zijden toegevoegd en er zijn zeven nul-elementen in ${\displaystyle L}$ beneden de diagonaal.

Alternatief 2

Neem ter vergelijking volgende matrix met als nul/niet-nul-structuur:

${\displaystyle A_{1}={\begin{pmatrix}\ *&0&*&*&0&0\\\ 0&*&0&0&0&*\\\ *&0&*&0&0&0\\\ *&0&0&*&*&0\\\ 0&0&0&*&*&*\\\ 0&*&0&0&*&*\end{pmatrix}}}$

Matrix ${\displaystyle A_{1}}$ is afgeleid uit ${\displaystyle A}$' door rijen en kolommen 3 te verwisselen met 4, en 2 met 6. De corresponderende graaf is:

De eerste eliminatiestap is dezelfde als in alternatief 1. We verwijderen knoop 1 en de zijden 1-3 en 1-4:

Hier moeten we ook een nieuwe zijde 3-4 toevoegen. Dit geeft de volgende eliminatiegraaf:

Knoop 2 heeft een enkele buur, 6. De eliminatie van knoop 2 vereist dus geen opvulling:

Hetzelfde geldt voor knoop 3, met als enige buur 4:

De eliminatie van knopen 4 en 5 verloopt verder analoog als in alternatief 1, zonder dat er opvulling nodig is:

en ten slotte

De nul/niet-nul-structuur van de Choleskyfactor van ${\displaystyle A_{1}}$ is:

${\displaystyle L_{1}={\begin{pmatrix}\ *&0&0&0&0&0\\\ 0&*&0&0&0&0\\\ *&0&*&0&0&0\\\ *&0&*&*&0&0\\\ 0&0&0&*&*&0\\\ 0&*&0&0&*&*\end{pmatrix}}}$

Er is in dit geval slechts 1 "fill-in"-zijde toegevoegd en er zijn negen nul-elementen in ${\displaystyle L,}$ tegenover zeven in het eerste alternatief.

De grafentheoretische analyse van de Cholesky-decompositie werd voor het eerst gedetailleerd uitgevoerd door Donald J. Rose in zijn doctoraatsthesis aan de Harvard-universiteit uit 1970.^[2] Het bepalen van een optimale volgorde van de vergelijkingen om de fill-in te minimaliseren (het minimum-fill-problem) is een NP-volledig probleem^[3] waarvoor verschillende algoritmen ontwikkeld zijn.

Men spreekt over een perfecte eliminatie-volgorde wanneer er geen enkele nieuwe kant moet toegevoegd worden in de eliminatie. Een graaf heeft een perfecte eliminatie-volgorde als en slechts als het een chordale graaf is (dit is een getriangulariseerde graaf). In een chordale graaf heeft elke cyclus van lengte 4 of meer een koorde, dit is een kant tussen twee niet-opeenvolgende knopen in de cyclus. Het minimum fill-in problem is equivalent aan het vinden van het minimumaantal kanten dat men aan een graaf moet toevoegen om er een chordale graaf van te maken.

De matrix ${\displaystyle A}$ uit het bovenstaande voorbeeld en zijn bijhorende graaf konden zonder opvulling en met een perfecte eliminatie-volgorde verwerkt door in plaats van de volgorde 1-2-3-4-5-6, de volgorde 4-1-3-5-2-6 te nemen.

Voetnoten ↑ Claude Brezinski. "La méthode de Cholesky." Revue d'Histoire des Mathématiques (2005), vol. 11 nr. 2, blz. 205-238. Gearchiveerd op 17 maart 2014. ↑ Vermeld in: J. Alan George. "Computer implementation of the finite element method." Report STAN-CS-71-208, Stanford University, Februari 1971. ↑ Mihalis Yannakakis. "Computing the minimum fill-in is NP-complete." SIAM Journal on Algebraic Discrete Methods (1981), vol. 2 nr. 1, blz. 77-79.