Remez exchange for a minimax polynomial

What it fills

• Polynomial approximation

  Given a target function, a domain and an error $\varepsilon$, return a polynomial of definite parity, bounded on $[-1,1]$, that is $\varepsilon$-close to the target on that domain, with an explicit degree.

When it applies

Remez is a classical approximation-theory algorithm, not a quantum one. The citation here is the QSP paper that employs it, alongside a Fourier-Chebyshev expansion, and should not be read as its origin. In that use the coefficients are solved from a reference set of degree-plus-two sampled points, which is then adjusted.

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

given  the target function f, a degree d, and the interval

choose a reference set of d + 2 sample points on the interval

repeat:
    solve for the degree-d polynomial P and the equioscillation error E
        that fit f at the d + 2 reference points
    move the reference set to the extrema of  f - P
until the reference set stops moving

return P

# the genuine minimax polynomial of degree d, where truncating a Chebyshev
# series is only near-optimal -- a tighter polynomial of the same degree
# handed on to the phase-factor stage

# classical approximation theory, not a quantum algorithm: the citation on
# this record is the QSP paper that employs it, not its origin

Cost, as the source states it

Dong, Meng, Whaley and Lin state no complexity for the Remez exchange itself. Their stated benefit is comparative: for approximating $1/x$, the minimax polynomial reaches the same accuracy at a degree smaller by a factor of 2–3 than the Fourier–Chebyshev truncation — their Table III has degrees 303–1519 (odd parity) and 280–1400 (even) against 759–4035 for the truncation, at $\varepsilon_0 = 10^{-14}$ over $\kappa = 10$–$50$. Stated in the full text; the abstract carries no formula.

Implementations

• QSPPACK's Remez solver (Solvers/JuliaSolver/Remez.ipynb)

  • Efficient phase-factor evaluation in quantum signal processing

  QSPPACK is the phase-factor toolbox released by authors of the paper this record's own citation names. Its readme lists two kinds of phase-factor solvers and then adds a third component: "The package also contains an implementation of the Remez algorithm for finding polynomial approximation." The notebook implementing it is dated "02/2020" and signed "Author: X. Meng" — a coauthor of the cited paper, and the same month the paper was first posted to arXiv — though the notebook itself carries no execution record tying it to a specific figure or table in the paper.

  The `Remez` function solves, at each iteration, a $(d{+}1)\times(d{+}1)$ linear system in arbitrary precision (Julia `BigFloat`) for the Chebyshev-basis coefficients of the odd, even or unconstrained-parity approximant plus a sign-alternating deviation term, evaluated at the current reference points; finds the sign-change roots of the residual with Brent's method; searches each bracketed sub-interval by dense sampling for the point of largest residual magnitude to move the reference to; and repeats until the deviation and the achieved $L^\infty$ error agree to within its stopping tolerance or `maxiter` is reached. The notebook states its own caveat directly: "We note that our implementation does not strictly follow the reference. The algorithm may crush [sic] when the degree is very large or the problem is ill-conditioned," and cites its reference as E. W. Cheney's *Introduction to Approximation Theory* — the same classical-approximation-theory framing this record's own `conditions` field gives Remez.

  Inputs are a Julia-callable target function, a requested `degree`, a `parity` flag (0 even, 1 odd, 2 unconstrained), the interval endpoints, and either an initial reference (the roots of a Chebyshev polynomial by default) or a user-supplied one. The notebook's own worked example approximates $1/(4\kappa x)$ — a rescaled $1/x$ — over $[1/\kappa,1]$ at $\kappa=10$, even parity, degree 60, target accuracy $10^{-6}$, using 512-bit precision.

  Repository `qsppack/QSPPACK`, path `Solvers/JuliaSolver/Remez.ipynb`, function `Remez(targetf, parity, degree, xapp, lef, rig, maxiter, sample_size, eps)`. Its output — the Chebyshev-basis coefficients — is, per the example cell's own comment, meant to be written to a `.mat` file so the caller can "solve for corresponding phase factors via optimization method": the hand-off to the phase-factor stage this record's `summary` describes. That optimization-based phase-factor solver lives separately, in `Solvers/Optimization`, and is Dong, Meng, Whaley and Lin's own method from the same paper.

