The Remez algorithm starts with the function ${\displaystyle f}$  to be approximated and a set ${\displaystyle X}$  of ${\displaystyle n+2}$  sample points ${\displaystyle x_{1},x_{2},...,x_{n+2}}$  in the approximation interval, usually the extrema of Chebyshev polynomial linearly mapped to the interval. The steps are:

• Solve the linear system of equations

${\displaystyle b_{0}+b_{1}x_{i}+...+b_{n}x_{i}^{n}+(-1)^{i}E=f(x_{i})}$  (where ${\displaystyle i=1,2,...n+2}$ ),

for the unknowns ${\displaystyle b_{0},b_{1}...b_{n}}$  and E.

• Use the ${\displaystyle b_{i}}$  as coefficients to form a polynomial ${\displaystyle P_{n}}$ .
• Find the set ${\displaystyle M}$  of points of local maximum error ${\displaystyle |P_{n}(x)-f(x)|}$ .
• If the errors at every ${\displaystyle m\in M}$  are of equal magnitude and alternate in sign, then ${\displaystyle P_{n}}$  is the minimax approximation polynomial. If not, replace ${\displaystyle X}$  with ${\displaystyle M}$  and repeat the steps above.

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.^[2]

Choice of initialization edit

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

${\displaystyle \lVert f-L_{n}(f)\rVert _{\infty }\leq (1+\lVert L_{n}\rVert _{\infty })\inf _{p\in P_{n}}\lVert f-p\rVert }$ 

with the norm or Lebesgue constant of the Lagrange interpolation operator L_n of the nodes (t₁, ..., t_n + 1) being

${\displaystyle \lVert L_{n}\rVert _{\infty }={\overline {\Lambda }}_{n}(T)=\max _{-1\leq x\leq 1}\lambda _{n}(T;x),}$ 

T being the zeros of the Chebyshev polynomials, and the Lebesgue functions being

${\displaystyle \lambda _{n}(T;x)=\sum _{j=1}^{n+1}\left|l_{j}(x)\right|,\quad l_{j}(x)=\prod _{\stackrel {i=1}{i\neq j}}^{n+1}{\frac {(x-t_{i})}{(t_{j}-t_{i})}}.}$ 

Theodore A. Kilgore,^[3] Carl de Boor, and Allan Pinkus^[4] proved that there exists a unique t_i for each L_n, although not known explicitly for (ordinary) polynomials. Similarly, ${\displaystyle {\underline {\Lambda }}_{n}(T)=\min _{-1\leq x\leq 1}\lambda _{n}(T;x)}$ , and the optimality of a choice of nodes can be expressed as ${\displaystyle {\overline {\Lambda }}_{n}-{\underline {\Lambda }}_{n}\geq 0.}$ 

For Chebyshev nodes, which provides a suboptimal, but analytically explicit choice, the asymptotic behavior is known as^[5]

${\displaystyle {\overline {\Lambda }}_{n}(T)={\frac {2}{\pi }}\log(n+1)+{\frac {2}{\pi }}\left(\gamma +\log {\frac {8}{\pi }}\right)+\alpha _{n+1}}$ 

(γ being the Euler–Mascheroni constant) with

${\displaystyle 0<\alpha _{n}<{\frac {\pi }{72n^{2}}}}$  for ${\displaystyle n\geq 1,}$ 

and upper bound^[6]

${\displaystyle {\overline {\Lambda }}_{n}(T)\leq {\frac {2}{\pi }}\log(n+1)+1}$ 

Lev Brutman^[7] obtained the bound for ${\displaystyle n\geq 3}$ , and ${\displaystyle {\hat {T}}}$  being the zeros of the expanded Chebyshev polynomials:

${\displaystyle {\overline {\Lambda }}_{n}({\hat {T}})-{\underline {\Lambda }}_{n}({\hat {T}})<{\overline {\Lambda }}_{3}-{\frac {1}{6}}\cot {\frac {\pi }{8}}+{\frac {\pi }{64}}{\frac {1}{\sin ^{2}(3\pi /16)}}-{\frac {2}{\pi }}(\gamma -\log \pi )\approx 0.201.}$ 

Rüdiger Günttner^[8] obtained from a sharper estimate for ${\displaystyle n\geq 40}$ 

${\displaystyle {\overline {\Lambda }}_{n}({\hat {T}})-{\underline {\Lambda }}_{n}({\hat {T}})<0.0196.}$ 

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 ${\displaystyle x_{0},x_{1},...x_{n+1}}$ , solve the linear system of n+2 equations

${\displaystyle b_{0}+b_{1}x_{i}+...+b_{n}x_{i}^{n}+(-1)^{i}E=f(x_{i})}$  (where ${\displaystyle i=0,1,...n+1}$ ),

for the unknowns ${\displaystyle b_{0},b_{1},...b_{n}}$  and E.

It should be clear that ${\displaystyle (-1)^{i}E}$  in this equation makes sense only if the nodes ${\displaystyle x_{0},...,x_{n+1}}$  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 ${\displaystyle O(n^{2})}$  arithmetic operations while a standard solver from the library would take ${\displaystyle O(n^{3})}$  operations. Here is the simple proof:

