Using minimaxApprox
What is minimax approximation?
Least-squares fitting minimizes the average squared error;
interpolation forces the fit through chosen points and lets the error do
whatever it likes in between. Minimax approximation asks a different
question: what is the polynomial (or rational function) of a given
degree whose largest error anywhere on the interval is
as small as possible?
The answer has a distinctive fingerprint: at the optimum, the error
curve oscillates between its extreme positive and negative values at
least degree + 2 times (degree1 + degree2 + 2
for rational approximation), a property called equioscillation.
Below, a degree-5 minimax fit to sin(x) on
[0, 1], plotted with plot.minimaxApprox:
fit <- minimaxApprox(sin, 0, 1, 5)
plot(fit)
The error touches its bounds seven times, alternating sign every
time—no other degree-5 polynomial can do better in the worst case.
The exchange algorithm used to find such a fit dates to Remez (1934;
the translation cited in ?minimaxApprox is Remez, 1962),
refined for rational approximation by Fraser and Hart (1962) and Cody,
Fraser, and Hart (1968) into the form this package follows most closely.
A newer barycentric formulation—used here as an additional
basis option—represents the trial function by its values rather than by
coefficients, following Pachón and Trefethen (2009) for polynomials and
Filip, Nakatsukasa, Trefethen, and Beckermann (2018) for rational
functions; see Choosing a basis below.
Quick start
fit <- minimaxApprox(exp, 0, 1, 5)
fit
#> $`Polynomial Basis`
#> [1] "Chebyshev"
#>
#> $a
#> [1] 1.753388e+00 8.503917e-01 1.052087e-01 8.722105e-03 5.434367e-04
#> [6] 2.715572e-05
#>
#> $aMono
#> [1] 0.99999887 1.00007946 0.49909610 0.17040197 0.03480057 0.01390373
#>
#> $ExpectedAbsError
#> [1] 1.12957e-06
#>
#> $ObservedAbsError
#> [1] 1.12957e-06
#>
#> $Ratio
#> [1] 1
#>
#> $Difference
#> [1] 3.895547e-16
#>
#> $Warnings
#> [1] FALSE
coef() extracts the coefficient vectors:
coef(fit)
#> $a
#> [1] 1.753388e+00 8.503917e-01 1.052087e-01 8.722105e-03 5.434367e-04
#> [6] 2.715572e-05
#>
#> $aMono
#> [1] 0.99999887 1.00007946 0.49909610 0.17040197 0.03480057 0.01390373
a holds the coefficients in whatever basis was used to
fit (Chebyshev by default); aMono is always the
monomial-basis (raw x^k) equivalent, useful when you need
ordinary polynomial coefficients regardless of how the fit was computed.
minimaxEval() evaluates the fit at new points, using the
same basis as the fit unless told otherwise:
x <- seq(0, 1, length.out = 11)
minimaxEval(x, fit)
#> [1] 0.9999989 1.1051718 1.2214020 1.3498579 1.4918250 1.6487224 1.8221193
#> [8] 2.0137519 2.2255400 2.4596039 2.7182807
In the printed summary, ExpectedAbsError is the leveled
error the Remez iteration converged to; ObservedAbsError is
the largest error actually seen at the final reference points.
Ratio and Difference are the two ways the two
are compared for convergence (see Options below); for a
clean converged result they will be at or near 1 and 0 respectively.
Choosing a basis
Three bases are available via the basis argument,
matched case-insensitively and by abbreviation ("m",
"c", "b").
Monomial (x^k) is the most familiar and
the worst conditioned: the Vandermonde-style linear solve it relies on
degrades quickly with degree or with intervals far from
[-1, 1].
Chebyshev (the default) evaluates T_k
on the range mapped to [-1, 1]; off that
interval, a are coefficients of T_k in the
mapped variable, not raw x (a breaking change as
of 0.6.0; see NEWS.Rd). aMono is always in raw
x regardless of basis, recovered by composing the
mapped-Chebyshev coefficients with the affine map; that composition is
itself a monomial-basis operation, so aMono inherits some
of the conditioning cost of the monomial basis at high degree on wide,
far-off-range intervals.
Barycentric represents the trial function by its
values at the reference points rather than by coefficients, so there is
no linear solve to go singular. It is the most numerically robust choice
and converges in cases the classical bases cannot reach. For example, on
[10, 20] at degree 30:
tryCatch(minimaxApprox(exp, 10, 20, 30, basis = "monomial"),
error = function(e) message(conditionMessage(e)))
#> The algorithm neither converged when looking for a polynomial of degree 30 nor when looking for a polynomial of degree 31.
bfit <- minimaxApprox(exp, 10, 20, 30, basis = "barycentric")
bfit$ObsErr
#> [1] 3.72529e-09
bfit$Warning
#> [1] FALSE
Monomial hard-errors; barycentric converges cleanly. (Chebyshev also
manages this case, but only via the machine-precision fallback described
in Trust and diagnostics below, at a coarser error than
barycentric achieves here.)
For best evaluation accuracy, use minimaxEval() on the
barycentric object directly—it evaluates the stored barycentric
representation, which is more accurate than the derived monomial
coefficients. The one potentially lossy step is the aMono
conversion; convResid reports how much accuracy that
conversion cost:
bfit$convResid
#> [1] 2.384186e-07
Rational approximation
Rational approximation fits a ratio of two polynomials, requested as
degree = c(numerator, denominator). It is most useful when
the function has a nearby singularity (in the real or complex plane)
that a polynomial can only chase by adding more and more degree. Compare
a polynomial and a rational fit to a function with a pole just outside
[-1, 1]:
f <- function(x) 1 / (x + 1.05)
poly8 <- minimaxApprox(f, -1, 1, 8)
poly8$ObsErr # degree 8, still poor
#> [1] 0.7854444
rat21 <- minimaxApprox(f, -1, 1, c(2, 1)) # classical (Cody) rational
#> Warning in minimaxApprox(f, -1, 1, c(2, 1)): Convergence to requested ratio and tolerance not achieved in 30 iterations.
#> 30 successive calculated solutions were too close to each other to warrant further iterations.
#> The ratio is Inf times expected and the difference is 7.99361e-15 from the expected.
rat21$ObsErr
#> [1] 7.993606e-15
The degree-(2,1) rational is many orders of magnitude more accurate
than the degree-8 polynomial. The barycentric basis handles the same
request via the adaptive algorithm of Filip, Nakatsukasa, Trefethen, and
Beckermann (2018):
rat21b <- minimaxApprox(f, -1, 1, c(2, 1), basis = "barycentric")
rat21b$ObsErr
#> [1] 3.996803e-15
rat21b$Warning
#> [1] FALSE
The barycentric rational implementation in this release is a
deliberate subset of the published algorithm, with three limitations. It
supports absolute error only; requesting
relErr = TRUE with basis = "b" is an error.
Degenerate or defective requests—for example, an even numerator
and denominator degree for an even function, which has fewer effective
parameters than the requested degree pair—are detected and reported as
an error naming a workaround, rather than repaired by degree
reduction:
g <- function(x) 1 / (x ^ 2 + 0.01)
tryCatch(minimaxApprox(g, -1, 1, c(3, 3), basis = "barycentric"),
error = function(e) message(conditionMessage(e)))
#> The barycentric rational approximation of degree c(3, 3) appears to be degenerate or defective (every solution of the leveled-error equations has a pole inside the approximation interval). Handling degenerate/defective rational approximations (degree reduction) is not implemented for the barycentric basis; try reducing both degrees, e.g. c(2, 2), or use the Chebyshev or monomial basis.
And a trial denominator with a zero inside the approximation interval
is likewise reported as an error, mirroring the classical rational
path.
Relative error
By default, minimaxApprox() minimizes absolute
error. Setting relErr = TRUE instead minimizes relative
error, (fn(x) - approx(x)) / fn(x), useful when the
function’s magnitude varies enormously over the interval and a uniform
absolute tolerance would waste accuracy where the function is small and
starve it where the function is large:
relfit <- minimaxApprox(exp, 0.5, 2, 6, relErr = TRUE)
relfit$ObsErr # this is now a maximum *relative* error, not absolute
#> [1] 4.075071e-07
Relative error is undefined wherever fn is exactly 0, so
an interval whose endpoints (or interior, more generally) include a zero
of fn cannot be fit with relErr = TRUE:
tryCatch(minimaxApprox(sin, -1, 1, 5, relErr = TRUE),
error = function(e) message(conditionMessage(e)))
#> The algorithm neither converged when looking for a polynomial of degree 5 nor when looking for a polynomial of degree 6.
Trust and diagnostics
Every returned object carries a Warning flag.
Warning = TRUE means one of two things happened: the
iteration did not converge to the requested tolerance within
maxiter (see Options below), or—on every
approximation path, as of 0.6.0—the returned result’s true maximum error
was checked against a dense evaluation grid and found to materially
exceed the leveled error the algorithm reported. In the second case,
ExpErr is still a valid lower bound on the true
minimax error, not a certificate of it.
Because of this, a non-converged result’s ObsErr need
not bound the true maximum error—the final reference points may not
track the approximation’s actual extrema. Verify with
minimaxErr() on a dense grid before relying on such a
result:
chk <- minimaxApprox(exp, 5, 6, 10)
#> Warning in minimaxApprox(exp, 5, 6, 10): Convergence to requested ratio and tolerance not achieved in 100 iterations.
#> The ratio is 1.02208 times expected and the difference is 6.50762e-14 from the expected.
chk$Warning
#> [1] TRUE
chk$ObsErr # what the object reports
#> [1] 3.012701e-12
grid <- seq(5, 6, length.out = 100001)
max(abs(minimaxErr(grid, chk))) # what a dense grid shows
#> [1] 3.069545e-12
Here the two are close but not identical—exactly the kind of gap the
Warning flag exists to surface.
Options
Most calls need no opts at all. Two are worth
knowing:
maxiter (default 100) caps the number of Remez
iterations. If a fit hits this cap and warns of non-convergence,
increasing maxiter occasionally resolves it, though a
persistent warning at high maxiter more often signals a
genuinely hard case (near machine precision, or non-normal for the
function—see Trust and diagnostics) than one that
simply needs more iterations.
convrat and tol control how closely the
observed and expected error must agree before the algorithm declares
convergence—by ratio (convrat, default
1 + 1e-9) or by absolute difference (tol,
default 1e-14), whichever check is more informative for the
current error magnitude. Loosening either lets marginal cases converge
sooner, at the cost of a less tightly verified result.
tailtol and ztol also exist but do not
apply to basis = "b", which has no coefficient solve to
restart or prune.
References
Remez, E. I. (1962) General computational methods of Chebyshev
approximation: The problems with linear real parameters. US Atomic
Energy Commission, Division of Technical Information. AEC-tr-4491
Fraser, W. and Hart, J. F. (1962) “On the computation of rational
approximations to continuous functions”, Communications of the
ACM, 5(7), 401-403.
Cody, W. J., Fraser, W. and Hart, J. F. (1968) “Rational Chebyshev
approximation using linear equations”, Numerische Mathematik,
12, 242-251.
Pachón, R. and Trefethen, L. N. (2009) “Barycentric-Remez algorithms
for best polynomial approximation in the chebfun system”, BIT
Numerical Mathematics, 49(4), 721-741.
Berrut, J.-P. and Trefethen, L. N. (2004) “Barycentric Lagrange
interpolation”, SIAM Review, 46(3), 501-517.
Higham, N. J. (2004) “The numerical stability of barycentric Lagrange
interpolation”, IMA Journal of Numerical Analysis, 24(4),
547-556.
Filip, S.-I., Nakatsukasa, Y., Trefethen, L. N. and Beckermann, B.
(2018) “Rational minimax approximation via adaptive barycentric
representations”, SIAM Journal on Scientific Computing, 40(4),
A2427-A2455.
Nakatsukasa, Y., Sète, O. and Trefethen, L. N. (2018) “The AAA
algorithm for rational approximation”, SIAM Journal on Scientific
Computing, 40(3), A1494-A1522.