# Mathematica 10.2 is slow to NIntegrate highly oscillatory integrals

I noticed that Mathematica V10.2 is taking much longer to run my code when compared with V9. One function that seems affected is NIntegrate. For example, for a highly oscillatory integrand like

f = Interpolation[Table[{x, Cos[x 100]}, {x, 0, π, 10^-6}]]


evaluating the numerical integration

NIntegrate[f[x], {x, 0, 30 π/400},
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]//AbsoluteTiming


took 0.00286sec (on V9) and 0.48sec (on V10.2) -- a factor of 160.

Would you know why Mathematica 10.2 is taking so much longer to calculate some functions?

• Hi there:) I edited your question for you and added code formatting. To learn how to properly format your question yourself please have a look here Jan 22 '16 at 19:41
• NIntegrate isn't every supposed to be ran on an Interpolation. So they probably didn't care too much about making sure it runs fast. They've upgraded a lot of things with Interpolation since version 9. Perhaps version 9 had a special shortcut it used directly on the InterpolationObject instead of calling it like a numerical function? Who knows. Jan 22 '16 at 22:26
• Actually, judging from the results, it's probably the opposite. I bet version 10 tries to do something smart behind the scenes with the Interpolation instead of treating it as a function. Jan 22 '16 at 22:33
• That was my first guess, but the rule "SymbolicProcessing" -> 0 should avoid MM of doing smart things when performing numerical calculations. Jan 23 '16 at 13:21

Having put in some time trying to see what's going, I've found a few things, but I don't have a perfectly clear picture. I believe the issue is with the large InterpolatingFunction in the integrand and not with NIntegrate per se. The time it takes for NIntegrate to set up the integration is much longer in V10, but the integration itself runs in about the same amount of time in V10 as in V9. One can get around the issue with ?NumericQ.

What is unclear still is how the size of the InterpolatingFunction comes to be an issue inside NIntegrate. Is it something that NIntegrate does? Or is it something Mathematica does that is outside the control of NIntegrate? Why in V10 but not in V9? I can't be sure. I suspect it is outside NIntegrate's control. It's probably worth reporting to WRI, even if it's not exactly a bug. They might be able to fix the problem in a future version.

One can modify Anton Antonov's tracing code posted in his answer to Determining which rule NIntegrate selects automatically to construct a rudimentary profiling tool. Anton expresses some reservations about his code, so you should read his post before using the code below. It generates a moderate amount of output, but it would make this a very large answer to include it all. What I will do is summarize the evidence and give the code to generate it.

## Issues with the integrand

The OP's interpolating function is an exact symbolic interpolating function. It takes several seconds to create it. Perhaps it is needed in that form, but if machine-precision numerical evaluation is the only purpose, it would be better (from the point of view of speed) to create it as a machine-precision function. Apparently there is only a speed issue in V10. One can get around it in V10 by protecting the interpolating function with ?NumericQ. One recovers timings similar to V9 this way, which suggests that it is not the oscillatory behavior of the integrand that is causing the slow evaluation.

f = Interpolation[Table[{x, Cos[x 100]}, {x, 0, π, 10^-6}]]; // AbsoluteTiming
fN = Interpolation[Table[{x, Cos[x 100]}, {x, 0., π, 10^-6}]]; // AbsoluteTiming
fNQ[x_?NumericQ] := f[x];
fNNQ[x_?NumericQ] := fN[x];
(*
{21.0106, Null}
{9.84533, Null}
*)

ByteCount[f]
ByteCount[fN]
(* V9
904756768  -- f   0.9 GB
75399504  -- fN
V10.3.1
703705096  -- f   0.7 GB
75399608  -- fN
*)


