# Huge computation error in Float32?

**URL:** <https://discourse.julialang.org/t/huge-computation-error-in-float32/106746>\
**Category:** Numerics\
**Created:** [November 26, 2023, 2:47pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746 "2023-11-26T14:47:50Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [November 26, 2023, 2:47pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/1 "2023-11-26T14:47:50Z")

</div>

the following shows that `s1` was hugely mis-computed. It happens when doing calculations in `Float32`. Everything is fine in `Float64`.

```julia
Random.seed!(1);
n = 1000000;
a = rand(Float32, n);
b = rand(Float32, 100, n);
s1 = 0f0;
s2 = 0f0;
si = zeros(Float32, n);
for i in 1:n
    for j in 1:100
        s1 += a[i] * b[j, i]
        si[i] += a[i] * b[j, i]
    end
    s2 += si[i]
end

julia> s1
1.6777216f7
julia> s2
2.4987638f7

aF64 = Float64.(a);
bF64 = Float64.(b);
s1F64 = 0.0;
s2F64 = 0.0;
siF64 = zeros(n);
for i in 1:n
    for j in 1:100
        s1F64 += aF64[i] * bF64[j, i]
        siF64[i] += aF64[i] * bF64[j, i]
    end
    s2F64 += siF64[i]
end

julia> s1F64
2.4989047580531508e7
julia> s2F64
2.4989047580524024e7

```

Is it a big problem??? using v1.9.1. Thanks for your attention.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [November 26, 2023, 2:54pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/2 "2023-11-26T14:54:33Z")

</div>

> [@tomtom](#):
>
> Is it a big problem???

It’s definitely a problem if you were expecting an accurate answer. But it’s also not unexpected. Everything you’re seeing here is documented, expected rounding behaviour of floating point arithmetic.

---

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [November 26, 2023, 3:01pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/3 "2023-11-26T15:01:06Z")

</div>

It’s a _ **big** _ problem if u see the difference between s1 and s2.

Given all the rounding errors in floating point calculations, s1 and s2 should be identical, so at least very close. Now they aren’t.

The results in Float64 are much more reasonable. So why this _ **huge** _ discrepancy happens only in Float32???

Or, in other words, is there any explanation that the computation of s2 is so much more accurate than s1???

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [November 26, 2023, 3:05pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/4 "2023-11-26T15:05:09Z")

</div>

You’re accumulating error at every step of the calculation, and your calculation has on the other of `1000000 * 100` steps. That’s a **lot** of accumulated error, and `Float32` has very low precision.

---

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [November 26, 2023, 3:06pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/5 "2023-11-26T15:06:39Z")

</div>

The problem is s2 is just getting the _ **same** _ times of error accumulation! Why the difference?

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [November 26, 2023, 3:16pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/6 "2023-11-26T15:16:13Z")

</div>

it’s because `s2` is accumulated with a more numerically stable pattern. `s1` just grows and grows, and the larger it becomes, the larger the rounding errors will be.

Once `s1` is of order `1f7`, the epsilon between floats is of order `1f0`:

```julia
julia> eps(1f7)
1.0f0

```

I strongly recommend reading about floating point error. It’s a subtle and complicated topic.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [November 26, 2023, 3:16pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/7 "2023-11-26T15:16:56Z")

</div>

I assume Mason is right: this is an example of accumulated rounding error. If Mason’s suggestion is correct, this has _nothing_ to do with Julia and whether it is v1.9.1 or not: you would then get similar error in MATLAB or Python or whatever, as long as you use Float32 and use comparable algorithms.

---

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [November 26, 2023, 3:18pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/8 "2023-11-26T15:18:59Z")

</div>

thanks

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [November 26, 2023, 3:20pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/9 "2023-11-26T15:20:51Z")

</div>

Here’s a perhaps simpler example that might help you understand:

```julia
julia> let v = rand(Float32, 1000000 * 100)
           s1 = 0f0
           s2 = 0f0
           for vi ∈ v
               s1 += vi
           end
           chunks = Iterators.partition(v, 100)
           for chunk ∈ chunks
               s3 = 0f0
               for vi ∈ chunk
                   s3 += vi
               end
               s2 += s3
           end
           @info "" s1 s2 sum(v)
       end
┌ Info: 
│ s1 = 1.6777216f7
│ s2 = 5.0003064f7
└ sum(v) = 5.0002704f7

```

Splitting the sum into chunks that we sum over individually is a good way to reduce errors.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [November 26, 2023, 3:32pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/10 "2023-11-26T15:32:00Z")

</div>

Maybe something I should spell out even more clearly @tomtom, is that once your accumulator hits `1f7`, adding a number less than `1` won’t do anything to it:

```julia
julia> v1 = 1f7
1.0f7

julia> (v1 + 0.5f0) === 1f7
true

```

This is a catastrophic loss of precision. The benefit of summing up the numbers in chunks is that you’re more likely to add up numbers that are similar in magnitude, and thus avoid situations where you round away the entire answer.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [November 26, 2023, 3:42pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/11 "2023-11-26T15:42:45Z")

</div>

Just for reference, adding the obligatory paper here: [What every computer scientist should know about floating-point arithmetic](https://dl.acm.org/doi/10.1145/103162.103163)  
Mathematically, floating-point arithmetic is broken in being neither associative nor commutative.

---

<div class="post-metadata">

**Author:** ![and](https://avatars.discourse-cdn.com/v4/letter/a/2acd7d/32.png) [@and](https://discourse.julialang.org/u/and)\
**Post date:** [November 26, 2023, 5:53pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/12 "2023-11-26T17:53:39Z")

</div>

Here the topic is explained in a very accessible and gentle manner using Julia: [Fundamentals of Numerical Computation. Floating-point numbers](https://tobydriscoll.net/fnc-julia/intro/floating-point.html)

---

<div class="post-metadata">

**Author:** ![arch.d.robison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arch.d.robison/32/17699_2.png) [@arch.d.robison](https://discourse.julialang.org/u/arch.d.robison)\
**Post date:** [November 27, 2023, 2:37pm UTC](https://discourse.julialang.org/t/huge-computation-error-in-float32/106746/13 "2023-11-27T14:37:39Z")

</div>

> [@bertschi](#):
>
> floating-point arithmetic is broken in being neither associative nor commutative

Clarification: Except when both operands are NaN, IEEE floating-point addition and multiplication _are commutative_. The commutativity is required because the IEEE standard requires that the result of addition or multiplication be the representable value closest to the exact result. In other words, the result might not be exact, but the order of the two operands does not matter.