[std-interval] More on interval computations as proofs

[Guillaume Melquiond]

Le mercredi 04 octobre 2006 à 16:29 -0700, Lawrence Crowl a écrit :

> My understanding is that interval computations tend to require
> some rewriting anyway.  That is, straightforward translations of
> plain floating-point computations tend to have wide intervals.
> Am I mistaken?

Let's take a concrete example: Newton's iteration for finding the roots
of a function f (its derivative will be d).

template<typename T> T f(T) { /* a complicated function */ }
template<typename T> T d(T) { /* an even more complicated derivative */ }

Starting from a point x, the following iteration might converge
quadratically towards a root:

x = x - f(x) / d(x);

Now if you translate the algorithm to intervals by replacing the
floating-point variable x by an interval X:

X = X - f(X) / d(X);

you will get an iteration that does not converge: interval X just grows
larger until containing whole. This is a typical example where
straightforward translation is a bad use of interval arithmetic.

Let's consider interval Newton's iteration instead:

X = intersect(X, midpoint(X) - div_rel(f(interval(midpoint(X))), d(X)));

The important point is that f and d haven't been changed; they really
are straightforward translations of floating-point computations to
interval computations. Moreover, if the computed interval before the
intersection is strictly contained in X and if the discontinuity flag is
not raised: you know that X will not vanish later on, as there exists a
root in it.

So, to qualify your point on straightforward translation, I would say
that you are right for resolution algorithms (finding roots, solving
differential equations, and so on), as they simply won't work. But for
evaluation functions, you don't have to rewrite them. Just change
"double" to "interval<double>", and the interval resolution algorithms
will take care of these evaluation functions.

Best regards,

Guillaume

PS: obviously, it's always possible to build a floating-point function
that will behave wrong once translated to intervals, for example by
inserting (1. / (1e-300 + (x-x)) - 1. / 1e-300) in it.