7

In Julia it appears that one picks up some error terms when doing finite differences using matrix multiplication versus shifts and addition.

julia> n = 1000
1000

julia> hessM = circshift(eye(n),-1) + circshift(eye(n),1) - 2* eye(n);

julia> function hess(x::Vector)
         return circshift(x,-1) + circshift(x,1) - 2*x
       end
hess (generic function with 1 method)

julia> x = rand(n);

julia> maximum(hessM * x - hess(x))
0.0

julia> x = rand(n)/300;

julia> maximum(hessM * x - hess(x))
8.673617379884035e-19

I understand that this is floating point error, since if I do x = rand(n) / 256 (or any power of 2 for that matter) the error goes away, but if we divide by something as simple as 3 the error crops up.

But where exactly is the floating point error coming in when running the above code? And why is it sensitive to the division?

Edit: in fact, I have localized the difference to 

julia> x = rand(n)/3;

julia> maximum(hessM * x + 2 * eye(n) * x - circshift(eye(n),-1) * x - circshift(eye(n),1)*x)
1.1102230246251565e-16

julia> x = rand(n);

julia> maximum(hessM * x + 2 * eye(n) * x - circshift(eye(n),-1) * x - circshift(eye(n),1)*x)
0.0

julia> x = rand(n) * 3;

julia> maximum(hessM * x + 2 * eye(n) * x - circshift(eye(n),-1) * x - circshift(eye(n),1)*x)
8.881784197001252e-16

What gives?

Willie Wong
  • 677
  • 2
  • 6
  • 15

1 Answers1

6

Edit (July 2021): it appears that the behavior will be changed as a side effect of the move of the default PRNG from Mersenne Twister to Xoshiro in the 1.7 release of Julia. See comments below.

It seems that this is tied to how Julia generates random numbers; I've opened a discussion on the Julia Language site. The current implementation of Julia's random number generator for the default range [0,1) for floats (in other words, calling simply rand()) always produces a 0 in the least significant bit for some reason or another (unlike MATLAB, for example). A side effect of this is that floating point errors are suppressed (since these are often errors coming from rounding/cancellation to the lowest bit). Multiplying/dividing by a power of 2 do not change the significand in the floating point representation, however multiplying/dividing by a non-power does. So generically after multiplying/dividing by some non-power, the least significant bit can be either 1 or zero, and now floating point error occurs as usually expected.

Willie Wong
  • 677
  • 2
  • 6
  • 15
  • 3
    Note that this will change soon since for Google Summer of Code there is a project to create RNG.jl which implements a bunch of RNGs, and use the best one to replace the current RNG. – Chris Rackauckas Jul 20 '16 at 16:03
  • @ChrisRackauckas Any idea what the status of this is now? Running the code in Julia 1.6 (replacing eye with an appropriate diagm call), I get the same results as in the question. Is there a place where there's either ongoing discussion or documentation, regarding this? – Sundar R Jul 01 '21 at 18:08
  • 1
    There's https://github.com/JuliaRandom/RandomNumbers.jl, https://github.com/JuliaRandom/StableRNGs.jl, and Julia Base just changed its RNG in Julia v1.7: https://github.com/JuliaLang/julia/pull/40546 – Chris Rackauckas Jul 01 '21 at 19:00
  • 1
    By chance, I just finished installed Julia 1.7-beta2. Testing in it, the error with rand(n)/300 returns the same value, and now the error with just rand(n) is also non-zero: 2.220446049250313e-16. Serendipitous that I came across this five-year old question right around the time the fix is to be released. :) – Sundar R Jul 02 '21 at 01:15