ENH: signal: add array API support / delegation to lfilter et al #21713

ev-br · 2024-10-15T18:43:13Z

~~Is on top of #20772~~

Add array API / delegation support to more routines from signaltools: lfilter, sosfilt, filtfilt, sosfiltfilt. Also lfiltic and lfilt_zi. While at it, also add deconvolve since it's only a small thing on top of lfilter anyway.

tylerjereddy

I noticed 213 new failures on GPU testing vs. main via SCIPY_DEVICE=cuda python dev.py test -j 8 -b all.

Many can be circumvented by skipping GPU tests with torch as usual, but I commented on a sample of one of the more interesting cases.

tylerjereddy · 2024-10-15T23:52:47Z

scipy/signal/tests/test_signaltools.py

        zi = sosfilt_zi(sos)

-        y, zf = sosfilt(sos, np.ones(40, dt), zi=zi)
-        assert_allclose_cast(zf, zi, rtol=1e-13)
+        y, zf = sosfilt(sos, xp.ones(40, dtype=dt), zi=zi)


Something interesting is happening with the test failure here via i.e., python dev.py test -t scipy/signal/tests/test_signaltools.py::TestSOSFilt::test_sosfilt_zi -b cupy

I noticed that zf in CuPy is the transpose of the result provided by NumPy and torch; i.e., this causes the above incantation to pass:

--- a/scipy/signal/tests/test_signaltools.py +++ b/scipy/signal/tests/test_signaltools.py @@ -3954,6 +3954,7 @@ class TestSOSFilt: with pytest.raises(ValueError, match='Invalid zi shape'): sosfilt(sos, x, zi=zi, axis=1) + @skip_xp_backends(cpu_only=True, exceptions=['cupy']) def test_sosfilt_zi(self, dt, xp): dt = getattr(xp, dt) sos = signal.butter(6, 0.2, output='sos') @@ -3961,7 +3962,7 @@ class TestSOSFilt: zi = sosfilt_zi(sos) y, zf = sosfilt(sos, xp.ones(40, dtype=dt), zi=zi) - xp_assert_close(zf, zi, rtol=1e-13, check_dtype=False) + xp_assert_close(zf.T, zi, rtol=1e-13, check_dtype=False) # Expected steady state value of the step response of this filter: ss = xp.prod(xp.sum(sos[:, :3], axis=-1) / xp.sum(sos[:, 3:], axis=-1))

We obviously can't/shouldn't use a patch like that, filtered to CuPy case. At first glance, I can't imagine this would be an intentional change on the CuPy side based on the docs. I do see that the zi shape convention is a bit different between SciPy and CuPy, but upon closer inspection I wasn't entirely convinced that was the cause since sosfilt_zi should compensate by preserving the larger zi convention in CuPy anyway.

ev-br · 2024-10-16T17:03:13Z

Hm. Apparently, cupyx.scipy.signal.sosfilt_zi is not a drop-in replacement for scipy.signal.sosfilt_zi:

sosfilt_zi returns zi of different shapes in scipy and cupyx: in SciPy it's (n_sections, 2) and in cupyx.scipy it's (n_sections, 4). Both are documented but the docs do not explain what these 2 and 4 are.
given identical sos inputs, the outputs of scipy.signal.sosfilt_zi and cupyx.scipy.signal.sosfilt_zi are numerically very different:

scipy:

(Pdb) p zi
array([[ 0.00525033, -0.00120091],
       [ 0.07729171, -0.02862319],
       [ 0.91711742, -0.65303262]])
(Pdb) a    
sos = array([[ 3.40537653e-04,  6.81075305e-04,  3.40537653e-04,
         1.00000000e+00, -1.03206941e+00,  2.75707942e-01],
       [ 1.00000000e+00,  2.00000000e+00,  1.00000000e+00,
         1.00000000e+00, -1.14298050e+00,  4.12801598e-01],
       [ 1.00000000e+00,  2.00000000e+00,  1.00000000e+00,
         1.00000000e+00, -1.40438489e+00,  7.35915191e-01]])

and cupyx.scipy:

(Pdb) p zi
array([[1.        , 1.        , 0.00559087, 0.00559087],
       [0.00559087, 0.00559087, 0.08288258, 0.08288258],
       [0.08288258, 0.08288258, 1.        , 1.        ]])
