Repository: scipy/scipy
Pull request: MAINT: Use `sparse.diags` instead of `spdiags` internally.
Modifies `sparse.eye()` to use `sparse.diags` internally instead of `sparse.spdiags`. This modifies the behavior in the case where the diagonal offset(s) are greater than the maximum shape of the array. Previously, this would result in a sparse array with 0 stored elements and, in the case of the DIA format, negative values for the  diagonals:

```python
# Current behavior
In [1]: import scipy as sp

In [2]: a = sp.sparse.eye(3, 3, -5)

In [3]: a
Out[3]: 
<3x3 sparse matrix of type '<class 'numpy.float64'>'
        with -2 stored elements (1 diagonals) in DIAgonal format>

In [4]: a.toarray()
Out[4]: 
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
```

With this PR, this scenario will instead raise a `ValueError` indicating that the resulting array would be empty due to the offsets being OOB:

```
# This PR
In [1]: import scipy as sp

In [2]: a = sp.sparse.eye(3, 3, -5)
Traceback (most recent call last)
   ...
ValueError: Offset -5 (index 0) out of bounds
```

#### Reference issue
Closes #18555 

#### Additional information
This proposal is predicated on the assumption that users generally don't intend to create empty sparse arrays with `eye()`/`identity()`. However, if this extra input validation is not desired, then an alternative approach would be to call the `dia_array` constructor directly instead of using `spdiags`. This would maintain the current behavior while getting rid of the dis-recommended usage of `spdiags` internally.

Removing the internal usage of `spdiags` will reduce the amount of duplication necessary to add new sparse array creation functions (e.g. `eye_array`).

Problem statement:
BUG: `spdiags`, `eye`, `identity` corner case when offset greater than array shape
### Describe your issue.

The `spdiags` creation function allows the diagonal offsets to be greater than the matrix shape, which results in an empty array in `dia` format with a negative number of non-zero values:

```python
>>> m = n = 3  # Square, 3x3
>>> k = -5  # diagonal offset
>>> data = np.ones(m)
>>> a = sp.sparse.spdiags(data, diags=k, m=m, n=n)
>>> a
<3x3 sparse matrix of type '<class 'numpy.float64'>'
        with -2 stored elements (1 diagonals) in DIAgonal format>
>>> a.toarray()  # 3x3 array of zeros
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
```

Both `identity` and `eye` are implemented using `spdiags` (at least for some formats), which means that this behavior extends to those functions as well.

This behavior is not strictly incorrect, but is likely not what users want nor expect. `diags`, which the documentation notes is preferred over `spdiags`, explicitly validates that the diagonal offset is less than the minimum dimension of the requested shape:

```python
>>> sp.sparse.diags([data], offsets=[k], shape=(m, n))
Traceback (most recent call last)
   ...
ValueError: Offset -5 (index 0) out of bounds
```

### Proposed change

I propose to modify `eye` and `identity` to use `diags` instead of `spdiags`. This will change behavior so that `eye` and `identity` will raise when the diagonal offset is greater than the shape.

### Reproducing Code Example

```python
See above
```


### Error message

```shell
N/A
```


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

```shell
1.11.0.dev0+2204.675b17e 1.24.3 sys.version_info(major=3, minor=11, micro=3, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/include
    lib directory: /usr/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH= DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK=0
      NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=1 HASWELL MAX_THREADS=8
    pc file directory: /usr/lib/pkgconfig
    version: 0.3.21
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/include
    lib directory: /usr/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH= DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK=0
      NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=1 HASWELL MAX_THREADS=8
    pc file directory: /usr/lib/pkgconfig
    version: 0.3.21
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.10.4
Compilers:
  c:
    commands: ccache, cc
    linker: ld.bfd
    name: gcc
    version: 13.1.1
  c++:
    commands: ccache, c++
    linker: ld.bfd
    name: gcc
    version: 13.1.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 0.29.35
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 13.1.1
  pythran:
    include directory: ../../../.virtualenvs/sp-dev/lib/python3.11/site-packages/pythran
    version: 0.13.1
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /home/ross/.virtualenvs/sp-dev/bin/python
  version: '3.11'
```

Hints from issue discussion:
+1, this would be a nice improvement.

