Got it. Yes, I’m very interested in this, too. How natural op(xs...) == reduce(op, xs) really is? If we implement optimization for op(...f.(xs)) (or another syntax), are people going to use it? Is it less “scary” and more intuitive than reduce?
Isn’t it enough to lower op(...f.(xs)) to reduce(op, broadcasted(f, xs)) so that broadcasting and reduction are fused? Can we do more than this?
Yes. Discussion is here: stop being clever with reduction result types · Issue #20560 · JuliaLang/julia · GitHub
Right: the reduction syntax acts as a barrier for broadcast and removes the call to materialize. I’m not sure anything else is required. We already discussed the dims keyword on the issue and it seemed clear that was an orthogonal problem.
For me:
julia> a = fill(rand(1000), 1000)
julia> @time reduce(vcat, a');
1.385112 seconds (4.99 k allocations: 3.729 GiB, 16.95% gc time)
julia> @time vcat(a'...);
0.150527 seconds (361.89 k allocations: 31.163 MiB, 1.50% gc time, 69.42% compilation tim
e)
For hcat, simply don’t use the ellipsis
Instead of reduce(vcat, a') it seems faster to do reduce(hcat,a)':
using BenchmarkTools
a = fill(rand(1000), 1000)
@btime vcat($a'...) # 3.374 ms (5017 allocations: 7.80 MiB)
@btime reduce(hcat,$a)' # 1.058 ms (2 allocations: 7.63 MiB)
@btime reduce(vcat, $a') # 1.298 s (4991 allocations: 3.73 GiB)
for the concatenation along the other direction, the following two seem to be equivalent:
@btime vcat($a...) # 1.046 ms (3 allocations: 7.64 MiB)
@btime reduce(vcat, $a) # 1.051 ms (3 allocations: 7.64 MiB)
Julia defines the identity function if you would like to avoid the x->x.
I’m just curious if my understanding is correct: Do these discrepancies arise because of Julia’s column-major ordering for N-D arrays?
No, it’s because there is an optimised version of reduce(hcat, a), but reduce(vcat, a') goes to the fallback routine since a' is not an AbstractVector.