[R] Modified Cholesky decomposition for sparse matrices

[Michael Braun]

I am trying to estimate a covariance matrix from the Hessian of a posterior mode.  However, this Hessian is indefinite (possibly because of numerical/roundoff issues), and thus, the Cholesky decomposition does not exist.  So, I want to use a modified Cholesky algorithm to estimate a Cholesky of a pseudovariance that is reasonably close to the original matrix.  I know that there are R packages that contain code for Gill-Murray and Schnabel-Eskow algorithms for standard, dense, base-R matrices.  But my Matrix is large (k=30000), and sparse (block-arrow structure, stored as a dsCMatrix class from the Matrix package).  

Is anyone aware of existing code (or perhaps an algorithm that is easy to adapt) that would perform a modified Cholesky decomposition on a large, sparse indefinite matrix, preferably working on sparseMatrix classes?  Alternatively, is there a way I could compute a sparse LDL' decomposition from an existing R function, and quickly modify the output? 

Thanks,

Michael
 


-------------------------------------------
Michael Braun
Associate Professor of Management Science
MIT Sloan School of Management
100 Main St.., E62-535
Cambridge, MA 02139
braunm at mit.edu
617-253-3436