rounding errors with real numbers.

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
15 messages Options
Reply | Threaded
Open this post in threaded view
|

rounding errors with real numbers.

Matthias Fischmann


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
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe

signature.asc (196 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

haskell-2
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
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Lennart Augustsson
In reply to this post by Matthias Fischmann
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
> [hidden email]
> http://www.haskell.org/mailman/listinfo/haskell-cafe

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

[Solved] rounding errors with real numbers.

Matthias Fischmann
In reply to this post by haskell-2

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

> To: Matthias Fischmann <[hidden email]>
> Cc: [hidden email]
> From: Chris Kuklewicz <[hidden email]>
> 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

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe

signature.asc (196 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Jared Updike
In reply to this post by Lennart Augustsson
> 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 ")-:"
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Brian Hulley
In reply to this post by Matthias Fischmann
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.

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Henning Thielemann
In reply to this post by Matthias Fischmann

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
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Matthias Fischmann

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

> To: Matthias Fischmann <[hidden email]>
> cc: [hidden email]
> From: Henning Thielemann <[hidden email]>
> 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

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe

signature.asc (196 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Henning Thielemann

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.
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Ben Rudiak-Gould
In reply to this post by Henning Thielemann
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

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Matthias Fischmann
In reply to this post by Henning Thielemann

> 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

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe

signature.asc (196 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Henning Thielemann

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

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Matthias Fischmann


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

> To: Matthias Fischmann <[hidden email]>
> cc: [hidden email]
> From: Henning Thielemann <[hidden email]>
> 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.

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe

signature.asc (196 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Henning Thielemann

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.

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: rounding errors with real numbers.

Matthias Fischmann

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 <[hidden email]>
> cc: [hidden email]
> From: Henning Thielemann <[hidden email]>
> 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.

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe

signature.asc (196 bytes) Download Attachment