Document our restrictions on the `init` kwarg

[christiangnrd]

Also reported at JuliaGPU/AcceleratedKernels.jl#28.

Examples:

julia> using OhMyThreads

julia> scheduler = StaticScheduler(; chunksize = 2, split = :batch)
StaticScheduler
├ Chunking: fixed size (2), split :consecutive
└ Threadpool: default

julia> sumvec = [1, 1, 1, 1, 1];

julia> treduce(+, sumvec; scheduler, init=0.0) # Correct
5.0

julia> reduce(+, sumvec; init=0.0) # Correct
5.0

julia> treduce(+, sumvec; scheduler, init=1.0) # Incorrect
8.0

julia> reduce(+, sumvec; init=1.0)
6.0

julia> vcatvec = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0], [7.0, 8.0]];

julia> treduce(vcat, vcatvec; scheduler, init=Float64[]) # Correct
8-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0
 5.0
 6.0
 7.0
 8.0

julia> reduce(vcat, vcatvec; init=Float64[]) # Correct
8-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0
 5.0
 6.0
 7.0
 8.0

julia> treduce(vcat, vcatvec; scheduler, init=Float64[0.0]) # Incorrect
10-element Vector{Float64}:
 0.0
 1.0
 2.0
 3.0
 4.0
 0.0 <- WRONG
 5.0
 6.0
 7.0
 8.0

julia> reduce(vcat, vcatvec; init=Float64[0.0]) # Correct
9-element Vector{Float64}:
 0.0
 1.0
 2.0
 3.0
 4.0
 5.0
 6.0
 7.0
 8.0

[MasonProtter]

This is an incorrect usage of init. The init kwarg in OhMyThreads must be a neutral element for whatever the reducing operator is you've provided.

I thought that I had spelled this out more clearly in the docstrings, but looking at them, I can see why someone would miss this. Currently the docstring just says

• init: initial value of the reduction. Will be forwarded to mapreduce for the task-local sequential parts of the calculation.

What this is saying is that when you do treduce(op, [1,2,3,4]; init=1) we do something like

task1 = @spawn reduce(op, [1, 2]; init=1)
task2 = @spawn reduce(op, [3, 4]; init=1)
result = mapreduce(fetch, op, [task1, task2]; init=1)

This means that for a reduction to be valid, init must be such that op(init, x) ≈ x up to whatever tolerance you care about.

Thanks for bringing this to my attention, it needs to be documented better

[christiangnrd]

If a PR analogous to JuliaGPU/AcceleratedKernels.jl#20 for this package were opened, would it be merged or is it too breaking?

[MasonProtter]

Probably not, no. We follow the same semantics that Base.reduce follows:

help?> reduce
search: reduce treduce mapreduce tmapreduce treducemap replace

  reduce(op, itr; [init])

  Reduce the given collection itr with the given binary
  operator op. If provided, the initial value init must be a
  neutral element for op that will be returned for empty
  collections. It is unspecified whether init is used for
  non-empty collections.

So even if OhMyThreads didn't insert a init value in the middle of your reduction here, Base.reduce reserves the right to do so still, and we're calling (map)reduce under the hood.

It appears that AcceleratedKernels is just incorrectly using init and providing false promises for what it does.

[christiangnrd]

Thanks! I hadn't thought to check the Base docstrings to see what behaviour is defined and what isn't.

The second Base docstring for reduce is quite ambiguous and could probably be more clear, but I can see how it's trying to convey what the behaviour is here.

[christiangnrd]

Seems like this conversation has been had before JuliaLang/julia#49042

[anicusan]

We follow the same semantics that Base.reduce follows:

Not really - while the docs are ambiguous (and you find arguments either way, e.g. Jeff Bezanson in #49042: "pretend it's added to the collection once, so that's probably what the spec should be."), using it you see:

julia> reduce(+, [1,2,3], init=10)
16

As it stands, treduce cannot be used as a drop-in replacement of Base.reduce.

Some other arguments for it:

• ?accumulate says "An initial value init may optionally be provided by a keyword argument" - it's odd that init would mean initial value in accumulate, but op neutral in reduce.
• Across the JuliaGPU ecosystem there are tests to check for the same results between the Base and GPU versions, including that init is not used as op neutral (example and AcceleratedKernels issue).
• To make it work with independent parallel processing - including for accumulate which is apparently more clearly specified - we needed to split init and neutral anyways (GPUArrays mapreduce example).

[MasonProtter]

Not really - while the docs are ambiguous (and you find arguments either way, e.g. Jeff Bezanson in JuliaLang/julia#49042: "pretend it's added to the collection once, so that's probably what the spec should be."), using it you see:

You're conflating the implementation of reduce with the API of reduce. Currently, most of the methods implementing reduce typically only uses the init argument once, but semantically the function allows itself to use init as many times as it likes. Changes are currently being made (slowly) to the mapreduce functionality to change how init is handled, and to document how we will use it more clearly.

Just because (map)reduce currently does one thing does not mean that can be relied upon into the future.

See JuliaLang/julia#53945 which is more recent than Jeff's comment and is the mostly agreed upon direction forward decided by the triage team.

This is actually vital for the pairwise reductions being well defined (which is very important for the parallel reductions we do in OhMyThreads). The fact that mapreduce doesn't always re-use init actually does cause some problems with mapreduce on certain collections and operations, and this is something that JuliaLang/julia#55318 is working to address by using the init multiple times.

---

?accumulate says "An initial value init may optionally be provided by a keyword argument" - it's odd that init would mean initial value in accumulate, but op neutral in reduce.

That's not odd to me.

• First, my understanding is that init means the value used to initialize a (sub-)reduction. That's fully consistent with what they say
• Second, accumulate is a sequential operation in Base. i.e. foldl-based, not reduce-based. If you don't do any re-associations, you never have this problem. This is why foldl also doesn't say anything about init needing to be neutral either. (yes, I'm aware it's possible to parallelize accumulate operations, but this is not something that Base does. It has a separate internal function for a accumulate_pairwise).

---

Across the JuliaGPU ecosystem there are tests to check for the same results between the Base and GPU versions, including that init is not used as op neutral (example and JuliaGPU/AcceleratedKernels.jl#16).

The GPU ecosystem has more freedom to give different guarantees here if desired. They control the reduction all the way from top to bottom and completely reimplement it how they see fit. OhMyThreads works by splitting up data into chunks and then passing those chunks to Base.mapreduce. This means we are beholden to respecting the interface that Base tells us to respect.

---

To make it work with independent parallel processing - including for accumulate which is apparently more clearly specified - we needed to split `init and neutral anyways (GPUArrays mapreduce example).

I haven't tried to implement a parallelized accumulate yet for OhMyThreads, so I haven't really looked into this in too much detail. If it's truly desprately needed, we could perhaps have a kwarg like leftmost_value or firstvalue or something that does what you call init.