# How to reduce the InterpolatingFunction building overhead?

I want a linear interpolation from the following example list:

list = {{0.0005023, 22.24}, {0.01457, 21.47}, {0.04922, 19.79},
{0.07484, 18.7}, {0.104, 17.55}, {0.1331, 16.52}, {0.1632, 15.49},
{0.1888, 14.52}, {0.2215, 13.31}, {0.2506, 12.16}, {0.3024, 10.01},
{0.3435, 8.304}, {0.3943, 6.036}, {0.4098, 5.329}, {0.4726, 2.384}};


The easiest way is to use:

Interpolation[list, InterpolationOrder -> 1]


but my list will be changing a lot, and the InterpolatingFunction takes a lot of time to build:

Timing[
Table[Interpolation[list, InterpolationOrder -> 1][q], {q,
0.0006, 0.4, 0.00001}];]


is 10× slower than:

test=Interpolation[list, InterpolationOrder -> 1];
Timing[Table[test[q], {q, 0.0006, 0.4, 0.00001}];]


How can I remove the overhead?

EDIT (following JxB comment)

This compiled version is 5 times faster than the original version, but I don't think Partition is compiling (it appears between all the Lists when I use FullForm); and there's also a CopyTensor that doesn't look good:

Compile[{{list, _Real, 2}, {value, _Real, 0}},
Module[{temp},
temp = Select[
Partition[list, 2, 1], #[[1, 1]] <= value && #[[2, 1]] > value &][[1]
];
temp[[1, 2]] +
(value - temp[[1, 1]])/(temp[[2, 1]] - temp[[1, 1]])*(temp[[2, 2]] - temp[[1, 2]])
]
]


