6
$\begingroup$

Is there a trick to speed up the following calculations?

modelLN2 = Compile[{{AN, _Real}, {tN, _Real}, {tr, _Real}, {ALN, _Real}, {xmf,_Real}, {w, _Real}, {y0, _Real}},Table[AN*Exp[-(i*0.128 - 0.128)/tN] + Exp[-(i*0.128 - 0.128)/tr]*NIntegrate[ALN*Exp[-(1/w^2)*(Log[k/xmf])^2]*Exp[-k*(i*0.128 - 0.128)], {k,0.0001,4.0}, Method -> {"TrapezoidalRule", "SymbolicProcessing" -> False}, AccuracyGoal -> 6], {i, 1.0, 5000.0, 1.0}], RuntimeOptions -> "Speed"]

I was hoping that 'Compile" would be helpful, but it does not change a lot. On my computer it takes approximately 2.5 s to do the above calculations, and I have to call modelLN2 function hundreds of times. Any ideas?

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ Please provide a full usage example for this function. Otherwise we cannot reproduce what exactly too 2.5 seconds. $\endgroup$ Commented Jul 27 at 16:05
  • $\begingroup$ Well, modelLN2 returns a list which is convolved with another list using ListConvolve. After that the resulting list is compared with experimental data to make a fit. In other words, modelLN2 is a theoretical model which is used to fit experimental data via the FindMinimum function - this function calls modelLN2 many times, and it takes about 30 minutes to do the fit to data. The 2.5 s is the time my computer needs to call modelLN2 only once. $\endgroup$ Commented Jul 27 at 17:21
  • 1
    $\begingroup$ This may not apply here, but in some cases e.g. the upper bound of the integral is the variable parameter, the integral can be converted to a differential equation, which can be solved efficiently with NDSolve $\endgroup$ Commented Jul 27 at 17:28
  • 1
    $\begingroup$ With "usage example" I meant an example that shows what the parameters of the function typically look like. I mean, real numeric parameters. This is particularly important when exponential functions are involved. $\endgroup$ Commented Jul 27 at 18:36

4 Answers 4

Reset to default
11
$\begingroup$

A lot of the time in NIntegrate is spent on adaptive subdivision and checking errors. If you need to compute many similar integrals, then a fixed quadrature might make sense. (But there is always the danger that the quadrature rule is badly adapted to some of the functions). Here is a version of your function that takes a list of quadature points points and a list of weights.

integrate = Compile[{{AN, _Real}, {tN, _Real}, {tr, _Real}, {ALN, _Real}, {xmf, _Real}, {w, _Real}, {y0, _Real}, {points, _Real, 1}, {weights, _Real, 1} }, Module[{n, t, point, weight, sum}, n = Min[Length[points], Length[weights]]; Table[ sum = 0.; Do[ point = Compile`GetElement[points, k]; weight = Compile`GetElement[weights, k]; sum += weight Exp[-(1/w^2) (Log[point/xmf])^2] Exp[-point (i 0.128 - 0.128)] , {k, 1, n}]; AN*Exp[-(i*0.128 - 0.128)/tN] + Exp[-(i 0.128 - 0.128)/tr] ALN sum , {i, 1.0, 5000.0, 1.0}] ], CompilationTarget -> "C", RuntimeOptions -> "Speed" ];

Here a usage example:

points = Subdivide[0.0001, 4.0, 1000]; weights = ConstantArray[0., Length[points]]; weights[[1 ;; -2]] += 0.5 Differences[points]; weights[[2 ;; -1]] += 0.5 Differences[points]; aa = modelLN2[1., 1., 1., 1., 1., 1., 1.]; // AbsoluteTiming // First bb = integrate[1., 1., 1., 1., 1., 1., 1., points, weights]; // AbsoluteTiming // First Max[aa - bb]

1.68816

0.037391

1.35239*10^-7

For smooth functions a Gauss quadrature might be more efficient.

