rounding errors with real numbers.

newer
trying to understand runProcess...

Matthias Fischmann

26 Feb 2006 26 Feb '06

7:33 a.m.

hi, I think this is the well-known issue of using real numbers in decimal representation on a machine that thinks binary, but I don't know what to do with it, and some of you maybe do. I want to shift+stretch a list of doubles into a given interval. example: | x1 = [2, 3, 4, 5, 10] | y1 = normInterval x1 0 1 | => y1 = [0.0,0.125,0.25,0.375,1.0] The function that does this looks something like this: | normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + lower) ps | where | oldLower = head ps | oldUpper = last ps | stretch = (upper - lower) / (oldUpper - oldLower) However, with bigger numbers I get rounding errors: | x2 = [0.0,1.9569471624266144e-3,5.870841487279843e-3,1.5655577299412915e-2,3.913894324853229e-2,9.393346379647749e-2,0.2191780821917808,0.5009784735812133,1.1272015655577299,2.504892367906066] | y2 = normInterval x2 0 1 | => y2 = [0.0,7.8125e-4,2.3437500000000003e-3,6.25e-3,1.5625000000000003e-2,3.7500000000000006e-2,8.750000000000001e-2,0.2,0.45000000000000007,0.9999999999999999] The solution that pops to my mind is very simple: | normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = repair (map (\ x -> (x - oldLower) * stretch + lower) ps) | where | oldLower = head ps | oldUpper = last ps | stretch = (upper - lower) / (oldUpper - oldLower) | | -- fix rounding error: | repair [i] = [upper] | repair (h:t) = h : repair t It works, but it's ugly. Is there a better way to do this? Thanks, Matthias

Attachments:

signature.asc (application/pgp-signature — 189 bytes)

Show replies by date

Chris Kuklewicz

26 Feb 26 Feb

8 a.m.

