The stability of extended Floater–Hormann interpolants

Abstract

We present a new analysis of the stability of extended Floater–Hormann interpolants, in which both noisy data and rounding errors are considered. Contrary to what is claimed in the current literature, we show that the Lebesgue constant of these interpolants can grow exponentially with the parameters that define them, and we emphasize the importance of using the proper interpretation of the Lebesgue constant in order to estimate correctly the effects of noise and rounding errors. We also present a simple condition that implies the backward instability of the barycentric formula used to implement extended interpolants. Our experiments show that extended interpolants mentioned in the literature satisfy this condition and, therefore, the formula used to implement them is not backward stable. Finally, we explain that the extrapolation step is a significant source of numerical instability for extended interpolants based on extrapolation.

What can we prove about the stability of extended interpolants

This appendix illustrates the difficulties in building a general and realistic stability theory for extended interpolants, a theory which would take into account the errors introduced by the current implementations of floating point arithmetic. We explain that the stability of extended interpolants is sensitive to the way we implement the extrapolation step, and that the accuracy of this step depends on the cancellation of the rounding errors. In fact, the errors incurred by extended interpolants can be enormous when we use the extrapolation formula proposed in [2] and [11], and compute $\tilde{\mathbf {y}}$ according to the following procedure:

• (a) If i is even, set the rounding mode upward and evaluate $\tilde{y}_i$ as in (6)–(8).

• (b) If i is odd, set the rounding mode downward and evaluate $\tilde{y}_i$ as in (6)–(8).

In this scenario the overall effect of rounding errors can be much larger than what one would expect from the already large Lebesgue constants, as illustrated in Figs. 13 and 14. In the plots corresponding to $\tilde{\mathbf {y}}$ evaluated as in (10)–(11) in these figures, $\tilde{\mathbf {y}}$ was obtained by matrix multiplication, with $a_{ij}$ and $b_{ij}$ computed in multiple precision and then rounded to double precision i.e., with $a_{ij}$ and $b_{ij}$ as accurate as possible.

Fig. 13

Extended Floater–Hormann. Log10 of the forward error for ${f}\left( t \right) = {\sin }\left( 20t \right)$ for $t \in [-1,1]$ and $n = 200$. By “$\tilde{y}$ by Taylor” we mean $\tilde{y}$ computed by Taylor series as in [2] and (6)–(8), and by “$\tilde{y}$ by matrix mult.” we mean $\tilde{y}$ computed by matrix multiplication of $\mathbf {y}$ by the matrices with entries $a_{ij}$ and $b_{ij}$, as in (10)–(11)

Fig. 14

Extended Floater–Hormann. Log10 of (forward error divided by the Lebesgue constant) for ${f}\left( t \right) = {\sin }\left( 20t \right)$, $t \in [-1,1]$ and $n = 200$

We emphasize that the choices of rounding modes in the steps (a) and (b) above are not frivolous. Their purpose is to help us understand what can be proved about the numerical stability of extended interpolants, so that we do not try to prove something that cannot be proved. It is unlikely that $\tilde{\mathbf {y}}$ will be evaluated as in the steps (a) and (b) when rounding to nearest. In this mode, there is a 50  % chance of rounding up in each flop and, under the questionable hypothesis of independence of the rounding errors, there would be a minuscule probability of $2^{-2\left( d+1 \right) d}$ of having all the intermediate results in the evaluation of $\tilde{\mathbf {y}}$ rounded up when rounding to nearest. Since the set of floating point numbers is finite, such a coincidence may be impossible. However, our experiments indicate that it is difficult to build a realistic theory on the effects of rounding errors on extended interpolants, because the rounding errors induced by our changes of rounding modes would be allowed by usual models of floating point arithmetic, with $\epsilon$ replaced by $2 \epsilon$. More precisely: when evaluating $\tilde{\mathbf {y}}$, with our choices of rounding modes, we monitored the relative errors

$$\begin{aligned} \left| \frac{{\mathrm {fl}}\left( x + y \right) - \left( x + y \right) }{\left( x + y \right) }\right| \quad \mathrm {and} \quad \left| \frac{{\mathrm {fl}}\left( x * y \right) - \left( x * y \right) }{\left( x * y \right) }\right| \end{aligned}$$

for each operation we performed, and found all of them to be smaller than $1.97 \epsilon$.

In other words, a stability theory for extended interpolants based on the usual models of floating point arithmetic would need to cover the changes of rounding modes in steps (a) and (b) above and, as a result, its predictions would be too pessimistic. Therefore, a realistic stability theory for extended interpolants will require additional hypothesis regarding the floating point arithmetic. By contrast, under the usual models of floating point arithmetic [10], we already have realistic theories bounding the rounding errors in terms of $\epsilon$, n and the Lebesgue constant for other barycentric interpolation schemes, as in [6, 9, 13–15].