Re: [Help-gsl] Recovering from failed Cholesky decomp.

[Brian Gough]

At Wed, 24 Mar 2010 12:38:53 -0600,
Daniel Neilson wrote:
>   I have an application in which I'm using gsl_linalg_cholesky_decomp() 
> to decompose a symmetric positive-definite matrix, A. However, the 
> decomp fails in a very small percentage of the matrices that I feed to 
> this function due to round-off error and all of that fun stuff.
> 
>   Does anyone happen to know of a good method that will still give me 
> the L L^T decomposition for A when one of these errors is detected? For 
> example, is there a method whereby I can add a small-magnitude diagonal 
> matrix to A, Cholesky decomp that, and modify the result to obtain the 
> decomp as though I'd done it with A instead of the modified A?

In linalg/cholesky.c there is a function which handles all the square roots

    double
    quiet_sqrt (double x)  
         /* avoids runtime error, for checking matrix for positive definiteness 
*/
    {
      return (x >= 0) ? sqrt(x) : GSL_NAN;
    }

You make a copy of cholesky.c in your application and modify this to
return 0 if x is negative but "close enough" to zero.

-- 
Brian Gough

GNU Scientific Library -
http://www.gnu.org/software/gsl/