Remez algorithm

The Remez algorithm starts with the function $f$ to be approximated and a set $X$ of $n\; +\; 2$ sample points $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

:$b\_0\; +\; b\_1\; x\_i+\; ...\; +b\_n\; x\_i\; ^\; n\; +\; (-1)^\; i\; E\; =\; f(x\_i)$ (where $i=1,\; 2,\; ...\; n+2$),

:for the unknowns $b\_0,\; b\_1...b\_n$ and E.

• Use the $b\_i$ as coefficients to form a polynomial $P\_n$.
• Find the set $M$ of points of local maximum error $|P\_n(x)\; -\; f(x)|$.
• If the errors at every $m\; \backslash in\; M$ are of equal magnitude and alternate in sign, then $P\_n$ is the minimax approximation polynomial. If not, replace $X$ with $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.{{cite journal |doi=10.1145/321281.321282 |first=W. |last=Fraser |title=A Survey of Methods of Computing Minimax and Near-Minimax Polynomial Approximations for Functions of a Single Independent Variable |journal=J. ACM |volume=12 |pages=295–314 |year=1965 |issue=3 |s2cid=2736060 |doi-access=free }}

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

:$\backslash lVert\; f\; -\; L\_n(f)\backslash rVert\_\backslash infty\; \backslash le\; (1\; +\; \backslash lVert\; L\_n\backslash rVert\_\backslash infty)\; \backslash inf\_\{p\; \backslash in\; P\_n\}\; \backslash lVert\; f\; -\; p\backslash rVert$

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

:$\backslash lVert\; L\_n\backslash rVert\_\backslash infty\; =\; \backslash overline\{\backslash Lambda\}\_n(T)\; =\; \backslash max\_\{-1\; \backslash le\; x\; \backslash le\; 1\}\; \backslash lambda\_n(T;\; x),$

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

:$\backslash lambda\_n(T;\; x)\; =\; \backslash sum\_\{j\; =\; 1\}^\{n\; +\; 1\}\; \backslash left|\; l\_j(x)\; \backslash right|,\; \backslash quad\; l\_j(x)\; =\; \backslash prod\_\{\backslash stackrel\{i\; =\; 1\}\{i\; \backslash ne\; j\}\}^\{n\; +\; 1\}\; \backslash frac\{(x\; -\; t\_i)\}\{(t\_j\; -\; t\_i)\}.$

Theodore A. Kilgore,{{cite journal |doi=10.1016/0021-9045(78)90013-8 |first=T. A. |last=Kilgore |title=A characterization of the Lagrange interpolating projection with minimal Tchebycheff norm |journal=J. Approx. Theory |volume=24 |pages=273–288 |year=1978 |issue=4 |doi-access= }} Carl de Boor, and Allan Pinkus{{cite journal |doi=10.1016/0021-9045(78)90014-X |first1=C. |last1=de Boor |first2=A. |last2=Pinkus |title=Proof of the conjectures of Bernstein and Erdös concerning the optimal nodes for polynomial interpolation |journal=Journal of Approximation Theory |volume=24 |pages=289–303 |year=1978 |issue=4 |doi-access=free }} proved that there exists a unique t_i for each L_n, although not known explicitly for (ordinary) polynomials. Similarly, $\backslash underline\{\backslash Lambda\}\_n(T)\; =\; \backslash min\_\{-1\; \backslash le\; x\; \backslash le\; 1\}\; \backslash lambda\_n(T;\; x)$, and the optimality of a choice of nodes can be expressed as $\backslash overline\{\backslash Lambda\}\_n\; -\; \backslash underline\{\backslash Lambda\}\_n\; \backslash ge\; 0.$

For Chebyshev nodes, which provides a suboptimal, but analytically explicit choice, the asymptotic behavior is known as{{cite journal |first1=F. W. |last1=Luttmann |first2=T. J. |last2=Rivlin |title=Some numerical experiments in the theory of polynomial interpolation |journal=IBM J. Res. Dev. |volume=9 |pages=187–191 |year=1965 |issue=3 |doi= 10.1147/rd.93.0187}}

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

({{math|γ}} being the Euler–Mascheroni constant) with

: for $n\; \backslash ge\; 1,$

