# How can I speed up my compiled function that is returned by another function?

Overview

I have a function that acts like the Interpolation on sparse $n$-dimensional data using a simple implementation of RBF interpolation method. I want my function to return a compiled function that will run fast. What I get works but it is much slower that I think it should be.

My code

Clear[RBFInterpolation]

Options[RBFInterpolation] = {
"DistanceFunction" -> (Norm[#1 - #2] &),
"RadialBasisFunction" -> (Sqrt[#1^2 + #2^2/4] &),
"RadialScale" -> Automatic, "Debug" -> False, "Compile" -> False};

RBFInterpolation[cptab_, opts : OptionsPattern[RBFInterpolation]] :=
Module[
{ro, xpts, fundata, Φ, disfun, λ, RBF, x},

xpts = #[[1]] & /@ cptab;
fundata = #[[2]] & /@ cptab;
disfun = OptionValue["DistanceFunction"];

Φ =
Table[disfun[xpts[[i]], xpts[[j]]], {i, 1, Length[xpts]},{j,1,Length[xpts]}];

Which[
ro = Median[
Flatten[Table[
Drop[Φ[[i]], {i}], {i, 1,
Length[Φ]}]]],

True,
OptionValue["RadialScale"], " So I'm going to make it up"];
ro = Median[
Flatten[Table[
Drop[Φ[[i]], {i}], {i, 1, Length[Φ]}]]]
];

If[OptionValue["Debug"], Print["ro=", ro]];
If[OptionValue["Debug"],
Print["Distance function on first two points"];
Print["point 1 ->", xpts[[1]]];
Print["point 2 ->", xpts[[2]]];
Print["Distance ->", disfun[xpts[[1]], xpts[[2]]]];
Print["Radial Basis Function on Distance ->",
RBF[disfun[xpts[[1]], xpts[[2]]], ro]]
];

Φ = Map[RBF[#, ro] &, Φ, {2}];

If[OptionValue["Debug"],
Print["Element of Φ[[1,1]]=", Φ[[1,1]]]];

λ = Inverse[Φ].fundata;

If[OptionValue["Debug"],
Print["First element of λ[[1]]=", λ[[i]]]];

If[OptionValue["Compile"],
Return[
With[{xi = x, λi = λ, xptsi = xpts, roi = ro},
Compile[{{xi, _Real, 1}},
Sum[λi[[i]] RBF[disfun[xi, xptsi[[i]]], roi], {i, 1,
Length[λ]}]]]],

Return[
Function[x,
Sum[λ[[i]] RBF[disfun[x, xpts[[i]]], ro], {i, 1,
Length[λ]}]]]
]
];


Most of this function is not of interest to my question. I think the key is where I Return[] the compiled function.

Return[
With[{xi = x, λi = λ, xptsi = xpts, roi = ro},
Compile[{{xi, _Real, 1}},
Sum[λi[[i]] RBF[disfun[xi, xptsi[[i]]], roi], {i, 1,
Length[λ]}]]]]


Testing the Function

The following code can be used to run and test the timing of the returned function.

First make a "Truth" function to sample then interpolate

Clear[truth]
truth[x_] := Product[Sin[x[[i]]], {i, 1, Length[x]}];


Make up some data

n = 100;
d = 5;
cpts = RandomReal[{-π/2, π/2}, {n, d}];
cptab = {#, truth[#]} & /@ cpts;
xpts = #[[1]] & /@ cptab;
fundata = #[[2]] & /@ cptab;


Test the speed of the returned functions

Print["Normal Function:"];
Timing[funFun = RBFInterpolation[cptab, "Compile" -> False];]
Timing[funFun[#] & /@ xpts;]

Print["Compile Function:"];
Timing[funFunc = RBFInterpolation[cptab, "Compile" -> True];]
Timing[funFunc[#] & /@ xpts;]
i = 1;
Print["Normal function: ", funFun[xpts[[i]]]];
Print["Complie function: ", funFunc[xpts[[i]]]];
Print["The real right answer: ", fundata[[i]]];


I get results like this:

Normal Function:
{0.080987,Null}
{0.123981,Null}

Compile Function:
{0.092986,Null}
{0.156977,Null}

Normal function: -0.0182901
Complie function: -0.0182901


So as you can see it works but it is not faster. How do I make this faster?

Simpler test that is faster!?

The code:

n = 10;
a = RandomReal[{-1, 1}, n];
f = Table[2 π i, {i, 1, n}];
ϕ = RandomReal[{0, 2 π}, n];
Clear[Nfun]
Nfun[t_] := Sum[a[[i]] Cos[f[[i]] t + ϕ[[i]]], {i, 1, n}];
Nfunc = Compile[{{t, _Real}},
Evaluate[Sum[a[[i]] Cos[f[[i]] t + ϕ[[i]]], {i, 1, n}]]];

Clear[makeNfunc]
makeNfunc[a_, f_, ϕ_] := Module[{n},
n = Length[a];
Return[
Compile[{{t, _Real}},
Evaluate[Sum[a[[i]] Cos[f[[i]] t + ϕ[[i]]], {i, 1, n}]]]]
];
NfuncR = makeNfunc[a, f, ϕ];


Run the code:

npts = 10000;
data = RandomReal[{0, 10}, npts];
Timing[Nfun[#] & /@ data;]
Timing[Nfunc[#] & /@ data;]
Timing[NfuncR[#] & /@ data;]


The output:

{0.585911, Null}
{0.012998, Null}
{0.012998, Null}


So in this simple case the compiled code is about 45 times faster for both the function compiled inline Nfunc and the function that was returned by the makeNfunc, NfuncR

So the question is what is the problem with my original function above?

-
I just want to add a link to an excellent package that provides a lot of $n$-space support Obtuse Package (I just need some more speed!) –  c186282 Aug 29 '13 at 22:45
Compile has been improved since version 7, which I use, but when I run the line Timing[funFunc[#] & /@ xpts;] I get CompiledFunction::cfte: Compiled expression 0. should be a rank 1 tensor of machine-size real numbers. >> CompiledFunction::cfex: Could not complete external evaluation at instruction 20; proceeding with uncompiled evaluation. >> Do you see similar errors? –  Mr.Wizard Aug 30 '13 at 0:59
In the code as shown in the post I do not get any errors. Early on I was playing with the rank which gave me similar errors. I then found that 1 worked which makes since because xi is a vector. I get no indication that what is returned is not compiled, like the error text proceeding with uncompiled evaluation. implies. If I just evaluate funFunc I get CompiledFunction[...stuff..] as expected. Everything works it is just no faster. –  c186282 Aug 30 '13 at 1:13
I added an answer. If I am correct using CompilationOptions -> {"InlineExternalDefinitions" -> True} will fix your problem. –  Mr.Wizard Aug 30 '13 at 1:18

I had to modify your code to get it to work without error in version 7. Once I did it that appears to be working correctly and faster than the non-compiled code.

I needed to inject the values of RBF and disfun into the Compile using With:

With[{iRBF = RBF, idisfun = disfun},
If[OptionValue["Compile"],
Return[With[{xi = x, λi = λ, xptsi = xpts, roi = ro},
Compile[{{xi, _Real, 1}},
Sum[λi[[i]] iRBF[idisfun[xi, xptsi[[i]]], roi], {i, 1,
Length[λ]}]]]],
Return[Function[x, Sum[λ[[i]] iRBF[idisfun[x, xpts[[i]]], ro], {i, 1, Length[λ]}]]]]
]


I believe that in later versions this can be done with:

CompilationOptions -> {"InlineExternalDefinitions" -> True}


n = 300;
d = 5;
cpts = RandomReal[{-\[Pi]/2, \[Pi]/2}, {n, d}];
cptab = {#, truth[#]} & /@ cpts;
xpts = #[[1]] & /@ cptab;
fundata = #[[2]] & /@ cptab;

Print["Normal Function:"];
Timing[funFun = RBFInterpolation[cptab, "Compile" -> False];]
Timing[funFun /@ xpts;]

Print["Compile Function:"];
Timing[funFunc = RBFInterpolation[cptab, "Compile" -> True];]
Timing[funFunc /@ xpts;]
i = 1;
Print["Normal function: ", funFun[xpts[[i]]]];
Print["Complie function: ", funFunc[xpts[[i]]]];
Print["The real right answer: ", fundata[[i]]];


Normal Function:

{0.514, Null}

{0.546, Null}

Compile Function:

{0.515, Null}

{0.094, Null}

Normal function: 0.000268092

Complie function: 0.000268092

Thank you I had tried many combos with With but I missed that one thank you. –  c186282 Aug 30 '13 at 1:20
Thank you for the advice I will toggle the answer for 24 hours. Also thank you for the edit. How did you turn \Pi` in the real symbol? –  c186282 Aug 30 '13 at 1:25