Numerieke Wiskunde II:
LU-decompositie van een matrix

Het oplossen van een stelsels vergelijkingen met twee of drie variabelen is, met behulp van de Regel van Kramer, relatief eenvoudig. Met de berekening van drie of vier determinanten en twee of drie delingen heb je de oplossing.
Als het stelsel vier of (veel) meer variabelen heeft, wordt het een ander verhaal. Het berekenen van de determinanten wordt lastiger naarmate er meer onbekenden zijn. In de huidige stand van de techniek (lees: toepassing van numerieke wiskunde op computers) zijn stelsels met een paar honderdduizend onbekenden niet uitzonderlijk.
Het mag duidelijk zijn dat de Regel van Kramer hiervoor niet geschikt is.

Gauss-eliminatie zou een alternatief kunnen zijn, maar alleen voor niet te grote stelsels vergelijkingen. Persoonlijk vind ik dat zes vergelijkingen wel het maximum is voor een oplossing op papier. In een computerprogramma heb je in principe geen beperkingen, maar Gauss-eliminatie kan erg veel rekenkracht vergen en daardoor traag zijn. Aan Gauss-eliminatie kleeft het nadeel dat een efficiënte aanpak creativiteit vereist. Dat is lastig te programmeren.
Bij Gauss-eliminatie zoek je vergelijkingen die veel op elkaar lijken. Die trek je, na vermenigvuldiging met een factor, van elkaar af, net zolang tot er een boven- of onderdriehoeksmatrix ontstaat. Die is gemakkelijk oplosbaar.

Een stelsel dat kan worden beschreven met een bovendriehoeksmatrix of een onderdriehoeksmatrix is eenvoudig oplosbaar, zie onderstaand voorbeeld:
 
Beschouw het stelsel:

y1	=	a11x1	+	a12x2	+	a13x3
y2	=			a22x2	+	a23x3
y3	=					a33x3

De matrix-coëfficiënten aij en de vector in het linker lid y zijn bekend. In het rechter lid is de vector x de onbekende die moet worden opgelost.
 
In matrix-vector vorm wordt het stelsel geschreven als:

y1	=	a11	a12	a13	·	x1
y2		0	a22	a23		x2
y3		0	0	a33		x3

Eenvoudig is in te zien dat:
x3 = y3 / a33
x2 = (y2 - a23x3) / a22
x1 = (y1 - a12x2 - a13x3) / a11
Omdat dit een bovendriehoeksmatrix is, wordt er van beneden naar boven opgelost. De oplossing van een onderdriehoeksmatrix gaat op dezelfde manier, maar dan van boven naar beneden.
Merk op dat er steeds wordt gedeeld door de elementen van de hoofddiagonaal. Dat komt later nog van pas.

De oplossing met LU-decompositie
De gedachte is om de systeemmatrix [S] op te splitsen in een onderdriehoeksmatrix [L] (van het Engelse woord 'lower') en een bovendriehoeksmatrix [U] (van 'upper'), die elk eenvoudig oplosbaar zijn. De stappen zijn:
 
Stel: y = [S]·x en [S] = [L]·[U].
 
Dat geeft: y = [L]·[U]·x.
 
Stel: δ = [U]·x. Dit wordt 'substitution' genoemd.
 
Los δ op uit y = [L]·δ.
 
Los x op uit δ = [U]·x. Dit wordt 'backsubstitution' genoemd.

De aanpak bij het splitsen van [S] in [L]·[U]
We gaan uit van:

s11	s12	s13	=	l11	l12	l13	·	u11	u12	u13
s21	s22	s23		l21	l22	l23		u21	u22	u23
s31	s32	s33		l31	l32	l33		u31	u32	u33

Om van [L] een lower matrix te maken moeten l12, l13 en l23 gelijk zijn aan nul.
Om van [U] een upper matrix te maken moeten u21, u31 en u32 gelijk zijn aan nul.
Er ontstaat dan:

s11	s12	s13	=	l11	0	0	·	u11	u12	u13
s21	s22	s23		l21	l22	0		0	u22	u23
s31	s32	s33		l31	l32	l33		0	0	u33

Eenvoudig is te zien dat s11 = l11· u11.

Een axioma in de rekenkunde is dat elk getal is te schrijven als het product van twee andere getallen. Als je er één kiest, ligt de andere vast. Bij de LU-decompositie nemen we voor de getallen op de hoofddiagonaal van [L} altijd 1. We hebben dus nu:
u11 = s11. Zie onderstaand.

s11	s12	s13	=	1	0	0	·	s11	u12	u13
s21	s22	s23		l21	1	0		0	u22	u23
s31	s32	s33		l31	l32	1		0	0	u33

Als we nu verder gaan met de eerste kolom, dan zien we dat:
l21 = s21/s11, en dat
l31 = s31/s11. We hebben nu:

s11	s12	s13	=	1	0	0	·	s11	u12	u13
s21	s22	s23		s21/s11	1	0		0	u22	u23
s31	s32	s33		s31/s11	l32	1		0	0	u33

Nu gaan we de eerste (bovenste) rij afmaken. We zien dat:
u12 = s12 en u13 = s13. Dit geeft:

