5
$\begingroup$

I am using Mathematica for numerical integrals in optics. By accident, I found out that NIntegrate becomes faster by an order of magnitude for a function I'm integrating, if I take the following steps:

  1. Evaluate[] the function;
  2. Compile the function;
  3. Define a new function that takes only numerical arguments and applies them to the compiled function.

I haven't been able to reproduce this with other functions; only this one so far. I do know that this integrand is highly oscillatory. Here's my minimal working example:

ClearAll["Global`*"];
ClearSystemCache[];

gaussians[x_] :=
  0.1 Sum[
    Exp[-34000. ((2 m)*0.00004)^2]*1/((Pi*0.00001^2)^0.25)*
     Exp[-(x - (2 m)*0.00004)^2/(2*0.00001^2)],
    {m, -200, 200}
    ];
integrand[x_, kx_] := gaussians[x] Exp[-I*kx*x];
numericIntegrand[x_?NumericQ, kx_?NumericQ] := integrand[x, kx]
compiledIntegrand = 
  Compile[{x,kx}, integrand[x, kx]];
numericalCompiledIntegrand[x_?NumericQ, kx_?NumericQ] := 
  compiledIntegrand[x, kx];
compiledEvaluatedIntegrand = 
  Compile[{x, kx}, Evaluate[integrand[x, kx]]];
numericalCompiledEvaluatedIntegrand[x_?NumericQ, kx_?NumericQ] := 
  compiledEvaluatedIntegrand[x, kx];
integral[f_, kx_] := 
  NIntegrate[f[x, kx], {x, -0.008, 0.008}, MaxRecursion -> 15(*,
   AccuracyGoal\[Rule]12,PrecisionGoal\[Rule]9*)];

ParallelTable[TimeConstrained[
  AbsoluteTiming[integral[function, 150000.]],
  60, "Exceeded time limit"
  ],
 {function, {integrand, numericIntegrand, compiledIntegrand, 
   numericalCompiledIntegrand, compiledEvaluatedIntegrand, 
   numericalCompiledEvaluatedIntegrand}}
 ]

And the output:

CompiledFunction::cfsa :  Argument x at position 1 should be a machine-size real number.

CompiledFunction::cfsa :  Argument x at position 1 should be a machine-size real number.

NIntegrate::slwcon :  Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::slwcon :  Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::slwcon :  Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::slwcon :  Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::slwcon :  Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::slwcon :  Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

{{55.2211, -3.03556*10^-6 + 
   2.87568*10^-19 I}, "Exceeded time limit", {53.982, -3.03556*10^-6 \
+ 2.87568*10^-19 I}, "Exceeded time limit", {55.0044, -3.03556*10^-6 \
+ 2.87568*10^-19 I}, {4.56948, -3.03556*10^-6 + 2.87568*10^-19 I}}

As you can see from the table output, only this specific sequence of steps gives the speedup I'm referring to. integrand takes the integrand as is, numericIntegrand uses a "numerical" version of the integrand (i.e. one that takes only numerical values), compiledIntegrand uses a compiled version of the integrand, numericalCompiledIntegrand uses a numerical version of the compiled integrand, compiledEvaluatedIntegrand uses a compiled version of the integrand with Evaluate[] applied before compilation, and finally, numericalCompiledEvaluatedIntegrand applies the steps I outlined above. All of them take over 50 seconds, except for numericalCompiledEvaluatedIntegrand, which takes less than 5 seconds. What could be the cause of this?

$\endgroup$
5
  • 2
    $\begingroup$ If you evaluate CompiledFunctionTools`CompilePrint[compiledIntegrand] you will see MainEvaluate[ Hold[integrand][ R0, R1]] which means that the integrand is not correctly compiled withour Evaluate. This is why you cannot get any speed-up from that. $\endgroup$ Mar 18, 2021 at 17:49
  • $\begingroup$ When I evaluate that command, I get CompiledFunctionTools`CompilePrint[CompiledFunction[Argument count: 2 Argument types: {_Real,_Real} ]]. Doesn't that mean it is evaluated in my version (12.1.1.0, Student Edition, 64-bit Windows)? $\endgroup$
    – KBS
    Mar 18, 2021 at 18:03
  • 1
    $\begingroup$ Hm, that's odd. Try to execute Needs["CCodeGenerator"]` before you call CompiledFunctionTools`CompilePrint. $\endgroup$ Mar 18, 2021 at 18:25
  • 2
    $\begingroup$ If you use immediate instead of delayed assignments, a few of your issues will go away already. $\endgroup$
    – Roman
    Mar 18, 2021 at 18:33
  • 1
    $\begingroup$ @HenrikSchumacher That does it (should be Needs["CCodeGenerator`"], but the markup language messed up the command, so you need to put two backticks before and two after the code). Thanks. $\endgroup$
    – KBS
    Mar 18, 2021 at 18:34