Your solution works, but is slightly wasteful with (repair) traversing the whole list again. Here is a slightly more efficient expression: -- Precondition: The first parameter (xs) is sorted (ascending) : -- assert (all (zipWith (<=) (xs, tail xs))) -- low' < high' -- low < high normInterval :: [Double] -> Double -> Double -> [Double] normInterval xs low high = let low' = head xs; high' = last xs; scale = (high-low)/(high'-low') middle = init (tail xs) in (low : [scale*(x-low')+low) | x <-middle])++[high] -- Chris

Matthias Fischmann

11:31 a.m.

New subject: [Solved] rounding errors with real numbers.

On Sun, Feb 26, 2006 at 01:00:54PM +0000, Chris Kuklewicz wrote:

...

To: Matthias Fischmann Cc: haskell-cafe@haskell.org From: Chris Kuklewicz Date: Sun, 26 Feb 2006 13:00:54 +0000 Subject: Re: [Haskell-cafe] rounding errors with real numbers.

Your solution works, but is slightly wasteful with (repair) traversing the whole list again. Here is a slightly more efficient expression:

-- Precondition: The first parameter (xs) is sorted (ascending) : -- assert (all (zipWith (<=) (xs, tail xs))) -- low' < high' -- low < high normInterval :: [Double] -> Double -> Double -> [Double] normInterval xs low high = let low' = head xs; high' = last xs; scale = (high-low)/(high'-low') middle = init (tail xs) in (low : [scale*(x-low')+low) | x <-middle])++[high]

hi chris, i like this code better, too. slightly better still, because of fewer typos: normInterval :: [Double] -> Double -> Double -> [Double] normInterval ps lower upper = assert check ([lower] ++ [ stretch * (p - oldLower) + lower | p <- middle ] ++ [upper]) where check = lower < upper && oldLower < oldUpper && (and (zipWith (<=) ps (tail ps))) oldLower = head ps oldUpper = last ps stretch = (upper - lower) / (oldUpper - oldLower) middle = init (tail ps) but then i'll just take this as the best solution. thanks, matthias

Lennart Augustsson

10:28 a.m.

Well, if you are relying on exact results from floating point arithmetic you're in trouble no matter what you do. I would just ignore the slight error and when finally printing the results do some rounding. Trying to fudge things is just going to bite you somewhere else. (BTW, I much prefer the name "floating point" to "real". I think the latter should be reserved for when you actually have real numbers.) -- Lennart Matthias Fischmann wrote:

...

hi,

I think this is the well-known issue of using real numbers in decimal representation on a machine that thinks binary, but I don't know what to do with it, and some of you maybe do.

I want to shift+stretch a list of doubles into a given interval. example:

| x1 = [2, 3, 4, 5, 10] | y1 = normInterval x1 0 1 | => y1 = [0.0,0.125,0.25,0.375,1.0]

The function that does this looks something like this:

| normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + lower) ps | where | oldLower = head ps | oldUpper = last ps | stretch = (upper - lower) / (oldUpper - oldLower)

However, with bigger numbers I get rounding errors:

| x2 = [0.0,1.9569471624266144e-3,5.870841487279843e-3,1.5655577299412915e-2,3.913894324853229e-2,9.393346379647749e-2,0.2191780821917808,0.5009784735812133,1.1272015655577299,2.504892367906066] | y2 = normInterval x2 0 1 | => y2 = [0.0,7.8125e-4,2.3437500000000003e-3,6.25e-3,1.5625000000000003e-2,3.7500000000000006e-2,8.750000000000001e-2,0.2,0.45000000000000007,0.9999999999999999]

The solution that pops to my mind is very simple:

| normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = repair (map (\ x -> (x - oldLower) * stretch + lower) ps) | where | oldLower = head ps | oldUpper = last ps | stretch = (upper - lower) / (oldUpper - oldLower) | | -- fix rounding error: | repair [i] = [upper] | repair (h:t) = h : repair t

It works, but it's ugly. Is there a better way to do this?

Thanks, Matthias

------------------------------------------------------------------------

_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

Jared Updike

11:40 a.m.

...

Well, if you are relying on exact results from floating point arithmetic you're in trouble no matter what you do.

As long as you don't do anything irrational (exp, sin, sqrt, etc.), you should be able to get away with using Rational. Number constants with decimals are not automatically constructors for floating point numbers in Haskell; they are exact (fromRational) until you make them Doubles or some other floating point value. If you have to use Doubles for other reasons (performance/memory, interfacing with other code, etc.) this won't help you... Just my 2c. Jared. -- http://www.updike.org/~jared/ reverse ")-:"

Brian Hulley

12:06 p.m.

Matthias Fischmann wrote:

...

| -- fix rounding error: | repair [i] = [upper] | repair (h:t) = h : repair t

Just to point out that this only fixes the last element of the list, so inputs like [1,2,10.8,10.8] would not be handled properly if you require the same input values to map to the same output values (I assume such inputs don't arise in the context you're using but in a general context the above wouldn't be a solution). Another thing is that when using floating point numbers, is there really any difference between 1.0 and 0.9999999 anyway? It's usually not recommended to ever test floats for equality since, depending on the architecture, the "same" float can end up being represented differently depending on what optimizations are happening eg an implementation could conceivably be making use of two different fp units if values are passed between different concurrently executing threads in a multi-processor or distributed processing environment... Regards, Brian.

Henning Thielemann

27 Feb 27 Feb

9:11 a.m.

On Sun, 26 Feb 2006, Matthias Fischmann wrote:

...

I think this is the well-known issue of using real numbers in decimal representation on a machine that thinks binary, but I don't know what to do with it, and some of you maybe do.

I want to shift+stretch a list of doubles into a given interval. example:

| x1 = [2, 3, 4, 5, 10] | y1 = normInterval x1 0 1 | => y1 = [0.0,0.125,0.25,0.375,1.0]

The function that does this looks something like this:

| normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + lower) ps

Is there --------------------------------------------------------------^ a cancellation problem? Maybe you should use a kind of convex combination, that is (x-oldLower)*a + (oldUpper-x)*b

Matthias Fischmann

1:51 p.m.

On Mon, Feb 27, 2006 at 03:11:35PM +0100, Henning Thielemann wrote:

...

To: Matthias Fischmann cc: haskell-cafe@haskell.org From: Henning Thielemann Date: Mon, 27 Feb 2006 15:11:35 +0100 (MET) Subject: Re: [Haskell-cafe] rounding errors with real numbers.

On Sun, 26 Feb 2006, Matthias Fischmann wrote:

...
I think this is the well-known issue of using real numbers in decimal representation on a machine that thinks binary, but I don't know what to do with it, and some of you maybe do.

I want to shift+stretch a list of doubles into a given interval. example:

| x1 = [2, 3, 4, 5, 10] | y1 = normInterval x1 0 1 | => y1 = [0.0,0.125,0.25,0.375,1.0]

The function that does this looks something like this:

| normInterval :: [Double] -> Double -> Double -> [Double] | normInterval ps lower upper = map (\ x -> (x - oldLower) * stretch + lower) ps

Is there --------------------------------------------------------------^ a cancellation problem?

what's a cancellation problem?

...

Maybe you should use a kind of convex combination, that is

(x-oldLower)*a + (oldUpper-x)*b

i don't quite understand this either. is 'x' an old element in my input list and your expression is the corresponding new element? then how does the resulting curve relate to the original one? does this keep the ratios between distances between elements in the list intact? (this is the property that i am interested in.) but it sounds intriguing. perhaps i should play with this a little and find out myself. thanks, matthias

Henning Thielemann

2:48 p.m.

On Mon, 27 Feb 2006, Matthias Fischmann wrote:

...

On Mon, Feb 27, 2006 at 03:11:35PM +0100, Henning Thielemann wrote:

...
Is there a cancellation problem?

what's a cancellation problem?

zu deutsch: Auslöschung cancellation happens for instance here: 1 + 1e-50 - 1 == 0

...

...
Maybe you should use a kind of convex combination, that is

(x-oldLower)*a + (oldUpper-x)*b

i don't quite understand this either.

You have to adjust 'a' and 'b' in order to obtain 0 on x==oldLower and 1 on x==oldUpper. The expression still depends linearly on x (plus an absolute term). Maybe it reduces the cancellations if the lower and the upper bound differ much in magnitude.

Matthias Fischmann

2 Mar 2 Mar

12:42 p.m.

...

cancellation happens for instance here: 1 + 1e-50 - 1 == 0

the function again (in the wasteful original form, for clarity). where do you think cancellation may take place? isn't what you call canellation a generic rounding error? ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ normInterval :: [Double] -> Double -> Double -> [Double] normInterval ps lower upper = repair (map (\ x -> (x - oldLower) * stretch + lower) ps) where oldLower = head ps oldUpper = last ps stretch = (upper - lower) / (oldUpper - oldLower) -- fix rounding error: repair [i] = [upper] repair (h:t) = h : repair t ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

...

...
...
(x-oldLower)*a + (oldUpper-x)*b

i got this into my head though. neat. thanks! i will rewrite the function right now. cheers, matthias

Henning Thielemann

3 Mar 3 Mar

7:29 a.m.

On Thu, 2 Mar 2006, Matthias Fischmann wrote:

...

...
cancellation happens for instance here: 1 + 1e-50 - 1 == 0

the function again (in the wasteful original form, for clarity). where do you think cancellation may take place? isn't what you call canellation a generic rounding error?

Cancellation is a special kind of rounding error. Rounding errors appear everywhere, in (*), sin, exp and so on, but cancellations only arise on differences. They are especially bad, because as the example above shows, even if all numbers are given in double precision in the computation a+b-a, no digit of the result is correct, that is 100% rounding error! The danger of cancellation is everywhere where you subtract numbers of similar magnitude. In your example: If lower=-100000, upper=1, x approximately oldUpper then you get the effect at the '+' in ((x - oldLower) * stretch + lower).

Matthias Fischmann

9:02 a.m.

On Fri, Mar 03, 2006 at 01:29:06PM +0100, Henning Thielemann wrote:

...

To: Matthias Fischmann cc: haskell-cafe@haskell.org From: Henning Thielemann Date: Fri, 3 Mar 2006 13:29:06 +0100 (MET) Subject: Re: [Haskell-cafe] rounding errors with real numbers.

On Thu, 2 Mar 2006, Matthias Fischmann wrote:

...
...
cancellation happens for instance here: 1 + 1e-50 - 1 == 0

the function again (in the wasteful original form, for clarity). where do you think cancellation may take place? isn't what you call canellation a generic rounding error?

Cancellation is a special kind of rounding error. Rounding errors appear everywhere, in (*), sin, exp and so on, but cancellations only arise on differences. They are especially bad, because as the example above shows, even if all numbers are given in double precision in the computation a+b-a, no digit of the result is correct, that is 100% rounding error! The danger of cancellation is everywhere where you subtract numbers of similar magnitude.

1 + epsilon - 1 == epsilon, which is zero except for a very small rounding error somewhere deep in the e-minus-somethings. how is the error getting worse than that, for which numbers? m.

Henning Thielemann

9:19 a.m.

On Fri, 3 Mar 2006, Matthias Fischmann wrote:

...

1 + epsilon - 1 == epsilon, which is zero except for a very small rounding error somewhere deep in the e-minus-somethings. how is the error getting worse than that, for which numbers?

I meant the relative error. epsilon should be the result, but the computer says 0, so the absolute error is epsilon and the relative error is 1, that is the number of reliable digits in the mantissa of the result is 0.

Matthias Fischmann

9:28 a.m.

ok, now i am intimidated enough to give up. (-:. thanks for trying, though. m. On Fri, Mar 03, 2006 at 03:19:44PM +0100, Henning Thielemann wrote:

...

To: Matthias Fischmann cc: haskell-cafe@haskell.org From: Henning Thielemann Date: Fri, 3 Mar 2006 15:19:44 +0100 (MET) Subject: Re: [Haskell-cafe] rounding errors with real numbers.

On Fri, 3 Mar 2006, Matthias Fischmann wrote:

...
1 + epsilon - 1 == epsilon, which is zero except for a very small rounding error somewhere deep in the e-minus-somethings. how is the error getting worse than that, for which numbers?

I meant the relative error. epsilon should be the result, but the computer says 0, so the absolute error is epsilon and the relative error is 1, that is the number of reliable digits in the mantissa of the result is 0.

Ben Rudiak-Gould

28 Feb 28 Feb

9:31 a.m.

Henning Thielemann wrote:

...

Maybe you should use a kind of convex combination, that is

(x-oldLower)*a + (oldUpper-x)*b

Maybe lower*(1-z) + upper*z where z = (x-oldLower) / (oldUpper-oldLower). I think this will always map oldLower and oldUpper to lower and upper exactly, but I'm not sure it's monotonic. -- Ben

7018

Age (days ago)

7023

Last active (days ago)

List overview

Download

14 comments

7 participants

participants (7)

Ben Rudiak-Gould
Brian Hulley
Chris Kuklewicz
Henning Thielemann
Jared Updike
Lennart Augustsson
Matthias Fischmann