Repository: scipy/scipy
Pull request: ENH: Update isspmatrix behavior
Closes #15790

Updates `isspmatrix` so that it is no longer an alias of `issparse`. After #18440, it became possible to define the functions so that `issparse` returns `True` for both sparse arrays and sparse matrices, while `isspmatrix` returns True *only* for sparse matrices and not sparse arrays. This makes the behavior of the functions more accurately reflect their name.

Note that this is a breaking change necessitating updating internal usage of `isspmatrix`. 

Problem statement:
BUG: `isspmatrix` doesn't account for sparse arrays
### Describe your issue.

There are two functions in `scipy.sparse` which determine whether or not an object is a sparse container: `sp.sparse.isspmatrix` and `sp.sparse.issparse`. Based on the names of these functions, one might expect the former to return `True` for the sparse matrix containers (e.g. `csr_matrix`) and the latter to return true for either sparse matrices or sparse arrays.

It turns out that `issparse` is just an alias of `isspmatrix`, and `isspmatrix` returns True for sparse arrays as well as matrices, which might be a misleading.

Given the addition of the sparse array containers, I would propose modifying the behavior of these functions to account for the distinction between sparse arrays and matrices, i.e. something like:

```python
def issparse(x):
    # Note: the sparse array containers are currently subclasses of
    # spmatrix, so this will be true for both arrays and matrices
    return isinstance(x, spmatrix)

def isspmatrix(x):
    return isinstance(x, spmatrix) and not isinstance(x, _sparray)
```

### Reproducing Code Example

```python
>>> import scipy as sp
>>> a = sp.sparse.csr_array(np.eye(3))
>>> m = sp.sparse.csr_matrix(np.eye(3))
>>> sp.sparse.isspmatrix(m)  # Expect True
True
>>> sp.sparse.issparse(m)  # Expect True
True
>>> sp.sparse.issparse(a)  # Expect True
True
>>> sp.sparse.isspmatrix(a)  # Expect False: problem
True
```


### Error message

```shell
N/A
```


### SciPy/NumPy/Python version information

1.9.0.dev0+1669.b9427e1 1.22.3 sys.version_info(major=3, minor=10, micro=2, releaselevel='final', serial=0)

Hints from issue discussion:
Thanks @rossbar. I agree with the proposed semantics for `issparse` and `ismatrix`.

This shows an unfortunate issue though:
```python
def isspmatrix(x):
    return isinstance(x, spmatrix) and not isinstance(x, _sparray)
```
The sparse array classes are subclasses of the matrix classes.  That means existing `isinstance(x, spmatrix)` in the wild are not quite doing the right thing. Perhaps that was discussed in @stefanv's PR, I can't remember. It's not easy to fix, so at least it should be pointed out with a warning not to use `isinstance` near the top of http://scipy.github.io/devdocs/reference/sparse.html I think.
> The sparse array classes are subclasses of the matrix classes. That means existing isinstance(x, spmatrix) in the wild are not quite doing the right thing

Agreed, though fortunately this shouldn't really be breaking existing usages of `isinstance(x, spmatrix)`, since the sparse array classes were only recently added in 1.8.

> It's not easy to fix,

I actually don't think it would be that bad. The main issue is that the `spmatrix` class is what has the concrete implementation of almost all of the methods. One idea would be to introduce a new class, say `class _SparseBase` that both `spmatrix` and `_sparray` derive from. That would clean up the inheritance model and thus resolving this issue as well, e.g.

```python
def issparse(x):
    return isinstance(x, _SparseBase)  # True for both spmatrix and sparray

def isspmatrix(x):
    return isinstance(x, spmatrix)  # Now only True for spmatrix
```
That sounds like a good plan.
I'm not sure we need to introduce _SparseBase to accomplish this, but I agree with the proposed semantics. We should fix this sooner than later to reduce the possibility of breakage.

One side note: it's possible that this will hamper adoption of the array-style classes, as older code written with `isspmatrix` will now no longer work as intended. But maybe that's a good thing, as it's likely the older code would have subtle bugs related to slice shapes, etc.
> The sparse array classes are subclasses of the matrix classes.

I was very emotional after seeing my attempts to explicitly block sparse matrix argument failed, because:
```
isinstance(sparse_array, sp.coo_matrix) is True!
```

This breaks LSP in a bad way.

From practical side, I believe that (co-)inheritance should be avoided here. 

I see main scenarios as these:
- code operates with matrices (will check for `isspmatrix`)
- code operates with arrays (will check for `issparray`)
  - I think code shouldn't work with both at the same time: many differences that are easy to overlook
- code operates with neither (better not providing any check here)
- code accepts 'sparse whatever' only to convert it to dense or their format of choice.

Previously written code used `issparse` to pass by only scipy sparse matrices (as @perimosocordiae noted), and better just deprecate that function.

@arogozhnikov You are right, of course. This is because the current sparse array implementation is a bit of a hack, building on top of matrices. We plan on addressing this in the upcoming months.

/cc @dschult @rossbar

Commit message: MAINT: Rm isspmatrix_fmt usage from csgraph pyx files.

Patch to review:
diff --git a/scipy/sparse/_base.py b/scipy/sparse/_base.py
index 1b3a07e510af..1f95ffa07905 100644
--- a/scipy/sparse/_base.py
+++ b/scipy/sparse/_base.py
@@ -7,6 +7,8 @@
                        get_sum_dtype, isdense, isscalarlike,
                        matrix, validateaxis,)
 
+from ._matrix import spmatrix
+
 __all__ = ['isspmatrix', 'issparse',
            'SparseWarning', 'SparseEfficiencyWarning']
 
@@ -1442,7 +1444,7 @@ def getrow(self, i):
 
 
 def issparse(x):
-    """Is x of a sparse array type?
+    """Is `x` of a sparse array type?
 
     Parameters
     ----------
@@ -1452,23 +1454,46 @@ def issparse(x):
     Returns
     -------
     bool
-        True if x is a sparse array, False otherwise
-
-    Notes
-    -----
-    issparse and isspmatrix are aliases for the same function.
+        True if `x` is a sparse array, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import csr_array, issparse
-    >>> issparse(csr_array([[5]]))
+    >>> from scipy.sparse import csr_array, csr_matrix, issparse
+    >>> issparse(csr_matrix([[5]]))
     True
