Cholesky-decompositie

Uit Wikipedia, de vrije encyclopedie
Ga naar: navigatie, zoeken

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:

 A = LL^T

waarin L een benedendriehoeksmatrix is. L^T is de getransponeerde matrix van L. L noemt men de Choleskyfactor van 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]

Bepaling van de Cholesky-decompositie[bewerken]

De Cholesky-factor van een n×n-matrix A kan met het volgende algoritme recursief berekend worden in n stappen. Elke stap levert een kolom van L op.

De vergelijking  A = LL^T wordt uitgeschreven als:


\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 L bepaald:

 l_{1,1} = \sqrt{a_{1,1}} (het element op de diagonaal)
L_{2,1} = \frac{A_{2,1}}{l_{1,1}} (dit is de rest van de kolom).

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

A_{2,2} - L_{2,1}L_{2,1}^T = L_{2,2}L_{2,2}^T

(L_{2,1}L_{2,1}^T is het uitwendig product van de kolomvectorL_{2,1} en de rijvector L_{2,1}^T)

Dit is opnieuw een vergelijking van de vorm  A = L L^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 L_{2,2} te vinden, wat de tweede kolom van 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 P bij, zodat de decompositie wordt uitgedrukt als:

PAP^T = LL^T

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

Voorbeeld[bewerken]

We zoeken de decompositie van de matrix

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

Eerste stap[bewerken]


\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 L wordt:

l_{1,1} = \sqrt 4 = 2
l_{2,1} = 12/2 = 6
l_{3,1} = -16/2 = -8

De 2x2-matrix voor de volgende stap is


\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[bewerken]

Analoog aan de eerste stap, levert de tweede stap:

l_{2,2} = \sqrt 1 = 1
l_{3,2} = 5


Blijft over de "matrixvergelijking" met 1x1-matrices:


\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[bewerken]

Uit de vorige vergelijking volgt eenvoudigweg dat

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

De Cholesky-decompositie van A is dus

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[bewerken]

De symmetrische vorm van de Cholesky-decompositie is:

 A = L_1DL_1^T

Hierin is L_1 een benedendriehoeksmatrix met 1 op de hoofddiagonaal, en D is een diagonaalmatrix.

Uit

LL^T = L_1DL_1^T

volgt

 L = L_1D^{1/2}.

De elementen in D zijn de kwadraten van de elementen op de hoofddiagonaal van L. De kolommen van L_1 zijn gelijk aan de kolommen van L gedeeld door de elementen op de hoofddiagonaal van L.

In het bovenstaande voorbeeld wordt dit

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}
     4 &  0 & 0 \\
    0   &  1 &  0 \\
    0  &  0  &  9 \\
  \end{pmatrix}

Dus

 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[bewerken]

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

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

Inverse matrix[bewerken]

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

De inverse matrix van A = LL^T wordt dan gegeven door:

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

Grafentheoretische interpretatie[bewerken]

Elke symmetrische vierkante n-op-n-matrix A kan beschouwd worden als de bogenmatrix van een graaf G met n knopen: er is een kant tussen knopen i en j als het matrixelement aij verschilt van nul.

Merk op: de elementen op de diagonaal van 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 i zijn de knopen die door een kant verbonden zijn met knoop i.

De Cholesky-decompositie kan men nu beschrijven als een opeenvolging van eliminaties van een knoop uit G zodat men een reeks eliminatiegrafen bekomt met telkens een knoop minder:

G = G0, G1, G2, ..., Gn-1

De i-de stap in de decompositie correspondeert met de afleiding van Gi uit Gi-1 door verwijdering van de knoop i en van de kanten die in die knoop samenkomen. Indien nodig moeten nieuwe kanten worden toegevoegd zodat de buren van knoop i paarsgewijs verbonden zijn. Deze nieuwe kanten noemt men "opvulling" (fill-in).

De "nul/niet-nul"-structuur van de Choleskyfactor L kan men uit de opeenvolgende eliminatiegrafen afleiden: de elementen in de i-de kolom van L die verschillen van nul, staan op de rijen met de nummers van de buren van knoop i in graaf Gi-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 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 A geschikt te rangschikken vooraleer de Cholesky-decompositie te maken, zoals het volgende voorbeeld duidelijk maakt:

Voorbeeld[bewerken]

Alternatief 1[bewerken]

Stel de matrix A heeft de volgende nul/niet-nul-structuur

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:

GraphA1.png

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

GraphA1cut.png

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

GraphA2.png

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

GraphA3.png

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

GraphA4.png

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

GraphA5.png

en tenslotte

GraphA6.png

De nul/niet-nul-structuur van de Choleskyfactor van A is bijgevolg:

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" kanten toegevoegd en er zijn zeven nul-elementen in L beneden de diagonaal.

Alternatief 2[bewerken]

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

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 A1 is afgeleid uit A door rijen en kolommen 3 te verwisselen met 4, en 2 met 6. De corresponderende graaf is:

GraphB1.png

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

GraphB1cut.png

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

GraphB2.png

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

GraphB3.png

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

GraphB4.png

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

GraphA5.png

en tenslotte

GraphA6.png

De nul/niet-nul-structuur van de Choleskyfactor van A1 is:

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"-kant toegevoegd en er zijn negen nul-elementen in 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 minimum aantal kanten dat men aan een graaf moet toevoegen om er een chordale graaf van te maken.

De graaf 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.

Zie ook[bewerken]

Bronnen, noten en/of referenties