# Combining two interpolating functions via Piecewise or If is very slow

For a given set of parameters $$(i,l,m)$$ I have a function which has been evaluated on a grid $$x \in (-2,2)$$. Since this function is not smooth (only continuous) at the "interface" $$x=0$$, I am interpolating separately the values defined for $$x<0$$ and $$x>0$$. I then combine them into a single function using Piecewise or If.

xmin = -2; xInterface = 0; xmax = 2;
xgrid =  Range[xmin, xmax, 1/100];

Clear[myPolyLeft,myPolyRight,leftInterpol, rightInterpol, fullInterpol]
Monitor[Do[
(* set up the table left and right of the interface *)
interfaceValue = RandomReal[{1, 10}];
myPolyLeft[x_] = RandomReal[{1, 3}] x^2 + RandomReal[{1, 10}] x + interfaceValue;
myPolyRight[x_] = RandomReal[{1, 3}] x^2 + RandomReal[{1, 10}] x + interfaceValue;
leftTable = Table[myPolyLeft[xmin + (xInterface - xmin) i/200], {i, 0, 200}];
rightTable = Table[myPolyRight[xInterface + (xmax - xInterface) i/200], {i, 0, 200}];

(* interpolate and stitch the solutions together *)
leftInterpol[i][l, m] = Interpolation[Transpose[{xgrid[[;; 201]], leftTable}]];
rightInterpol[i][l, m] = Interpolation[Transpose[{xgrid[[201 ;;]], rightTable}]];
fullInterpol[i][l, m][x_] =  If[x < xInterface, Evaluate[leftInterpol[i][l, m][x]], Evaluate[rightInterpol[i][l, m][x]]]
, {i, 1, 10}, {l, 2, 40}, {m, -l, l}]
, {i, l, m}]


What I am surpised to see is that evaluating fullInterpol on a grid takes about 1000x longer than evaluating the indvidual "left" and "right" interpolating functions:

(* Do using i:1-10,l:2-40 *)
Table[rightInterpol[10][3, 2][xInterface + i (xmax - xInterface)/100], {i, 0, 100}]; // AbsoluteTiming
Table[leftInterpol[10][3, 2][xmin + i (xInterface - xmin)/100], {i, 0, 100}]; // AbsoluteTiming
Table[fullInterpol[10][3, 2][xmin + i (xmax - xmin)/100], {i, 0, 100}]; // AbsoluteTiming
(* Timing Output:0.000275,0.000258,0.312824 *)


Even more perplexing, as @DanielHuber pointed out, if I then set

fun[10][3, 2][x_] = fullInterpol[10][3, 2][x];


and do a timing test on fun:

Table[fun[10][3, 2][xmin + i (xmax - xmin)/100], {i, 0,  100}]; // AbsoluteTiming
(* Timing Output: 0.000327 *)


Furthermore, the time to evaluate fullInterpol seems to be affected by the range of parameter $$(i,l,m)$$. For example, choosing to only apply the Do loop for i=10, the timings are then:

(* Do using i:10-10, l:2-40 *)
Table[rightInterpol[10][3, 2][xInterface + i (xmax - xInterface)/100], {i, 0, 100}]; // AbsoluteTiming
Table[leftInterpol[10][3, 2][xmin + i (xInterface - xmin)/100], {i, 0, 100}]; // AbsoluteTiming
Table[fullInterpol[10][3, 2][xmin + i (xmax - xmin)/100], {i, 0, 100}]; // AbsoluteTiming
(*Timing Output: 0.000277, 0.00026, 0.000985 *)


The evaluation of the two interpolating functions are unaffected as expected, but fullInterpol is now much faster. Something similar can be seen by reducing the range of $$l$$ instead.

What is going on?

Perhaps avoid SubValues and even DownValues. The problem is set up for simple table look up, and with this approach, it is just as fast as evaluating the individual interpolating functions:

Clear[fullInterpol2]
Monitor[interpData = Table[
{leftInterpol[i][l, m], rightInterpol[i][l, m]}, {i, 1, 10}, {l,
2, 40}, {m, -l, l}], {i, l, m}];
fullInterpol2[i_, l_, m_, x_] :=
interpData[[i, l - 1, m + l + 1, 1 + UnitStep[x]]][x];

Table[rightInterpol[10][3, 2][
xInterface + i (xmax - xInterface)/100], {i, 0,
100}]; // AbsoluteTiming
Table[leftInterpol[10][3, 2][xmin + i (xInterface - xmin)/100], {i, 0,
100}]; // AbsoluteTiming
foo1 = Table[
fullInterpol[10][3, 2][xmin + i (xmax - xmin)/100], {i, 0,
100}]; // AbsoluteTiming
foo2 = Table[
fullInterpol2[10, 3, 2, xmin + i (xmax - xmin)/100], {i, 0,
100}]; // AbsoluteTiming
(*
{0.000353, Null}
{0.000332, Null}
{0.135156, Null}
{0.000336, Null}
*)

(* Check *)
foo1 == foo2
(*  True  *)


You have a bug in the definition of fullInterpol. "If" has the attribute "HoldRest". This makes the subsequent evaluation symbolically. To se this, consider:

l = 3; m = 2;
fullInterpol[i][l, m][x_] =
If[x < xInterface, leftInterpol[i][l, m][x],
rightInterpol[i][l, m][x]]


this returns:

If[x < 0, leftInterpol[i][l, m][x], rightInterpol[i][l, m][x]]


l and m have no value.

To fix this, you must force evaluation like:

l = 3; m = 2;
fullInterpol[i][l, m][x_] =
If[x < xInterface, Evaluate[leftInterpol[i][l, m][x]],
Evaluate[rightInterpol[i][l, m][x]]]


This now returns:

If[x < 0, leftInterpol[i][3, 2][x], rightInterpol[i][3, 2][x]]

• Hi Daniel. Thank you. I edited the question to fix this bug. I still obtain the same kind of timing tests though Feb 17, 2023 at 14:04
• If you set: fun[x_] = fullInterpol[10][3, 2][x]; and do the timing on fun, then it is fast. I do not understand this. Maybe if you can simplify it it is worth a report to [email protected]. Feb 17, 2023 at 14:47
• I updated the question. This is quite perpexing... Feb 17, 2023 at 16:48