Aside: It's unfortunate that the DIA `nnz` semantics don't really match the rest of the formats, where `nnz` reflects the number of stored entries (including duplicates). If you're using `nnz` to keep track of the storage cost of a sparse array/matrix, then the DIA result isn't what you want.

Commit message: ENH: Prefer sparse.diags over spdiags for eye().

Switches the implementation of sparse.eye to use sparse.diags
instead of spdiags. This eliminates a corner-case where spdiags
would result in dia arrays with negative diagonals when the
diagonal offset was greater than the array shape.

Patch to review:
diff --git a/scipy/sparse/_construct.py b/scipy/sparse/_construct.py
index 95ac89cb1a01..c685c11c7741 100644
--- a/scipy/sparse/_construct.py
+++ b/scipy/sparse/_construct.py
@@ -348,8 +348,8 @@ def eye(m, n=None, k=0, dtype=float, format=None):
             data = np.ones(n, dtype=dtype)
             return coo_matrix((data, (row, col)), (n, n))
 
-    diags = np.ones((1, max(0, min(m + k, n))), dtype=dtype)
-    return spdiags(diags, k, m, n).asformat(format)
+    data = np.ones((1, max(0, min(m + k, n))), dtype=dtype)
+    return diags(data, offsets=[k], shape=(m, n), dtype=dtype).asformat(format)
 
 
 def kron(A, B, format=None):
diff --git a/scipy/sparse/tests/test_construct.py b/scipy/sparse/tests/test_construct.py
index 37c5543acd20..a7400374c269 100644
--- a/scipy/sparse/tests/test_construct.py
+++ b/scipy/sparse/tests/test_construct.py
@@ -248,9 +248,27 @@ def test_eye(self):
         for m in [3, 5]:
             for n in [3, 5]:
                 for k in range(-5,6):
-                    assert_equal(construct.eye(m, n, k=k).toarray(), np.eye(m, n, k=k))
-                    if m == n:
-                        assert_equal(construct.eye(m, k=k).toarray(), np.eye(m, n, k=k))
+                    # scipy.sparse.eye deviates from np.eye here. np.eye will
+                    # create arrays of all 0's when the diagonal offset is
+                    # greater than the size of the array. For sparse arrays
+                    # this makes less sense, especially as it results in dia
+                    # arrays with negative diagonals. Therefore sp.sparse.eye
+                    # validates that diagonal offsets fall within the shape of
+                    # the array. See gh-18555.
+                    if abs(k) > m or k > n:
+                        pytest.raises(
+                            ValueError, match="Offset.*out of bounds"
+                        )
+                    else:
+                        assert_equal(
+                            construct.eye(m, n, k=k).toarray(),
+                            np.eye(m, n, k=k)
+                        )
+                        if m == n:
+                            assert_equal(
+                                construct.eye(m, k=k).toarray(),
+                                np.eye(m, n, k=k)
+                            )
 
     def test_eye_one(self):
         assert_equal(construct.eye(1).toarray(), [[1]])


Code review comments:
diff_hunk: @@ -248,9 +248,27 @@ def test_eye(self):
         for m in [3, 5]:
             for n in [3, 5]:
                 for k in range(-5,6):
-                    assert_equal(construct.eye(m, n, k=k).toarray(), np.eye(m, n, k=k))
-                    if m == n:
-                        assert_equal(construct.eye(m, k=k).toarray(), np.eye(m, n, k=k))
+                    # scipy.sparse.eye deviates from np.eye here. np.eye will
+                    # create arrays of all 0's when the diagonal offset is
+                    # greater than the size of the array. For sparse arrays
+                    # this makes less sense, especially as it results in dia
+                    # arrays with negative diagonals. Therefore sp.sparse.eye
+                    # validates that diagonal offsets fall within the shape of
+                    # the array. See gh-18555.
+                    if abs(k) > m or k > n:
+                        pytest.raises(
+                            ValueError, match="Offset.*out of bounds"
+                        )
path: scipy/sparse/tests/test_construct.py
text: @user1:
I think you need to fun the function to get the ValueError

@author:
Of course - nice catch!

I should've been suspicious when it didn't raise on the `match` locally, since I never get regex's right on my first attempt :). Fixed in ead4a42 by using `pytest.raises` as a context manager. This also highlighted a bit of tightening that was necessary on the conditional for determining when a diagonal is OOB. Thanks @user1 !