Nonlinear model fitting in Mathematica: The quest for more speed

March 29th, 2010 | Categories: math software, mathematica, programming | Tags:

A Mathematica user recently asked me if I could help him speed up the following Mathematica code. He has some data and he wants to model it using the following model

dFIO2 = 0.8;
PBH2O = 732;
dCVO2[t_] := 0.017 Exp[-4 Exp[-3 t]];
dPWO2[t_?NumericQ, v_?NumericQ, qL_?NumericQ, q_?NumericQ] := 
 PBH2O*((v/(v + qL))*dFIO2*(1 - E^((-(v + qL))*t)) + 
    q*NIntegrate[E^((v + qL)*(\[Tau] - t))*dCVO2[\[Tau]], {\[Tau], 0, t}])

Before starting to apply this to his actual dataset, he wanted to check out the performance of the NonlinearModelFit[] command. So, he created a completely noise free dataset directly from the model itself

data = Table[{t, dPWO2[t, 33, 25, 55]}, {t, 0, 6, 6./75}];

and then, pretending that he didn’t know the parameters that formed that dataset, he asked NonlinearModelFit[] to find them:

fittedModel = NonlinearModelFit[data, dPWO2[t, v, qL, q], {v, qL, q}, t]; // Timing
params = fittedModel["BestFitParameters"]

On my machine this calculation completes in 11.21 seconds and gives the result we expect:

{11.21, Null} 
{v -> 33., qL-> 25., q -> 55.} 

11.21 seconds is too slow though since he wants to evaluate this on real data many hundreds of times. How can we make it any faster?

I started off by trying to make NIntegrate go faster since this will be evaluated potentially hundreds of times by the optimizer. Turning off Symbolic pre-processing did the trick nicely in this case. To turn off Symbolic pre-processing you just add

Method -> {Automatic, "SymbolicProcessing" -> 0}

to the end of the NIntegrate command. So, we re-define his model function as

dPWO2[t_?NumericQ, v_?NumericQ, qL_?NumericQ, q_?NumericQ] := 
 PBH2O*((v/(v + qL))*dFIO2*(1 - E^((-(v + qL))*t)) + 
    q*NIntegrate[E^((v + qL)*(\[Tau] - t))*dCVO2[\[Tau]], {\[Tau], 0, t}
,Method -> {Automatic, "SymbolicProcessing" -> 0}])

This speeds things up almost by a factor of 2.

fittedModel = NonlinearModelFit[data, dPWO2[t, v, qL, q], {v, qL, q}, t]; // Timing
params = fittedModel["BestFitParameters"]

{6.13, Null}
{v -> 33., qL -> 25., q -> 55.}

Not bad but can we do better?

Giving Mathematica a helping hand.

Like most optimizers, NonlinearModelFit[] will work a lot faster if you give it a decent starting guess. For example if we provide a close-ish guess at the starting parameters then we get a little more speed

fittedModel = 
   NonlinearModelFit[data, 
    dPWO2[t, v, qL, q], {{v, 25}, {qL, 20}, {q, 40}}, 
    t]; // Timing
params = fittedModel["BestFitParameters"]

{5.78, Null}
{v -> 33., qL -> 25., q -> 55.}

give an even better starting guess and we get yet more speed

fittedModel = 
   NonlinearModelFit[data, 
    dPWO2[t, v, qL, q], {{v, 31}, {qL, 24}, {q, 54}}, 
    t]; // Timing
params = fittedModel["BestFitParameters"]

{4.27, Null}
{v -> 33., qL -> 25., q -> 55.}

Ask the audience

So, that’s where we are up to so far. Neither of us have much experience with this particular part of Mathematica so we are hoping that someone out there can speed this calculation up even further. Comments are welcomed.

  1. March 29th, 2010 at 18:12
    Reply | Quote | #1

    I’m quite surprised it does it this fast; Your ‘fit-model’ is really nasty (the model includes integrations of exp of exp of something!)

    I don’t have a lot of experience with the NonlinearModelFit command; however I’ve seen my code speed up by giving all the constraint for the variables (for example that they have to be positive, and below a certain value). And as you already said specifying a nice guess can help quite a bit.

    Have you tried using Compile ?

    Good luck with this tough one!

  2. March 31st, 2010 at 10:18
    Reply | Quote | #2

    Hi Sander

    I agree – a few seconds isn’t bad considering the complexity of the model. Compile doesn’t help apparently (tried by the user, not me).

    Cheers,
    Mike

  3. Graeme
    February 17th, 2014 at 23:50
    Reply | Quote | #3

    Hi,

    I recently started doing something similar (running many bootstrap iterations using NLMF) and I was wondering: do you know of a way to speed up the NMinimize method? I’ve done transforms so that the parameter space is the real line, so the most exhaustive method is preferable (automatic is much faster but yields different results most of the time). I assume turning off symbolic wouldn’t help much here?

    Thanks,
    -Graeme