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?
Lx+Ly+Lz == 1
there is no need to calculate all three integrals. CalculateLx
andLy
then setLz=1-Lx-Ly
$\endgroup$NIntegrate
. This should speed things up considerably. This may also reduce numerical noise (originating in subdivision strategies ofNIntegrate
) and helpNonlinearModelFit
to estimate derivates. The method choice here is probably the Gauss-Newton method (or variants of it) which requires first derivatives with respect to the unknown parameters. However, I am not familiar with the internals ofNonlinearModelFit
... $\endgroup$NIntegrate
. $\endgroup$