The problem

Consider some pre-generated list with coordinates x1,x2,x3,x4 and values of some function func at these points. For fixed x1 belonging to the x1 range from list but generally not equal to any of the grid points, I need to interpolate it somehow and integrate it with some pre-factor function prefactor[x1,x2,x3,x4] over x2, x3, x4.

The problem is that a brute-force integral using Interpolation is very slow, caused by interpolation built in Mathematica. I attempted to improve the integration speed by ~5-10 times (depending on the values of the parameter and x1) by keeping an accuracy of ~10% - by "mapping" the data onto a grid if assuming the given value of x1 and summing it over x2,x3,x4 (the code is given below). However, the implementation is very ugly - I am sure it is not optimized and contains rough approximations, while the speedup may be higher.

Could you please help me speed up the integration in a proper way?

Test data

This is the test list:

func[x1_, x2_, x3_, x4_] = If[x3 > x1, Exp[-(Sqrt[x1^2 + x3^2]/5)]*Sin[x2]*(Cos[20*x2])^x4*((100 - x4)/50)^2, 10^-90];
gridx1 = Table[x, {x, 0.05, 5.1, (5.1 - 0.05)/30}];
gridx2 = Table[10^x, {x, -5, Log10[0.05], (Log10[0.05] + 5)/20}];
gridx3 = Table[10^x, {x, Log10[0.051], Log10[350], (Log10[350] - Log10[0.051])/50}];
gridx4 = Table[x, {x, 38., 88., 2.}];
GridIn1 = {gridx1, gridx2, gridx3, gridx4};
GridIn = Tuples[GridIn1];
list = Join[GridIn, Partition[func @@@ GridIn, 1], 2];

Below, there are the interpolation funcInt, the pre-factor prefactor and the integral intv:

funcInt[x1_,x2_,x3_,x4_] = 10^(Interpolation[Log10[list],InterpolationOrder->1][Log10[x1],Log10[x2],Log10[x3],Log10[x4]]);
prefactor[x1_, x2_, x3_, x4_, ct_] = Exp[-x4/(Cos[x2]*ct*x3/x1)]/(Cos[x2]*ct*x3/x1);
intv[x1_, ct_] := NIntegrate[funcInt[x1, x2, x3, x4]*prefactor[x1, x2, x3, x4, ct], {x2, Min[gridx2], Max[gridx2]}, {x3, Min[gridx3], Max[gridx3]}, {x4, Min[gridx4], Max[gridx4]}, Method -> "AdaptiveMonteCarlo"]

Here, ct is a positive parameter with meaningful values between 0.01 and 10^8. Like everything that I do in Mathematica, it is very slow:

intv[4, 30] // AbsoluteTiming

{0.708241, 0.0000166047}

I interpolated over the logarithmized data to accurately integrate into the region where the integrand drops exponentially (the case of small ct).

The potential speedup may be (compared to the integration of the initial function func) 10 times:

intv0[x1_, ct_] := 
  func[x1, x2, x3, x4]*prefactor[x1, x2, x3, x4, ct], {x2, 
   Min[gridx2], Max[gridx2]}, {x3, Min[gridx3], Max[gridx3]}, {x4, 
   Min[gridx4], Max[gridx4]}, Method -> "AdaptiveMonteCarlo"]
intv0[4, 30] // AbsoluteTiming

{0.0670157, 0.0000184466}

It may also be imlroved further, depending of the performance of the built-in integration methods in Mathematica.

My attempt

I made the following attempt. For the given value of x1, I "map" Log10[list] onto a grid {{x1},gridx21,gridx31, gridx41}, where gridx21,gridx31,gridx41 are denser than the initial grid gridx1,.... Then, having the grid I compute the size of the intervals Delta x, multiply the mapped list with the intervals and the prefactor, and then sum over it.

It gives me a $\sim 5$x speedup per list of ct, with ~10% accuracy compared to intv.

However, I am pretty sure that the implementation is very dumb. First, it is not defined correctly mathematically (there are some rough approximations) - instead of using, e.g., Darboux sum $\sum_{i}\text{sup}(f_{x\in \{x_{i}\}})\Delta x_{i}$ I just sum $\sum_{i}f_{x_{i}}\Delta x_{i}$, where $f(x_{i})$ is just a value of the function at a point $x_{i}$, with no adequate relation between $x_{i}$ and $\Delta x_{i}$. Second, in particular, it becomes very inconvenient if there are several lists (with different grids) that I need to interpolate and integrate. I would appreciate it if you comment on this implementation.

