Re: [Help-gsl] strange NLS behaviour

[Patrick Alken]

Hello,

  Glad you're using the new interface - it is much more powerful than 
the old one. I will try to address your questions, though you didn't 
give any source code so it is difficult to diagnose exactly what's going on.

1. Since you're using a finite-difference approximation to the Jacobian, 
the output should say that there are no Jacobian evaluations. The 
counter is incremented only when your _df function is called. One thing 
you might try is to use centered differences to compute J instead of 
forward. This will result in a more accurate J (the old 
multifit_fdfsolver interface only allowed forward differences). It will 
cost you about twice as many function evaluations, but it may converge 
in fewer iterations due to a more accurate Jacobian matrix. Based on 
your "time elapsed" it looks like you can compute your residuals pretty 
quickly, so doubling the number of function calls may not matter much. 
You can set this via:

fdf_params.fdtype = GSL_MULTIFIT_NLINEAR_CTRDIFF

2. chisq/dof - I personally don't use this statistic since I find it 
difficult to interpret. chisq is the sum of squares of your residuals, 
which often have units. Therefore the chisq/dof depends on the units of 
the problem, so you can always scale your units in a way which will make 
chisq/dof as small or as big as you want. Perhaps it is more useful if 
you happen to have dimensionless residuals.

I think its better to plot the residuals themselves to determine the 
quality of the fit. If you see any systematic structure in the residuals 
then there could be problems with the fit (or with the underlying 
model). But you're obtaining parameter values very close to the true 
values so it seems ok to me.

3. Condition number of J - just curious why you think its too big? Do 
you have reason to suspect it should be close to 1? Since you're using a 
finite-difference Jacobian, that will probably make the conditioning 
worse - I think the conditioning of a finite difference operator goes 
something like O(1/h^2), so the conditioning gets worse as you reduce 
the step size. But in any case a condition number of 1e6 is reasonable 
and shouldn't affect the accuracy of the method too much, especially if 
you're using the QR (default) solver.

4. Error estimates - I don't see where you list the error estimates on 
your parameters, just a correlation matrix. Are you using the same error 
estimation as in the "Exponential Fit" example program? That program 
uses the chisq/dof to estimate the parameter error. In your case, chisq 
is tiny so its not surprising you end up with tiny errors. If your 
residuals are all really small, then it might be true that the parameter 
errors are also small.

So, it seems to me your program is working correctly. Are you getting 
different results when using the old interface (fdfsolver)?

Another thing I'd recommend is to run your program with 
Levenberg-Marquardt and then one of the dogleg methods. You might find 
that one method signficantly outperforms the other. For example, in a 
inverse problem I'm working on, I found that LM takes over 20 iterations 
to converge, while dogleg only takes 8, so its worth trying different 
methods.

Patrick

