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

[David Komanek]

Hi Brian,

I use the weights both in r[i] and Jacobians. What I realized after
hours of playing with my code is that the coefficient values are quite
different than the known solution from a different software (on the same
input data). So, there is probably some other problem in my code and
re-weighting only makes it more important and visible. Now I have a
different set of input data for which re-weighting does not trigger the
solver error. In this case, re-weighting makes final values of
coefficients slightly different with a better R-squared value than with
fixed weights. But still it is far away from coefficients computed by
another software.

I really tried find the problem in my code, but with no success yet. The
only thing, which I think is really different than in the example in the
gsl documentation (with static weights, of course), it seems to be the
"data" structure, which is used as a second parameter of func_f and
func_df functions. I hope thiese structure data are used only in the
code of these functions and not by any other internal solver procedure,
so this structure can have all the members I need (not only n, y and
sigma) in any order. Is it right ? If yes, then I really do not know
what could be wrong, but it would mean the problem is somewhere in my
head :-)

Thanks again. Sincerely,

  David


Dne 10.11.2010 6:44, Brian Hawkins napsal(a):
> 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
>>>>
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl