Skip to content

Polynomials from repeated linear interpolation

August 5, 2012

It’s fairly well-known that Bezier curves can be evaluated using repeated linear interpolation – de Casteljau’s algorithm. It’s also fairly well-known that a generalization of this algorithm can be used to evaluate B-Splines: de Boor’s algorithm. What’s not as well known is that it’s easy to construct interpolating polynomials using a very similar approach, leading to an algorithm that is, in a sense, halfway between the two.

In the following, I’ll write l(\alpha,x,y) := (1 - \alpha)x + \alpha y for linear interpolation. I’ll stick with quadratic curves since they are the lowest-order curves to show “interesting” behavior for the purposes of this article; everything generalizes to higher degrees in the obvious way.

De Casteljau’s algorithm

De Casteljau’s algorithm is a well-known algorithm to evaluate Bezier Curves. There’s plenty of material on this elsewhere, so as usual, I’ll keep it brief. Assume we have three control points x_0, x_1, x_2. In the first stage, we construct three constant (degree-0) polynomials for the three control points:

p_0^{[0]}(t) = x_0
p_1^{[0]}(t) = x_1
p_2^{[0]}(t) = x_2

These are then linearly interpolated to yield two linear (degree-1) polynomials:

p_0^{[1]}(t) = l(t, p_0^{[0]}(t), p_1^{[0]}(t))
p_1^{[1]}(t) = l(t, p_1^{[0]}(t), p_2^{[0]}(t))

which we then interpolate linearly again to give the final result

p_0^{[2]}(t) = l(t, p_0^{[1]}(t), p_1^{[1]}(t)) = B(t).

Note I give the construction of the full polynomials here; the actual de Casteljau algorithm gets rid of them immediately by evaluating all of them as soon as they appear (only ever doing linear interpolations). Anyway, the general construction rule we’ve been following is this:

p_i^{[d]}(t) = l(t, p_i^{[d-1]}(t), p_{i+1}^{[d-1]}(t))

De Boor’s algorithm

De Boor’s algorithm is the equivalent to de Casteljau’s algorithm for B-Splines. Again, there’s plenty of material out there on it, so I’ll keep it brief: We again start with constant functions for our data points. This time, the exact formulas depend on the degree of the spline we’ll be using. I’ll be using the degree k=2 (quadratic) here. We’ll also need a knot vector (t_i) which determines where our knots are; knots are (very roughly) the t’s corresponding to the control points. I’ll be using slightly different indexing from what’s normally used to make the similarities more visible, and ignore issues such as picking the right set of control points to interpolate from:

p_0^{[0]}(t) = x_0
p_1^{[0]}(t) = x_1
p_2^{[0]}(t) = x_2

Then we linearly interpolate based on the knot vector:

p_0^{[1]}(t) = l((t - t_0) / (t_2 - t_0), p_0^{[0]}(t), p_1^{[0]}(t))
p_1^{[1]}(t) = l((t - t_1) / (t_3 - t_1), p_1^{[0]}(t), p_2^{[0]}(t))

and interpolate again one more time to get the results:

p_0^{[2]}(t) = l((t - t_1) / (t_2 - t_1), p_0^{[1]}(t), p_1^{[1]}(t))

The general recursion formula for de Boor’s algorithm (with this indexing convention, which is non-standard, so do not use this for reference!) is this:

p_i^{[d]}(t) = l((t - t_{i+d-1}) / (t_{i+k} - t_{i+d-1}), p_i^{[d-1]}(t), p_{i+1}^{[d-1]}(t))

Interpolating polynomials from linear interpolation

There’s multiple constructions for interpolating polynomials; the best-known are probably the Lagrange polynomials (which form a basis for the interpolating polynomials of degree n for a given set of nodes t_i) and the Newton polynomials (since polynomial interpolation has unique solutions, these give the same results, but the Newton formulation is more suitable for incremental evaluation).

What’s less well known is that interpolating polynomials also obey a simple triangular scheme based on repeated linear interpolation: Again, we start the same way with constant polynomials

p_0^{[0]}(t) = x_0
p_1^{[0]}(t) = x_1
p_2^{[0]}(t) = x_2

and this time we have associated nodes t_0, t_1, t_2 and want to find the interpolating polynomial I(t) such that I(t_0)=x_0, I(t_1)=x_1, I(t_2)=x_2. Same as before, we first try to find linear functions that solve part of the problem. A reasonable choice is:

p_0^{[1]}(t) = l((t - t_0) / (t_1 - t_0), p_0^{[0]}(t), p_1^{[0]}(t))
p_1^{[1]}(t) = l((t - t_1) / (t_2 - t_1), p_1^{[0]}(t), p_2^{[0]}(t))

Note the construction here. p_0^{[1]} is a linear polynomial that interpolates the data points (t_0,x_0), (t_1,x_1), and we get it by interpolating between two simpler (degree-0) polynomials that interpolate only (t_0,x_0) and (t_1,x_1), respectively: we simply make sure that at t=t_0, we use p_0^{[0]}, and at t=t_1, we use p_1^{[0]}. All of this is easiest to visualize when t_0 \le t_1 \le t_2, but it in facts works with them in any order. p_1^{[1]} is constructed the same way.

To construct our final interpolating polynomial, we use the same trick again:

p_0^{[2]}(t) = l((t - t_0) / (t_2 - t_0), p_0^{[1]}(t), p_1^{[1]}(t)) = I(t).

Note this one is a bit subtle. We linearly interpolate between two polynomials that both in turn interpolate (t_1,x_1); this means we already know that the result will also pass through this point. So t_1 is taken care of, and we only need to worry about t_0 and t_2 – and for each of the two, one of our two polynomials does the job, so we can do the linear interpolation trick again. The generalization of this approach to higher degrees requires that we make sure that both of our input polynomials at every step interpolate all of the middle points, so we only need to fix up the ends. But this is easy to arrange – the general pattern should be clear from the construction above. This gives us our recursive construction rule:

p_i^{[d]}(t) = l((t - t_i) / (t_{i+d} - t_i), p_i^{[d-1]}(t), p_{i+1}^{[d-1]}(t))

All of this is, of course, not new; in fact, this is just Neville’s algorithm. But in typical presentations, this is derived purely algebraically from the properties of Newton interpolation and divided differences, and it’s not pointed out that the linear combination in the recurrence is, in fact, a linear interpolation – which at least to me makes everything much easier to visualize.

The punchline

The really interesting bit to me though is that, starting from the exact same initial conditions, we get three different important interpolation / approximation algorithms that differ only in how they choose their interpolation factors:

de Casteljau: p_i^{[d]}(t) = l(t, p_i^{[d-1]}(t), p_{i+1}^{[d-1]}(t))
Neville: p_i^{[d]}(t) = l((t - t_i) / (t_{i+d} - t_i), p_i^{[d-1]}(t), p_{i+1}^{[d-1]}(t))
de Boor: p_i^{[d]}(t) = l((t - t_{i+d-1}) / (t_{i+k} - t_{i+d-1}), p_i^{[d-1]}(t), p_{i+1}^{[d-1]}(t))

I think this quite pretty. B-Splines with the right knot vector (e.g. [0,0,0,1,1,1] for the quadratic curves we’ve been using) are just Bezier Curves, that bit is well known. But what’s less well known is that Neville’s Algorithm (and hence regular polynomial interpolation) is just another triangular linear interpolation scheme that fits inbetween the two.

From → Maths

4 Comments
  1. Robin Green permalink

    OK, so convert each method to weighted sums of basis functions. Basis functions can be converted between each other, and Fourier analysis is just projecting functions to a sinusoidal basis so there’s Spectral Analysis, which leads into Generalized Sampling and Shannon’s Information theory. Basis functions also lead to Wavelets and, if you relax the self-similarity preconditions, you’re into the world of Frames. It’s all part of the same thing!

    Reply
    • fgiesen permalink

      The connection from anything polynomial to Fourier Analysis is far more direct though: Just go from real polynomials p(x) to the complex domain p(z) and substitute z=exp(-it). Presto trigonometric polynomials (which are really just finite Fourier series).

      Reply
      • Per Vognsen permalink

        You didn’t quite spell it out. But the real punchline is that the substitution establishes a relationship between complex power series on the unit circle and Fourier series. With only a little work, it gives you an extremely weak but still practically useful part of the convergence theory for Fourier series. Namely, if the complex-valued function defined on the unit circle has an analytic continuation to the entire unit disk, then we know from complex analysis 101 that this extended function has a convergent power series expansion around the origin, and hence the original function on the unit circle has a convergent Fourier expansion.

  2. Per Vognsen permalink

    By the way, another way you see this relationship between complex power series and Fourier series is in Cauchy’s integral formula for calculating derivatives (and hence coefficients of power series). It integrates around a circle and it should look really familiar–it’s completely analogous to the formula for calculating Fourier coefficients where you integrate over a sinusoidal period.

    Reply

Leave a comment

«
»