Compute the standard n-th degree interpolant ${\displaystyle p_{1}(x)}$  to ${\displaystyle f(x)}$  at the first n+1 nodes and also the standard n-th degree interpolant ${\displaystyle p_{2}(x)}$  to the ordinates ${\displaystyle (-1)^{i}}$ 

${\displaystyle p_{1}(x_{i})=f(x_{i}),p_{2}(x_{i})=(-1)^{i},i=0,...,n.}$ 

To this end, use each time Newton's interpolation formula with the divided differences of order ${\displaystyle 0,...,n}$  and ${\displaystyle O(n^{2})}$  arithmetic operations.

The polynomial ${\displaystyle p_{2}(x)}$  has its i-th zero between ${\displaystyle x_{i-1}}$  and ${\displaystyle x_{i},\ i=1,...,n}$ , and thus no further zeroes between ${\displaystyle x_{n}}$  and ${\displaystyle x_{n+1}}$ : ${\displaystyle p_{2}(x_{n})}$  and ${\displaystyle p_{2}(x_{n+1})}$  have the same sign ${\displaystyle (-1)^{n}}$ .

The linear combination ${\displaystyle p(x):=p_{1}(x)-p_{2}(x)\!\cdot \!E}$  is also a polynomial of degree n and

${\displaystyle p(x_{i})=p_{1}(x_{i})-p_{2}(x_{i})\!\cdot \!E\ =\ f(x_{i})-(-1)^{i}E,\ \ \ \ i=0,\ldots ,n.}$ 

This is the same as the equation above for ${\displaystyle i=0,...,n}$  and for any choice of E. The same equation for i = n+1 is

${\displaystyle p(x_{n+1})\ =\ p_{1}(x_{n+1})-p_{2}(x_{n+1})\!\cdot \!E\ =\ f(x_{n+1})-(-1)^{n+1}E}$  and needs special reasoning: solved for the variable E, it is the definition of E:

${\displaystyle E\ :=\ {\frac {p_{1}(x_{n+1})-f(x_{n+1})}{p_{2}(x_{n+1})+(-1)^{n}}}.}$ 

As mentioned above, the two terms in the denominator have same sign: E and thus ${\displaystyle p(x)\equiv b_{0}+b_{1}x+\ldots +b_{n}x^{n}}$  are always well-defined.

The error at the given n+2 ordered nodes is positive and negative in turn because

${\displaystyle p(x_{i})-f(x_{i})\ =\ -(-1)^{i}E,\ \ i=0,...,n\!+\!1.}$ 

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 ${\displaystyle {\tilde {p}}(x)}$ , then the difference ${\displaystyle p(x)-{\tilde {p}}(x)=(p(x)-f(x))-({\tilde {p}}(x)-f(x))}$  would still be positive/negative at the n+2 nodes ${\displaystyle x_{i}}$  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 ${\displaystyle b_{0}+b_{1}x+...+b_{n}x^{n}}$  to ${\displaystyle p(x)}$ .

Step 3 improves upon the input nodes ${\displaystyle x_{0},...,x_{n+1}}$  and their errors ${\displaystyle \pm E}$  as follows.

In each P-region, the current node ${\displaystyle x_{i}}$  is replaced with the local maximizer ${\displaystyle {\bar {x}}_{i}}$  and in each N-region ${\displaystyle x_{i}}$  is replaced with the local minimizer. (Expect ${\displaystyle {\bar {x}}_{0}}$  at A, the ${\displaystyle {\bar {x}}_{i}}$  near ${\displaystyle x_{i}}$ , and ${\displaystyle {\bar {x}}_{n+1}}$  at B.) No high precision is required here, the standard line search with a couple of quadratic fits should suffice. (See ^[9])

Let ${\displaystyle z_{i}:=p({\bar {x}}_{i})-f({\bar {x}}_{i})}$ . Each amplitude ${\displaystyle |z_{i}|}$  is greater than or equal to E. The Theorem of de La Vallée Poussin and its proof also apply to ${\displaystyle z_{0},...,z_{n+1}}$  with ${\displaystyle \min\{|z_{i}|\}\geq E}$  as the new lower bound for the best error possible with polynomials of degree n.

Moreover, ${\displaystyle \max\{|z_{i}|\}}$  comes in handy as an obvious upper bound for that best possible error.

Step 4: With ${\displaystyle \min \,\{|z_{i}|\}}$  and ${\displaystyle \max \,\{|z_{i}|\}}$  as lower and upper bound for the best possible approximation error, one has a reliable stopping criterion: repeat the steps until ${\displaystyle \max\{|z_{i}|\}-\min\{|z_{i}|\}}$  is sufficiently small or no longer decreases. These bounds indicate the progress.

Some modifications of the algorithm are present on the literature.^[10] These include:

• Replacing more than one sample point with the locations of nearby maximum absolute differences.^{[citation needed]}
• Replacing all of the sample points with in a single iteration with the locations of all, alternating sign, maximum differences.^[11]
• Using the relative error 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;
• Including zero-error point constraints.^[11]
• The Fraser-Hart variant, used to determine the best rational Chebyshev approximation.^[12]

• Approximation theory

• Minimax Approximations and the Remez Algorithm, background chapter in the Boost Math Tools documentation, with link to an implementation in C++
• Intro to DSP
• Aarts, Ronald M.; Bond, Charles; Mendelsohn, Phil & Weisstein, Eric W. "Remez Algorithm". MathWorld.