s11	s12	s13	=	1	0	0	·	s11	s12	s13
s21	s22	s23		s21/s11	1	0		0	u22	u23
s31	s32	s33		s31/s11	l32	1		0	0	u33

Nu lossen we de tweede kolom verder op:
s22 = s12· s21/s11 + u22. Hieruit volgt u22.
s32 = s12· s31/s11 + l32· u22. Hieruit volgt l32.

Nu zijn alleen u23 en u33 nog niet bekend. Die lossen we op uit:
s23 = s13· s21/s11 + u23. Hieruit volgt u23.
Tenslotte: s33 = s13· s31/s11 + u23·l32 + u33. Hieruit volgt u33.

Een uitgewerkt voorbeeld

Decomponeer de matrix [S] =	1	2	3	=	1	0	0	·	u11	u12	u13
	4	5	6		l21	1	0		0	u22	u23
	7	8	9		l31	l32	1		0	0	u33

Oplossing:

• u11 = s11, dus u11 = 1.
• l21 = s21/s11, dus l21 = 4 / 1 = 4.
• l31 = s31/s11, dus l31 = 7 / 1 = 7.
• u12 = s12, dus u12 = 2.
• u13 = s13, dus u13 = 3.
• u22 = s22 - s12· s21/s11, dus u22 = 5 - 2×4/1 = -3.
• l32 = (s32 - s12· s31/s11) / u22, dus l32 = (8 - 2×7/1) / -3 = -6 / -3 = 2.
• u23 = s23 - s13· s21/s11 = 6 - 3×4/1 = -6.
• u33 = s33 - s13· s31/s11 - u23· l32 = 9 - 3×7/1 - 2×-6 = 9 - 21 + 12 = 0.

We hebben nu:

[S] =	1	2	3	=	1	0	0	·	1	2	3
	4	5	6		4	1	0		0	-3	-6
	7	8	9		7	2	1		0	0	0

Bespreking:

• Als je de determinant van de matrix [S] zou uitrekenen, komt daar nul uit. Dat betekent dat [S] singulier is. Dat blijkt uit de oplossing van [U]: Daar is de determinant ook nul van (alleen nullen op de onderste rij). Je kunt de determinant van [U] ook berekenen door de elementen op de hoofddiagonaal met elkaar te vermenigvuldigen. Als er een nul in de hoofddiagonaal van [U] staat is het stelsel [S] niet oplosbaar.
• Merk op dat een slecht geconditioneerde matrix (determinant is bijna nul) of een singuliere matrix (determinant is gelijk aan nul) probleemloos kan worden ontbonden met LU-decompositie. Dat is in lijn met de eigenschap van matrices dat elke matrix kan worden geschreven als het product van twee andere matrices.
• De determinant van [L] is gelijk aan 1. Dat betekent dat δ altijd kan worden opgelost uit y = [L]·δ.
• De determinant van [S] is gelijk aan de determinant van [L] maal de determinant van [U]. Omdat de determinant van [L] gelijk is aan 1, zijn de determinanten van [S] en [U] gelijk aan elkaar.
• Als de matrix [S] problemen heeft, komen die pas te voorschijn bij de backsubstitution. Grote rekenpakketten kijken, voordat er überhaupt met oplossen van [U] en [L] wordt begonnen, even de [U]-matrix na. Dan wordt ook gekeken naar de verhouding tussen de absolute waarden van de getallen tussen de elementen op de hoofddiagonaal van [S] en de overeenkomende elementen op de hoofddiagonaal van [U]. Als die verhouding te groot is, ontstaat het risico van underflow en wordt de oplossing onnauwkeurig.

Om wat te experimenteren met LU-decompositie is een applicatie gemaakt, klik op de knop hiernaast.
De applicatie kan matrices ontbinden met afmetingen van minimaal 2 × 2. Er is geen restrictie aan de grootte van de matrices, behalve dan het werkgeheugen en de rekenkracht van de computer. De applicatie berekent ook de determinant van [S].

De applicatie heeft een debug-mode. Als je in de file ludecomp.js de waarde van debugAndCheck op true zet, wordt de oorspronkelijke matrix [S] terug gerekend uit [L]·[U]. Het resultaat verschijnt in het console. De rijen van de matrix worden achter elkaar gezet als een enkele reeks getallen. Voor een 3 × 3-matrix zijn er dus 9 getallen.

De code van de applicatie kun je downloaden om zelf aan door te ontwikkelen.

Bron:
Boek:
Kammer, ir. R; Numerieke methoden voor Technici, § 8.4 en 8.5. 2^e druk 1977. Uitg. Agon-Elsevier, ISBN 90 10 01840 7.

Downloaden:
 
Druk op de knop:  File: voorb650.zip, 3722 bytes.

Opmerking:
 
De hier beschreven aanpak voor LU-decompositie is bekend als de methode van Doolittle. Hierbij zijn alle elementen op de hoofddiagonaal van [L] gelijk aan 1. Er is ook de –minder bekende– methode van Crout, waar de enen juist op de hoofddiagonaal van [U] staan. De methode van Cholesky is speciaal voor de decompositie van symmetrische matrices, waarbij de matrix-elementen aij gelijk zijn aan de elementen aji.
Symmetrische matrices komen veel voor in technisch-wetenschappelijk rekenwerk aan o.a. staalconstructies.