# Arbitrary precision spline interpolation

Current implementation of Interpolation does not allow arbitrary precision spline interpolation. Yu-Sung Chang say 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 are the names of such spline and of this kind of parametrization?). Also Ｊ. Ｍ. gives in this answer 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? –  belisarius Jan 10 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 at 12:51

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!

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.

-

The following is a shameful plug of Ｊ. Ｍ. 's answer 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],  MovingAverage[ArrayPad[tvals, -1], order],
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}
*)


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 at 16:28