Any suggestions? (I don't want to compile to C.)

• Perhaps Timing[Table[ Evaluate@Interpolation[list, InterpolationOrder -> 1][q], {q, 0.0006, 0.4, 0.00001}];] ? May 2, 2012 at 12:46
• Since you are using linear interpolation, it might be straightforward to build your own compiled version of an interpolating function.
– JxB
May 2, 2012 at 12:53
• Do your x-value remain the same from list to list? May 2, 2012 at 14:50
• @rcollyer my x-value also changes May 2, 2012 at 16:21
• Opps... from the answers and comments I see that apparently I didn't explain myself correctly. The Table was there just to raise the timing to readable measurements (I should have used Do...). The calls to "test" will be made almost for one q at a time, and in-between list may change. May 2, 2012 at 16:30

You can use binary search with Compile. I failed inlining (Compile was complaining endlessly about types mismatch), so I included a binary search directly into Compile-d function. The code for binary search itself corresponds to the bsearchMin function from this answer.

Clear[linterp];
linterp =
Compile[{{lst, _Real, 2}, {pt, _Real}},
Module[{pos  = -1 , x = lst[[All, 1]], y = lst[[All, 2]], n0 = 1,
n1 = Length[lst], m = 0},
While[n0 <= n1, m = Floor[(n0 + n1)/2];
If[x[[m]] == pt,
While[x[[m]] == pt  && m < Length[lst], m++];
pos = If[m == Length[lst], m, m - 1];
Break[];
];
If[x[[m]] < pt, n0 = m + 1, n1 = m - 1]
];
If[pos == -1, pos = If[x[[m]] < pt, m, m - 1]];
Which[
pos == 0,
y[[1]],
pos == Length[x],
y[[-1]],
True,
y[[pos]] + (y[[pos + 1]] - y[[pos]])/(x[[pos + 1]] -
x[[pos]])*(pt - x[[pos]])
]],
CompilationTarget -> "C"];


This is about 20 times faster, on my benchamrks:

AbsoluteTiming[
Table[Interpolation[list,InterpolationOrder->1][q],{q,0.0006,0.4,0.00001}];
]


{1.453,Null}

AbsoluteTiming[
Table[linterp[list,q],{q,0.0006,0.4,0.00001}];
]


{0.063,Null}

• I just realized the OP specifically said "I don't want to compile to C." What timings do you get on v8 using WVM? May 2, 2012 at 19:18
• @Mr.Wizard It's about twice slower. The difference is not so dramatic, because the complexity itself is only logarithmic. May 2, 2012 at 19:48
• @P.Fonseca It is, but it is a hassle. You will need to do this for all platforms (since shared libs are platform-dependent), and you will have to dispatch to a right one depending on the platform. Nothing too complicated though, and sounds like a very good project (since this is a general problem). I think, a general module for serialization of compiled functions would be very useful. The problem would be, that to use it on a single machine, one would need all cross-compilers to other platforms installed. Macs seem best equipped for that, but it would still be some work to set that up. May 2, 2012 at 21:03
• @P.Fonseca if you are precompiling, then sure, I can understand only using WVM. However, if you load CCompilerDriver  the function DefaultCCompiler[] will generate the message CreateLibrary::nocomp and return $Failed if no default compiler was set. Personally, I'd Quiet the message and test for $Failed directly, as there is less chance of the message escaping. Using that info, you can choose your compilation target. May 3, 2012 at 1:23
• @rcollyer "that may be counter productive with the OPs purpose" - this is exactly the reason why I did not go there. Inlining a list would make sense when the list doesn't change, which is not the case OP is interested in. May 3, 2012 at 9:47

Combinatorica functions are often not well optimized, so there may very well be a faster binary search algorithm. If that can be found, this might be effective:

Needs["Combinatorica"]

f[{{a_, b_}, {c_, d_}}][x_] := b + (d - b)/(c - a) (x - a)

list[[Floor@{#, # + 1}]] & @ BinarySearch[list[[All, 1]], 0.33]

f[%][0.33]

8.86436


Check:

Interpolation[list, InterpolationOrder -> 1][0.33]

8.86436

• I was fighting with inlining in Compile (and lost), and was beaten to it by you, as a result (so lost again). +1. May 2, 2012 at 18:57
• @Leonid if I had taken the time to look up the faster binary search (which I believe you posted before) you'd have been first. +1 to you as well. :-) May 2, 2012 at 19:11
• It seems you don't need to load Combinatorica anymore; you can use GeometricFunctionsBinarySearch[] instead. May 10, 2013 at 14:03

Here is a linear interpolation routine that uses binary search with a few refinements (in particular, the binary search is skipped in the case of equispaced abscissas), as well as a stabilized version of the linear interpolation formula:

lerp = Compile[{{dat, _Real, 2}, {x, _Real}},
Module[{n = Length[dat], k = 1, l, m, r, xa, ya},
{xa, ya} = Transpose[dat];
l = Min[Max[2, 1 + Quotient[x - First[xa],
(Last[xa] - First[xa])/(n - 1)]], n - 1];

If[xa[[l]] <= x,
r = l + 1;
While[r < n && xa[[r]] <= x,
l = r; k *= 2; r = Min[l + k, n]],
{l, r} = {l - 1, l};
While[1 < l && x < xa[[l]],
r = l; k *= 2; l = Max[1, r - k]]];

While[r - l > 1,
m = Quotient[l + r, 2];
If[x < xa[[m]], r = m, l = m]];

({xa[[r]] - x, x - xa[[l]]}/(xa[[r]] - xa[[l]])).ya[[{l, r}]]],
RuntimeOptions -> "Speed"]


Even without the compilation to C, the method is quite fast on my box:

AbsoluteTiming[Table[Interpolation[data, InterpolationOrder -> 1][q],
{q, 0.0006, 0.4, 0.00001}];][[1]]
15.206078

AbsoluteTiming[Table[lerp[data, q], {q, 0.0006, 0.4, 0.00001}];][[1]]
0.693506


Something with memory ?

myTest[alist_] :=  myTest[alist] = Interpolation[alist, InterpolationOrder -> 1]

Timing[Table[myTest[list][q], {q, 0.0006, 0.4, 0.00001}];]

(* {0.187,Null} *)

test=Interpolation[list,InterpolationOrder->1];
Timing[Table[test[q],{q,0.0006,0.4,0.00001}];]

(* {0.172,Null} *)

• This might work. It builds the InterpolatingFunction` just once for each new list. I don't think I will have more than a couple thousands lists during a session, and so, I believe still manageable. Do you see a way to purge the memorization (memory or quantity) if it passes over a certain value (obvious, without too much overhead…)? May 2, 2012 at 16:37
• @P.Fonseca Do you mean remove the previous definitions if they exceed a certain memory usage ? If so, sorry I don't. May 2, 2012 at 18:38
• Yes. If the total myTest definitions exceed a certain quantity or memory usage May 2, 2012 at 18:41