# How to accurately decompose positive semidefinite matrix and use the lower triangular part in linear equations

I have $n$ arbitrary $p\times 1$ vectors $x_i$, and $p\times k$ matrices $A_i$, and $n$ $p \times p$ positive semidefinite matrices $S_i$, where some (often most) of the $S_i$'s are same (for example only two different $S$ matrices, one positive definite which applies to $i=1,\ldots n-1$ and semidefinite $S$ for $i=n$). I want to make following linear transformations for all $x_i$ and $A_i$:

$x_i^* = L_i^{-1}x_i$ and $A_i^* = L_i^{-1}A_i$, where

$L_i$ is a lower triangular matrix so that $L_i L'_i=S_i$ (or even better $L_i D_i L'_i=S_i$ where $D_i$ is diagonal matrix.)

$n$ can be range from few to thousands or even more, $p$ is usually less than 10 and $k$ is usually less than 100. $S$ can contain rows and columns with zeros, and it is originally formed as $BB'$ where $B$ is a lower tridiagonal matrix with nonnegative diagonal elements.These $B$'s are not available though.

I'm wondering what would be the optimal way of accomplishing this in terms of speed and accuracy, with more weight on accuracy?

Currently I'm using my self written LDL decomposition function for different $S$, invert $L$ and compute the transformations, as I have beem dealing with cases with small number of different $S$ and large $n$. If I have understood correctly, it would be wiser to just solve the linear equations without explicitly inverting the matrix in terms of accuracy, but for speed it seems to be better option?

I'm using Fortran, and I would rather use some LAPACK function for the factorization instead of my own LDL decomposition which I cannot guarantee to be stable etc., but I'm not sure which would be the good way.

edit: This may be better suited for stackoverflow, so this can be removed by whoever has the power to do so.

-
this page may have the information you need on factorization schemes in LAPACK. As far as which site it should be on, this question is well suited for this one, since this is pretty much exactly what we do. –  Godric Seer Dec 13 '12 at 15:13
If one of your $S_i$'s is truly singular, then you will need to rethink your algorithm, as the inverse does not exist (in fact, one or more diagonal entries of the $D$ in your $LDL^T$ factorization will be zero). Are you sure that the $S_i$'s are not indefinite? –  Jack Poulson Dec 13 '12 at 17:13
Ok, now that I think of it, I don't actually need the zero rows/columns at all, so I can work on a subsets of $S$, $A$ and $x$. But what about semidefinitess of $S$, if $S$ is build as $B'B$, can this result at least numerically semidefinite $S$ instead of definite $S$? $B$ is build as diag(B)=exp(0.5*param), and lower triangular part of $B$ is filled without constraints. –  Hemmo Dec 14 '12 at 6:57
Is $B$ lower-triangular? If so, $B$ is a (non-conventional) Cholesky factor for $S$ and you don't need to perform any factorizations. In addition, it looks like the diagonal of $B$ is guaranteed to be positive, and so $B$ cannot be singular. –  Jack Poulson Dec 20 '12 at 16:35