This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Bug in Levenbeg-Marquardt Algorithm & brute-force patch
- From: Hans Ekkehard Plesser <hans dot plesser at itf dot nlh dot no>
- To: gsl-discuss at sources dot redhat dot com
- Date: Mon, 25 Feb 2002 14:38:56 +0100 (CET)
- Subject: Bug in Levenbeg-Marquardt Algorithm & brute-force patch
- Reply-to: hans dot plesser at itf dot nlh dot no
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);
--------------------------------------------------------------------