From mboxrd@z Thu Jan 1 00:00:00 1970 Received: (from majordomo@localhost) by pauillac.inria.fr (8.7.6/8.7.3) id XAA09546; Mon, 10 Nov 2003 23:10:21 +0100 (MET) X-Authentication-Warning: pauillac.inria.fr: majordomo set sender to owner-caml-list@pauillac.inria.fr using -f Received: from concorde.inria.fr (concorde.inria.fr [192.93.2.39]) by pauillac.inria.fr (8.7.6/8.7.3) with ESMTP id XAA09571 for ; Mon, 10 Nov 2003 23:10:20 +0100 (MET) Received: from herd.plethora.net (herd.plethora.net [205.166.146.1]) by concorde.inria.fr (8.11.1/8.11.1) with ESMTP id hAAMAJ112101 for ; Mon, 10 Nov 2003 23:10:19 +0100 (MET) Received: from bhurt.plethora.net (bhurt.plethora.net [205.166.146.49]) by herd.plethora.net (8.11.6/8.10.1) with ESMTP id hAAMAEC11453; Mon, 10 Nov 2003 16:10:15 -0600 (CST) Date: Mon, 10 Nov 2003 17:09:13 -0600 (CST) From: Brian Hurt X-X-Sender: bhurt@localhost.localdomain To: Eric Dahlman cc: Ocaml Mailing List Subject: Re: [Caml-list] Rounding mode In-Reply-To: <3FAFF6AD.4090009@atcorp.com> Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII X-Loop: caml-list@inria.fr X-Spam: no; 0.00; caml-list:01 rounding:01 simplistic:01 rounding:01 horrified:99 unexpected:01 replacements:01 numerically:01 numerically:01 analysts:99 figuring:01 differing:01 additive:01 statistic:99 randomized:01 Sender: owner-caml-list@pauillac.inria.fr Precedence: bulk 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