Numerical Continued Fraction Interpolation

We show that highly accurate approximations can often be obtained from constructing Thiele interpolating continued fractions by a Greedy selection of the interpolation points together with an early termination condition. The obtained results are comparable with the outcome from state-of-the-art rational interpolation techniques based on the barycentric form.


1.
Introduction. An alternative title to this note could be The unreasonable effectiveness of Thiele continued fractions. Indeed, it is known [2,4] and further detailed in Section 3.5 that the construction of Thiele interpolating continued fractions can suffer from numerical instabilities. As argued in [4], a careful selection of the ordering of the interpolation points is therefore needed. Here we propose a Greedy strategy, one which takes the next point where the error is maximal. It is shown in Theorem 3.2 that, at least from an existence point of view, this strategy is well motivated. We cannot prove that it is the overall best strategy, but we rather give ample numerical evidence that the approach can lead to results that are competitive with other well-known approaches [12,9,3] for univariate rational interpolation and that the observed behavior is typical.
Once the existence problem is out of the way, we show how also best approximations can be obtained using the approach from [5] and how poles and zeros can be calculated directly from a generalised eigenvalue problem in Theorem (B.1). Combined, these cover most of the important tools required from a numerical rational interpolation method.
2. Thiele continued fraction rational interpolation. Consider a sequence of distinct complex points (x i ) i∈N together with function valuations f (x i ) = f i of the complex function f (x). The continued fraction generates rational interpolants when considering its successive convergents. The nth convergent has the property that C n (x i ) = f i for i = 0, . . . , n, provided that the inverse differences ϕ i [x 0 , . . . , x i ] = ∞ exist and none of the tails vanish at x i−1 for i > 0. When T i,n (x i−1 ) = 0, then x i−1 is unattainable, meaning that it is a common zero of the numerator and denominator of C n (x) and in addition 3. Existence, sampling and numerical stability.
110] that the numerator and denominator degrees of the convergents C n (x) are at most n/2 and n/2 respectively. If desired, more off-diagonal interpolants can be obtained by considering f (x) − p (x) instead of f (x), with p (x) an interpolating polynomial of degree at most . The introduction of such a polynomial could be useful to dampen functions f (x) with a wide varying range. However, in what follows our interest is purely in Thiele continued fractions, hence we do not explore further potential benefits of non-diagonal interpolants.
3.2. Unattainability. While unattainability is undesirable, it is rarely encountered in practice except for pathological cases which are mostly related to symmetry of f (x). Therefore we only mention that for i = 1, . . . , n, it can be checked a posteriori if desired: either by checking for zeros of the tails T i,n (x) at x i−1 or from direct evaluation of the convergents C n (x) at the interpolation points. The next theorem says that a Thiele continued fraction is irreducible, up to possible unattainable points.
Theorem 3.1. Provided that the convergents C n (x) exist, its only common factors are of the form (x − x i ) with x i ∈ {x 0 , . . . , x n−1 } an interpolation point. In such a case the x i are unattainable and are characterized by vanishing of the tails T i+1,n (x i ) = 0.
Proof. Elaborate proofs can for instance be found in [18,4].