1 Answer 1

7
$\begingroup$

The reason for Evaluate is that Compile is surprisingly bad at recognizing the expression that it is suppose to optimize. One can use Evaluate as a workaround, but I prefer to use With as follows (Block is supposed to shield x and kx from external definitions):

cf = Block[{x, kx},
   With[{code = integrand[x, kx]},
    Compile[{{x, _Real}, {kx, _Real}},
     code,
     CompilationTarget -> "C"
     ]
    ]
   ];

The above also uses CompilationTarget -> "C" which will attempt to compile the code if your system has a C compiler installed. That will make the execution even faster.

Now we wrap the compiled code by a "numeric" function:

f[x_?NumericQ, kx_?NumericQ] := compiledEvaluatedIntegrand[x, kx];

The point of this is that NIntegrate attempts to analyze every analytic input to figure out which integration strategy is the most promising and to compile the expression if possible. For some reason, it does not realize that we have compiled our function already and needlessly inserts symbolic values into it which kicks off a cascade of error messages. The pattern _?NumericQ just prevents this symbolic evaluation pass.

Finally we can do the integration. Since the integrand is highly oscillatory, "LocalAdaptive" is quite a reasonable choice for the integration method.

kx = 150000.;
NIntegrate[f[x, kx], {x, -0.008, 0.008}, Method -> "LocalAdaptive"]

On my machine, this is more then 4 times faster than numericalCompiledEvaluatedIntegrand with Method -> "LocalAdaptive" and 13 times faster then numericalCompiledEvaluatedIntegrand without Method -> "LocalAdaptive".

Edit

Just to point out the commentbelow by Michael E2: The following combination makes the use of the pattern _?NumericQ obsolete:

cf2 = Block[{x, kx},
   With[{code = integrand[x, kx]},
    Compile[{{x, _Real}, {kx, _Real}},
     code,
     CompilationTarget -> "C",
     RuntimeOptions -> {"EvaluateSymbolically" -> False}
     ]
    ]
   ];
NIntegrate[cf2[x, kx], {x, -0.008, 0.008}, Method -> {"LocalAdaptive", "SymbolicProcessing" -> False}]

And it is even another 20% faster!

$\endgroup$
6
  • 4
    $\begingroup$ You can also make a compiled function be "numeric" with RuntimeOptions -> {"EvaluateSymbolically" -> False}. $\endgroup$
    – Michael E2
    Mar 18, 2021 at 19:19
  • 2
    $\begingroup$ Hm. I have just tried that but NIntegrate still throws two(!) errors "CompiledFunction::cfsa: Argument -x at position 1 should be a machine-size real number.". =/ $\endgroup$ Mar 18, 2021 at 19:26
  • 2
    $\begingroup$ Odd that cf[-x, 150000.] gives that message but cf[x, 150000.] does not. Nintegrate is checking symmetry. With Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, which makes sense for a compiled function, the message goes away. $\endgroup$
    – Michael E2
    Mar 18, 2021 at 19:31
  • 2
    $\begingroup$ Ah! Great! I had tried "SymbolicProcessing" -> False before but to no avail. It seems to work only in combination with RuntimeOptions -> {"EvaluateSymbolically" -> False}. $\endgroup$ Mar 18, 2021 at 19:36
  • 1
    $\begingroup$ Don't you think the message with cf2[-x, 150000.] is a bug? I can't quite figure out the pattern. I thought it was the -1 factor in Times[-1, x]. But I get the error for cf2[Exp[x], 150000.] but not for cf2[Sin[x], 150000.]. Maybe the trigger is whether the head is Plus, Times, or Power. Apparently, it is not using NumericQ. $\endgroup$
    – Michael E2
    Mar 19, 2021 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.