Termination conditions

[ivan-pi]

Is there any interest to have custom termination conditions? This could be implemented with a user-provided callback (procedure pointer). If nothing else it could be useful for comparison with other root-finders in terms of speed or number of function evaluations needed for a given accuracy.

Below is a summary of what the other libraries have (in addition to the classic maximum iterations criterion).

---

Boost

The default condition used by Boost is

when the relative distance between a and b is less than four times the machine epsilon for T, or 2^{1-bits}, whichever is the larger. In other words, you set bits to the number of bits of precision you want in the result. The minimal tolerance of four times the machine epsilon of type T is required to ensure that we get back a bracketing interval, since this must clearly be at greater than one epsilon in size. While in theory a maximum distance of twice machine epsilon is possible to achieve, in practice this results in a great deal of "thrashing" given that the function whose root is being found can only ever be accurate to 1 epsilon at best.

To get "maximum" precision you would do something like this:

  int digits = std::numeric_limits<T>::digits;  // Maximum possible binary digits accuracy for type T.
  // Some fraction of digits is used to control how accurate to try to make the result.
  int get_digits = digits - 3;                  // We have to have a non-zero interval at each step, so
                                                // maximum accuracy is digits - 1.  But we also have to
                                                // allow for inaccuracy in f(x), otherwise the last few
                                                // iterations just thrash around.

In Fortran the equivalent would be

real(dp) :: x, f
integer :: get_digits
get_digits = digits(x) - 3

Boost also gives the possibility to specify a custom termination functor for the root bracket:

template <class T>
struct eps_custom
{
   eps_custom();
   bool operator()(const T& a, const T& b) const; // user-defined
};

It provides also a few predefined termination condition functors to stop at the nearest integer of the true root, or at the integer ceiling/floor of the true root. I don't really know what is the usage case for these.

SciPy

The methods in Scipy generally use a combination of absolute and relative tolerance with small variations picking the factor multiplying the relative tolerance:

• Bisection

    dm = xb - xa;
 // ...
        dm *= .5;
        if (fm == 0 || fabs(dm) < xtol + rtol*fabs(xm)) {
            solver_stats->error_num = CONVERGED;
            return xm;
        }

• Brent's method (same as bisection other just multiplied by 1/2 on both sides)

        delta = (xtol + rtol*fabs(xcur))/2;
        sbis = (xblk - xcur)/2;
        if (fcur == 0 || fabs(sbis) < delta) {
            solver_stats->error_num = CONVERGED;
            return xcur;
        }

• Ridder's method (uses full interval instead of half the interval, tolerance calculated with new estimate)

        tol = xtol + rtol*xn;
        if (fn == 0.0 || fabs(xb - xa) < tol) {
            solver_stats->error_num = CONVERGED;
            return xn;
        }

• TOMS748 (uses np.isclose(a,b) which is defined as absolute(a - b) <= (atol + rtol * absolute(b)))

        a, b = self.ab[:2]
        if np.isclose(a, b, rtol=self.rtol, atol=self.xtol):
            return _ECONVERGED, sum(self.ab) / 2.0

Roots.jl

For most algorithms, convergence is decided when

The value |f(x_n)| <= tol with tol = max(atol, abs(x_n)*rtol), or

the values x_n ≈ x_{n-1} with tolerances xatol and xrtolandf(x_n) ≈ 0with a _relaxed_ tolerance based onatolandrtol`.

GSL

GSL provides three search stopping criteria:

• gsl_root_test_interval - $| a - b| < epsabs + epsrel \min(|a|,|b|)$, with some special handling when $[a,b]$ contains the origin
• gsl_root_test_delta - $|x_1 - x_0| < \hbox{\it epsabs} + \hbox{\it epsrel}\,|x_1|$
• gsl_root_test_residual - $|f| < epsabs$

Since the GSL solver uses a reverse-communication interface, a user can easily supply his own stopping condition too.