Re: [Caml-list] Rounding mode

[Brian Hurt <bhurt@spnz.org>]

On Mon, 10 Nov 2003, Eric Dahlman wrote:

> Brian Hurt wrote:
> 
> >
> >Simplistic discussion of why changing the rounding modes is usefull for 
> >numeric programming.  If you know numerical analysis, stop reading now 
> >:-).
> >  
> >
> I read it anyway ;-)

The reason I asked you to stop reading was that I was going to be 
massively simplifying a subject that is in effect it's own branch of PhD 
level mathematics.  The result will be somewhat like those books on 
relativity aimed at grade school kids.  Physicists who read them are 
regularly horrified that important subtlities are glossed over, hand waved 
past, or simply blatantly ignored.  The reason is, of course, I was 
writting an email not a book :-).

Also, I'm far from an expert in this field myself (although I know a 
couple).  I know just enough to know how much I don't know.

[Snip happens]

> Now I don't know that it is this simple.  

You're right, it's not at all that simple.  For any given two rounding 
modes, it's possible to construct algorithms for which both rounding modes 
will error in the same direction- i.e. the two results will not bracket 
the "real" answer, and will also misdirect you as to how large your error 
bounds really are.  Also, doing this trick only tells you about the data 
you are using.  Algorithms can have unexpected failure modes where they 
lose a lot more presision.  A classic example of this is subtraction.  You 
may know that x = 1.00000000000002 to 15 digits of precision, and y = 
1.00000000000001 to 15 digits of precision, but you only know x-y = 
0.00000000000001 to 1 digit of precision.  Most of your precision just got 
canceled.  Opps.

Interval arithmetic has similiar problems.  There are no mechanistic 
replacements for numerically analyzing your algorithms.  Both are usefull 
as experimental evidence that your analysis is correct.  But numerically 
analyzing your algorithms is a subtle and tricky- which is why my advice 
is to leave to the numerical analysts, and only use algorithms they have 
blessed (if possible).

Hmm.  Actually, I wonder if type systems could be extended to detect 
numerical instabilities?

> This is incorrect because you are thinking about real numbers while 
> working with floating point numbers, this is especially dangerous when 
> you are operating near or below the epsilon for the particular floating 
> point representation like we are now.  If you were to look at the set of 
> numbers which are precisely represented by a given floating point 
> representation you would see that it is concentrated near 0 and that it 
> tends to thin out as you move away from there.   The distance  between 
> any  two adjacent floating point numbers is highly variable and not 
> uniform or "nice" by almost any measure.  I really am having a hard time 
> figuring out what you can say about the tiny e terms in this case as we 
> move between regions with differing "floating point resolutions".  

Most of the numerical proofs I have seen assume error is a scale factor,
not an additive factor.  I.e., it's x * (1 + e), not x + e.  Then x*y with
error terms becomes x*(1+e_x)*y*(1+e_y) = x*y*(1 + e_x + e_y + e_x*e_y).  
This way I sidestep any assumptions about the "sizes" of x and y.  You
still have to deal with the corner cases where x and y become too small,
and you start hitting denormalized numbers.  This is more accurate, but
more complex math.  Thus the simplification to (x + e), while making sure
the conclusions are still more or less correct.

> As for a numerically stable algorithm being able to produce a result 
> which is more accurate than the original inputs that is total hogwash.  
> You cannot ever have more accuracy than was present in your original 
> input, ever. Nihil ex nihilo.

Yep.  But you can have less, even none what so ever.  Stable algorithms
are stable precisely because the errors they introduce tend to cancel each
other out instead of accumulating.  Which was the point I was trying to
make- I'll admit I was speaking inaccurately when I implied that stable
algorithms somehow magically added precision that wasn't there already.

You can make both type I and type II errors (warping statistic
terminology), and both are incorrect.  Type I errors would here be
assuming that because the computer printed out 15 digits of output, all 15
digits are to be beleived (despite the fact that the inputs are only
accurate to 3 digits).  Type II errors are assuming that errors only
accumulate.  This leads to demands for, for example, 2048-bit floating
point types.  The hope here is that if you can make the errors small
enough, you may still have some precision when you're done.

In my experience, type II errors are most often made by programmers who 
have recently "discovered" that floating point introduces errors into 
calculations.  Which is why I was addressing it early- the solution isn't 
to throw more precision at the problem, the solution is to make the errors 
introduced cancel each other.  And in so doing, unwittingly played to 
making type I errors.

[ snip: randomized rounding as the default ]

> As an oldby I would be disconcerted by this as well.  All you are doing
> is trading one sort of numerical noise for another which you guess
> should be better behaved but in the absence of a neat proof to the
> contrary I don't believe is.  There are just *SO* many sources of huger
> error terms in numerical calculations that the tiny rounding errors are
> really lost in the noise.  I even tried to come up with a case where it
> would matter but that got so farcical that I gave up on it.

All stable algorithms I know of remain stable in the face of randomized 
rounding.  There are some algorithms which are instable with any fixed 
rounding scheme, but which are stable with randomized rounding.

Consider the following case: adding a list of three numbers together-
1.0, 0.5e-15, and 0.5e-15.  Adding them in that order is unstable no 
matter what fixed rounding you use.  The correct answer is, obviously, 
1.000000000000001.  If you always round up, you get 1.000000000000002, 
while if you always round down, you get 1.000000000000000.  With 
randomized rounding, you have a 50% shot of getting the right answer- more 
likely than any other answer.  Now, replace that three-element list with 
1.0 followed by 0.5e-15 repeated a million times.

The stable algorithm here is to sort your long list of integers to add 
first, and always add from smallest to largest.  If you do this with the 
three element list, you'll add 0.5e-15 to itself, getting 1e-15, a number 
large enough not get totally rounded off when added to 1.0.  Adding 
0.5e-15 to itself a million times gets you 0.5e-9, which when added to 1.0 
always gets the right answer, with very little dependence on what rounding 
happens.  Replacing an instable algorithm with a stable algorithm is 
always a good idea.  Of course, in this case, you're replacing an O(n) 
algorithm with an O(n log n) algorithm.

There are reasons why people who know what they're doing may want to play 
with rounding modes.  And while randomized rounding may in some sense be 
"better", I agree that making it the default it a bad idea.  I'm not even 
sure if all architectures offer it.

I've got to admit to a certain fondness for the writtings of Prof. Kahan:
http://www.cs.berkeley.edu/~wkahan/

In fact, I'll go so far as saying the following should be required 
reading for the Ocaml developers:
http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf

(Well, except for operator overloading, I'll disagree with Kahan on that. 
 Although for what Kahan wants, I wonder if Ocaml's ability to define new
infix operators isn't enough, or even better.)

> 
> What I was really wondering about when I first asked this question was 
> if the other person had exactly that kind of example where they were 
> using some clever trick to usefully misapply floating point calculations 
> in their problem domain.

I agree- this is an incredibly bad idea.

-- 
"Usenet is like a herd of performing elephants with diarrhea -- massive,
difficult to redirect, awe-inspiring, entertaining, and a source of
mind-boggling amounts of excrement when you least expect it."
                                - Gene Spafford 
Brian


-------------------
To unsubscribe, mail caml-list-request@inria.fr Archives: http://caml.inria.fr
Bug reports: http://caml.inria.fr/bin/caml-bugs FAQ: http://caml.inria.fr/FAQ/
Beginner's list: http://groups.yahoo.com/group/ocaml_beginners