Crout-decompositie

De Crout-decompositie is een algoritme voor de LU-decompositie van een vierkante niet-singuliere matrix in een benedendriehoeksmatrix ${\displaystyle L}$ en een bovendriehoeksmatrix ${\displaystyle U.}$ In de matrix ${\displaystyle U}$ zijn de elementen op de hoofddiagonaal gelijk aan 1.

De methode is genoemd naar Prescott Durand Crout, een wiskundige aan het Massachusetts Institute of Technology die de methode in 1941 heeft beschreven.[1]

Beschrijving

Via Gauss-eliminatie

Zij ${\displaystyle A}$ een ${\displaystyle n\times n}$-matrix met elementen ${\displaystyle a_{i,j}.}$ De Crout-decompositie is formeel gezien een Gauss-eliminatie. Het algoritme begint dan met als ${\displaystyle L}$ de ${\displaystyle n\times n}$-eenheidsmatrix en ${\displaystyle U}$ gelijk aan ${\displaystyle A.}$ In ${\displaystyle n-1}$ stappen wordt ${\displaystyle U}$ herleid tot een bovendriehoeksmatrix en ${\displaystyle L}$ tot een benedendriehoeksmatrix, zodanig dat het matrixproduct ${\displaystyle LU}$ steeds gelijk blijft aan ${\displaystyle A.}$

In de ${\displaystyle k}$-de stap worden het elementen op de hoofddiagonaal in de ${\displaystyle k}$-de kolom van ${\displaystyle U}$ op een gebracht en de elementen eronder op nul. Dit gebeurt door rij ${\displaystyle k}$ te delen door het pivot-element ${\displaystyle u_{k,k}.}$ Deze waarde kopiëren we in de corresponderende positie in matrix ${\displaystyle L.}$ Daarna vermenigvuldigen we de ${\displaystyle k}$-de rij met ${\displaystyle u_{j,k}}$ en trekken ze dan af van rij ${\displaystyle j,}$ wat een nieuwe rij ${\displaystyle j}$ geeft met nul in de ${\displaystyle k}$-de kolom (${\displaystyle j>k}$). De factoren ${\displaystyle u_{j,k}}$ kopiëren we eveneens in de corresponderende kolom van ${\displaystyle L.}$

Indien ${\displaystyle U}$ een element ${\displaystyle u_{k,k}}$ op de hoofddiagonaal heeft dat nul is, gaat dit natuurlijk niet op. In dat geval moet de rij ${\displaystyle k}$ verwisseld worden met een onderliggende rij ${\displaystyle j}$ waarin ${\displaystyle u_{j,k}}$ niet nul is. De verwisseling gebeurt ook in ${\displaystyle L.}$ De verwisselingen kunnen bijgehouden worden in een permutatiematrix ${\displaystyle P.}$ In het begin is ${\displaystyle P}$ de eenheidsmatrix. Wanneer twee rijen in ${\displaystyle U}$ en ${\displaystyle L}$ verwisseld worden, worden dezelfde rijen in ${\displaystyle P}$ verwisseld. De decompositie is dan niet langer ${\displaystyle A=LU}$ maar ${\displaystyle A=PLU.}$

Rechtstreekse berekening

Indien men rekening houdt met de specifieke structuur van ${\displaystyle L}$ en ${\displaystyle U,}$ kan men de coëfficiënten ervan een voor een berekenen. Immers als men het matrixproduct ${\displaystyle A=LU}$ uitschrijft, verkrijgt men een stelsel van ${\displaystyle n\times n}$ lineaire vergelijkingen in de ${\displaystyle n\times n}$ onbekende coëfficiënten van ${\displaystyle L}$ en ${\displaystyle U.}$ Deze vergelijkingen kan men zo rangschikken dat elke volgende vergelijking een nieuwe coëfficiënt geeft; bijvoorbeeld:

${\displaystyle a_{1,1}=l_{1,1},a_{2,1}=l_{2,1}}$ enz. (dit is de eerste kolom van ${\displaystyle L}$)

${\displaystyle a_{1,2}=l_{1,1}u_{1,2}}$ of ${\displaystyle u_{1,2}={\frac {a_{1,2}}{a_{1,1}}}}$

