Repository: cupy/cupy
Pull request: Fix sparse min/max to return sparse matrix
Close #3504.

This PR fixes sparse min/max to return a value as a sparse matrix. I will parameterize the test in another PR.

Problem statement:
Sparse min/max have less SciPy-compatibility
The outptus of sparse min/max are less compatible to those of Scipy.

```Python
import numpy
import scipy

import cupy
from cupy import testing
import cupyx.scipy

def _make(xp, sp):
    m = testing.shaped_random((10, 20), xp=xp, dtype='f', scale=1)
    m[m < 0.95] = 0
    return sp.sparse.csr_matrix(m)

for (xp, sp) in [(numpy, scipy), (cupy, cupyx.scipy)]:
    m = _make(xp, sp).min(axis=0)
    print(m.dtype, m.shape, type(m))
```

**Output**
```Python
float32 (1, 20) <class 'scipy.sparse.coo.coo_matrix'>
float64 (20,) <class 'cupy.core.core.ndarray'>
```

Hints from issue discussion:

Commit message: Return sparse matrix

Patch to review:
diff --git a/cupyx/scipy/sparse/data.py b/cupyx/scipy/sparse/data.py
index 891dc0af470..b33fe84443d 100644
--- a/cupyx/scipy/sparse/data.py
+++ b/cupyx/scipy/sparse/data.py
@@ -1,5 +1,6 @@
 import cupy
 from cupyx.scipy.sparse import base
+from cupyx.scipy.sparse import coo
 from cupyx.scipy.sparse import util
 
 
@@ -139,6 +140,7 @@ def _min_or_max_axis(self, axis, min_or_max, sum_duplicates, nonzero):
         N = self.shape[axis]
         if N == 0:
             raise ValueError("zero-size array to reduction operation")
+        M = self.shape[1 - axis]
 
         mat = self.tocsc() if axis == 0 else self.tocsr()
         if sum_duplicates:
@@ -146,8 +148,20 @@ def _min_or_max_axis(self, axis, min_or_max, sum_duplicates, nonzero):
 
         # Do the reudction
         value = mat._minor_reduce(min_or_max, axis, nonzero)
+        major_index = cupy.arange(M)
 
-        return value
+        mask = value != 0
+        major_index = cupy.compress(mask, major_index)
+        value = cupy.compress(mask, value)
+
+        if axis == 0:
+            return coo.coo_matrix(
+                (value, (cupy.zeros(len(value)), major_index)),
+                dtype=self.dtype, shape=(1, M))
+        else:
+            return coo.coo_matrix(
+                (value, (major_index, cupy.zeros(len(value)))),
+                dtype=self.dtype, shape=(M, 1))
 
     def _min_or_max(self, axis, out, min_or_max, sum_duplicates, non_zero):
         if out is not None:
diff --git a/tests/cupyx_tests/scipy_tests/sparse_tests/test_csc.py b/tests/cupyx_tests/scipy_tests/sparse_tests/test_csc.py
index 36b21e3edc7..cd9dfbea315 100644
--- a/tests/cupyx_tests/scipy_tests/sparse_tests/test_csc.py
+++ b/tests/cupyx_tests/scipy_tests/sparse_tests/test_csc.py
@@ -915,10 +915,9 @@ def test_min_sparse_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=0))
-        da_scipy_values = numpy.array(dm_data.min(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=0)
+        da_scipy_values = dm_data.min(axis=0)
+        testing.assert_array_equal(da_scipy_values.A, da_scipy_values.A)
 
     def test_min_dense_axis_0(self):
         dm_data = numpy.random.random((10, 20))
@@ -929,10 +928,9 @@ def test_min_dense_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=0))
-        da_scipy_values = numpy.array(dm_data.min(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=0)
+        da_scipy_values = dm_data.min(axis=0)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_min_axis_0_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -943,10 +941,10 @@ def test_min_axis_0_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=0, nonzero=True))
-        da_numpy_values = numpy.array([10, 1, 2, 3, 4,
-                                       5, 6, 7, 8, 9]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.min(axis=0, nonzero=True)
+        da_numpy_values = numpy.array([[10, 1, 2, 3, 4,
+                                        5, 6, 7, 8, 9]]).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_min_sparse_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -958,10 +956,9 @@ def test_min_sparse_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=1))
-        da_scipy_values = numpy.array(dm_data.min(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=1)
+        da_scipy_values = dm_data.min(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_min_dense_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -972,10 +969,9 @@ def test_min_dense_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=1))
-        da_scipy_values = numpy.array(dm_data.min(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=1)
+        da_scipy_values = dm_data.min(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_min_axis_1_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -986,10 +982,11 @@ def test_min_axis_1_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=1, nonzero=True))
-        da_numpy_values = numpy.array([1, 10, 20, 30, 40,
-                                       50, 60, 70, 80, 90]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.min(axis=1, nonzero=True)
+        da_numpy_values = numpy.array(
+            [1, 10, 20, 30, 40, 50, 60, 70, 80, 90]
+        ).reshape(10, 1).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_max_sparse_axis_0(self):
         dm_data = numpy.random.random((10, 20))
@@ -1001,10 +998,9 @@ def test_max_sparse_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=0))
-        da_scipy_values = numpy.array(dm_data.max(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=0)
+        da_scipy_values = dm_data.max(axis=0)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_dense_axis_0(self):
         dm_data = numpy.random.random((10, 20))
@@ -1015,10 +1011,9 @@ def test_max_dense_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=0))
-        da_scipy_values = numpy.array(dm_data.max(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=0)
+        da_scipy_values = dm_data.max(axis=0)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_axis_0_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -1029,10 +1024,10 @@ def test_max_axis_0_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=0, nonzero=True))
-        da_numpy_values = numpy.array([90, 91, 92, 93, 94,
-                                       95, 96, 97, 98, 99]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.max(axis=0, nonzero=True)
+        da_numpy_values = numpy.array([[90, 91, 92, 93, 94,
+                                        95, 96, 97, 98, 99]]).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_max_sparse_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -1044,10 +1039,9 @@ def test_max_sparse_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=1))
-        da_scipy_values = numpy.array(dm_data.max(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=1)
+        da_scipy_values = dm_data.max(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_dense_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -1058,10 +1052,9 @@ def test_max_dense_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=1))
-        da_scipy_values = numpy.array(dm_data.max(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=1)
+        da_scipy_values = dm_data.max(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_axis_1_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -1072,10 +1065,11 @@ def test_max_axis_1_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=1, nonzero=True))
-        da_numpy_values = numpy.array([9, 19, 29, 39, 49,
-                                       59, 69, 79, 89, 99]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.max(axis=1, nonzero=True)
+        da_numpy_values = numpy.array(
+            [9, 19, 29, 39, 49, 59, 69, 79, 89, 99]
+        ).reshape(10, 1).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_argmin_sparse_axis_0(self):
         dm_data = numpy.random.random((10, 20))
diff --git a/tests/cupyx_tests/scipy_tests/sparse_tests/test_csr.py b/tests/cupyx_tests/scipy_tests/sparse_tests/test_csr.py
index 5a2bbc1c30b..dde701b79f7 100644
--- a/tests/cupyx_tests/scipy_tests/sparse_tests/test_csr.py
+++ b/tests/cupyx_tests/scipy_tests/sparse_tests/test_csr.py
@@ -1036,10 +1036,9 @@ def test_min_sparse_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=0))
-        da_scipy_values = numpy.array(dm_data.min(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=0)
+        da_scipy_values = dm_data.min(axis=0)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_min_dense_axis_0(self):
         dm_data = numpy.random.random((10, 20))
@@ -1050,10 +1049,9 @@ def test_min_dense_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=0))
-        da_scipy_values = numpy.array(dm_data.min(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=0)
+        da_scipy_values = dm_data.min(axis=0)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_min_axis_0_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -1064,10 +1062,10 @@ def test_min_axis_0_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=0, nonzero=True))
-        da_numpy_values = numpy.array([10, 1, 2, 3, 4,
-                                       5, 6, 7, 8, 9]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.min(axis=0, nonzero=True)
+        da_numpy_values = numpy.array([[10, 1, 2, 3, 4,
+                                        5, 6, 7, 8, 9]]).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_min_sparse_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -1079,10 +1077,9 @@ def test_min_sparse_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=1))
-        da_scipy_values = numpy.array(dm_data.min(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=1)
+        da_scipy_values = dm_data.min(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_min_dense_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -1093,10 +1090,9 @@ def test_min_dense_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=1))
-        da_scipy_values = numpy.array(dm_data.min(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=1)
+        da_scipy_values = dm_data.min(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_min_axis_1_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -1107,10 +1103,11 @@ def test_min_axis_1_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=1, nonzero=True))
-        da_numpy_values = numpy.array([1, 10, 20, 30, 40,
-                                       50, 60, 70, 80, 90]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.min(axis=1, nonzero=True)
+        da_numpy_values = numpy.array(
+            [1, 10, 20, 30, 40, 50, 60, 70, 80, 90]
+        ).reshape(10, 1).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_max_sparse_axis_0(self):
         dm_data = numpy.random.random((10, 20))
@@ -1122,10 +1119,9 @@ def test_max_sparse_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=0))
-        da_scipy_values = numpy.array(dm_data.max(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=0)
+        da_scipy_values = dm_data.max(axis=0)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_dense_axis_0(self):
         dm_data = numpy.random.random((10, 20))
@@ -1136,10 +1132,9 @@ def test_max_dense_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=0))
-        da_scipy_values = numpy.array(dm_data.max(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=0)
+        da_scipy_values = dm_data.max(axis=0)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_axis_0_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -1150,10 +1145,10 @@ def test_max_axis_0_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=0, nonzero=True))
-        da_numpy_values = numpy.array([90, 91, 92, 93, 94,
-                                       95, 96, 97, 98, 99]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.max(axis=0, nonzero=True)
+        da_numpy_values = numpy.array([[90, 91, 92, 93, 94,
+                                        95, 96, 97, 98, 99]]).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_max_sparse_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -1165,10 +1160,9 @@ def test_max_sparse_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=1))
-        da_scipy_values = numpy.array(dm_data.max(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=1)
+        da_scipy_values = dm_data.max(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_dense_axis_1(self):
         dm_data = numpy.random.random((10, 20))
@@ -1179,10 +1173,9 @@ def test_max_dense_axis_1(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=1))
-        da_scipy_values = numpy.array(dm_data.max(axis=1).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.max(axis=1)
+        da_scipy_values = dm_data.max(axis=1)
+        testing.assert_array_equal(da_cupy_values.A, da_scipy_values.A)
 
     def test_max_axis_1_nonzero(self):
         dm_data = numpy.arange(0, 100, 1).reshape((10, 10)).astype(float)
@@ -1193,10 +1186,11 @@ def test_max_axis_1_nonzero(self):
                                        cupy.array(dm_sparse.indptr)),
                                       shape=(10, 10))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.max(axis=1, nonzero=True))