So, first, let us define the "mapping" function - one may use any of the beautiful solutions presented in this post:

nd = 4;
cf4 = Module[{xgvars = Unique["xg"] & /@ slist @@ Range@nd, 
     igvars = Unique["ig"] & /@ slist @@ Range@nd, 
     tgvars = Unique["tg"] & /@ slist @@ Range@nd, 
     ivars = Unique["i"] & /@ slist @@ Range@nd, 
     svars = Unique["s"] & /@ slist @@ Range@nd, 
     tvars = Unique["t"] & /@ slist @@ Range@nd, 
     jvars = Unique["j"] & /@ slist @@ Range@nd}, 
       Compile[{seq@{xgvars, _Real, 1}, {y, _Real, nd}}, 
        Module[{seq@igvars, seq@tgvars, seq@ivars, seq@svars, 
         seq[igvars = 
           Floor[xgvars] - UnitStep[xgvars - indexed@Dimensions@y]];
         seq[tgvars = xgvars - igvars];
         Table[seq[tvars = Compile`GetElement[tgvars, jvars]];
          seq[ivars = Compile`GetElement[igvars, jvars]];
          seq[svars = 1. - tvars];
            Times @@ #[[All, 1]] Compile`GetElement[y, 
                Sequence @@ #[[All, 2]]] & /@ 
             Tuples@Transpose[{{svars, ivars}, {tvars, 
                 ivars + 1}}, {2, 3, 1}]], 
          seq@{jvars, Length@igvars}]], CompilationTarget -> "C", 
        RuntimeAttributes -> {Listable}, Parallelization -> True, 
        RuntimeOptions -> "Speed"], Except[seq | eval | indexed]] /. 
      seq[expr_] :> 
       RuleCondition@(Sequence @@ 
             Except[slist | indexed]] /. {l_slist :> l[[i]], 
             indexed[l_] :> Inactive[Compile`GetElement][l, i]}, {i, 
            nd}]) /. eval@expr_ :> 
      RuleCondition[Activate[expr /. slist -> List]]] // Activate;

Next, let us define the out grid and the intervals:

(*Slightly different maximal value of x2max since cf4 gets stuck for x2max=0.05*)
gridx21 = Table[x, {x, 10^-5., 0.05, (0.05 - 10^-5.)/40}];
(*Out grid for x3. Will be computed once x1 is fixed*)
gridx31temp[x1_] := Table[x, {x, x1, 350., (350. - x1)/200}];
gridx41 = Table[x, {x, 38., 88., 0.5}];
DxVals[list_] := 
    list[[i]] - list[[i - 1]], {i, 2, Length[list], 1}], {list[[-1]] -
(*Delta x2 and Delta x4*)
Dx2vals = DxVals[gridx21];
Dx4vals = DxVals[gridx41];
(*Reshaped array of logarithmized values of the function from list. Needs for mapping*)
vals = ArrayReshape[
   Log10[list[[All, 5]]], {Length[gridx1], Length[gridx2], 
    Length[gridx3], Length[gridx4]}];
(*Compiled code which computes the product of Delta x*)
DxValsTotComp = 
  Hold@Compile[{{Dx2vals, _Real, 1}, {Dx3vals, _Real, 
       1}, {Dx4vals, _Real, 1}}, 
     Table[{Dx2vals[[i]]*Dx3vals[[j]]*Dx4vals[[k]]}, {i, 1, 
       Length[Dx2vals]}, {j, 1, Length[Dx3vals]}, {k, 1, 
       Length[Dx4vals]}], CompilationTarget -> "C", 
     RuntimeOptions -> "Speed"] // ReleaseHold;

Next, let us define the block returning the mapped list:

(*Block producing the mapped list*)
BlockDiscreteIntegral[x1_] := Block[{},
  gridx31 = gridx31temp[x1];
  Dx3vals = DxVals[gridx31];
  GridOut1 = Log10[{{x1 // N}, gridx21, gridx31, gridx41}];
  GridOut = 10^Tuples[{GridOut1[[2]], GridOut1[[3]], GridOut1[[4]]}];
  Dxvals = 
   Flatten[DxValsTotComp[Dx2vals, Dx3vals, Dx4vals], {1, 2, 3}];
  (*Mapping routine*)
  xig = MapThread[
    Interpolation[Transpose@{#, Range@Length@#}, 
       InterpolationOrder -> 1][#2] &, {Log10[GridIn1], GridOut1}];
  Join[GridOut, 10^Partition[Flatten@cf4[Sequence @@ xig, vals], 1], 
   Dxvals, 2]

Finally, this is the block returning the integral for the list ctlist:

(*Compiled summed integral*)
IntegralCompiled = 
 Hold@Compile[{{tab, _Real, 2}, {ct, _Real}, {x1, _Real}}, 
         Cos[#[[1]]]*ct*#[[2]]/x1) #[[4]]*#[[5]]) & /@ tab], 
    CompilationTarget -> "C", RuntimeOptions -> "Speed"] // ReleaseHold
BlockIntegralList[x1_, ctlist_] := Block[{},
  dat = BlockDiscreteIntegral[x1];
  Table[{ctlist[[i]], IntegralCompiled[dat, ctlist[[i]], x1]}, {i, 1, 
    Length[ctlist], 1}]

Let us compare the integral and the discrete integral:

mtest = 4;
ctlist = {0.5, 1, 2, 3, 4, 5, 6, 10, 50, 100, 1000, 10^4} // N;
t1 = BlockIntegralList[mtest, ctlist]; // AbsoluteTiming // First
t2 = Table[{ctlist[[i]], intv[mtest, ctlist[[i]]]}, {i, 1, 
      Length[ctlist], 1}]; // AbsoluteTiming // First
t2[[All, 2]]^-1 (t1[[All, 2]] - t2[[All, 2]])



{0.0916109, 0.0606622, 0.0491823, 0.0590111, 0.0686524, 0.113182,
0.0597432, 0.0117492, -0.0875845, -0.0977589, -0.032049, -0.0811554}

Another approach would be to generate values of the interpolated function for random points inside the integration domain and sum them. In principle, I may use cf4 for this by calling it for any particular point, but the problem is that it is not optimized for this purpose.


2 Answers 2


By inspection the Transpose of your list consists of an immense data set of periodic step functions.

First impression

In[32]:= Dimensions[list\[Transpose]]
Out[32]= {5,863226}

In[49]:= Mean/@(list\[Transpose])

Out[49]= {2.575,0.00686475,42.3706,63.,0.0000545084}

In[71]:= GraphicsGrid[{(

enter image description here

Its not possible to find an interpolation by splines, of course. A blind interpolation will end up as a more complicated nested call of the original construction. Then its probably better to dress up the construction function with a linear interpolation between adjacent points in the grid.

Just take the Transpose, get the step parts and interpolate them as pure functions. Then one can define a function distributing the arguments at will over the grid after of pure functions of one variable.

If the periodicity is exact, one may reduce the variables with a Mod to the single representative interval.

For construction work with toy dimensions of small integers and try Fourier and Tensor symmetry checks and to find the most efficient storage and access method.

Since any integral is in fact an orderless sum over values times the volumes sum dx occupied by points of the same value y (Riemann sum), for very irregular functions the Lebesgue integral is helpful:

Sort the values in intervals of value in (y,y+dy), and add the volumes Gathered in the parameter space V=f^-1 (y,y,+dy).

The great advantage of the discrete Lebesgue integral approximation:

There i no local rounding to zero of inexact values or unlucky sampling; instead positive and negative values are gathered separately and cancellation is an absolute effect convergent (or not) in the limit dy-> 0

  • $\begingroup$ Thanks! What do you think about my attempt to speedup the integration? $\endgroup$ Apr 18, 2023 at 11:32

On top of what @John Taylor said, I would parallize this part:

t2 = Table[{ctlist[[i]], intv[mtest, ctlist[[i]]]}, {i, 1, 
      Length[ctlist], 1}]; // AbsoluteTiming // First

to this:

t2 = ParallelTable[{ctlist[[i]], intv[mtest, ctlist[[i]]]}, {i, 1, 
      Length[ctlist], 1}]; // AbsoluteTiming // First

I doubt there is more of a generaliseable approach that can be taken with the integration itself. My attempts at diving deeper into tweaking the accuracy were in vain.


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.