remez.qbk

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

[section:remez The Remez Method]

The [@http://en.wikipedia.org/wiki/Remez_algorithm Remez algorithm]
is a methodology for locating the minimax rational approximation
to a function.  This short article gives a brief overview of the method, but
it should not be regarded as a thorough theoretical treatment, for that you
should consult your favorite textbook.

Imagine that you want to approximate some function f(x) by way of a rational
function R(x), where R(x) may be either a polynomial P(x) or a ratio of two
polynomials P(x)/Q(x) (a rational function).  Initially we'll concentrate on the 
polynomial case, as it's by far the easier to deal with, later we'll extend 
to the full rational function case.  

We want to find the "best" rational approximation, where
"best" is defined to be the approximation that has the least deviation
from f(x).  We can measure the deviation by way of an error function:

E[sub abs](x) = f(x) - R(x)

which is expressed in terms of absolute error, but we can equally use
relative error:

E[sub rel](x) = (f(x) - R(x)) / |f(x)|

And indeed in general we can scale the error function in any way we want, it
makes no difference to the maths, although the two forms above cover almost
every practical case that you're likely to encounter.

The minimax rational function R(x) is then defined to be the function that
yields the smallest maximal value of the error function.  Chebyshev showed
that there is a unique minimax solution for R(x) that has the following
properties:

* If R(x) is a polynomial of degree N, then there are N+2 unknowns:
the N+1 coefficients of the polynomial, and maximal value of the error
function.
* The error function has N+1 roots, and N+2 extrema (minima and maxima).
* The extrema alternate in sign, and all have the same magnitude.

That means that if we know the location of the extrema of the error function
then we can write N+2 simultaneous equations:

R(x[sub i]) + (-1)[super i]E = f(x[sub i])

where E is the maximal error term, and x[sub i] are the abscissa values of the
N+2 extrema of the error function.  It is then trivial to solve the simultaneous
equations to obtain the polynomial coefficients and the error term.

['Unfortunately we don't know where the extrema of the error function are located!]

[h4 The Remez Method]

The Remez method is an iterative technique which, given a broad range of
assumptions, will converge on the extrema of the error function, and therefore
the minimax solution.

In the following discussion we'll use a concrete example to illustrate
the Remez method: an approximation to the function e[super x][space] over
the range \[-1, 1\].

Before we can begin the Remez method, we must obtain an initial value
for the location of the extrema of the error function.  We could "guess"
these, but a much closer first approximation can be obtained by first  
constructing an interpolated polynomial approximation to f(x).

In order to obtain the N+1 coefficients of the interpolated polynomial
we need N+1 points (x[sub 0]...x[sub N]): with our interpolated form 
passing through each of those points
that yields N+1 simultaneous equations:

f(x[sub i]) = P(x[sub i]) = c[sub 0] + c[sub 1]x[sub i] ... + c[sub N]x[sub i][super N]

Which can be solved for the coefficients c[sub 0]...c[sub N] in P(x).

Obviously this is not a minimax solution, indeed our only guarantee is that f(x) and 
P(x) touch at N+1 locations, away from those points the error may be arbitrarily
large.  However, we would clearly like this initial approximation to be as close to
f(x) as possible, and it turns out that using the zeros of an orthogonal polynomial
as the initial interpolation points is a good choice.  In our example we'll use the 
zeros of a Chebyshev polynomial as these are particularly easy to calculate, 
interpolating for a polynomial of degree 4, and measuring /relative error/
we get the following error function:

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

Which has a peak relative error of 1.2x10[super -3].

While this is a pretty good approximation already, judging by the 
shape of the error function we can clearly do better.  Before starting
on the Remez method propper, we have one more step to perform: locate
all the extrema of the error function, and store
these locations as our initial ['Chebyshev control points].

[note
In the simple case of a polynomial approximation, by interpolating through
the roots of a Chebyshev polynomial we have in fact created a ['Chebyshev
approximation] to the function: in terms of /absolute error/
this is the best a priori choice for the interpolated form we can
achieve, and typically is very close to the minimax solution.

However, if we want to optimise for /relative error/, or if the approximation
is a rational function, then the initial Chebyshev solution can be quite far
from the ideal minimax solution.  

A more technical discussion of the theory involved can be found in this
[@http://math.fullerton.edu/mathews/n2003/ChebyshevPolyMod.html online course].]

[h4 Remez Step 1]

The first step in the Remez method, given our current set of
N+2 Chebyshev control points x[sub i], is to solve the N+2 simultaneous
equations:

P(x[sub i]) + (-1)[super i]E = f(x[sub i])

To obtain the error term E, and the coefficients of the polynomial P(x).

This gives us a new approximation to f(x) that has the same error /E/ at
each of the control points, and whose error function ['alternates in sign]
at the control points.  This is still not necessarily the minimax 
solution though: since the control points may not be at the extrema of the error
function.  After this first step here's what our approximation's error
function looks like:

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

Clearly this is still not the minimax solution since the control points
are not located at the extrema, but the maximum relative error has now
dropped to 5.6x10[super -4].

[h4 Remez Step 2]

The second step is to locate the extrema of the new approximation, which we do 
in two stages:  first, since the error function changes sign at each
control point, we must have N+1 roots of the error function located between
each pair of N+2 control points.  Once these roots are found by standard root finding 
techniques, we know that N extrema are bracketed between each pair of
roots, plus two more between the endpoints of the range and the first and last roots.
The N+2 extrema can then be found using standard function minimisation techniques.

We now have a choice: multi-point exchange, or single point exchange.

In single point exchange, we move the control point nearest to the largest extrema to
the absissa value of the extrema.

In multi-point exchange we swap all the current control points, for the locations
of the extrema.

In our example we perform multi-point exchange.

[h4 Iteration]

The Remez method then performs steps 1 and 2 above iteratively until the control
points are located at the extrema of the error function: this is then
the minimax solution.

For our current example, two more iterations converges on a minimax
solution with a peak relative error of
5x10[super -4] and an error function that looks like:

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

[h4 Rational Approximations]

If we wish to extend the Remez method to a rational approximation of the form

f(x) = R(x) = P(x) / Q(x)

where P(x) and Q(x) are polynomials, then we proceed as before, except that now
we have N+M+2 unknowns if P(x) is of order N and Q(x) is of order M.  This assumes
that Q(x) is normalised so that it's leading coefficient is 1, giving
N+M+1 polynomial coefficients in total, plus the error term E.

The simultaneous equations to be solved are now:

P(x[sub i]) / Q(x[sub i]) + (-1)[super i]E = f(x[sub i])

Evaluated at the N+M+2 control points x[sub i].

Unfortunately these equations are non-linear in the error term E: we can only
solve them if we know E, and yet E is one of the unknowns!

The method usually adopted to solve these equations is an iterative one: we guess the
value of E, solve the equations to obtain a new value for E (as well as the polynomial
coefficients), then use the new value of E as the next guess.  The method is
repeated until E converges on a stable value.

These complications extend the running time required for the development