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?