An extended Prony’s interpolation scheme on an equispaced grid

Abstract An interpolation scheme on an equispaced grid based on the concept of the minimal order of the linear recurrent sequence is proposed in this paper. This interpolation scheme is exact when the number of nodes corresponds to the order of the linear recurrent function. It is shown that it is still possible to construct a nearest mimicking algebraic interpolant if the order of the linear recurrent function does not exist. The proposed interpolation technique can be considered as the extension of the Prony method and can be useful for describing noisy and defected signals.


Introduction
It is well known that the n-point polynomial interpolants in equally spaced points to a function f on [-1, 1] do not necessarily converge as n ! 1, even if f is analytic. Instead one may see wild oscillations near the endpoints, an effect known as the Runge phenomenon [9]. Moreover, the interpolation process becomes exponentially illconditioned, as shown first by Turetskii [13] and later independently by Schonhage [11]. This ill-conditioning means that even if the interpolants converge in theory, they will diverge in floating point arithmetic, at least for values of x near the endpoints, because of exponential amplification of rounding errors.
On the other hand, the polynomial interpolation in Chebyshev points is numerically stable since the associated Lebesgue constants are of size O .log .n// [3]. Is is shown in [4] that the Chebyshev interpolant can be evaluated in floating point arithmetic by Salzer's barycentric formula [10]. Moreover, Chebyshev interpolants are used in the Chebfun software where polynomials in degrees of tens of thousands are routinely used for practical computations [7,12].
The main objective of this paper is to propose an extended Prony-type algebraic interpolation scheme on an equispaced grid. Our goal is to develop a strategy for finding a nearest algebraic interpolant to an analytic function. It will be demonstrated that an approach based on the nearest algebraic interpolant does not suppress the Runge phenomenon, but interpolation errors produced by this method are much lower if compared to the classic schemes on equispaced grids. Moreover, the proposed algebraic interpolation method can be effectively exploited for analytic interpolation of noisy and/or defected signals. This paper is organized as follows. Linear recurring functions and linear recurring sequences are discussed in Section 2; computational examples of linear recurring sequences versus Prony decomposition are shown in Section 3; Extended Prony's functions and their properties are given in Section 4; algorithm of the extended Pronytype interpolation is presented in Section 5; the interpolation of the real time series is investigated in Section 6; concluding remarks are given in Section 7.

Preliminaries
The definition of the linear recurring function, the linear recurring sequence and their properties will be presented in this section.

A short overview of the Prony method
Fourier expansions are successfully exploited for the approximation of a function f .x/ in a variety of theoretical and practical applications. And even though higher terms of the Fourier series are usually infinitesimal, that may cause substantial complications. In contrast, Prony's method computes an approximation to function f .x/ by using only a finite number of damped complex exponentials [14]: where m 2 N and c n ; n 2 C. Prony-type methods are successfully exploited for effective and accurate approximation of different functions and signals [15][16][17][18] (though the determination of m remains problematic in general). As mentioned previously, the main objective of this paper is to use the extended Prony-type scheme to interpolate a function on an equispaced grid. The extended Prony scheme comprises complex exponentials and polynomials [1,2,6]: where m 2 N; n 2 C and Q n .x/ are polynomials in x with complex coefficients and non-negative integer powers of x.

The minimal order of linear recurring sequence
Let us consider a sequence: p 0 ; p 1 ; p 2 ; ::: WD .p j I j 2 Z 0 / where elements p j can be real or complex numbers. Then, a sequence of Hankel matrices reads: H n WD .p i Cj 2 / 1Äi;j Än D 2 6 6 6 4 p 0 p 1 ::: p n 1 p 1 p 2 ::: p n ::: p n 1 p n ::: p 2n 2 3 7 7 7 5 I n D 1; 2; ::: : The Hankel transform (the sequence of determinants of Hankel matrices) .d n I n 2 N / reads: Definition 2.1. The minimal order of the recurring sequence .
Note that rk j k The opposite statement holds also. If (4) holds true, then rank.p j I j 2 Z 0 / D m 1 C m 2 C ::: C m l : Rigorous proof of this theorem is given in [5].
Definition 2.5. A sequence .p j I j 2 Z 0 / is a linear recurring sequence (LRS) if elements of that sequence can be expressed in the form of (4).
For example, Thus, the sequence q j D P l rD1 P m r 1 kD0 O rk j k Á j r ; j D 0; 1; 2; : : : is an LRS and its minimal order of LRS is equal to rank.q j I j 2 Z 0 / D m 1 C m 2 C ::: C m l because O rm r 1 D rm r 1 ¤ 0.
This system of linear equations has a unique solution [5].
Remark 2.8. Let rank.p j I j 2 Z 0 / D m and the first 2m elements of that series are known. Then it is possible to use (3), (6) and (4) to calculate all elements of that sequence.
Corollary 2.9. Let the sequence .p j I j 2 Z 0 / be an LRS and its minimal order be equal to m. Then the following linear recurrence relation holds true [5]: B 0 p j C B 1 p j C1 C ::: C B m 1 p j Cm 1 D p j Cm I j D 0; 1; ::: where constants B 0 ; B 1 ; :::; B m 1 2 C exist and do not depend on j .
On the other hand (according to (2)), the characteristic polynomial yields one root 1 D 1; its multiplicity index is m 1 D 3. Therefore, according to (4), Coefficients 10 ; 11 ; 12 are determined in order to fit the recurrence: 10 D 0I 11 D 1I 12 D 2. Thus, finally, Remark 2.11. The determination of the minimal order of sequence fp j g by using Definition 2.1 has a cost of

