Why does NIntegrate become much faster if I evaluate a function, then compile it, then make it numerical?

Question

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?

Answer 1

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!