[SciPy-user] best way of matrix G => L, D: G = L D L^T

[Dominique Orban]

On 8/24/07, dmitrey <openopt at ukr.net> wrote:

> I need best way of getting from matrix G matrices L, D: G = L D L^T (L -
> unit lower triangular matrix, D - diagonal)
> (IIRC it's related to LU-decomposition)

The decomposition that you mention is called a Cholesky factorization.
It is only applicable to *symmetric* and *positive definite* matrices.
Assuming you are not using sparse matrices, you can obtain L and D
from numpy.linalg.cholesky:

In [4]: A = numpy.matrix( [ [4, 2, 1], [2, 4, 1], [1, 1, 4] ], dtype=float )
In [5]: A
Out[5]:
matrix([[ 4.,  2.,  1.],
        [ 2.,  4.,  1.],
        [ 1.,  1.,  4.]])

In [6]: numpy.linalg.cholesky(A)
Out[6]:
matrix([[ 2.        ,  0.        ,  0.        ],
        [ 1.        ,  1.73205081,  0.        ],
        [ 0.5       ,  0.28867513,  1.91485422]])

Here you get L and D combined (there is no need for the function to
return the transpose of L... it is trivial to write a function that
solves triangular systems with L and its transpose). This
decomposition is usually denoted L L^t.

If the decomposition that you want is L D L^t, the L you are looking
for is the above matrix with diagonal elements replaced by 1, and D is
the diagonal matrix whose diagonal elements are the squares of the
diagonal elements in the above matrix.

In SciPy, there is also scipy.linalg.cholesky, scipy.linalg.cho_factor
and scipy.linalg.cho_solve. I don't know if those work with sparse
matrices though.

> (it's for implementing Gill-Murray Stable Newton's Method)
> then I need to solve the equation
>
> L^T d = e[t]
> for d , where e[t] is a unit vector with the t-th component of e[t] being 1
>
> What's the best way of doing it via numpy/scipy?
>
> Has using something like colamd before LU any sense here?

Approximate minimum degree ordering is really relevant for large and
sparse matrices. If that's what you're working with, you'll need to
have a sparse Cholesky hooked up. I don't know if there is one in
SciPy. For instance, Choldmod (from the same author as Umfpack) is a
great code. There are also a couple simpler ones on the same author's
webpage (Tim Davis at University of Florida, Gainseville) called "ldl"
and another one in his "Csparse" library.

Typically, a sparse factorization code will do some ordering
beforehand. However, ldl doesn't (because it is meant to be as concise
as possible).

Dominique