Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ENH: improve scipy.special.log_softmax accuracy in edge cases by a factor of 2**126 to 2**1022 #19549

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
+63 −5
Open

ENH: improve scipy.special.log_softmax accuracy in edge cases by a factor of 2**126 to 2**1022 #19549

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e573855
ENH: improve `scipy.special.log_softmax` accuracy
JasonGross Nov 17, 2023
1c0dca0
ENH: `log_softmax` supports tuple `axis` again
JasonGross Jan 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
60 changes: 55 additions & 5 deletions scipy/special/_logsumexp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,24 +284,74 @@ def log_softmax(x, axis=None):
>>> y
array([ 0., -inf])

>>> subnormal32 = np.finfo(np.float32).smallest_subnormal
>>> x = np.array([0, np.log(subnormal32)], dtype=np.float32)
>>> y = log_softmax(x)
>>> y
array([-1.40130e-45, -1.03279e+02], dtype=float32)

>>> with np.errstate(divide='ignore'):
... y = np.log(softmax(x))
...
>>> y
array([ 0. , -103.27893], dtype=float32)

>>> subnormal64 = np.finfo(np.float64).smallest_subnormal
>>> x = np.array([0, np.log(subnormal64)], dtype=np.float64)
>>> y = log_softmax(x)
>>> y
array([-4.9407e-324, -7.4444e+002])

>>> with np.errstate(divide='ignore'):
... y = np.log(softmax(x))
...
>>> y
array([ 0. , -744.44007])
"""

x = _asarray_validated(x, check_finite=False)

x_max = np.amax(x, axis=axis, keepdims=True)
# work around https://github.com/numpy/numpy/issues/25623
if isinstance(axis, tuple):
x_max = np.amax(x, axis=axis, keepdims=True)
else:
x_argmax = np.argmax(x, axis=axis, keepdims=True)
# work around https://github.com/numpy/numpy/issues/25622
if axis is None:
x_argmax = x_argmax.flatten()
x_max = np.take_along_axis(x, x_argmax, axis=axis)

finite_max_mask = np.isfinite(x_max)

if x_max.ndim > 0:
x_max[~np.isfinite(x_max)] = 0
elif not np.isfinite(x_max):
x_max = 0
x_max[~finite_max_mask] = 0
elif not finite_max_mask:
x_max = np.zeros_like(x_max)

tmp = x - x_max
exp_tmp = np.exp(tmp)

# work around https://github.com/numpy/numpy/issues/25623
if isinstance(axis, tuple):
pass
else:
# we know that exp_tmp at the location of the max is either 1 or infinite,
# depending on finite_max_mask, so we can set it to zero and use log1p
if exp_tmp.ndim > 0:
exp_tmp_max = np.take_along_axis(exp_tmp, x_argmax, axis=axis)
exp_tmp_max[finite_max_mask] = 0
np.put_along_axis(exp_tmp, x_argmax, exp_tmp_max, axis=axis)
elif finite_max_mask:
exp_tmp = np.zeros_like(exp_tmp)

# suppress warnings about log of zero
with np.errstate(divide='ignore'):
s = np.sum(exp_tmp, axis=axis, keepdims=True)
out = np.log(s)
# work around https://github.com/numpy/numpy/issues/25623
if isinstance(axis, tuple):
out = np.log(s)
else:
out = np.log1p(s)

out = tmp - out
return out
8 changes: 8 additions & 0 deletions scipy/special/tests/test_log_softmax.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,14 @@

@pytest.mark.parametrize('x, expected', [
(np.array([1000, 1]), np.array([0, -999])),
# we shouldn't return zero on the smallest subnormal input
(np.array([-np.log(np.finfo(np.float32).smallest_subnormal), 0], dtype=np.float32),
np.array([float.fromhex('-0x1.00000p-149'), float.fromhex('-0x1.9d1dap+6')],
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This test (I believe) is the one test that fails on Windows with

Not equal to tolerance rtol=1e-13, atol=0

Mismatched elements: 1 / 2 (50%)
Max absolute difference among violations: 5.e-324
Max relative difference among violations: 1.
 ACTUAL: array([   0.      , -744.440072])
 DESIRED: array([-4.940656e-324, -7.444401e+002])

This seems like a bug in the windows implementation of floating point arithmetic. What should we do? Is there a way to relax this test only on Windows?

dtype=np.float32)),
(np.array([-np.log(np.finfo(np.float64).smallest_subnormal), 0], dtype=np.float64),
np.array([float.fromhex('-0x0.0000000000001p-1022'),
float.fromhex('-0x1.74385446d71c3p+9')],
dtype=np.float64)),
# Expected value computed using mpmath (with mpmath.mp.dps = 200) and then
# converted to float.
Expand Down