10

4

**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"];
RBF = OptionValue["RadialBasisFunction"];
Φ =
Table[disfun[xpts[[i]], xpts[[j]]], {i, 1, Length[xpts]},{j,1,Length[xpts]}];
Which[
OptionValue["RadialScale"] == Automatic,
ro = Median[
Flatten[Table[
Drop[Φ[[i]], {i}], {i, 1,
Length[Φ]}]]],
NumberQ[OptionValue["RadialScale"]],
ro = OptionValue["RadialScale"],
True,
Print["I cannot understand \"RadialScale\"->",
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
The real right answer: -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?

The routine given here for thin plate splines might be useful to you; you will only need to change the underlying RBF.

– J. M.'s ennui – 2015-07-25T14:36:42.017I 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 – 2013-08-29T22:45:40.630`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 – 2013-08-30T00:59:36.567

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 – 2013-08-30T01:13:49.357I added an answer. If I am correct using

`CompilationOptions -> {"InlineExternalDefinitions" -> True}`

will fix your problem. – Mr.Wizard – 2013-08-30T01:18:33.053