(Pdb) a
sos = array([[ 3.40537653e-04,  6.81075305e-04,  3.40537653e-04,
         1.00000000e+00, -1.03206941e+00,  2.75707942e-01],
       [ 1.00000000e+00,  2.00000000e+00,  1.00000000e+00,
         1.00000000e+00, -1.14298050e+00,  4.12801598e-01],
       [ 1.00000000e+00,  2.00000000e+00,  1.00000000e+00,
         1.00000000e+00, -1.40438489e+00,  7.35915191e-01]])

Before opening a cupy issue, let me ping @andfoy : as the author of cupy's sosfilt and sosfilt_zi, could you weigh in on the differences between scipy's and cupy's implementations --- is this something to fix on the cupy side maybe?

andfoy · 2024-10-16T17:41:39Z

The difference is that the state is not being represented as a sum of delays:

Given IIR coefficients $[1, a, b]$ and previous states $[k_{-2}, k_{-1}]$, SciPy states are represented as:

z_{-1} = k_{-1} * a + k_{-2} * b 
z_0 =  k_{-1} * b

Such that when computing the output y[0] and y[1] over input x it yields the following values

y[0] = x[0] + z_{-1}
y[1] = x[1] + y[0] * a  + z_{0}

This kind of sum carrying is possible given that SciPy filter computation occurs serially.

When computing the filter output in GPU, the previous state values are required as they are, that is [k_{-2}, k_{-1}] since such representation simplifies the computation process outlined in https://github.com/cupy/cupy/blob/fb6ffd71d3da200d4300c5760327a278cd8d9ee3/cupyx/scipy/signal/_signaltools.py#L696, while CuPy lfilter/sosfilt previous states can be converted into SciPy format using lfiltic, the reverse would imply solving a linear equation of the form:

$$\begin{bmatrix} c_0 & c_1 & \cdots & c_{n -1} & c_{n-2} \\\ c_1 & c_2 & \cdots & c_{n - 2} & 0 \\\ \vdots & \vdots & \ddots & \vdots & \vdots \\\ c_{n - 2} & 0 & \cdots & 0 & 0 \end{bmatrix} \begin{bmatrix} k_{-1} \\\ k_{-2} \\\ \vdots \\\ k_{-N} \end{bmatrix} = \begin{bmatrix} z_{-N} \\\ z_{-N + 1} \\\ \vdots \\\ z_{0} \end{bmatrix}$$

ev-br · 2024-10-16T18:20:27Z

Thanks Edgar! This should probably be great to document, for both scipy and cupy.

That said, this does not explain the need to transpose noted in #21713 (comment) though?

andfoy · 2024-10-16T19:18:57Z

Now that you mention it, maybe this line comes to mind w.r.t output shape manipulation? https://github.com/cupy/cupy/blob/fb6ffd71d3da200d4300c5760327a278cd8d9ee3/cupyx/scipy/signal/_iir_utils.py#L730

andfoy · 2024-10-17T17:15:17Z

I've just opened a PR in CuPy that fixes the issue at hand: cupy/cupy#8677

ev-br · 2024-10-19T15:20:51Z

Added a commit with a slew of skips for CuPy, and a separate commit to avoid delegating for functions which deliberately deviate in scipy and cupy.
I believe skipping tests is the right thing to do today; when Edgar's CuPy PR lands, we will be able to undo most of the skips, hopefully.

Note: Some CuPy functions deliberately deviate from scipy: {lfilter,sosfilt}_zi return arrays of different shapes. Skip dispatching for these cases; a user will have to be explicit on whether they want the 'original' SciPy output (then call lfilter_zi with NumPy array arguments) or the CuPy output (then call cupyx.scipy.signal.lfilter_zi with CuPy array arguments).

scipy/signal/_signaltools.py

lucascolley · 2024-12-09T21:34:16Z

scipy/signal/_signaltools.py

+        b = xp.concat((b, xp.zeros(n - b.shape[0], dtype=b.dtype)))
+
+    dt = xp.result_type(a, b)
+    IminusA = np.eye(n - 1) - linalg.companion(a).T


does this work, passing a to linalg.companion for non-NumPy arrays?

I suppose this goes through __array__. As long as it works, I'd keep it as is, until at some point scipy.linalg.companion gets array api compatible, and we'll fix this in one go.

scipy/signal/_signaltools.py

The replacement for `np.atleast_1d(x)` is `xpx.atleast_nd(xp.asarray(x), ndim=1, xp=xp)`, no less.

ev-br · 2024-12-10T09:48:06Z

Updated to address review comments, add some array-like tests, improve a link to a CuPy fix to investigate after a CuPy makes a 14.0 release, and added enough skips for tests to pass locally with SCIPY_DEVICE=cuda.