LRS versus Prony decomposition -two computational examples
Let us consider a sequence p j ; j D 0; :::; N where Q s; N 2 N, 3 Ä Q s Ä N , Q s is an upper bound of the number of exponentials and the bounds 0;1 are positive. An LRS of a given sequence could be found using the Prony interpolation method [21]. The main steps of this algorithm read: 1. Determine the smallest singular value of the rectangular Hankel matrix (H WD .p oCg / N Q s; Q s oD1;gD0 and use singular value decomposition to find the related right singular vector u D .u l / Q s lD0 . 2. Compute all zeros of polynomial P Q s lD0 u l z l and determine all zeros Q z t , t D 1; :: :::; Q M / as least squares solution of the overdetermined linear Vandermonde-type system: Example 3.1. Let us consider a sequence p j WD j; j D 0; :::; N 1, N D 21. Let us find a LRS of a given sequence using two alternative methods: a) the concept of the minimal order of LRS; b) the Prony interpolation method [21].
a) The minimal order of LRS of a given sequence is 2 because the sequence of determinants of Hankel matrices reads .0; 1; 0; 0; :::/. Then the characteristic polynomial reads: The roots are: 1 ; 2 D 1. Now, the linear algebraic system of equations (6) yields: : Then the LRS of the given sequence reads: where Q r .x/ D P m r 1 k r D0 a rk r x k r ; m r 1; a rk r 2 C; a r;.m r 1/ ¤ 0; r D 0; 1; : : : ; n and x; f .x/ 2 R: It is important to note that n is finite in (8).   is a LRS (x 0 ; h 2 R are fixed parameters).

Extended Prony's functions and their properties
Proof. Let the function y D f .x/ be an EPF. Then, the following equalities hold for all x 0 ; h 2 R: where coefficients b rk can be expressed in terms of coefficients a r0 ; a r1 ; :::; a rm r 1 . It can be noted that the index k, parameters x 0 and h do not depend on j. The introduction of the symbol reduces (10) into the LRS: r I j D 0; 1; 2; : : : : Theorem 4.2. Let .y j I j 2 Z 0 / be an LRS. Then the following inequalities hold true: 0 Ä rank.y j I j 2 Z 0 / Ä m 1 C m 2 C ::: C m r : Proof. (5) can be used to express terms in (12) (coefficients b rk ¤ 0). Combining like terms (some roots r may coincide) yields the estimate in (13).
Proof. (11) and the definition of the complex logarithm yield: r h WD Ln r D ln j r j C i.arg r C 2 k r /I k r D 0;˙1;˙2; :::I r D 1; l: It is important to select such coefficients k r that different r would correspond to different r for all r. That ensures the equality of the EPF.
It can be noted that the selection of the parameter h may be a non-trivial task in the general case. Therefore, it is important to define the concept of the sufficiently small step h. .ln j r j C i arg r / I where 0 Ä jarg r j < and r D 1; l. Proof. (14) yields: a r h D jln r j; b r h D arg r C 2 k r .h/ where r D a r C i b r and 0 < jb r j < C1; a 2 R. Thus, b r h D arg r while 0 < jb r j h < . Therefore, k r .h/ D 0 and 0 < h < jbr j . The selection of h 0 WD min r jbr j finalizes the proof. It can be noted that h 0 would be unbounded if max Proof. Let i.e. r WD r .h/ is a function of h; 0 < h < C1. Then, b r h D arg r .h/ C 2 k r .h/: Example 4.13. It can be noted that the step h 0 D is a sufficiently small step for f .x/ D cos x. Then parameters k and l must be set to 0.

