remez.qbk

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

of rational approximations quite considerably. It is often desirable
to obtain a rational rather than polynomial approximation none the less:
rational approximations will often match more difficult to approximate
functions, to greater accuracy, and with greater efficiency, than their
polynomial alternatives.  For example, if we takes our previous example
of an approximation to e[super x], we obtained 5x10[super -4] accuracy
with an order 4 polynomial.  If we move two of the unknowns into the denominator
to give a pair of order 2 polynomials, and re-minimise, then the peak relative error drops
to 8.7x10[super -5].  That's a 5 fold increase in accuracy, for the same number 
of terms overall.

[h4 Practical Considerations]

Most treatises on approximation theory stop at this point.  However, from
a practical point of view, most of the work involves finding the right
approximating form, and then persuading the Remez method to converge
on a solution.

So far we have used a direct approximation:

f(x) = R(x)

But this will converge to a useful approximation only if f(x) is smooth.  In
addition round-off errors when evaluating the rational form mean that this
will never get closer than within a few epsilon of machine precision.  
Therefore this form of direct approximation is often reserved for situations
where we want efficiency, rather than accuracy.

The first step in improving the situation is generally to split f(x) into
a dominant part that we can compute accurately by another method, and a 
slowly changing remainder which can be approximated by a rational approximation.
We might be tempted to write:

f(x) = g(x) + R(x)

where g(x) is the dominant part of f(x), but if f(x)\/g(x) is approximately
constant over the interval of interest then:

f(x) = g(x)(c + R(x))

Will yield a much better solution: here /c/ is a constant that is the approximate
value of f(x)\/g(x) and R(x) is typically tiny compared to /c/.  In this situation
if R(x) is optimised for absolute error, then as long as its error is small compared
to the constant /c/, that error will effectively get wiped out when R(x) is added to
/c/.

The difficult part is obviously finding the right g(x) to extract from your
function: often the asymptotic behaviour of the function will give a clue, so
for example the function __erfc becomes proportional to 
e[super -x[super 2]]\/x as x becomes large.  Therefore using:

erfc(z) = (C + R(x)) e[super -x[super 2]]/x

as the approximating form seems like an obvious thing to try, and does indeed
yield a useful approximation.

However, the difficulty then becomes one of converging the minimax solution.
Unfortunately, it is known that for some functions the Remez method can lead
to divergent behaviour, even when the initial starting approximation is quite good.
Furthermore, it is not uncommon for the solution obtained in the first Remez step
above to be a bad one: the equations to be solved are generally "stiff", often
very close to being singular, and assuming a solution is found at all, round-off
errors and a rapidly changing error function, can lead to a situation where the
error function does not in fact change sign at each control point as required.
If this occurs, it is fatal to the Remez method.  It is also possible to
obtain solutions that are perfectly valid mathematically, but which are
quite useless computationally: either because there is an unavoidable amount
of roundoff error in the computation of the rational function, or because
the denominator has one or more roots over the interval of the approximation.
In the latter case while the approximation may have the correct limiting value at
the roots, the approximation is nonetheless useless.

Assuming that the approximation does not have any fatal errors, and that the only
issue is converging adequately on the minimax solution, the aim is to
get as close as possible to the minimax solution before beginning the Remez method.
Using the zeros of a Chebyshev polynomial for the initial interpolation is a 
good start, but may not be ideal when dealing with relative errors and\/or
rational (rather than polynomial) approximations.  One approach is to skew
the initial interpolation points to one end: for example if we raise the
roots of the Chebyshev polynomial to a positive power greater than 1 
then the roots will be skewed towards the middle of the \[-1,1\] interval, 
while a positive power less than one
will skew them towards either end.  More usefully, if we initially rescale the
points over \[0,1\] and then raise to a positive power, we can skew them to the left 
or right.  Returning to our example of e[super x][space] over \[-1,1\], the initial
interpolated form was some way from the minimax solution:

[$../graphs/remez-2.png]