${\displaystyle a_{1,3}=l_{1,1}u_{1,3}}$ of ${\displaystyle u_{1,3}={\frac {a_{1,3}}{a_{1,1}}}}$ enz. Dit maakt de eerste rij van ${\displaystyle U}$ volledig.

${\displaystyle a_{2,2}=l_{2,1}u_{1,2}+l_{2,2}}$ of ${\displaystyle l_{2,2}=a_{2,2}-l_{2,1}u_{1,2}}$

${\displaystyle a_{3,2}=l_{3,1}u_{1,2}+l_{3,2}}$ of ${\displaystyle l_{3,2}=a_{3,2}-l_{3,1}u_{1,2}}$ enz. (dit vervolledigt de tweede kolom van ${\displaystyle L}$)

Op deze manier kan men afwisselend een kolom van ${\displaystyle L}$ en een rij van ${\displaystyle U}$ rechtstreeks berekenen. Dit is de "compacte" vorm van de Crout-decompositie.

Voorbeeld

Gegeven is de matrix:

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

Neem als matrix ${\displaystyle L}$ de eenheidsmatrix en ${\displaystyle U}$ gelijk aan ${\displaystyle A}$.

${\displaystyle LU={\begin{pmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\\\end{pmatrix}}{\begin{pmatrix}3&-7&-2&2\\-3&5&1&0\\6&-4&0&-5\\-9&5&-5&12\\\end{pmatrix}}}$

Eerste stap

Deel de eerste rij van ${\displaystyle U}$ door 3, en vermenigvuldig de eerste rij achtereenvolgens met −3, 6, 9 en trek ze af van de tweede, derde en vierde rij zodat de elementen onder de diagonaal in de eerste kolom van ${\displaystyle U}$ nul worden. De getallen 3, −3, 6 en 9 kopiëren we naar de eerste kolom van ${\displaystyle L,}$ die dus een kopie wordt van de eerste kolom van ${\displaystyle U.}$ Dit geeft als tussenstand:

${\displaystyle LU={\begin{pmatrix}3&0&0&0\\-3&1&0&0\\6&0&1&0\\-9&0&0&1\\\end{pmatrix}}{\begin{pmatrix}1&-7/3&-2/3&2/3\\0&-2&-1&2\\0&10&4&-9\\0&-16&-11&18\\\end{pmatrix}}}$

Men kan nagaan dat het product van deze twee matrices de oorspronkelijk matrix ${\displaystyle A}$ is.

Tweede stap

Om het element op de hoofddiagonaal in de tweede kolom van ${\displaystyle U}$ op een te krijgen, moet de tweede rij van ${\displaystyle U}$ gedeeld worden door −2. Om de rest van de tweede kolom van ${\displaystyle U}$ op nul te krijgen, moet de tweede rij van ${\displaystyle U}$ achtereenvolgens vermenigvuldig worden met 10 en −16, en afgetrokken van de derde, respectievelijk de vierde rij. De waarden −2, 10 en −16 gaan naar de tweede kolom van ${\displaystyle L:}$, die is dus een kopie van de tweede kolom van ${\displaystyle U}$ vanaf de tweede rij. Dit geeft:

${\displaystyle LU={\begin{pmatrix}3&0&0&0\\-3&-2&0&0\\6&10&1&0\\-9&-16&0&1\\\end{pmatrix}}{\begin{pmatrix}1&-7/3&-2/3&2/3\\0&1&1/2&-1\\0&0&-1&1\\0&0&-3&2\\\end{pmatrix}}}$

Derde stap

In de derde stap wordt de derde rij van ${\displaystyle U}$ gedeeld door −1, waarna de derde rij vermenigvuldig wordt met −3 en afgetrokken van de vierde rij. De factoren −1 en −3 gaan naar de derde kolom van ${\displaystyle L}$:

${\displaystyle LU={\begin{pmatrix}3&0&0&0\\-3&-2&0&0\\6&10&-1&0\\-9&-16&-3&1\\\end{pmatrix}}{\begin{pmatrix}1&-7/3&-2/3&2/3\\0&1&1/2&-1\\0&0&1&-1\\0&0&0&-1\\\end{pmatrix}}}$

Vierde stap

In de vierde en laatste stap moet alleen nog het element in de linker benedenhoek van ${\displaystyle U}$ op een gebracht worden, door het te delen door −1, en die waarde te kopiëren naar ${\displaystyle L.}$ Dit geeft als eindresultaat:

${\displaystyle LU={\begin{pmatrix}3&0&0&0\\-3&-2&0&0\\6&10&-1&0\\-9&-16&-3&-1\\\end{pmatrix}}{\begin{pmatrix}1&-7/3&-2/3&2/3\\0&1&1/2&-1\\0&0&1&-1\\0&0&0&1\\\end{pmatrix}}}$

Daarmee is de decompositie van ${\displaystyle A}$ in ${\displaystyle LU}$ beëindigd. In dit geval was er geen verwisseling van rijen nodig.

Opmerkingen

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

1. Vorm de matrices ${\displaystyle L}$ en ${\displaystyle U}$: ${\displaystyle LUx=B}$. We noemen ${\displaystyle y}$ het (voorlopig onbekende) product ${\displaystyle Ux}$.
2. Los ${\displaystyle Ly=B}$ op door voorwaartse substitutie.
3. Los ${\displaystyle Ux=y}$ op door achterwaartse substitutie.

Als bijvoorbeeld ${\displaystyle B={\begin{pmatrix}-36&20&2&-34\end{pmatrix}}^{T},}$ vinden we achtereenvolgens:

• ${\displaystyle 3y_{1}=-36}$ of ${\displaystyle y_{1}=-12}$
• ${\displaystyle -3y_{1}-2y_{2}=20}$ waaruit ${\displaystyle y_{2}=8}$
• ${\displaystyle 6y_{1}+10y_{2}-y_{3}=2}$ waaruit ${\displaystyle y_{3}=6}$
• ${\displaystyle -9y_{1}-16y_{2}-3y_{3}-y_{4}=-34}$ waaruit ${\displaystyle y_{4}=-4}$

en:

• ${\displaystyle x_{4}=y_{4}=-4}$
• ${\displaystyle x_{3}-x_{4}=y_{3}=6}$ waaruit ${\displaystyle x_{3}=2}$
• ${\displaystyle x_{2}+{\tfrac {x_{3}}{2}}-x_{4}=y_{2}=8}$ waaruit ${\displaystyle x_{2}=3}$
• ${\displaystyle x_{1}-{\tfrac {7x_{2}}{3}}-{\tfrac {2x_{3}}{3}}+{\tfrac {2x_{4}}{3}}=y_{1}=-12}$ waaruit ${\displaystyle x_{1}=-1}$

Indien deze methode in een computerprogramma wordt geïmplementeerd, is het niet nodig om twee afzonderlijke matrices voor ${\displaystyle L}$ en ${\displaystyle U}$ te gebruiken. Ze kunnen compact in een enkele matrix bewaard worden (daarvoor kan zelfs de matrix ${\displaystyle A}$ dienen als deze mag overschreden worden). Het is immers niet nodig om de elementen op de hoofddiagonaal van ${\displaystyle U}$ te bewaren, die per definitie gelijk zijn aan 1. De rij-permutaties kunnen bijgehouden worden in een permutatievector van ${\displaystyle n}$ elementen 1 tot en met ${\displaystyle n;}$ als twee rijen verwisseld worden worden de corresponderende elementen in de permutatievector verwisseld.

De Crout-decompositie gelijkt sterk op de Doolittle-decompositie. Dit is ook een LU-decompositie met dit verschil dat in een Doolittle-decompositie de elementen op de hoofddiagonaal van de benedendriehoeksmatrix ${\displaystyle L}$ gelijk zijn aan 1.

Websites

• Cholesky, Doolittle and Crout Factorization

Voetnoten

1. Prescott D. Crout. "A Short Method for Evaluating Determinants and Solving Systems of Linear Equations With Real or Complex Coefficients " (1941), Transactions of the American Institute of Electrical Engineers, vol. 60 nr. 12, blz. 1235–1240. DOI:10.1109/T-AIEE.1941.5058258