.
Lagrange Interpolation§
Before diving into splines, let’s have a look at an arguably simpler interpolation method using polynomials: Lagrange interpolation.
This is easy to implement, but as we will see, it has quite severe limitations, which will motivate us to look into splines later.
[1]:
import matplotlib.pyplot as plt
import numpy as np
One-dimensional Example§
Assume we have \(N\) time instants \(t_i\), with \(0 \le i < N\):
[2]:
ts = -1.5, 0.5, 1.7, 3, 4
… and for each time instant we are given an associated value \(x_i\):
[3]:
xs = 2, -1, 1.3, 3.14, 1
Our task is now to find a function that yields the given \(x_i\) values for the given times \(t_i\) and some “reasonable” interpolated values when evaluated at time values in-between.
The idea of Lagrange interpolation is to create a separate polynomial for each of the \(N\) given time instants, which will be weighted by the associated \(x\). The final interpolation function is the weighted sum of these \(N\) polynomials.
In order for this to actually work, the polynomials must fulfill the following requirements:
Each polynomial must yield \(1\) when evaluated at its associated time \(t_i\).
Each polynomial must yield \(0\) at all other instances in the set of given times.
To satisfy the second point, let’s create a product with a term for each of the relevant times and make each of those factors vanish when evaluated at their associated time. As an example we look at the basis for \(t_3 = 3\):
[4]:
def maybe_polynomial_3(t):
t = np.asarray(t)
return (t - (-1.5)) * (t - 0.5) * (t - 1.7) * (t - 4)
[5]:
maybe_polynomial_3(ts)
[5]:
array([ -0. , 0. , -0. , -14.625, 0. ])
As we can see, this indeed fulfills the second requirement. Note that we were given 5 time instants, but we need only 4 product terms (corresponding to the 4 roots of the polynomial).
Now, for the first requirement, we can divide each term to yield \(1\) when evaluated at \(t = 3\) (luckily, this will not violate the second requirement). If each term is \(1\), the whole product will also be \(1\):
[6]:
def polynomial_3(t):
t = np.asarray(t)
return (
(t - (-1.5)) / (3 - (-1.5)) *
(t - 0.5) / (3 - 0.5) *
(t - 1.7) / (3 - 1.7) *
(t - 4) / (3 - 4))
[7]:
polynomial_3(ts)
[7]:
array([ 0., -0., 0., 1., -0.])
That’s it!
To get a better idea what’s going on between the given time instances, let’s plot this polynomial (with a little help from helper.py):
[8]:
from helper import grid_lines
[9]:
plot_times = np.linspace(ts[0], ts[-1], 100)
[10]:
plt.plot(plot_times, polynomial_3(plot_times))
grid_lines(ts, [0, 1])
We can see from its shape that this is a polynomial of degree 4, which makes sense because the product we are using has 4 terms containing one \(t\) each. We can also see that it has the value \(0\) at each of the initially provided time instances \(t_i\), except for \(t_3 = 3\), where it has the value \(1\).
The above calculation can be easily generalized to be able to get any one of the set of polynomials defined by an arbitrary list of time instants:
[11]:
def lagrange_polynomial(times, i, t):
"""i-th Lagrange polynomial for the given time values, evaluated at t."""
t = np.asarray(t)
product = np.multiply.reduce
return product([
(t - times[j]) / (times[i] - times[j])
for j in range(len(times))
if i != j
])
Now we can calculate and visualize all 5 polynomials for our 5 given time instants:
[12]:
polys = np.column_stack([lagrange_polynomial(ts, i, plot_times)
for i in range(len(ts))])
[13]:
plt.plot(plot_times, polys)
grid_lines(ts, [0, 1])
Finally, the interpolated values can be obtained by applying the given \(x_i\) values as weights to the polynomials and summing everything together:
[14]:
weighted_polys = polys * xs
[15]:
interpolated = np.sum(weighted_polys, axis=-1)
[16]:
plt.plot(plot_times, weighted_polys)
plt.plot(plot_times, interpolated, color='black', linestyle='dashed')
plt.scatter(ts, xs, color='black')
grid_lines(ts)
Neville’s Algorithm§
An alternative way to calculate interpolated values is Neville’s algorithm (see also [BG88], figure 2). We mention this algorithm mainly because it is referenced in the derivation of non-uniform Catmull–Rom splines and the description of the Barry–Goldman algorithm.
As main building block, we need a linear interpolation between two values in a given time interval:
[17]:
def lerp(xs, ts, t):
"""Linear interpolation.
Returns the interpolated value at time *t*,
given the two values *xs* at times *ts*.
"""
x_begin, x_end = xs
t_begin, t_end = ts
return (x_begin * (t_end - t) + x_end * (t - t_begin)) / (t_end - t_begin)
In each stage of the algorithm, linear interpolation is used to interpolate between adjacent values, leading to one fewer value than in the stage before. The new values are used as input to the next stage and so on. When there is only one value left, this value is the result.
The only tricky part is to choose the appropriate time interval for each interpolation. In the first stage, the intervals between the given time values are used. In the second stage, each time interval is combined with the following one, leading to one fewer time intervals in total. In the third stage, each time interval is combined with the following two intervals, and so on until the last stage, where all time intervals are combined into a single large interval.
[18]:
def neville(xs, ts, t):
"""Lagrange interpolation using Neville's algorithm.
Returns the interpolated value at time(s) *t*,
given the values *xs* at times *ts*.
"""
assert len(xs) == len(ts)
if not np.isscalar(t):
return np.array([neville(xs, ts, time) for time in t])
while len(xs) > 1:
step = len(ts) - len(xs) + 1
xs = [
lerp(*args, t)
for args in zip(zip(xs, xs[1:]), zip(ts, ts[step:]))]
return xs[0]
[19]:
plt.plot(plot_times, neville(xs, ts, plot_times))
plt.scatter(ts, xs)
grid_lines(ts)
Two-dimensional Example§
Lagrange interpolation can of course also be used in higher-dimensional spaces. To show this, let’s create a class:
[20]:
class Lagrange:
def __init__(self, vertices, grid):
assert len(vertices) == len(grid)
self.vertices = np.array(vertices)
self.grid = list(grid)
def evaluate(self, t):
# Alternatively, we could simply use this one-liner:
#return neville(self.vertices, self.grid, t)
if not np.isscalar(t):
return np.array([self.evaluate(time) for time in t])
polys = [lagrange_polynomial(self.grid, i, t)
for i in range(len(self.grid))]
weighted_polys = self.vertices.T * polys
return np.sum(weighted_polys, axis=-1)
Since this class has the same interface as the splines that are discussed in the following sections, we can use a spline helper function from helper.py for plotting:
[21]:
from helper import plot_spline_2d
This time, we have a list of two-dimensional vectors and the same list of associated times as before:
[22]:
l1 = Lagrange([(2, -3), (-1, 0), (1.3, 1), (3.14, 0), (1, -1)], ts)
[23]:
plot_spline_2d(l1)
Runge’s Phenomenon§
This seems to work quite well, but as indicated above, Lagrange implementation has a severe limitation. This limitation gets more apparent when using more vertices, which leads to a higher-degree polynomial.
[24]:
vertices = [
(1, 0),
(1, 2),
(3, 0),
(2, -1),
(2.5, 1.5),
(5, 2),
(6, 1),
(5, 0),
(6, -2),
(7, 2),
(4, 4),
]
times = range(len(vertices))
[25]:
l2 = Lagrange(vertices, times)
plot_spline_2d(l2)
Here we see a severe overshooting effect, most pronounced at the beginning and the end of the curve. This effect is called Runge’s phenomenon.
Long story short, Lagrange interpolation is typically not usable for drawing curves. For comparison, let’s use the same positions and time values and create a Catmull–Rom spline:
[26]:
import splines
[27]:
cr_spline = splines.CatmullRom(vertices, times)
[28]:
plot_spline_2d(cr_spline)
This clearly doesn’t have the overshooting problem we saw above.
Note
The splines.CatmullRom class uses “natural” end conditions by default.