The extended Prony-type interpolation method
Let us assume that a function f .x/ is not necessarily EPF. Then the reconstruction of the closest EPF to the function f .x/ in the interval a Ä x Ä b would be an important practical problem which is discussed in details in this section. First of all the step h of the regular grid in the interval OEaI b should be selected. The selection of the step h is directly related to the order of the EPF F .x/ which will be used to mimic the original function f .x/. Let the order of F .x/ is m (we wish to mimic f .x/ using an EPF with the order equal to m). Then, the step h reads:  . 1 C 2 j C 3 j 2 /e j 1 C . 4 C 5 j C 6 j 2 /e j 2 C 7 e j 3 C 8 e j 4 D f .j h/; j D 0; 7: Solutions of the linear algebraic system of equation now reads: 1;2;4;5 D 0; 3;6 D 0:0667i ; 7;8 D 1 2 . Finally, the expression of F 8 .x/ reads: The algebraic interpolant is illustrated in Figure 4.

Algebraic interpolation of the real time series
The ability of algebraic interpolation on a regular grid (for the nearest algebraic interpolant) suggests interesting possibilities for application of this approximation scheme for the real time series. We use a time series of monthly PMI Composite Index. A PMI reading above 50 percent indicates that the manufacturing economy is generally expanding and below 50 percent that it is generally declining. The unfiltered PMI data S k in interval k D 1; 2; : : : ; 788 is shown in Figure 5. Circles denote nodes of the equispaced grid; the zoomed part of the algebraic interpolant is illustrated in part (b).
We will use the algebraic interpolation scheme in equally spaced grid on interval [1; 788]. First, we preselect the order of the EPF which will be used to mimic the experimental data. For the order of the EPF equal to m, the step h D 787 2m 1 and p n D S j ; j D 1 C nh; n D 0; 1; : : : ; 2m 1. We again perform a number of computational experiments for different values of m; m D 1; 2; : : : ; 90. The algorithm of algebraic interpolation produces 90 different EPF F m .x/. RMSE of algebraic interpolation is now computed as the square root of the sum of squared differences between values of the given data and the values of the EPF F m .x/ at all sampling points in the interval [1; 788]. Computational experiments show that the best result (the minimal RMSE = 3.17) is achieved at m D 56 (Figure 6). The graph of F 56 .x/ is shown in Figure 5. As noted previously, the extended Prony-type interpolation scheme outperforms the Lagrange polynomial on equispaced grinds. The Lagrangian interpolant is illustrated in Figure 7 (the nodes are the same as in Figure 5). Runge's effect prevents Lagrange interpolation to be a reasonable approximation (the maximum absolute value of the Lagrange interpolant in Figure 7 is equal to 3:1634 10 32 ), yet the proposed extended Prony-type interpolation does not have that problem. As mentioned in the Introduction, the polynomial interpolation in Chebyshev points is numerically stable even for high polynomial degrees. However, real-world time series are usually recorded using a constant sampling rate. Thus, it would be complicated to find values of this time series at Chebyshev points without transforming the scale (Chebyshev interpolation is straightforward for continuous functions of course). F 56 .x/ in Figure 5 is a good example of such alternative interpolation.

Concluding remarks
An algebraic interpolation scheme on an equispaced grid is presented in this paper. It is demonstrated that the proposed scheme can be used for the identification of the nearest algebraic interpolant for a given function. This computational effect can be explained by the fact that though all nodes are located inside the bounded interval, the interpolant is reconstructed in the global domain. Such interpolation scheme can be extended to the extrapolation scheme what can be successfully exploited in time series prediction applications [8]. On the other hand, this advantageous feature can be successfully exploited for the analytic approximation of noisy and/or defected realworld signals.
Numerical experiments have shown that optimal interpolants produced by the proposed algebraic technique based on the order of LRS outperform classical Lagrange polynomial interpolants on equispaced grids. This effect can be explained by the fact that the functional base used to construct algebraic interpolants is wider compared to polynomial interpolants. Really, (8) would represent a polynomial interpolant if all indexes r I r D 0; 1; : : : ; n would be equal to zero.
The main drawback of the proposed interpolation scheme is that rather complex computations are required as the number of nodes becomes large. In the first place this is associated with the necessity to find all roots of the characteristic Hankel equation. Thus, the proposed scheme of interpolation looses its aptitude when the number of nodes becomes higher than one hundred. Nevertheless, the scheme preserves interesting potential of practical applications at lower number of nodes. The explicit error bound of the interpolation and the possibility of using adaptive grids remains a definite objective of future research.
Finally, it can be noted that the proposed scheme can be considered as an effective numerical tool for the identification of nearest algebraic "skeleton" functions and extends the applicability of classical interpolation schemes to real-world data contaminated with the inevitable noise.