Needs["NumericalDifferentialEquationAnalysis`"]; {gpoints, gweights} = NumericalDifferentialEquationAnalysis`GaussianQuadratureWeights[9, 0, 1]\[Transpose]; partition = Subdivide[0.0001, 4.0, 10]; points = Flatten[Plus[ KroneckerProduct[Most[partition], 1. - gpoints], KroneckerProduct[Rest[partition], gpoints] ]]; weights = Flatten[KroneckerProduct[Differences[partition], gweights]]; cc = integrate[1., 1., 1., 1., 1., 1., 1., points, weights]; // AbsoluteTiming // First Max[aa - cc]

0.006521

8.88444*10^-9

That's about 250 times faster than OP's code.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Very nice! Here's +1 $\endgroup$ Commented Jul 27 at 17:49
  • $\begingroup$ Thank you very much! =) $\endgroup$ Commented Jul 27 at 18:34
  • $\begingroup$ Thank you very much for your help. The first part of your code speeds up my calculations five times. The calculation takes 2.6416 s for my code and 0.528483 s for your code. But I do not have a C compiler, unfortunately. Anyway, this really helps! As for the Gauss quadrature, perhaps it is even faster, but when I use it the result is incorrect. In this case, in 0.0642595 seconds I obtained Max[aa-cc]=1.85338. I simply copy-pasted your code. Do you know why my result is different? Do I need to install some Package? $\endgroup$ Commented Jul 27 at 18:41
  • $\begingroup$ Oops. I had accidentally delete the redefinition of the weights in a former edit. Please try again. $\endgroup$ Commented Jul 27 at 19:30
  • 2
    $\begingroup$ Wow, now it works, and it is very fast! I got the result in 0.0497891 s without a C compiler. Impressive. This is precisely what I needed. Thank you very much for the code and for your advice; you helped me a lot! And I will definitively look for a C compiler. $\endgroup$ Commented Jul 27 at 19:47
5
$\begingroup$

The point here is that one can approach compiled speeds with thoughtful Mma code.

For reference:

  • I used @Henrik's example with all parameters set to 1..
  • A high-precision version opHP of the OP's integrals are computed.
  • The OP's modelLN2[1., 1., 1., 1., 1., 1., 1.] took 1.79s on my machine.
  • @Henrik's cc (Gaussian quadrature) example took 0.00359s on my machine.
  • All timings below are under 0.002s.

High-precision version of OP's code (it took almost 27s but timing is irrelevant):

opHP = Block[{AN = 1, tN = 1, tr = 1, ALN = 1, xmf = 1, w = 1, y0 = 1}, Table[ AN*Exp[-(i*0.128`32 - 0.128`32)/tN] + Exp[-(i*0.128`32 - 0.128`32)/tr]* NIntegrate[ ALN*Exp[-(1/w^2)*(Log[k/xmf])^2]*Exp[-k*(i*0.128 - 0.128)] // Rationalize , {k, 0.0001, 4.0} , Method -> {"GaussKronrodRule", "Points" -> 9, "SymbolicProcessing" -> 0}, MinRecursion -> 2, WorkingPrecision -> 32] , {i, 1.0, 5000.0, 1.0}]]; // AbsoluteTiming

A fast, uncompiled version:

{pts, wts} = Most@NIntegrate`GaussRuleData[31, MachinePrecision]; dd = Block[{ AN = 1., tN = 1., tr = 1., ALN = 1., xmf = 1., w = 1., y0 = 1. , k = Rescale[pts, {0, 1}, {0.0001, 4.}] , wts = (4. - 0.0001) wts}, Block[{iRange = Subdivide[0., -4999.*0.128, 4999]}, AN*Exp[iRange/tN] + ALN*(Exp[iRange/tr]*Exp[ Table[(-1/w^2)*(Log[k/xmf])^2 + k*i2 , {i2, iRange/tr}] ] . wts) ] ]; // RepeatedTiming (* {0.00190383, Null} *) MinMax[(opHP - dd)/opHP] (* relative precision *) (* {-8.649*10^-9, 9.59972*10^-9} *)

Compare with compiled code with targets "C" and "WVM" (the default):

int2 = Compile[{ {AN, _Real}, {tN, _Real}, {tr, _Real}, {ALN, _Real}, {xmf, _Real}, {w, _Real}, {y0, _Real} , {k, _Real, 1}, {wts, _Real, 1}}, Block[{ii = Table[(1 - i)*0.128, {i, 1.0, 5000.0, 1.0}]}, AN*Exp[ii/tN] + ALN*(Exp[ii/tr]*Exp[ Table[(-1/w^2)*(Log[k/xmf])^2 + k*i2, {i2, ii/tr}] ] . wts) ] (*,CompilationTarget->"C"*) (*optional*) , RuntimeOptions -> "Speed" ]; (* "C" *) {pts, wts} = Most@NIntegrate`GaussRuleData[31, MachinePrecision]; ee = int2[1., 1., 1., 1., 1., 1., 1., Rescale[pts, {0, 1}, {0.0001, 4.}], (4 - 0.0001) wts]; // RepeatedTiming (* {0.001502, Null} *) (* "WVM" *) {pts, wts} = Most@NIntegrate`GaussRuleData[31, MachinePrecision]; ee = int2[1., 1., 1., 1., 1., 1., 1., Rescale[pts, {0, 1}, {0.0001, 4.}], (4 - 0.0001) wts]; // RepeatedTiming (* {0.00189007, Null} *) MinMax[(opHP - ee)/opHP] (* relative precision *) (* {-8.649*10^-9, 9.59972*10^-9} *)