-
-    >>> from scipy.sparse import issparse
+    >>> issparse(csr_array([[5]])
+    True
+    >>> issparse(np.array([[5]]))
+    False
     >>> issparse(5)
     False
     """
     return isinstance(x, _sparray)
 
 
-isspmatrix = issparse
+def isspmatrix(x):
+    """Is `x` of a sparse matrix type?
+
+    Parameters
+    ----------
+    x
+        object to check for being a sparse matrix
+
+    Returns
+    -------
+    bool
+        True if `x` is a sparse matrix, False otherwise
+
+    Examples
+    --------
+    >>> from scipy.sparse import csr_array, csr_matrix, isspmatrix
+    >>> isspmatrix(csr_matrix([[5]]))
+    True
+    >>> isspmatrix(csr_array([[5]]))
+    False
+    >>> isspmatrix(np.array([[5]]))
+    False
+    >>> isspmatrix(5)
+    False
+    """
+    return isinstance(x, spmatrix)
diff --git a/scipy/sparse/_bsr.py b/scipy/sparse/_bsr.py
index c6be241752ef..c6728f2e8a84 100644
--- a/scipy/sparse/_bsr.py
+++ b/scipy/sparse/_bsr.py
@@ -11,7 +11,7 @@
 from ._matrix import spmatrix, _array_doc_to_matrix
 from ._data import _data_matrix, _minmax_mixin
 from ._compressed import _cs_matrix
-from ._base import isspmatrix, _formats, _sparray
+from ._base import issparse, _formats, _sparray
 from ._sputils import (isshape, getdtype, getdata, to_native, upcast,
                        check_shape)
 from . import _sparsetools
@@ -122,8 +122,8 @@ class bsr_array(_cs_matrix, _minmax_mixin):
     def __init__(self, arg1, shape=None, dtype=None, copy=False, blocksize=None):
         _data_matrix.__init__(self)
 
-        if isspmatrix(arg1):
-            if isspmatrix_bsr(arg1) and copy:
+        if issparse(arg1):
+            if arg1.format == self.format and copy:
                 arg1 = arg1.copy()
             else:
                 arg1 = arg1.tobsr(blocksize=blocksize)
@@ -372,14 +372,12 @@ def _mul_sparse_matrix(self, other):
         R,n = self.blocksize
 
         # convert to this format
-        if isspmatrix_bsr(other):
+        if issparse(other) and other.format == "bsr":
             C = other.blocksize[1]
         else:
             C = 1
 
-        from ._csr import isspmatrix_csr
-
-        if isspmatrix_csr(other) and n == 1:
+        if other.format == "csr" and n == 1:
             other = other.tobsr(blocksize=(n,C), copy=False)  # lightweight conversion
         else:
             other = other.tobsr(blocksize=(n,C))
@@ -697,7 +695,7 @@ def _with_data(self,data,copy=True):
 
 
 def isspmatrix_bsr(x):
-    """Is x of a bsr_array type?
+    """Is `x` of a bsr_matrix type?
 
     Parameters
     ----------
@@ -707,19 +705,19 @@ def isspmatrix_bsr(x):
     Returns
     -------
     bool
-        True if x is a bsr matrix, False otherwise
+        True if `x` is a bsr matrix, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import bsr_array, isspmatrix_bsr
-    >>> isspmatrix_bsr(bsr_array([[5]]))
+    >>> from scipy.sparse import bsr_array, bsr_matrix, csr_matrix, isspmatrix_bsr
+    >>> isspmatrix_bsr(bsr_matrix([[5]])
     True
-
-    >>> from scipy.sparse import bsr_array, csr_matrix, isspmatrix_bsr
+    >>> isspmatrix_bsr(bsr_array([[5]]))
+    False
     >>> isspmatrix_bsr(csr_matrix([[5]]))
     False
     """
-    return isinstance(x, bsr_matrix) or isinstance(x, bsr_array)
+    return isinstance(x, bsr_matrix)
 
 
 class bsr_matrix(spmatrix, bsr_array):
diff --git a/scipy/sparse/_compressed.py b/scipy/sparse/_compressed.py
index eccd9b6a867d..49a17e5e8719 100644
--- a/scipy/sparse/_compressed.py
+++ b/scipy/sparse/_compressed.py
@@ -7,7 +7,7 @@
 import numpy as np
 from scipy._lib._util import _prune_array
 
-from ._base import _sparray, isspmatrix, SparseEfficiencyWarning
+from ._base import _sparray, issparse, SparseEfficiencyWarning
 from ._data import _data_matrix, _minmax_mixin
 from . import _sparsetools
 from ._sparsetools import (get_csr_submatrix, csr_sample_offsets, csr_todense,
@@ -25,7 +25,7 @@ class _cs_matrix(_data_matrix, _minmax_mixin, IndexMixin):
     def __init__(self, arg1, shape=None, dtype=None, copy=False):
         _data_matrix.__init__(self)
 
-        if isspmatrix(arg1):
+        if issparse(arg1):
             if arg1.format == self.format and copy:
                 arg1 = arg1.copy()
             else:
@@ -233,7 +233,7 @@ def __eq__(self, other):
         elif is_pydata_spmatrix(other):
             return NotImplemented
         # Sparse other.
-        elif isspmatrix(other):
+        elif issparse(other):
             warn("Comparing sparse matrices using == is inefficient, try using"
                  " != instead.", SparseEfficiencyWarning, stacklevel=3)
             # TODO sparse broadcasting
@@ -271,7 +271,7 @@ def __ne__(self, other):
         elif is_pydata_spmatrix(other):
             return NotImplemented
         # Sparse other.
-        elif isspmatrix(other):
+        elif issparse(other):
             # TODO sparse broadcasting
             if self.shape != other.shape:
                 return True
@@ -298,7 +298,7 @@ def _inequality(self, other, op, op_name, bad_scalar_msg):
         elif isdense(other):
             return op(self.todense(), other)
         # Sparse other.
-        elif isspmatrix(other):
+        elif issparse(other):
             # TODO sparse broadcasting
             if self.shape != other.shape:
                 raise ValueError("inconsistent shapes")
@@ -369,7 +369,7 @@ def multiply(self, other):
         if isscalarlike(other):
             return self._mul_scalar(other)
         # Sparse matrix or vector.
-        if isspmatrix(other):
+        if issparse(other):
             if self.shape == other.shape:
                 other = self.__class__(other)
                 return self._binopt(other, '_elmul_')
@@ -574,7 +574,7 @@ def _maximum_minimum(self, other, npop, op_name, dense_check):
                 return mat
         elif isdense(other):
             return npop(self.todense(), other)
-        elif isspmatrix(other):
+        elif issparse(other):
             return self._binopt(other, op_name)
         else:
             raise ValueError("Operands not compatible.")
diff --git a/scipy/sparse/_coo.py b/scipy/sparse/_coo.py
index 114ad3ec92a1..a9f146290a0e 100644
--- a/scipy/sparse/_coo.py
+++ b/scipy/sparse/_coo.py
@@ -10,7 +10,7 @@
 
 from ._matrix import spmatrix, _array_doc_to_matrix
 from ._sparsetools import coo_tocsr, coo_todense, coo_matvec
-from ._base import isspmatrix, SparseEfficiencyWarning, _sparray
+from ._base import issparse, SparseEfficiencyWarning, _sparray
 from ._data import _data_matrix, _minmax_mixin
 from ._sputils import (upcast, upcast_char, to_native, isshape, getdtype,
                        getdata, downcast_intp_index,
@@ -161,8 +161,8 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False):
                 self.data = getdata(obj, copy=copy, dtype=dtype)
                 self.has_canonical_format = False
         else:
-            if isspmatrix(arg1):
-                if isspmatrix_coo(arg1) and copy:
+            if issparse(arg1):
+                if arg1.format == self.format and copy:
                     self.row = arg1.row.copy()
                     self.col = arg1.col.copy()
                     self.data = arg1.data.copy()
@@ -595,7 +595,7 @@ def _mul_multivector(self, other):
 
 
 def isspmatrix_coo(x):
-    """Is x of coo_array type?
+    """Is `x` of coo_matrix type?
 
     Parameters
     ----------
@@ -605,19 +605,19 @@ def isspmatrix_coo(x):
     Returns
     -------
     bool
-        True if x is a coo matrix, False otherwise
+        True if `x` is a coo matrix, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import coo_array, isspmatrix_coo
-    >>> isspmatrix_coo(coo_array([[5]]))
+    >>> from scipy.sparse import coo_array, coo_matrix, csr_matrix, isspmatrix_coo
+    >>> isspmatrix_coo(coo_matrix([[5]]))
     True
-
-    >>> from scipy.sparse import coo_array, csr_matrix, isspmatrix_coo
+    >>> isspmatrix_coo(coo_array([[5]]))
+    False
     >>> isspmatrix_coo(csr_matrix([[5]]))
     False
     """
-    return isinstance(x, coo_matrix) or isinstance(x, coo_array)
+    return isinstance(x, coo_matrix)
 
 
 class coo_matrix(spmatrix, coo_array):
diff --git a/scipy/sparse/_csc.py b/scipy/sparse/_csc.py
index 3c739695e2fb..622583ad8c8d 100644
--- a/scipy/sparse/_csc.py
+++ b/scipy/sparse/_csc.py
@@ -235,7 +235,7 @@ def _swap(self, x):
 
 
 def isspmatrix_csc(x):
-    """Is x of csc_array type?
+    """Is `x` of csc_matrix type?
 
     Parameters
     ----------
@@ -245,19 +245,19 @@ def isspmatrix_csc(x):
     Returns
     -------
     bool
-        True if x is a csc matrix, False otherwise
+        True if `x` is a csc matrix, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import csc_array, isspmatrix_csc
-    >>> isspmatrix_csc(csc_array([[5]]))
+    >>> from scipy.sparse import csc_array, csc_matrix, coo_matrix, isspmatrix_csc
+    >>> isspmatrix_csc(csc_matrix([[5]]))
     True
-
-    >>> from scipy.sparse import csc_array, csr_matrix, isspmatrix_csc
+    >>> isspmatrix_csc(csc_array([[5]]))
+    False
     >>> isspmatrix_csc(csr_matrix([[5]]))
     False
     """
-    return isinstance(x, csc_matrix) or isinstance(x, csc_array)
+    return isinstance(x, csc_matrix)
 
 
 class csc_matrix(spmatrix, csc_array):
diff --git a/scipy/sparse/_csr.py b/scipy/sparse/_csr.py
index 53d48aaaae39..9d95ccb82941 100644
--- a/scipy/sparse/_csr.py
+++ b/scipy/sparse/_csr.py
@@ -332,7 +332,7 @@ def _get_arrayXslice(self, row, col):
 
 
 def isspmatrix_csr(x):
-    """Is x of csr_array type?
+    """Is `x` of csr_matrix type?
 
     Parameters
     ----------
@@ -342,19 +342,19 @@ def isspmatrix_csr(x):
     Returns
     -------
     bool
-        True if x is a csr matrix, False otherwise
+        True if `x` is a csr matrix, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import csr_array, isspmatrix_csr
-    >>> isspmatrix_csr(csr_array([[5]]))
+    >>> from scipy.sparse import csr_array, csr_matrix, coo_matrix, isspmatrix_csr
+    >>> isspmatrix_csr(csr_matrix([[5]]))
     True
-
-    >>> from scipy.sparse import csc_matrix, csr_array, isspmatrix_csc
-    >>> isspmatrix_csr(csc_matrix([[5]]))
+    >>> isspmatrix_csr(csr_array([[5]]))
+    False
+    >>> isspmatrix_csr(coo_matrix([[5]]))
     False
     """
-    return isinstance(x, csr_matrix) or isinstance(x, csr_array)
+    return isinstance(x, csr_matrix)
 
 
 class csr_matrix(spmatrix, csr_array):
diff --git a/scipy/sparse/_dia.py b/scipy/sparse/_dia.py
index 142eda34a335..046afcb29382 100644
--- a/scipy/sparse/_dia.py
+++ b/scipy/sparse/_dia.py
@@ -7,7 +7,7 @@
 import numpy as np
 
 from ._matrix import spmatrix, _array_doc_to_matrix
-from ._base import isspmatrix, _formats, _sparray
+from ._base import issparse, _formats, _sparray
 from ._data import _data_matrix
 from ._sputils import (isshape, upcast_char, getdtype, get_sum_dtype, validateaxis, check_shape)
 from ._sparsetools import dia_matvec
@@ -89,20 +89,21 @@ class dia_array(_data_matrix):
     def __init__(self, arg1, shape=None, dtype=None, copy=False):
         _data_matrix.__init__(self)
 
-        if isspmatrix_dia(arg1):
-            if copy:
-                arg1 = arg1.copy()
-            self.data = arg1.data
-            self.offsets = arg1.offsets
-            self._shape = check_shape(arg1.shape)
-        elif isspmatrix(arg1):
-            if isspmatrix_dia(arg1) and copy:
-                A = arg1.copy()
+        if issparse(arg1):
+            if arg1.format == "dia":
+                if copy:
+                    arg1 = arg1.copy()
+                self.data = arg1.data
+                self.offsets = arg1.offsets
+                self._shape = check_shape(arg1.shape)
             else:
-                A = arg1.todia()
-            self.data = A.data
-            self.offsets = A.offsets
-            self._shape = check_shape(A.shape)
+                if arg1.format == self.format and copy:
+                    A = arg1.copy()
+                else:
+                    A = arg1.todia()
+                self.data = A.data
+                self.offsets = A.offsets
+                self._shape = check_shape(A.shape)
         elif isinstance(arg1, tuple):
             if isshape(arg1):
                 # It's a tuple of matrix dimensions (M, N)
@@ -444,7 +445,7 @@ def resize(self, *shape):
 
 
 def isspmatrix_dia(x):
-    """Is x of dia_array type?
+    """Is `x` of dia_matrix type?
 
     Parameters
     ----------
@@ -454,19 +455,19 @@ def isspmatrix_dia(x):
     Returns
     -------
     bool
-        True if x is a dia matrix, False otherwise
+        True if `x` is a dia matrix, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import dia_array, isspmatrix_dia
-    >>> isspmatrix_dia(dia_array([[5]]))
+    >>> from scipy.sparse import dia_array, dia_matrix, coo_matrix, isspmatrix_dia
+    >>> isspmatrix_dia(dia_matrix([[5]]))
     True
-
-    >>> from scipy.sparse import dia_array, csr_matrix, isspmatrix_dia
-    >>> isspmatrix_dia(csr_matrix([[5]]))
+    >>> isspmatrix_dia(dia_array([[5]]))
+    False
+    >>> isspmatrix_dia(coo_matrix([[5]]))
     False
     """
-    return isinstance(x, dia_matrix) or isinstance(x, dia_array)
+    return isinstance(x, dia_matrix)
 
 
 class dia_matrix(spmatrix, dia_array):
diff --git a/scipy/sparse/_dok.py b/scipy/sparse/_dok.py
index 69f8a073e8ab..4449cf8aba8e 100644
--- a/scipy/sparse/_dok.py
+++ b/scipy/sparse/_dok.py
@@ -8,7 +8,7 @@
 import numpy as np
 
 from ._matrix import spmatrix, _array_doc_to_matrix
-from ._base import _sparray, isspmatrix
+from ._base import _sparray, issparse
 from ._index import IndexMixin
 from ._sputils import (isdense, getdtype, isshape, isintlike, isscalarlike,
                        upcast, upcast_scalar, check_shape)
@@ -80,8 +80,8 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False):
         if isinstance(arg1, tuple) and isshape(arg1):  # (M,N)
             M, N = arg1
             self._shape = check_shape((M, N))
-        elif isspmatrix(arg1):  # Sparse ctor
-            if isspmatrix_dok(arg1) and copy:
+        elif issparse(arg1):  # Sparse ctor
+            if arg1.format == self.format and copy:
                 arg1 = arg1.copy()
             else:
                 arg1 = arg1.todok()
@@ -246,20 +246,21 @@ def __add__(self, other):
                 if aij:
                     new[key] = aij
             # new.dtype.char = self.dtype.char
-        elif isspmatrix_dok(other):
-            if other.shape != self.shape:
-                raise ValueError("Matrix dimensions are not equal.")
-            # We could alternatively set the dimensions to the largest of
-            # the two matrices to be summed.  Would this be a good idea?
-            res_dtype = upcast(self.dtype, other.dtype)
-            new = self._dok_container(self.shape, dtype=res_dtype)
-            dict.update(new, self)
-            with np.errstate(over='ignore'):
-                dict.update(new,
-                           ((k, new[k] + other[k]) for k in other.keys()))
-        elif isspmatrix(other):
-            csc = self.tocsc()
-            new = csc + other
+        elif issparse(other):
+            if other.format == "dok":
+                if other.shape != self.shape:
+                    raise ValueError("Matrix dimensions are not equal.")
+                # We could alternatively set the dimensions to the largest of
+                # the two matrices to be summed.  Would this be a good idea?
+                res_dtype = upcast(self.dtype, other.dtype)
+                new = self._dok_container(self.shape, dtype=res_dtype)
+                dict.update(new, self)
+                with np.errstate(over='ignore'):
+                    dict.update(new,
+                               ((k, new[k] + other[k]) for k in other.keys()))
+            else:
+                csc = self.tocsc()
+                new = csc + other
         elif isdense(other):
             new = self.todense() + other
         else:
@@ -274,16 +275,17 @@ def __radd__(self, other):
                 aij = dict.get(self, (key), 0) + other
                 if aij:
                     new[key] = aij
-        elif isspmatrix_dok(other):
-            if other.shape != self.shape:
-                raise ValueError("Matrix dimensions are not equal.")
-            new = self._dok_container(self.shape, dtype=self.dtype)
-            dict.update(new, self)
-            dict.update(new,
-                       ((k, self[k] + other[k]) for k in other.keys()))
-        elif isspmatrix(other):
-            csc = self.tocsc()
-            new = csc + other
+        elif issparse(other):
+            if other.format == "dok":
+                if other.shape != self.shape:
+                    raise ValueError("Matrix dimensions are not equal.")
+                new = self._dok_container(self.shape, dtype=self.dtype)
+                dict.update(new, self)
+                dict.update(new,
+                           ((k, self[k] + other[k]) for k in other.keys()))
+            else:
+                csc = self.tocsc()
+                new = csc + other
         elif isdense(other):
             new = other + self.todense()
         else:
@@ -423,7 +425,7 @@ def resize(self, *shape):
 
 
 def isspmatrix_dok(x):
-    """Is x of dok_array type?
+    """Is `x` of dok_array type?
 
     Parameters
     ----------
@@ -433,19 +435,19 @@ def isspmatrix_dok(x):
     Returns
     -------
     bool
-        True if x is a dok matrix, False otherwise
+        True if `x` is a dok matrix, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import dok_array, isspmatrix_dok
-    >>> isspmatrix_dok(dok_array([[5]]))
+    >>> from scipy.sparse import dok_array, dok_matrix, coo_matrix, isspmatrix_dok
+    >>> isspmatrix_dok(dok_matrix([[5]]))
     True
-
-    >>> from scipy.sparse import dok_array, csr_matrix, isspmatrix_dok
-    >>> isspmatrix_dok(csr_matrix([[5]]))
+    >>> isspmatrix_dok(dok_array([[5]]))
+    False
+    >>> isspmatrix_dok(coo_matrix([[5]]))
     False
     """
-    return isinstance(x, dok_matrix) or isinstance(x, dok_array)
+    return isinstance(x, dok_matrix)
 
 
 class dok_matrix(spmatrix, dok_array):
diff --git a/scipy/sparse/_index.py b/scipy/sparse/_index.py
index dc2e30d56425..030d05d7c0db 100644
--- a/scipy/sparse/_index.py
+++ b/scipy/sparse/_index.py
@@ -116,8 +116,8 @@ def __setitem__(self, key, x):
         if i.shape != j.shape:
             raise IndexError('number of row and column indices differ')
 
-        from ._base import isspmatrix
-        if isspmatrix(x):
+        from ._base import issparse
+        if issparse(x):
             if i.ndim == 1:
                 # Inner indexing, so treat them like row vectors.
                 i = i[None]
@@ -266,7 +266,7 @@ def _unpack_index(index):
     Valid type for row/col is integer, slice, or array of integers.
     """
     # First, check if indexing with single boolean matrix.
-    from ._base import _sparray, isspmatrix
+    from ._base import _sparray, issparse
     if (isinstance(index, (_sparray, np.ndarray)) and
             index.ndim == 2 and index.dtype.kind == 'b'):
         return index.nonzero()
@@ -291,7 +291,7 @@ def _unpack_index(index):
         elif idx.ndim == 2:
             return idx.nonzero()
     # Next, check for validity and transform the index as needed.
-    if isspmatrix(row) or isspmatrix(col):
+    if issparse(row) or issparse(col):
         # Supporting sparse boolean indexing with both row and col does
         # not work because spmatrix.ndim is always 2.
         raise IndexError(
diff --git a/scipy/sparse/_lil.py b/scipy/sparse/_lil.py
index bf0f26a988b3..68a273345682 100644
--- a/scipy/sparse/_lil.py
+++ b/scipy/sparse/_lil.py
@@ -10,7 +10,7 @@
 import numpy as np
 
 from ._matrix import spmatrix, _array_doc_to_matrix
-from ._base import _sparray, isspmatrix
+from ._base import _sparray, issparse
 from ._index import IndexMixin, INT_TYPES, _broadcast_arrays
 from ._sputils import (getdtype, isshape, isscalarlike, upcast_scalar,
                        check_shape, check_reshape_kwargs)
@@ -86,8 +86,8 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False):
         self.dtype = getdtype(dtype, arg1, default=float)
 
         # First get the shape
-        if isspmatrix(arg1):
-            if isspmatrix_lil(arg1) and copy:
+        if issparse(arg1):
+            if arg1.format == "lil" and copy:
                 A = arg1.copy()
             else:
                 A = arg1.tolil()
@@ -323,7 +323,7 @@ def __setitem__(self, key, x):
             # Fast path for full-matrix sparse assignment.
             if (isinstance(row, slice) and isinstance(col, slice) and
                     row == slice(None) and col == slice(None) and
-                    isspmatrix(x) and x.shape == self.shape):
+                    issparse(x) and x.shape == self.shape):
                 x = self._lil_container(x, dtype=self.dtype)
                 self.rows = x.rows
                 self.data = x.data
