NIntegrate over a list of functions

Question

This question is the result of these other two questions. Question 1 and 2. I thought it would be better to ask a new question rather than deleting previous one.

I think When NIntegrate is used on lists it becomes slow. I give the integrand a list of data points and hence integrand is a list of functions. When the number of elements are below 20, Mathematica is either as fast as or faster than Matlab. Yet, when the number of elements increases Mathematica gets slow:

SetSystemOptions["CatchMachineUnderflow" -> False]
shaxis = Table[1.0*i, {i, 1, 20}];
Total[NIntegrate[shaxis/(x^3 + 10), {x, 0, Infinity}, Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}]] // AbsoluteTiming

{0.039000, 54.7079717948704}

laxis = Table[1.0*i, {i, 1, 2046}];
Total[NIntegrate[laxis/(x^3 + 10), {x, 0, Infinity}, Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}]] // AbsoluteTiming

{4.058000, 545537.7346865428}

Both Matlab and Mathematica, in this case, use global adaptive quadrature method. The performance as a function of number of elements of the list depends on the integrand too. In general, it seems it is the length of the list which is the factor making Mathematica slow.

Matlab's documentation about integral:

q = integral(fun,xmin,xmax) numerically integrates function fun from xmin to xmax using global adaptive quadrature and default error tolerances.

Edit

I think the problem is that Matlab uses vectorization and hence uses all of the CPU cores available, but Mathematica won't do the same. If I use a very long list in Matlab, it uses all the cores but not Mathematica.

Matlab's code:

clc
clear all;
f = @(x,t)  t./((x.^3) + 10);
t = 1:0.001:2000;
tic
sum(integral(@(x) f(x,t),0,Inf,'ArrayValued',true))
toc

In Matlab code I used 'ArrayValued' to enable the integrand to receive a vector. I have the same problem, not using all the cores, in this question too. There again Matlab uses vectorization.

Answer 1

NIntegrate does each integral separately

This has been observed before: NIntegrate piecewise vector function, Nested NIntegrate of vector function. It is also clear from the following BenchmarkPlot:

int[n_] := Block[{shaxis},
   shaxis = Table[1.0*i, {i, 1, n}];
   NIntegrate[shaxis/(x^3 + 10), {x, 0, Infinity}, 
    Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}]];
