BUG: fix pow method for sparrays with power zero (scipy#19335)

* fix pow method for sparrays with power 0

A**0 means every entry is 1, not just the originally non-zero entries.
So, the A**x function changes the nonzero structure to dense in a
discontinuous fashion near x=0. There is a question then as to whether
this method should return a sparse array with all entries 1, or return
a dense array, or raise an exception with a workaround described.

This PR implements raising an exception for `n` a zero valued scalar.
The message (and doc_string description of `n`) says if they want to
compute A**x for x==0, they should special case that value to use
`np.ones(A.shape, dtype=A.dtype)`

* Update scipy/sparse/tests/test_array_api.py

Co-authored-by: Jake Bowhay <[email protected]>

* change from ValueError to NotImplementedError and update test

* Thanks for these suggestions:)

Co-authored-by: Julien Jerphanion <[email protected]>

* fix error message

---------

Co-authored-by: Jake Bowhay <[email protected]>
Co-authored-by: Julien Jerphanion <[email protected]>

Author: 3 people

scipy/sparse/_data.py

@@ -107,12 +107,25 @@ def power(self, n, dtype=None):
 
         Parameters
         ----------
-        n : n is a scalar
+        n : scalar
+            n is a non-zero scalar (nonzero avoids dense ones creation)
+            If zero power is desired, special case it to use `np.ones`
 
         dtype : If dtype is not specified, the current dtype will be preserved.
+
+        Raises
+        ------
+        NotImplementedError : if n is a zero scalar
+            If zero power is desired, special case it to use
+            `np.ones(A.shape, dtype=A.dtype)`
         """
         if not isscalarlike(n):
             raise NotImplementedError("input is not scalar")
+        if not n:
+            raise NotImplementedError(
+                "zero power is not supported as it would densify the matrix.\n"
+                "Use `np.ones(A.shape, dtype=A.dtype)` for this case."
+            )
 
         data = self._deduped_data()
         if dtype is not None:

scipy/sparse/tests/test_array_api.py

@@ -140,10 +140,16 @@ def test_matmul(A):
     assert np.all((A @ A.T).todense() == A.dot(A.T).todense())
 
 
-@parametrize_square_sparrays
-def test_pow(B):
-    assert isinstance((B**0), scipy.sparse.sparray), "Expected array, got matrix"
-    assert isinstance((B**2), scipy.sparse.sparray), "Expected array, got matrix"
+@parametrize_sparrays
+def test_power_operator(A):
+    assert isinstance((A**2), scipy.sparse.sparray), "Expected array, got matrix"
+
+    # https://github.com/scipy/scipy/issues/15948
+    npt.assert_equal((A**2).todense(), (A.todense())**2)
+
+    # power of zero is all ones (dense) so helpful msg exception
+    with pytest.raises(NotImplementedError, match="zero power"):
+        A**0
 
 
 @parametrize_sparrays
@@ -354,12 +360,6 @@ def test_spilu():
     )
 
 
-@parametrize_sparrays
-def test_power_operator(A):
-    # https://github.com/scipy/scipy/issues/15948
-    npt.assert_equal((A**2).todense(), (A.todense())**2)
-
-
 @pytest.mark.parametrize(
     "cls,indices_attrs",
     [