# Multidimensional arbitrary precision spline interpolation on the grid

This question is a generalization of the previous one for multiple dimensions. In the answer to that question an implementation for the clamped spline interpolation for 1D case and arbitrary degree of spline is given. How can it be extended for multidimensional case?

• I've figured out how to handle two and three variables, but I'm having difficulty extending it to an arbitrary number of dimensions. If it's okay with you, I can post the incomplete solution. Commented Jul 11, 2015 at 20:59
• What does the arbitrary precision mean? Could you give the difinition or a demo to show arbitrary precision interpolation? Or I would like to know which case did you need to arbitrary precision interpolation?
– xyz
Commented Oct 19, 2015 at 3:25
• @Shutao You could easily find the answer in the Documentation: "Arbitrary‐Precision Numbers," "Arbitrary‐Precision Calculations" or in Wikipedia: "Arbitrary-precision arithmetic." An example of arbitrary precision interpolation with short discussion is given in previous answer by J.M.. Commented Oct 19, 2015 at 7:53
• In fact, @AlexeyPopkov the MachinePrecsion is enough in my work.
– xyz
Commented Oct 19, 2015 at 9:31
• Here's one possible reason to use arbitrary precision: if you're using Gröbner basis methods to find intersections of B-spline surfaces (they are, after all, piecewise polynomials!), it is useful to have a representation that can work with exact or arbitrary-precision arithmetic, since these methods can sometimes be sensitive. Commented Oct 20, 2015 at 2:29

I still haven't figured out how to write a routine for arbitrary dimension, but I'm posting my (incomplete!) solution in case people might have ideas on extending what I have.

## Bivariate interpolant

Here is a random bivariate polynomial, which we'll use for generating test data:

f[x_, y_] := -2 + 4 x^2 + 4 x^3 - 3 y - 5 x^2 y + 5 y^2 + 5 x y^2 + y^3


Here's some test data from f[x, y], sampled at non-equispaced points:

da = Flatten[Table[{{x, y}, f[x, y]},
{x, {-2, -4/3, -1/5, 1/9, 3/4, 1}},
{y, {0, 1/6, 3/8, 9/5, 2}}], 1];


Here's the reference interpolating function:

{p, q} = {2, 3}; (* spline degrees in the two variables *)
ipf = Interpolation[da, InterpolationOrder -> {p, q}, Method -> "Spline"];


Some preliminary processing to separate out independent and dependent variables:

