Fixes initial value for reduce in logsumexp

[itsdfish]

This PR request fixes the initial value for reduce in logsumexp. This error was discovered in the following issue in Turing (see final comment).

[mohamed82008]

This is not right, X can be all negative.

[itsdfish]

Thanks for catching the error. Do you know what the correct initial value should be?

[mohamed82008]

I think we should not have init at all. Just check if X is empty and return log(zero(T)) then. Otherwise, keep going with the reduce.

[itsdfish]

Good call. I made the changes you suggested and the results hold for both negative and positive values. Please review if you can.

[mohamed82008]

Would be nice to add a ForwardDiff test to catch such error if it happens again in the future.

[nalimilan]

This isn't correct when dims != :. Can you describe the problem you're trying to fix (not knowing anything about what ForwardDiff is trying to do in that issue)? Anyway we'll need to add a test for it.

[mohamed82008]

The problem is that log(ForwardDiff.Dual{Nothing}(0.0, 0.0)) is ForwardDiff.Dual{Nothing}(-Inf, NaN). The NaN then propagates to u when doing the reduce.

[mohamed82008]

ForwardDiff.gradient(logsumexp, rand(3)) now gives a vector of NaN.

[nalimilan]

I guess we could use another way to get -Inf of the right type, e.g. convert(T, -Inf). Would that work? Or is there another solution which is OK for ForwardDiff?

[mohamed82008]

Yup init=convert(T, -Inf) works like a charm!

[mohamed82008]

I am curious though, why is the solution here not correct when dims != :? Could you please give an example?

[nalimilan]

Yup init=convert(T, -Inf) works like a charm!

Cool. Actually this won't work for integer inputs, so better use something like oftype(log(zero(T)), -Inf). Would you update the PR to use that? It would also be nice to have a test or at least a comment explaining why we don't use the simplest syntax.

I am curious though, why is the solution here not correct when dims != :? Could you please give an example?

julia> reduce(max, zeros(0, 2), dims=1, init=log(zero(Int)))
1×2 Array{Float64,2}:
 -Inf  -Inf

julia> reduce(max, zeros(0, 2), dims=1)
ERROR: ArgumentError: reducing over an empty collection is not allowed

[mohamed82008]

@nalimilan the reduce won't be called if isempty(X) is true which it is in case of zeros(0, 2). Oh but the shape will still be wrong and it is not type stable.

[mohamed82008]

I don't have write access to this PR, so @itsdfish please use:

u = reduce(max, X, dims=dims, init=oftype(log(zero(T)), -Inf))

instead of

u = reduce(max, X, dims=dims, init=log(zero(T)))

[itsdfish]

The updates seem to work on my end.

[nalimilan]

Thanks. One last check: is it really expected/correct that log(ForwardDiff.Dual{Nothing}(0.0, 0.0)) gives ForwardDiff.Dual{Nothing}(-Inf, NaN)? Where does that NaN come from? Maybe we should discuss this with ForwardDiff developers?

Anyway, if we make this change, we should add a comment to ensure somebody doesn't revert this without knowing the reason why it exists.

EDIT: I've found http://www.juliadiff.org/ForwardDiff.jl/stable/user/advanced/#Fixing-NaN/Inf-Issues-1, so this seems to be expected.

[itsdfish]

To make things easier, I added you both as collaborators in case you want to make changes.

[mohamed82008]

Perhaps merge this?

[nalimilan]

#83