Repository: cupy/cupy
Pull request: Support sparse pointwise division by vectors or matrices
Closes https://github.com/cupy/cupy/issues/3830

This PR is to support sparse matrix point-wise division by another matrix or vector.

Problem statement:
Sparse pointwise division
This is a very important feature as right now there is no way to perform point wise division of vectors, dense or sparse matrices with other sparse matrices. As a workaround, we have been multiplying by the reciprocal, but this is not immediately obvious to most users and adds complexity.

Pointwise multiplication was recently added to sparse matrices. It shouldn't be too hard to repurpose the multiplication to a division operation. 

Tagging @anaruse

Hints from issue discussion:

Commit message: Support sparse pointwise division by vectors or matrices

Patch to review:
diff --git a/cupyx/scipy/sparse/csr.py b/cupyx/scipy/sparse/csr.py
index 51806930a8f..7290c28f5ee 100644
--- a/cupyx/scipy/sparse/csr.py
+++ b/cupyx/scipy/sparse/csr.py
@@ -183,7 +183,7 @@ def __rdiv__(self, other):
         raise NotImplementedError
 
     def __truediv__(self, other):
-        """Point-wise division by scalar"""
+        """Point-wise division by another matrix, vector or scalar"""
         if util.isscalarlike(other):
             if self.dtype == numpy.complex64:
                 # Note: This is a work-around to make the output dtype the same
@@ -196,7 +196,16 @@ def __truediv__(self, other):
                     dtype = numpy.promote_types(numpy.float64, other.dtype)
             d = cupy.array(1. / other, dtype=dtype)
             return multiply_by_scalar(self, d)
-        # TODO(anaruse): Implement divide by dense or sparse matrix
+        elif util.isdense(other):
+            other = cupy.atleast_2d(other)
+            check_shape_for_pointwise_op(self.shape, other.shape)
+            return self.todense() / other
+        elif base.isspmatrix(other):
+            check_shape_for_pointwise_op(self.shape, other.shape)
+            dtype = numpy.promote_types(self.dtype, other.dtype)
+            if dtype.char not in 'FD':
+                dtype = numpy.promote_types(numpy.float64, dtype)
+            return self.todense().astype(dtype) / other.todense().astype(dtype)
         raise NotImplementedError
 
     def __rtruediv__(self, other):
@@ -477,6 +486,15 @@ def isspmatrix_csr(x):
     return isinstance(x, csr_matrix)
 
 
+def check_shape_for_pointwise_op(a_shape, b_shape):
+    a_m, a_n = a_shape
+    b_m, b_n = b_shape
+    if not (a_m == b_m or a_m == 1 or b_m == 1):
+        raise ValueError('inconsistent shape')
+    if not (a_n == b_n or a_n == 1 or b_n == 1):
+        raise ValueError('inconsistent shape')
+
+
 def multiply_by_scalar(sp, a):
     data = sp.data * a
     indices = sp.indices.copy()
@@ -485,12 +503,9 @@ def multiply_by_scalar(sp, a):
 
 
 def multiply_by_dense(sp, dn):
+    check_shape_for_pointwise_op(sp.shape, dn.shape)
     sp_m, sp_n = sp.shape
     dn_m, dn_n = dn.shape
