Editing Approximation theory (section)

==Remez's algorithm==
The [[Remez algorithm]] (sometimes spelled Remes) is used to produce an optimal polynomial ''P''(''x'') approximating a given function ''f''(''x'') over a given interval. It is an iterative algorithm that converges to a polynomial that has an error function with ''N''+2 level extrema. By the theorem above, that polynomial is optimal.

Remez's algorithm uses the fact that one can construct an ''N''<span style="padding-left:0.1em;">th</span>-degree polynomial that leads to level and alternating error values, given ''N''+2 test points.

Given ''N''+2 test points <math>x_1</math>, <math>x_2</math>, ... <math>x_{N+2}</math> (where <math>x_1</math> and <math>x_{N+2}</math> are presumably the end points of the interval of approximation), these equations need to be solved:

:<math>\begin{align}
          P(x_1) - f(x_1) &= +\varepsilon \\
          P(x_2) - f(x_2) &= -\varepsilon \\
          P(x_3) - f(x_3) &= +\varepsilon \\
                          &\ \ \vdots \\
  P(x_{N+2}) - f(x_{N+2}) &= \pm\varepsilon.
\end{align}</math>

The right-hand sides alternate in sign.

That is,

:<math>\begin{align}
  P_0 + P_1 x_1 + P_2 x_1^2 + P_3 x_1^3 + \dots + P_N x_1^N - f(x_1) &= +\varepsilon \\
  P_0 + P_1 x_2 + P_2 x_2^2 + P_3 x_2^3 + \dots + P_N x_2^N - f(x_2) &= -\varepsilon \\
                                                                     &\ \ \vdots
\end{align}</math>

Since <math>x_1</math>, ..., <math>x_{N+2}</math> were given, all of their powers are known, and <math>f(x_1)</math>, ..., <math>f(x_{N+2})</math> are also known. That means that the above equations are just ''N''+2 linear equations in the ''N''+2 variables <math>P_0</math>, <math>P_1</math>, ..., <math>P_N</math>, and <math>\varepsilon</math>. Given the test points <math>x_1</math>, ..., <math>x_{N+2}</math>, one can solve this system to get the polynomial ''P'' and the number <math>\varepsilon</math>.

The graph below shows an example of this, producing a fourth-degree polynomial approximating <math>e^x</math> over [−1, 1]. The test points were set at
−1, −0.7, −0.1, +0.4, +0.9, and 1. Those values are shown in green. The resultant value of <math>\varepsilon</math> is 4.43 × 10<sup>−4</sup>

[[Image:Remesdemo.png|thumb|center|300px|Error of the polynomial produced by the first step of Remez's algorithm, approximating e<sup>x</sup> over the interval [−1, 1]. Vertical divisions are 10<sup>−4</sup>.]]

The error graph does indeed take on the values <math>\pm \varepsilon</math> at the six test points, including the end points, but that those points are not extrema. If the four interior test points had been extrema (that is, the function ''P''(''x'')''f''(''x'') had maxima or minima there), the polynomial would be optimal.

The second step of Remez's algorithm consists of moving the test points to the approximate locations where the error function had its actual local maxima or minima. For example, one can tell from looking at the graph that the point at −0.1 should have been at about −0.28. The way to do this in the algorithm is to use a single round of [[Newton's method]]. Since one knows the first and second derivatives of {{nowrap|''P''(''x'') − ''f''(''x'')}}, one can calculate approximately how far a test point has to be moved so that the derivative will be zero.

Calculating the derivatives of a polynomial is straightforward. One must also be able to calculate the first and second derivatives of ''f''(''x''). Remez's algorithm requires an ability to calculate <math>f(x)\,</math>, <math>f'(x)\,</math>, and <math>f''(x)\,</math> to extremely high precision. The entire algorithm must be carried out to higher precision than the desired precision of the result.

After moving the test points, the linear equation part is repeated, getting a new polynomial, and Newton's method is used again to move the test points again. This sequence is continued until the result converges to the desired accuracy. The algorithm converges very rapidly. Convergence is quadratic for well-behaved functions—if the test points are within <math>10^{-15}</math> of the correct result, they will be approximately within <math>10^{-30}</math> of the correct result after the next round.

Remez's algorithm is typically started by choosing the extrema of the Chebyshev polynomial <math>T_{N+1}</math> as the initial points, since the final error function will be similar to that polynomial.