However, if we first skew the interpolation points to the left (rescale them
to \[0, 1\], raise to the power 1.3, and then rescale back to \[-1,1\]) we
reduce the error from 1.3x10[super -3][space]to 6x10[super -4]:

[$../graphs/remez-5.png]

It's clearly still not ideal, but it is only a few percent away from
our desired minimax solution (5x10[super -4]).

[h4 Remez Method Checklist]

The following lists some of the things to check if the Remez method goes wrong, 
it is by no means an exhaustive list, but is provided in the hopes that it will
prove useful.

* Is the function smooth enough?  Can it be better separated into
a rapidly changing part, and an asymptotic part?
* Does the function being approximated have any "blips" in it?  Check
for problems as the function changes computation method, or
if a root, or an infinity has been divided out.  The telltale
sign is if there is a narrow region where the Remez method will
not converge.
* Check you have enough accuracy in your calculations: remember that
the Remez method works on the difference between the approximation
and the function being approximated: so you must have more digits of
precision available than the precision of the approximation
being constructed.  So for example at double precision, you
shouldn't expect to be able to get better than a float precision
approximation.
* Try skewing the initial interpolated approximation to minimise the
error before you begin the Remez steps.
* If the approximation won't converge or is ill-conditioned from one starting
location, try starting from a different location.
* If a rational function won't converge, one can minimise a polynomial
(which presents no problems), then rotate one term from the numerator to
the denominator and minimise again.  In theory one can continue moving
terms one at a time from numerator to denominator, and then re-minimising, 
retaining the last set of control points at each stage.
* Try using a smaller interval.  It may also be possible to optimise over
one (small) interval, rescale the control points over a larger interval,
and then re-minimise.
* Keep absissa values small: use a change of variable to keep the abscissa
over, say \[0, b\], for some smallish value /b/.

[h4 References]

The original references for the Remez Method and it's extension
to rational functions are unfortunately in Russian:

Remez, E.Ya., ['Fundamentals of numerical methods for Chebyshev approximations], 
"Naukova Dumka", Kiev, 1969.

Remez, E.Ya., Gavrilyuk, V.T., ['Computer development of certain approaches 
to the approximate construction of solutions of Chebyshev problems 
nonlinearly depending on parameters], Ukr. Mat. Zh. 12 (1960), 324-338.

Gavrilyuk, V.T., ['Generalization of the first polynomial algorithm of 
E.Ya.Remez for the problem of constructing rational-fractional 
Chebyshev approximations], Ukr. Mat. Zh. 16 (1961), 575-585.

Some English language sources include:

Fraser, W., Hart, J.F., ['On the computation of rational approximations 
to continuous functions], Comm. of the ACM 5 (1962), 401-403, 414.

Ralston, A., ['Rational Chebyshev approximation by Remes' algorithms], 
Numer.Math. 7 (1965), no. 4, 322-330.

A. Ralston, ['Rational Chebyshev approximation, Mathematical 
Methods for Digital Computers v. 2] (Ralston A., Wilf H., eds.), 
Wiley, New York, 1967, pp. 264-284.

Hart, J.F. e.a., ['Computer approximations], Wiley, New York a.o., 1968.

Cody, W.J., Fraser, W., Hart, J.F., ['Rational Chebyshev approximation 
using linear equations], Numer.Math. 12 (1968), 242-251.

Cody, W.J., ['A survey of practical rational and polynomial 
approximation of functions], SIAM Review 12 (1970), no. 3, 400-423.

Barrar, R.B., Loeb, H.J., ['On the Remez algorithm for non-linear 
families], Numer.Math. 15 (1970), 382-391.

Dunham, Ch.B., ['Convergence of the Fraser-Hart algorithm for rational 
Chebyshev approximation], Math. Comp. 29 (1975), no. 132, 1078-1082.

G. L. Litvinov, ['Approximate construction of rational
approximations and the effect of error autocorrection],
Russian Journal of Mathematical Physics, vol.1, No. 3, 1994.

[endsect][/section:remez The Remez Method]

[/ 
  Copyright 2006 John Maddock and Paul A. Bristow.
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]