# Precision issue: (A.+R.-A.-R) do not vanish ?!
**URL:** https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615
**Category:** Performance
**Tags:** question
**Created:** [July 23, 2018, 8:24pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615 "2018-07-23T20:24:53Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)
#### Post date: [July 23, 2018, 8:24pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/1 "2018-07-23T20:24:53Z")
```julia
N = 1023
A = rand(1000:20000,(N,N));
R = rand(ComplexF64,N,N);
println(sum(abs.(A.+R.-A.-R)))
println(sum(abs.(R.+A.-R.-A)))
```
On my machine it gives
```julia
4.53665288446814e-7
5.0590642786119133e-11
```
What is happening?
---
### Author:  [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)
#### Post date: [July 23, 2018, 8:29pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/2 "2018-07-23T20:29:24Z")
Machine \epsilon due to order of operations?
For
```julia
sum(abs.(R .- R .+ A .- A))
```
I indeed get `0.0`.
---
### Author:  [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)
#### Post date: [July 23, 2018, 8:34pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/3 "2018-07-23T20:34:38Z")
Yes, `R.-R` and `A.-A` actually yield zero matrices.
Python3 with numpy library also have this issue.
I came across this issue by testing my code and I got desperate ! I don’t know how many such precision issues are there due to my undisciplined coding habit !!!
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [July 23, 2018, 8:35pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/4 "2018-07-23T20:35:24Z")
> [@Lian\_Yunlong](#):
>
> Python3 with numpy library also have this issue.
Floating-point arithmetic itself has this issue. Unless the two quantities happen to have exactly the same exponents, there’s going to be some precision loss when adding.
---
### Author:  [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)
#### Post date: [July 23, 2018, 8:37pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/5 "2018-07-23T20:37:29Z")
But are there generally some principles / common practices that minimizes this uncertainty?
---
### Author:  [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)
#### Post date: [July 23, 2018, 8:39pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/6 "2018-07-23T20:39:09Z")
See [here](https://en.wikipedia.org/wiki/Machine_epsilon) for an introduction to the concept.
Also see the first section of [Numerical Recipes](http://e-maxx.ru/bookz/files/numerical_recipes.pdf).
---
### Author:  [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)
#### Post date: [July 23, 2018, 8:40pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/7 "2018-07-23T20:40:12Z")
Oh ! That is quite an education …
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [July 23, 2018, 8:42pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/8 "2018-07-23T20:42:39Z")
It’s a whole field of computational study. In general, don’t expect floating-point arithmetic to be precisely associative. If you expect something to be zero, don’t be too surprised when it’s not quite zero (avoid comparisons to zero whenever possible). It’s better to check approximate equality of two non-zero quantities than to check that their difference is exactly zero. The non-zero quantities give a decent sense of scale and thus when an approximate equality is ok; by the time you’ve subtracted them, you’ve lost all sense of scale so it’s hard to tell if a non-zero residual is acceptable. At the end of a long computation, if half of your bits are correct, you’re doing pretty well. This is why having 64 bits of precision is so important—so that you have a decent chance that the first dozen of those bits are actually right.
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [July 23, 2018, 8:43pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/9 "2018-07-23T20:43:52Z")
But beware of copying any of the recipes in _Numerical Recipes_—they are often quite flawed (there may be legal and copyright issues as well). The general principles it explains are good to understand though.
---
### Author:  [@geoalchimista](https://discourse.julialang.org/u/geoalchimista)
#### Post date: [July 23, 2018, 8:44pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/10 "2018-07-23T20:44:55Z")
> [@Lian\_Yunlong](#):
>
> R = rand(ComplexF64,N,N);
~~Also, your Complex64 is made of Float32 types. That is the primary reason you get a result bigger than `eps(Float64)`.~~
Edit: My apologies, didn’t see the 0.7.0 change. Then indeed this issue may be a roundoff issue, because `A` is an `Int64` array populated with values between 1000 and 20000, while the modulus of any element in `R` is smaller than one.
Edit2: If you do `println(sum(abs.(A.-A.+R.-R)))`, you’d always get `0.0`, which confirms that _floating point math is not associative_.
Edit3: On the other hand, you were summing up 1023 \* 1023 elements, so the total roundoff error of `4.5e-7` is not surprising.
---
### Author:  [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)
#### Post date: [July 23, 2018, 8:45pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/11 "2018-07-23T20:45:51Z")
> [@StefanKarpinski](#):
>
> that their difference is exactly zero. The quantities give a decent sense of when an approximate equality is o
I never copy recipes. I write my own code and test, and … fail 😃
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [July 23, 2018, 8:47pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/12 "2018-07-23T20:47:04Z")
The `ComplexF64` type is actually `Complex{Float64}`; we changed it from the old `Complex64` which was a 64-bit complex type with `Float32` parts because this kept tripping people up.
---
### Author:  [@geoalchimista](https://discourse.julialang.org/u/geoalchimista)
#### Post date: [July 23, 2018, 8:50pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/13 "2018-07-23T20:50:16Z")
I see. My apologies for not catching up with the 0.7.0 changes.
---
### Author:  [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)
#### Post date: [July 23, 2018, 8:50pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/14 "2018-07-23T20:50:27Z")
> [@StefanKarpinski](#):
>
> But beware of copying any of the recipes in _Numerical Recipes_ —they are often quite flawed. The general principles it explains are good to understand though.
I was referring to the explanations in the text, not the code. There are dubious copyright issues with the code itself. I knew one of the authors, so fortunately I got the catharsis of complaining about this in person. 😄
> [@geoalchimista](#):
>
> Also, your Complex64 is made of Float32 types.
He was using `ComplexF64` (note the `F`) which was introduced in 0.7 as an alias for `Complex{Float64}`.
Indeed you should _ **never** _ do equality comparisons with floating point values. This principle was deeply ingrained in me from an early age, but lately, as I’ve come to understand a bit better how code works especially in Julia, I’ve started to sometimes when I know it should be ok (such as `sum(R .- R)` in the example above). I should probably immediately abandon that habit though, it’s bad practice. Never do equality checks with floats.
---
### Author:  [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)
#### Post date: [July 23, 2018, 8:52pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/15 "2018-07-23T20:52:23Z")
Stefan, may I ask about how to defend my work if some referee says "Your work is ultimately meaningless because the uncontrolled / unknown precision issues in JuliaLang could easily ruin the calculation? " Are there any published papers to cite if I use Julia for scientific computing? do I have to cite anything?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 23, 2018, 8:59pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/16 "2018-07-23T20:59:15Z")
> [@Lian\_Yunlong](#):
>
> Stefan, may I ask about how to defend my work if some referee says "Your work is ultimately meaningless because the uncontrolled / unknown precision issues in JuliaLang could easily ruin the calculation? " Are there any published papers to cite if I use Julia for scientific computing? do I have to cite anything?
Julia is just using floating point arithmetic. Fortran, C, C++, Python, MATLAB, R, etc. all have exactly the same issue with exactly the same precision. This is built into your computer. Everyone knows about it. If a reviewer says you have precision issues, it would be due to using a numerically unstable algorithm and not anything to do with a programming language.
---
### Author:  [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)
#### Post date: [July 23, 2018, 9:11pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/17 "2018-07-23T21:11:31Z")
I do think the same. Gonna rethink all the calculations I’ve done tonight …
By the way, the problem actually came from my effort to test a column reduction function I wrote with `CUDAnative.jl`. The function sum the first axis of `array2d_in` and write results to `array1d_out`. The code is ugly because I explicitly use `Uint32`, in order to make the code portable. Many GPUs prefer 32bit integers.
```julia
using CUDAdrv
using CUDAnative
Curry = CuArray
#NOTE [:,b] ==> [b]
function ReduceC!(
dims::NTuple{2,UInt32},
array1d_out,
array2d_in)
# -------------------------------------------------
(dim_to_reduce::UInt32, dim1::UInt32) = dims
tid::UInt32 = UInt32(threadIdx().x)
bid::UInt32 = UInt32(blockIdx().x)
sdata = @cuStaticSharedMem(ComplexF64,1024) #NOTE doubled shared memory
sdata[tid] = ComplexF64(0)
if bid<=dim1
# ------------------------------ PART 1 ------------------------------
# accumulate on the shared memory in each block
tl::UInt32 = tid
while tl <= dim_to_reduce
@inbounds sdata[tid] += array2d_in[tl,bid];
tl += 0x00000400
end
# ------------------------------ PART 2 ------------------------------
# folding the shared memory
sync_threads()
if tid<=0x00000200
@inbounds sdata[tid] += sdata[tid + 0x00000200]
end
sync_threads()
if tid<=0x00000100
@inbounds sdata[tid] += sdata[tid + 0x00000100]
end
sync_threads()
if tid<=0x00000080
@inbounds sdata[tid] += sdata[tid + 0x00000080]
end
sync_threads()
if tid<=0x00000040
@inbounds sdata[tid] += sdata[tid + 0x00000040]
end
sync_threads()
if tid<=0x00000020
@inbounds sdata[tid] += sdata[tid + 0x00000020]; sync_warp()
@inbounds sdata[tid] += sdata[tid + 0x00000010]; sync_warp()
@inbounds sdata[tid] += sdata[tid + 0x00000008]; sync_warp()
@inbounds sdata[tid] += sdata[tid + 0x00000004]; sync_warp()
@inbounds sdata[tid] += sdata[tid + 0x00000002]; sync_warp()
@inbounds sdata[tid] += sdata[tid + 0x00000001]; sync_warp()
end
sync_threads()
if tid==0x00000001
@inbounds array1d_out[bid] = sdata[0x00000001]
end
end
end
```
To test the code, I wote
```julia
N = 1023
A = rand(1000:20000,(N,N))
R = rand(ComplexF64,N,N)
cuAF = Curry( ComplexF64.(A) .+ R )
cuBF = Curry( fill(ComplexF64(0),N) )
@cuda threads=1024 blocks=N ReduceC!(UInt32.((N,N)),cuBF,cuAF)
REDUCTION_CPU = sum(Array(cuAF),dims=1)[1,:]
REDUCTION_GPU = Array(cuBF)
println( sum(abs.(REDUCTION_CPU.-REDUCTION_GPU)) )
```
The result (different from `0.0`) confused me and I suspected that GPU handles `Float64` differently … I am still sorting out from all online materials.
I am not very familiar with GPU programming, and unable to source the numeric errors …
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 23, 2018, 10:23pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/18 "2018-07-23T22:23:15Z")
> [@Lian\_Yunlong](#):
>
> I am not very familiar with GPU programming, and unable to source the numeric errors …
GPUs compute in parallel. That plus the lack of associativity means it can be far harder to control such floating point errors and they could change each run due to how things parallelize. In this case, the reduction ordering can be random due to the parallelism, and so exactly the difference from this test case causes numerical differences (i.e. the next value summed cannot be assumed to be the next value in the array, but changing order changes the floating point errors).
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [July 23, 2018, 10:30pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/19 "2018-07-23T22:30:37Z")
Julia is generally more careful with numerical stability, gives the user more control over it, and offers many more ways to validate numerical accuracy of computations than ANY other programming language. If you want to investigate the numerical stability of your algorithm, there are many high quality options in Julia, including:
- [big floats](https://docs.julialang.org/en/latest/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic-1)
- [double floats](https://github.com/JuliaMath/DoubleFloats.jl)
- [interval arithmetic](https://github.com/JuliaIntervals/IntervalArithmetic.jl)
- [validated numerics](https://github.com/JuliaIntervals/ValidatedNumerics.jl)
None of these are an option in Fortran, C, C++, Python, Matlab or R. The real question is: if you’re not using Julia, then how do you know that your algorithm is accurate and numerically stable?
---
### Author:  [@ExpandingMan](https://discourse.julialang.org/u/ExpandingMan)
#### Post date: [July 23, 2018, 11:13pm UTC](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615/20 "2018-07-23T23:13:00Z")
> [@Lian\_Yunlong](#):
>
> may I ask about how to defend my work if some referee says "Your work is ultimately meaningless because the uncontrolled / unknown precision issues in JuliaLang could easily ruin the calculation? "
Woah, this is just an uncertainty much like the others. In publications it should be quoted along with your other uncertainties. Like other uncertainties, it is bounded. I don’t know what journal this is, but any referee with any experience whatsoever in working on computational problems should understand that this uncertainty exists. In practice, usually uncertainty due to floating point arithmetic with `Float64` is negligible relative to statistical and systematic uncertainties, even in most fields of physics.
[Next page](https://discourse.julialang.org/t/precision-issue-a-r-a-r-do-not-vanish/12615.md?page=2)