lucascolley

looks pretty good, thanks!

scipy/signal/tests/test_signaltools.py

lucascolley

thanks Evgeni!

ev-br requested review from larsoner, andyfaff, person142, steppi and ilayn as code owners October 15, 2024 18:43

ev-br added the array types Items related to array API support and input array validation (see gh-18286) label Oct 15, 2024

ev-br removed enhancement A new feature or improvement scipy.special scipy._lib CI Items related to the CI tools such as CircleCI, GitHub Actions or Azure Meson Items related to the introduction of Meson as the new build system for SciPy labels Oct 15, 2024

lucascolley added the enhancement A new feature or improvement label Oct 15, 2024

tylerjereddy reviewed Oct 15, 2024

View reviewed changes

andfoy mentioned this pull request Oct 17, 2024

Fix sosfilt state output shape when ndim < 2 cupy/cupy#8677

Merged

ev-br force-pushed the sigtools_convolve_cupy_lfilt branch from 46eac28 to 58f4b2b Compare October 19, 2024 15:07

ev-br mentioned this pull request Oct 19, 2024

MAINT: signal deprecation cleanups #21724

Closed

ev-br mentioned this pull request Oct 30, 2024

ENH: signal.windows: add array API support (take 2) #21783

Merged

ev-br mentioned this pull request Dec 1, 2024

BUG: signal: actually reject objects in correlate/convolve #21984

Merged

ev-br force-pushed the sigtools_convolve_cupy_lfilt branch from e9ad517 to 5a4876e Compare December 9, 2024 20:28

ev-br mentioned this pull request Dec 9, 2024

ENH: array types, signal: delegate to CuPy and JAX for correlations and convolutions #20772

Merged

lucascolley changed the title ~~ENH: array types, signal: add array API support / delegation to lfilter et al~~ ENH: signal: add array API support / delegation to lfilter et al Dec 9, 2024

lucascolley reviewed Dec 9, 2024

View reviewed changes

BUG: signal: restore (and test for) array-like list inputs

bb0cf02

The replacement for `np.atleast_1d(x)` is `xpx.atleast_nd(xp.asarray(x), ndim=1, xp=xp)`, no less.

ev-br force-pushed the sigtools_convolve_cupy_lfilt branch from 9b61c03 to bb0cf02 Compare December 10, 2024 09:23

MAINT: signal: link to a cupy fix

35cf20d

lucascolley reviewed Dec 10, 2024

View reviewed changes

scipy/signal/tests/test_signaltools.py Show resolved Hide resolved

scipy/signal/tests/test_signaltools.py Outdated Show resolved Hide resolved

scipy/signal/tests/test_signaltools.py Outdated Show resolved Hide resolved

scipy/signal/tests/test_signaltools.py Show resolved Hide resolved

lucascolley added this to the 1.16.0 milestone Dec 10, 2024

ev-br force-pushed the sigtools_convolve_cupy_lfilt branch from 9351b41 to b062440 Compare December 11, 2024 06:33

address review comments

40a626e

ev-br force-pushed the sigtools_convolve_cupy_lfilt branch from b062440 to 40a626e Compare December 11, 2024 06:51

lucascolley approved these changes Dec 11, 2024

View reviewed changes

lucascolley merged commit e29b9a4 into scipy:main Dec 11, 2024
37 checks passed

Uh oh!

ENH: signal: add array API support / delegation to lfilter et al #21713

ENH: signal: add array API support / delegation to lfilter et al #21713

Uh oh!

Conversation

ev-br commented Oct 15, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

tylerjereddy left a comment

Choose a reason for hiding this comment

Uh oh!

tylerjereddy Oct 15, 2024

Choose a reason for hiding this comment

Uh oh!

ev-br commented Oct 16, 2024

Uh oh!

andfoy commented Oct 16, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ev-br commented Oct 16, 2024

Uh oh!

andfoy commented Oct 16, 2024

Uh oh!

andfoy commented Oct 17, 2024

Uh oh!

ev-br commented Oct 19, 2024

Uh oh!

Uh oh!

Uh oh!

lucascolley Dec 9, 2024

Choose a reason for hiding this comment

Uh oh!

ev-br Dec 10, 2024

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

ev-br commented Dec 10, 2024

Uh oh!

lucascolley left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

lucascolley left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

ev-br commented Oct 15, 2024 •

edited

Loading

andfoy commented Oct 16, 2024 •

edited

Loading