@@ -522,7 +522,7 @@ def _prepare_index_for_memoryview(i, j, x=None):
 
 
 def isspmatrix_lil(x):
-    """Is x of lil_array type?
+    """Is `x` of lil_matrix type?
 
     Parameters
     ----------
@@ -532,19 +532,19 @@ def isspmatrix_lil(x):
     Returns
     -------
     bool
-        True if x is a lil matrix, False otherwise
+        True if `x` is a lil matrix, False otherwise
 
     Examples
     --------
-    >>> from scipy.sparse import lil_array, isspmatrix_lil
-    >>> isspmatrix_lil(lil_array([[5]]))
+    >>> from scipy.sparse import lil_array, lil_matrix, coo_matrix, isspmatrix_lil
+    >>> isspmatrix_lil(lil_matrix([[5]]))
     True
-
-    >>> from scipy.sparse import lil_array, csr_matrix, isspmatrix_lil
-    >>> isspmatrix_lil(csr_matrix([[5]]))
+    >>> isspmatrix_lil(lil_array([[5]]))
+    False
+    >>> isspmatrix_lil(coo_matrix([[5]]))
     False
     """
-    return isinstance(x, lil_matrix) or isinstance(x, lil_array)
+    return isinstance(x, lil_matrix)
 
 
 class lil_matrix(spmatrix, lil_array):
diff --git a/scipy/sparse/_spfuncs.py b/scipy/sparse/_spfuncs.py
index b6ec302724d5..8e9b0abcede6 100644
--- a/scipy/sparse/_spfuncs.py
+++ b/scipy/sparse/_spfuncs.py
@@ -3,8 +3,8 @@
 
 __all__ = ['count_blocks','estimate_blocksize']
 
-from ._csr import isspmatrix_csr, csr_array
-from ._csc import isspmatrix_csc
+from ._base import issparse
+from ._csr import csr_array
 from ._sparsetools import csr_count_blocks
 
 
@@ -14,7 +14,7 @@ def estimate_blocksize(A,efficiency=0.7):
     Returns a blocksize=(r,c) such that
         - A.nnz / A.tobsr( (r,c) ).nnz > efficiency
     """
