This page was generated from doc/euclidean/lagrange.ipynb. Interactive online version: .

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.

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\) …

ts = -1.5, 0.5, 1.7, 3.5, 4

… and for each time instant we are given an associated value \(x_i\):

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 \(\ell_i(t)\) for each of the \(N\) given time instants, which will be weighted by the associated \(x_i\). The final interpolation function is the weighted sum of these \(N\) polynomials:

\begin{equation*} L(t) = \sum_{i=0}^{N-1} x_i \ell_i(t) \end{equation*}

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. For example, let’s look at the basis for \(i = 3\):

def maybe_polynomial_3(t):
    t = np.asarray(t)
    return (
        (t - (-1.5)) *
        (t - 0.5) *
        (t - 1.7) *
        (t - 4))

maybe_polynomial_3(ts)

array([ -0. ,   0. ,  -0. , -13.5,   0. ])

As we can see, this indeed fulfills the second requirement. Note that we were given 5 time instants, but we only need 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 = t_3 = 3.5\) (luckily, this will not violate the second requirement). If each term is \(1\), the whole product will also be \(1\):

def polynomial_3(t):
    t = np.asarray(t)
    return (
        (t - (-1.5)) / (3.5 - (-1.5)) *
        (t - 0.5) / (3.5 - 0.5) *
        (t - 1.7) / (3.5 - 1.7) *
        (t - 4) / (3.5 - 4))

polynomial_3(ts)

array([ 0., -0.,  0.,  1., -0.])

That’s it!

To get a better idea what’s going on between the given time instances \(t_i\), let’s plot this polynomial (with a little help from helper.py):

from helper import grid_lines

plot_times = np.linspace(ts[0], ts[-1], 100)

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.5\), 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:

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
    ])

Putting this in mathematic notation, Lagrange basis polynomials can be written as

\begin{equation*} \ell_i(t) = \prod_{\substack{j=0\\i \ne j}}^{N-1} \frac{t - t_j}{t_i - t_j}. \end{equation*}

Now we can calculate and visualize all 5 basis polynomials for our 5 given time instants:

polys = np.column_stack(
    [lagrange_polynomial(ts, i, plot_times) for i in range(len(ts))])

plt.plot(plot_times, polys)
grid_lines(ts, [0, 1])

Finally, the interpolated values \(L(t)\) can be obtained by applying the given \(x_i\) values as weights to the polynomials \(\ell_i(t)\) and summing everything up together:

weighted_polys = polys * xs

interpolated = np.sum(weighted_polys, axis=-1)

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

def lerp(xs, ts, t):
    """Linear intERPolation.

    Returns the interpolated value(s) at time(s) *t*,
    given two values/vertices *xs* at times *ts*.

    The two x-values can be scalars or vectors,
    or even higher-dimensional arrays
    (as long as the shape of *t* is compatible).

    """
    x_begin, x_end = map(np.asarray, xs)
    t_begin, t_end = ts
    if not np.isscalar(t):
        # This allows using an array of *t* values:
        t = np.expand_dims(t, axis=-1)
    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 less 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 less 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.

Barry and Goldman [BG88] show (in figure 2) the cubic case, which looks something like this:

\begin{equation*} \def\negspace{\!\!\!\!\!\!} \begin{array}{ccccccccccccc} &&&&&& \boldsymbol{p}_{0,1,2,3} &&&&&& \\ &&&&& \negspace \frac{t_3 - t}{t_3 - t_0} \negspace && \negspace \frac{t - t_0}{t_3 - t_0} \negspace &&&&& \\ &&&& \boldsymbol{p}_{0,1,2} &&&& \boldsymbol{p}_{1,2,3} &&&& \\ && & \negspace \frac{t_2 - t}{t_2 - t_0} \negspace && \negspace \frac{t - t_0}{t_2 - t_0} \negspace & & \negspace \frac{t_3 - t}{t_3 - t_1} \negspace && \negspace \frac{t - t_1}{t_3 - t_1} \negspace & && \\ && \boldsymbol{p}_{0,1} &&&& \boldsymbol{p}_{1,2} &&&& \boldsymbol{p}_{2,3} && \\ & \negspace \frac{t_1 - t}{t_1 - t_0} \negspace && \negspace \frac{t - t_0}{t_1 - t_0} \negspace & & \negspace \frac{t_2 - t}{t_2 - t_1} \negspace && \negspace \frac{t - t_1}{t_2 - t_1} \negspace & & \negspace \frac{t_3 - t}{t_3 - t_2} \negspace && \negspace \frac{t - t_2}{t_3 - t_2} \negspace & \\ \boldsymbol{x}_0 &&&& \boldsymbol{x}_1 &&&& \boldsymbol{x}_2 &&&& \boldsymbol{x}_3 \end{array} \end{equation*}

