Round-off error?

68 views
Skip to first unread message

D Yong

unread,
Oct 30, 2013, 11:35:53 PM10/30/13
to
Hi everyone,

The following lines of code are stumping me.

ClearAll[y]
dt = 0.15;
Do[y[i*dt] = "hi there", {i, 0, 20}]
Table[{i, y[(i + 1)*dt], y[i*dt + dt]}, {i, 0, 19}]

Here is the output:

0 hi there hi there
1 hi there hi there
2 hi there hi there
3 hi there hi there
4 hi there hi there
5 hi there y(0.9)
6 hi there y(1.05)
7 hi there hi there
8 hi there hi there
9 hi there hi there
10 hi there hi there
11 hi there hi there
12 hi there y(1.95)
13 hi there hi there
14 hi there hi there
15 hi there hi there
16 hi there hi there
17 hi there hi there
18 hi there y(2.85)
19 hi there hi there


It seems that y[i*dt+dt] does not match y[(i+1)*dt], perhaps due to numerical round-off? If you switch to dt=15/100, the problem goes away. You might also try different values for dt to see what happens.


If it is numerical round-off, then why does this evaluate to true?

6*.15 == 5*.15 + .15
True


Thanks in advance for your input!
Darryl

Itai Seggev

unread,
Nov 2, 2013, 2:20:54 AM11/2/13
to
The first thing to do whenever Mathematica does something unexpected is to look
at the FullForm.

In[1]:= 6*.15 //FullForm
Out[1]//FullForm= 0.8999999999999999`

In[2]:= 5*.15 + .15 //FullForm
Out[2]//FullForm= 0.9`

So clearly these are two different numbers. And the answer to your last
question is easy--Equal by default ignores the last seven bits in comparison,
pricesly to avoid too many False due to machine precision error. You can
control this using the variable Internal`$EqualTolerance:

In[3]:= Block[{Internal`$EqualTolerance = 0}, 6*.15 == 5*.15 + .15]

Out[3]= False

Note that this variable is specified using decimal digits, which is why its
default value is 2.10721 == N[Log[10, 2^7].

It's also easy to understand why it goes away with 15/100--you're now using
exact numbers and you end up with the exact result 9/10 in both cases.

However, you have an added twist, which is that you're using a function
definition / pattern. (And this is something I didn't fully understand until I
started poking at this, so thanks for the good question!) All patterns,
including literal patterns like 0.9`, are hashed for efficiency reasons. On
the other hand Equal is no longer strickly transitive (since two numbers which
different from a third by 6 bits need not differ from each other by less than
7). So machine numbers are effectively binned in the hash process to provide
something fast and transitive. You are suffereing because sometimes the two
computations are in the same bin and sometimes they are not.

In[31]:= Internal`HashSameQ[.15*6, .15*5 + .15]
Out[31]= False

In[32]:= Internal`HashSameQ[.15*5, .15*4 + .15]
Out[32]= True

This binning only affects patterns/definitions tied to specific approximate
numbers. If you want to avoid the binning, you could rewrite your definition
in the form

y[x_] /; x == num = "Hi there"

This will make sure that two numbers which are Equal trigger the same
definition, but at performance cost. How to balance those is dependent on your
needs and preferences.

Here I show that this change actually works. Notice that I'm using With in the
Do loop to make sure that num is evaluated before the definition of y is
created, which is necessary when generating definitions programmatically.

In[27]:= ClearAll[y]
dt=0.15;
Do[With[{num=dt*i}, y[x_] /; x==num="hi there"],{i,0,20}]
Table[{i,y[(i+1)*dt],y[i*dt+dt]},{i,0,19}]//TableForm
Out[30]//TableForm=
0 hi there hi there
1 hi there hi there
2 hi there hi there
3 hi there hi there
4 hi there hi there
5 hi there hi there
6 hi there hi there
7 hi there hi there
8 hi there hi there
9 hi there hi there
10 hi there hi there
11 hi there hi there
12 hi there hi there
13 hi there hi there
14 hi there hi there
15 hi there hi there
16 hi there hi there
17 hi there hi there
18 hi there hi there
19 hi there hi there


--
Itai

Itai Seggev
Mathematica Algorithms R&D

Bill Rowe

unread,
Nov 2, 2013, 2:21:15 AM11/2/13
to
On 10/30/13 at 11:47 PM, darry...@gmail.com (D Yong) wrote:

>The following lines of code are stumping me.

>ClearAll[y] dt = 0.15; Do[y[i*dt] = "hi there", {i, 0, 20}]
>Table[{i, y[(i + 1)*dt], y[i*dt + dt]}, {i, 0, 19}]

>Here is the output:

<output snipped>

>It seems that y[i*dt+dt] does not match y[(i+1)*dt], perhaps due to
>numerical round-off?

The problem is machine arithmetic, the order in which things are
done and the fact 0.15 can not be represented exactly in a
finite number of binary digits. You can verify the issue is
caused by this by doing:

In[3]:= dt = .15;
n = 5;

In[5]:= RealDigits[(1 + n) dt] === RealDigits[n dt + dt]

Out[5]= False

>If you switch to dt=15/100, the problem goes away. You might also try
>different values for dt to see what happens.

Right. By using 15/100 you are using exact arithmetic rather
than machine precision values and there is no round off.

>If it is numerical round-off, then why does this evaluate to true?

>6*.15 == 5*.15 + .15 True

Look at the documentation for Equal and note under details

Approximate numbers with machine precision or higher are
considered equal if they differ in at most their last seven
binary digits (roughly their last two decimal digits).

And

In[6]:= Pick[Range[Length@RealDigits[(1 + n) dt, 2][[1]]],
Equal @@@
Transpose[{RealDigits[(n + 1) dt, 2][[1]],
RealDigits[n dt + dt, 2][[1]]}], False]

Out[6]= {53}

shows the two machine precision numbers differ by 1 bit and
consequently are considered equal by Equal as documented.

This type of issue is inherent in floating point arithmetic done
on any modern computer.


Richard Fateman

unread,
Nov 4, 2013, 11:11:24 PM11/4/13
to
On 11/1/2013 11:21 PM, Bill Rowe wrote:

>
> shows the two machine precision numbers differ by 1 bit and
> consequently are considered equal by Equal as documented.
>
> This type of issue is inherent in floating point arithmetic done
> on any modern computer.
>
This is not at all true. Mathematica chooses to implement a fuzzy
version of equality that does not correspond to any hardware
floating-point operation on any modern computer that I am aware of.


Mathematica's definition of Equal is a misfeature that leads to
many contradictions to ordinary arithmetic results. This has
been pointed out numerous times and presumably will continue to
provoke questions indefinitely.

RJF



>


Bill Rowe

unread,
Nov 6, 2013, 12:19:44 AM11/6/13
to
On 11/4/13 at 11:17 PM, fat...@cs.berkeley.edu (Richard Fateman)
wrote:

>On 11/1/2013 11:21 PM, Bill Rowe wrote:

>>shows the two machine precision numbers differ by 1 bit and
>>consequently are considered equal by Equal as documented.

>>This type of issue is inherent in floating point arithmetic done on
>>any modern computer.

>This is not at all true. Mathematica chooses to implement a fuzzy
>version of equality that does not correspond to any hardware
>floating-point operation on any modern computer that I am aware of.

Just to clarify my comment. My intended meaning was (n+1)*dt
will differ from n*dt+dt for some values of n and dt, a
characteristic of any modern computer and a fundamental
characteristic of floating point arithmetic. I had not intended
to imply Mathematica's design for Equal is universal.


Reply all
Reply to author
Forward
0 new messages