• Chebfun's `minimax` (barycentric Remez exchange)

  Chebfun's general-purpose best-approximation routine, in `minimax.m` at the root of the `chebfun/chebfun` repository. Its own header states plainly: "P = MINIMAX(F, M) computes the minimax polynomial approximation of degree M to the real function F using the Remez algorithm," where F may be a chebfun, a function handle, or a string, over a domain that defaults to $[-1,1]$ — a general target function and interval, not a filter's frequency response. The header adds that the file "supersedes REMEZ," an earlier, separate implementation in the same toolbox.

  In the two-argument form `MINIMAX(F, M)` this record's `about` field quotes — Chebfun's ordinary polynomial-approximation call — the reference set is seeded once at $N+2$ Chebyshev extrema (`xk = chebpts(N + 2, ...)`) and handed to a single call of the internal `minimaxKernel`, which on each round calls a local subroutine named `exchange`, whose own header states the algorithm outright: "EXCHANGE Modify an equioscillation reference using the Remez algorithm," returning "a FLAG indicating whether there were at least N+2 alternating extrema of the error to form the next reference." This polynomial path makes no symmetry assumption (`symFlag` is hardcoded to 0) and, unlike QSPPACK's solver above, takes no parity argument either; automatic detection of odd/even symmetry via `adjustDegreesForSymmetries`, and the fallback chain that retries a failed reference construction — first with a Caratheodory–Fejer-based initial reference, then an AAA-Lawson-based one, then two cumulative-distribution-function-based ones, before raising `CHEBFUN:CHEBFUN:minimax:failure` — both belong only to `minimax.m`'s separate rational-approximation branch, reached when the caller supplies an explicit third argument.

  A target `f` (chebfun, function handle, or string), a requested degree `m`, an optional domain `[a, b]` (default $[-1,1]$), and an optional user-supplied starting reference `xk`; a tolerance on the equioscillation error and a maximum iteration count are also accepted as name-value options.

  Repository `chebfun/chebfun`, file `minimax.m`; the reference-exchange step is the local function `exchange` in the same file, and the Chebyshev-point initialization for the polynomial case is the call to `chebpts` inside `minimax.m`'s own top-level branch for $n=0$.

• SciPy's `scipy.signal.remez` (Parks–McClellan FIR design)

  SciPy's public entry point for the Remez exchange algorithm, but scoped to one long-standing application rather than to a general target function: designing finite-impulse-response filter taps against a piecewise-constant frequency-response specification. Its docstring states directly: "Calculate the minimax optimal filter using the Remez exchange algorithm," computing "filter-coefficients for the finite impulse response (FIR) filter whose transfer function minimizes the maximum error between the desired gain and the realized gain in the specified frequency bands." Read honestly against this record: the caller supplies band edges and one desired gain per band, not an arbitrary function handle and an arbitrary interval, so this is the filter-design specialization of Remez rather than the general-purpose construction the record describes.

  The Python wrapper in `scipy/signal/_fir_filter_design.py` validates its arguments and forwards to a compiled extension function, `_sigtools._remez` (defined in `scipy/signal/_sigtoolsmodule.cc`, registered under the name `_remez`). That file credits its algorithm directly, in a header comment: "Code taken from remez.c by Erik Kvaleberg which was converted from an original FORTRAN by... JAMES H. MCCLELLAN... THOMAS W. PARKS... LAWRENCE R. RABINER." Its own `remez()` subroutine comment states the algorithm in the same terms this record's own theory hop uses: "THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS FUNCTION WITH A SUM OF COSINES," working from "A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE EXTREMAL FREQUENCIES," and it "MINIMIZES THE CHEBYSHEV ERROR BY DETERMINING THE BSMINEST [sic] LOCATION OF THE EXTREMAL FREQUENCIES... AND THEN CALCULATES THE COEFFICIENTS OF THE BEST APPROXIMATION." A sum of cosines in the frequency variable is, under $x=\cos\omega$, a linear combination of Chebyshev polynomials, so the routine is the same exchange algorithm on an algebraic polynomial in disguise.

  `numtaps` (filter length), `bands` (a monotonic sequence of frequency band edges), `desired` (one gain value per band), an optional per-band `weight`, a filter `type` ('bandpass', 'differentiator' or 'hilbert'), `maxiter` (default 25) and `grid_density` (default 16, so the dense grid used internally has `(numtaps + 1) * grid_density` points).

  Public function `scipy.signal.remez`, defined in `scipy/signal/_fir_filter_design.py`; it calls the C extension function registered as `_sigtools._remez` in `scipy/signal/_sigtoolsmodule.cc`, which in turn calls the file's own internal `pre_remez()` and `remez()` routines.

What it needs

Nothing below this — it bottoms out here.

Other ways to fill the same slot

Different approaches

• Truncated Chebyshev expansion

  Expand the target in Chebyshev polynomials and truncate once the coefficients have fallen below the error budget. For $1/x$ the expansion is taken of the odd function $f(x) = (1 − (1 − x²)^b)/x$; for $e^{-ixt}$ the Jacobi-Anger identity supplies Bessel coefficients that decay super-exponentially once the order passes about $t$.

In the Atlas

No record in the Atlas covers this yet. The catalogue is circuits and primitives; this part of the literature is not in it.

Sources

• Efficient phase-factor evaluation in quantum signal processing↗Yulong Dong, Xiang Meng, K. Birgitta Whaley, Lin Lin · 2020theory yes · simulation yes · hardware no