The polynomial \(\boldsymbol{p}_{0,1,2,3}(t)\) at the apex can be evaluated for \(t_0 \le t \le t_3\). For a detailed explanation of this triangular scheme, see the notebook about the Barry–Goldman algorithm. Neville’s algorithm can be implemented for arbitrary degree:

def neville(xs, ts, t):
    """Lagrange interpolation using Neville's algorithm.

    Returns the interpolated value(s) at time(s) *t*,
    given the values *xs* at times *ts*.

    """
    if len(xs) != len(ts):
        raise ValueError('xs and ts must have the same length')
    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]

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 little class:

class Lagrange:

    def __init__(self, vertices, grid):
        assert len(vertices) == len(grid)
        self.vertices = vertices
        self.grid = grid

    def evaluate(self, t):
        return neville(self.vertices, self.grid, t)

Since this class has the same interface as the splines that will be discussed in the following sections, we can use a spline helper function from helper.py for plotting:

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:

l1 = Lagrange([(2, -2), (-1, 0), (0.3, 0.5), (3.14, -1), (1, -1)], ts)

plot_spline_2d(l1)

Runge’s Phenomenon

This seems to work to some degree, 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.

vertices1 = [
    (-2, 3),
    (1, 1),
    (3, -1),
    (2, -1),
    (2.5, 1.5),
    (5, 2),
    (6, 1),
    (5, 0),
    (6.5, -1),
    (7, 0),
    (6, 3),
]

l2 = Lagrange(vertices1, range(len(vertices1)))
plot_spline_2d(l2)

Here we see a severe overshooting effect, most pronounced at the beginning and the end of the curve. Moving some vertices can make this even worse. This effect is called Runge’s phenomenon. A possible mitigation for this overshooting is to use so-called Chebyshev nodes as time instances:

def chebyshev_nodes(a, b, n):
    k = np.arange(n) + 1
    nodes = np.cos(np.pi * (2 * k - 1) / (2 * n))
    return (a + b) / 2 - (b - a) * nodes / 2

l3 = Lagrange(vertices1, chebyshev_nodes(0, len(vertices1) - 1, len(vertices1)))
plot_spline_2d(l3)

This is definitely better. But it gets worse again when we move a few of the vertices.

vertices2 = [
    (0, -1),
    (1, 1),
    (3, -1),
    (2.5, 1.5),
    (5, 2),
    (6, 0.5),
    (6, 0),
    (4, -1),
    (6.5, -1),
    (7, 2),
    (8, 0),
]

l4 = Lagrange(vertices2, chebyshev_nodes(0, len(vertices2) - 1, len(vertices2)))
plot_spline_2d(l4)

Long story short, Lagrange interpolation is typically not suitable for drawing curves. For comparison, and as a teaser for the following sections, let’s use the same vertices to create a uniform Catmull–Rom spline:

import splines

cr_spline = splines.CatmullRom(vertices2)

plot_spline_2d(cr_spline)

And to get an even better fit, we can try a centripetal Catmull–Rom spline:

cr_centripetal_spline = splines.CatmullRom(vertices2, alpha=0.5)

plot_spline_2d(cr_centripetal_spline)

Note

The class splines.CatmullRom uses “natural” end conditions by default.