-        da_numpy_values = numpy.array([9, 19, 29, 39, 49,
-                                       59, 69, 79, 89, 99]).astype(float)
-        assert numpy.array_equal(da_cupy_values, da_numpy_values)
+        da_cupy_values = cp_matrix.max(axis=1, nonzero=True)
+        da_numpy_values = numpy.array(
+            [9, 19, 29, 39, 49, 59, 69, 79, 89, 99]
+        ).reshape(10, 1).astype(float)
+        testing.assert_array_equal(da_cupy_values.A, da_numpy_values)
 
     def test_argmin_sparse_axis_0(self):
         dm_data = numpy.random.random((10, 20))


Code review comments:
diff_hunk: @@ -915,10 +915,9 @@ def test_min_sparse_axis_0(self):
                                        cupy.array(dm_data.indptr)),
                                       shape=(10, 20))
 
-        da_cupy_values = cupy.asnumpy(cp_matrix.min(axis=0))
-        da_scipy_values = numpy.array(dm_data.min(axis=0).todense().ravel())
-        da_scipy_values = da_scipy_values[0, :]
-        assert numpy.array_equal(da_cupy_values, da_scipy_values)
+        da_cupy_values = cp_matrix.min(axis=0)
+        da_scipy_values = dm_data.min(axis=0)
+        testing.assert_array_equal(da_scipy_values.A, da_scipy_values.A)
path: tests/cupyx_tests/scipy_tests/sparse_tests/test_csc.py
text: @user1:
`testing.assert_array_equal` does not support sparse, while `testing.numpy_cupy_array_equal` and `testing.numpy_cupy_allclose` do.

@author:
I make the tests use `testing.numpy_cupy_arrray_equal`. Would you check them? To test `nonzero` parameter, which SciPy does not have yet https://github.com/scipy/scipy/pull/11899 and will be marked as experimental in https://github.com/cupy/cupy/pull/3583, I supply a bit messy SciPy inputs.