BenchmarkPlot[int, # &, 2^Range[4, 11], "IncludeFits" -> True, 
 TimeConstraint -> 20.`]

I would not know why there has been no interest in developing such functionality for NIntegrate. NDSolve has it.

Using NDSolve instead

NDSolve uses a local estimate of the error to control the step size. It is not going to be as robust as the global adaptive strategy of NIntegrate, but it can do a pretty good job. On the OP's test example, it does quite well, but those functions are similar and relatively easy to integrate.

We need to remap the infinite interval {x, 0, Infinity} to a finite one. The substitution x -> t/(1 - t) pulls it back to {t, 0, 1}. As it is not as good at estimating error, I had to fiddle with PrecisionGoal and AccuracyGoal to get as accurate a solution as possible. This is a drawback to this method. If you set them higher, NDSolveValue complains and stops integration before reaching the end. Below we'll compare it with the default NDSolve settings, NIntegrate timings and the exact solution.

shaxis = Table[1.0*i, {i, 1, 2046}];
ClearAll[df];
df[t_] = shaxis/(x^3 + 10) Dt[x] /. x -> t/(1 - t) /. Dt[t] -> 1 // Simplify;
ndvals = NDSolveValue[
    {a'[t] == df[t], a[0] == ConstantArray[0, Length@shaxis]},
    a[1], {t, 0, 1},
    PrecisionGoal -> 16, AccuracyGoal -> 11]; // AbsoluteTiming
(*   {0.349559, Null}   *)

ndvals0 = NDSolveValue[
    {a'[t] == df[t], a[0] == ConstantArray[0, Length@shaxis]},
    a[1], {t, 0, 1}]; // AbsoluteTiming
(*   {0.240336, Null}   *)

nivals = NIntegrate[
    shaxis/(x^3 + 10), {x, 0, Infinity}]; // AbsoluteTiming
(*   {7.49675, Null}   *)

(* xzczd's solution *)
NIntegrate[shaxis/(x^3 + 10), {x, 0, Infinity}, 
    Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", 
      "SymbolicProcessing" -> 0}]; // AbsoluteTiming
(*   {2.29383, Null}   *)

exact = shaxis Integrate[1/(x^3 + 10), {x, 0, Infinity}];

With the special precision/accuracy settings, the NDSolve solutions are a little better than the default NIntegrate ones. (Again, though, the integrands are not numerically difficult.)

ListPlot[Log10@{Abs[(ndvals - exact)/exact], Abs[(nivals - exact)/exact]}]

The default NDSolve has 8 digits of precision, which is the default goal:

{Norm[(ndvals0 - exact)]/Norm[exact],
 Norm[(ndvals - exact)]/Norm[exact],
 Norm[(nivals - exact)]/Norm[exact]}
(*  {6.57234*10^-8, 1.2117*10^-13, 4.61459*10^-13}  *)

Vector version of NIntegrate

We can easily adapt the implementation of a global adaptive strategy in tutorial/NIntegrateIntegrationStrategies to handle vector integrands. With a little refactoring, we can make take advantage of the built-in vectorization of Mathematica functions. In this implementation, the integrand needs to be Listable. On the integrand df above, it uses four cores (of the 4-core/8-virtual-core i7 on my Mac) typical of Mathematica's vectorized functions.

Of course this circumvents all the other nice features of NIntegrate, but it might be comparable to the Matlab function (of which I'm ignorant).

res = IStrategyGlobalAdaptive[df, {0, 1}, 10^-8]; // AbsoluteTiming
Norm[(res - exact)]/Norm[exact]
(*
  {0.952552, Null}
  1.28022*10^-15
*)

Code:

{absc, weights, errweights} = NIntegrate`GaussKronrodRuleData[5, MachinePrecision];

ClearAll[IRuleEstimate, IStrategyGlobalAdaptive];
IRuleEstimate[f_, a_, b_] :=
  (b - a) {weights, errweights}.Transpose[f[Rescale[absc, {0, 1}, {a, b}]]];

IStrategyGlobalAdaptive[f_, {aArg_, bArg_}, tol_] /; 
   MatrixQ[f[N@{aArg, bArg}], NumericQ] :=
  Module[{integral, error, regions, r1, r2, a = N[aArg], b = N[bArg], c},

   {integral, error} = IRuleEstimate[f, a, b];

   (* boundaries, integral, error *)
   regions = {{a, b, integral, error}};

   While[Norm[error] >= tol*Norm[integral],
    (* splitting of the region with the largest error *)
    {a, b} = regions[[-1, 1 ;; 2]]; c = (a + b)/2;

    (* integration of the left region *)
    {integral, error} = IRuleEstimate[f, a, c];
    r1 = {a, c, integral, error};

    (* integration of the right region *)
    {integral, error} = IRuleEstimate[f, c, b];
    r2 = {c, b, integral, error};

    (* sort the regions: the largest error is last *)
    regions = Join[Most[regions], {r1, r2}];
    regions = SortBy[regions, Norm@*Last];

    (* global integral and error vectors *)
    {integral, error} = Total[regions[[All, -2 ;;]]];
    ];

   integral
   ];

Answer 2

Adding "SymbolicProcessing" -> 0 (it's probably the "default setting" of Matlab, right?) and making use of parallelism gives me a 3X speedup on my dual-core old laptop:

laxis = ParallelTable[1.0 i, {i, 1, 2046}];
Total[ParallelMap[
   NIntegrate[#/(x^3 + 10), {x, 0, Infinity}, 
     Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", 
       "SymbolicProcessing" -> 0}] &, laxis]] // AbsoluteTiming