and upper bound{{cite book |first=T.J. |last=Rivlin |chapter=The lebesgue constants for polynomial interpolation |chapter-url=https://link.springer.com/chapter/10.1007/BFb0063594 |doi=10.1007/BFb0063594 |series=Lecture Notes in Mathematics |volume=399 |editor-last=Garnir |editor-first=H.G. |editor2-last=Unni |editor2-first=K.R. |editor3-last=Williamson |editor3-first=J.H. |title=Functional Analysis and its Applications |publisher=Springer |date=1974 |isbn=978-3-540-37827-3 |pages=422–437 }}

:$\backslash overline\{\backslash Lambda\}\_n(T)\; \backslash le\; \backslash frac\{2\}\{\backslash pi\}\; \backslash log(n\; +\; 1)\; +\; 1$

Lev Brutman{{cite journal |doi=10.1137/0715046 |first=L. |last=Brutman |title=On the Lebesgue Function for Polynomial Interpolation |journal=SIAM J. Numer. Anal. |volume=15 |pages=694–704 |year=1978 |issue=4 |bibcode=1978SJNA...15..694B }} obtained the bound for $n\; \backslash ge\; 3$, and $\backslash hat\{T\}$ being the zeros of the expanded Chebyshev polynomials:

:

Rüdiger Günttner{{cite journal |doi=10.1137/0717043 |first=R. |last=Günttner |title=Evaluation of Lebesgue Constants |journal=SIAM J. Numer. Anal. |volume=17 |pages=512–520 |year=1980 |issue=4 |bibcode=1980SJNA...17..512G }} obtained from a sharper estimate for $n\; \backslash ge\; 40$

:

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 $x\_0,\; x\_1,\; ...\; x\_\{n+1\}$, solve the linear system of n+2 equations

:$b\_0\; +\; b\_1\; x\_i+\; ...\; +b\_n\; x\_i\; ^\; n\; +\; (-1)\; ^\; i\; E\; =\; f(x\_i)$ (where $i=0,\; 1,\; ...\; n+1$),

:for the unknowns $b\_0,\; b\_1,\; ...b\_n$ and E.

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

Compute the standard n-th degree interpolant $p\_1(x)$ to $f(x)$ at the first n+1 nodes and also the standard n-th degree interpolant

$p\_2(x)$ to the ordinates $(-1)^i$

:$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 $0,\; ...,n$ and $O(n^2)$ arithmetic operations.

The polynomial $p\_2(x)$ has its i-th zero between $x\_\{i-1\}$ and $x\_i,\backslash \; i=1,\; ...,n$, and thus no further zeroes between $x\_n$ and $x\_\{n+1\}$: $p\_2(x\_n)$ and $p\_2(x\_\{n+1\})$ have the same sign $(-1)^n$.

The linear combination

$p(x)\; :=\; p\_1\; (x)\; -\; p\_2(x)\backslash !\backslash cdot\backslash !E$ is also a polynomial of degree n and

:$p(x\_i)\; =\; p\_1(x\_i)\; -\; p\_2(x\_i)\backslash !\backslash cdot\backslash !\; E\; \backslash \; =\; \backslash \; f(x\_i)\; -\; (-1)^i\; E,\backslash \; \backslash \; \backslash \; \backslash \; i\; =0,\; \backslash ldots,\; n.$

This is the same as the equation above for $i\; =\; 0,\; ...\; ,n$ and for any choice of E.

The same equation for i = n+1 is

:$p(x\_\{n+1\})\; \backslash \; =\; \backslash \; p\_1(x\_\{n+1\})\; -\; p\_2(x\_\{n+1\})\backslash !\backslash cdot\backslash !E\; \backslash \; =\; \backslash \; f(x\_\{n+1\})\; -\; (-1)^\{n+1\}\; E$ and needs special reasoning: solved for the variable E, it is the definition of E:

:$E\; \backslash \; :=\; \backslash \; \backslash 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 $p(x)\; \backslash equiv\; b\_0\; +\; b\_1x\; +\; \backslash ldots\; +\; b\_nx^n$ are always well-defined.

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

:$p(x\_i)\; -\; f(x\_i)\; \backslash \; =\; \backslash \; -(-1)^i\; E,\backslash \; \backslash \; i\; =\; 0,\; ...\; ,\; n\backslash !+\backslash !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 $\backslash tilde\; p(x)$, then the difference

$p(x)-\backslash tilde\; p(x)\; =\; (p(x)\; -\; f(x))\; -\; (\backslash tilde\; p(x)\; -\; f(x))$ would still be positive/negative at the n+2 nodes $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

$b\_0\; +\; b\_1x\; +\; ...\; +\; b\_nx^n$ to $p(x)$.

Step 3 improves upon the input nodes $x\_0,\; ...,\; x\_\{n+1\}$ and their errors $\backslash pm\; E$ as follows.

In each P-region, the current node $x\_i$ is replaced with the local maximizer $\backslash bar\{x\}\_i$ and in each N-region $x\_i$ is replaced with the local minimizer. (Expect $\backslash bar\{x\}\_0$ at A, the $\backslash bar\; \{x\}\_i$ near $x\_i$, and $\backslash 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 {{cite book |last1=Luenberger |first1=D.G. |last2=Ye |first2=Y. |chapter=Basic Descent Methods |chapter-url=https://link.springer.com/chapter/10.1007/978-0-387-74503-9_8 |title=Linear and Nonlinear Programming |publisher=Springer |edition=3rd |series=International Series in Operations Research & Management Science |volume=116 |date=2008 |isbn=978-0-387-74503-9 |pages=215–262 |doi=10.1007/978-0-387-74503-9_8}})

Let $z\_i\; :=\; p(\backslash bar\{x\}\_i)\; -\; f(\backslash bar\{x\}\_i)$. Each amplitude $|z\_i|$ is greater than or equal to E. The Theorem of de La Vallée Poussin and its proof also

apply to $z\_0,\; ...\; ,z\_\{n+1\}$ with $\backslash min\backslash \{|z\_i|\backslash \}\; \backslash geq\; E$ as the new

lower bound for the best error possible with polynomials of degree n.

Moreover, $\backslash max\backslash \{|z\_i|\backslash \}$ comes in handy as an obvious upper bound for that best possible error.

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

Some modifications of the algorithm are present on the literature.{{Citation |last1=Egidi |first1=Nadaniela |title=A New Remez-Type Algorithm for Best Polynomial Approximation |date=2020 |url=http://link.springer.com/10.1007/978-3-030-39081-5_7 |work=Numerical Computations: Theory and Algorithms |volume=11973 |pages=56–69 |editor-last=Sergeyev |editor-first=Yaroslav D. |place=Cham |publisher=Springer |doi=10.1007/978-3-030-39081-5_7 |isbn=978-3-030-39080-8 |last2=Fatone |first2=Lorella |last3=Misici |first3=Luciano |s2cid=211159177 |editor2-last=Kvasov |editor2-first=Dmitri E.}} These include:

• Replacing more than one sample point with the locations of nearby maximum absolute differences.{{Citation needed|date=March 2022}}
• Replacing all of the sample points with in a single iteration with the locations of all, alternating sign, maximum differences.{{cite journal |last1=Temes |first1=G.C. |last2=Barcilon |first2=V. |last3=Marshall |first3=F.C. |title=The optimization of bandlimited systems |journal=Proceedings of the IEEE |volume=61 |issue=2 |pages=196–234 |date=1973 |doi=10.1109/PROC.1973.9004 |issn=0018-9219}}
• 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.
• The Fraser-Hart variant, used to determine the best rational Chebyshev approximation.{{Cite journal |last=Dunham |first=Charles B. |date=1975 |title=Convergence of the Fraser-Hart algorithm for rational Chebyshev approximation |url=https://www.ams.org/mcom/1975-29-132/S0025-5718-1975-0388732-9/ |journal=Mathematics of Computation |language=en |volume=29 |issue=132 |pages=1078–1082 |doi=10.1090/S0025-5718-1975-0388732-9 |issn=0025-5718|doi-access=free }}

{{Portal|Mathematics}}

• {{annotated link|Hadamard's lemma}}
• {{annotated link|Laurent series}}
• {{annotated link|Padé approximant}}
• {{annotated link|Newton series}}
• {{annotated link|Approximation theory}}
• {{annotated link|Function approximation}}

{{Reflist}}

• [https://www.boost.org/doc/libs/1_47_0/libs/math/doc/sf_and_dist/html/math_toolkit/toolkit/internals2/minimax.html Minimax Approximations and the Remez Algorithm], background chapter in the Boost Math Tools documentation, with link to an implementation in C++
• [http://www.bores.com/courses/intro/filters/4_equi.htm Intro to DSP]
• {{MathWorld|urlname=RemezAlgorithm|title=Remez Algorithm|author1-link=Ronald Aarts|author=Aarts, Ronald M.|author2=Bond, Charles|author3=Mendelsohn, Phil|author4= Weisstein, Eric W.|name-list-style=amp}}

Category:Polynomials

Category:Approximation theory

Category:Numerical analysis