On 09/13/2016 05:59 AM, viktor drobot wrote:
> Hi, community!
>
> I model kinetics of enzymatic reactions. There are problems for which it's
> necessary to find suitable values of kinetic parameters so the real kinetic
> data could be described by our model.
>
> To facilitate my test runs I generate some "real kinetic data" by solving
> the system of ODEs with known kinetic parameters and then add some Gaussian
> noise. This is the reference curves. After that I try to solve the inverse
> problem - find out the most appropriate values of these parameters with
> initial approximation that is close enough to the true values. I use GLS
> NLS routines (v 2.2.1) for that.
>
> First of all I solve ODEs with another approximation and then find the
> differences (Y[t] - y[t]) between reference points Y[t] (generated as
> stated above) and the new calculated y[t], where 't' is the time point.
> Usually I process several kinetic curves at once. I store these differences
> in 'f' vector of function to be minimized. Because of the nature of the
> problem there is no analytical Jacobian possible so I set corresponding
> field of gsl_multifit_nlinear_fdf structure to NULL. In case of unweighted
> NLS I multiply the diagonal elements of the covariance matrix by the
> variance of the residuals. Then the NLS solver does all work for me and I
> get some strange results.
>
> The output says that there was no Jacobian evaluations. Kinetic parameters
> obtained seems to be quiet right but other statistics weird enough. Look at
> the example below.
>
> InverseSolver v0.1.0
> ********************
> Run ended at Thu Sep  8 17:13:27 2016
>
>
> ====================
> Optimized parameters
> --------------------
> Kcat = 2.4849549376350680e+01 +/- 2.3301006770134422e-01 (true value is
> 2.5e+01)
> Km   = 2.5962066490476204e-05 +/- 1.2805709675629306e-06 (true value is
> 2.5e-05)
> Kp   = 5.6728528418318680e-05 +/- 3.3698382714302584e-06 (true value is
> 5.0e-05)
>
>
> ====================
> Minimization details
> --------------------
> The fit converged after 13 iteration(s)
> Reason for stopping: small step size
> Objective function evaluations = 69
> Number of data points          = 380
> Model parameters               = 3
> dof                            = 377
> chi                            = 2.7339109755226917e-05
> chisq                          = 7.4742692220834363e-10
> chisq/dof                      = 1.9825647803934845e-12
> Time elapsed: 0.05957 secs
>
>
> ====================================
> Correlation matrix of fit parameters
> ------------------------------------
>       Kcat        Km          Kp
> Kcat  1.0000e+00  8.5928e-01  6.5048e-01
> Km    8.5928e-01  1.0000e+00  9.4124e-01
> Kp    6.5048e-01  9.4124e-01  1.0000e+00
>
>
> ==================
> Iterations history
> ------------------
> Iter # Kcat                    Km
> Kp                      cond(J)                 chi
> 0       2.0000000000000000e+01  2.0000000000000002e-05
> 2.0000000000000002e-05                     inf  3.1567955122657597e-04
> 1       2.7173704932869438e+01  2.6087621471921223e-05
> 3.0462571586587205e-05  2.9025496731552132e+06  8.3486976629633576e-05
> 2       2.4543231868632770e+01  2.5481466984423840e-05
> 4.6317587677535904e-05  3.0659026221618010e+06  4.9189276130628142e-05
> 3       2.4799130142472123e+01  2.5333806359622357e-05
> 5.3688047773966246e-05  2.1896626012014118e+06  2.7592054291596016e-05
> 4       2.4830032619351048e+01  2.5841007936620306e-05
> 5.6370153753761361e-05  2.0525957706617371e+06  2.7340238653871502e-05
> 5       2.4848466472445963e+01  2.5956371903529538e-05
> 5.6716387131234307e-05  2.0040307592909317e+06  2.7339110581465998e-05
> 6       2.4849535469371190e+01  2.5961998929997054e-05
> 5.6728382061723916e-05  1.9990772957287871e+06  2.7339109755370303e-05
> 7       2.4849547296667978e+01  2.5962054330720485e-05
> 5.6728499949112989e-05  1.9989372038075817e+06  2.7339109755230966e-05
> 8       2.4849547539278458e+01  2.5962056427776209e-05
> 5.6728505891823464e-05  1.9989391076721402e+06  2.7339109755230180e-05
> 9       2.4849549814327389e+01  2.5962073280806688e-05
> 5.6728551401861829e-05  1.9989361463323494e+06  2.7339109755228415e-05
> 10      2.4849550175380344e+01  2.5962072172440512e-05
> 5.6728543239809625e-05  1.9989390874822082e+06  2.7339109755227405e-05
> 11      2.4849549199258245e+01  2.5962065484919628e-05
> 5.6728525792610483e-05  1.9989382142781706e+06  2.7339109755227290e-05
> 12      2.4849549376350691e+01  2.5962066490476238e-05
> 5.6728528418318701e-05  1.9989384807737938e+06  2.7339109755226951e-05
> 13      2.4849549376350680e+01  2.5962066490476204e-05
> 5.6728528418318680e-05  1.9989391375842493e+06  2.7339109755226917e-05
>
>
> =============================
> Per-file Fischer's statistics
> -----------------------------
> [/home/viktor/science/enzyme_kinetics/data/test/172-186/0172_1s-0075.dat]
>                    dof = 51
> Number of data points = 54
>           F-statistics = 3.2567463900742809e+04
>
> [/home/viktor/science/enzyme_kinetics/data/test/172-186/0173_1s-0100.dat]
>                    dof = 62
> Number of data points = 65
>           F-statistics = 1.0932769130444922e+04
>
> [/home/viktor/science/enzyme_kinetics/data/test/172-186/0174_1s-0125.dat]
>                    dof = 72
> Number of data points = 75
>           F-statistics = 2.2345670716799857e+04
>
> [/home/viktor/science/enzyme_kinetics/data/test/172-186/0175_1s-0150.dat]
>                    dof = 84
> Number of data points = 87
>           F-statistics = 4.4888120828034567e+04
>
> [/home/viktor/science/enzyme_kinetics/data/test/172-186/0176_1s-0175.dat]
>                    dof = 96
> Number of data points = 99
>           F-statistics = 3.9235545459615321e+04
>
> chisq/dof value is too small. It's much smaller than 1. Condition number of
> Jacobian is too big, much bigger than 1. Error estimates are too small
> despite of some Gaussian noise in my generated data. Seems like the
> calculation is unstable.
>
> I completely stuck with these strange results. What I'm doing wrong? Any
> help is greatly appreciated!