Coefficients of a Polynomial Approximation in Minimax Sense

Question

I am working on a Uniform Random Number Generator using a IEEE paper, and I got stuck with the coefficients for a Piecewise Polynomial Approximation using Horner's Rule :

$$ y = ((C_d x + C_{d-1})x + \ldots)x + C_0 $$

where $x$ is the input, $d$ is the polynomial degree, and $C_i$ are the polynomial coefficients.

According to the paper the polynomial coefficients are found in a mini-max sense, that minimizes the maximum absolute error.

Can anyone help with how I can generate these coefficients?

ref:<Link> http://www.ee.usyd.edu.au/people/philip.leong/UserFiles/File/papers/bm_tc06.pdf

Answer 1

An iterative method for finding the best "maximum norm" approximation by polynomial of degree at most $d$ to a given smooth function $f(x)$ on a bounded interval $I = [a,b]$ was proposed by Evgeny Yakovlevich Remez in 1934, and has come to be known as the Remez exchange algorithm.

The underlying idea is a characterization of the "mini-max" polynomial approximation called the equioscillation theorem attributed to Chebyshev:

Among all polynomials of degree $d$ or less, the polynomial $p(x)$ minimizes the max-norm of the error $||f - p||_\infty$ on $I = [a,b]$ if and only if there exist $d+2$ points $x_0 \lt x_1 \lt \ldots \lt x_{d+1}$ in $I$ such that the successive errors $f(x_i) - p(x_i)$ are equal in absolute value to the maximum $||f - p||_\infty$ and alternate in signs: $$ f(x_i) - p(x_i) = p(x_{i+1}) - f(x_{i+1}) $$ for $i=0,1,\ldots,d$.

Various implementations of the Remez exchange algorithm are available, and if you are specifically interested in the Matlab environment, see this package by Sherif Tawfik (2005). More about a "built-in" Matlab implementation is discussed under this DSP.SE Question.

The specific problem of approximating $f(x) = \sqrt x$ on $[1,4]$ using a first degree polynomial can be solved almost by inspection. Because of the concavity of this function, the best fit line will have two of its equioscillation points at the ends of the interval, and must necessarily be parallel to the secant line passing through the endpoints $(1,1)$ and $(4,2)$, namely $y = (x+2)/3$.

The only computation we need here is to find where the maximum of:

$$\sqrt{x} - \frac{1}{3}(x + 2) $$

occurs (the critical point) and split that "absolute error" difference with the endpoints.

A standard calculus approach locates the critical point at $x_1 = 9/4$, where the height of curve $\sqrt{x}$ above the secant line is $1/12$. This is illustrated in the graph below by the dotted vertical red line:

Fig. 1 The graph of $y=\sqrt{x}$ for $1\le x\le 4$ and chord through endpoints

Together with endpoints $x_0=1$ and $x_2=4$, we achieve the equioscillation criteria with:

$$ p(x) = \frac{1}{3}(x + 2) + \frac{1}{24} $$

The interested Reader may wish to simplify this polynomial further.