Repository: scipy/scipy
Pull request: BUG: sparse: enable sparray for sparse.linalg.norm
Fixes #21690 

This allow spasre arrays (sparray) in `sparse.linalg.norm`. As reported in #21690 this wasn't working because the indexing of the 1D results from (max/sum) reductions didn't align with later reductions.

Note that tests already existed for the case when `axis` was specified as an integer. This bug only arises for `axis` as a 2-tuple or None.

Added tests to test this case.   

It is pretty clear that more testing of sparse.linalg and sparse.csgraph is needed for sparray. But that's a different PR.

Problem statement:
BUG: sparse.linalg: calling `norm` on `sparray` does not work for 1-norms
### Describe your issue.

I'm implementing code for sparse matrices and decided to use `sparray` instead of `spmatrix` (trying to preempt the shift towards `sparray` here). I want to compute the 1-norm of a sparse matrix, however I get the error shown below. Similarly, the norm orders `inf`, `-inf`, `0`, `-1` fail to compute.

It seems to me that the reduction in `mat.sum(axis=0)` from a 2D shape to a 1D shape messes up the call to `max(axis=1)`.

### Reproducing Code Example

```python
import scipy.sparse as sps

mat = sps.random_array((100, 100), density=0.1)
sps.linalg.norm(mat, ord=1)
```


### Error message

```shell
Traceback (most recent call last):
  File "[path-to-my-code]/scipy_norm_sparray_bug.py", line 5, in <module>
    sps.linalg.norm(mat, ord=1)
  File "[path-to-my-venv]/lib/python3.10/site-packages/scipy/sparse/linalg/_norm.py", line 151, in norm
    return abs(x).sum(axis=row_axis).max(axis=col_axis)[0,0]
  File "[path-to-my-venv]/lib/python3.10/site-packages/numpy/core/_methods.py", line 41, in _amax
    return umr_maximum(a, axis, None, out, keepdims, initial, where)
numpy.exceptions.AxisError: axis 1 is out of bounds for array of dimension 1
```


### SciPy/NumPy/Python version and system information

```shell
1.14.1 1.26.4 sys.version_info(major=3, minor=10, micro=11, releaselevel='final', serial=0)
[path-to-my-venv]/lib/python3.10/site-packages/scipy/__config__.py:154: UserWarning: Install `pyyaml` for better output
  warnings.warn("Install `pyyaml` for better output", stacklevel=1)
{
  "Compilers": {
    "c": {
      "name": "clang",
      "linker": "ld64",
      "version": "15.0.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.0.11",
      "commands": "cython"
    },
    "c++": {
      "name": "clang",
      "linker": "ld64",
      "version": "15.0.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld64",
      "version": "13.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.16.1",
      "include directory": "../../../../../../private/var/folders/hw/1f0gcr8d6kn9ms0_wn0_57qc0000gn/T/pip-build-env-mmy0d4vk/overlay/lib/python3.10/site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    },
    "build": {
      "cpu": "aarch64",
      "family": "aarch64",
      "endian": "little",
      "system": "darwin"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "accelerate",
      "found": true,
      "version": "unknown",
      "detection method": "extraframeworks",
      "include directory": "unknown",
      "lib directory": "unknown",
      "openblas configuration": "unknown",
      "pc file directory": "unknown"
    },
    "lapack": {
      "name": "accelerate",
      "found": true,
      "version": "unknown",
      "detection method": "extraframeworks",
      "include directory": "unknown",
      "lib directory": "unknown",
      "openblas configuration": "unknown",
      "pc file directory": "unknown"
    },
    "pybind11": {
      "name": "pybind11",
      "version": "2.12.0",
      "detection method": "config-tool",
      "include directory": "unknown"
    }
  },
  "Python Information": {
    "path": "/private/var/folders/hw/1f0gcr8d6kn9ms0_wn0_57qc0000gn/T/cibw-run-d9imboit/cp310-macosx_arm64/build/venv/bin/python",
    "version": "3.10"
  }
}
```

Hints from issue discussion:
I can verify this bug.  And It looks like there is a fairly easy fix because the axis for the second function in each of these types of norms doesn't need to be specified.  
For example, for `ord=1`,
```python
abs(x).sum(axis=row_axis).max().item()
```
works for both spmatrix and sparray.

This shows that we have work to do to get `sparse.linalg` and `sparse.csgraph` compatible with the sparse array interface.
Thank you @peoe !!

Commit message: handle sparse.linalg.norm rows/cols better for sparray

Patch to review:
diff --git a/scipy/_lib/boost_math b/scipy/_lib/boost_math
index 260f7600e95d..54e5acf7faea 160000
--- a/scipy/_lib/boost_math
+++ b/scipy/_lib/boost_math
@@ -1 +1 @@
-Subproject commit 260f7600e95dcf87b228d46a4a37cf594df5437b
+Subproject commit 54e5acf7faeab1400031d84bbfb2c46e81d82442
diff --git a/scipy/sparse/linalg/_norm.py b/scipy/sparse/linalg/_norm.py
index 4cd7fd6f50dc..021d3fd9afbd 100644
--- a/scipy/sparse/linalg/_norm.py
+++ b/scipy/sparse/linalg/_norm.py
@@ -123,7 +123,7 @@ def norm(x, ord=None, axis=None):
     x = x.tocsr()
 
     if axis is None:
