Missing `reduce(-, x,dims = 2)` implementation

[miguelraz]

Found via tweet:
https://twitter.com/code_report/status/1582766635873886209?s=20&t=8MmrCCKan8fryH3G-HZJNw

Couldn't find a related issue.

Maybe just a case of adding it to these cases?

julia/base/reduce.jl

Line 346 in 35191ad

reduce_empty(::typeof(add_sum), ::Type{Union{}}) = _empty_reduce_error(add_sum, Union{})

[ribeiro-juliano]

Same happens with /. If you don't feed in an init, it returns a similar error:

ERROR: MethodError: no method matching reducedim_init(::typeof(identity), ::typeof(/), ::Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, ::Int64)
Closest candidates are:
  reducedim_init(::Any, ::typeof(min), ::AbstractArray, ::Any) at reducedim.jl:128
  reducedim_init(::Any, ::typeof(max), ::AbstractArray, ::Any) at reducedim.jl:128
  reducedim_init(::Any, ::typeof(&), ::Union{Base.AbstractBroadcasted, AbstractArray}, ::Any) at reducedim.jl:206
  ...
Stacktrace:
 [1] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}, dims::Int64)
   @ Base ./reducedim.jl:371
 [2] #mapreduce#765
   @ ./reducedim.jl:357 [inlined]
 [3] #reduce#767
   @ ./reducedim.jl:406 [inlined]
 [4] top-level scope
   @ REPL[9]:1

But it works with init=1.0 and reduce(/, x) also shows expected behavior if init=1.0 was considered the default.

[Moelf]

right, - is inverse of + and the identity is the same, zero(T).

for * and / the identity is one(T)

julia> reduce(*, a; dims=2)
2×1 Matrix{Float64}:
 0.3107919222529072
 0.05811077660665645

julia> reduce(/, a; dims=2)
ERROR: MethodError: no method matching reducedim_init(::typeof(identity), ::typeof(/), ::Matrix{Float64}, ::Int64)
Closest candidates are:
  reducedim_init(::Any, ::typeof(min), ::AbstractArray, ::Any) at reducedim.jl:128
  reducedim_init(::Any, ::typeof(max), ::AbstractArray, ::Any) at reducedim.jl:128
  reducedim_init(::Any, ::typeof(&), ::Union{Base.AbstractBroadcasted, AbstractArray}, ::Any) at reducedim.jl:206

[ribeiro-juliano]

one(T) would still display the error in my ill-conceived code I jumbled together just to test this because the array I generated was of integers, meaning one(T)=1, which still causes an error:

julia>x = reshape(1:15, (5,3))
5×3 reshape(::UnitRange{Int64}, 5, 3) with eltype Int64:
 1   6  11
 2   7  12
 3   8  13
 4   9  14
 5  10  15
julia> reduce(/,x,dims=2,init=1)
ERROR: InexactError: Int64(0.5)
Stacktrace:
  [1] Int64
    @ ./float.jl:788 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:968 [inlined]
  [4] setindex!
    @ ./multidimensional.jl:674 [inlined]
  [5] macro expansion
    @ ./reducedim.jl:317 [inlined]
  [6] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [7] _mapreducedim!(f::typeof(identity), op::typeof(/), R::Matrix{Int64}, A::Matrix{Float64})
    @ Base ./reducedim.jl:316
  [8] mapreducedim!(f::Function, op::Function, R::Matrix{Int64}, A::Matrix{Float64})
    @ Base ./reducedim.jl:324
  [9] _mapreduce_dim
    @ ./reducedim.jl:368 [inlined]
 [10] #mapreduce#765
    @ ./reducedim.jl:357 [inlined]
 [11] #reduce#767
    @ ./reducedim.jl:406 [inlined]
 [12] top-level scope
    @ REPL[28]:1

This doesn't seem to be expected behavior since performing reduce(/,x) still works even with an array of integers. Either way the Float case causing problems I think would warrant the same treatment as - , though.

[mkitti]

Here is where the default is defined for + and * on current master.

julia/base/reducedim.jl

Lines 107 to 112 in 0d52506

	function reducedim_init(f, op::Union{typeof(+),typeof(add_sum)}, A::AbstractArray, region)
	_reducedim_init(f, op, zero, sum, A, region)
	end
	function reducedim_init(f, op::Union{typeof(*),typeof(mul_prod)}, A::AbstractArray, region)
	_reducedim_init(f, op, one, prod, A, region)
	end

Is there a potential issue between signed and unsigned values here?

[andrewjradcliffe]

I can give a partial explanation. Only a few binary operators have a sensible neural element. Recall also that a binary operator need not be associative; -, /,^ are some examples.

Now, we can simplify some cases, e.g. when reducing over an entire array (treating it as a vector), as we can make a reasonable claim as to what the first element should be, thereby alleviating the need for a neutral element. This is why expressions like reduce(op, A) will typically work -- the first element is assumed to be first(A)

However, in the multidimensional case, it is ambiguous as to what the first element should be given that our operator may not be associative -- if I say, for example, for A ∈ ℝᴺˣᴹˣᴾ, reduce(op, A, dims=(1,2)), do I intend that a first element be found on the first row of the first column of each slice? Perhaps. Given a binary operator without the associative property, however, this implies that I am okay with starting the computation with whatever value happens to the first element -- a seemingly random choice in the general case (it may make sense in a specific case where the values held in the array are intentionally ordered).
In essence, unless a binary operator is associative, there is no way to sensibly begin a computation, hence, the need for a user to supply an initial value.

See also the reduce docstring.
As for the neutrality of 0 w.r.t. -:

# with a neutral element
julia> a = [1, 2, 3];

julia> v0 = 0
0

julia> v1 = -(v0, a[1])
-1

julia> v2 = -(v1, a[2])
-3

julia> v3 = -(v2, a[3])
-6

# without a neutral element
julia> v1 = a[1]
1

julia> v2 = -(v1, a[2])
-1

julia> v3 = -(v2, a[3])
-4