help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Brian Gough
Subject: Re: [Help-gsl] Recovering from failed Cholesky decomp.
Date: Tue, 06 Apr 2010 10:38:51 +0100
User-agent: Wanderlust/2.14.0 (Africa) Emacs/22.2 Mule/5.0 (SAKAKI)

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/




reply via email to

[Prev in Thread] Current Thread [Next in Thread]