LU decomposition using two Cholesky decompositions (part two)

Alec Jacobson

March 04, 2012

weblog/

I previously posted about how often you can compute a LU decomposition of a matrix Q using two Cholesky decompositions. Specifically when your matrix is of the form: 
Q = / A  C \
    \ CT * /
In that post I assumed you had a subroutine that computes the Cholesky factorization of a symmetric positive-definite matrix A: 
LLT=A
Where L is a lower-triangular matrix. Often, routines exist to give you an even better factorization which includes a permutation matrix S: 
LLT=STAS
This often results in a more sparse L. Working through the LU decomposition using this type of Cholesky factorization is a bit tedious and requires some careful book keeping. This time I want to fill in the missing parts of the following system: 
/ ? * \ / A  C \ / ? * \ = / ? * \ / ? ? \
\ * ? / \ CT * / \ * ? /   \ ? ? / \ * ? /
Since A is symmetric positive definite we use our factorization STAS = LLT which gives us: 
/ ST * \ / A  C \ / S * \ = / L * \ / LT ? \
\ *  ? / \ CT * / \ * ? /   \ ? ? / \ *  ? /
But this means that know we need that STC = L ?. But L is lower triangular so it’s easy to compute: 
M = L-1STC
We can then place M into our LU decomposition: 
/ ST * \ / A  C \ / S * \ = / L  * \ / LT M \
\ *  ? / \ CT * / \ * ? /   \ MT ? / \ *  ? /
Now at this point we could use cholesky factorization without the permutation matrix to complete the missing parts, but we might as well use our better decomposition again. To do this we must be a bit careful. We want to create the zeros in the lower right corner of the system. To do this we need a factorization of MTM. But our new routine gives us: 
KKT = TTMTMT
where K another lower triangular matrix and T another permutation matrix. So we replace M with MT in our system and plug in K. This gives us: 
/ ST * \ / A  C \ / S * \ = / L    * \ / LT  MT \
\ *  ? / \ CT * / \ * ? /   \ TTMT K / \ *  -KT /
Finally, compute the missing pieces by evaluating the right hand side and we indeed get our original system conjugated with a permutation matrix composed of the permutation matrices from both of the cholesky factorizations: 
/ ST *  \ / A  C \ / S * \ = / L    * \ / LT  MT \
\ *  TT / \ CT * / \ * T /   \ TTMT K / \ *  -KT /