Re: Floating point inaccuracies

[Dirk Laurie <dpl@...>]

On Wed, Dec 01, 2010 at 03:03:12AM +0200, Andrew Lentvorski wrote:
> On 11/30/10 8:48 AM, steve donovan wrote:
> > On Tue, Nov 30, 2010 at 6:33 PM, Roberto Ierusalimschy
> > <roberto@inf.puc-rio.br>  wrote:
> >> Something like linterval?
> >> (http://www.tecgraf.puc-rio.br/~lhf/ftp/lua/#linterval)
> >
> > Yes, that's exactly it. The source of all wisdom has some very
> > readable paragraphs:
> >
> > http://en.wikipedia.org/wiki/Interval_arithmetic
> 
> What it fails to mention is that intervals tend to blow up faster than 
> actual errors for most complex computations.
> 
> So, for example, if you have to invert a matrix, the interval result 
> winds up with gigantic values and huge intervals while the standard 
> binary floating point result winds up within delta of the actual value 
> most of the time.
> 
We're getting seriously off-topic, but these statements, like
all statements containing phrases like "tend to" and "most of 
the time", are only true because those terms are user-definable.

The wiki article fails to mention over-pessimistic blowup of
intervals, because that is a simplistic and dated view.

Facts:
 
1. If N arithmetic operations each committing an absolute error 
   bounded by eps are made, the absolute error will be bounded
   by N*eps, and interval arithmetic does just that.
2. If these errors are uncorrelated, the expected value of the
   actual error is however a small constant multiple of 
   sqrt(N)*eps.  You have actually lost say three digits of 
   precision, but interval arithmetic says correctly that 
   you could be unlucky and lose six digits.
   If "N*eps" is "huge" and "sqrt(N)*eps" is "delta", how large
   must N be?  I can't figure that out.
3. There is a certain class of computational algorithm, called
   "backward stable", in which the error caused by the algorithm
   is not much larger than the error caused by forcing the data
   into your floating-point system, i.e. the algorithm gives
   the exact solution to a problem close to the one you actually
   set out to solve.  Standard methods for matrix inversion have
   this property, which you can recognize by the fact that 
   A*roughinverse(A) is close to I, even though roughinverse(A)
   can be quite far from inverse(A).
4. Intervals can't represent the set of matrices B with the property
   "A*B is close to I".  For an ill-conditioned matrix, a matrix
   interval that contains both inverse(A) and roughinverse(A)
   also contains matrices B such that A*B is not close to I. In this 
   sense, yes, interval methods give pessimistic results, but only 
   when something is rotten in the state of the original problem.
5. There are specialized interval matrix inversion routines
   that produce realistic-sized intervals, but they are much
   more sophisticated that just doing ordinary Gaussian elimination
   on interval-valued objects.  For example, the interval Newton
   method described in the Wiki article can be used to refine the
   intervals.

Dirk