# How to produce the following table faster?

Question

Consider a toy function func[x1,x2,x3,x4] and a toy grid gridTot of the coordinates x1,x2,x3,x4:

func[x1_, x2_, x3_, x4_] =
If[x3 > x1,
1/(x1^2 + 3^2) Exp[-(x1^2 + x3^2)^(1/3.)]*x3/(x3^2 + 20^2) Sin[x2]*
Cos[10 x2]*x4/(x4 + 50), 0];
gridx1 = Table[x1, {x1, 0.05, 5.1, (5.1 - 0.05)/30}];
gridx2 = Table[20^x2, {x2, -5., Log10[0.05], (Log10[0.05] + 5)/20}];
gridx3 = Table[10^
x3, {x3, Log10[0.05], Log10[350.], (Log10[350] - Log10[0.05])/50}];
gridx4 = Table[x4, {x4, 38, 88, 2.}];
gridTot = Tuples[{gridx1, gridx2, gridx3, gridx4}];


I may produce the table x1,x2,x3,x4,func[x1,x2,x3,x4], where x1,x2,x3,x4 belong to gridTot fast using Compile:

compiled =
Hold@Compile[{{gridTot, _Real, 2}},
Table[{gridTot[[i]][[1]], gridTot[[i]][[2]], gridTot[[i]][[3]],
gridTot[[i]][[4]],
func[gridTot[[i]][[1]], gridTot[[i]][[2]], gridTot[[i]][[3]],
gridTot[[i]][[4]]]}, {i, 1, Length[gridTot], 1}],
CompilationTarget -> "C", RuntimeOptions -> "Speed"] /.
DownValues@func // ReleaseHold;
functab = compiled[gridTot]; // AbsoluteTiming


{0.134968,Null}

Let us now assume that we have the table functab instead of func, and we want to interpolate it (obtaining funcInt[x1,x2,x3,x4], to produce the table of values x1,x2,x3,x4,funcint[x1,x2,x3,x4] for a different grid gridTot1:

funcInt[x1_, x2_, x3_, x4_] =
Interpolation[functab, InterpolationOrder -> 1][x1, x2, x3, x4];
gridx11 = Table[x1, {x1, 0.051, 5.09, (5.09 - 0.051)/30}];
gridx21 =
Table[20^x2, {x2, -4.9, Log10[0.049], (Log10[0.049] + 4.9)/20}];
gridx31 =
Table[10^x3, {x3, Log10[0.051], Log10[349.], (
Log10[349.] - Log10[0.051])/50}];
gridx41 = Table[x4, {x4, 38.1, 79.9, 2.}];
gridTot1 =
Flatten[Table[{gridx11[[i]], gridx21[[j]], gridx31[[k]],
gridx41[[l]]}, {i, 1, Length[gridx11]}, {j, 1,
Length[gridx21]}, {k, 1, Length[gridx31]}, {l, 1,
Length[gridx41]}], {1, 2, 3, 4}];


Now, Compile does not work since we are dealing with the interpolation:

compiled1 =
Hold@Compile[{{gridTot, _Real, 2}},
Table[funcInt[gridTot[[i]][[1]], gridTot[[i]][[2]],
gridTot[[i]][[3]], gridTot[[i]][[4]]], {i, 1, Length[gridTot],
1}], CompilationTarget -> "C", RuntimeOptions -> "Speed"] /.
DownValues@funcInt // ReleaseHold;
functab1 = compiled1[gridTot1]; // AbsoluteTiming


just gets stuck. If removing , CompilationTarget -> "C", RuntimeOptions -> "Speed", it works, but very slow:

functab1 = compiled1[gridTot1]; // AbsoluteTiming


{6.4564,Null}

This means that Compile does not handle the interpolation properly. Actually, there is almost no improvement compared to the ordinary code:

functab1 =
Table[funcInt[gridTot[[i]][[1]], gridTot[[i]][[2]],
gridTot[[i]][[3]], gridTot[[i]][[4]]], {i, 1, Length[gridTot],
1}]; // AbsoluteTiming


{9.9496,Null}

Question

Could you please tell me whether there are efficient ways to speedup producing functab1? I need to speedup this procedure since in a real case I am dealing with a much larger grid and a much slower interpolated function.

Update

I have tried to use cf2 function for my example:

vals = functab[[All, 5]];
xiset1[list_] := Sort[RandomReal[MinMax[list], Length[list]]]
x1set1 = xiset1[gridx1];
x2set1 = xiset1[gridx2];
x3set1 = xiset1[gridx3];
x4set1 = xiset1[gridx4];
xsettot =
Flatten[Table[{x1set1[[i]], x2set1[[j]], x3set1[[k]],
x4set1[[l]]}, {i, 1, Length[x1set1]}, {j, 1, Length[x2set1]}, {k,
1, Length[x3set1]}, {l, 1, Length[x4set1]}], {1, 2, 3, 4}];


and

threadCount =
"ParallelThreadNumber" /. ("ParallelOptions" /.
SystemOptions["ParallelOptions"]);
cresult2 =
Join @@ cf2[gridx1, gridx2, gridx3, gridx4, vals, xsettot,


But it got stuck.

Update 2

In case of cf2, if I want to make the same function for 3D grid, should I just replace Plus[...] with

Plus[s1 s2 s3 CompileGetElement[y, i1, i2, i3],
s1 s2 t3 CompileGetElement[y, i1, i2, 1 + i3],
s1 t2 s3 CompileGetElement[y, i1, 1 + i2, i3],
t1 s2 s3 CompileGetElement[y, 1 + i1, i2, i3],
s1 t2 t3 CompileGetElement[y, i1, 1 + i2, 1 + i3],
t1 s2 t3 CompileGetElement[y, 1 + i1, i2, 1 + i3],
t1 t2 s3 CompileGetElement[y, 1 + i1, 1 + i2, i3],
t1 t2 t3 CompileGetElement[y, 1 + i1, 1 + i2,
1 + i3]]/((s1 + t1) (s2 + t2) (s3 + t3))
, {k, begin, end}]


and remove the block

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g4];
a = 1;
b = n;
ga = CompileGetElement[g4, a];
gb = CompileGetElement[g4, b];
z = CompileGetElement[x, k, 4];
If[z <= ga,
i4 = 1;
t4 = ga;
s4 = CompileGetElement[g4, i4 + 1];
,
If[z >= gb,
i4 = n - 1;
t4 = CompileGetElement[g4, i4];
s4 = gb;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g4, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i4 = a; t4 = z - ga; s4 = gb - z;
];
];
];


