Interpolation

Interpolation#

Interpolation means to fill in a function between known values. The data for interpolation are a set of points x and a set of function values y, and the result is a function f from some function class so that f(x) = y. Typically this function class is something simple, like Polynomials of bounded degree, piecewise constant functions, or splines.

Let’s look at a few examples

import numpy as np
import matplotlib.pyplot as plt

Here are a couple of functions we’ll look at:

fs = [
    lambda x : (x - 3) * (x + 3) * x, # cubic
    lambda x : np.exp(-x**2 / 2),     # gaussian
    lambda x : np.sin(3*x) / (3*x),       # sinc function
    lambda x : 1 / (np.exp(-2*x) + 1)   # logistic
]

x = np.linspace(-4,4,500)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
plt.show(fig)

../_images/b0461333103739b68f4aec734e396099a516fdf5b8f539d9038e90a3c6bb602e.png

If we sample a couple of points on each function, we would see something like the following:

x = np.linspace(-4,4,500)
x6 = np.linspace(-4,4,6)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
    y = f(x6)
    ax[i].scatter(x6, y)
    
plt.show(fig)

../_images/cd33f8e3c86bd90ca393cb3645aa5a0aa02eec711ed0eeeb15d1c28b6f578c24.png

An interpolation just uses the sampled points and function values to try to reconstruct the original function.

Scipy provides a high-level interface for doing this with scipy.interpolate.interp1d. You can pass in a variety of arguments, and the output is an interpolation object, which implements the __call__ method, so you can evaluate the interpolation.

from scipy.interpolate import interp1d

x = np.linspace(-4,4,500)
xk = np.linspace(-4,4,6)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = interp1d(xk, yk, kind='linear')
    ax[i].plot(x, interp(x), '--')
    
    
    
plt.show(fig)

../_images/86767dce448de63a8104cf2e483edb8651c8d332a722a07fb7af69bb35a139b3.png

If we increase the number of points used for interpolation, then we can get better approximations to the function.

x = np.linspace(-4,4,500)
xk = np.linspace(-4,4,10)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = interp1d(xk, yk, kind='linear')
    ax[i].plot(x, interp(x), '--')
    
    
    
plt.show(fig)

../_images/06af013e395530ee13b98e2b162fd3af8b80fa324233a168be66caaa5ba30bfe.png

There are a variety of options for the kind keyword.

kind='linear': piecewise linear interpolation
kind='nearest': nearest-neighbor interpolation (piecewise constant)
kind='quadratic': quadratic spline (piecewise quadratic)
kind='cubic': cubic spline (piecewise cubic)

and more. See the documentation for details

kinds = ('linear', 'nearest', 'quadratic', 'cubic')
f = fs[2]

x = np.linspace(-4,4,500)
xk = np.linspace(-4,4,10)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, kind in enumerate(kinds):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = interp1d(xk, yk, kind=kind)
    ax[i].plot(x, interp(x), '--')
    ax[i].set_title(kind)
    
    
    
plt.show(fig)

../_images/ee77ba0b043d682554864a41c8145789814aced37e83d5aef3753383f0de5ff1.png

Barycentric Interpolation#

There are also several interpolation methods not covered by interp1d, or which offer more customization. These include:

BarycentricInterpolator
KroghInterpolator
PchipInterpolator

and so on. See the interpolation documentation for more.

Barycentric interpolation can be done using BarycentricInterpolator, which constructs a Lagrange Polynomial of minimum degree which interpolates the points.

See Berrut and Trefethen, 2004 for details.

from scipy.interpolate import BarycentricInterpolator

x = np.linspace(-4,4,500)
xk = np.linspace(-4,4,10)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = BarycentricInterpolator(xk, yk)
    ax[i].plot(x, interp(x), '--')
    
    
    
plt.show(fig)

../_images/1a958b761f72890fd8f38f002a2cf58d794b2e2cd17caf4ab1facaebb1f2d062.png

a big difference between a Barycentric interpolator and the sorts of interpolants you contruct using interp1d is that it is not defined piecewise - it is a degree k-1 polynomial which fits each of the k interpolation points exactly. If the function being interpolated is not a lower-degree polynomial, this can lead to large fluctuations away from the origin, which is called the Runge phenomenon.

You can minimize the effects of the Runge phenomenon by using Chebyshev interpolation points: cos(i*np.pi/(k-1)), i=0,...,k-1

def chebyshev(k, scale=1):
    """
    return k Chebyshev interpolation points in the range [-scale, scale]
    """
    return scale*np.cos(np.arange(k) * np.pi / (k-1))

xk = chebyshev(5, scale=4)
xk

array([ 4.00000000e+00,  2.82842712e+00,  2.44929360e-16, -2.82842712e+00,
       -4.00000000e+00])

x = np.linspace(-4,4,500)
xk = chebyshev(10, scale=4)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = BarycentricInterpolator(xk, yk)
    ax[i].plot(x, interp(x), '--')
    
    
plt.show(fig)

../_images/dc87e553b8440c98d5ab23c673911537c05f87d6e65870f2207fda3552b550dc.png

Higher-Dimensional Interpolation#

Higher dimensional interpolations can be computed using either scipy.interpolate.interp2d (for 2-dimensions), or scipy.interpolate.griddata (for any number of dimensions). There are also a variety of more specific options which you can find in the documentation.

The following function is found in the SciPy interpolation tutorial

f = lambda x, y : x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2

We’ll want to interpolate onto a grid covering the unit square. The following two cells generate equivalent data:

gx, gy = np.mgrid[0:1:100j, 0:1:100j]

gx, gy = np.meshgrid(np.linspace(0,1,100), np.linspace(0,1,100))

We’ll sample the function at some random points

from scipy.interpolate import griddata

n = 500
x = np.random.rand(n, 2)
z = f(x[:,0], x[:,1])

fig, ax = plt.subplots(2, 2, figsize=(10,10))
ax = ax.flatten()

# keyword arguments for imshow
imkw = {'extent': (0,1,0,1), 'origin': 'lower'}

ax[0].scatter(x[:,0], x[:,1], marker='x', color='k')
ax[0].imshow(f(gx, gy), **imkw)
ax[0].set_title('samples')

for i, method in enumerate(['nearest', 'linear', 'cubic']):
    grid_interpolation = griddata(x, z, (gx, gy), method=method, fill_value=0)
    ax[i+1].imshow(grid_interpolation, **imkw)
    ax[i+1].set_title(method)

    
    
plt.show(fig)

../_images/85c6446a40f6209962b98749467e84fd1e5120b964b87eeb6b41fd578d68bd94.png

interp2d should be used with a grid of points both when fitting the interpolant and when evaluating

from scipy.interpolate import interp2d

gx, gy = np.mgrid[0:1:20j, 0:1:20j]
z = f(gx, gy)

fig, ax = plt.subplots(2, 2, figsize=(10,10))
ax = ax.flatten()

# keyword arguments for imshow
imkw = {'extent': (0,1,0,1), 'origin': 'lower'}

ax[0].scatter(gx.flatten(), gy.flatten(), marker='x', color='k')
ax[0].imshow(z.T, **imkw)
ax[0].set_title('samples')

for i, method in enumerate(['linear', 'cubic', 'quintic']):
    interp = interp2d(gx, gy, z, kind=method, fill_value=0)
    grid_interpolation = interp(np.linspace(0,1,100), np.linspace(0,1,100))
    ax[i+1].imshow(grid_interpolation, **imkw)
    ax[i+1].set_title(method)

    
    
plt.show(fig)

/tmp/ipykernel_1229/3204843088.py:17: DeprecationWarning: `interp2d` is deprecated!
`interp2d` is deprecated in SciPy 1.10 and will be removed in SciPy 1.13.0.

For legacy code, nearly bug-for-bug compatible replacements are
`RectBivariateSpline` on regular grids, and `bisplrep`/`bisplev` for
scattered 2D data.

In new code, for regular grids use `RegularGridInterpolator` instead.
For scattered data, prefer `LinearNDInterpolator` or
`CloughTocher2DInterpolator`.

For more details see
`https://scipy.github.io/devdocs/notebooks/interp_transition_guide.html`

  interp = interp2d(gx, gy, z, kind=method, fill_value=0)
/tmp/ipykernel_1229/3204843088.py:18: DeprecationWarning:         `interp2d` is deprecated!
        `interp2d` is deprecated in SciPy 1.10 and will be removed in SciPy 1.13.0.

        For legacy code, nearly bug-for-bug compatible replacements are
        `RectBivariateSpline` on regular grids, and `bisplrep`/`bisplev` for
        scattered 2D data.

        In new code, for regular grids use `RegularGridInterpolator` instead.
        For scattered data, prefer `LinearNDInterpolator` or
        `CloughTocher2DInterpolator`.

        For more details see
        `https://scipy.github.io/devdocs/notebooks/interp_transition_guide.html`

  grid_interpolation = interp(np.linspace(0,1,100), np.linspace(0,1,100))

../_images/d0a7a4b7b084d5f270f85017f8ca298e1fbb74787a237489536fb1277d3063cb.png

Extrapolate#

After SciPy version 0.17.0, scipy.interpolate.interp1d allows extrapolation, by setting fill_value='extrapolate'. Please see documentations here for reference.

x = np.linspace(-10,10,500)
xk = np.linspace(-4,4,10)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = interp1d(xk, yk, kind='linear', fill_value='extrapolate')
    ax[i].plot(x, interp(x), '--')
    
    
    
plt.show(fig)

../_images/1928fa8b47d71ea260f1a71c866983f340014fa77de4c5a1d851a1c409103549.png

kinds = ('linear', 'nearest', 'quadratic', 'cubic')
f = fs[2]

x = np.linspace(-10,10,500)
xk = np.linspace(-4,4,10)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, kind in enumerate(kinds):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = interp1d(xk, yk, kind=kind, fill_value='extrapolate')
    ax[i].plot(x, interp(x), '--')
    ax[i].set_title(kind)
    
    
plt.show(fig)

../_images/a130726dd97db5c048c9088e6531b1c1589131836aed007113e4b1914cec6c95.png

Barycentric interpolation with BarycentricInterpolator natively supports extrapolation.

x = np.linspace(-8,8,500)
xk = chebyshev(10, scale=4)

fig, ax = plt.subplots(2, 2, figsize=(8,8))
ax = ax.flatten()

for i, f in enumerate(fs):
    y = f(x)
    ax[i].plot(x, y)
    
    yk = f(xk)
    ax[i].scatter(xk, yk)
    
    interp = BarycentricInterpolator(xk, yk)
    ax[i].plot(x, interp(x), '--')
    
    
plt.show(fig)

../_images/9fd1f85427f2d307dbc1416e01da07baab13f2a8b457dc6e6af2e2e75d1d93a2.png

Interpolation

Contents

Interpolation#

Barycentric Interpolation#

Higher-Dimensional Interpolation#

Extrapolate#