# Arbitrary precision spline interpolation

The current implementation of Interpolation does not allow arbitrary precision spline interpolation. Yu-Sung Chang says here that "it is not hard to implement it manually for arbitrary precision using BSplineBasis."

I am new to splines and not familiar with BSplineBasis. What is the best way to implement arbitrary precision spline interpolation in Mathematica?

I am interested primarily in an analog of the Method -> "Spline" of Interpolation which I investigated in this answer (BTW, what is the name of such a spline and of this kind of parametrization?). Also, J. M. gives in this answer an interesting implementation of a spline with centripetal parametrization, but it is not clear how to use it for interpolation.

P.S. I need an implementation where the degree of a spline is arbitrary. At the moment I am more interested in even-degree splines.

• Yu-Sung's pointer is wolfram.com/xid/0c0rpcn4bku6m-dj3q8d, where there are two examples of such a thing (I haven't tested them). What are those example lacking? – Dr. belisarius Jan 10 '14 at 12:45
• That is an example of cubic B-spline interpolation. I need an implementation where the degree of spline is arbitrary. – Alexey Popkov Jan 10 '14 at 12:51
• Please see my question about B-spline Basis Function – xyz May 28 '15 at 8:12

The OP linked to an answer of mine for interpolating over general point sets; for constructing a single interpolating function, a slight modification of my procedure is needed. (In particular, you don't need centripetal or chord-length parametrization in this case.)

data = {{0, 0}, {1/10, 3/10}, {1/2, 3/5}, {1, -1/5}, {2, 3}, {3, -6/5}};


To review the problems with using the built-in Interpolation[], let's try building an interpolating quartic spline from data:

m = 4; (* spline order *)
sp1 = Interpolation[data, InterpolationOrder -> m, Method -> "Spline"];


Now, evaluate:

sp1[3/2]
0.26753591659167364


Herein lies the problem: even though the contents of data and the argument fed to sp1 are exact, the output is a machine precision number. Yeesh!

So, we go back to basics and build the interpolating spline ourselves from BSplineBasis[]. First, we generate the knots. The clamped condition is the most appropriate for interpolation, so we generate the appropriate knot sequence like so:

n = Length[data];
{xa, ya} = Transpose[data]
knots = Join[ConstantArray[xa[[1]], m + 1],
If[m + 2 <= n, MovingAverage[ArrayPad[xa, -1], m], {}],
ConstantArray[xa[[-1]], m + 1]];


To generate the control points for the B-spline, we construct the relevant linear system and solve it with LinearSolve[]:

cp = LinearSolve[Outer[BSplineBasis[{m, knots}, #2, #1] &,
xa, Range[0, Length[data] - 1]], ya];


Note that all the entries of cp are exact.

At this juncture, you might think that we can use BSplineFunction[] on the control points and the knots, but alas, the function also only evaluates at machine precision. We are left with no choice but to use BSplineBasis[] once more:

sp2[x_] = cp.Table[BSplineBasis[{m, knots}, j - 1, x], {j, n}];


Try it out:

sp2[3/2]
25251843/94386740

N[%]
0.2675359165916738


As a verification that we were able to reproduce the control points used internally by Interpolation[], let's do a comparison:

sp1bs = Cases[sp1, _BSplineFunction, ∞][[1]];
sp1bs["ControlPoints"]
{0., 0.777830630128766, 1.3470257035045392, -6.670028631670085, 13.100332135107115, -1.2}

N[cp]
{0., 0.777830630128766, 1.347025703504539, -6.6700286316700845, 13.100332135107113, -1.2}


Certainly, I can't end this post without at least one picture, so:

Plot[{sp1[x], sp2[x]}, {x, 0, 3}]


Plot[sp1[x] - sp2[x], {x, 0, 3}, PlotRange -> All]


I leave the task of bundling everything in a routine as an exercise. If you want to learn more, Piegl and Tiller's The NURBS Book is the canonical reference.

• It is good to have you back. – Mr.Wizard May 24 '15 at 10:31
• I was able to borrow a friend's computer with a Mathematica installation, so I went for it and answered some questions. It was fun while it lasted… :) – J. M.'s technical difficulties May 24 '15 at 10:45
• Did you try contacting the Stack Exchange team about the matter? I understand if you do not wish to pursue it but I think they might help. – Mr.Wizard May 24 '15 at 10:48
• Maybe I'll ask them sometime… thanks for the concern. :) – J. M.'s technical difficulties May 24 '15 at 10:51
• @J.M. May I ask you also about a generalization of the spline interpolation (at least) for the 2D case? In the linked answer it is said that "Multi-dimension is simply tensor product version." I do not see this as a simple thing but it would be very helpful to have such a generalization. May be a separate thread would be the best place for this. – Alexey Popkov May 26 '15 at 14:01

Here is my implementation of quadric spline arbitrary precision interpolation which is defined exactly as in Interpolation with options Method->"Spline", InterpolationOrder->2.

## Theoretical background

Quadric spline interpolation for n datapoints is defined as a Piecewise function which consists of n-2 parabolas which are spliced together in the middles of successive datapoints (with exceptions for first and last intervals) by two conditions: successive parabolas at these middle points must be equal and their first derivatives must be equal too. So we have 2(n-3) such conditions.

Let us designate jth datapoint as {x[j],y[j]} and define jth parabola as

s[j_, xx_] = y[j + 1] + b[j] (xx - x[j + 1]) + c[j] (xx - x[j + 1])^2;


This definition automatically makes jth parabola equal to y[j + 1] at x[j + 1]. We have to add also 2 boundary conditions:

s[1, x[1]] == y[1]
s[n - 2, x[n]] == y[n]


And 2(n-3) splicing conditions:

Table[{
s[j, (x[j + 1] + x[j + 2])/2] == s[j + 1, (x[j + 1] + x[j + 2])/2],
ds[j, (x[j + 1] + x[j + 2])/2] == ds[j + 1, (x[j + 1] + x[j + 2])/2]
}, {j, 1, n - 3}]


where ds is first derivative:

ds[j_, xx_] = D[s[j, xx], xx];


Now we combine everything together and convert into matrix form:

eqs = Flatten[{
s[1, x[1]] == y[1],
Table[{
s[j, (x[j + 1] + x[j + 2])/2] == s[j + 1, (x[j + 1] + x[j + 2])/2],
ds[j, (x[j + 1] + x[j + 2])/2] == ds[j + 1, (x[j + 1] + x[j + 2])/2]
}, {j, 1, n - 3}],
s[n - 2, x[n]] == y[n]}];
arrs = Simplify@
Normal@CoefficientArrays[eqs, Flatten@Array[{b[#], c[#]} &, n - 2]];
MatrixForm /@ (4 arrs)


Looking at the matrices it is easy to see that they have periodic structure and contain only elements of the forms x[j+1]-x[j] and y[j+1]-y[j] with some numerical coefficients. This periodic structures can be expressed as SparceArrays.

## Implementation

Assuming that data contains an array of {x[j],y[j]}, we define

Δx[i_] := Subtract @@ data[[{i + 1, i}, 1]];
Δy[i_] := Subtract @@ data[[{i + 1, i}, 2]];


Two matrices in the arrs can be expressed as follows:

bB = SparseArray[{1 -> -4 Δy[1], -1 -> 4 Δy[n - 1], i_ /; OddQ[i] :> 0, i_ /; EvenQ[i] :> 4 Δy[i/2 + 1]}, 2 (n - 2)];
mM = SparseArray[
Join[{{1, 1} -> -4 Δx[1], {1, 2} -> 4 Δx[1]^2, {-1, -1} -> 4 Δx[n - 1]^2, {-1, -2} -> 4 Δx[n - 1]},
Table[Band[{2 i, 2 i - 1}] -> {{2 Δx[i + 1], Δx[i + 1]^2, 2 Δx[i + 1], -Δx[i + 1]^2}, {4, 4 Δx[i + 1], -4, 4 Δx[i + 1]}}, {i, 1, n - 3}]], {2 (n - 2), 2 (n - 2)}];


Now for obtaining coefficients b[i], c[i] we can use LinearSolve:

bcM = LinearSolve[mM, bB];


Grouping coefficients with identical indexes together:

bcM = Partition[bcM, 2];


The intervals at which the parabolas are defined:

intervList =
Partition[
Join[{data[[1, 1]]},
MovingAverage[data[[2 ;; -2, 1]], 2], {data[[-1, 1]]}], 2, 1];


Now we compile all the spline data in one object:

splineData = Table[{data[[i + 1]], bcM[[i]], intervList[[i]]}, {i, n - 2}];


splineData contains everything that is needed to construct the spline. It contains elements of the form:

{{x[i + 1], y[i + 1]}, {b[i], c[i]}, {xmin, xmax}}


It is very handy to use splineData for exact numerical integration (see here).

To produce the explicit piecewise function we define the constructor (HornerForm is already applied and gives 27% speedup):

makeSpline[splineData_List, x_Symbol] :=
Piecewise[
Append[Table[{d[[1, 2]] - d[[1, 1]] d[[2, 1]] +
d[[1, 1]]^2 d[[2, 2]] + x (d[[2, 1]] + x d[[2, 2]] -
2 d[[1, 1]] d[[2, 2]]), #1 <= x <= #2 & @@ d[[3]]}, {d, splineData}],
{Indeterminate, True}]]


It can be used as follows:

spline[\[FormalX]_] = makeSpline[splineData, \[FormalX]];


Compiling the code for creation of splineData into one function:

ClearAll[toSplineData, makeSpline];
Options[toSplineData] = {Method -> Automatic};
toSplineData[data_, OptionsPattern[]] /; MatrixQ[data, NumberQ] :=
Module[{Δx, Δy, bB, mM, bcM, intervList,
n = Length[data], dataS = Sort[data]},
Δx[i_] :=
Subtract @@ dataS[[{i + 1, i}, 1]]; Δy[i_] :=
Subtract @@ dataS[[{i + 1, i}, 2]];
bB = SparseArray[{1 -> -4 Δy[1], -1 ->
4 Δy[n - 1], i_ /; OddQ[i] :> 0,
i_ /; EvenQ[i] :> 4 Δy[i/2 + 1]}, 2 (n - 2)];
mM = SparseArray[
Join[{{1, 1} -> -4 Δx[1], {1, 2} ->
4 Δx[1]^2, {-1, -1} ->
4 Δx[n - 1]^2, {-1, -2} ->
4 Δx[n - 1]},
Table[Band[{2 i,
2 i - 1}] -> {{2 Δx[i + 1], Δx[
i + 1]^2,
2 Δx[i + 1], -Δx[i + 1]^2}, {4,
4 Δx[i + 1], -4,
4 Δx[i + 1]}}, {i, 1, n - 3}]], {2 (n - 2),
2 (n - 2)}];
bcM = LinearSolve[mM, bB, Method -> OptionValue[Method]];
bcM = Partition[bcM, 2];
intervList =
Partition[
Join[{dataS[[1, 1]]},
MovingAverage[dataS[[2 ;; -2, 1]], 2], {dataS[[-1, 1]]}], 2, 1];
Table[{dataS[[i + 1]], bcM[[i]], intervList[[i]]}, {i, n - 2}]];


Now to construct the spline from data we can evaluate

spline[\[FormalX]_] = makeSpline[toSplineData[data], \[FormalX]];


Any suggestions and improvements are welcome!

## What is so special about quadric spline interpolation?

Quadric spline interpolation with splicing points in the middle of successive data points is locally invariant and weakly depends on boundary condition. (But the same is not true for quadric spline interpolation with splicings in data points.) It gives much lesser artifacts than "usual" cubic spline interpolation. Moreover, it seems that splines of higher degree do not give any advantages over such quadric splines. As opposed to cubic spline, data points are not special points at which polynomials are spliced, and it makes much easier exact integration in these points.

• Could you give a demo of your solution? THANKS:) – xyz Oct 10 '15 at 10:38
• @Shutao As you can see from the definition for the parabolas (the first formula in my answer) from which the spline is formed, my solution is non-parametric and so (as well as Interpolation) it is not intended for interpolating an arbitrary dataset (like yours). – Alexey Popkov Oct 10 '15 at 11:40
• @Shutao My solution does not use BSplineBasis and is not intended to be compatible with BSplineFunction and BSplineCurve. But it has an advantage of being very simple and efficient, mathematically clear and very handy when one needs to perform exact integration of interpolated data (without using NIntegrate). Actually I developed my solution exactly for the latter kind of applications. – Alexey Popkov Oct 10 '15 at 11:40
• @Shutao While I appreciate your solution very much I should note that your answer (as well as belisarius' answer) does not answer the original question because the question is explicitly about creation of arbitrary-precision version of Interpolate (preferably with richer set of spline interpolation methods), not about construction of BSplineCurve from a set of points. – Alexey Popkov Oct 10 '15 at 11:40
• @Shutao It would be great if you add a non-parametric version of your solution as another answer: something like BSplineInterpolate[pts][x]: BSplineInterpolate would accept a list of points as the first argument and return something line InterpolatingFunction (in my answer makeSpline returns a Piecewise function object), the latter can accept the abscissa x as its argument and return the interpolated value of the ordinate y. – Alexey Popkov Oct 10 '15 at 11:41

The following is a shameful plug of J. M.'s answer that you already linked to show how to get the parametrization by using BSplineBasis[] with arbitrary precision and packed into a function:

splineInterp[data_, order_, prec_] :=
Module[{parametrizeCurve, tvals, bas, ctrlpts, knots},
parametrizeCurve[pts_ /; MatrixQ[pts, NumericQ], a : (_?NumericQ) : 1/2] :=
FoldList[
Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]];
tvals = parametrizeCurve[data];
(*knots for interpolating B-spline*)
knots =
Join[
ConstantArray[0, order + 1],
ConstantArray[1, order + 1]];
(*basis function matrix*)
bas =
Table[
BSplineBasis[{order, knots}, j - 1, N[tvals[[i]], prec]],
{i, Length[data]}, {j, Length[data]}];
ctrlpts = LinearSolve[bas, data];
Return[
Sum[
ctrlpts[[i + 1]] BSplineBasis[{order, knots}, i, #],
{i, 0, Length[testData] - 1}] &]
]


Usage

SeedRandom[42];
testData = RandomReal[{0, 1}, {10, 2}, WorkingPrecision -> prec];
f = splineInterp[testData, 5, 50];
f[1/10]
ParametricPlot[
f[t], {t, 0, 1},
Axes -> None,
Epilog ->
{Directive[Green, PointSize[Large]], Point[testData]}, Frame -> True
]

(*
{0.06025513038208479326777055464103084790562316, \
0.5927955887343273319037301352154452480202510}
*)


Please refer to J. M.'s answer for the details.

### Explanation

The centripetal distribution:

The function parametrizeCurve complete this algorithm

parametrizeCurve[pts_ /; MatrixQ[pts, NumericQ], a : (_?NumericQ) : 1/2] :=
FoldList[
Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]]

• I do not see in your answer how can I use this as interpolation in a regular way: interpolation[x] where x is not a parametrization variable t but the actual independent variable on which the dependent variable y depends by y[x]. – Alexey Popkov Jan 10 '14 at 16:28
• @AlexeyPopkov, Please see my explanation – xyz May 29 '15 at 13:04
• @ShutaoTang Thank you for the answer to my comment, now I see that the second argument of parametrizeCurve may be arbitrary but must be lesser or equal to $1$. BTW it would be better for others if you would post this explanation under your own answer. :^) – Alexey Popkov May 29 '15 at 13:22
• @AlexeyPopkov, My pleasue. I used the chord-length method to generate the parameters in my answer, so I think it is suitable to place this explanation in belisarius's answer:-) – xyz May 29 '15 at 13:34