-    if not (isspmatrix_csr(A) or isspmatrix_csc(A)):
+    if not (issparse(A) and A.format in ("csc", "csr")):
         A = csr_array(A)
 
     if A.nnz == 0:
@@ -67,10 +67,10 @@ def count_blocks(A,blocksize):
     if r < 1 or c < 1:
         raise ValueError('r and c must be positive')
 
-    if isspmatrix_csr(A):
-        M,N = A.shape
-        return csr_count_blocks(M,N,r,c,A.indptr,A.indices)
-    elif isspmatrix_csc(A):
-        return count_blocks(A.T,(c,r))
-    else:
-        return count_blocks(csr_array(A),blocksize)
+    if issparse(A):
+        if A.format == "csr":
+            M,N = A.shape
+            return csr_count_blocks(M,N,r,c,A.indptr,A.indices)
+        elif A.format == "csc":
+            return count_blocks(A.T,(c,r))
+    return count_blocks(csr_array(A),blocksize)
diff --git a/scipy/sparse/csgraph/_flow.pyx b/scipy/sparse/csgraph/_flow.pyx
index b3b24da416e8..d6dd88c365a9 100644
--- a/scipy/sparse/csgraph/_flow.pyx
+++ b/scipy/sparse/csgraph/_flow.pyx
@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from scipy.sparse import csr_matrix, isspmatrix_csr
+from scipy.sparse import csr_matrix, issparse
 
 cimport numpy as np
 
@@ -224,7 +224,7 @@ def maximum_flow(csgraph, source, sink, *, method='dinic'):
     modifying the capacities of the new graph appropriately.
 
     """
-    if not isspmatrix_csr(csgraph):
+    if not (issparse(csgraph) and csgraph.format == "csr"):
         raise TypeError("graph must be in CSR format")
     if not issubclass(csgraph.dtype.type, np.integer):
         raise ValueError("graph capacities must be integers")
diff --git a/scipy/sparse/csgraph/_laplacian.py b/scipy/sparse/csgraph/_laplacian.py
index 4a47ab665da8..3ac929239a5c 100644
--- a/scipy/sparse/csgraph/_laplacian.py
+++ b/scipy/sparse/csgraph/_laplacian.py
@@ -3,7 +3,7 @@
 """
 
 import numpy as np
-from scipy.sparse import isspmatrix
+from scipy.sparse import issparse
 from scipy.sparse.linalg import LinearOperator
 
 
@@ -340,12 +340,12 @@ def laplacian(
 
     if form == "array":
         create_lap = (
-            _laplacian_sparse if isspmatrix(csgraph) else _laplacian_dense
+            _laplacian_sparse if issparse(csgraph) else _laplacian_dense
         )
     else:
         create_lap = (
             _laplacian_sparse_flo
-            if isspmatrix(csgraph)
+            if issparse(csgraph)
             else _laplacian_dense_flo
         )
 
diff --git a/scipy/sparse/csgraph/_matching.pyx b/scipy/sparse/csgraph/_matching.pyx
index da16a3a856fd..1a4f582a74e1 100644
--- a/scipy/sparse/csgraph/_matching.pyx
+++ b/scipy/sparse/csgraph/_matching.pyx
@@ -6,7 +6,7 @@ cimport numpy as np
 from libc.math cimport INFINITY
 
 
-from scipy.sparse import isspmatrix_coo, isspmatrix_csc, isspmatrix_csr
+from scipy.sparse import issparse
 
 np.import_array()
 
@@ -130,10 +130,11 @@ def maximum_bipartite_matching(graph, perm_type='row'):
      [2 0 0 3]]
 
     """
-    if isspmatrix_csc(graph) or isspmatrix_coo(graph):
-        graph = graph.tocsr()
-    elif not isspmatrix_csr(graph):
+    if not issparse(graph):
+        raise TypeError("graph must be sparse")
+    if graph.format not in ("csr", "csc", "coo"):
         raise TypeError("graph must be in CSC, CSR, or COO format.")
+    graph = graph.tocsr()
     i, j = graph.shape
     x, y = _hopcroft_karp(graph.indices, graph.indptr, i, j)
     return np.asarray(x if perm_type == 'column' else y)
@@ -443,9 +444,9 @@ def min_weight_full_bipartite_matching(biadjacency_matrix, maximize=False):
     28.0
 
     """