-    if not (sp_m == dn_m or sp_m == 1 or dn_m == 1):
-        raise ValueError('inconsistent shape')
-    if not (sp_n == dn_n or sp_n == 1 or dn_n == 1):
-        raise ValueError('inconsistent shape')
     m, n = max(sp_m, dn_m), max(sp_n, dn_n)
     nnz = sp.nnz * (m // sp_m) * (n // sp_n)
     dtype = numpy.promote_types(sp.dtype, dn.dtype)
@@ -567,12 +582,9 @@ def cupy_multiply_by_dense():
 
 
 def multiply_by_csr(a, b):
+    check_shape_for_pointwise_op(a.shape, b.shape)
     a_m, a_n = a.shape
     b_m, b_n = b.shape
-    if not (a_m == b_m or a_m == 1 or b_m == 1):
-        raise ValueError('inconsistent shape')
-    if not (a_n == b_n or a_n == 1 or b_n == 1):
-        raise ValueError('inconsistent shape')
     m, n = max(a_m, b_m), max(a_n, b_n)
     a_nnz = a.nnz * (m // a_m) * (n // a_n)
     b_nnz = b.nnz * (m // b_m) * (n // b_n)
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 7f3c6d2d30b..b9872f29a67 100644
--- a/tests/cupyx_tests/scipy_tests/sparse_tests/test_csr.py
+++ b/tests/cupyx_tests/scipy_tests/sparse_tests/test_csr.py
@@ -124,7 +124,7 @@ def _make_col(xp, sp, dtype):
 
 
 def _make_shape(xp, sp, dtype):
-    return sp.csr_matrix((3, 4))
+    return sp.csr_matrix((3, 4), dtype=dtype)
 
 
 @testing.parameterize(*testing.product({
@@ -1106,6 +1106,30 @@ def test_divide_scalarlike(self, xp, sp):
         y = m / xp.array(2)
         return y.toarray()
 
+    @testing.numpy_cupy_allclose(sp_name='sp')
+    def test_divide_dense_row(self, xp, sp):
+        m = self.make(xp, sp, self.dtype)
+        x = xp.arange(4, dtype=self.dtype)
+        return m / x
+
+    @testing.numpy_cupy_allclose(sp_name='sp')
+    def test_divide_dense_col(self, xp, sp):
+        m = self.make(xp, sp, self.dtype)
+        x = xp.arange(3, dtype=self.dtype).reshape(3, 1)
+        return m / x
+
+    @testing.numpy_cupy_allclose(sp_name='sp')
+    def test_divide_dense_matrix(self, xp, sp):
+        m = self.make(xp, sp, self.dtype)
+        x = xp.arange(12, dtype=self.dtype).reshape(3, 4)
+        return m / x
+
+    @testing.numpy_cupy_allclose(sp_name='sp')
+    def test_divide_csr_matrix(self, xp, sp):
+        m = self.make(xp, sp, self.dtype)
+        x = _make4(xp, sp, self.dtype)
+        return m / x
+
 
 @testing.parameterize(*testing.product({
     'dtype': [numpy.float32, numpy.float64, numpy.complex64, numpy.complex128],


Code review comments:
diff_hunk: @@ -477,6 +486,15 @@ def isspmatrix_csr(x):
     return isinstance(x, csr_matrix)
 
 
+def check_shape_for_pointwise_op(a_shape, b_shape):
path: cupyx/scipy/sparse/csr.py
text: @user1:
I noticed SciPy doesn't check the input shapes. Can't we just let the division raise the error if the shapes are incompatible?

@user1:
(I do see this was implemented elsewhere and this is just a refactoring, so I'm simply wondering and not requesting a change 😅)

@author:
The input shapre is checked here.
https://github.com/scipy/scipy/blob/3af12341854f9f9a6ef0a3a9c853e80bcc5e8002/scipy/sparse/compressed.py#L1249-L1250
The difference is that this PR allows broadcasting (like point-wise multiply), while the SciPy does not..

@user2:
> The difference is that this PR allows broadcasting (like point-wise multiply), while the SciPy does not..

SciPy does not allow broadcasting on `sparse / sparse`, while `sparse.multiply(sparse)` allows.  Could you make the PR compatible with SciPy?

@author:
OK, I will fix it to raise an exception when broadcasting is needed for compatibility with SciPy.

diff_hunk: @@ -196,7 +196,16 @@ def __truediv__(self, other):
                     dtype = numpy.promote_types(numpy.float64, other.dtype)
             d = cupy.array(1. / other, dtype=dtype)
             return multiply_by_scalar(self, d)
-        # TODO(anaruse): Implement divide by dense or sparse matrix
+        elif util.isdense(other):
+            other = cupy.atleast_2d(other)
+            check_shape_for_pointwise_op(self.shape, other.shape)
+            return self.todense() / other
+        elif base.isspmatrix(other):
+            check_shape_for_pointwise_op(self.shape, other.shape)
+            dtype = numpy.promote_types(self.dtype, other.dtype)
+            if dtype.char not in 'FD':
+                dtype = numpy.promote_types(numpy.float64, dtype)
+            return self.todense().astype(dtype) / other.todense().astype(dtype)
path: cupyx/scipy/sparse/csr.py
text: @user1:
Could you add a comment for calling `todense()`? You can also link to https://github.com/cupy/cupy/issues/3830#issuecomment-678928623 if you like.

@user1:
Looks like SciPy calls a specialized routine when `other` is also a sparse matrix:
https://github.com/scipy/scipy/blob/3af12341854f9f9a6ef0a3a9c853e80bcc5e8002/scipy/sparse/compressed.py#L1245
Maybe we can leave a TODO here?

@author:
That's right, SciPy uses a specialized routine to create the intermediate result in a sparse matrix if `other` is sparse matrix. However, if the data type is a floating point, it is converted to a dense matrix to show that the result of division by zero is NaN. Therefore, if you convert the sparse matrix to a dense matrix beforehand and then perform the division (as this PR does), I think the result will not change.

@user2:
Does the following reduce memory usage?
```python
dense = self.todense().astype(dtype, copy=False)
dense /= other.todense()
return dense
```