Repository: scipy/scipy
Pull request: BUG: fix pow method for sparrays with power zero
Even after the fix #16430 fixed #15948 so the sparse array ** operator matches numpy **, they don't match when the power is `0`. The result of `A ** 0` is `np.ones(A.shape)` for a ndarray, but sparse arrays set `A ** 0` to ones only for originally non-zero values.  Thus 
```python
A = sp.sparse.csr_array([[0,2,-3]])
result = (A ** 0).toarray()         # => np.array([[0, 1, 1]])

np.array([[0, 2, -3]] ** 0  #  =>  np.array([[1, 1, 1]])
```

`A**0` should mean every entry is `1`, not just the originally non-zero entries. So the current return values are incorrect for sparse arrays (BUG).

Furthermore, this means the nonzero structure of `A**x` changes 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 switch to a dense array, or raise an exception with a workaround described.

This PR implements raising a `ValueError` when `n` is a zero valued scalar. The message (and doc_string description of `n`) says if they want to compute `A**n` for `n==0`, they should special case that value to use `np.ones(A.shape, dtype=A.dtype)`

If they really meant to do this, I'm thinking they could use something like:
```python
result = A ** x if x else np.ones(A.shape, dtype=A.dtype)
```
And we could make our function do that too, but I think it might surprise people who don't expect conversion to a dense matrix.  Either way, the current return value of a matrix with some 1 and some 0 values is wrong and doesn't match numpy results.


Problem statement:
BUG: Power operator behaviour for sparse arrays
### Describe your issue.

Currently, applying the `**` operator to a sparse array implements the `numpy.linalg.matrix_power` function. This is not compatible with the `np.array` implementation of `**`, which does element-wise power. This should be fixed as part of the migration from sparse "matrices" to sparse "arrays". Additionally, the `matrix_power` function should then also be included in the `sparse.linalg` package.

### Reproducing Code Example

```python
from scipy.sparse import csr_array, random

A = csr_array(random(5,5,0.1))
print(A.A)  
print((A**2).A)  #This does matrix_power, eventhough it is initialised as a sparse array (not a matrix)
print((A.A)**2)  #This does element-wise power, since it is a numpy array
```


### Error message

```shell
N/A
```


### SciPy/NumPy/Python version information

1.8.0 1.20.3 sys.version_info(major=3, minor=9, micro=7, releaselevel='final', serial=0)

Hints from issue discussion:
My duplicate issue (#15950) is closed, but there is a slight difference which might be important for you, so I just want to comment here. 
This issue is for sparse arrays, and the same problem occurs also with the sparse matrices (which, I guess, is no longer supported ?)

```python
import scipy.sparse as sp

someMatrix = sp.csc_matrix(([1, -1], ([0, 1], [1, 0])), shape=(2, 2))

print(someMatrix.A)                 # print the matrix
print((someMatrix**2).A)        # raises the matrix to a power
print((someMatrix.A)**2)        # element-wise power   
```
These are indeed not duplicates for that reason, but #15950 is also not a bug. For matrix objects, `**` is indeed designed to be matrix power, and for arrays (e.g. `someMatrix.A`), `**` is designed to be elementwise.
Thanks for filing this, @jpdehollain. I agree that the new sparse arrays should treat `**` as an elementwise operation like a regular ndarray.

Fixing this should be a simple matter of adding more overrides to the `*` ones here: https://github.com/scipy/scipy/blob/64dd617cd2c7d1bb34a6e6dfc2ca37bad6ac063b/scipy/sparse/_arrays.py#L45-L50
@perimosocordiae Could this please be assigned to me? Happy to take a look - looks like it'll also need tests added to make sure the overrides work 

Also can I please check if it's possible to build from source with python 3.10? When running ` python runtests.py -v` from root folder with an isolated download of Scipy 1.8 I get `numpy.distutils.system_info.NotFoundError: No BLAS/LAPACK libraries found.`
Go for it! Please tag me when your PR is ready to review.
Hey! Is this issue still up? I would like to take it. Thanks.

Commit message: 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)`

Patch to review:
diff --git a/scipy/sparse/_data.py b/scipy/sparse/_data.py
index 87e5f6534e38..8d137efdefb2 100644
--- a/scipy/sparse/_data.py
+++ b/scipy/sparse/_data.py
@@ -107,12 +107,23 @@ 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
+        ------
+        ValueError : 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 ValueError("zero power makes all ones (dense).\n"
+                "Use `np.ones(A.shape, dtype=A.dtype)` for this case.")
 
         data = self._deduped_data()
         if dtype is not None:
diff --git a/scipy/sparse/tests/test_array_api.py b/scipy/sparse/tests/test_array_api.py
index d92924bb2d79..9b03cb699016 100644
--- a/scipy/sparse/tests/test_array_api.py
+++ b/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(ValueError):
+        A**0
 
 
 @parametrize_sparrays
@@ -346,12 +352,6 @@ def test_spilu():
     npt.assert_allclose(LU.solve(np.array([1, 2, 3, 4])), [1, 0, 0, 0])
 
 
-@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",
     [


Code review comments:
diff_hunk: @@ -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(ValueError):
path: scipy/sparse/tests/test_array_api.py
text: @user1:
```suggestion
    with pytest.raises(ValueError, match="zero power"):
```