NonlinearModelFit runs too slow when fitting data to a function that does multiple numerical integrations

Question

I have a function of several parameters defined in a module which does multiple numerical integrations. Then I try to fit that function to some sample data via with NonlinearModelFit, varying several parameters, using differential evolution as the minimization method, and applying constraints. The problem I have run into is that it takes approximately an hour for NonlinearModelFit to find a solution. Since playing with initial guesses is necessary to get a decent fit, this makes it too slow to be useful.

Below, fTest is a a simplified version of the model function I describe above. This function depends only on one parameter, D. It contains the same numerical integrals as my original function but I have removed all the other operations for simplicity. When I try to use it for fitting with NonlinearModelFit, it now takes approximately 7-10 minutes to find a solution.

particleSize = {a2 -> 3, b2 -> 2, c2 -> 1};

fTest[f_, D_?NumericQ] := 
   Module[{a, b, c, prefactor, Lx, Ly, Lz},
     a = a2 /. particleSize;
     b = b2 /. particleSize;
     c = c2 /. particleSize;
     prefactor = (Sqrt[a^2 + D] Sqrt[b^2 + D] Sqrt[c^2 + D])/2;
     Lx = 
       Re[
         prefactor 
           NIntegrate[
             1/((s + a^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), 
             {s, 0, ∞}, 
             WorkingPrecision -> 10]];
     Ly = 
       Re[
         prefactor 
           NIntegrate[
             1/((s + b^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), 
             {s, 0, ∞}, 
             WorkingPrecision -> 10]];
     Lz = 
       Re[
         prefactor 
           NIntegrate[
             1/((s + c^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), 
             {s, 0, ∞}, 
             WorkingPrecision -> 10]];
     (*We know that Lx+Ly+Lz = 1, the below is for testing*)
     f (Lx + Ly + Lz)]

Testing

fTest[10, 0.1]

Generate some fake data to fit to the above function:

dataTest = Table[{n, n + RandomReal[{-0.6, 0.6}]}, {n, 1, 10}];
errorDataTest = 1 & /@ Range[10];

Then I do some non-linear fitting:

AbsoluteTiming[
  nlm = 
    NonlinearModelFit[
      dataTest,{fTest[f, Delta], 0 < Delta < 1}, 
      {{Delta, 0.085}}, f, 
      Weights -> errorDataTest , 
      Method -> {NMinimize, Method -> {"DifferentialEvolution"}}]]

The non-linear fitting takes more than 7-10 min to run, depending of the noise of the data.

Reducing WorkingPrecision to 5 instead of 10 introduces small errors and does not speed the calculation up.

To look at the fit results:

 Show[
   ListPlot[dataTest], 
   Plot[nlm[f], {f, 1, 10}, PlotStyle -> Orange], 
   PlotRange -> All]

 nlm["ParameterTable"]

 Grid[Transpose[{#, nlm[#]} &[{"AdjustedRSquared", "AIC", "BIC", "RSquared"}]], 
   Alignment -> Left]

I have also tested doing a fit varying two parameters instead of one. Defining a second test function, fTest2, with two parameters:

fTest2[f_, D_?NumericQ, offset_] := 
   Module[{a, b, c, prefactor, Lx, Ly, Lz},
     a = a2 /. particleSize;
     b = b2 /. particleSize;
     c = c2 /. particleSize;
     prefactor = (Sqrt[a^2 + D] Sqrt[b^2 + D] Sqrt[c^2 + D])/2;
     Lx = Re[prefactor NIntegrate[1/((s + a^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), {s, 0, ∞}, WorkingPrecision -> 5]];
     Ly = Re[prefactor NIntegrate[1/((s + b^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), {s, 0, ∞}, WorkingPrecision -> 5]];
     Lz = Re[prefactor NIntegrate[1/((s + c^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), {s, 0, ∞}, WorkingPrecision -> 5]];
    (*We know that Lx+Ly+Lz = 1, the below is for testing*)
    offset + f (Lx + Ly + Lz)]

And now trying a non-linear fit varying two parameters instead of one:

AbsoluteTiming[
  nlm = 
    NonlinearModelFit[
      dataTest, {fTest2[f, Delta, offset1], 0 < Delta < 1, 0 < offset1 < 1},
      {{Delta, 0.085}, {offset1, 0}}, f,
      Weights -> errorDataTest ,
      Method -> {NMinimize, Method -> {"DifferentialEvolution"}}]]

The above takes approx. 25 minutes.

For the actual calculations I need to do (as opposed to the above test functions) the original function depends on various parameters and when I try the non-linear fitting varying four parameters, it takes around an hour. Since the fitting is very sensitive to initial guesses, I would need to play a bit with the guesses and also ideally with the parameters of the DifferentialEvolution method ("ScalingFactor", "CrossProbability"), but I cannot do this if it takes approximately one hour per attempt.

Is there any way to improve the speed or is this it?

---

After the comments I now have a faster 2-parameter test function (fTest4 below) that does Lz=1-Lx-Ly, and also I have turned SymbolicProcessing off in NIntegrate. The two changes make the two-parameter fitting (see below) take approx. 6 minutes instead of approx. 25 minutes for fTest2 above. So that's already a substantial improvement:

fTest4[f_, D_?NumericQ, offset_] := Module[{a, b, c, prefactor, Lx, Ly, Lz},
 a = a2 /. particleSize;
 b = b2 /. particleSize;
 c = c2 /. particleSize;
 prefactor = (Sqrt[a^2 + D] Sqrt[b^2 + D] Sqrt[c^2 + D])/2;
 Lx = Re[prefactor NIntegrate[1/((s + a^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), {s, 0, \[Infinity]}, WorkingPrecision -> 10, 
  Method -> {Automatic, "SymbolicProcessing" -> 0}]];
 Ly = Re[prefactor NIntegrate[1/((s + b^2 + D) Sqrt[(s + a^2 + D) (s + b^2 + D) (s + c^2 + D)]), {s, 0, \[Infinity]}, WorkingPrecision -> 10, 
  Method -> {Automatic, "SymbolicProcessing" -> 0}]];
 Lz = 1 - Lx - Ly;
 (*We know that Lx+Ly+Lz = 1, the below is for testing*)
 offset + f (Lx + Ly + Lz)
]

Two-parameter fitting:

AbsoluteTiming[nlm = NonlinearModelFit[dataTest,
 {fTest4[f, Delta, offset1], 0 < Delta < 1, 0 < offset1 < 1},
 {{Delta, 0.085}, {offset1, 0}},
 f,
 Weights -> errorDataTest ,
 Method -> {NMinimize, Method -> {"DifferentialEvolution"}}
]]

Any more suggestions for speeding this up?

Answer 1

This is sort of an expansion of my comment above.

Function to integrate:

F = {s, a, b, c, d} \[Function] (Sqrt[a^2 + d] Sqrt[b^2 + d] Sqrt[c^2 + d])/2/((s + a^2 + d) Sqrt[(s + a^2 + d) (s + b^2 + d) (s + c^2 + d)]);

Transformation of the unit interval to $[0,\infty]$:

\[Phi] = t \[Function] Evaluate[Simplify[1/(1 - t) - 1]];

Computation of evaluation points and integration weights by using Gauss quadrature and transformation formula for integrals (incorporating @J.M.'s comment):

m = 3;
n = 100;
{p, \[Omega]} = Developer`ToPackedArray[N[Most[NIntegrate`GaussRuleData[m, 16]]]]/n;
(* quadrature data for the unit interval *)
pts0 = Join @@ Table[p + i/n, {i, 0, n - 1}];
weights0 = Join @@ ConstantArray[\[Omega], n];
(* quadrature data for the interval from 0 to Infinity *)
pts = (\[Phi] /@ pts0);
weights = (\[Phi]' /@ pts0) weights0;

Compilation of function:

Block[{s, a, b, c, d},
  With[{code = F[s, a, b, c, d]},
   cF = Compile[{{s, _Real}, {a, _Real}, {b, _Real}, {c, _Real}, {d, _Real}}, code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True
     ]
   ]
];

Test example:

a = 3.; b = 2.; c = 1.; d = 1.;
bla = Simplify[Integrate[F[s, a, b, c, d], s]];
int = F[s, a, b, c, d];
trueLx = Limit[bla, s -> \[Infinity]] - Limit[bla, s -> 0]; //AbsoluteTiming // First
Lx = Re[ NIntegrate[int, {s, 0, \[Infinity]}, WorkingPrecision -> 5]];// AbsoluteTiming // First
myLx = weights.cF[pts, a, b, c, d]; // AbsoluteTiming // First
trueLx - myLx
trueLx - Lx

(* 0.001997 *)
(* 0.003963 *)
(* 0.000082 *)
(* -0.0000127242 *)
(* -0.0000624892 *)

Thus, with comparable error, we have a 50x speed-up. Accuracy can be increased by increasing m and n. See also the documentation of NumericalDifferentialEquationAnalysis`GaussianQuadratureWeights (this function is another way to compute the Gauss quadrature points and weights).

However, this method relies on the assumption that the integrant is several times differentiable (2 m-fold) and has no singularities in $[0,\infty]$. Moreover, it is hard to obtain more that single precision accuracy with this method.