LU decomposition using two Cholesky factorizations (part one)

Alec Jacobson

November 03, 2011

weblog/

My colleague showed me a neat linear algebra trick to compute the LU decomposition of a special type of matrix using two Cholesky factorizations. This is typical much faster. We'd like to factorize Q = LU where L is lower triangular and U is upper triangular. If Q is symmetric positive-definite we can use Cholesky factorization to compute Q = LL^T. Often this unfortunately is not the case and we typically fall back on the slower LU decomposition. But we can use a trick if Q is a special matrix of the following block form:

Q = / A  C \
    \ C^T * /

where A is an n by n symmetric, positive definite matrix, C is a full row rank m by n matrix and * is all zeros (here m by m). These types of systems are actually very common. They arrise often during energy minimization with linear equality constraints. For example I might have the quadratic energy:

minimize x^T A x + x^T B
subject to C^T x= D

Using lagrange multiplers I can transform this minimization into finding the sandle point of

[x^T λ^T] / A  C \ /x\ + [x^T λ^T] /  B  \
        \ C^T * / \λ/           \ -2D /

And to find the sandle point we take derivatives with respect to each x and λ resulting in

/ A  C \ /x\ = / -B/2 \
\ C^T * / \λ/   \  D   /

Then we can solve for x and λ by computing an LU factorization of

/ A  C \
\ C^T * /

Just like our matrix is divided into 4 blocks, we can divide our target L and U into 4 blocks:

/ A  C \ = / ?  * \ / ?  ? \
\ C^T * /   \ ?  ? / \ *  ? /

Since we are assuming that A is symmetric positive definite we can compute a Cholesky factorization, namely let A = LL^T. Then we can place L into the appropriate corners of the lower and triangular parts:

/ A  C \ = / L  * \ / L^T  ? \
\ C^T * /   \ ?  ? / \ *   ? /

But this means that know we need that C = L ?. But L is lower triangular so it's easy to compute:

 
M = L^-1C

We can then place M into our LU decomposition:

/ A  C \ = / L  * \ / L^T  M \
\ C^T * /   \ M^T ? / \ *  ? /

We're almost done. All that's left is reproduce the zeros in the bottom right corner on the left hand side. This means we need the two missing pieces when multiplied to equal -M^TM. But again this is made easy via Cholesky facotrization. M^TM is clearly symmetric so we can Choleksy factorize it into:

M^TM = KK^T

where K is lower triangular. So filling in with K we have our finally LU decomposition:

/ A  C \ = / L   * \ / L^T  M \
\ C^T * /   \ M^T -K / \ *  K^T /

We only needed two cholesky facotrizaions. And further each Cholesky factorization can be completed according to the characteristics of A (our original quadratic coefficients) and C our constraints. So if A is sparse then we can use sparse cholesky factorization to compute L. But then say if C on the other hand is dense, M will come out dense and we can use dense Cholesky factorization on M^TM to compute K. Update: See my following post about taken advantage of cholesky decomposition which give you permutation matrices.