This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Bug in Levenbeg-Marquardt Algorithm & brute-force patch



Hi!

I think I came across a bug in the Levenberg-Marquardt nonlinear least
squares algorithm.  

The problem lies in function lmpar, file multifit/lmpar.c:

In this file, matrix r contains the upper rectangular part of the QR
decomposition of the Jacobian J.  The problem arises from the
following piece of code (ll. 265-271)

 compute_newton_bound (r, newton, dxnorm, perm, diag, w);

  {
    double wnorm = enorm (w);
    double phider = wnorm * wnorm;
    par_lower = fp / (delta * phider);
  }

If matrix r is rank deficient, compute_newton_bound sets vector w to
zero, so that par_lower becomes Inf.  As a consequence, the
LM-parameter par becomes NaN, and the next trial vector contains only
NaNs. 

As a brute-force fix, I added additional checks gsl_finite(par_lower),
gsl_finite(par_upper) around any statements using par_lower or
par_upper (probably unneccessary for par_upper, but just for safety).
The resulting code seems to work, but I think someone with better
understanding of the LM-code should provide a more sensible patch.

My patch is appended as a diff file.  

Please reply directly, I am not subscribed to the mailing list.

Best regards,
Hans

---------------------------------------------------------------------
Dr. Hans Ekkehard Plesser             Tel.  :           +47 6494 8832
Physics Section / ITF                 Fax   :           +47 6494 8810
Agricultural University of Norway     e-mail: hans.plesser@itf.nlh.no
N-1432 Ås, Norway		      WWW   :    arken.nlh.no/~itfhep
---------------------------------------------------------------------

--------------------------------------------------------------------
298c298
<   if (par > par_upper)
---
>   if ( gsl_finite(par_upper) && par > par_upper)
306c306
<   else if (par < par_lower)
---
>   else if ( gsl_finite(par_lower) && par < par_lower)
419c419
<       if (par > par_lower)
---
>       if ( gsl_finite(par_lower) && par > par_lower)
430c430
<       if (par < par_upper)
---
>       if ( gsl_finite(par_upper) && par < par_upper)
445c445,446
<   par = GSL_MAX_DBL (par_lower, par + par_c);
---
>   if ( gsl_finite(par_lower) )    
>     par = GSL_MAX_DBL (par_lower, par + par_c);
--------------------------------------------------------------------


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]