Re: [Help-gsl] GSL 2.1 nonlinear least squares question

[Patrick Alken]

I routinely solve problems with millions of residuals (rows) and thousands of 
parameters (columns). In fact for the next release I'm developing a new LM 
solver designed for large systems which can take advantage of sparse structure 
in the Jacobian.

For now you're better off using the 'lmniel' algorithm which will perform 
faster than lmsder for large jacobians


On Feb 18, 2016 4:34 AM, viktor drobot <address@hidden> wrote:
>
> Thank you for explanation!
> But I still doubt is it normal to pass vectors with hundreds of elements 
> (and, so, Jacobians with hudnreds of rows)? My usual computational task 
> requires about 1000 experimental points but sometimes the number of points 
> grows drastically. Does the GSL manages some special way of storing such huge 
> objects?
>
> 2016-02-18 1:22 GMT+03:00 Patrick Alken <address@hidden>:
>>
>> Well in principle you could pass SS as your only residual, but you won't 
>> really benefit from the Levenberg-Marquardt algorithm. The benefit really 
>> comes when you can individually compute the f_i, in your example,
>>
>> f_i = exp(t_i) - model(t_i,c)
>> with 1 <= i <= n
>>
>> where c is a p-vector of your model parameters. The Jacobian will then be a 
>> n-by-p matrix where n is the number of residuals (ie: number of timestamps).
>>
>> What you are trying to do is define a single residual:
>>
>> g_i = sqrt [ sum_i f_i^2 ]
>>
>> with the f_i defined above. In this case your Jacobian would be a 1-by-p 
>> matrix, which I don't think is even allowed in the LM algorithm (I can't 
>> remember exactly where right now but I believe LM requires that the number 
>> of residuals is >= the number of model parameters).
>>
>> So if you insist on a single residual you're probably better off with 
>> another minimization algorithm.
>>
>> The nonlinear least squares solver can numerically approximate the Jacobian 
>> with finite differences, but only if you supply it the full vector of f_i 
>> values.
>>
>> Hope this helps,
>> Patrick
>>
>>
>> On 02/17/2016 03:08 PM, viktor drobot wrote:
>>>
>>> Hello!
>>>
>>>
>>> I'm working on application which determines some kinetic constants from
>>> experimental data graph. Currently I use Nelder-Mead multidimensional
>>> minimization algorithm to find the minimum of least squares functional. My
>>> program computes the sum of squares directly from experimental data:
>>>
>>> SS = sum (exp(time_i) - calc(time_i))^2
>>>
>>> the calc() function solves the initial value problem for some predefined
>>> ODE system which corresponds to the model under investigation and returns
>>> the value of substance concentration at time point time_i. With this
>>> approach I don't need Jacobians and this pretty good because for my case I
>>> can't find the analytical form of Jacobi matrix.
>>>
>>> Recently I've found that GSL 2.1 offers the Levenberg-Marquardt algorithm
>>> with Tikhonov's regularization to solve the nonlinear least squares
>>> problem. Also I've found that these routines can compute the Jacobian
>>> internally with finite difference formulae without external manipulations.
>>>
>>> But I've stucked with gsl_multifit_function_fdf routine. As I understood
>>> the internal implementation of such routines require vector with
>>> experimental data points (in GSL they are named as f_i). This drives me a
>>> bit crazy because I not sure if I can pass the only one function f_i - the
>>> previous SS. I don't want to pass over the whole vector of curve points
>>> because sometimes I have enormous number of input files and so the final
>>> vector size will be too large.
>>>
>>> So I'm asking your, dear experts, what should I do to correctly implement
>>> nonlinear least squares problem solver as described above. Can I pass only
>>> the final sum of squares as the only f_i function or no?
>>>
>>
>>
>
>
>
> -- 
> С уважением,
> Дробот Виктор