# 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.

```
>>> import numpy as np
>>> 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
<Compressed Sparse Row sparse array of dtype 'float64'
with 12 stored elements and shape (3, 7)>
```

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. ]]
```