Re: [Caml-list] Rounding mode

[Eric Dahlman <edahlman@atcorp.com>]

Brian Hurt wrote:

>On Mon, 10 Nov 2003, Eric Dahlman wrote:
>
>  
>
>>Christophe Raffalli wrote:
>>
>>    
>>
>>>I saw a previous discussion about rounding mode for OCaml (in 2000).
>>>
>>>As anyone implemented the functions to change the rounding mode for 
>>>floating point from OCaml (possibly as a patch to add primitives and 
>>>save a C function call) ?
>>>
>>>This is necessary to implement interval arithmetic ...
>>>      
>>>
>>Somewhat off topic but why is this necessary from a numerical math type 
>>of perspective.  I am honestly curious as I don't see how this would 
>>interact with the calculation in a meaningful way.
>>
>>    
>>
>
>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 ;-)

[Snip a description of how floating point works.]

>Consider, for example, multiplying two numbers- x and y.  But you don't 
>have x and y exactly, you have (x + e) and (y + e) for some error amount 
>e.  You multiply them, and you get (x + e)*(y + e) = x*y + x*e + y*e + 
>e*e.  Now, e*e is vanishingly small, we can ignore it.  And, assuming x 
>and y are both about 1.0, our answer is x*y + ~2*e- in other words, we've 
>about doubled the amount of error.  Even starting at 1e-15, doubling it 
>every operation and you eat up all your precision surprisingly fast.
>  
>
>How can we tell how much error we have?  Well, you can either do a whole
>bunch of graduate-level mathematics, or you can run the same algorithm
>with different rounding modes.  If you always round towards positive
>infinity, you'll always be adding an error term- i.e. it'll be (x + e) and
>(y + e) with a result of x*y + 2e.  If you always round to negative
>infinity, you'll be subtracting the error term- (x - e) and (y - e) with a 
>result x*Y - 2e.  Comparing the results gives you both a good idea of 
>where the answer really is, and how big your errors are.
>  
>
Now I don't know that it is this simple.  I understand that by always 
rounding in one direction you are in effect  calculating a bound on the 
error in that direction but that will only work in the most simplistic 
of cases.  For instance if half of your values had their sign changed 
before being summed the sign of half  the associated error terms would 
also change.  The result would not give you any information about the 
effective lower/upper bound you thought you were calculating.  It seems 
to me that this approach is only applicable to problems which would be 
trivially solvable with graduate level mathematics and of absolutely no 
value in any interesting cases.

>Note, all is not necessarily lost.  If instead we have (x + e) and (y - e)  
>a little bit of quick algebra shows that the result is actually *more*
>accurate than either of the original inputs.  The errors tend to cancel
>instead of accumulate.  Numerical analysts call algorithms whose errors
>tend to cancel "stable", while algorithsm whose errors tend to accumulate
>are "unstable".  A good textbook on numerical methods will generally only 
>discuss stable algorithms- use 'em.
>  
>
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".  

Just to illustrate this the different "resolutions" the smallest 
positive double float which can be represented is  2.225 x 10^(-308) 
while the epsilon for doubles in general is 1.11 x 10^(-16) and can be 
thought of at the coarsest resolution.  That is a big difference and is 
because floats are packed in tightly around 0.

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.

>The effect can also be sort of acheived by using a random rounding method.  
>For each operation, the CPU picks one of it's rounding methods at random 
>to use (well, pseudo-random anyways).  Now, while this will in general 
>increase the accuracy of floating point programs, it also means that every 
>time you run a numerical program, even over the same data, you will get a 
>slightly different answer.  Newbies wondering why their program keeps 
>changing it's answer will probably rival newbies wondering why their stack 
>overflows pretty quickly.  Which is why I'm hesitant to advocate doing 
>this to all programs.  I still would like to be able to select the 
>rounding mode- either at compile time or at run time, which ever is 
>easier.  Actually, run time is probably easier to implement and more 
>usefull (I can switch rounding modes for different parts of the program). 
>  
>
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.

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.

-Eric

-------------------
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