Multidimensional arbitrary precision spline interpolation on the grid

Question

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?

Answer 1

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