> There are many techniques to interpolate between a given set of points. Polynomial interpolation can perfectly fit N points with an N-1 degree polynomial, but this approach can be problematic for large a N; high-degree polynomials tend to overfit their data, and suffer from other numerical issues like Runge's phenomenon.
This is a misconception that's often repeated. High degree polynomial interpolations are problematic if you use equispaced points. If you use Chebyshev points, they are highly accurate and in general perform much better than cubic splines. See myth 1 from Lloyd Trefethen's paper: https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf
I don’t see how “If you use Chebyshev points, polynomial interpolation isn’t problematic” is a refutation of the claim “Polynomial interpolation between a given set of points can be problematic”. That given set isn’t necessarily, and typically isn’t, a set of Chebyshev points.
That’s like saying “Integer factorization isn’t hard. If you pick powers of ten, it’s easy”.
> In fact, polynomial interpolants in Chebyshev points are problem-free when evaluated by the barycentric interpolation formula. They have the same behavior as discrete Fourier series for period functions, whose reliability nobody worries about. The introduction of splines is a red herring: the true advantage of splines is not that they converge where polynomials fail to do so, but that they are more easility adapted to irregular point distributions and more localized.
In some cases cubic splines are a convenient and good enough tool, computationally cheap because the tridiagonal system involved can be solved in linear time.
You don't need to use O(n^3) Gaussian elimination to find natural cubic splines: you can write the system of equations as a tridiagonal matrix, which can be solved in linear time.
You don’t need to solve a system of equations at all, and it’s unusual to do so unless you have very specific requirements, right? All of the reasons the article used to justify going that direction can be met (more or less, depending on more nuance than the article used) with existing spline types. Plus you can trade some of the overshoot properties of the article’s solution for just a little bit of extra curvature, with a Catmull-Rom spline for example, which might be preferable for many people for multiple reasons.
Spline interpolation is a rich and extensively studied field, dating back to at least the 1940s. Here are a few comments about related work. The interpolation algorithm itself can be adjusted in numerous ways. One effective algorithm for reducing overshooting artifacts is Akima spline interpolation, but there are many others. Cubic splines have noteworthy theoretical properties; you can search for "energy minimizer cubic spline" to find out more. As already mentioned in @creata's comment, instead of using Gaussian elimination, you can take advantage of the fact that the matrix is banded and solve it using a banded solve algorithm (search for DGBSV in LAPACK). There is also research that explores the relationship between shift-invariant spaces and functions that can be expressed as a linear combination of splines. Finally, multivariate interpolation is an extension of the 1-d theory to n-d. This is used in image and surface interpolation, for example.
I was hoping someone would mention the variational properties of splines, so I'm glad to see you bring up their energy minimization. In Strang's Intro to Applied Mathematics, he briefly mentions (at the end of the section on splines) that a spline interpolation can be found by minimum principles, since it's essentially passing a beam in bending through the control points. I was curious if anyone has used that practically, eg writing a finite element spline solver?
I am unaware of any effort in that direction, but I did find a 1987 paper by Höllig [1] that seems to do a kind of converse (approximates FEM solutions using splines. I haven't read the paper, but that's my quick take). Recovering splines using a FEM might be an interesting learning exercise, since the analytic solution is known and the model seems simple. The motivating physical model is described by flat splines [2], and Carl de Boor has photographs of physical splines in action [3].
Polynomials, sampling theory, signal processing, Fourier series, etc. are all so closely related that learning even a little bit of one of these things opens up such a wealth of other things to think about and explore. Splines like this basically answer the question "what polynomial fits these points the best?" which, it turns out, is the same question you need to pose to begin formulating gaussian quadrature. And gaussian quadrature is super cool because it approximates integrals about a billion times better than riemann summation[1] for the same amount of computational effort.
I'm taking a numerical algorithm class and have been learning about gaussian quadrature for the first time. It's super cool stuff. I haven't been able to wrap my head around why it works entirely and have found readings about legendre polynomials really difficult to understand. Do you have any intuition you can share?
It's by no means rigorous, but the gist of it is that
1. You're integrating the polynomial that best approximates the function over the interval. This is why quadrature is exact for polynomials up to some order: the best approximating polynomial is the polynomial itself.
2. Polynomials are good approximations when the function is smooth. Most useful functions are.
3. Because of the smoothness, the behavior in the middle of the interval is largely affected by the behavior at the edges, so you need more densely sample the edges. Think of it like wiggling a string. You're only allowed to wiggle the end, but that very specifically defines the behavior in the middle.
4. There are lots of polynomial sets - Legendre, Chebyshev, Hermite, etc. They're each useful because they're orthogonal to a special weight function, and Legendre polynomials are kind of the default set because they have the simplest weight function: 1.
Did you know you can do spline interpolation on a piece of paper with the ruler and pencil? It's a nice demonstration that you can do to high school students, and maybe even to middle schoolers.
Here's an example, with fairly round numbers, because we need to use text instead of paper.
Let's say you have a function f, and its values at 1,2 and 3. f(1) = 1, f(2) = 3, f(3) = 2. We plot these 3 points, we call them A,B,C. What is the spline interpolation value at 1.5 ?
You start with the linear interpolation between A (the point (1,1)) and B (the point (2,3)). You get the value 2. Let's call the point (1.5, 2) as M.
Then you do the linear extrapolation of the segment BC. The new point is N with coordinates (1.5, 3.5).
Now move M horizontally until it hits the vertical line y=1, call that point M' (it has coordinates (1,2)). Move N to the vertical line y=3, it has coordinates (3,3.5). Draw the line M'N'. Where it intersects the vertical line y=1.5, that's where the spline interpolation is. It is the point (1.5, 2.375).
Of course, this is not only useful for paper and pencil stuff. That's how splines are actually implemented in lots of packages. You can easily implement it in Excel, if you need to, all with only cell formulas.
Stineman interpolation [0] is an alternative that is quick and "dirty" imputer for missing values. Of course there are monotonic-preserving versions [1] of splines that are excellent as well.
For my uses, cubic interpolation is simply not of value when I have noisy points and/or unpredictable behaviour between points --- Kalman smoothers or Gaussian process/Kriging gives me both a good mean estimate between points and a sense of the error associated with the interpolation.
I second this. In quantitative finance cubic splines are often, simply, out of the question. However, we end up using montonic cubic splines quite often.
Somewhat related: Something that helped me understand splines more intuitively was to see them animated in a way that they can be played with. https://www.jasondavies.com/animated-bezier/
That's a cool page but it's about Bézier curves, while this article is about cubic splines. They are not the same!
The most notable difference is that a Bézier curve is defined in terms of control points that do not (generally) pass through the curve, while a polynomial spline is defined by the points it passes through; there are no additional control points.
Coincidentally, I recently put up a profile coffee roasting planner built on cubic splines. The constraints of the cubic spline match up nicely with the behavior of coffee in a roasting machine so doing it this way gets you plans that are easy to follow and fast to put together or edit. This one is just done as a web site, but prior to that I'd been using a prototype in C++ (unreleased because no UI) for several years. Try not to look at the code too closely, not a web dev. https://crucs.net.
The Chaikin curve is a quadratic B-spline, which is an approximating spline rather than an interpolating spline. Chaikin is very useful for lots of applications, but the reason it doesn’t fit the stated constraints in the article is because it doesn’t interpolate the control points, and doesn’t “naturally” stop at the endpoints of the curve either (though it’s obviously very easy to make your Chaikin curve hit endpoints if you want to). Playing with interpolating quadratic splines is a fun exercise!
This is a misconception that's often repeated. High degree polynomial interpolations are problematic if you use equispaced points. If you use Chebyshev points, they are highly accurate and in general perform much better than cubic splines. See myth 1 from Lloyd Trefethen's paper: https://people.maths.ox.ac.uk/trefethen/mythspaper.pdf