Polynomial Interpolation methods/algorithms

[juha_p]

Been testing several interpolation algorithms (actually >60 already):

Code: Select all

Function: Sin(x) 
Range: [-pi/4;pi/4]

4-point interpolation (equally spaced points) of function Sin(Sqrt(x))/Sqrt(x)
Test: fast test with 10000 samples (in Range)

AnalyticInt	: rel.err:~  7.7601e-06  @ pos:~ -7.2360e-01 - Piecewise analytic interpolation
BaryInv		: rel.err:~  1.0390e-02  @ pos:~ -7.7952e-01 - Barycentric interpolation with inverse distance weighting
Barylag		: rel.err:~  1.2966e-08  @ pos:~ -3.1796e-01 - Barycentric lagrange interpolation
BulirschStoer	: rel.err:~  3.0480e-08  @ pos:~  3.1698e-01 - Rational Bulirsch-Stoer interpolation
Cakima		: rel.err:~  9.3585e-07  @ pos:~  2.7874e-01 - Piecewise cubic Akima interpolation
CBezier		: rel.err:~  6.1144e-05  @ pos:~  3.2204e-01 - Piecewise cubic Bezier spline interpolation
Chermite	: rel.err:~  1.3788e-06  @ pos:~ -3.3441e-01 - Piecewise cubic Hermite spline interpolation (Finite difference)
Chermite	: rel.err:~  8.5600e-05  @ pos:~ -2.9732e-01 - Piecewise cubic Hermite spline interpolation (Catmull-Rom spline)
Chermite	: rel.err:~  8.5600e-05  @ pos:~ -2.9732e-01 - Piecewise cubic Hermite spline interpolation (Monotone interpolation)
Chermite	: rel.err:~  6.6922e-07  @ pos:~ -7.8529e-01 - Piecewise cubic Hermite spline interpolation (Monotone with Lam harmonic mean)
ConstInt	: rel.err:~  2.2558e-02  @ pos:~  6.2665e-01 - Piecewise constant interpolation
CosInt		: rel.err:~  4.8498e-03  @ pos:~  7.6286e-01 - Piecewise cosine interpolation
CSPline		: rel.err:~  6.1144e-05  @ pos:~  3.2204e-01 - Piecewise cubic spline interpolation
CubicFarrow	: rel.err:~  9.4201e-02  @ pos:~  5.1166e-01 - Cubic Farrow interpolation
CubiConv	: rel.err:~  1.3788e-06  @ pos:~ -3.3441e-01 - Cubic convolution interpolation
CubInt		: rel.err:~  1.1613e+01  @ pos:~ -7.8540e-01 - Piecewise cubic interpolation (Finite difference)
CubInt		: rel.err:~  1.1498e+01  @ pos:~ -7.8540e-01 - Piecewise cubic interpolation (Catmull-Rom spline)
CubInt		: rel.err:~  1.1498e+01  @ pos:~ -7.8540e-01 - Piecewise cubic interpolation (Monotone interpolation)
CubInt		: rel.err:~  1.1613e+01  @ pos:~ -7.8540e-01 - Piecewise cubic interpolation (Monotone with Lam harmonic mean)
ExpInt		: rel.err:~ -9.9120e-05  @ pos:~  7.8539e-01 - Piecewise exponential interpolation
FloaterHormann	: rel.err:~  1.2966e-08  @ pos:~  3.1795e-01 - Rational interpolation using the Floater-Hormann Method
FractInt  	: - self-affine fractal interpolation for specific interpolation points
HermInt	  	: - piecewise Hermite interpolation
HilbInt   	: rel.err:~ -7.3502e-01  @ pos:~  7.8539e-01 - Piecewise Hilbert transform interpolation (Rectangular window)
HilbInt		: rel.err:~ -9.9586e-01  @ pos:~  7.8539e-01 - Piecewise Hilbert transform interpolation (Lanczos window)
HilbInt		: rel.err:~ -9.5708e-01  @ pos:~  7.8539e-01 - Piecewise Hilbert transform interpolation (Hamming window)
HilbInt		: rel.err:~ -9.6904e-01  @ pos:~  7.8539e-01 - Piecewise Hilbert transform interpolation (Squared cosine window)
HilbInt		: rel.err:~ -8.9266e-01  @ pos:~  7.8539e-01 - Piecewise Hilbert transform interpolation (Kaiser window (alpha = 3))
HilbInt		: rel.err:~ -9.0682e-01  @ pos:~  7.8539e-01 - Piecewise Hilbert transform interpolation (Blackman window)
Krogh     	: rel.err:~ -1.0313e-08  @ pos:~  7.8530e-01 - Interpolation using Krogh algorithm
LagInt	  	: rel.err:~ -1.1928e-06  @ pos:~  7.8539e-01 - Piecewise Lagrange interpolation
LinInt    	: rel.err:~  1.4836e-04  @ pos:~  6.2766e-01 - Piecewise linear interpolation
MinimaxPoly	: - Optimal polynomial via the Remez algorithm
MinimaxRat	: - Optimal rational polynomial via the Remez algorithm
MQSpline	: rel.err:~  6.6922e-07  @ pos:~ -7.8529e-01 - Piecewise monotone quadratic spline interpolation (Lam harmonic mean)
MQSpline	: rel.err:~  5.8454e-05  @ pos:~ -3.2521e-01 - Piecewise monotone quadratic spline interpolation (Schumaker estimation)
Neville   	: rel.err:~  1.2966e-08  @ pos:~ -3.1795e-01 - Interpolation using Neville's Method
NewtInt		: rel.err:~  1.2966e-08  @ pos:~  3.1795e-01 - Interpolation of equally-spaced points (Forwards method)
NewtInt		: rel.err:~  1.2969e-08  @ pos:~ -3.1795e-01 - Interpolation of equally-spaced points (Backwards method)
NewtInt		: rel.err:~  1.1149e-02  @ pos:~ -3.6180e-01 - Interpolation of equally-spaced points (Stiring's method)
NewtInt		: rel.err:~  7.4651e-05  @ pos:~  3.3456e-01 - Interpolation of equally-spaced points (Everett's method)
PolyCof   	: rel.err:~  5.0620e-05  @ pos:~ -7.8540e-01 - Non-Vandermonde polynomial interpolation coefficients
QHermite	: rel.err:~  1.4051e-06  @ pos:~  3.2608e-01 - Piecewise quintic Hermite interpolation (Finite difference)
QHermite  	: rel.err:~  1.0433e-04  @ pos:~  2.9340e-01 - Piecewise quintic Hermite interpolation (Catmull-Rom spline)
QSpline   	: rel.err:~  1.8161e-02  @ pos:~  6.2782e-01 - Piecewise quadratic spline interpolation
RatInt    	: rel.err:~ -8.1855e-07  @ pos:~  7.8539e-01 - Simple rational polynomial interpolation
RCHermite 	: rel.err:~ -1.1989e-04  @ pos:~  7.8539e-01 - Piecewise rational cubic Hermite spline interpolation
Said		: rel.err:~  1.9109e-01  @ pos:~  3.4258e-01 - Piecewise Said interpolation (Sinc)
Said		: rel.err:~  1.0942e-01  @ pos:~  3.0902e-01 - Piecewise Said interpolation (Lanczos, M=2 )
Said		: rel.err:~  1.2912e-01  @ pos:~  3.1993e-01 - Piecewise Said interpolation (Lanczos, M=3 )
Said 		: rel.err:~  1.4596e-01  @ pos:~ -3.3161e-01 - Piecewise Said interpolation (Lanczos, M=4 )
Said		: rel.err:~  1.5936e-01  @ pos:~ -3.3681e-01 - Piecewise Said interpolation (Lanczos, M=5 )
Said		: rel.err:~  9.2656e-02  @ pos:~ -2.9740e-01 - Piecewise Said interpolation (Blackman-Harris, N=6)
Said		: rel.err:~  1.1773e-01  @ pos:~ -3.1600e-01 - Piecewise Said interpolation (Cubic B-Spline)
Said		: rel.err:~  5.1442e-02  @ pos:~  2.5192e-01 - Piecewise Said interpolation (Mitchell-Netravali, B=C=1/3)
Schwerner 	: rel.err:~  3.0480e-08  @ pos:~ -3.1698e-01 - Rational interpolation using the Schneider-Werner Method
SHermite  	: rel.err:~  1.4296e-06  @ pos:~  3.2123e-01 - Piecewise septic Hermite interpolation
SincDInt  	: rel.err:~  1.9484e-02  @ pos:~ -3.1358e-01 - Piecewise discrete sinc interpolation
SincInt   	: rel.err:~  8.7148e-02  @ pos:~ -2.9989e-01 - Piecewise sinc interpolation
Steffen   	: rel.err:~  1.3788e-06  @ pos:~ -3.3441e-01 - Monotonic Steffen interpolation
Stineman  	: rel.err:~  1.2843e-06  @ pos:~ -3.3236e-01 - Monotonic Stineman interpolation
Thiele    	: rel.err:~ -7.3597e-09  @ pos:~  7.8539e-01 - Rational Thiele interpolation
TrigInt   	: rel.err:~ -1.6538e-04  @ pos:~  7.8539e-01 - Piecewise trigonometric interpolation

Tested methods comes from Joe Henning (2022). 
Interpolation Utilities (https://www.mathworks.com/matlabcentral/fileexchange/36800-interpolation-utilities), MATLAB Central File Exchange. Retrieved May 9, 2022. 
Krogh's algorithm is from scipy library which can be fround from https://programtalk.com/vs2/?source=python/11627/scipy/scipy/

but, I have not found usable source code for Krogh's interpolation algorithms, which ones, by many sources, might be a bit faster than those commonly mentioned.

I've found few papers as like these:
Krogh: Efficient Algorithms for Polynomial Interpolation (Semantic Scholar)
Macleod: A comparison of algorithms for polynomial interpolation (Science Direct):

but, as those mathematical equations there are out of my math skills, these are not very helpful for me.

KroghInterpolator (Python implementation) I've bumped to can be found in scipy library (Program Talk) but, I'm not familiar with python language ... and the implementation looks (by my skills) quite complex to implement in C/C++. Managed to found out what's happening there in main-loop but, (A*) there are tasks like initalization, reshaping arrays and axis rolling which I did not got solved and also, how/what those c and Vk arrays acts/holds in in calculations (two dimensional in definition, one dimensional in calculations).

<code snipped deleted>

So, in short, looking a quick path

Q: Which one of Krogh's algorithms they used in scipy implementation?
Q: Any suggestions for A*
Q: Any open source (C/C++/Matlab) Krogh algorithm implementations available?
---

BTW, It looks like (at least) GCC and Clang can handle well some of those calculations involved, as like it does as for an example, in case of Lagrange interpolation:

Desmos 3-point 4-point
Compiler Explorer -
Quick C++ benchmark Clang GCC

[2DaT]

I never had a problem with newton interpolating polynomial.

https://en.wikipedia.org/wiki/Newton_polynomial

For small n (n<10) the source of error is usually the bad choice of nodes. It's well known that Chebyshev nodes are good a priori choice for polynomial interpolation.

https://en.wikipedia.org/wiki/Chebyshev_nodes

These nodes are great for polynomial approximation if you don't want to bother with remez algorithm. Though if you need to optimize for relative error, I don't know good one step method (and I don't need to, because remez exists).

I kinda have a "hack" for generating somewhat good polynomials that interpolate in relative sense. Instead of building an interpolating polynomials on N cheb nodes, build weighted equioscillaton polynomial on N+1 in the form of P(x) +(-1)^n*f(x)*E. This is first step of the remez algorithm and you can find how to construct it here:
https://en.wikipedia.org/wiki/Remez_algorithm

The results are rather nice: https://www.desmos.com/calculator/jc42z9iblm

But if you care that much, might aswell do full Remez algorithm. It's rather simple, especially on grid. Instead of treating the interval as continuity, use like 100k points in it. This doesn't give you "true" minmax, but it's very close for well behaved functions.


If you want to approximate something on the fly then high order polynomial interpolation is usually a bad choice. If it's something graphics related obvious choice would be splines, if it's audio it would be some type of sinc interpolation.

[mystran]

What are we trying to do?

There is one unique interpolating (=passes through points) polynomial of minimum degree (=n-1) for a given set of n points. This can be computed in the Newton form, Lagrange form (straight or barycentric) and probably a bunch of other ways, but all of these lead to the same polynomial.

That said, polynomial interpolation typically doesn't make a whole lot of sense, except as a building block for either polynomial approximation or spline interpolation (and really, even for splines consider Hermite-splines instead).

[juha_p]

What are we trying to do?
...

Well, I have for years used LibreOffice Calc (Trendline -tool) for to get polynomials for to represent certain functions / curves. It has working otherwise well but, at least in Ubuntu, Calc tends to freeze the computer when doing certain tasks and also one have to copy equations from Calc chart to LibreOffice Draw program, break them there, copy then to a ASCII editor, make changes (replace minus signs with correct one, change powers representation format, etc.) before polynomials are usable in programming IDEs or places like Desmos, etc.) This has been cumbersome task to do especially when freezing happens and computer needs reboot ... . .
By learning/implementing various interpolation methods/algorithms, this task become simpler.

...

There is one unique interpolating (=passes through points) polynomial of minimum degree (=n-1) for a given set of n points. This can be computed in the Newton form, Lagrange form (straight or barycentric) and probably a bunch of other ways, but all of these lead to the same polynomial.

...

Maybe not exactly same polynomials.

Here are results for few (most accurate from my testing) 4-point interpolations of function sin(sqrt(x))/sqrt(x) (shape is almost straight line):

Code: Select all

Barylag		: rel.err:~  1.2966e-08  @ pos:~ -3.1796e-01 - Barycentric lagrange interpolation
BulirschStoer	: rel.err:~  3.0480e-08  @ pos:~  3.1698e-01 - Rational Bulirsch-Stoer interpolation
FloaterHormann	: rel.err:~  1.2966e-08  @ pos:~  3.1795e-01 - Rational interpolation using the Floater-Hormann Method
Krogh     	: rel.err:~ -1.0313e-08  @ pos:~  7.8530e-01 - Interpolation using Krogh algorithm
LagInt	  	: rel.err:~ -1.1928e-06  @ pos:~  7.8539e-01 - Piecewise Lagrange interpolation
Neville   	: rel.err:~  1.2966e-08  @ pos:~ -3.1795e-01 - Interpolation using Neville's Method
NewtInt		: rel.err:~  1.2966e-08  @ pos:~  3.1795e-01 - Interpolation of equally-spaced points (Forwards method)
NewtInt		: rel.err:~  1.2969e-08  @ pos:~ -3.1795e-01 - Interpolation of equally-spaced points (Backwards method)
Schwerner 	: rel.err:~  3.0480e-08  @ pos:~ -3.1698e-01 - Rational interpolation using the Schneider-Werner Method
Thiele    	: rel.err:~ -7.3597e-09  @ pos:~  7.8539e-01 - Rational Thiele interpolation (P=1)
Thiele		: rel.err:~ -7.3597e-09  @ pos:~  7.8539e-01 - Rational Thiele interpolation (P=2)

I have not tried to tweak anything nor do speed benchmark for these.

I kinda have a "hack" for generating somewhat good polynomials that interpolate in relative sense. Instead of building an interpolating polynomials on N cheb nodes, build weighted equioscillaton polynomial on N+1 in the form of P(x) +(-1)^n*f(x)*E. This is first step of the remez algorithm and you can find how to construct it here:
https://en.wikipedia.org/wiki/Remez_algorithm

The results are rather nice: https://www.desmos.com/calculator/jc42z9iblm

Yes, your results looks nice.

As a comparison, I made 4-point Lagrange, Neville and Thiele Interpolation methods (sin(pi/2*x) as in your example): https://www.desmos.com/calculator/axv5y9obr1
which shows similar error response (dunno why Octave reports much lower error level: Thiele : rel.err:~ 6.5770e-13 @ pos:~ 1.0000e+00 (I did not scan near zero values as Desmos plot shows that area being OK). By combining Thiele and Neville methods one can get even better interpolation accuracy (a6(x) in desmos plot).
I have not tried to find best possible positions for interpolation points and the equations can also be shortened quite a bit to get calculation faster (as for an example, speed for Thiele can be dropped from (rdtsc/val) 92.6 to 2.9 by combining and forming (tried Horner)).

[mystran]

...

There is one unique interpolating (=passes through points) polynomial of minimum degree (=n-1) for a given set of n points. This can be computed in the Newton form, Lagrange form (straight or barycentric) and probably a bunch of other ways, but all of these lead to the same polynomial.

...

Maybe not exactly same polynomials.

Given the conditions of passing a polynomial of degree N-1 through N points, that's literally all the degrees of freedom we have so there is exactly one such polynomial.

Note that rational functions will give different results, choosing a different set of points will give different results, doing piecewise interpolation (=splines) will give different results... but if you're just doing straight polynomial interpolation and you get different results from different algorithms, that's probably just numerical errors.

edit: In general, when playing with polynomials (or rational functions) it's worthwhile to try to understand what is going on to the point where you can mentally differentiate between choosing what polynomial (or rational function) you actually want vs. choosing an algorithm to compute that polynomial (or rational function).

Counting degrees of freedom generally works well: passing through a point takes one degree of freedom, forcing a specific derivative on a specific point takes one degree of freedom, forcing even or odd symmetry drops all odd or even order terms... and once you run out of degrees of freedom you'll have your polynomial (or rational function for a given choice of numerator/denominator degrees).

edit2: If you're trying to approximate a given function by fitting an interpolating polynomial through some chosen set of points, then the method you use to compute that polynomial absolutely does not matter except in terms of numerical errors. The only thing that really matters in this case is the choice of points. If you get to choose the points, then choosing something like Chebychev nodes is typically much better than choosing uniform points. If you're stuck with uniform points, consider splines (eg. Hermite, Catmull-Rom, something like that) instead.

[juha_p]

<accidentally quoted my own post>

[mystran]

Some of your interpolation methods say "rational" while many (all the polynomial ones?) have (effectively) the exact same errors. If you plot them, I predict those with exact same errors are exact same polynomials up to some numerical tolerances due to different computation methods.

Code: Select all

Barylag		: rel.err:~  1.2966e-08  @ pos:~ -3.1796e-01 - Barycentric lagrange interpolation
FloaterHormann	: rel.err:~  1.2966e-08  @ pos:~  3.1795e-01 - Rational interpolation using the Floater-Hormann Method
Neville   	: rel.err:~  1.2966e-08  @ pos:~ -3.1795e-01 - Interpolation using Neville's Method
NewtInt		: rel.err:~  1.2966e-08  @ pos:~  3.1795e-01 - Interpolation of equally-spaced points (Forwards method)
NewtInt		: rel.err:~  1.2969e-08  @ pos:~ -3.1795e-01 - Interpolation of equally-spaced points (Backwards method)

Look at the numbers from this reduced list. They are same, except Newton's backwards method probably running into some slight numerical precision issue. While FloaterHormann says "Rational" in the description, I'd guess based on the numbers that we probably have denominator degree zero here?

Now, if we consider Trendlines that's computing "regressions" (gotta love the statistics terminology) what it's probably doing is a least-squares fit. As it turns out, the problem of finding a polynomial can be written as a standard system of equations Ax=b where A is some NxN matrix. If x has N dimensions (polynomial of degree N-1) there is generally exactly one solution which you can find solving the system directly. If x has less than N dimensions (over-determined system), we can do an ordinary least-squares fit (eg. premultiply both sides by At the transpose of A, so that we end up with AtAx=Atb which can then be solved directly). Again, in the case where there is one exact fit (with polynomials; you can do regressions with all kinds of functions), you'll get the same polynomial as any other method.

[juha_p]

What interpolation method Mathematica uses in InterpolatingPolynomial()?

Noticed that it accepts non-integer n values... so I ran this code:

Code: Select all

<< FunctionApproximations`
SetOptions[SelectedNotebook[], PrintPrecision -> 12]

f[x_] = Sin[Sqrt[x]]/Sqrt[x];
a = -0.00000000001;
b = Pi/4;
n = 3.5; 
		   b - a            b - a
XY = N[Table[{a + ------- k, f[a + ------- k]}, {k, 0, n}]];
		      n		      n
P[x_] = N[Expand[InterpolatingPolynomial[XY, x]], 9];

HornerForm[N[Expand[x*P[x*x]], 9]]]

which resulted:

Code: Select all

P[x] = x (1. +x^2 (-0.166666482107+x^2 (0.00833182380375 -0.000194733764456 x^2)))

Here's comparison with n=3, n=3.5 and n= 3.9 against Remez method:
https://www.desmos.com/calculator/g45upysrf2

[2DaT]

Here's comparison with n=3, n=3.5 and n= 3.9 against Remez method:
https://www.desmos.com/calculator/8u7xlfmq0d

What's the point of using anything else if you have Remez available ? It's best by definition.

[mystran]

Here's comparison with n=3, n=3.5 and n= 3.9 against Remez method:
https://www.desmos.com/calculator/8u7xlfmq0d

What's the point of using anything else if you have Remez available ? It's best by definition.

Well, just like there is one unique polynomial (of minimum degree) that goes through a set of points, there is also one unique (up to symmetries in some cases?) polynomial (of a given degree) that minimizes the maximum error with respect to some continuous function to be approximated. There are several ways to find this "minimax" polynomial and one of these is indeed Remez which essentially combines the above two concepts and iteratively searches for the set of points that result in the minimax polynomial when you compute an interpolation polynomial.

One could probably argue whether the Remez algorithm specifically is always the "best" way to find the minimax poly (not that I'm trying to claim that Remez wasn't usually the algorithm that you want), but if you want to minimize the maximum error then the minimax polynomial is certainly the "best" polynomial by definition.

[juha_p]

Well, just like there is one unique polynomial (of minimum degree) that goes through a set of points, there is also one unique (up to symmetries in some cases?) polynomial (of a given degree) that minimizes the maximum error with respect to some continuous function to be approximated. There are several ways to find this "minimax" polynomial and one of these is indeed Remez which essentially combines the above two concepts and iteratively searches for the set of points that result in the minimax polynomial when you compute an interpolation polynomial.

One could probably argue whether the Remez algorithm specifically is always the "best" way to find the minimax poly (not that I'm trying to claim that Remez wasn't usually the algorithm that you want), but if you want to minimize the maximum error then the minimax polynomial is certainly the "best" polynomial by definition.

I did compare InterpolatingPolynomial, remez (Sollya) and MiniMaxApproximation (in this sine(x) case used as an example) and 1st mentioned resulted smallest error (~ 4.6e-9), Remez error was (~8.05e-9) and MiniMax had little more error.

[Borogove]

Take a look at Olli Niemitalo's "Elephant Paper" -- a lot of analysis on resampling interpolators with an eye toward computational efficiency tradeoffs.

[Kraku]

If by chance you are trying to interpolate audio data, then you want to check this out:
http://yehar.com/blog/?p=197

[mystran]

If by chance you are trying to interpolate audio data, then you want to check this out:
http://yehar.com/blog/?p=197

Wasn't it already linked in the post just before yours?

[Kraku]

If by chance you are trying to interpolate audio data, then you want to check this out:
http://yehar.com/blog/?p=197

Wasn't it already linked in the post just before yours?

I should learn to read the whole thread through before answering.