## Posts Tagged ‘factorization’

### LU decomposition using two Cholesky decompositions (part two)

Sunday, March 4th, 2012

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 /
``````

### LU decomposition using two Cholesky factorizations (part one)

Thursday, November 3rd, 2011

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 = LLT`. 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 \
\ CT * /
``````

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 xT A x + xT B
subject to CT x= D
``````

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

``````
[xT λT] / A  C \ /x\ + [xT λT] /  B  \
\ CT * / \λ/           \ -2D /
``````

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

``````
/ A  C \ /x\ = / -B/2 \
\ CT * / \λ/   \  D   /
``````

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

``````
/ A  C \
\ CT * /
``````

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

``````
/ A  C \ = / ?  * \ / ?  ? \
\ CT * /   \ ?  ? / \ *  ? /
``````

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

``````
/ A  C \ = / L  * \ / LT  ? \
\ CT * /   \ ?  ? / \ *   ? /
``````

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  * \ / LT  M \
\ CT * /   \ MT ? / \ *  ? /
``````

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 `-MTM`. But again this is made easy via Cholesky facotrization. `MTM` is clearly symmetric so we can Choleksy factorize it into:

``````
MTM = KKT
``````

where K is lower triangular.

So filling in with `K` we have our finally LU decomposition:

``````
/ A  C \ = / L   * \ / LT  M \
\ CT * /   \ MT -K / \ *  KT /
``````

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 `MTM` to compute `K`.

Update: See my following post about taken advantage of cholesky decomposition which give you permutation matrices.