Stable log-space ops in SciPy: logsumexp, log_expit, logaddexp, and the all--inf softmax NaN
Measured on NumPy 2.5.4 / SciPy 1.18.1. Naive log(sum(exp(x))) for x=[1000,1001,1002] gives inf and softmax gives [nan nan nan]; scipy.special.logsumexp gives 1002.4076 and scipy.special.softmax [0.090 0.245 0.665]. log(1/(1+exp(-z))) is -inf at z=-800 and exactly 0 at z=40; scipy.special.log_expit gives -800 and -4.25e-18. Binary cross-entropy for a confident wrong logit: -log(1-expit(40)) = inf, while np.logaddexp(0, 40) = 40.0 (use logaddexp(0, z) - y*z). Log of a difference: logsumexp(a, b=[1,-1], return_sign=True) gives log(e-1) correctly. Gotcha: a fully masked row (all -inf) still gives NaN from scipy softmax, with only a RuntimeWarning 'invalid value encountered in subtract'; logsumexp of it is -inf, which is fine. Handle all-masked rows explicitly.
Every number reproduced on NumPy 2.5.4 / SciPy 1.18.1 (CPython 3.12.16): inf and [nan nan nan] naive vs 1002.4076 and [0.090 0.245 0.665]; -inf/0.0 naive vs log_expit -800.0 and -4.248e-18; inf vs logaddexp(0, 40) = 40.0. logsumexp([1, 0], b=[1, -1], return_sign=True) gives (0.5413, 1.0) = log(e-1). The all -inf row gives softmax [nan nan] with 'invalid value encountered in subtract', and logsumexp gives -inf.