4

I tried the following:

spline= interpolate.InterpolatedUnivariateSpline(X, Y, k=3)
coefs= spline.get_coeffs()

With five values in each of X and Y, I ended up with coefs also having five values. Given that five data points implies four spline sections, and that a cubic polynomial has four coefficients, I would have expected to get four times four= 16 coefficients. Does anyone know how to interpret the values that are returned by the get_coeffs method? Is there any place where this is documented?

  • Somewhat relevant: http://stackoverflow.com/questions/22488637/getting-spline-equation-from-univariatespline-object/25330648#25330648 – Jonathan March Feb 18 '15 at 00:49
  • and http://stackoverflow.com/questions/13384859/coefficients-of-spline-interpolation-in-scipy – Jonathan March Feb 18 '15 at 00:53
  • and http://scipy-user.10969.n7.nabble.com/help-interpreting-univariate-spline-td12857.html – Jonathan March Feb 18 '15 at 00:54
  • They are the B-spline coefficients, i.e., `spline(x) == sum(coefs[k]*B(knots,k,x) for k in range(5))`. The `B` functions are the B-spline basis set corresponding to your set of knots; the boundary conditions are IIRC handled by inserting coinciding knots at the edges as necessary. What B-splines are is explained in the books listed http://docs.scipy.org/doc/scipy-dev/reference/generated/scipy.interpolate.splrep.html – pv. Feb 18 '15 at 01:05

1 Answers1

2

These are not the coefficients of x, x**2, and so forth: such monomials are ill-suited to representing splines. Rather, they are coefficients of B-splines which are computed for the specific grid on which interpolation is done. The number of B-splines involved is equal to the number of data points, and so is the number of coefficients. As an example, suppose we want to interpolate these data:

xv = [0, 1, 2, 3, 4, 5]
yv = [3, 6, 5, 7, 9, 1]

Begin with the simpler case of degree k=1 (piecewise linear spline). Then the B-splines are these "triangular hat" functions:

hats

There are 6 of them. Each of them is equal to 1 at "its" grid point, and 0 at all other grid points. This makes it really easy to write the interpolating spline: it is y[0]*b[0] + ... + y[5]*b[5]. And indeed, get_coeffs shows the coefficients are the the y-values themselves.

InterpolatedUnivariateSpline(xv, yv, k=1).get_coeffs()  # [ 3.,  6.,  5.,  7.,  9.,  1.]

Cubic splines

Now it gets hairy, because we need "hats" that are smooth, rather than pointy as those above. The smoothness requirement forces them to be wider, so each B-spline has nonzero values on several grid points. (Technicality: a cubic B-spline has nonzero values at up to 3 knots, but on the chart below, 1 and 4, despite being grid points, are not knots due to so-called "not a knot" condition. Never mind this.) Here are the B-splines for our grid:

b-splines

To get these, I used older methods splrep and splev of scipy.interpolate, which call the same fitpack routines under the hood. The coefficient vector here is the second entry of the tuple tck; I modify it to have one 1 and the rest 0, thus creating a basis spline (b-spline).

k = 3
tck = splrep(xv, yv, s=0, k=k)
xx = np.linspace(min(xv), max(xv), 500)
bsplines = []
for j in range(len(xv)):
    tck_mod = (tck[0], np.arange(len(xv)+2*k-2) == j, k)
    bsplines.append(splev(xx, tck_mod))
    plt.plot(xx, bsplines[-1])

Now that we have a list bsplines, we can use the coefficients returned by get_coeffs to put them together ourselves into an interpolating spline:

coeffs = InterpolatedUnivariateSpline(xv, yv, k=3).get_coeffs()
interp_spline = sum([coeff*bspline for coeff, bspline in zip(coeffs, bsplines)])
plt.plot(xx, interp_spline)

interpolation

If you want a formula for the pieces of these B-splines, the Cox-de Boor recursion formula on B-splines can help but these are a chore to compute by hand.

SymPy can give formulas for B-splines, but there is a little twist. One should pass in a padded set of knots, by repeating the end values like

[0, 0, 0, 0, 2, 3, 5, 5, 5, 5]

This is because at 0 and 5 all four coefficients change values, while at 1 and 4 none of them do, so they are omitted ("not a knot"). (Also, the current version of SymPy (1.1.1) has an issue with repeated knots, but this will be fixed in the next version.)

from sympy import symbols, bspline_basis_set, plot
x = symbols('x')
xv_padded = [0, 0, 0, 0, 2, 3, 5, 5, 5, 5] 
bs = bspline_basis_set(3, xv_padded, x)

Now bs is an array of scary-looking piecewise formulas:

[Piecewise((-x**3/8 + 3*x**2/4 - 3*x/2 + 1, (x >= 0) & (x <= 2)), (0, True)), 
Piecewise((19*x**3/72 - 5*x**2/4 + 3*x/2, (x >= 0) & (x <= 2)), (-x**3/9 + x**2 - 3*x + 3, (x >= 2) & (x <= 3)), (0, True)), 
Piecewise((-31*x**3/180 + x**2/2, (x >= 0) & (x <= 2)), (11*x**3/45 - 2*x**2 + 5*x - 10/3, (x >= 2) & (x <= 3)), (-x**3/30 + x**2/2 - 5*x/2 + 25/6, (x >= 3) & (x <= 5)), (0, True)), 
Piecewise((x**3/30, (x >= 0) & (x <= 2)), (-11*x**3/45 + 5*x**2/3 - 10*x/3 + 20/9, (x >= 2) & (x <= 3)), (31*x**3/180 - 25*x**2/12 + 95*x/12 - 325/36, (x >= 3) & (x <= 5)), (0, True)), 
Piecewise((x**3/9 - 2*x**2/3 + 4*x/3 - 8/9, (x >= 2) & (x <= 3)), (-19*x**3/72 + 65*x**2/24 - 211*x/24 + 665/72, (x >= 3) & (x <= 5)), (0, True)), 
Piecewise((x**3/8 - 9*x**2/8 + 27*x/8 - 27/8, (x >= 3) & (x <= 5)), (0, True))]