apart from removing obvious variables g4 etc.?

Update 3

I have changed cf2 to interpolate over 3D arrays:

D3toD3 = Compile[{{g1, _Real, 1}, {g2, _Real, 1}, {g3, _Real,
1}, {y, _Real, 3}, {x, _Real,
Module[{jobCount, begin, end, q, r, i1, i2, i3, t1, t2, t3, s1, s2,
s3}, jobCount = Length[x];
q = Quotient[jobCount, threadCount];
r = Mod[jobCount, threadCount];
begin =
q (threadID - 1) + Quotient[r (threadID - 1), threadCount] + 1;
Table[Block[{a, ga, b, gb, c, gc, z, n}, n = Length[g1];
a = 1;
b = n;
ga = CompileGetElement[g1, a];
gb = CompileGetElement[g1, b];
z = CompileGetElement[x, k, 1];
If[z <= ga, i1 = 1;
t1 = ga;
s1 = 0.;, If[z >= gb, i1 = n - 1;
t1 = 0.;
s1 = gb;, While[b - a > 1, c = Quotient[a + b, 2];
gc = CompileGetElement[g1, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];];
i1 = a; t1 = z - ga; s1 = gb - z;];];];
Block[{a, ga, b, gb, c, gc, z, n}, n = Length[g2];
a = 1;
b = n;
ga = CompileGetElement[g2, a];
gb = CompileGetElement[g2, b];
z = CompileGetElement[x, k, 2];
If[z <= ga, i2 = 1;
t2 = ga;
s2 = 0.;, If[z >= gb, i2 = n - 1;
t2 = 0.;
s2 = gb;, While[b - a > 1, c = Quotient[a + b, 2];
gc = CompileGetElement[g2, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];];
i2 = a; t2 = z - ga; s2 = gb - z;];];];
Block[{a, ga, b, gb, c, gc, z, n}, n = Length[g3];
a = 1;
b = n;
ga = CompileGetElement[g3, a];
gb = CompileGetElement[g3, b];
z = CompileGetElement[x, k, 3];
If[z <= ga, i3 = 1;
t3 = ga;
s3 = 0.;, If[z >= gb, i3 = n - 1;
t3 = 0.;
s3 = gb;, While[b - a > 1, c = Quotient[a + b, 2];
gc = CompileGetElement[g3, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];];
i3 = a; t3 = z - ga; s3 = gb - z;];];];
Plus[s1 s2 s3 CompileGetElement[y, i1, i2, i3],
s1 s2 t3 CompileGetElement[y, i1, i2, 1 + i3],
s1 t2 s3 CompileGetElement[y, i1, 1 + i2, i3],
t1 s2 s3 CompileGetElement[y, 1 + i1, i2, i3],
s1 t2 t3 CompileGetElement[y, i1, 1 + i2, 1 + i3],
t1 s2 t3 CompileGetElement[y, 1 + i1, i2, 1 + i3],
t1 t2 s3 CompileGetElement[y, 1 + i1, 1 + i2, i3],
t1 t2 t3 CompileGetElement[y, 1 + i1, 1 + i2,
1 + i3]]/((s1 + t1) (s2 + t2) (s3 + t3))
, {k, begin, end}]], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"];


Consider a test data temp_distr.dat. I define in/out grids InGrid, OutGrid (the points of the latter belongs to the domain covered by the former) and built-in interpolation:

SetDirectory[NotebookDirectory[]];
data = Import["temp-distr.dat", "Table"];
InGridx1 = DeleteDuplicates[data[[All, 1]]];
InGridx2 = DeleteDuplicates[data[[All, 2]]];
InGridx3 = DeleteDuplicates[data[[All, 3]]];
vals = ArrayReshape[
data[[All, 4]], {Length[InGridx1], Length[InGridx2],
Length[InGridx3]}];
OutGridi[list_, n_] :=
Table[x, {x, Min[list], Max[list], (Max[list] - Min[list])/n}]
OutGridx1 = OutGridi[InGridx1, 50];
OutGridx2 = OutGridi[InGridx2, 50];
OutGridx3 = OutGridi[InGridx3, 130];
OutGrid = Tuples[{OutGridx1, OutGridx2, OutGridx3}];
FuncInt[x1_, x2_, x3_] =
Interpolation[data, InterpolationOrder -> 1][x1, x2, x3];


Finally, I compare the grids produced with FuncInt and D3toD3:

threadCount =
"ParallelThreadNumber" /. ("ParallelOptions" /.
SystemOptions["ParallelOptions"]);
outlist1 =
Join @@ D3toD3[InGridx1, InGridx2, InGridx3, vals, OutGrid,
outlist2 = FuncInt @@@ OutGrid; // AbsoluteTiming


The absolute difference may be very large:

Max[Abs[outlist1 - outlist2]]


0.0122018

Edit

Henrik Schumacher has developed the fastest approach, but it requires additional software and hence cannot be used "from box". The approach developed by Lukas Lang is slower but does not require additional soft. For different people, any of them may be preferable.

• Did you check which step is slowing it down? Is it the call to Table? Or is it Flatten? Or is it func? Apr 14, 2023 at 17:57
• @Zatrapilla : it is func. Apr 14, 2023 at 18:03
• I think vals has the wrong shape. It should be an array of size {Length[gridx1],Length[gridx2],Length[gridx3],Length[gridx4]}. Apr 15, 2023 at 10:34
• We have gridTot == Tuples[{gridx1, gridx2, gridx3, gridx4}], right? Then vals = ArrayReshape[functab[[All,5]],{Length[gridx1],Length[gridx2],Length[gridx3],Length[gridx4]}] should work. Apr 15, 2023 at 10:47
• Or vals = Outer[func, gridx1, gridx2, gridx3, gridx4];. Apr 15, 2023 at 10:49

Here's an attempt to improve further upon @HenrikSchumacher's excellent answer:

nd = 4;
cf3 = 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},
Inactivate[
Compile[{seq@{xgvars, _Real, 1}, {y, _Real, nd}},
Module[{seq@igvars, seq@tgvars, seq@ivars, seq@svars,
seq@tvars},
seq[igvars =
Floor[xgvars] - UnitStep[xgvars - indexed@Dimensions@y]];
seq[tgvars = xgvars - igvars];
Table[seq[tvars = CompileGetElement[tgvars, jvars]];
seq[ivars = CompileGetElement[igvars, jvars]];
seq[svars = 1. - tvars];

eval@Total[
Times @@ #[[All, 1]] CompileGetElement[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 @@
Table[Inactivate[expr,
Except[slist | indexed]] /. {l_slist :> l[[i]],
indexed[l_] :> Inactive[CompileGetElement][l, i]}, {i,
nd}]) /. eval@expr_ :>
RuleCondition[Activate[expr /. slist -> List]]] // Activate;

dims = {30, 20, 50, 30};
g = {Sort[RandomReal[{-1, 1}, dims[[1]]]],
Sort[RandomReal[{-1, 1}, dims[[2]]]],
Sort[RandomReal[{-1, 1}, dims[[3]]]],
Sort[RandomReal[{-1, 1}, dims[[4]]]]};
vals = RandomReal[{-1, 1}, dims];

f = Interpolation[Transpose[{Tuples[g], Flatten[vals]}],
InterpolationOrder -> 1];

n = 25;
xg = Sort@RandomReal[{#[[1]], #[[-1]]}, n] & /@ g;

result = f @@@ Tuples@xg; // AbsoluteTiming // First
cresult = cf[Sequence @@ g, vals, Tuples@xg]; //
AbsoluteTiming // First
"ParallelThreadNumber" /. ("ParallelOptions" /.
SystemOptions["ParallelOptions"]);
c2result =
Join @@ cf2[Sequence @@ g, vals, Tuples@xg, Range[threadCount],
threadCount]; // AbsoluteTiming // First
Max[Abs[result - c2result]]
Interpolation[
Transpose@{#, Range@Length@#},
InterpolationOrder -> 1
][#2] &,
{g, xg}
];
c3result = Flatten@cf3[
Sequence @@ xig,
vals
]; // AbsoluteTiming // First
Max[Abs[result - c3result]]
(* 5.06924 *)
(* 0.0531016 *)
(* 0.0504067 *)
(* 5.71765*10^-15 *)
(* 0.0152314 *)
(* 5.71765*10^-15 *)


As you can see, cf3 is even faster than cf and cf2 from @HenrikSchumacher's answer. It also has the advantage that changing the number of dimensions of the interpolation is simply done by changing nd (the code takes care to generate the appropriate expression to be compiled). This version epxloits the fact that if we want to interpolate a regular grid g onto a regular grid xg, we can do the step of finding the appropriate fractional grid indices xig only once for each grid axis. The compiled function then gets the lists of fractional grid indices and generates the table inside.

• Thanks! But my attempt in the question body produces the result within the same time. Apr 14, 2023 at 17:12
• @JohnTaylor Hmm, I would have thought it should have made more of a difference... I'll try to see whether I can think of something faster. In the meantime: Why exactly do you need it to be faster? Do you need to run a similar computation repeatedly? If yes, what changes from run to run? Apr 14, 2023 at 21:14
• @JohnTaylor Please see the updated answer for an improved version of HenrikSchumacher's answer Apr 15, 2023 at 12:03
• @JohnTaylor I have updated the answer so that cf3 now works properly with e.g. xg=g. The issue was caused by the interpolation trying to access the grid point after the last, which is now avoided Apr 20, 2023 at 12:21
• For those kinds of grids, @HenrikSchumacher's solution is probably better: I'd expect the speedup due to the compiled index computation to outweigh the slowdown due to the repeated computation of the first index. You could of course try to combine the two strategies for maximum speed Apr 22, 2023 at 11:03

tl;dr: As always with my post, find the most recent implementation at the bottom of the post.

InterpolationFunctions are indeed much less efficient than they could/should be. Here is a compiled implementation of a scalar interpolation function in four variables with linear interpolation in each direction. It is designed for random access and for rectangular, but irregularly spaced grids. It finds the correct grid point via binary search in each dimension. Moreover, it is Listable and parallelized.

cf = Compile[{{g1, _Real, 1}, {g2, _Real, 1}, {g3, _Real, 1}, {g4, _Real, 1}, {y, _Real, 4}, {x, _Real, 1}},
Module[{i1, i2, i3, i4, t1, t2, t3, t4, s1, s2, s3, s4},

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g1];
a = 1; ga = CompileGetElement[g1, a];
b = n; gb = CompileGetElement[g1, b];
z = CompileGetElement[x, 1];
If[z <= ga,
i1 = 1; t1 = 0.;
,
If[z >= gb,
i1 = n - 1; t1 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g1, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i1 = a; t1 = (z - ga)/(gb - ga);
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g2];
a = 1; ga = CompileGetElement[g2, a];
b = n; gb = CompileGetElement[g2, b];
z = CompileGetElement[x, 2];
If[z <= ga,
i2 = 1; t2 = 0.;
,
If[z >= gb,
i2 = n - 1; t2 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g2, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i2 = a; t2 = (z - ga)/(gb - ga);
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g3];
a = 1; ga = CompileGetElement[g3, a];
b = n; gb = CompileGetElement[g3, b];
z = CompileGetElement[x, 3];
If[z <= ga,
i3 = 1; t3 = 0.;
,
If[z >= gb,
i3 = n - 1; t3 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g3, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i3 = a; t3 = (z - ga)/(gb - ga);
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g4];
a = 1; ga = CompileGetElement[g4, a];
b = n; gb = CompileGetElement[g4, b];
z = CompileGetElement[x, 4];
If[z <= ga,
i4 = 1; t4 = 0.;
,
If[z >= gb,
i4 = n - 1; t4 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g4, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i4 = a; t4 = (z - ga)/(gb - ga);
];
];
];

s1 = (1. - t1);
s2 = (1. - t2);
s3 = (1. - t3);
s4 = (1. - t4);

Plus[
s1 s2 s3 s4 CompileGetElement[y, i1, i2, i3, i4],
s1 s2 s3 t4 CompileGetElement[y, i1, i2, i3, 1 + i4],
s1 s2 t3 s4 CompileGetElement[y, i1, i2, 1 + i3, i4],
s1 s2 t3 t4 CompileGetElement[y, i1, i2, 1 + i3, 1 + i4],
s1 t2 s3 s4 CompileGetElement[y, i1, 1 + i2, i3, i4],
s1 t2 s3 t4 CompileGetElement[y, i1, 1 + i2, i3, 1 + i4],
s1 t2 t3 s4 CompileGetElement[y, i1, 1 + i2, 1 + i3, i4],
s1 t2 t3 t4 CompileGetElement[y, i1, 1 + i2, 1 + i3, 1 + i4],
t1 s2 s3 s4 CompileGetElement[y, 1 + i1, i2, i3, i4],
t1 s2 s3 t4 CompileGetElement[y, 1 + i1, i2, i3, 1 + i4],
t1 s2 t3 s4 CompileGetElement[y, 1 + i1, i2, 1 + i3, i4],
t1 s2 t3 t4 CompileGetElement[y, 1 + i1, i2, 1 + i3, 1 + i4],
t1 t2 s3 s4 CompileGetElement[y, 1 + i1, 1 + i2, i3, i4],
t1 t2 s3 t4 CompileGetElement[y, 1 + i1, 1 + i2, i3, 1 + i4],
t1 t2 t3 s4 CompileGetElement[y, 1 + i1, 1 + i2, 1 + i3, i4],
t1 t2 t3 t4 CompileGetElement[y, 1 + i1, 1 + i2, 1 + i3, 1 + i4]
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];


Here is a usage and test example:

dims = {30, 20, 50, 30};
g = Sort[RandomReal[{-1, 1}, #]] & /@ dims;
vals = RandomReal[{-1, 1}, dims];

f = Interpolation[
Transpose[{Tuples[g, Flatten[vals]}],
InterpolationOrder -> 1];

n = 1000000;
x = Transpose[RandomReal[{#[[1]], #[[-1]]}, n] & /@ g];

result = f @@@ x; // AbsoluteTiming // First
cresult = cf[Sequence @@ g, vals, x]; // AbsoluteTiming // First
Max[Abs[result - cresult]]

8.33928

0.100488

7.77156*10^-16


So it seems to produce essentially the correct results with an 80-fold speedup.

This is certainly not optimal. For one, the fine grained parallelization probably leads to some avoidable overhead. Moreover, as OP wants to apply this interpolation to a presorted list of points, one could first search a neighborhood of the previously found grid point and do the binary search only if this is not successful. Quite likely, this would be way more cache friendly and thus faster.

Edit

Here is an improved version of the above. It gets rid of 3 out of 4 divisions (which does not seem to matter much) and has improved parallelism.

cf2 = Compile[{{g1, _Real, 1}, {g2, _Real, 1}, {g3, _Real, 1}, {g4, _Real, 1}, {y, _Real, 4}, {x, _Real, 2}, {threadID, _Integer}, {threadCount, _Integer}},
Module[{jobCount, begin, end, q, r, i1, i2, i3, i4, t1, t2, t3, t4, s1, s2, s3, s4},

jobCount = Length[x];
q = Quotient[jobCount, threadCount];
r = Mod[jobCount, threadCount];
begin = q (threadID - 1) + Quotient[r (threadID - 1), threadCount] + 1;

Table[

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g1];
a = 1;
b = n;
ga = CompileGetElement[g1, a];
gb = CompileGetElement[g1, b];
z = CompileGetElement[x, k, 1];
If[z <= ga,
i1 = a;
s1 = 1.;
t1 = 0.;
,
If[z >= gb,
i1 = b - 1;
s1 = 0.;
t1 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g1, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i1 = a; t1 = z - ga; s1 = gb - z;
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g2];
a = 1;
b = n;
ga = CompileGetElement[g2, a];
gb = CompileGetElement[g2, b];
z = CompileGetElement[x, k, 2];
If[z <= ga,
i2 = a;
s2 = 1.;
t2 = 0.;
,
If[z >= gb,
i2 = b - 1;
s2 = 0.;
t2 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g2, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i2 = a; t2 = z - ga; s2 = gb - z;
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g3];
a = 1;
b = n;
ga = CompileGetElement[g3, a];
gb = CompileGetElement[g3, b];
z = CompileGetElement[x, k, 3];
If[z <= ga,
i3 = a;
s3 = 1.;
t3 = 0.;
,
If[z >= gb,
i3 = b - 1;
s3 = 0.;
t3 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g3, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i3 = a; t3 = z - ga; s3 = gb - z;
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g4];
a = 1;
b = n;
ga = CompileGetElement[g4, a];
gb = CompileGetElement[g4, b];
z = CompileGetElement[x, k, 4];
If[z <= ga,
i4 = a;
s4 = 1.;
t4 = 0.;
,
If[z >= gb,
i4 = b - 1;
s4 = 0.;
t4 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g4, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i4 = a; t4 = z - ga; s4 = gb - z;
];
];
];
Plus[
s1 s2 s3 s4 CompileGetElement[y, i1, i2, i3, i4],
s1 s2 s3 t4 CompileGetElement[y, i1, i2, i3, 1 + i4],
s1 s2 t3 s4 CompileGetElement[y, i1, i2, 1 + i3, i4],
s1 s2 t3 t4 CompileGetElement[y, i1, i2, 1 + i3, 1 + i4],
s1 t2 s3 s4 CompileGetElement[y, i1, 1 + i2, i3, i4],
s1 t2 s3 t4 CompileGetElement[y, i1, 1 + i2, i3, 1 + i4],
s1 t2 t3 s4 CompileGetElement[y, i1, 1 + i2, 1 + i3, i4],
s1 t2 t3 t4 CompileGetElement[y, i1, 1 + i2, 1 + i3, 1 + i4],
t1 s2 s3 s4 CompileGetElement[y, 1 + i1, i2, i3, i4],
t1 s2 s3 t4 CompileGetElement[y, 1 + i1, i2, i3, 1 + i4],
t1 s2 t3 s4 CompileGetElement[y, 1 + i1, i2, 1 + i3, i4],
t1 s2 t3 t4 CompileGetElement[y, 1 + i1, i2, 1 + i3, 1 + i4],
t1 t2 s3 s4 CompileGetElement[y, 1 + i1, 1 + i2, i3, i4],
t1 t2 s3 t4 CompileGetElement[y, 1 + i1, 1 + i2, i3, 1 + i4],
t1 t2 t3 s4 CompileGetElement[y, 1 + i1, 1 + i2, 1 + i3, i4],
t1 t2 t3 t4 CompileGetElement[y, 1 + i1, 1 + i2, 1 + i3, 1 + i4]
]/((s1 + t1) (s2 + t2) (s3 + t3) (s4 + t4)),
{k, begin, end}
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];


It has a slightly different syntax; call it like this:

threadCount = "ParallelThreadNumber" /. ("ParallelOptions" /. SystemOptions["ParallelOptions"]);
cresult2 =  Join @@ cf2[Sequence @@ g, vals, x, Range[threadCount], threadCount];


It is roughly twice as fast as the preceding version, thus 160 times as fast as the InterpolationFunction. But still, 80% of runtime seem to go into the binary search...

Edit 3

This is probably over the top, but I was curious and implemented the above approach based on binary search directly in C++. Using LibraryLink as a jit-compiler, it can handle arbitrary-dimensional multivariate functions (and also vector-valued functions).

If you want to run this, then you might have to adjust the path dirOpenMP to your local installation of OpenMP and the settings for "CompileOptions".

Needs["CCompilerDriver"]

(*Path to local OpenMP installation. This should work for OpenMP installed via homebrew on Apple Silicon.*)
dirOpenMP = "/opt/homebrew/opt/libomp";

Quiet[Cases[DownValues[cInterpolate], HoldPattern[_ :> x_LibraryFunction] :> LibraryFunctionUnload[x]]];
ClearAll[cInterpolate];
cInterpolate[DomDim_Integer, AmbDim_Integer] :=
cInterpolate[DomDim, AmbDim] = Module[{lib, code, name},
name = "cInterpolate_" <> IntegerString[DomDim] <> "_" <>
IntegerString[AmbDim];
code = StringJoin["

#include \"WolframLibrary.h\"

#include <cstdint>
#include <vector>
#include <array>

#include \"WolframLibrary.h\"

template<typename Int, typename Int1, typename Int2>
Int JobPointer( const Int job_count, const Int1 thread_count, const Int2 k )
{
return job_count/static_cast<Int>(thread_count)*static_cast<Int>(k) + \
}

template< size_t DomDim, size_t AmbDim, typename Scal_in_, typename Scal_out_>
class Interpolator
{
// DomDim   - domain dimension
// AmbDim   - ambient dimension (dimension of the target space)
// Scal_in_ - scalar type for the inputs
// Scal_in_ - scalar type for the outputs

static_assert(DomDim >= 1);

static_assert(AmbDim >= 1);

public:

using Int           = size_t;
using Scal_in       = Scal_in_;
using Scal_out      = Scal_out_;

using Index_T       = std::array<Int,DomDim>;

using KnotVector_T  = const Scal_in *;
using Grid_T        = std::array<KnotVector_T,DomDim>;

using Input_T       = std::array<Scal_in, DomDim>;
using Output_T      = std::array<Scal_out,AmbDim>;

static constexpr Int corner_count = (1 << DomDim);

protected:

Index_T dims;
Grid_T  grid;

const   Scal_out * const      values;

mutable Index_T               grid_idx = {0};

mutable std::array<Input_T,2> t = {{{0},{1}}};
mutable Input_T               x = {0};
mutable Output_T              y = {0};

public:

Interpolator( const Grid_T & grid_, const Index_T & dims_, const Scal_out * const values_ )
:   grid(grid_)
,   dims(dims_)
,   values(values_)
{}

template< typename I>
Interpolator( const Scal_in * const * const grid_, const I * const dims_, const Scal_out * const values_ )
:   values(values_)
{
std::copy( &dims_[0], &dims_[DomDim], &dims[0] );
std::copy( &grid_[0], &grid_[DomDim], &grid[0] );
}

void Evaluate( const Scal_in * const x_, Scal_out * const y_ )
{
Find();
Eval();
Write(y_);
}

const Index_T & GridIndex() const
{
return grid_idx;
}

Int GlobalIndex( const Index_T & idx ) const
{
if constexpr( DomDim == 1 )
{
return idx[0];
}
else
{
Int global_idx = 0;

for( Int k = 0; k < DomDim - 1; ++k )
{
global_idx = (global_idx + idx[k]) * dims[k+1];
}

global_idx += idx[DomDim-1];

return global_idx;
}
}

void Read( const Scal_in * const x_ )
{
std::copy( &x_[0], &x_[DomDim], &x[0] );
}

void Write( Scal_out * const y_ ) const
{
std::copy( &y[0], &y[AmbDim], &y_[0] );
}

void Eval()
{
Scal_in volume = static_cast<Scal_in>(1);

for( Int k = 0; k < DomDim; ++k )
{
volume *= t[0][k] + t[1][k];
}

const Scal_in volume_inv = static_cast<Scal_in>(1) / volume;

std::fill( &y[0], &y[AmbDim], static_cast<Scal_out>(0) );

for( Int kdx = 0; kdx < corner_count; ++kdx )
{
// For each corner of the grid cuboid, we have to find the global position global_idx of the corner in the grid and its coeffcient coeff

Index_T idx = grid_idx;

Scal_in coeff = volume_inv;

for( Int k = 0; k < DomDim; ++k )
{
Int bit  = (kdx >> (DomDim-1-k)) & 1;
idx[k]  += bit;
coeff   *= t[bit][k];
}

const Int global_idx = AmbDim * GlobalIndex(idx);

auto coeff_converted = static_cast<Scal_out>(coeff);

for( Int k = 0; k < AmbDim; ++k )
{
y[k] += coeff_converted * values[global_idx+k];
}
}
}

void Find( const Scal_in z, const Int k )
{
KnotVector_T g = grid[k];

Int a = 0;
Int b = dims[k]-1;

Scal_in g_a = g[a];
Scal_in g_b = g[b];

if( z < g_a )
{
grid_idx[k] = a;
t[0][k]     = static_cast<Scal_in>(1);
t[1][k]     = static_cast<Scal_in>(0);

}
else if ( z > g_b )
{
grid_idx[k] = b-1;
t[0][k]     = static_cast<Scal_in>(0);
t[1][k]     = static_cast<Scal_in>(1);
}
else
{
// Binary search
while( b > a + 1 )
{
Int     c   = a + (b-a)/2;
Scal_in g_c = g[c];

if( z <= g_c )
{
b   = c;
g_b = g_c;
}
else
{
a   = c;
g_a = g_c;
}
}

grid_idx[k] = a;
t[0][k]     = g_b - z;
t[1][k]     = z - g_a;

}
}

void Find()
{
for( Int k = 0; k < DomDim; ++k )
{
Find( x[k], k );
}
}

}; // class Interpolator

EXTERN_C DLLEXPORT int " <> name <>
"(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
constexpr size_t DomDim = ", IntegerString[DomDim], ";
constexpr size_t AmbDim = ", IntegerString[AmbDim], ";

std::vector<mint>   dims (DomDim);
std::vector<mreal*> grid (DomDim);

for( size_t k = 0; k < DomDim; ++k )
{
MTensor g = MArgument_getMTensor(Args[k]);

dims[k] = libData->MTensor_getDimensions(g)[0];
grid[k] = libData->MTensor_getRealData(g);
}

MTensor vals = MArgument_getMTensor(Args[DomDim  ]);
MTensor x_   = MArgument_getMTensor(Args[DomDim+1]);
MTensor y_   = MArgument_getMTensor(Args[DomDim+2]);

size_t thread_count = MArgument_getInteger(Args[DomDim+3]);

const mreal * const x = libData->MTensor_getRealData(x_);
mreal * const y = libData->MTensor_getRealData(y_);

const size_t n = libData->MTensor_getDimensions(x_)[0];

#pragma omp parallel for num_threads( thread_count )
{
Interpolator<DomDim,AmbDim,mreal,mreal> f ( &grid[0], &dims[0], libData->MTensor_getRealData(vals) );

const size_t i_begin = JobPointer( n, thread_count, thread   );
const size_t i_end   = JobPointer( n, thread_count, thread+1 );

for( size_t i = i_begin; i < i_end; ++i )
{
f.Evaluate( &x[DomDim*i], &y[AmbDim*i] );
}
}

libData->MTensor_disown(y_);

return LIBRARY_NO_ERROR;
}"];
lib = CreateLibrary[code, name, "Language" -> "C++",
"ShellOutputFunction" -> Print,(*"ShellCommandFunction"->Print,*)
"CompileOptions" -> {" -Wall", "-Wextra",
"-Wno-unused-parameter", "-mmacosx-version-min=12.0",
"-std=c++20", "-Ofast", "-flto", "-Xpreprocessor -fopenmp",
"IncludeDirectories" -> {FileNameJoin[{dirOpenMP, "include"}]},
"LibraryDirectories" -> {FileNameJoin[{dirOpenMP, "lib"}]}
];
{
Sequence @@ ConstantArray[{Real, 1, "Constant"}, DomDim],
{Real, DomDim + Boole[AmbDim != 1], "Constant"},
{Real, 2, "Constant"},
If[AmbDim == 1, {Real, 1, "Shared"}, {Real, 2, "Shared"}],
Integer
},
"Void"
]];

Interpolate[grids_, vals_, x_] :=
Module[{threadCount, DomDim, AmbDim, y},
"ParallelThreadNumber" /. ("ParallelOptions" /.
SystemOptions["ParallelOptions"]);
DomDim = Length[grids];
y = ConstantArray[0., {Length[x]}];
AmbDim =
If[TensorRank[vals] == DomDim, 1,
Dimensions[vals][[DomDim + 1]]];
cInterpolate[DomDim, AmbDim][Sequence @@ grids, vals, x, y,
y
];


Here a brief test:

y = Interpolate[g, vals, x]; // RepeatedTiming // First
Max[Abs[result - y]]


0.00973586

6.66134*10^-16

This is more than nine times faster than cf2 and more than 800 times faster than the InterpolationFunction.

Edit 4

Here is a 3D-version of cf2.

D3toD3 = Compile[{{g1, _Real, 1}, {g2, _Real, 1}, {g3, _Real, 1}, {y, _Real, 3}, {x, _Real, 2}, {threadID, _Integer}, {threadCount, _Integer}},
Module[{jobCount, begin, end, q, r, i1, i2, i3, t1, t2, t3, s1, s2, s3},

jobCount = Length[x];

q = Quotient[jobCount, threadCount];
r = Mod[jobCount, threadCount];
begin = q (threadID - 1) + Quotient[r (threadID - 1), threadCount] + 1;

Table[

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g1];
a = 1;
b = n;
ga = CompileGetElement[g1, a];
gb = CompileGetElement[g1, b];
z = CompileGetElement[x, k, 1];
If[z <= ga,
i1 = a;
s1 = 1.;
t1 = 0.;
,
If[z >= gb,
i1 = b - 1;
s1 = 0.;
t1 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g1, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i1 = a; s1 = gb - z; t1 = z - ga;
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g2];
a = 1;
b = n;
ga = CompileGetElement[g2, a];
gb = CompileGetElement[g2, b];
z = CompileGetElement[x, k, 2];
If[z <= ga,
i2 = a;
s2 = 1.;
t2 = 0.;
,
If[z >= gb,
i2 = b - 1;
s2 = 0.;
t2 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g2, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i2 = a; s2 = gb - z; t2 = z - ga;
];
];
];

Block[{a, ga, b, gb, c, gc, z, n},
n = Length[g3];
a = 1;
b = n;
ga = CompileGetElement[g3, a];
gb = CompileGetElement[g3, b];
z = CompileGetElement[x, k, 3];
If[z <= ga,
i3 = a;
s3 = 1.;
t3 = 0.;
,
If[z >= gb,
i3 = b - 1;
s3 = 0.;
t3 = 1.;
,
While[b - a > 1,
c = Quotient[a + b, 2];
gc = CompileGetElement[g3, c];
If[z <= gc, b = c; gb = gc;, a = c; ga = gc];
];
i3 = a; s3 = gb - z; t3 = z - ga;
];
];
];

Plus[
s1 s2 s3 CompileGetElement[y, i1, i2, i3],
s1 s2 t3 CompileGetElement[y, i1, i2, 1 + i3],
s1 t2 s3 CompileGetElement[y, i1, 1 + i2, i3],
s1 t2 t3 CompileGetElement[y, i1, 1 + i2, 1 + i3],
t1 s2 s3 CompileGetElement[y, 1 + i1, i2, i3],
t1 s2 t3 CompileGetElement[y, 1 + i1, i2, 1 + i3],
t1 t2 s3 CompileGetElement[y, 1 + i1, 1 + i2, i3],
t1 t2 t3 CompileGetElement[y, 1 + i1, 1 + i2, 1 + i3]
]/((s1 + t1) (s2 + t2) (s3 + t3))
, {k, begin, end}]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"
]

• Thank you very much for the answer! First, a side question: is it possible to use the interpolation you have implemented for the numerical integration? If yes, then should we expect a speedup of the integration compared to the standard Interpolation? Apr 15, 2023 at 9:52
• Maybe. My functions are optimized for bulk evaluation, but NIntegrate prefers to evaluate functions one at a time for adaptive quadrature rules. Nonetheless, cf or a wrapper like F = {x1, x2, x3, x4} |-> cf[g1, g2, g3, g4, vals, {x1, x2, x3, x4}]; might work quite fine. Apr 15, 2023 at 10:04
• Thanks again. I have faced a problem: when trying to use cf2 for my example from the question, it got stuck (please see an update to the question). The same is the case for cf. What do you think, what may be the problem? Apr 15, 2023 at 10:23
• Your solution helps amazingly! Could you please tell me one thing? In your last implementation, what is g? Apr 16, 2023 at 19:47
• Ooops. g = {g1, g2, g3, g4}`. (I replaced it.) Apr 16, 2023 at 20:17