The compiled-to-C version is significantly faster, between 20% and 25%. The WVM version is just a little faster. (Note: There is variation in both AbsoluteTiming[] and RepeatedTiming[], with less in RepeatedTiming[]. The WVM was usually faster than the uncompiled version. I think it is correct to say it is faster.)

Improve this answer
$\endgroup$
0
5
$\begingroup$

My favorite trick for this kind of thing is to use ParametricNDSolve instead of NIntegrate. It's the easiest way to pre-process the numeric integration:

modelLN3 = With[{ sol = ParametricNDSolveValue[ { y'[k] == ALN*Exp[-(1/w^2)*(Log[k/xmf])^2]*Exp[-k*(i*0.128 - 0.128)], y[0.0001] == 0. }, y[4.0], {k, 0.0001, 4.0}, {ALN, w, xmf, i} ] }, Function[{AN, tN, tr, ALN, xmf, w, y0}, Table[ AN*Exp[-(i*0.128 - 0.128)/tN] + Exp[-(i*0.128 - 0.128)/tr]*sol[ALN, w, xmf, i], {i, 1.0, 5000.0, 1.0}] ] ]; (res1 = modelLN2[1., 1., 1., 1., 1., 1., 1.]); // AbsoluteTiming (res2 = modelLN3[1., 1., 1., 1., 1., 1., 1.]); // AbsoluteTiming res1 - res2 // MinMax

As you can see, it achieves a 5x speed-up without even having to compromise on AccuracyGoal:

Timing comparison

Edit

To explain this trick in a little more detail: if you have an integral of the form:

int[p1_, p2_, ...] := NIntegrate[f[x, p1, p2, ...], {x, x0, x1}]

you can re-cast this as:

int = ParametricNDSolveValue[ {y'[x] == f[x, p1, p2, ...], y[x0] == 0}, y[x1], {x, x0, x1}, {p1, p2, ...} ]

If x1 is not a fixed value, you can add it to the list of parameters as well.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Shouldn't your "favorit trick" use y'[i] inside NDSolve? $\endgroup$ Commented Jul 28 at 14:14
  • $\begingroup$ @UlrichNeumann Explain? The integral is over k, right? The i is for the Table outside the integral. $\endgroup$ Commented Jul 28 at 14:40
  • $\begingroup$ I only know this trick if the upper limit of the integral is the argument of the integral. Sorry for not understanding your approach in detail. $\endgroup$ Commented Jul 28 at 15:09
  • 1
    $\begingroup$ I got it, thanks. How did you get the initial condition y[.001]==.0001? $\endgroup$ Commented Jul 28 at 15:44
  • $\begingroup$ With y[0.0001] == 0. your model modelLN3 gives much better result res1 - res2 // MinMax (*{-9.92388*10^-8, 2.89796*10^-8}*) $\endgroup$ Commented Jul 28 at 20:42
4
$\begingroup$

Just an addition to @HenrikSchumacher`s nice answer.

Without compilation (incl. fast C-compiler) Mathematica(v12.2) & Gaussian quadrature evaluates ~25times faster than QP's result aa . Not so bad I think.

aa = modelLN2[1., 1., 1., 1., 1., 1., 1.]; // AbsoluteTiming (* 3.81s *) Needs["NumericalDifferentialEquationAnalysis`"]; xi = Subdivide[0. , 4.0, 10 ]; xwGauss =Flatten[Map[GaussianQuadratureWeights[9(*5*), #[[1]], #[[2]]] & ,Partition[xi, 2, 1]], 1] ; int = Module[{AN = 1., tN = 1., tr = 1., ALN = 1. , xmf = 1., w = 1., y0 = 1. }, Table[AN*Exp[-(i*0.128 - 0.128)/tN] + Exp[-(i*0.128 - 0.128)/tr]* Total@Map[(ALN*Exp[-(1/w^2)*(Log[#[[1]]/xmf])^2]* Exp[-#[[1]]*(i*0.128 - 0.128)]) #[[2]] &, xwGauss], {i, 1.0, 5000.0, 1.0}] ] (* 0.14s *)
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.