-    if not isspmatrix_csr(biadjacency_matrix)\
-            and not isspmatrix_csc(biadjacency_matrix)\
-            and not isspmatrix_coo(biadjacency_matrix):
+    if not issparse(biadjacency_matrix):
+        raise TypeError("graph must be sparse")
+    if biadjacency_matrix.format not in ("csr", "csc", "coo"):
         raise TypeError("graph must be in CSC, CSR, or COO format.")
 
     if not (np.issubdtype(biadjacency_matrix.dtype, np.number) or
@@ -473,7 +474,7 @@ def min_weight_full_bipartite_matching(biadjacency_matrix, maximize=False):
     # The algorithm expects more columns than rows in the graph.
     if j < i:
         biadjacency_matrix_t = biadjacency_matrix.T
-        if not isspmatrix_csr(biadjacency_matrix_t):
+        if biadjacency_matrix_t.format != "csr":
             biadjacency_matrix_t = biadjacency_matrix_t.tocsr()
         b = np.asarray(_lapjvsp(biadjacency_matrix_t.indptr,
                                 biadjacency_matrix_t.indices,
@@ -482,7 +483,7 @@ def min_weight_full_bipartite_matching(biadjacency_matrix, maximize=False):
         indices = np.argsort(b)
         return (b[indices], a[indices])
     else:
-        if not isspmatrix_csr(biadjacency_matrix):
+        if biadjacency_matrix.format != "csr":
             biadjacency_matrix = biadjacency_matrix.tocsr()
         b = np.asarray(_lapjvsp(biadjacency_matrix.indptr,
                                 biadjacency_matrix.indices,
diff --git a/scipy/sparse/csgraph/_reordering.pyx b/scipy/sparse/csgraph/_reordering.pyx
index 304777c4a0f6..f45dd4344505 100644
--- a/scipy/sparse/csgraph/_reordering.pyx
+++ b/scipy/sparse/csgraph/_reordering.pyx
@@ -5,9 +5,7 @@
 import numpy as np
 cimport numpy as np
 from warnings import warn
-from scipy.sparse import (csr_matrix, isspmatrix, isspmatrix_coo,
-                          isspmatrix_csc, isspmatrix_csr,
-                          SparseEfficiencyWarning)
+from scipy.sparse import csr_matrix, issparse, SparseEfficiencyWarning
 from . import maximum_bipartite_matching
 
 np.import_array()
@@ -70,7 +68,9 @@ def reverse_cuthill_mckee(graph, symmetric_mode=False):
     array([3, 2, 1, 0], dtype=int32)
     
     """
-    if not (isspmatrix_csc(graph) or isspmatrix_csr(graph)):
+    if not issparse(graph):
+        raise TypeError("Input graph must be sparse")
+    if graph.format not in ("csc", "csr"):
         raise TypeError('Input must be in CSC or CSR sparse matrix format.')
     nrows = graph.shape[0]
     if not symmetric_mode:
@@ -226,10 +226,10 @@ def structural_rank(graph):
     4
 
     """
-    if not isspmatrix:
+    if not issparse:
         raise TypeError('Input must be a sparse matrix')
-    if not isspmatrix_csr(graph):
-        if not (isspmatrix_csc(graph) or isspmatrix_coo(graph)):
+    if graph.format != "csr":
+        if graph.format not in ("csc", "coo"):
             warn('Input matrix should be in CSC, CSR, or COO matrix format',
                     SparseEfficiencyWarning)
         graph = csr_matrix(graph)
diff --git a/scipy/sparse/csgraph/_shortest_path.pyx b/scipy/sparse/csgraph/_shortest_path.pyx
index b356521888ab..01502ee19d0a 100644
--- a/scipy/sparse/csgraph/_shortest_path.pyx
+++ b/scipy/sparse/csgraph/_shortest_path.pyx
@@ -14,7 +14,7 @@ import warnings
 import numpy as np
 cimport numpy as np
 
-from scipy.sparse import csr_matrix, isspmatrix
+from scipy.sparse import csr_matrix, issparse
 from scipy.sparse.csgraph._validation import validate_graph
 
 cimport cython
@@ -157,14 +157,14 @@ def shortest_path(csgraph, method='auto',
                    copy_if_dense=(not overwrite),
                    copy_if_sparse=(not overwrite))
 
-    cdef bint issparse
+    cdef bint is_sparse
     cdef ssize_t N      # XXX cdef ssize_t Nk fails in Python 3 (?)
 
     if method == 'auto':
         # guess fastest method based on number of nodes and edges
         N = csgraph.shape[0]
-        issparse = isspmatrix(csgraph)
-        if issparse:
+        is_sparse = issparse(csgraph)
+        if is_sparse:
             Nk = csgraph.nnz
             if csgraph.format in ('csr', 'csc', 'coo'):
                 edges = csgraph.data
@@ -301,7 +301,7 @@ def floyd_warshall(csgraph, directed=True,
     dist_matrix = validate_graph(csgraph, directed, DTYPE,
                                  csr_output=False,
                                  copy_if_dense=not overwrite)
-    if not isspmatrix(csgraph):
+    if not issparse(csgraph):
         # for dense array input, zero entries represent non-edge
         dist_matrix[dist_matrix == 0] = INFINITY
 
diff --git a/scipy/sparse/csgraph/_tools.pyx b/scipy/sparse/csgraph/_tools.pyx
index dde34fe78de0..bfe27971022f 100644
--- a/scipy/sparse/csgraph/_tools.pyx
+++ b/scipy/sparse/csgraph/_tools.pyx
@@ -8,8 +8,7 @@ Tools and utilities for working with compressed sparse graphs
 import numpy as np
 cimport numpy as np
 
-from scipy.sparse import csr_matrix, isspmatrix,\
-    isspmatrix_csr, isspmatrix_csc, isspmatrix_lil
+from scipy.sparse import csr_matrix, issparse
 
 np.import_array()
 
@@ -315,10 +314,11 @@ def csgraph_to_dense(csgraph, null_value=0):
     """
     # Allow only csr, lil and csc matrices: other formats when converted to csr
     # combine duplicated edges: we don't want this to happen in the background.
-    if isspmatrix_csc(csgraph) or isspmatrix_lil(csgraph):
-        csgraph = csgraph.tocsr()
-    elif not isspmatrix_csr(csgraph):
+    if not issparse(csgraph):
+        raise ValueError("csgraph must be sparse")
+    if csgraph.format not in ("lil", "csc", "csr"):
         raise ValueError("csgraph must be lil, csr, or csc format")
+    csgraph = csgraph.tocsr()
 
     N = csgraph.shape[0]
     if csgraph.shape[1] != N:
@@ -483,13 +483,13 @@ def reconstruct_path(csgraph, predecessors, directed=True):
     # Fix issue #4018:
     # If `pind` and `indices` are empty arrays, `data` is a sparse matrix
     # (it is a numpy.matrix otherwise); handle this case separately.
-    if isspmatrix(data):
+    if issparse(data):
         data = data.todense()
     data = data.getA1()
 
     if not directed:
         data2 = csgraph[indices, pind]
-        if isspmatrix(data2):
+        if issparse(data2):
             data2 = data2.todense()
         data2 = data2.getA1()
         data[data == 0] = np.inf
diff --git a/scipy/sparse/csgraph/_validation.py b/scipy/sparse/csgraph/_validation.py
index d93060c9b7b7..2f111d4b2dcd 100644
--- a/scipy/sparse/csgraph/_validation.py
+++ b/scipy/sparse/csgraph/_validation.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy.sparse import csr_matrix, isspmatrix, isspmatrix_csc
+from scipy.sparse import csr_matrix, issparse
 from ._tools import csgraph_to_dense, csgraph_from_dense,\
     csgraph_masked_from_dense, csgraph_from_masked
 
@@ -17,10 +17,10 @@ def validate_graph(csgraph, directed, dtype=DTYPE,
 
     # if undirected and csc storage, then transposing in-place
     # is quicker than later converting to csr.
-    if (not directed) and isspmatrix_csc(csgraph):
+    if (not directed) and issparse(csgraph) and csgraph.format == "csc":
         csgraph = csgraph.T
 
-    if isspmatrix(csgraph):
+    if issparse(csgraph):
         if csr_output:
             csgraph = csr_matrix(csgraph, dtype=DTYPE, copy=copy_if_sparse)
         else:
diff --git a/scipy/sparse/linalg/_dsolve/linsolve.py b/scipy/sparse/linalg/_dsolve/linsolve.py
index 939cbd04a277..0c408f737266 100644
--- a/scipy/sparse/linalg/_dsolve/linsolve.py
+++ b/scipy/sparse/linalg/_dsolve/linsolve.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 from numpy import asarray
-from scipy.sparse import (isspmatrix_csc, isspmatrix_csr, isspmatrix,
+from scipy.sparse import (issparse,
                           SparseEfficiencyWarning, csc_matrix, csr_matrix)
 from scipy.sparse._sputils import is_pydata_spmatrix
 from scipy.linalg import LinAlgError
@@ -209,13 +209,13 @@ def spsolve(A, b, permc_spec=None, use_umfpack=True):
     if is_pydata_spmatrix(A):
         A = A.to_scipy_sparse().tocsc()
 
-    if not (isspmatrix_csc(A) or isspmatrix_csr(A)):
+    if not (issparse(A) and A.format in ("csc", "csr")):
         A = csc_matrix(A)
         warn('spsolve requires A be CSC or CSR matrix format',
                 SparseEfficiencyWarning)
 
     # b is a vector only if b have shape (n,) or (n, 1)
-    b_is_sparse = isspmatrix(b) or is_pydata_spmatrix(b)
+    b_is_sparse = issparse(b) or is_pydata_spmatrix(b)
     if not b_is_sparse:
         b = asarray(b)
     b_is_vector = ((b.ndim == 1) or (b.ndim == 2 and b.shape[1] == 1))
@@ -264,7 +264,7 @@ def spsolve(A, b, permc_spec=None, use_umfpack=True):
             b_is_sparse = False
 
         if not b_is_sparse:
-            if isspmatrix_csc(A):
+            if A.format == "csc":
                 flag = 1  # CSC format
             else:
                 flag = 0  # CSR format
@@ -281,7 +281,7 @@ def spsolve(A, b, permc_spec=None, use_umfpack=True):
             # b is sparse
             Afactsolve = factorized(A)
 
-            if not (isspmatrix_csc(b) or is_pydata_spmatrix(b)):
+            if not (b.format == "csc" or is_pydata_spmatrix(b)):
                 warn('spsolve is more efficient when sparse b '
                      'is in the CSC matrix format', SparseEfficiencyWarning)
                 b = csc_matrix(b)
@@ -390,7 +390,7 @@ def csc_construct_func(*a, cls=type(A)):
     else:
         csc_construct_func = csc_matrix
 
-    if not isspmatrix_csc(A):
+    if not (issparse(A) and A.format == "csc"):
         A = csc_matrix(A)
         warn('splu converted its input to CSC format', SparseEfficiencyWarning)
 
@@ -482,7 +482,7 @@ def csc_construct_func(*a, cls=type(A)):
     else:
         csc_construct_func = csc_matrix
 
-    if not isspmatrix_csc(A):
+    if not (issparse(A) and A.format == "csc"):
         A = csc_matrix(A)
         warn('spilu converted its input to CSC format',
              SparseEfficiencyWarning)
@@ -547,7 +547,7 @@ def factorized(A):
         if noScikit:
             raise RuntimeError('Scikits.umfpack not installed.')
 
-        if not isspmatrix_csc(A):
+        if not (issparse(A) and A.format == "csc"):
             A = csc_matrix(A)
             warn('splu converted its input to CSC format',
                  SparseEfficiencyWarning)
@@ -638,7 +638,7 @@ def spsolve_triangular(A, b, lower=True, overwrite_A=False, overwrite_b=False,
         A = A.to_scipy_sparse().tocsr()
 
     # Check the input for correct type and format.
-    if not isspmatrix_csr(A):
+    if not (issparse(A) and A.format == "csr"):
         warn('CSR matrix format is required. Converting to CSR matrix.',
              SparseEfficiencyWarning)
         A = csr_matrix(A)
diff --git a/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py b/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
index 64b66ddae268..4740685a3dd7 100644
--- a/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
+++ b/scipy/sparse/linalg/_dsolve/tests/test_linsolve.py
@@ -14,7 +14,7 @@
 import scipy.linalg
 from scipy.linalg import norm, inv
 from scipy.sparse import (spdiags, SparseEfficiencyWarning, csc_matrix,
-        csr_matrix, identity, isspmatrix, dok_matrix, lil_matrix, bsr_matrix)
+        csr_matrix, identity, issparse, dok_matrix, lil_matrix, bsr_matrix)
 from scipy.sparse.linalg import SuperLU
 from scipy.sparse.linalg._dsolve import (spsolve, use_solver, splu, spilu,
         MatrixRankWarning, _superlu, spsolve_triangular, factorized)
@@ -35,7 +35,7 @@
     has_umfpack = False
 
 def toarray(a):
-    if isspmatrix(a):
+    if issparse(a):
         return a.toarray()
     else:
         return a
@@ -328,9 +328,9 @@ def test_shape_compatibility(self):
                 assert_array_almost_equal(toarray(x2), x, err_msg=repr((b, spmattype, 2)))
 
                 # dense vs. sparse output  ("vectors" are always dense)
-                if isspmatrix(b) and x.ndim > 1:
-                    assert_(isspmatrix(x1), repr((b, spmattype, 1)))
-                    assert_(isspmatrix(x2), repr((b, spmattype, 2)))
+                if issparse(b) and x.ndim > 1:
+                    assert_(issparse(x1), repr((b, spmattype, 1)))
+                    assert_(issparse(x2), repr((b, spmattype, 2)))
                 else:
                     assert_(isinstance(x1, np.ndarray), repr((b, spmattype, 1)))
                     assert_(isinstance(x2, np.ndarray), repr((b, spmattype, 2)))
diff --git a/scipy/sparse/linalg/_eigen/arpack/arpack.py b/scipy/sparse/linalg/_eigen/arpack/arpack.py
index 1221d50daba9..8ab5ed34dccf 100644
--- a/scipy/sparse/linalg/_eigen/arpack/arpack.py
+++ b/scipy/sparse/linalg/_eigen/arpack/arpack.py
@@ -38,7 +38,7 @@
 import numpy as np
 import warnings
 from scipy.sparse.linalg._interface import aslinearoperator, LinearOperator
-from scipy.sparse import eye, issparse, isspmatrix, isspmatrix_csr
+from scipy.sparse import eye, issparse
 from scipy.linalg import eig, eigh, lu_factor, lu_solve
 from scipy.sparse._sputils import isdense, is_pydata_spmatrix
 from scipy.sparse.linalg import gmres, splu
@@ -1038,7 +1038,7 @@ def dtype(self):
 
 def _fast_spmatrix_to_csc(A, hermitian=False):
     """Convert sparse matrix to CSC (by transposing, if possible)"""
-    if (isspmatrix_csr(A) and hermitian
+    if (A.format == "csr" and hermitian
             and not np.issubdtype(A.dtype, np.complexfloating)):
         return A.T
     elif is_pydata_spmatrix(A):
@@ -1051,7 +1051,7 @@ def _fast_spmatrix_to_csc(A, hermitian=False):
 def get_inv_matvec(M, hermitian=False, tol=0):
     if isdense(M):
         return LuInv(M).matvec
-    elif isspmatrix(M) or is_pydata_spmatrix(M):
+    elif issparse(M) or is_pydata_spmatrix(M):
         M = _fast_spmatrix_to_csc(M, hermitian=hermitian)
         return SpLuInv(M).matvec
     else:
@@ -1072,7 +1072,7 @@ def get_OPinv_matvec(A, M, sigma, hermitian=False, tol=0):
                 A = A + 0j
             A.flat[::A.shape[1] + 1] -= sigma
             return LuInv(A).matvec
-        elif isspmatrix(A) or is_pydata_spmatrix(A):
+        elif issparse(A) or is_pydata_spmatrix(A):
             A = A - sigma * eye(A.shape[0])
             A = _fast_spmatrix_to_csc(A, hermitian=hermitian)
             return SpLuInv(A).matvec
@@ -1080,8 +1080,8 @@ def get_OPinv_matvec(A, M, sigma, hermitian=False, tol=0):
             return IterOpInv(_aslinearoperator_with_dtype(A),
                              M, sigma, tol=tol).matvec
     else:
-        if ((not isdense(A) and not isspmatrix(A) and not is_pydata_spmatrix(A)) or
-                (not isdense(M) and not isspmatrix(M) and not is_pydata_spmatrix(A))):
+        if ((not isdense(A) and not issparse(A) and not is_pydata_spmatrix(A)) or
+                (not isdense(M) and not issparse(M) and not is_pydata_spmatrix(A))):
             return IterOpInv(_aslinearoperator_with_dtype(A),
                              _aslinearoperator_with_dtype(M),
                              sigma, tol=tol).matvec
diff --git a/scipy/sparse/linalg/_eigen/lobpcg/lobpcg.py b/scipy/sparse/linalg/_eigen/lobpcg/lobpcg.py
index 8e032f71ddcd..a46ce5ceb412 100644
--- a/scipy/sparse/linalg/_eigen/lobpcg/lobpcg.py
+++ b/scipy/sparse/linalg/_eigen/lobpcg/lobpcg.py
@@ -22,7 +22,7 @@
 from scipy.linalg import (inv, eigh, cho_factor, cho_solve,
                           cholesky, LinAlgError)
 from scipy.sparse.linalg import LinearOperator
-from scipy.sparse import isspmatrix
+from scipy.sparse import issparse
 from numpy import block as bmat
 
 __all__ = ["lobpcg"]
@@ -493,7 +493,7 @@ def lobpcg(
                         f"The shape {A.shape} of the primary matrix\n"
                         f"defined by a callable object is wrong.\n"
                     )
-            elif isspmatrix(A):
+            elif issparse(A):
                 A = A.toarray()
             else:
                 A = np.asarray(A)
@@ -513,7 +513,7 @@ def lobpcg(
                             f"The shape {B.shape} of the secondary matrix\n"
                             f"defined by a callable object is wrong.\n"
                         )
-                elif isspmatrix(B):
+                elif issparse(B):
                     B = B.toarray()
                 else:
                     B = np.asarray(B)
diff --git a/scipy/sparse/linalg/_eigen/tests/test_svds.py b/scipy/sparse/linalg/_eigen/tests/test_svds.py
index 7782d1b61b8d..3b6bfead16cd 100644
--- a/scipy/sparse/linalg/_eigen/tests/test_svds.py
+++ b/scipy/sparse/linalg/_eigen/tests/test_svds.py
@@ -7,7 +7,7 @@
 import pytest
 
 from scipy.linalg import svd, null_space
-from scipy.sparse import csc_matrix, isspmatrix, spdiags, random
+from scipy.sparse import csc_matrix, issparse, spdiags, random
 from scipy.sparse.linalg import LinearOperator, aslinearoperator
 if os.environ.get("SCIPY_USE_PROPACK"):
     has_propack = True
@@ -23,7 +23,7 @@
 def sorted_svd(m, k, which='LM'):
     # Compute svd of a dense matrix m, and return singular vectors/values
     # sorted.
-    if isspmatrix(m):
+    if issparse(m):
         m = m.toarray()
     u, s, vh = svd(m)
     if which == 'LM':
diff --git a/scipy/sparse/linalg/_expm_multiply.py b/scipy/sparse/linalg/_expm_multiply.py
index c71d1f50fc11..721c2ef34205 100644
--- a/scipy/sparse/linalg/_expm_multiply.py
+++ b/scipy/sparse/linalg/_expm_multiply.py
@@ -16,7 +16,7 @@
 
 def _exact_inf_norm(A):
     # A compatibility function which should eventually disappear.
-    if scipy.sparse.isspmatrix(A):
+    if scipy.sparse.issparse(A):
         return max(abs(A).sum(axis=1).flat)
     elif is_pydata_spmatrix(A):
         return max(abs(A).sum(axis=1))
@@ -26,7 +26,7 @@ def _exact_inf_norm(A):
 
 def _exact_1_norm(A):
     # A compatibility function which should eventually disappear.
-    if scipy.sparse.isspmatrix(A):
+    if scipy.sparse.issparse(A):
         return max(abs(A).sum(axis=0).flat)
     elif is_pydata_spmatrix(A):
         return max(abs(A).sum(axis=0))
@@ -91,7 +91,7 @@ def traceest(A, m3, seed=None):
 
 def _ident_like(A):
     # A compatibility function which should eventually disappear.
-    if scipy.sparse.isspmatrix(A):
+    if scipy.sparse.issparse(A):
         return scipy.sparse._construct.eye(A.shape[0], A.shape[1],
                 dtype=A.dtype, format=A.format)
     elif is_pydata_spmatrix(A):
diff --git a/scipy/sparse/linalg/_interface.py b/scipy/sparse/linalg/_interface.py
index b143eab00a8f..71e267bb7a87 100644
--- a/scipy/sparse/linalg/_interface.py
+++ b/scipy/sparse/linalg/_interface.py
@@ -44,7 +44,7 @@
 
 import numpy as np
 
-from scipy.sparse import isspmatrix
+from scipy.sparse import issparse
 from scipy.sparse._sputils import isshape, isintlike, asmatrix, is_pydata_spmatrix
 
 __all__ = ['LinearOperator', 'aslinearoperator']
@@ -325,7 +325,7 @@ def matmat(self, X):
         _matmat method to ensure that y has the correct type.
 
         """
-        if not (isspmatrix(X) or is_pydata_spmatrix(X)):
+        if not (issparse(X) or is_pydata_spmatrix(X)):
             X = np.asanyarray(X)
 
         if X.ndim != 2:
@@ -337,7 +337,7 @@ def matmat(self, X):
         try:
             Y = self._matmat(X)
         except Exception as e:
-            if isspmatrix(X) or is_pydata_spmatrix(X):
+            if issparse(X) or is_pydata_spmatrix(X):
                 raise TypeError(
                     "Unable to multiply a LinearOperator with a sparse matrix."
                     " Wrap the matrix in aslinearoperator first."
@@ -371,7 +371,7 @@ def rmatmat(self, X):
         This rmatmat wraps the user-specified rmatmat routine.
 
         """
-        if not (isspmatrix(X) or is_pydata_spmatrix(X)):
+        if not (issparse(X) or is_pydata_spmatrix(X)):
             X = np.asanyarray(X)
 
         if X.ndim != 2:
@@ -384,7 +384,7 @@ def rmatmat(self, X):
         try:
             Y = self._rmatmat(X)
         except Exception as e:
-            if isspmatrix(X) or is_pydata_spmatrix(X):
+            if issparse(X) or is_pydata_spmatrix(X):
                 raise TypeError(
                     "Unable to multiply a LinearOperator with a sparse matrix."
                     " Wrap the matrix in aslinearoperator() first."
@@ -434,7 +434,7 @@ def dot(self, x):
         elif np.isscalar(x):
             return _ScaledLinearOperator(self, x)
         else:
-            if not isspmatrix(x) and not is_pydata_spmatrix(x):
+            if not issparse(x) and not is_pydata_spmatrix(x):
                 # Sparse matrices shouldn't be converted to numpy arrays.
                 x = np.asarray(x)
 
@@ -487,7 +487,7 @@ def _rdot(self, x):
         elif np.isscalar(x):
             return _ScaledLinearOperator(self, x)
         else:
-            if not isspmatrix(x) and not is_pydata_spmatrix(x):
+            if not issparse(x) and not is_pydata_spmatrix(x):
                 # Sparse matrices shouldn't be converted to numpy arrays.
                 x = np.asarray(x)
 
@@ -872,7 +872,7 @@ def aslinearoperator(A):
         A = np.atleast_2d(np.asarray(A))
         return MatrixLinearOperator(A)
 
-    elif isspmatrix(A) or is_pydata_spmatrix(A):
+    elif issparse(A) or is_pydata_spmatrix(A):
         return MatrixLinearOperator(A)
 
     else:
diff --git a/scipy/sparse/linalg/_matfuncs.py b/scipy/sparse/linalg/_matfuncs.py
index 1f99cf7f60dc..bad9d75d0f9b 100644
--- a/scipy/sparse/linalg/_matfuncs.py
+++ b/scipy/sparse/linalg/_matfuncs.py
@@ -13,7 +13,7 @@
 import numpy as np
 from scipy.linalg._basic import solve, solve_triangular
 
-from scipy.sparse._base import isspmatrix
+from scipy.sparse._base import issparse
 from scipy.sparse.linalg import spsolve
 from scipy.sparse._sputils import is_pydata_spmatrix
 
@@ -67,7 +67,7 @@ def inv(A):
 
     """
     # Check input
-    if not (scipy.sparse.isspmatrix(A) or is_pydata_spmatrix(A)):
+    if not (scipy.sparse.issparse(A) or is_pydata_spmatrix(A)):
         raise TypeError('Input must be a sparse matrix')
 
     # Use sparse direct solver to solve "AX = I" accurately
@@ -111,7 +111,7 @@ def _onenorm_matrix_power_nnm(A, p):
 
 def _is_upper_triangular(A):
     # This function could possibly be of wider interest.
-    if isspmatrix(A):
+    if issparse(A):
         lower_part = scipy.sparse.tril(A, -1)
         # Check structural upper triangularity,
         # then coincidental upper triangularity if needed.
@@ -152,7 +152,7 @@ def _smart_matrix_product(A, B, alpha=None, structure=None):
         raise ValueError('expected B to be a rectangular matrix')
     f = None
     if structure == UPPER_TRIANGULAR:
-        if (not isspmatrix(A) and not isspmatrix(B)
+        if (not issparse(A) and not issparse(B)
                 and not is_pydata_spmatrix(A) and not is_pydata_spmatrix(B)):
             f, = scipy.linalg.get_blas_funcs(('trmm',), (A, B))
     if f is not None:
@@ -604,7 +604,7 @@ def _expm(A, use_exact_onenorm):
     # carefully handling sparse scenario
     if A.shape == (0, 0):
         out = np.zeros([0, 0], dtype=A.dtype)
-        if isspmatrix(A) or is_pydata_spmatrix(A):
+        if issparse(A) or is_pydata_spmatrix(A):
             return A.__class__(out)
         return out
 
@@ -614,13 +614,13 @@ def _expm(A, use_exact_onenorm):
 
         # Avoid indiscriminate casting to ndarray to
         # allow for sparse or other strange arrays
-        if isspmatrix(A) or is_pydata_spmatrix(A):
+        if issparse(A) or is_pydata_spmatrix(A):
             return A.__class__(out)
 
         return np.array(out)
 
     # Ensure input is of float type, to avoid integer overflows etc.
-    if ((isinstance(A, np.ndarray) or isspmatrix(A) or is_pydata_spmatrix(A))
+    if ((isinstance(A, np.ndarray) or issparse(A) or is_pydata_spmatrix(A))
             and not np.issubdtype(A.dtype, np.inexact)):
         A = A.astype(float)
 
@@ -702,7 +702,7 @@ def _solve_P_Q(U, V, structure=None):
     """
     P = U + V
     Q = -U + V
-    if isspmatrix(U) or is_pydata_spmatrix(U):
+    if issparse(U) or is_pydata_spmatrix(U):
         return spsolve(Q, P)
     elif structure is None:
         return solve(Q, P)
diff --git a/scipy/sparse/tests/test_array_api.py b/scipy/sparse/tests/test_array_api.py
index 7d24f34009b0..acd6aabde928 100644
--- a/scipy/sparse/tests/test_array_api.py
+++ b/scipy/sparse/tests/test_array_api.py
@@ -488,3 +488,60 @@ def test_format_property():
         assert M._format == fmt
         with pytest.raises(AttributeError):
             M.format = "qqq"
+
+
+def test_issparse():
+    m = scipy.sparse.eye(3)
+    a = scipy.sparse.csr_array(m)
+    assert not m._is_array
+    assert a._is_array
+
+    # Both sparse arrays and sparse matrices should be sparse
+    assert scipy.sparse.issparse(a)
+    assert scipy.sparse.issparse(m)
+
+    # ndarray and array_likes are not sparse
+    assert not scipy.sparse.issparse(a.todense())
+    assert not scipy.sparse.issparse(m.todense())
+
+
+def test_isspmatrix():
+    m = scipy.sparse.eye(3)
+    a = scipy.sparse.csr_array(m)
+    assert not m._is_array
+    assert a._is_array
+
+    # Should only be true for sparse matrices, not sparse arrays
+    assert not scipy.sparse.isspmatrix(a)
+    assert scipy.sparse.isspmatrix(m)
+
+    # ndarray and array_likes are not sparse
+    assert not scipy.sparse.isspmatrix(a.todense())
+    assert not scipy.sparse.isspmatrix(m.todense())
+
+
+@pytest.mark.parametrize(
+    ("fmt", "fn"),
+    (
+        ("bsr", scipy.sparse.isspmatrix_bsr),
+        ("coo", scipy.sparse.isspmatrix_coo),
+        ("csc", scipy.sparse.isspmatrix_csc),
+        ("csr", scipy.sparse.isspmatrix_csr),
+        ("dia", scipy.sparse.isspmatrix_dia),
+        ("dok", scipy.sparse.isspmatrix_dok),
+        ("lil", scipy.sparse.isspmatrix_lil),
+    ),
+)
+def test_isspmatrix_format(fmt, fn):
+    m = scipy.sparse.eye(3, format=fmt)
+    a = scipy.sparse.csr_array(m).asformat(fmt)
+    assert not m._is_array
+    assert a._is_array
+
+    # Should only be true for sparse matrices, not sparse arrays
+    assert not fn(a)
+    assert fn(m)
+
+    # ndarray and array_likes are not sparse
+    assert not fn(a.todense())
+    assert not fn(m.todense())
diff --git a/scipy/sparse/tests/test_base.py b/scipy/sparse/tests/test_base.py
index 102fde58badd..0ee61dd02e0e 100644
--- a/scipy/sparse/tests/test_base.py
+++ b/scipy/sparse/tests/test_base.py
@@ -33,7 +33,7 @@ class for generic tests" section.
 import scipy.sparse as sparse
 from scipy.sparse import (csc_matrix, csr_matrix, dok_matrix,
         coo_matrix, lil_matrix, dia_matrix, bsr_matrix,
-        eye, isspmatrix, SparseEfficiencyWarning)
+        eye, issparse, SparseEfficiencyWarning)
 from scipy.sparse._sputils import (supported_dtypes, isscalarlike,
                                    get_index_dtype, asmatrix, matrix)
 from scipy.sparse.linalg import splu, expm, inv
@@ -1528,7 +1528,7 @@ def test_elementwise_multiply_broadcast(self):
                     assert_raises(ValueError, i.multiply, j)
                     continue
                 sp_mult = i.multiply(j)
-                if isspmatrix(sp_mult):
+                if issparse(sp_mult):
                     assert_almost_equal(sp_mult.toarray(), dense_mult)
                 else:
                     assert_almost_equal(sp_mult, dense_mult)
@@ -1596,7 +1596,7 @@ def test_small_multiplication(self):
         # test that A*x works for x with shape () (1,) (1,1) and (1,0)
         A = self.spcreator([[1],[2],[3]])
 
-        assert_(isspmatrix(A * array(1)))
+        assert_(issparse(A * array(1)))
         assert_equal((A * array(1)).toarray(), [[1], [2], [3]])
 
         assert_equal(A @ array([1]), array([1, 2, 3]))
@@ -2894,7 +2894,7 @@ def test_fancy_indexing_randomized(self):
         S = self.spcreator(D)
 
         SIJ = S[I,J]
-        if isspmatrix(SIJ):
+        if issparse(SIJ):
             SIJ = SIJ.toarray()
         assert_equal(SIJ, D[I,J])
 
@@ -3142,7 +3142,7 @@ def test_fancy_indexing_ndarray(self):
             S = self.spcreator(D)
 
             SIJ = S[I,J]
-            if isspmatrix(SIJ):
+            if issparse(SIJ):
                 SIJ = SIJ.toarray()
             assert_equal(SIJ, D[I,J])
 
@@ -3800,7 +3800,7 @@ def test_fancy_indexing_broadcast(self):
         S = self.spcreator(D)
 
         SIJ = S[I,J]
-        if isspmatrix(SIJ):
+        if issparse(SIJ):
             SIJ = SIJ.toarray()
         assert_equal(SIJ, D[I,J])
 
@@ -4012,7 +4012,7 @@ def test_fancy_indexing_broadcast(self):
         S = self.spcreator(D)
 
         SIJ = S[I,J]
-        if isspmatrix(SIJ):
+        if issparse(SIJ):
             SIJ = SIJ.toarray()
         assert_equal(SIJ, D[I,J])
 


Code review comments:
diff_hunk: @@ -1452,23 +1454,46 @@ def issparse(x):
     Returns
     -------
     bool
-        True if x is a sparse array, False otherwise
-
-    Notes
-    -----
-    issparse and isspmatrix are aliases for the same function.
+        True if `x` is a sparse array, False otherwise
path: scipy/sparse/_base.py
text: @user1:
Must sparse matrices be mentioned as well?
```suggestion
        True if `x` is a sparse array or a sparse matrix, False otherwise
```

diff_hunk: @@ -372,14 +372,12 @@ def _mul_sparse_matrix(self, other):
         R,n = self.blocksize
 
         # convert to this format
-        if isspmatrix_bsr(other):
+        if issparse(other) and other.format == "bsr":
path: scipy/sparse/_bsr.py
text: @user1:
Technically we don't need to call issparse here because we're in _mul_sparse_matrix

diff_hunk: @@ -91,7 +91,7 @@ def traceest(A, m3, seed=None):
 
 def _ident_like(A):
     # A compatibility function which should eventually disappear.
-    if scipy.sparse.isspmatrix(A):
+    if scipy.sparse.issparse(A):
         return scipy.sparse._construct.eye(A.shape[0], A.shape[1],
path: scipy/sparse/linalg/_expm_multiply.py
text: @user1:
If given a sparse array, this will return a spmatrix. We should wrap this such that it returns a sparse array when `A` is a sparse array.