# Imprecision of mean over iterators of large vectors

**URL:** <https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759>\
**Category:** General Usage\
**Created:** [December 17, 2018, 6:43pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759 "2018-12-17T18:43:34Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![cstjean](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cstjean/32/1444_2.png) [@cstjean](https://discourse.julialang.org/u/cstjean)\
**Post date:** [December 17, 2018, 6:43pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/1 "2018-12-17T18:43:35Z")

</div>

To wit:

```julia
julia> N = 100000000; aa = rand(Float32, N);

julia> mean((x for x in aa))
0.16777216f0

julia> mean(aa)
0.500059f0

```

I know what’s happening: Float32 has limited accuracy, and `mean` over an iterator computes a “naive” sum, which saturates at some point. However, it caught me by surprise in my real-world scenario:

```julia
julia> mean(skipmissing(Umat))
1.0638367f0 V

julia> mean(filter(!ismissing, Umat))
3.1320891f0 V

```

That’s a bit scary. Is that just Something I Must Watch Out For? It feels like there should be a reasonable algorithm that gives good accuracy without sacrificing speed. Something like

```julia
function sum(iter)
    total = zero(eltype(iter))
    minitotal = zero(eltype(iter))
    i = 0
    for x in iter
        if (i+=1) % 100000 == 0
            total += minitotal
            minitotal = zero(eltype(iter))
        end
        minitotal += x
    end
    return total + minitotal
end

```

I’m not particularly good with numerical methods, but wouldn’t that be similar in efficiency, with much better accuracy?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 17, 2018, 7:02pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/2 "2018-12-17T19:02:00Z")

</div>

> [@cstjean](#):
>
> Is that just Something I Must Watch Out For?

Yes, this is a big topic. For example, see [Kahan summation algorithm - Wikipedia](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) . The Julia standard library has alternative summation implementation rules for this.

---

<div class="post-metadata">

**Author:** ![cstjean](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cstjean/32/1444_2.png) [@cstjean](https://discourse.julialang.org/u/cstjean)\
**Post date:** [December 17, 2018, 7:16pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/3 "2018-12-17T19:16:51Z")

</div>

That’s the algorithm used by `sum(::Vector)`, right? Why can’t it be used with an iterable/generator?

---

<div class="post-metadata">

**Author:** ![simonbyrne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonbyrne/32/19_2.png) [@simonbyrne](https://discourse.julialang.org/u/simonbyrne)\
**Post date:** [December 17, 2018, 7:20pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/4 "2018-12-17T19:20:42Z")

</div>

That’s a great example of the problem of sequential summation.

For arrays, Julia uses a hybrid sequential/pairwise approach for summation (using sequential for small chunks, then accumulating these chunks in a pairwise manner), which is why `sum(aa)` avoids this problem.

However there’s no reason why we couldn’t do the same thing for general iterators: could you open an issue?

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [December 17, 2018, 8:17pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/5 "2018-12-17T20:17:27Z")

</div>

> [@simonbyrne](#):
>
> However there’s no reason why we couldn’t do the same thing for general iterators: could you open an issue?

Not sure about general iterables (the algorithm would have to be adapted since we can’t access items at arbitrary positions like for arrays), but for `skipmissing` we could definitely do that.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 17, 2018, 8:42pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/6 "2018-12-17T20:42:55Z")

</div>

> [@simonbyrne](#):
>
> However there’s no reason why we couldn’t do the same thing for general iterators: could you open an issue?

A direct implementation of pairwise summation requires random access, which is we we don’t use the pairwise algorithm for general iterators.

However, I suppose it might be possible to implement pairwise summation on a sequential iterator by allocating an explicit “stack” of O(log₂n) accumulators, and doing accumulation on the correct one as we go along.

---

<div class="post-metadata">

**Author:** ![simonbyrne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonbyrne/32/19_2.png) [@simonbyrne](https://discourse.julialang.org/u/simonbyrne)\
**Post date:** [December 17, 2018, 9:47pm UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/7 "2018-12-17T21:47:19Z")

</div>

The current implementation does access them in order, it just uses the call stack as the storage stack. However you do need to know the length ahead of time (so you know how deep the stack needs to be).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [December 18, 2018, 7:48am UTC](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759/8 "2018-12-18T07:48:21Z")

</div>

> [@simonbyrne](#):
>
> you do need to know the length ahead of time (so you know how deep the stack needs to be)

Not necessarily; you can just expand as necessary. Eg after 10 levels of depth are used up, add another 10. Someone just needs to code this, which can be tricky to do in a compiler-friendly way if the iterator has unknown eltype.