Cubic spline interpolation | Hacker News Reader

Cubic spline interpolation

eli.thegreenplace.net

92 pointssigna1128 comments

> 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

> in general perform much better than cubic splines

Source? (I couldn't find that in the paper you linked.)

The paper gives the rate of convergence you get with Polynomial interpolants in Chebyshev nodes:

> If f has v derivatives, with the vth derivative being of bounded variation V, then ||f - p_n|| = O(V n^{-v}) as n -> ∞

and

> If f is analytic, the convergence is geometric, with ||f - p_n|| = O(p^{-n}) for some p > 1

You will not get that good of a rate of convergence with cubic splines. See https://www.researchgate.net/publication/243095286_On_the_Or...

This is further explained in Trefethen's book https://www.amazon.com/Approximation-Theory-Practice-Applied...

Quoting from Ch 14

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

You can see also the software package https://www.chebfun.org/ for Chebyshev interpolations with Matlab and https://github.com/rnburn/bbai for Chebyshev interpolation of arbitrary dimension functions with sparse grids for Python. And here is a quick notebook for an experiment you can run that will compare Chebyshev interpolants to cubic splines: https://github.com/rnburn/bbai/blob/master/example/13-sparse...

Thank you!

However, in many circumstances, you have equispaced nodes (not Chebyshev nodes), for which global polynomial interpolation is terrible and it is literally impossible to make a global method which is without pitfalls. Trefethen has some other papers about this, esp. http://people.maths.ox.ac.uk/~trefethen/impossibility.pdf, but also a new paper in 2023 with a different recommended method ("AAA") https://link.springer.com/article/10.1007/s10543-023-00959-x

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.

That's a great paper, thanks for the pointer.

> This is a misconception that's often repeated

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

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.

https://splines.readthedocs.io/en/latest/euclidean/natural-u...

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.

Here's a write-up I made a couple years ago which I think does a bit better describing the math in question:

https://observablehq.com/@jrus/cubic-spline

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

[1] https://www.semanticscholar.org/paper/Finite-element-methods...

[2] https://en.wikipedia.org/wiki/Flat_spline

[3] https://pages.cs.wisc.edu/~deboor/draftspline.html

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.

[1]: https://www.youtube.com/watch?v=k-yUdqRXijo

That post should come with a "Danger: Rabbit Hole" warning. What a great lecturer!

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.

Best resource on this topic:

https://youtu.be/jvPPXbo87ds

https://youtu.be/aVwxzDHniEw

Freya has great videos on other maths subjects, and makes them quite approachable! Definitely recommend checking her channel

These videos are about a completely different topic than the linked post.

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.

[0] https://archive.org/details/creativecomputing-1980-07/page/n...

[1] https://en.wikipedia.org/wiki/Monotone_cubic_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.

Would split-tweak (chaikins, jareks...) be applicable here? It's easy.

(It's weirdly difficult to find an example to link)

Here's chaikins

https://sighack.com/post/chaikin-curves

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!

Reply on news.ycombinator.com