Existence.
Non-existence of inverse differences on the other hand is more inconvenient as it renders contributions from subsequent tails meaningless. The existence problem depends entirely on the ordering of the (x i ) i∈N . Specifically, it is required for two consecutive convergents C i−1 (x) and C i (x) to be different (i > 0) in order for the inverse differences ϕ i [x 0 , . . . , x i ] to exist.
Theorem 3.2. If the points (x i ) i∈N are ordered such that every two consecutive convergents of the continued fraction (2.1) are different, then ϕ i [x 0 , . . . , The proof is detailed in Appendix A, where it is shown that the inverse differences can essentially be interpreted as the ratio of two (linearized) residuals of successive convergents. This interpretation is used next to motivate a Greedy selection of the interpolation points.
3.4. Adaptive Greedy selection. Given a finite sequence of points (x i ) 0≤i≤n , the condition of Theorem 3.2 where successive convergents should be different gives us a heuristic way to choose the next point x i+1 in the construction of C i+1 (x) with 0 < i < n. Given C i (x) and (x 0 , . . . , x i ), reorder the remaining points (x i+1 , . . . , x n ) to determine C i+1 (x) such that |C i (x i+1 ) − f (x i+1 )| is maximal. In this way, C n (x) is ultimate constructed in an adaptive Greedy way by choosing in each step the point where for 0 < i < n the error between C i (x) and f (x) is maximal. For the first point x 0 one can take a point where |f (x 0 )| is minimum. As such, at least one zero of f (x) is accurately represented when present in the data.
This Greedy selection ensures the existence of the inverse differences (see Appendix A). A similar strategy is adopted by the AAA approach [9], although there the choice is motivated by numerical considerations rather than existence. We cannot prove that the Greedy strategy is best for numerical stability in the context of Thiele continued fractions (discussed next), but the examples in this note suggest that the observed stable behaviour is typical.
3.5. Numerical stability. It is known [2] that the computation of the inverse differences and successive convergents in a Thiele interpolating continued fraction can suffer from numerical instabilities and worst case exponential loss of precision may occur in the calculation of the inverse differences [4]. Nevertheless, the backward evaluation of the continued fraction often leads to near machine precision (roughly 2e-16) magnitudes for the residuals f (x i ) − C n (x i ) in practice [2].
We do not attempt an analysis like in [7] because any error analysis would necessarily depend on both the function f (x) and the chosen interpolation points x i . For instance, a very strong stability result is recently obtained in [1] when restricting to Markov functions. Instead, we provide compelling numerical evidence in Sections 4 and 5 that accurate approximations are often obtained, even on hard numerical problems, when the Greedy selection strategy of Section 3.4 is employed.
3.6. Early termination. From the discussion in the previous section, it is important to obtain the Thiele interpolating continued fractions in as few steps as possible. Besides the greedy selection, we add a stopping criterion when constructing C n (x). If the maximum absolute error in the remaining points is below a prescribed tolerance, say tol=5e-15, then we stop the recursion. The precise condition used is It essentially tells us that within numerical tolerance, the function to be approximated is rational. This can also be understood from the deeper connection to the shape of the blocks in the Walsh table [17], a topic we do not address further.

Example: interpolation of |x| and √
x. To illustrate that the Greedy selection strategy works well in practice, we start with the example of interpolation of f (x) = |x| in Newman [11] points for x ∈ [−1, 1]. This problem demonstrates the approximation power of rational functions as compared to polynomials which can only achieve O(n −1 ) accuracy at best. However, it has proven quite challenging numerically and has been used to assess the performance of several rational approximations schemes based on for instance the barycentric form [12,9,6].
Newman approximations are rational interpolants in the 2n + 1 points These points cluster around the origin, approaching it at an exponential rate as n increases. It is shown in [19] that the asymptotic rate of convergence of rational interpolants to |x| is root-exponential O(n −1/2 e − √ n ). Figure 1 shows the results of Thiele interpolation in Newman points and f (x) = |x| with n up to 50. All obtained interpolants use all interpolation points in their construction, meaning that for n = 50 we have 2n + 1 = 101 interpolation points and we construct C 2n (x) = C 100 (x). Mind that this is a challenging example because C 100 (x) is a rational function of degree n = 50 and n = 50 in numerator and denominator respectively. For instance the robust approach [12, see Fig   An equivalent, but computationally easier, problem is to approximate f (x) = √ x on x ∈ [0, 1]. For this interpolation problem we take the square of the Newmanpoints (4.1) leading to n + 1 unique points in [0, 1]: Note that these points cluster even closer near the origin than the original Newman points (4.1). Figure 2 shows the results of Thiele interpolation in squared Newman points and f (x) = √ x with n up to 400. The error in the interpolation points seems to decay logarithmically with n. However, this is merely a consequence of the fact that for n > 60 the interpolants are no longer of full degree. Beyond that point, not all given interpolation points are used in the continued fraction construction and their degree only rises slowly. In fact, for n = 400 only 116 points are used after which the early termination kicks in. Nevertheless, the overall approximation quality improves in a root-exponential manner. We note that our implementation of the AAA method [9] fails for this example when n > 50 as it introduces real poles inside the interpolation interval.  |f (x) − C n (x)|.
Recent advances in best rational approximation have been made in [3,10,5]. We refer the interested reader to those references for further details on best approximations. Provided that C * n (x) is of full degree, it is known that there exists a so-called alternant set consisting of n + 2 ordered nodes where |f (x) − C * n (x)| attains its global extremum over all x ∈ [a, b] with alternating signs: f (x j ) − C * n (x j ) = (−1) δ+j λ, j = 0, . . . , n + 1.
Here δ ∈ {0, 1} and λ = max x∈[a,b] |f (x) − C * n (x)|. The alternating set is usually the starting point for the iterative Remez algorithm. From an interpolation point of view, the alternant set is not immediately useful because the construction of C n (x) requires only n + 1 points and the value λ is unknown. However, we can use an alternative approach as discussed next.
Due to continuity, the error f (x) − C * n (x) must attain zero between each pair (x j ,x j+1 ) of neighboring points of the alternant set. This means that there must exist at least n + 1 points x 0 , . . . , x n (5.2) a ≤x 0 < x 0 <x 1 < x 1 < · · · < x n <x n+1 ≤ b, such that f (x i ) = C * n (x i ) for i = 0, . . . , n. This idea is exploited in the BRASIL [5] algorithm which, given an initial guess of n+1 interpolation points, iteratively rescales the interval widths between successive interpolation points with the goal of equilibrating the local errors. In each step, a new set of interpolation points is determined. Figure 3 illustrates the best rational aproximation for f (x) = sin(20x)/(1 + 25x 2 ) on x ∈ [−1, 2]. To obtain the initial interpolation points, a low accuracy continued fraction C 49 (x) is constructed from 100 Chebychev points of the first kind between −1 and 2. Then the BRASIL iteration is ran, which relocates the interpolation points in each step. The interpolation itself is done using the adaptive Thiele continued fractions approach rather than using the barycentric form as in [5]. Other parameters, such as step size are put equal to 0.01 and convergence acceleration was not implemented. The final error f (x) − C 49 (x) equioscillates n + 2 = 51 times on x ∈ [−1, 2] between the maximum error ≈ 1.76e-08.  For the last illustration, we turn to the example of [15]. In Figure 4, the results for f (x) = √ x on x ∈ [0, 1] are shown. We found that a convenient starting point was obtained from a low(er) accuracy approximation based on 1000 linearly spaced points between 0 and 1 raised to the power 6 where the endpoints are removed. The iteration is started with n + 1 = 81 points adaptively chosen from the aforementioned set leading to the full degree interpolant C * 80 (x) with maximum leveled error ≈ 4.39e-12. Compared to the interpolants of similar size (i.e. constructed with the same number of interpolation points) in Figure 2a, the best approximation is roughly twice as accurate.
It is worth noting that one easily obtains an approximation of the same quality for f (x) = |x| on x ∈ [−1, 1] by taking C * 80 (x 2 ). Such an approximation would be equivalent to the one obtained in [3, Fig. 7.1] of degree 80 in numerator and denominator. Mind that here the interpolation points range over more than 20 order of magnitudes and some of them are even below machine precision. A similar phenomenon is observed in [5, Fig. 3]. To be fair, the direct best approximation of f (x) = |x| on x ∈ [−1, 1] more often fails using the current approach, because the continued fraction representation struggles to maintain the exact symmetry. x Approximation error

Conclusions.
We have shown that highly accurate approximations can often be obtained with Thiele continued fractions when incorporating a Greedy selection of the interpolation points together with an early termination condition. The obtained results are comparable with state-of-the art rational interpolation techniques based on the barycentric form such as AAA [9] and can be used for rational minimax approximation [5]. Also the poles and zeros are relatively easy to obtain, making the approach attractive for practical purposes.
The main advantage of the Thiele continued fraction is the simplicity of its construction. It does not require specialized linear algebra implementations such as SVD, but rather relies on an elementary recursion (see Appendix C for an implementaion). Of course accumulation of rounding errors does occur, hence inevitably one can find examples where the approach will break down. The possibility of numerical breakdown is applicable for all rational interpolation approaches, including AAA [10, see §6 p.3172-3173].

Appendix A. Proof of Theorem 3.2.
For the proof of Theorem 3.2 we first relate the construction of the inverse difference to the (linearized) residuals of the convergents C n (x) of the continued fraction (2.1). To that end, let C n (x) = A n (x)/B n (x) where the nth numerator A n (x) and nth denominator B n (x) satisfy the recurrence relation [16] (A.1) Theorem A.1. If the points (x i ) i∈N are ordered such that every two consecutive convergents of the continued fraction (2.1) are different then for i = −1, 0, 1, . . .
where the linearized residuals R i (x) are defined as Proof. The proof is by induction.
And for i = 0 Assume that the hypothesis ϕ i [x 0 , . . . , holds. First note that from application of (A.1) we have Provided that R i−1 (x k ) = 0 and application of the induction hypothesis gives For the last step, the presumed ordering of the points is important. If two consecutive convergents C i−1 (x) and C i (x) are different, then the polynomial Due to the condition (A.2) it cannot be that both R i−1 (x) and R i (x) simultaneously vanish at other points.
Proof of Theorem 3.2. From Theorem A.1, we know that the inverse differences can be expressed as and the presumed ordering of the (x i ) i∈N ensures that the ratio is well-defined because its numerator and denominator cannot vanish simultaneously. In addition, ϕ i+1 [x 0 , . . . , x i , x i+1 ] = ∞ would imply R i (x i+1 ) = 0 which contradicts that consecutive convergents C i+1 (x) and C i (x) are different.
Appendix B. Poles, zeros and residues. The poles and zeros of a Thiele continued fraction can be extracted directly from application of a (generalised) eigenvalue problem. The key ingredient to that end is the continuant [14] representation of the partial numerator and denominator. Even though continuants are well-known in the continued fraction literature, to the best of our knowledge their use to extract poles and zeros has not been exploited before.

B.1. Poles.
Theorem B.1. Given distinct points x i , finite ϕ i [x 0 , . . . , x i ] = ∞ (i = 0, . . . , n) and its associated Thiele continued fraction (B.1) In case that n is odd, assume: Let the n × n matrices D and C respectively be defined as then the matrix pencil D − λC is regular and its finite eigenvalues coincide with the n/2 poles of C n (x) counting multiplicities.
Proof. Using the three term recurrence (A.1), it is well-known [14] that one can write the denominator B n (x) of C n (x) as the determinant of an n × n tridiagonal matrix, which is also called a continuant .
Clearly (B.2) can be written in the form det(D − xC) where the n × n matrices D and C are as defined above. Under the given assumptions, B n (x) is of exact degree n/2 . Hence, B n (x) = det(D − xC) ≡ 0 so that the matrix pencil D − λC is regular. After all, in case that n is even, then B n (x) is monic of exact degree n/2; regardless of the values ϕ i [x 0 , . . . , x i ] = ∞ (i = 0, . . . , n). In case that n is odd, then B n (x) is of exact degree (n − 1)/2 = n/2 with highest degree coefficient equal to the sum in (B.1).
Since the characteristic polynomial det(D − λC) = B n (λ), the poles of C n (x) coincide with the non-zero eigenvalues of D − λC. There are exactly n/2 such eigenvalues, because B n (λ) is of exact degree n/2 .

Remark.1
The condition (B.1) is seldom of practical importance. It is more of theoretical interest to prevent the trivial case B n (x) ≡ 0 and the occurrence of additional eigenvalues at infinity when B n (x) is not of exact degree n/2 . Remark 2. In practice, the poles of C n (x) can thus be found by solving the generalised eigenvalue problem Dv = λCv.
The eigenvalues at infinity can be discarded. There are at least n/2 infinite eigenvalues (also see Remark 1).

B.2. Zeros.
In an entirely similar fashion, also the roots of C n (x) can be determined from a generalised eigenvalue problem. Based on the continuant representation of the numerator of C n (x), one simply replaces the n × n matrices D and C by the (n + 1) × (n + 1) matrices, where D is now and swaps every occurrence of odd by even in condition (B.1) of Theorem (B.1).
Recall that C n (x) has at most n/2 zeros counting multiplicities.
B.3. Residues. Once the poles are determined, the residues can be found in a similar fashion as in [9]. Essentially, one applies a trapezoidal rule approximation to the contour integral of C n (x) over a small circle around each pole in the complex plane.