-        axis = (0, 1)
+        axis = tuple(range(x.ndim))
     elif not isinstance(axis, tuple):
         msg = "'axis' must be None, an integer or a tuple of integers"
         try:
@@ -134,7 +134,7 @@ def norm(x, ord=None, axis=None):
             raise TypeError(msg)
         axis = (int_axis,)
 
-    nd = 2
+    nd = x.ndim
     if len(axis) == 2:
         row_axis, col_axis = axis
         if not (-nd <= row_axis < nd and -nd <= col_axis < nd):
@@ -150,13 +150,13 @@ def norm(x, ord=None, axis=None):
             raise NotImplementedError
             #return _multi_svd_norm(x, row_axis, col_axis, amin)
         elif ord == 1:
-            return abs(x).sum(axis=row_axis).max(axis=col_axis)[0,0]
+            return abs(x).sum(axis=row_axis).max().item()
         elif ord == np.inf:
-            return abs(x).sum(axis=col_axis).max(axis=row_axis)[0,0]
+            return abs(x).sum(axis=col_axis).max().item()
         elif ord == -1:
-            return abs(x).sum(axis=row_axis).min(axis=col_axis)[0,0]
+            return abs(x).sum(axis=row_axis).min().item()
         elif ord == -np.inf:
-            return abs(x).sum(axis=col_axis).min(axis=row_axis)[0,0]
+            return abs(x).sum(axis=col_axis).min().item()
         elif ord in (None, 'f', 'fro'):
             # The axis order does not matter for this norm.
             return _sparse_frobenius_norm(x)
diff --git a/scipy/sparse/linalg/tests/test_norm.py b/scipy/sparse/linalg/tests/test_norm.py
index 96c2f65da75b..b5274b315fa9 100644
--- a/scipy/sparse/linalg/tests/test_norm.py
+++ b/scipy/sparse/linalg/tests/test_norm.py
@@ -12,16 +12,26 @@
 
 
 # https://github.com/scipy/scipy/issues/16031
+# https://github.com/scipy/scipy/issues/21690
 def test_sparray_norm():
     row = np.array([0, 0, 1, 1])
     col = np.array([0, 1, 2, 3])
     data = np.array([4, 5, 7, 9])
     test_arr = scipy.sparse.coo_array((data, (row, col)), shape=(2, 4))
     test_mat = scipy.sparse.coo_matrix((data, (row, col)), shape=(2, 4))
-    assert_equal(spnorm(test_arr, ord=1, axis=0), np.array([4, 5, 7, 9]))
-    assert_equal(spnorm(test_mat, ord=1, axis=0), np.array([4, 5, 7, 9]))
-    assert_equal(spnorm(test_arr, ord=1, axis=1), np.array([9, 16]))
-    assert_equal(spnorm(test_mat, ord=1, axis=1), np.array([9, 16]))
+    for ord in (1, np.inf, None):
+        for ax in [0, 1, None, (0, 1), (1, 0)]:
+            for A in (test_arr, test_mat):
+                expected = npnorm(A.toarray(), ord=ord, axis=ax)
+                assert_equal(spnorm(A, ord=ord, axis=ax), expected)
+    # test 1d
+    test_arr_1d = scipy.sparse.coo_array((data, (col,)), shape=(4,))
+    test_mat_col = scipy.sparse.coo_matrix((data, (col, [0, 0, 0, 0])), shape=(4, 1))
+    for ord in (1, np.inf, None):
+        for ax in [0, None]:
+            for A in (test_arr_1d, test_mat_col):
+                expected = npnorm(A.toarray(), ord=ord, axis=ax)
+                assert_equal(spnorm(A, ord=ord, axis=ax), expected)
 
 
 class TestNorm:
@@ -33,7 +43,7 @@ def setup_method(self):
     def test_matrix_norm(self):
 
         # Frobenius norm is the default
-        assert_allclose(spnorm(self.b), 7.745966692414834)        
+        assert_allclose(spnorm(self.b), 7.745966692414834)
         assert_allclose(spnorm(self.b, 'fro'), 7.745966692414834)
 
         assert_allclose(spnorm(self.b, np.inf), 9)
@@ -50,7 +60,7 @@ def test_matrix_norm(self):
 
     def test_matrix_norm_axis(self):
         for m, axis in ((self.b, None), (self.b, (0, 1)), (self.b.T, (1, 0))):
-            assert_allclose(spnorm(m, axis=axis), 7.745966692414834)        
+            assert_allclose(spnorm(m, axis=axis), 7.745966692414834)
             assert_allclose(spnorm(m, 'fro', axis=axis), 7.745966692414834)
             assert_allclose(spnorm(m, np.inf, axis=axis), 9)
             assert_allclose(spnorm(m, -np.inf, axis=axis), 2)


Code review comments:
path: scipy/_lib/boost_math
text: @user1:
Accidental submodule update here?

Some CI failures look real as well, but I'll leave that to the usual `sparse` experts.