Tangential remark? One difference in how V9 and V10 treats interpolating functions is how they are typeset. This shouldn't affect NIntegrate, but maybe the underlying causes are related. The difference is time is amazing. Create two input cells as below, select both cells, and execute them simultaneously (shift-return). Make sure $HistoryLength is at least 2. (* cell 1 *) SessionTime[] f[x] // AbsoluteTiming (* cell 2 *) SessionTime[] - %%  The second cell won't be executed until the output from the first cell is done typesetting. The timing from AbsoluteTiming is the same in V9 and V10, about 4.*10^-6. But the wall-clock time between executing the cells, roughly the time spent typesetting, is about 20 times longer in V10. (* Differences in session time *) 0.147116 -- V9 3.362884 -- V10.3.1  ### Timings of the four integrands Running the OP's NIntegrate command on the four functions {f, fN, fNQ, fNNQ} yields the following timings. In V10, the second time the timings are run, there is a speed-up of both f and fN. All subsequent timings are like the 2nd. The first timing for f can range from 4 to 10 seconds; the second is around 2 seconds. The results for fN seem to be consistently as shown. Note that if the order is changed to {fNQ, fNNQ, f, fN} in the Table, the timings do not change. All eight integrations evaluate the functions at 165 points (obtained via EvaluationMonitor). The numeric interpolating function fN is much faster than f, but V10 still lags way behind V9. The ?NumericQ functions fNQ and fNNQ are much more competitive with V9, although they are consistently slower. ### Dependence on the number of evaluations Until the number of evaluations grows large enough, the initialization time dominates. Note that the increases in time from the n = 515 to n = 1027 evaluations are about the same in V9 and V10. That suggests that the evaluation of the interpolating function is not the issue. (The results below were run with npts = 0; NIntegrate[f[x], {x, 0, 30 π/400}, Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, Method -> {"GaussKronrodRule", "Points" -> gkpts}}, MaxRecursion -> 0 EvaluationMonitor :> ++npts] // AbsoluteTiming  for gkpts in 1 + {8, 32, 128, 256, 512}.) ## Profiling the initialization of NIntegrate Using SessionTime[] in the code shown down below, we can print the elapsed time for each step in the initialization of the integration of the OP's f[x]. The timings below show that about 2 seconds are taken to start the integration in V10. In V9, it is similar in that the initialization takes up almost all the time it takes to evaluate the integral, namely about 0.13 seconds for initialization. (* V10 time : step 0.000696: UpValue call for: NIntegrateAutomaticStrategy :: {GlobalAdaptive,SymbolicProcessing->0} 2.108061: UpValue call for: NIntegrateGeneralRule ... 2.108530: UpValue call for: NIntegrateGaussKronrodRule ... 2.109095: UpValue call for: NIntegrateGlobalAdaptive ... *)  I don't know how to figure out what happens between the first two steps. One can omit the "SymbolicProcessing" -> 0 option to see that much, if not all, symbolic processing steps are omitted. Even with "SymbolicProcessing" -> 0, NIntegrate symbolically evaluated the integrand (for instance, to see if the integrand is a list, which it has to do even if "SymbolicProcessing" is turned off). ### Profiling code Note the profiling code below changes the definition of internal functions used by NIntegrate. Use Quit[] to reset. Do not execute more than once (nothing really bad will happen; Quit[] will still work and reset). Note the code uses the global symbol tmark to mark the time the computation starts. This has to be set by the user at the beginning of the computation like this: tmark = SessionTime[]; NIntegrate[...]  Profiling code symbNames = Names["NIntegrate*"]; symbNames = Append[Pick[symbNames, StringMatchQ[ symbNames, (__ ~~ "Rule") | (__ ~~ "Global" | "Local" | "MonteCarlo" | "Principal" | "Levin" | "Osc" ~~ ___)]], "NIntegrateAutomaticStrategy"]; symbs = ToExpression[#] & /@ symbNames; dvs = DownValues /@ symbs; uvs = UpValues /@ symbs; Unprotect /@ symbs; dvsNew = MapThread[ With[{s = #2}, DownValues[s] = ReplaceAll[#1, HoldPattern[ a_ :> b___] :> (a :> (Print[SessionTime[] - tmark, ": DownValue call for: ", Style[s, Red]]; b))]] &, {dvs, symbs, symbNames}]; uvsNew = MapThread[ With[{s = #2}, UpValues[s] = ReplaceAll[#1, HoldPattern[Block[vars_, CompoundExpression[b___]]] :> Block[{res = Block[vars, CompoundExpression[b]]}, Print[SessionTime[] - tmark, ": UpValue call for: ", Style[s, Blue], Style[" ::\n", Blue], res]; res] ] ] &, {uvs, symbs, symbNames}];  I used this little utility function to time the integrals for the various integrands. It sets the global symbol tmark to mark the time the computation starts. It is set up for this particular problem. For instance the value of the integral of the function that is being interpolated is -1/100, so the + 1/100 lets me see slight changes in the computed value, if any. ClearAll[timeNI]; SetAttributes[timeNI, HoldAll]; timeNI[fn_, opts___] := Module[{npts}, List[$Version,
tmark = SessionTime[];
npts = 0;
(NIntegrate[fn, {x, 0, 30 π/400}, opts,
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0},
EvaluationMonitor :> ++npts] + 1/100) // AbsoluteTiming,
Row[{npts, " points"}]
]
]

• Thank you very much for your deep analysis! I will check the hints and see how they work with my code. Jan 25 '16 at 16:47
• @Frederico You're welcome. Good luck! Jan 25 '16 at 20:48
• @MichaelE2 If you skip the command symbNames = Append[Pick[symbNames, ...  you might be able to see what is happening between NIntegrateAutomaticStrategy and NIntegrateGeneralRule. (Sorry, I can't spent time today to check it myself.) Impressive work! Jan 25 '16 at 22:22
• @AntonAntonov Thanks. I gave your suggestion a once-through and there's still a 2-sec. gap between NIntegrateAutomaticStrategy and the next call, NIntegrateGetMethodOptionValues. It is all quite fast from that point on. Jan 26 '16 at 1:25
• FWIW, AbsoluteTiming[InterpolatingFunction @@ f] takes about 1.5 sec. So if a reevaluation of f is triggered inside NIntegrate, it would account for at least 75% of the slowdown. Jan 26 '16 at 14:54