# Chebyshev rational approximation with Remez algorithm: how do I tell how many digits are accurate?

I have implemented the Remez algorithm for Chebyshev rational approximation using a multiple precision numerics library (mpmath via SymPy). Wikipedia has suggested using $\max\{|z_i|\} - \min\{|z_i|\}$ as a stopping criterion. Here $z_i$ are the values of the error function $E(x) = f(x) - r_i(x)$ evaluated at the local extrema, that is, $z_i = E(x_i)$. Other sources have suggested stopping when the alternates $x_i$ stop changing, but Wikipedia's suggestion seems more reasonable to me, as the theory says that the perfect minimax rational approximation has equioscillating error, and $\max\{|z_i|\} - \min\{|z_i|\}$ is a measure of closeness to that.

Empirically, when I perform the computations with $N$ digits of precision, the convergence happens when $\max\{|z_i|\} - \min\{|z_i|\}$ goes under $10^N$ (that is, iterating past that doesn't produce any further improvement), so I've set that as my stopping condition. If anyone knows why this might be the case (or might not actually be the case) that would be great to know.

My question is, once the approximation has converged using $N$ digits of precision and $\max\{|z_i|\} - \min\{|z_i|\} < \epsilon$, how many digits of the coefficients in final approximation $r(x)$ are correct? Obviously my end goal is to get at least enough digits for some double floating point calculations.

My only goal post thus far has been a paper with published coefficients.

My question is twofold:

1. How can I tell how many digits of the coefficients in my computed approximation are accurate?

2. How can I put a lower bound on the number of accurate digits (preferably by increasing the number of digits of precision in the algorithm, $N$, or by changing the convergence conditions)?

I only care about approximating $e^{-x}$ on $[0, \infty)$, although general answers would be nice. I've been using the technique of the above paper (translate to $[-1, 1]$ via $z = c\frac{1 + x}{1 - x}$ with suitably chosen $c$ to spread out the errors, then translate back with $x = \frac{z - c}{z + c}$).

So far, I've thought of the following:

Suppose we have a polynomial approximation instead of a rational approximation (i.e., $E(x) = f(x) - p(x) < \epsilon$ is the error for the best minimax approximation). If we perturb a coefficient of $p(x)$ by $\delta$, and let $$E^* = f(x) - p^*(x),$$ where $$p^*(x) = a_0 + \cdots + (a_k + \delta)x^k + \cdots + a_nx^n,$$ then $$|E^*| = |E - \delta x^k| <= \epsilon + \delta$$ because $x\in [-1, 1].$ But I'm not sure if this helps me. I also don't know how to apply it to the rational case (how does the denominator affect things?), and finally, I don't know how the inverse transformation $\frac{z - c}{z + c}$ will affect things.

• I assume the $z_{i}$ represent error at the extrema of $approx(x) - f(x)$, and that you iterate until the spread between minimum and maximum of the $\left|z_{i}\right|$ is less than some predetermined bound? That is a very common terminating condition in the Remez algorithm. As a rule of thumb, I stop the leveling process when $-log(spread) > -3*log(\epsilon_{machine})$, so bound of $10^{-50}$ for IEEE double precision. Intermediate computations needs very high precision, as parts of the computation is ill-conditioned, I routinely use 1000 bits even when computing double-precision coefficients – njuffa Dec 6 '16 at 5:41
• @njuffa yes you are right about the $z_i$ (I've explained it in the question). Rules of thumb are helpful of course, but I would really like to know how to tell how many correct digits I have. If I even knew that, I could use a large number of digits of precision and check that I am correct enough, but knowing a formula for how many to use would be helpful as well. But as I stand now, even if I use 1000 digits, how can I know that I have enough correct digits in the result (beyond checking against published digits, which I have to trust are all correct)? – asmeurer Dec 6 '16 at 5:48
• It's also worth noting that mpmath automatically increases the internal precision automatically to give results of the desired precision. So, for instance, when I solve the system $E = f(x_i) + r(x_i) + (-1)^i\varepsilon$, I know that the answer is correct up to $N$ digits of precision (it automatically verifies the solution it finds against an appropriate tolerance for $N$). – asmeurer Dec 6 '16 at 5:53
• Oops, sorry, I overlooked the definition of $z_{i}$ even though I read the question twice. $max{z_{i}}$ tells you the approximation error assuming the evaluation of the approximation is done in the high-precision floating-point currently used. The rounding of the coefficients to the desired target format (and evaluation of the approximation in that format) introduces new error, and I optimize that by heuristic searches. The first part of this issue is addressed by the fpminimax function of the Sollya tool, the second part an open research problem. – njuffa Dec 6 '16 at 6:17

For practical purposes you can compute the approximation for some large $K >> N$, say $K = N + 100$ (I've used $K=1000$). Then, although you do not know how many digits of the $K$-digit approximation are accurate, you can be sure that at least $N$ of them are, so you can then compare how many digits of the $N$-digit approximation agree with the first $N$ digits of the $K$ approximation.
This is a little hand-wavy, because it relies on knowing that the algorithm has properly "converged". I've actually had to adjust the $10^N$ condition I said above to match the published digits in that paper, so this does matter. However, (as far as I can tell) even partial convergence should have the first digits correct, so if $K$ is large enough, then this shouldn't be an issue.