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

[David Komanek]

Well, to reply myself with a bit more information:

I just had the oppoortunity to compare weighted sums of squares for my
computation by GSL and the other methods (SAS and 3rd party handmade by
one of our mathematicians) for the case weights are static through all
the iterations. All three methods have the same input data. SAS and 3rd
are producing same results with sum of 291,1819806 while GSL computes
slightly different coefficients by sum of 297,6829277 ? Probably my
question is dumb, but is there a possibility that this difference can be
caused by some limitation in implementation LM solvers ?

I'm using pure C++ with "doubles" as storage for the numbers. The
solution is for function

y = (a * exp (b * x)) / (1 + (a * exp (b * x)));
first derivatives expressions were authorized by a mathematician
count of [x;y] pairs is 26 and all are for 65.5 till 90.5 (step 1.0)

Thank you in advance for any opinion on this.

Sincerely,

  David
 

On 11/18/2010 12:15 AM, David Komanek wrote:
> Thank you, Brian.
>
> I just tried to get to the roots, so I took the example program from the
> manual and it works fine - I get exactly the same results as stated in
> the example. But with my own function, the fdf solver (both lmder and
> lmsder) finds another solution than SAS and 3rd party VBA routine (which
> gives the same results as SAS). It all is true for static weights only.
> Re-weighted ILS still does not work for me well. So I am giving up and
> will rewrite the VBA code to c++ to satisfy my colleagues. Maybe I
> should write my own solver function instead (more suitable to my
> particular type of problem) to still use GSL, but I have no idea how to
> do this, my math skills are poor, which is probably the reason for all
> my trouble.
>
> Thank you for your assistance.
>
> Best regards,
>
>   David
>
>
> Dne 17.11.2010 6:47, Brian Hawkins napsal(a):
>> Hi David,
>>
>> The parameters to the function are just a void pointer, so the solver
>> machinery has no way to know what's in there.  You're free to use it however
>> you like.
>>
>> Good luck,
>> Brian
>>
>> On Mon, Nov 15, 2010 at 9:01 AM, <address@hidden> wrote:
>>> Message: 1
>>> Date: Mon, 15 Nov 2010 02:20:05 +0100
>>> From: David Komanek <address@hidden>
>>> Subject: Re: [Help-gsl] Re: iteratively re-weighted least squares
>>> To: address@hidden
>>> Message-ID: <address@hidden>
>>> Content-Type: text/plain; charset=ISO-8859-2
>>>
>>> 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
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl