Mill Fighting - 1: Interpolation Splines
In this article, the lyrical hero challenges the optimal implementation of the classic polynomial Lagrange interpolator (Farrow), in the process of the battle accidentally opens and proves the trivial useless math magic spell with which he tries to press the opponent, but according to the results of all rounds of the battle, a draw is fixed by the judges .
Once a friend of mine asked me to simulate a simple and low-cost variant of resampling (changing the sampling frequency) of a digital signal. There were no special requirements for linear frequency response and phase response, and the desire to save time and memory was, therefore, some methods such as polyphase interpolating banks of FIR filters were immediately excluded from consideration, and the main attention was paid to local polynomial interpolation methods. Polynomials are quick and easy to read, storage of their coefficients requires little memory, and using a single set of coefficients, you can calculate any values within the original sampling interval.
It is known that through any N points (from here and further for definiteness we will consider a uniform grid - with the same step, although for many described techniques and methods this is not a prerequisite), so, through any N points with different x-coordinates, you can conduct a single polynomial of N-1 order uniquely determined by a set of N coefficients (read that LaTex on Habré is not directly, but through pictures that are heavy and short-lived, so I will only manage the upper-lower indices in html-tags):
')
f (t) = a n t n + a n-1 t n-1 + ..... + a 1 t + a 0
For convenience of calculations, we normalize our uniform argument grid to integer values, using the example of four points, choose a convenient interval of the argument (Farrow chose from [-2, 1] with central interval from [-1, 0]) and write the system of linear equations from the condition of passing the polynomial through our 4 points:
f (i) = y i
where i = {-2, -1, 0, 1}
You can decide as you like, even though Gauss can convert the matrix, but Farrow has optimized the number of operations to calculate the coefficients of this polynomial:
a 3 = (y 1 - y -2 ) / 6 + (y -1 - y 0 ) / 2
a 1 = (y 1 - y -1 ) / 2 - a 3
a 2 = y 1 - y 0 - a 1 - a 3
a 0 = y 0
Total operations:
This polynomial is used to calculate values that fall in the central interval of the original four points: [-1, 0] (with the corresponding rationing of the argument); for other intervals, the polynomials are calculated in the same way from the four points surrounding them. Very popular interpolator in narrow circles, deserved and revered. And for some reason I immediately wanted to come up with something that could surpass this implementation in a certain sense, taking into account the number of operations for the calculation.
You can see that both the sought polynom itself and all its derivatives are linear combinations of the unknown coefficients of the polynomial. Therefore, we can set the conditions for the equality of derivatives of any order at the nodes of the grid to the given values, and obtain a linear system of equations for the polynomial coefficients. Hermite's local cubic spline is based on four conditions: the values of the polynomial and its first derivative at the edges of the interval. And we can get the approximate values of these derivatives at the right points by differentiating the interpolating Lagrange polynomial passing through them. For example, we normalize the arguments of four points to the interval [-1, 2], take three points with arguments {-1, 0, 1} and draw a parabola through them, and calculate its coefficients. And since we only need the value of the first derivative of this parabola at the center point 0, it is enough to calculate its coefficient for the first degree of the argument, and it will be used as a derivative for the construction of the spline. The derivative at the right edge of the interval is calculated similarly. As a result of a simple calculation, we obtain the following system of equations for the spline coefficients:
f (i) = y i
f '(i) = (y i + 1 - y i-1 ) / 2
where i = {0, 1}
whose solution can be obtained in the following form:
y ' 0 = (y 1 - y -1 ) / 2
y ' 1 = (y 2 - y 0 ) / 2
a 0 = y 0
a 1 = y ' 1
a 3 = y ' 0 + y' 1 + 2 (y 0 - y 1 )
a 2 = y 1 - a 3 - a 1 - a 0
We can note that in the case of sequential (streaming) processing of input data, the derivative of the left boundary of the interval is the same as the derivative of the right boundary of the previous interval, which means that at each interval we need to calculate only one derivative - in the right boundary of the interval, and take the value in the left from memorized in the previous calculation in the right. Using this saving on matches, we obtain the following number of operations for calculating the polynomial coefficients:
This particular case of the Hermite spline is called Catmulla-Roma.
If we estimate the derivatives in the edges of the central interval not by three, but by five points, then their calculation will become a bit more complicated, but the other formulas will not change:
y ' 0 = ((y 1 - y -1 ) - (y 2 - y -2 ) / 8) 2/3
y ' 1 = ((y 2 - y 0 ) - (y 3 - y -1 ) / 8) 2/3
a 0 = y 0
a 1 = y ' 1
a 3 = y ' 0 + y' 1 + 2 (y 0 - y 1 )
a 2 = y 1 - a 3 - a 1 - a 0
Taking into account the previous remark about the need to calculate only one derivative at each interval, the number of operations:
By the number of operations, the spline of Catmulla-Roma significantly outperforms Farrow, Hermite 6 points slightly losing (by one addition). It is also worth noting that both splines have first order smoothness (continuous zero and first derivatives) due to their construction, unlike Farrow, which does not ensure the continuity of the derivative. But there remains the question of the accuracy of the approximation of functions by these splines. The sine was chosen as a test function, it was approximated by the above spline variants for a different grid step, which is normalized to the relative value of the number of points per period. The magnitude of the maximum deviation from the approximated function was measured. The results are presented on the graph. When choosing the logarithmic scale, the dependences on both axes are almost linear, converging to one point as the sampling frequency decreases and approaches the Nyquist frequency (2 points per approximated test sinus period).
The graph shows that the spline of Catmulla-Roma is inferior to Farrow in accuracy, and the latter is inferior to Hermite by 6 points, with almost the same number of operations required for the calculation. It was not possible to win a clear victory, but in principle, the result is not bad.
We can develop the general idea of the Hermite spline on derivatives of any order. For example, we construct a spline according to the following conditions: its values at the edges of the interval coincide with the values of the samples at these points, and the second derivatives at these points coincide with their estimates produced by parabolas through triples of points - that is, everything is like Catmull-Roma, only instead of the first, we will take the second derivatives - parabola, as a second-order polynomial, which will allow us to calculate them completely. And here we are surprised to find that the spline constructed in this way fully coincides with our esteemed Farrow, who is actually Lagrange polynomial. Those who wish can check this, I will even write a formula for evaluating the second derivative of a parabola on 3 points (although this is a trivial exercise):
y '' i = y i + 1 + y i-1 - 2 y i
Good match? Random because there are few points and small polynomial order? Numerically checking this fact on a different number of points and the order of polynomials, I was convinced that the spell works:
Suppose we have an interpolation polynomial constructed by some values on an even number of points of a uniform grid.
x -k , x - (k-1) , ... x 0 , x 1 , ... x k-1 , x k , x k + 1
Then all the even derivatives of this polynomial at the point x 0 coincide with the derivatives of the same orders of the interpolation polynomial constructed from the points
x -k , x - (k-1) , ... x 0 , x 1 , ... x k-1 , x k
and at the point x 1 - a polynomial constructed by the points
x - (k-1) , ... x 0 , x 1 , ... x k-1 , x k , x k + 1
Evidence. Consider a function on a uniform grid of an odd number of points — let's make an argument offset, giving the central node argument to 0
x -k , x - (k-1) , ... 0, x 1 , ... x k-1 , x k
There is its only interpolation polynomial P 0 (x). Add to this set of points one more x k + 1 (for the proof it does not matter, we add it to the right, left, or even to the middle and not to the nodes of the original grid). Interpolation polynomial on an extended set of points, provided a uniform grid can be represented in the form of Newton
P 1 (x) = P 0 (x) + A (xx -k ) (xx - (k-1) ) .... (x-0) .... (xx k-1 ) (xx k )
where A is a constant. Indeed, this additive additive requires zero values at the points of the original polynomial, and the constant A is determined from the equality condition P 1 (x k + 1 ) = y k + 1 - the required value of the polynomial at the added point. Due to the uniformity of the grid -x -i = xi , it means
P 1 (x) = P 0 (x) + A (x 2 -x k 2 ) (x 2 -x k-1 2 ) .... (x 2 -x 1 2 ) (x-0)
Obviously, the right side is the sum of P 0 (x) and a polynomial with nonzero coefficients only for odd powers of the argument. Hence, the coefficients for even degrees of the polynomials P 1 (x) and P 0 (x) coincide, and all even derivatives of any polynomial at zero are determined only by its corresponding coefficient with even degree of argument. Thus, all the even derivatives of P 1 (x) and P 0 (x) coincide.
The above reasoning is valid for any location of the added point x k + 1 , in particular, when adding it to the right or left of the original grid. Thus, the derivatives of even orders at the origin of the original polynomial and the polynomial with respect to the extended set of points coincide. PM
Since the interpolation polynomial of minimal order is unique, we can, taking into account the spell we have just proved, build it through conditions on the values of even derivatives at the edges of the central interval — for example, a 5th order polynomial through 6 points according to 6 conditions: values of counts and second and fourth derivatives calculated for every five points. In this case, use the same optimization described above - calculation of derivatives only in the right boundary of each interval. In the original formulation of passing a polynomial through all points, these invariants were not obvious. Apply this to our original Lagrange by 4 points and optimize the number of operations. I’ll give the Matlab code of my version right away, and the participant's nickname _Anatoliy from the electronic forum, which independently and independently defeated Farrow before my experiments:
As you can see, both codes contain fewer operations for calculating the polynomial coefficients, so in terms of sports, I believe that we both won. Although there are opinions that this is saving on matches, that nowadays, whenspace ships plow everywhere there are mathematical co-processors with integrated operations with floating-point numbers, the speed of operations is comparable to the speed of moving the contents of registers and especially reading / writing memory to save / recovery of calculated values - such optimizations do nothing, I still have a certain feeling of satisfaction with the results, despite their weak practical significance.
- Where do you see giants? - Sancho Panza asked.
- Yes, there they are, with huge hands, - answered his lord. - Some of them have arm lengths of almost two miles.
“Have mercy, señor,” Sancho said, “that what is seen there is not at all giants, but windmills; what you take for their hands are the wings: they whirl from the wind and set the millstones in motion.
“You can see an inexperienced adventure seeker,” said Don Quixote, “these are giants.” And if you are afraid, then drive off to the side and pray, and in the meantime I will join with them in a fierce and unequal battle ...
Once a friend of mine asked me to simulate a simple and low-cost variant of resampling (changing the sampling frequency) of a digital signal. There were no special requirements for linear frequency response and phase response, and the desire to save time and memory was, therefore, some methods such as polyphase interpolating banks of FIR filters were immediately excluded from consideration, and the main attention was paid to local polynomial interpolation methods. Polynomials are quick and easy to read, storage of their coefficients requires little memory, and using a single set of coefficients, you can calculate any values within the original sampling interval.
Meet the enemy - the implementation of Farrow
It is known that through any N points (from here and further for definiteness we will consider a uniform grid - with the same step, although for many described techniques and methods this is not a prerequisite), so, through any N points with different x-coordinates, you can conduct a single polynomial of N-1 order uniquely determined by a set of N coefficients (read that LaTex on Habré is not directly, but through pictures that are heavy and short-lived, so I will only manage the upper-lower indices in html-tags):
')
f (t) = a n t n + a n-1 t n-1 + ..... + a 1 t + a 0
For convenience of calculations, we normalize our uniform argument grid to integer values, using the example of four points, choose a convenient interval of the argument (Farrow chose from [-2, 1] with central interval from [-1, 0]) and write the system of linear equations from the condition of passing the polynomial through our 4 points:
f (i) = y i
where i = {-2, -1, 0, 1}
You can decide as you like, even though Gauss can convert the matrix, but Farrow has optimized the number of operations to calculate the coefficients of this polynomial:
a 3 = (y 1 - y -2 ) / 6 + (y -1 - y 0 ) / 2
a 1 = (y 1 - y -1 ) / 2 - a 3
a 2 = y 1 - y 0 - a 1 - a 3
a 0 = y 0
Total operations:
- 1 multiplication by a constant (1/6)
- 2 divisions by 2 (considered fast and in integer arithmetic and in floating point format)
- 8 addition / subtraction
This polynomial is used to calculate values that fall in the central interval of the original four points: [-1, 0] (with the corresponding rationing of the argument); for other intervals, the polynomials are calculated in the same way from the four points surrounding them. Very popular interpolator in narrow circles, deserved and revered. And for some reason I immediately wanted to come up with something that could surpass this implementation in a certain sense, taking into account the number of operations for the calculation.
First challenger - Catmulla-Roma spline
You can see that both the sought polynom itself and all its derivatives are linear combinations of the unknown coefficients of the polynomial. Therefore, we can set the conditions for the equality of derivatives of any order at the nodes of the grid to the given values, and obtain a linear system of equations for the polynomial coefficients. Hermite's local cubic spline is based on four conditions: the values of the polynomial and its first derivative at the edges of the interval. And we can get the approximate values of these derivatives at the right points by differentiating the interpolating Lagrange polynomial passing through them. For example, we normalize the arguments of four points to the interval [-1, 2], take three points with arguments {-1, 0, 1} and draw a parabola through them, and calculate its coefficients. And since we only need the value of the first derivative of this parabola at the center point 0, it is enough to calculate its coefficient for the first degree of the argument, and it will be used as a derivative for the construction of the spline. The derivative at the right edge of the interval is calculated similarly. As a result of a simple calculation, we obtain the following system of equations for the spline coefficients:
f (i) = y i
f '(i) = (y i + 1 - y i-1 ) / 2
where i = {0, 1}
whose solution can be obtained in the following form:
y ' 0 = (y 1 - y -1 ) / 2
y ' 1 = (y 2 - y 0 ) / 2
a 0 = y 0
a 1 = y ' 1
a 3 = y ' 0 + y' 1 + 2 (y 0 - y 1 )
a 2 = y 1 - a 3 - a 1 - a 0
We can note that in the case of sequential (streaming) processing of input data, the derivative of the left boundary of the interval is the same as the derivative of the right boundary of the previous interval, which means that at each interval we need to calculate only one derivative - in the right boundary of the interval, and take the value in the left from memorized in the previous calculation in the right. Using this saving on matches, we obtain the following number of operations for calculating the polynomial coefficients:
- 2 multiplication / division by 2
- 7 addition / subtraction
This particular case of the Hermite spline is called Catmulla-Roma.
The second challenger - Hermite spline at 6 points
If we estimate the derivatives in the edges of the central interval not by three, but by five points, then their calculation will become a bit more complicated, but the other formulas will not change:
y ' 0 = ((y 1 - y -1 ) - (y 2 - y -2 ) / 8) 2/3
y ' 1 = ((y 2 - y 0 ) - (y 3 - y -1 ) / 8) 2/3
a 0 = y 0
a 1 = y ' 1
a 3 = y ' 0 + y' 1 + 2 (y 0 - y 1 )
a 2 = y 1 - a 3 - a 1 - a 0
Taking into account the previous remark about the need to calculate only one derivative at each interval, the number of operations:
- 1 multiplication by a constant (2/3)
- 2 multiplication / division by grade 2
- 9 addition / subtraction
Comparison of options
By the number of operations, the spline of Catmulla-Roma significantly outperforms Farrow, Hermite 6 points slightly losing (by one addition). It is also worth noting that both splines have first order smoothness (continuous zero and first derivatives) due to their construction, unlike Farrow, which does not ensure the continuity of the derivative. But there remains the question of the accuracy of the approximation of functions by these splines. The sine was chosen as a test function, it was approximated by the above spline variants for a different grid step, which is normalized to the relative value of the number of points per period. The magnitude of the maximum deviation from the approximated function was measured. The results are presented on the graph. When choosing the logarithmic scale, the dependences on both axes are almost linear, converging to one point as the sampling frequency decreases and approaches the Nyquist frequency (2 points per approximated test sinus period).
The graph shows that the spline of Catmulla-Roma is inferior to Farrow in accuracy, and the latter is inferior to Hermite by 6 points, with almost the same number of operations required for the calculation. It was not possible to win a clear victory, but in principle, the result is not bad.
Spell
We can develop the general idea of the Hermite spline on derivatives of any order. For example, we construct a spline according to the following conditions: its values at the edges of the interval coincide with the values of the samples at these points, and the second derivatives at these points coincide with their estimates produced by parabolas through triples of points - that is, everything is like Catmull-Roma, only instead of the first, we will take the second derivatives - parabola, as a second-order polynomial, which will allow us to calculate them completely. And here we are surprised to find that the spline constructed in this way fully coincides with our esteemed Farrow, who is actually Lagrange polynomial. Those who wish can check this, I will even write a formula for evaluating the second derivative of a parabola on 3 points (although this is a trivial exercise):
y '' i = y i + 1 + y i-1 - 2 y i
Good match? Random because there are few points and small polynomial order? Numerically checking this fact on a different number of points and the order of polynomials, I was convinced that the spell works:
Suppose we have an interpolation polynomial constructed by some values on an even number of points of a uniform grid.
x -k , x - (k-1) , ... x 0 , x 1 , ... x k-1 , x k , x k + 1
Then all the even derivatives of this polynomial at the point x 0 coincide with the derivatives of the same orders of the interpolation polynomial constructed from the points
x -k , x - (k-1) , ... x 0 , x 1 , ... x k-1 , x k
and at the point x 1 - a polynomial constructed by the points
x - (k-1) , ... x 0 , x 1 , ... x k-1 , x k , x k + 1
Evidence. Consider a function on a uniform grid of an odd number of points — let's make an argument offset, giving the central node argument to 0
x -k , x - (k-1) , ... 0, x 1 , ... x k-1 , x k
There is its only interpolation polynomial P 0 (x). Add to this set of points one more x k + 1 (for the proof it does not matter, we add it to the right, left, or even to the middle and not to the nodes of the original grid). Interpolation polynomial on an extended set of points, provided a uniform grid can be represented in the form of Newton
P 1 (x) = P 0 (x) + A (xx -k ) (xx - (k-1) ) .... (x-0) .... (xx k-1 ) (xx k )
where A is a constant. Indeed, this additive additive requires zero values at the points of the original polynomial, and the constant A is determined from the equality condition P 1 (x k + 1 ) = y k + 1 - the required value of the polynomial at the added point. Due to the uniformity of the grid -x -i = xi , it means
P 1 (x) = P 0 (x) + A (x 2 -x k 2 ) (x 2 -x k-1 2 ) .... (x 2 -x 1 2 ) (x-0)
Obviously, the right side is the sum of P 0 (x) and a polynomial with nonzero coefficients only for odd powers of the argument. Hence, the coefficients for even degrees of the polynomials P 1 (x) and P 0 (x) coincide, and all even derivatives of any polynomial at zero are determined only by its corresponding coefficient with even degree of argument. Thus, all the even derivatives of P 1 (x) and P 0 (x) coincide.
The above reasoning is valid for any location of the added point x k + 1 , in particular, when adding it to the right or left of the original grid. Thus, the derivatives of even orders at the origin of the original polynomial and the polynomial with respect to the extended set of points coincide. PM
Last round - Lagrange after Farrow with Lagrange via derivatives
Since the interpolation polynomial of minimal order is unique, we can, taking into account the spell we have just proved, build it through conditions on the values of even derivatives at the edges of the central interval — for example, a 5th order polynomial through 6 points according to 6 conditions: values of counts and second and fourth derivatives calculated for every five points. In this case, use the same optimization described above - calculation of derivatives only in the right boundary of each interval. In the original formulation of passing a polynomial through all points, these invariants were not obvious. Apply this to our original Lagrange by 4 points and optimize the number of operations. I’ll give the Matlab code of my version right away, and the participant's nickname _Anatoliy from the electronic forum, which independently and independently defeated Farrow before my experiments:
clf reset N = 10; x = 1:N; y = rand(1, N); h = 0.01; g = y(3) - y(2); d = g - y(2) + y(1); plot(x, y, 'or', 'LineWidth', 2); hold on; grid on; axis on; title([' ', ... ' 3 , 2 .']); plot(x(1), y(1), 'b-', x(1), y(1), 'g:'); legend('', ' ', ' Anatoliy'); for k = 2:(N-2) % : 1 , 1 , 5 c = d; e = g; g = y(k+2) - y(k+1); d = g - e; b2 = c/2; b3 = (d - c)/6; b1 = e - b2 - b3; t = 0:h:1; f = b3.*t.^3 + b2.*t.^2 + b1.*t + y(k); plot(t+x(k), f, 'b-') % _Anatoliy: 1 , 2 , 6 p = y(k) - y(k+1); q = (y(k+2) - y(k))/2; a3 = ( (y(k+2) - y(k-1))/3 + p )/2; a2 = p + q; a1 = q - a3; t = -1:h:0; f = a3.*t.^3 + a2.*t.^2 + a1.*t + y(k+1); plot(t+1+x(k), f, 'g:') end As you can see, both codes contain fewer operations for calculating the polynomial coefficients, so in terms of sports, I believe that we both won. Although there are opinions that this is saving on matches, that nowadays, when
Source: https://habr.com/ru/post/282441/
All Articles