[Help-gsl] Re: iteratively re-weighted least squares

[Brian Hawkins]

Hi David,

I don't think there's anything wrong with the approach of re-computing the
weights.  The solver has no knowledge of the weights, only the functions you
provide.  Keep in mind that you should be supplying the residual function,
e.g.

r[i] = (y[i] - f(p, x[i])) * weight[i]

and similarly the rows of the Jacobian are scaled by the weights.  While I
think varying the weights between iterations is probably okay, you should
verify that they are well-behaved.  As a first cut, none are going to 0 or
inf relative to the other weights and/or machine precision.  I don't know to
what extent the GSL solvers are robust to poorly-scaled systems. Also if
your weights vary strongly between iterations, that would suggest you're
taking too-large a step or maybe the weight model is wrong.

The part of the GSL lmsder algorithm that is somewhat mysterious to me is
the "trust region" aspect.  In solvers I've written in the past, I
essentially hand-tuned the step size when necessary.  You might check to see
that you're taking reasonably sized steps (s->dx).  I don't know if there's
any interface for manual control of the step size (didn't see it in the
book).

Of course, try to make sure the parameters are observable in your data set.
 By that I mean you can see a pattern in the data that suggests you're
fitting an appropriate model.  I'm not sure offhand what your model or data
look like.  Maybe you're only testing with a subset of your data, and it's
not enough.

It's also possible that your model is just really tough to fit.  For
example, your performance could be very sensitive to the initial guess.
 What happens when you make your initial guess is very close to a known
solution (on known/simulated data)?

Hope that helps,
Brian

Message: 1
> Date: Sun, 07 Nov 2010 22:02:16 +0100
> From: David Komanek <address@hidden>
> Subject: Re: [Help-gsl] Re: iteratively re-weighted least squares
>        fitting
> To: address@hidden
> Message-ID: <address@hidden>
> Content-Type: text/plain; charset=ISO-8859-2
>
> Dear Brian,
>
> here is the main loop:
>
>    computeWeightsByYi(); // computes weights for the first round
>    do {
>        iter++;
>        status = gsl_multifit_fdfsolver_iterate(s);
>        if (status) break;
>        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
>        if (movingWeights) { // boolean to switch re-weighted and
> "normal" regression
>            gsl_vector *params = gsl_vector_alloc(s->x->size);
>            gsl_vector_memcpy(params, s->x);
>            setParameters(params); // copies current model parameters
> into the enclosing object member data field
>            computeWeightsByFx(); // computes new weights from the new
> set of parameters for the next iteration
>        }
>    }
>    while ((status == GSL_CONTINUE) && (iter < 10000));
>
> The solver is gsl_multifit_fdfsolver_lmsder.
>
> The function is f(x) = a * exp(b * (x+0.5)) / (1 + a * exp(b * (x+0.5)))
>
> With nearly certainty, the problem is on my side, but I wanted to be
> sure that re-weighting is o.k. with GSL or that I will need another
> approach. I am not a matematician, but I am pretty sure the partial
> derivatives etc. are o.k. in my case - my colleague uses another
> computing method with a success, but he has it as a set of VBA macros
> in Excel and we need to do it in C/C++. For some reason I think GSL will
> be better choice than blindly rewrite code lines from VBA to C++. What I
> like to know is if I am using the GSL right for this type of problem.
>
> Thank you. Kind regards,
>
>  David
>
> .
> Dne 6.11.2010 18:46, Brian Hawkins napsal(a):
> > David,
> >
> > What's your convergence criterion?  Is your system full-rank?  Have you
> had
> > success with this problem using a different solver?  I'm having a hard
> time
> > understanding why the GSL solver in particular would be giving you
> trouble.
> >
> > Regards,
> > Brian
> >
> > On Sat, Nov 6, 2010 at 9:01 AM, <address@hidden> wrote:
> >> Message: 1
> >> Date: Fri, 05 Nov 2010 23:09:54 +0100
> >> From: David Komanek <address@hidden>
> >> Subject: [Help-gsl] iteratively re-weighted least squares fitting
> >> To: address@hidden
> >> Message-ID: <address@hidden>
> >> Content-Type: text/plain; charset=ISO-8859-2
> >>
> >> Dear all,
> >>
> >> I ran into problems using weighted least squares fitting in GSL. For
> >> some reason I need to use IRLS modification of this method, so the
> >> weights are recomputed in every iteration. In the case the weights are
> >> computed at the beginning and being constant throug all the iterations,
> >> the procedure works fine. But when I adjust the weights in every
> >> iteration, this usually leads to an error:
> >>
> >> 27 iteration is not making progress towards solution
> >>
> >> I think it is because there are some internally tested conditions and
> >> some of them are not satisfied in this case. For example, in SAS, there
> >> is a special parameter to relax those conditions:
> >>
> >>
> >>
> http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_nlin_sect034.htm
> >>
> >> Is something like this possible with GSL ?
> >>
> >> Thank you in advance.
> >>
> >> David
> >>
>