{pts, vals} = Transpose[SplitBy[SortBy[da, First], #[[1, 1]] &], {3, 2, 1}];


Make the knot sequence for each independent variable:

makeKnots[list_?VectorQ, deg_Integer?NonNegative] :=
With[{n = Length[list]},
Join[ConstantArray[list[[1]], deg + 1],
If[deg + 2 <= n, MovingAverage[ArrayPad[list, -1], deg], {}],
ConstantArray[list[[-1]], deg + 1]]]

{u, v} = {pts[[1, All, 1]], pts[[All, 1, 2]]};
{uk, vk} = MapThread[makeKnots, {{u, v}, {p, q}}];


Build the control points:

{m, n} = {Length[u], Length[v]};
usol = LinearSolve[Outer[BSplineBasis[{p, uk}, #2, #1] &,
u, Range[0, m - 1], 1]];
vsol = LinearSolve[Outer[BSplineBasis[{q, vk}, #2, #1] &,
v, Range[0, n - 1], 1]];

cpts = vsol /@ Transpose[usol /@ vals];


Finally, the bivariate interpolating spline:

spf[x_, y_] = Fold[Dot, cpts, {Table[BSplineBasis[{q, vk}, k - 1, y], {k, n}],
Table[BSplineBasis[{p, uk}, k - 1, x], {k, m}]}];


Plot the two interpolants and the original function:

MapThread[Plot3D[#1[x, y], {x, -2, 1}, {y, 0, 2}, PlotLabel -> #2] &,
{{f, ipf, spf}, {"True", "InterpolatingFunction", "B-spline"}}]
// GraphicsRow


Evaluate ipf[] and spf[] at the same argument:

{ipf[-1, 1], spf[-1, 1]}
{-8.59478, -552429212/64275003}


Note that only the second function gave exact output.

The difference between ipf[] and spf[], showing good agreement:

Plot3D[ipf[x, y] - spf[x, y], {x, -2, 1}, {y, 0, 2}, PlotRange -> All]


If you want further confirmation, you can try recovering the control points and knots from ipf[], using a procedure similar to the one in my previous answer, and compare them with the control points and knots I generated here.

## Trivariate interpolant

Hopefully, you can see the similarities and differences between the previous example and this one:

(* random polynomial *)
mp[x_, y_, z_] := 3 - 7 x^2 + 2 x^3 - 2 y + 8 y^2 + 5 x y^2 + 8 y^3 + 4 x y^3 -
2 y^4 + 6 x z + 2 x^2 z - 6 x^3 z + 2 y z + 4 x y z -
2 x^2 y z - 3 y^2 z - 8 x y^2 z + 7 y^3 z - 5 z^2 + x z^2 +
2 x^2 z^2 + 6 y z^2 - 4 y^2 z^2 + 9 z^3 - 4 x z^3 -
3 y z^3 - 9 z^4
(* random data *)
da = Flatten[Table[{{x, y, z}, mp[x, y, z]}, {x, {-2, -3/2, -9/7, 7/9, 2, 3}},
{y, {1, 5/3, 5/2, 9/2, 5}}, {z, {-1, 5/7, 11/10, 2}}], 2];

{p, q, r} = {4, 3, 2}; (* B-spline degree for each variable *)
ipf = Interpolation[da, InterpolationOrder -> {p, q, r}, Method -> "Spline"];

{pts, vals} = Transpose[GatherBy[da, {#[[1, 1]] &, #[[1, 2]] &}], {4, 3, 2, 1}];

(* make knots *)
{u, v, w} = {pts[[1, 1, All, 1]], pts[[1, All, 1, 2]], pts[[All, 1, 1, 3]]};
{uk, vk, wk} = MapThread[makeKnots, {{u, v, w}, {p, q, r}}];

(* make control points *)
{l, m, n} = Length /@ {u, v, w};
usol = LinearSolve[Outer[BSplineBasis[{p, uk}, #2, #1] &, u, Range[0, l - 1], 1]];
vsol = LinearSolve[Outer[BSplineBasis[{q, vk}, #2, #1] &, v, Range[0, m - 1], 1]];
wsol = LinearSolve[Outer[BSplineBasis[{r, wk}, #2, #1] &, w, Range[0, n - 1], 1]];
cpts = Map[wsol, Transpose[Map[vsol,
Transpose[Map[usol, vals, {2}], {2, 3, 1}], {2}], {1, 3, 2}], {2}];

(* B-spline interpolant *)
spf[x_, y_, z_] = Fold[Dot, cpts, {Table[BSplineBasis[{r, wk}, k - 1, z], {k, n}],
Table[BSplineBasis[{q, vk}, k - 1, y], {k, m}],
Table[BSplineBasis[{p, uk}, k - 1, x], {k, l}]}];


Tests:

{ipf[2, 3, 1], spf[2, 3, 1]}
{383.531, 10447533501473/27240371550}

ContourPlot3D[#[x, y, z], {x, -2, 3}, {y, 1, 5}, {z, -1, 2},
BoxRatios -> Automatic, MaxRecursion -> 0] & /@
{ipf, spf} // GraphicsRow


• The program in Trivariate interpolant was very time-consuming.
– xyz
Commented Oct 19, 2015 at 7:40
• Agreed. For some reason, it seemed faster in version 8 than in version 10… Commented Oct 19, 2015 at 7:51
• Although the precision of usol is Infinity, when you use ContourPlot3D, the precison will change to the MachinePrecison. In addition, you use the property of B-spline basis. Namely, In generally, when $u_0 \in \left[u _i ,u _{i+1} \right)$, the B-spline function basis $N_{i-p,p},\cdots ,N_{i,p}$ are not equal to 0.
– xyz
Commented Oct 20, 2015 at 2:07
• That's correct. The only reason for showing the plot is to visually demonstrate that the explicitly-constructed B-spline is equivalent to the output of Interpolation[]. If you want to force the use of higher precision, the plotting functions all accept a WorkingPrecision setting. The explicitly constructed B-spline will keep up; the InterpolatingFunction[] will still be stuck at MachinePrecision. Commented Oct 20, 2015 at 2:19
• @Shutao, to repeat what I said: "look at the low order cases, try to determine a pattern"; sometimes, one just starts by guessing. :) As a simpler example: if you keep differentiating $x^k$ $n$ times for various small values of $k$ and $n$, you'll eventually notice the pattern $\frac{k!}{(k-n)!}x^{k-n}$, which you then proceed to prove inductively. A lot of work in combinatorics actually proceeds in this way. Commented Oct 20, 2015 at 3:25