LogSumExp of -inf returns nan instead of -inf

[ricardoV94]

x = at.vector('x')
res1 = at.log(at.sum(at.exp(x)))  # nan with sum
res2 = at.log(at.prod(at.exp(x)))  # but not with prod

fun = aesara.function([x], [res1, res2])
print(fun(np.array([-np.inf, -np.inf])))  # [array(nan), array(-inf)]

for reference in numpy we get -inf, together with a divide by zero encountered in log

print(np.log(np.sum(np.exp([-np.inf, -np.inf]))))  # -inf

This showed up in the pymc3 logsumexp function which returns nan with this input:

from pymc3.math import logsumexp

x = at.vector('x')
res = logsumexp(x)
fun = aesara.function([x], res) 

print(fun(np.array([-np.inf, -np.inf])))  # [nan]

Whereas the scipy reference works fine

print(scipy.special.logsumexp([-np.inf, -np.inf]))  # -inf

Weirdly it happens with addition, but not subtraction or multiplication, so maybe the problem is not the sum but the addition?

x = at.vector('x')
res1 = at.log(at.exp(x[0]) + at.exp(x[1]))
res2 = at.log(at.exp(x[0]) - at.exp(x[1]))
fun = aesara.function([x], [res1, res2])
print(fun(np.array([-np.inf, -np.inf])))  # [array(nan), array(-inf)]

[ricardoV94]

I think it has to do with an automatic optimization which is basically the same as what the pymc3 logsumexp attempts to do, but with this line missing:

https://github.com/pymc-devs/pymc3/blob/d7172c0a1a76301031d1b3b411d00643c416a0c4/pymc3/math.py#L192

x = at.vector('x')
res = at.log(at.exp(x[0]) + at.exp(x[1]))
fun = aesara.function([x], res)
aesara.dprint(fun)

Elemwise{Composite{(maximum(i0, i1) + scalar_softplus(((i0 - maximum(i0, i1)) + (i1 - maximum(i0, i1)))))}} [id A] ''   2
 |Subtensor{int64} [id B] ''   1
 | |x [id C]
 | |ScalarConstant{0} [id D]
 |Subtensor{int64} [id E] ''   0
   |x [id C]
   |ScalarConstant{1} [id F]

The softplus input in the graph above has a -inf - (-inf), yielding the nan in the output