Fixed the Remez exchange silently converging to a
reference-local fixed point that is not the global minimax,
returning a materially suboptimal approximation with no warning (for
example, atan on [0, 3] at degree 3 reported an expected
error 1.8 times smaller than the returned approximation's true maximum
error). Two structural defects were responsible: the root search only
ever examined the intervals between consecutive current reference
points, so a sign change between an endpoint and the outermost reference
point was invisible; and the maximize/minimize schedule was seeded at
the lower endpoint and flipped mechanically, so an error curve with one
more oscillation than the schedule assumed desynced the exchange into
returning roots instead of extrema. Both are fixed by the exchange
redesign described under “Changed”.
Fixed a crash (“segfault from C stack overflow”) evaluating a Chebyshev-basis approximation at several million or more points at once.
Added a clear error, instead of undefined behavior, when an internal Chebyshev evaluation would require more than 2^31 matrix cells.
Fixed lower > upper (or lower == upper) being silently
accepted and producing a badly suboptimal approximation with no warning;
this now raises a clear error suggesting the arguments be swapped.
Fixed the Remez-iteration stagnation check (isUnchanging)
reporting convergence as “unchanging” when errors were still
improving rapidly between iterations. The ratio test was checked in only
one direction. Now, it requires the ratio to be close to 1 in either
direction before treating the solution as stalled.
Fixed opts$ztol using the monomial-basis coefficient scale
even when the Chebyshev basis was requested, which could incorrectly zero
out genuine small-but-relevant Chebyshev coefficients on ranges
intersecting [-1, 1].
Fixed the degree-n+1 restart's “uppermost coefficient
is effectively 0” check silently accepting any negative top coefficient
regardless of magnitude (a missing abs()), which could return a
polynomial one degree lower than warranted.
Fixed severe numerical conditioning loss when fitting the
Chebyshev basis on ranges far from [-1, 1] (e.g. degree 10 on
[5, 6] had a condition number around 10^{24}), which could
cap the achievable accuracy or force convergence warnings well short of
the true minimax error (that [5, 6] case: observed error improved
from \sim 1.8 \times 10^{-10} to \sim 3 \times 10^{-12},
matching the true value). See the Breaking change below.
The Remez exchange was redesigned to be driven by the error curve itself rather than by the current reference's structure. Roots of the error equation are now located over the whole approximation interval from an oversampled grid. Each sign-region's search direction is taken from the error's own sign there instead of a mechanically alternated schedule, and the interval endpoints are always retained as exchange candidates. When the curve offers more alternating extrema than the reference holds, the new reference is the alternating window that contains the global extremum and maximizes its smallest error—the classical conditions for the leveled error to increase strictly toward the true minimax. Converged results for healthy inputs agree with version 0.5.0 at tolerance level but are not bitwise identical. Inputs that previously exited through a singular-matrix error or a stall may now exit through a different path as the iteration trajectory differs.
Breaking: Chebyshev-basis approximation now internally
maps the requested range [l, u] to [-1, 1] before evaluating
T_k, fixing the conditioning issue noted above. a (and
b) are therefore coefficients of T_k of the internally
mapped variable, not T_k(x), whenever the fitted range is not
[-1, 1] Results on [-1, 1] are unchanged. aMono and
bMono remain coefficients in the natural, raw x in every
case, and are the stable way to consume a Chebyshev-basis result without
needing to know about the internal map. See minimaxApprox
for a documented accuracy caveat on this conversion at high degree on
wide, far-off-[-1,1] ranges.
Added upfront validation of fn, lower, upper,
degree, and opts in minimaxApprox(), replacing
cryptic low-level failures (e.g. “missing value where TRUE/FALSE
needed”, “invalid argument type”) with clear, argument-specific
error messages. fn must be a function whose first argument is
x, although primitives such as exp and sin are still
accepted; lower/upper must be finite non-missing numeric
scalars; degree must be finite and non-missing;
opts$maxiter/miniter/conviter must be single
positive integers and opts$tol/convrat/tailtol/
ztol must be finite numeric scalars when supplied. The
tailtol and ztol options may still be NULL by design.
Corrected wording in the degree-n+1 restart messages, which
described the returned polynomial's length as its degree
(a degree-n polynomial has n+1 coefficients).
isOscil, used internally to test for the expected
alternating-sign error pattern, no longer propagates NA when passed
a NaN/NA error value. Previously this could have errored
deep inside the convergence check instead of failing cleanly.
The internal zero-basis perturbation, used in relative-error mode
when the approximated function is exactly 0 at a candidate extremum, now
escalates to a magnitude-scaled step when the fixed 1e-12 nudge
would be a no-op (points with |x| \gtrsim 4.5 \times 10^3). It
remains identical to the previous behavior at ordinary magnitudes.
Functions that the requested degree resolves to, or near, machine
floating-point precision—including functions that are exactly
polynomials, such as minimaxApprox(function(x) x^2, 0, 1, 2), and
well-behaved cases such as minimaxApprox(exp, -1, 1, 14)—no
longer hard-error or emit misleading non-convergence warnings. This case
can surface two ways, both now handled: (1) both the degree-n and
degree-n+1 Remez solves fail singular, previously
“neither converged”; or (2) on some BLAS/LAPACK
platforms the same case does not report singular but instead runs to
maxiter (or stalls “unchanging”) while genuinely wandering
at the machine-precision floor. In either case, when there is no minimax
problem left to solve, the degree-n interpolant through the
Chebyshev reference is returned instead, together with a warning that it
is not technically a Remez result but that its error is below the level
at which a distinct minimax solution can be discerned. Genuinely
non-representable cases, where no polynomial comes near the function,
still raise the original error or convergence warning.
Updated the package Title and Description to
reflect the expanded scope: minimax approximation by polynomials and
rational functions via both the classical and barycentric algorithms.
Added a reference-local certificate check to all
approximation paths. At a converged exit, the returned approximation's
maximum error on a dense grid is compared against its leveled error, and a
warning is raised when the former materially exceeds the latter. In that
case the reported expected error is a valid lower bound on the true
minimax error (de la Vallée Poussin) rather than a certificate or
gurantee of it. This can occur for degrees that are non-normal for the
function (e.g. even or odd functions at generic degrees) or when the
exchange stalls at a reference-local fixed point. The check is skipped
where its ratio would be noise: results at the machine-precision floor and
relative-error results with an exact zero of fn on the grid.
Added a "Barycentric" option to the basis argument,
implementing the barycentric-Remez algorithm for the best polynomial
approximation of Pachón and Trefethen (2009). It represents the
trial polynomial by its values at the reference points rather than by
coefficients in a fixed basis, avoiding the ill-conditioned linear solve
that limits the achievable degree in the classical bases. In consequence,
it converges at degrees and on ranges where the Monomial and Chebyshev
bases fail with a singular-matrix error—for example,
minimaxApprox(exp, -1, 1, 14) (and degree 50), and exp on
[5, 6] at degree 10—and, using capacity scaling, remains stable on
very wide or far-off-center intervals. Exactly-representable and
machine-precision-resolved functions are handled natively (leveled error
0, no rescue needed). Relative error is supported and a function
that is exactly zero at a reference point (typically an interval
endpoint), where the relative minimax does not exist, now raises a clear
error. The returned object carries the barycentric representation and a
coefficient-conversion residual diagnostic (convResid). See
minimaxApprox.
Added rational approximation support to the
"Barycentric" basis, implementing a subset of the adaptive
barycentric minimax algorithm of Filip, Nakatsukasa, Trefethen, and
Beckermann (2018), with the initial reference taken from the error
extrema of an adaptive Antoulas–Anderson (AAA) approximant
(Nakatsukasa, Sète, and Trefethen, 2018). The leveled error is
computed as an eigenvalue of a small symmetric matrix. Therefore, as for
the polynomial case, there is no ill-conditioned linear solve. Rather, the
classical smooth-function results are reproduced to convergence tolerance
while cases such as |x| converge to the certified minimax where
reference-based classical iterations can struggle. This release only
implements the algorithm using absolute error. Passing
relErr = TRUE will return an error. Degenerate or defective
requests (e.g. even degrees for an even function) are detected and
reported with a clear error rather than repaired by degree reduction. A
trial denominator with a zero in the interval is an error mirroring the
classical rational path. The returned object carries a and b
normalized as in the classical path, aMono and bMono, the
barycentric representation, and a named two-component convResid.
See minimaxApprox.
Added the package's first vignette, “Using minimaxApprox”,
covering the minimax criterion, the three bases and how to choose among
them, rational approximation and its barycentric-basis limitations,
relative error, the Warning/convResid diagnostics, and the
handful of opts most users need.
The xi argument to minimaxApprox is deprecated
and will be removed in the next release. The redesigned exchange
algorithm derives its reference from the error curve at every
iteration, making a user-supplied initial reference unnecessary.
Passing xi now emits a deprecation warning.
Fixed error in one package metadata unit test.
Breaking: Converted C code for Compensated Horner
method to use long double internally to eke out a bit more
precision. One unit test had to be updated as what was formerly an error
now converged.
Updated Github-related files and workflows.
Updated ‘WORDLIST’.
Reduced overwriting of base R function names.
Simplified code and nesting of various functions.
Package metadata tests now use vapply for type safety.
Efficiency tweaks to C code.
Use native CRAN DOI and update ‘README’.
Expanded package metadata testing.
If an explicitly passed conviter is larger than
miniter, then miniter is set to conviter.
Documentation made more consistent.
C code tweaked for clarity and efficiency.
Updated unit tests.
Set up Github workflow for testing with rhub.
One unit test adjusted to pass on Professor Ripley's OpenBLAS test machine which has different tolerances than the standard twelve on CRAN, six on Github, three on R Win Builder, and the one on R Mac Builder.
Breaking: Analysis now uses the Chebyshev polynomials of
the first kind, T_k, as the default basis instead of the monomials,
x^k. Monomials can still be used by passing the appropriate
basis parameter.
The internal functions underlying the Chebyshev polynomial calculations were ported to C for speed.
Breaking: Changed names in return object. What was
Basis is now Extrema. The basis parameter is now used
to capture whether monomials or Chebyshev polynomials should be/were used.
Breaking: User-facing functions now have basis
parameter to determine which basis to use; see appropriate help files.
Parameter positions have thus changed.
The print, plot, and coef methods were updated
to recognize or allow basis selection.
The print method now shows 14 digits for the convergence
ratio by default.
Internal code was streamlined and tolerances were tightened.
Updated documentation for new functionality.
The C files were split and header files created to keep Chebyshev and Horner code separate from each other.
Fixed potential division-by-zero problem in the “unchanging” test.
Fixed outdated parameter name in documentation.
Fixed incorrect ‘NEWS’ entries.
Breaking: More descriptive names in return object. See the
Value section of help("minimaxApprox") for details.
Breaking: When using relative error for convergence and one
of the basis points has a functional value of 0, instead of stopping with
an error, the algorithm will perturb that basis point by 1e-12 in
the appropriate direction and proceed with a warning.
Breaking: Forced convergence to look only at the
ratio if the difference between tested values is \le
.Machine$double.eps.
Breaking: Increased default for conviter from 10 to
30 as tests on some values now need more iterations to achieve ratio
convergence given a distance \le .Machine$double.eps.
Added unit tests for additional functionality.
Fixed a misnamed option value in some unit tests.
Refactored ‘README’ to split out security and contributions wordings to their own respective files and link the citation to the file on CRAN.
Changed some tests to be more acceptable for CRAN.
Added LTO_OPT to ‘Makevars’ and UseLTO to
‘DESCRIPTION’.
Test for calculating polynomial of degree n + 1 when one of
degree n fails, while passing both R Mac Builder and
Github's Mac platform, failed on CRAN's testbed.
So test is no longer run for Mac on CRAN.
When the polynomial algorithm fails with a singular matrix error
for degree n, it will now look for a polynomial of degree
n + 1. If found, and the contribution of that coefficient
is less than the new option tailtol (default 1e-10), it will
return the found polynomial without the uppermost coefficient—which is
effectively 0. It will message the user to that effect. Passing
NULL skips the test. (see
Issue 2 with
thanks to Hans W. Borchers).
New option ztol which acts similarly to
zapsmall in that it will zero out coefficient values—both
polynomial and rational—whose contributions are \le ztol
at each step of the iteration. The default, NULL, skips the test.
Ported Compensated Horner Scheme algorithms to C for speed. The algorithms are now roughly 10x faster.
Changed linear solution algorithm to use qr.solve when the
default solve fails with a singular matrix issue. This may
allow for higher degree polynomials in certain cases.
Breaking: Changed name of ‘round’ parameter
in print method to ‘digits’ to prevent any confusion with
base function round.
Updated unit tests.
Cleaned up and updated some documentation and various messages.
Tightened tolerances and cleaned up some internal functions.
Added ‘LICENSE’ to ‘.Rbuildignore’ and removed from ‘DESCRIPTION’ as suggested on R-pkg-devel.
Internal polynomial evaluations now uses Compensated Horner Scheme of Langlois et al. (2006) for increased precision and stability.
Clarified some error messages.
Remove one test which fails on M1 Mac (thanks to Professor Brian Ripley) and PowerPC (Issue 1; thanks to Sergey Fedorov) due to it returning a singular matrix where on Linux and Windows it does not.
Added file LICENSE with text of MPL-2.0.
Added copyright and license notices to source code where it was missing.
Added unit tests.
Correct documentation and make more consistent.
Provide better examples for minimaxApprox.
Updated README with live CRAN links.
Initial release.
So long as the package remains in its experimental development state—noted by a 0 major version—the API may change without incrementing the major version. Please read these development notes carefully. Breaking changes will be prefaced by “Breaking:”.