Remez algorithm

The Remez algorithm or Remez exchange algorithm, published by Evgeny Yakovlevich Remez in 1934, is an iterative algorithm used to find simple approximations to functions, specifically, approximations by functions in a Chebyshev space that are the best in the uniform norm "L"∞ sense.
A typical example of a Chebyshev space is the subspace of Chebyshev polynomials of order "n" in the space of real continuous functions on an interval, "C"["a", "b"]. The polynomial of best approximation within a given subspace is defined to be the one that minimizes the maximum absolute difference between the polynomial and the function. In this case, the form of the solution is precised by the equioscillation theorem.
Procedure.
The Remez algorithm starts with the function formula_1 to be approximated and a set formula_2 of formula_3 sample points formula_4 in the approximation interval, usually the extrema of Chebyshev polynomial linearly mapped to the interval. The steps are:
The result is called the polynomial of best approximation or the minimax approximation algorithm.
A review of technicalities in implementing the Remez algorithm is given by W. Fraser.
On the choice of initialization.
The Chebyshev nodes are a common choice for the initial approximation because of their role in the theory of polynomial interpolation. For the initialization of the optimization problem for function "f" by the Lagrange interpolant "L"n("f"), it can be shown that this initial approximation is bounded by
with the norm or Lebesgue constant of the Lagrange interpolation operator "L""n" of the nodes ("t"1, ..., "t""n" + 1) being
"T" being the zeros of the Chebyshev polynomials, and the Lebesgue functions being
Theodore A. Kilgore, Carl de Boor, and Allan Pinkus proved that there exists a unique "t""i" for each "L""n", although not known explicitly for (ordinary) polynomials. Similarly, formula_19, and the optimality of a choice of nodes can be expressed as formula_20
For Chebyshev nodes, which provides a suboptimal, but analytically explicit choice, the asymptotic behavior is known as
( being the Euler–Mascheroni constant) with
and upper bound
Lev Brutman obtained the bound for formula_25, and formula_26 being the zeros of the expanded Chebyshev polynomials:
Rüdiger Günttner obtained from a sharper estimate for formula_28
Detailed discussion.
This section provides more information on the steps outlined above. In this section, the index "i" runs from 0 to "n"+1.
Step 1: Given formula_30, solve the linear system of "n"+2 equations
It should be clear that formula_34 in this equation makes sense only if the nodes formula_35 are "ordered", either strictly increasing or strictly decreasing. Then this linear system has a unique solution. (As is well known, not every linear system has a solution.) Also, the solution can be obtained with only formula_36 arithmetic operations while a standard solver from the library would take formula_37 operations. Here is the simple proof:
Compute the standard "n"-th degree interpolant formula_38 to formula_39 at the first "n"+1 nodes and also the standard "n"-th degree interpolant
formula_40 to the ordinates formula_41
To this end, use each time Newton's interpolation formula with the divided
differences of order formula_43 and formula_36 arithmetic operations.
The polynomial formula_40 has its "i"-th zero between formula_46 and formula_47, and thus no further zeroes between formula_48 and formula_49: formula_50 and formula_51 have the same sign formula_52.
The linear combination
formula_53 is also a polynomial of degree "n" and
This is the same as the equation above for formula_55 and for any choice of "E".
The same equation for "i" = "n"+1 is
As mentioned above, the two terms in the denominator have same sign:
"E" and thus formula_58 are always well-defined.
The error at the given "n"+2 ordered nodes is positive and negative in turn because
The theorem of "de La Vallée Poussin" states that under this condition no polynomial of degree "n" exists with error less than "E". Indeed, if such a polynomial existed, call it formula_60, then the difference
formula_61 would still be positive/negative at the "n"+2 nodes formula_62 and therefore have at least "n"+1 zeros which is impossible for a polynomial of degree "n".
Thus, this "E" is a lower bound for the minimum error which can be achieved with polynomials of degree "n".
Step 2 changes the notation from
formula_63 to formula_64.
Step 3 improves upon the input nodes formula_65 and their errors formula_66 as follows.
In each P-region, the current node formula_62 is replaced with the local maximizer formula_68 and in each N-region formula_62 is replaced with the local minimizer. (Expect formula_70 at "A", the formula_71 near formula_62, and formula_73 at "B".) No high precision is required here,
the standard "line search" with a couple of "quadratic fits" should suffice. (See )
Let formula_74. Each amplitude formula_75 is greater than or equal to "E". The Theorem of "de La Vallée Poussin" and its proof also
apply to formula_76 with formula_77 as the new
lower bound for the best error possible with polynomials of degree "n".
Moreover, formula_78 comes in handy as an obvious upper bound for that best possible error.
Step 4: With formula_79 and formula_80 as lower and upper bound for the best possible approximation error, one has a reliable stopping criterion: repeat the steps until formula_81 is sufficiently small or no longer decreases. These bounds indicate the progress.
Variants.
Sometimes more than one sample point is replaced at the same time with the locations of nearby maximum absolute differences.
Sometimes all of the sample points are replaced in a single iteration with the locations of all, alternating sign, maximum differences.
Sometimes relative error is used to measure the difference between the approximation and the function, especially if the approximation will be used to compute the function on a computer which uses floating point arithmetic.
Sometimes zero-error point constraints are included in a Modified Remez Exchange Algorithm.