# 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. – J. M.'s discontentment Jul 11 '15 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 Oct 19 '15 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.. – Alexey Popkov Oct 19 '15 at 7:53
• In fact, @AlexeyPopkov the MachinePrecsion is enough in my work. – xyz Oct 19 '15 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. – J. M.'s discontentment Oct 20 '15 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[], 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 Oct 19 '15 at 7:40
• Agreed. For some reason, it seemed faster in version 8 than in version 10… – J. M.'s discontentment Oct 19 '15 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 Oct 20 '15 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. – J. M.'s discontentment Oct 20 '15 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. – J. M.'s discontentment Oct 20 '15 at 3:25