Piecewise polynomials and splines#
1D interpolation routines discussed in the previous section, work by constructing certain piecewise polynomials: the interpolation range is split into intervals by the so-called breakpoints, and there is a certain polynomial on each interval. These polynomial pieces then match at the breakpoints with a predefined smoothness: the second derivatives for cubic splines, the first derivatives for monotone interpolants and so on.
A polynomial of degree \(k\) can be thought of as a linear combination of
\(k+1\) monomial basis elements, \(1, x, x^2, \cdots, x^k\).
In some applications, it is useful to consider alternative (if formally
equivalent) bases. Two popular bases, implemented in scipy.interpolate are
B-splines (BSpline) and Bernstein polynomials (BPoly).
B-splines are often used for, for example, non-parametric regression problems,
and Bernstein polynomials are used for constructing Bezier curves.
PPoly objects represent piecewise polynomials in the ‘usual’ power basis.
This is the case for CubicSpline instances and monotone interpolants.
In general, PPoly objects can represent polynomials of
arbitrary orders, not only cubics. For the data array x, breakpoints are at
the data points, and the array of coefficients, c , define polynomials of
degree \(k\), such that c[i, j] is a coefficient for
(x - x[j])**(k-i) on the segment between x[j] and x[j+1] .
BSpline objects represent B-spline functions — linear combinations of
b-spline basis elements.
These objects can be instantiated directly or constructed from data with the
make_interp_spline factory function.
Finally, Bernstein polynomials are represented as instances of the BPoly class.
All these classes implement a (mostly) similar interface, PPoly being the most
feature-complete. We next consider the main features of this interface and
discuss some details of the alternative bases for piecewise polynomials.
Manipulating PPoly objects#
PPoly objects have convenient methods for constructing derivatives
and antiderivatives, computing integrals and root-finding. For example, we
tabulate the sine function and find the roots of its derivative.
>>> from scipy.interpolate import CubicSpline
>>> x = np.linspace(0, 10, 71)
>>> y = np.sin(x)
>>> spl = CubicSpline(x, y)
Now, differentiate the spline:
>>> dspl = spl.derivative()
Here dspl is a PPoly instance which represents a polynomial approximation
to the derivative of the original object, spl . Evaluating dspl at a
fixed argument is equivalent to evaluating the original spline with the nu=1
argument:
>>> dspl(1.1), spl(1.1, nu=1)
(0.45361436, 0.45361436)
Note that the second form above evaluates the derivative in place, while with
the dspl object, we can find the zeros of the derivative of spl:
>>> dspl.roots() / np.pi
array([-0.45480801, 0.50000034, 1.50000099, 2.5000016 , 3.46249993])
This agrees well with roots \(\pi/2 + \pi\,n\) of \(\cos(x) = \sin'(x)\). Note that by default it computed the roots extrapolated to the outside of the interpolation interval \(0 \leqslant x \leqslant 10\), and that the extrapolated results (the first and last values) are much less accurate. We can switch off the extrapolation and limit the root-finding to the interpolation interval:
>>> dspl.roots(extrapolate=False) / np.pi
array([0.50000034, 1.50000099, 2.5000016])
In fact, the root method is a special case of a more general solve
method which finds for a given constant \(y\) the solutions of the
equation \(f(x) = y\) , where \(f(x)\) is the piecewise polynomial:
>>> dspl.solve(0.5, extrapolate=False) / np.pi
array([0.33332755, 1.66667195, 2.3333271])
which agrees well with the expected values of \(\pm\arccos(1/2) + 2\pi\,n\).
Integrals of piecewise polynomials can be computed using the .integrate
method which accepts the lower and the upper limits of integration. As an
example, we compute an approximation to the complete elliptic integral
\(K(m) = \int_0^{\pi/2} [1 - m\sin^2 x]^{-1/2} dx\):
>>> from scipy.special import ellipk
>>> m = 0.5
>>> ellipk(m)
1.8540746773013719
To this end, we tabulate the integrand and interpolate it using the monotone
PCHIP interpolant (we could as well used a CubicSpline):
>>> from scipy.interpolate import PchipInterpolator
>>> x = np.linspace(0, np.pi/2, 70)
>>> y = (1 - m*np.sin(x)**2))**(-1/2)
>>> spl = PchipInterpolator(x, y)
and integrate
>>> spl.integrate(0, np.pi/2)
1.854074674965991
which is indeed close to the value computed by scipy.special.ellipk.
All piecewise polynomials can be constructed with N-dimensional y values.
If y.ndim > 1, it is understood as a stack of 1D y values, which are
arranged along the interpolation axis (with the default value of 0).
The latter is specified via the axis argument, and the invariant is that
len(x) == y.shape[axis]. As an example, we extend the elliptic integral
example above to compute the approximation for a range of m values, using
the NumPy broadcasting:
>>> from scipy.interpolate import PchipInterpolator
>>> m = np.linspace(0, 0.9, 11)
>>> x = np.linspace(0, np.pi/2, 70)
>>> y = 1 / np.sqrt(1 - m[:, None]*np.sin(x)**2)
Now the y array has the shape (11, 70), so that the values of y
for fixed value of m are along the second axis of the y array.
>>> spl = PchipInterpolator(x, y, axis=1) # the default is axis=0
>>> import matplotlib.pyplot as plt
>>> plt.plot(m, spl.integrate(0, np.pi/2), '--')
>>> from scipy.special import ellipk
>>> plt.plot(m, ellipk(m), 'o')
>>> plt.legend(['`ellipk`', 'integrated piecewise polynomial'])
>>> plt.show()
B-splines: knots and coefficients#
A b-spline function — for instance, constructed from data via a
make_interp_spline call — is defined by the so-called knots and coefficients.
As an illustration, let us again construct the interpolation of a sine function.
The knots are available as the t attribute of a BSpline instance:
>>> x = np.linspace(0, 3/2, 7)
>>> y = np.sin(np.pi*x)
>>> from scipy.interpolate import make_interp_spline
>>> bspl = make_interp_spline(x, y, k=3)
>>> print(bspl.t)
[0. 0. 0. 0. 0.5 0.75 1. 1.5 1.5 1.5 1.5 ]
>>> print(x)
[ 0. 0.25 0.5 0.75 1. 1.25 1.5 ]
We see that the knot vector by default is constructed from the input
array x: first, it is made \((k+1)\) -regular (it has k
repeated knots appended and prepended); then, the second and
second-to-last points of the input array are removed—this is the so-called
not-a-knot boundary condition.
In general, an interpolating spline of degree k needs
len(t) - len(x) - k - 1 boundary conditions. For cubic splines with
(k+1)-regular knot arrays this means two boundary conditions—or
removing two values from the x array. Various boundary conditions can be
requested using the optional bc_type argument of make_interp_spline.
The b-spline coefficients are accessed via the c attribute of a BSpline
object:
>>> len(bspl.c)
7
The convention is that for len(t) knots there are len(t) - k - 1
coefficients. Some routines (see the Smoothing splines section) zero-pad the c arrays so that
len(c) == len(t). These additional coefficients are ignored for evaluation.
We stress that the coefficients are given in the b-spline basis, not the power basis of \(1, x, \cdots, x^k\).
B-spline basis elements#
B-splines are piecewise polynomials, represented as linear combinations of b-spline basis elements — which themselves are certain linear combinations of usual monomials, \(x^m\) with \(m=0, 1, \dots, k\).
The b-spline basis is generally more computationally stable than the power basis and is useful for a variety of applications which include interpolation, regression and curve representation. The main feature is that these basis elements are localized and equal to zero outside of an interval defined by the knot array.
Specifically, a b-spline basis element of degree k (e.g. k=3 for cubics)
is defined by \(k+2\) knots and is zero outside of these knots.
To illustrate, plot a collection of non-zero basis elements on a certain
interval:
>>> k = 3 # cubic splines
>>> t = [0., 1.4, 2., 3.1, 5.] # internal knots
>>> t = np.r_[[0]*k, t, [5]*k] # add boundary knots
>>> from scipy.interpolate import BSpline
>>> import matplotlib.pyplot as plt
>>> for j in [-2, -1, 0, 1, 2]:
... a, b = t[k+j], t[-k+j-1]
... xx = np.linspace(a, b, 101)
... bspl = BSpline.basis_element(t[k+j:-k+j])
... plt.plot(xx, bspl(xx), label=f'j = {j}')
>>> plt.legend(loc='best')
>>> plt.show()
Here BSpline.basis_element is essentially a shorthand for constructing a spline
with only a single non-zero coefficient. For instance, the j=2 element in
the above example is equivalent to
>>> c = np.zeros(t.size - k - 1)
>>> c[-2] = 1
>>> b = BSpline(t, c, k)
>>> np.allclose(b(xx), bspl(xx))
True
If desired, a b-spline can be converted into a PPoly object using
PPoly.from_spline method which accepts a BSpline instance and returns a
PPoly instance. The reverse conversion is performed by the
BSpline.from_power_basis method. However, conversions between bases is best
avoided because it accumulates rounding errors.
Design matrices in the B-spline basis#
One common application of b-splines is in non-parametric regression. The reason is that the localized nature of the b-spline basis elements makes linear algebra banded. This is because at most \(k+1\) basis elements are non-zero at a given evaluation point, thus a design matrix built on b-splines has at most \(k+1\) diagonals.
As an illustration, we consider a toy example. Suppose our data are one-dimensional and are confined to an interval \([0, 6]\). We construct a 4-regular knot vector which corresponds to 7 data points and cubic, k=3, splines:
>>> t = [0., 0., 0., 0., 2., 3., 4., 6., 6., 6., 6.]
Next, take ‘observations’ to be
>>> xnew = [1, 2, 3]
and construct the design matrix in the sparse CSR format
>>> from scipy.interpolate import BSpline
>>> mat = BSpline.design_matrix(xnew, t, k=3)
>>> mat
<3x7 sparse array of type '<class 'numpy.float64'>'
with 12 stored elements in Compressed Sparse Row format>
Here each row of the design matrix corresponds to a value in the xnew array,
and a row has no more than k+1 = 4 non-zero elements; row j
contains basis elements evaluated at xnew[j]:
>>> with np.printoptions(precision=3):
... print(mat.toarray())
[[0.125 0.514 0.319 0.042 0. 0. 0. ]
[0. 0.111 0.556 0.333 0. 0. 0. ]
[0. 0. 0.125 0.75 0.125 0. 0. ]]