[Numpy-discussion] Different results from repeated calculation, part 2

[Alok Singhal]

On 14/08/08: 10:20, Keith Goodman wrote:
> A unit test is attached. It contains three tests:
> 
> In test1, I construct matrices x and y and then repeatedly calculate z
> = calc(x,y). The result z is the same every time. So this test passes.
> 
> In test2, I construct matrices x and y each time before calculating z
> = calc(x,y). Sometimes z is slightly different. But the x's test to be
> equal and so do the y's. This test fails (on Debian Lenny, Core 2 Duo,
> with libatlas3gf-sse2 but not with libatlas3gf-sse).
> 
> test3 is the same as test2 but I calculate z like this: z =
> calc(100*x,y) / (100 * 100). This test passes.
> 
> I get:
> 
> ======================================================================
> FAIL: repeatability #2
> ----------------------------------------------------------------------
> Traceback (most recent call last):
>   File "/home/[snip]/test/repeat_test.py", line 73, in test_repeat_2
>     self.assert_(result, msg)
> AssertionError: Max difference = 2.04946e-16

Could this be because of how the calculations are done?  If the
floating point numbers are stored in the cpu registers, in this case
(intel core duo), they are 80-bit values, whereas 'double' precision
is 64-bits.  Depending upon gcc's optimization settings, the amount of
automatic variables, etc., it is entirely possible that the numbers
are stored in registers only in some cases, and are in the RAM in
other cases.  Thus, in your tests, sometimes some numbers get stored
in the cpu registers, making the calculations with those values
different from the case if they were not stored in the registers.

See "The pitfalls of verifying floating-point computations" at
http://portal.acm.org/citation.cfm?doid=1353445.1353446 (or if that
needs subscription, you can download the PDF from
http://arxiv.org/abs/cs/0701192).  The paper has a lot of examples of
surprises like this.  Quote:

  We shall discuss the following myths, among others:

  ...

  - "Arithmetic operations are deterministic; that is, if I do z=x+y in
    two places in the same program and my program never touches x and y
    in the meantime, then the results should be the same."
  
  - A variant: "If x < 1 tests true at one point, then x < 1 stays true
    later if I never modify x."
  
  ...

-Alok

-- 
                                           *   *          
Alok Singhal                           *           *     *
http://www.astro.virginia.edu/~as8ca/                     
                                                   *    *