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.

Walking Randomly