soft/garrote thresholding returns NaN for zero-valued data at threshold 0

[pavelkomarov]

Summary

pywt.threshold(data, 0, mode='soft') returns NaN for any element of data that is exactly zero, instead of 0. Soft thresholding is mathematically well defined there — sign(x)·max(|x| − λ, 0) is 0 at x = 0 for every λ ≥ 0 — so this is purely an artifact of the multiplicative form used internally. mode='garrote' and threshold_firm have the same defect.

Versions: PyWavelets 1.8.0, NumPy 2.5.2, Python 3.13.

Reproducing

>>> import numpy as np, pywt
>>> pywt.threshold(np.array([0.0, 1.0, -2.0]), 0.0, mode='soft')
array([nan,  1., -2.])          # expected: array([ 0.,  1., -2.])

with a RuntimeWarning: invalid value encountered in divide from _thresholding.py:22.

Affected, all at value=0 on an input containing an exact zero:

mode / function	result	expected
soft	[nan, 1., -2.]	[0., 1., -2.]
garrote	[nan, 1., -2.]	[0., 1., -2.]
threshold_firm(d, 0.0, 0.0)	[nan, 1., -2.]	[0., 1., -2.]
soft, complex	[nan+nanj, 1.+1.j]	[0.+0.j, 1.+1.j]

hard, greater, and less are correct.

Cause

In soft (pywt/_thresholding.py):

with np.errstate(divide='ignore'):
    # divide by zero okay as np.inf values get clipped, so ignore warning.
    thresholded = (1 - value/magnitude)
    thresholded.clip(min=0, max=None, out=thresholded)
    thresholded = data * thresholded

The comment correctly anticipates data == 0 with value > 0: value/0 is inf, 1 - inf is -inf, that clips to 0, and 0 * 0 == 0. That path works and is silent.

What is missed is value == 0 and data == 0. Then value/magnitude is 0/0, which is not a divide-by-zero but an invalid operation yielding NaN; 1 - NaN is NaN, clip passes NaN through unchanged, and 0 * NaN is NaN. Because the errstate suppresses only divide and not invalid, the warning escapes to the user as well.

value=0 is a legitimate no-op request meaning "do not denoise", so callers reasonably pass it — and a single NaN detail coefficient then propagates through waverec to poison an entire reconstruction, which is how we hit this downstream in PyNumDiff.

Suggested fix

Zeroing where the magnitude is zero keeps the existing complex-data support (sign alone would not) and costs one masked write:

with np.errstate(divide='ignore', invalid='ignore'):
    thresholded = (1 - value/magnitude)
    thresholded.clip(min=0, max=None, out=thresholded)
    thresholded = data * thresholded
thresholded = np.where(magnitude > 0, thresholded, 0)

The same treatment applies to nn_garrote and to the value_high * (1 - value_low/magnitude) / vdiff expression in threshold_firm (_thresholding.py:241).

Happy to open a PR if you would like one.

---

Investigated and filed by Claude Opus 5 (Claude Code) on behalf of @pavelkomarov, PyNumDiff maintainer.