Repository: scipy/scipy
Pull request: ENH: sparse: nD cleanup and docs
Fixes #22112 

This cleanup was inspired from the stubs work in #22112 which showed that:
- CSC format could be created as an nD array.
- `tensordot` had no doc_string

While investigating I found CSR could also almost create nD...and the resulting error was not clear.
Also, the new function `dot` did not have a doc_string.  In addition, there is quite q bit of spurious whitespace from the nD PRs.  I also found and fixed an idx_dtype error in the new `_broadcast_to` method. And finally that `test_base.py` have tests for sum and mean that were copied and never changed the var names from mean to sum. So, there are a number of small changes here.  Each has it's own commit.

-remove spurious whitespace
-align exception messages for 3D input across formats (adds test too)
- fix idx_dtype error in `_broadcast_to`
- change test_base.py varnames to match quantity computed
- add doc_strings to dot and tensordot with examples.

This should be backported to the maintenance branch. 

Problem statement:
BUG/DOC: sparse: ND COO unexpected behaviour 1.15.0rc1
### Describe your issue.

`csc_array` accepts dense array-likes since 1.15.0rc1, and accepts other `csc_array` instances with `ndim > 2`, but raises a `ValueError` if passed a `coo_array` with `ndim > 2`.

### Reproducing Code Example

```pycon
>>> import numpy as np
>>> from scipy.sparse import coo_array, csc_array
>>> csc_array(np.zeros((2, 3, 4)))
<Compressed Sparse Column sparse array of dtype 'float64'
	with 0 stored elements and shape (2, 3, 4)>
>>> coo_array(np.zeros((2, 3, 4)))
<COOrdinate sparse array of dtype 'float64'
	with 0 stored elements and shape (2, 3, 4)>
>>> csc_array(_)
Traceback (most recent call last):
  File "<python-input-10>", line 1, in <module>
    csc_array(_)
    ~~~~~~~~~^^^
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_compressed.py", line 36, in __init__
    arg1 = arg1.asformat(self.format)
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_base.py", line 456, in asformat
    return convert_method(copy=copy)
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_coo.py", line 339, in tocsc
    raise ValueError(f'Cannot convert. CSC format must be 2D. Got {self.ndim}D')
ValueError: Cannot convert. CSC format must be 2D. Got 3D
```


### Error message

```shell
Traceback (most recent call last):
  File "<python-input-10>", line 1, in <module>
    csc_array(_)
    ~~~~~~~~~^^^
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_compressed.py", line 36, in __init__
    arg1 = arg1.asformat(self.format)
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_base.py", line 456, in asformat
    return convert_method(copy=copy)
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_coo.py", line 339, in tocsc
    raise ValueError(f'Cannot convert. CSC format must be 2D. Got {self.ndim}D')
ValueError: Cannot convert. CSC format must be 2D. Got 3D
```


### SciPy/NumPy/Python version and system information

```shell
1.15.0rc1 2.2.0 sys.version_info(major=3, minor=13, micro=1, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /opt/_internal/cpython-3.13.0/lib/python3.13/site-packages/scipy_openblas32/include
    lib directory: /opt/_internal/cpython-3.13.0/lib/python3.13/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.28 DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=64
    pc file directory: /project
    version: 0.3.28
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /opt/_internal/cpython-3.13.0/lib/python3.13/site-packages/scipy_openblas32/include
    lib directory: /opt/_internal/cpython-3.13.0/lib/python3.13/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.28 DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=64
    pc file directory: /project
    version: 0.3.28
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.13.6
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.11
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-fa6gfmf0/overlay/lib/python3.13/site-packages/pythran
    version: 0.17.0
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: /opt/python/cp313-cp313/bin/python
  version: '3.13'
```

Hints from issue discussion:
Is this a break? I thought the maximum `ndim` for COO arrays was `2` before this release?
> Is this a break? I thought the maximum `ndim` for COO arrays was `2` before this release?

It's mentioned in the [releaase notes](https://github.com/scipy/scipy/releases/tag/v1.15.0rc1):

> n-D sparse arrays in COO format can now be constructed, reshaped and used for basic arithmetic.

... although it doesn't mention the >=2-d support for CSC  🤔 
Oh and I just noticed that indexing a `csc_array` with `ndim > 2` also doesn't seem to work:

```pycon
>>> csc_array(np.ones((2, 3, 4)))[0, 0, 0]
Traceback (most recent call last):
  File "<python-input-17>", line 1, in <module>
    csc_array(np.ones((2, 3, 4)))[0, 0, 0]
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_index.py", line 57, in __getitem__
    row, col = index
    ^^^^^^^^
ValueError: too many values to unpack (expected 2)
```

```pycon
>>> csc_array(np.ones((2, 3, 4)))[0, 0]
Traceback (most recent call last):
  File "<python-input-18>", line 1, in <module>
    csc_array(np.ones((2, 3, 4)))[0, 0]
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_index.py", line 57, in __getitem__
    row, col = index
    ^^^^^^^^
ValueError: too many values to unpack (expected 2)
```

```pycon
>>> csc_array(np.ones((2, 3, 4)))[0]
Traceback (most recent call last):
  File "<python-input-19>", line 1, in <module>
    csc_array(np.ones((2, 3, 4)))[0]
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^
  File "/home/joren/.pyenv/versions/3.13.1/lib/python3.13/site-packages/scipy/sparse/_index.py", line 57, in __getitem__
    row, col = index
    ^^^^^^^^
ValueError: too many values to unpack (expected 2)
```

So, should I add shape-typing support to `csc_coo` in https://github.com/jorenham/scipy-stubs/pull/349, or will it be reverted to only 2-d? 
(I just finished the >2-d shape-typing support for `coo_array` 😅)
yeah, I don't think this is documented to work anywhere - is there a specific part of the release notes which seems to be missing?
Ohh, are these examples of things which worked in 1.14.1?
I haven't seen anything myself either. But there a similar "surprise features" in `scipy.sparse`: `coo_array.tensordot`, which doesn't even have a docstring (I did annotate it though) 🤷🏻

> Ohh, are these examples of things which worked in 1.14.1?

No `csc` and `coo` both didn't support >2-d then. `csc` only supported 2-d, and `coo` either 1-d or 2-d. The 1-d + 2-d support of the `dok`  and `csr` arrays didn't change (but I haven't seen any documentation on the 1-d support of those two either)
Okay, so what specifically is this issue about? Do you think we need to make any changes for 1.15.0?

> But there a similar "surprise features" in scipy.sparse: coo_array.tensordot, which doesn't even have a docstring

That does seem like an oversight, cc @dschult @perimosocordiae 
> Okay, so what specifically is this issue about? Do you think we need to make any changes for 1.15.0?


Either remove it before anyone notices (while nervously whistling, of course), or fix and document it before `1.15.0` final. 

I only quickly tried out the >2d `csc_matrix` and `coo_matrix.tensordot`. And of the 5 or so things I tried with `csc_matrix`, 2 turned out to be bugged. The `coo_matrix.tensordot` is missing a docstring, and raises a `TypeError` once you pass a 2-tuple or list to its `axes` parameter. The `axes=2` default also causes errors when "tensordotting" two 1-d arrays (sparse or dense doesn't matter), as that requires setting `axes=1`. 

Given that I found these bugs this quickly, then it's likely that more will pop up if looked at more thoroughly.

So maybe it's a bit too ambitious if both will need to be to documented, better tested, and bug-fixed, within the two weeks until the 1.15.0 final releaase? 
At least at this point no-one will notice that they won't be included in the final release, given the missing release notes 🤷🏻 
Are these bugs? Or just undocumented intended behaviour?

Commit message: add docstrings to dot and tensordot

Patch to review:
diff --git a/scipy/sparse/_compressed.py b/scipy/sparse/_compressed.py
index 1be38322e3ff..94825e674768 100644
--- a/scipy/sparse/_compressed.py
+++ b/scipy/sparse/_compressed.py
@@ -86,8 +86,11 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False, *, maxprint=None):
             except Exception as e:
                 raise ValueError(f"unrecognized {self.__class__.__name__} "
                                  f"constructor input: {arg1}") from e
-            if isinstance(self, sparray) and arg1.ndim < 2 and self.format == "csc":
+            if isinstance(self, sparray) and arg1.ndim != 2 and self.format == "csc":
                 raise ValueError(f"CSC arrays don't support {arg1.ndim}D input. Use 2D")
+            if arg1.ndim > 2:
+                raise ValueError(f"CSR arrays don't yet support {arg1.ndim}D.")
+
             coo = self._coo_container(arg1, dtype=dtype)
             arrays = coo._coo_to_compressed(self._swap)
             self.indptr, self.indices, self.data, self._shape = arrays
@@ -1415,17 +1418,17 @@ def _divide_sparse(self, other):
     def _broadcast_to(self, shape, copy=False):
         if self.shape == shape:
             return self.copy() if copy else self
-        
+
         shape = check_shape(shape, allow_nd=(self._allow_nd))
 
         if broadcast_shapes(self.shape, shape) != shape:
             raise ValueError("cannot be broadcast")
-        
+
         if len(self.shape) == 1 and len(shape) == 1:
             self.sum_duplicates()
             if self.nnz == 0: # array has no non zero elements
                 return self.__class__(shape, dtype=self.dtype, copy=False)
-            
+
             N = shape[0]
             data = np.full(N, self.data[0])
             indices = np.arange(0,N)
@@ -1434,14 +1437,14 @@ def _broadcast_to(self, shape, copy=False):
 
         # treat 1D as a 2D row
         old_shape = self._shape_as_2d
-            
+
         if len(shape) != 2:
             ndim = len(shape)
             raise ValueError(f'CSR/CSC broadcast_to cannot have shape >2D. Got {ndim}D')
-        
+
         if self.nnz == 0: # array has no non zero elements
             return self.__class__(shape, dtype=self.dtype, copy=False)
-        
+
         self.sum_duplicates()
         M, N = self._swap(shape)
         oM, oN = self._swap(old_shape)
diff --git a/scipy/sparse/_coo.py b/scipy/sparse/_coo.py
index 08d737428b30..efaba6cb9c68 100644
--- a/scipy/sparse/_coo.py
+++ b/scipy/sparse/_coo.py
@@ -667,24 +667,24 @@ def _rmatmul_dispatch(self, other):
                 o_ndim = other.ndim
             perm = tuple(range(o_ndim)[:-2]) + tuple(range(o_ndim)[-2:][::-1])
             tr = other.transpose(perm)
-            
+
             s_ndim = self.ndim
             perm = tuple(range(s_ndim)[:-2]) + tuple(range(s_ndim)[-2:][::-1])
             ret = self.transpose(perm)._matmul_dispatch(tr)
             if ret is NotImplemented:
                 return NotImplemented
-            
+
             if s_ndim == 1 or o_ndim == 1:
                 perm = range(ret.ndim)
             else:
                 perm = tuple(range(ret.ndim)[:-2]) + tuple(range(ret.ndim)[-2:][::-1])
             return ret.transpose(perm)
-        
+
 
     def _matmul_dispatch(self, other):
         if isscalarlike(other):
             return self.multiply(other)
-        
+
         if not (issparse(other) or isdense(other)):
             # If it's a list or whatever, treat it like an array
             other_a = np.asanyarray(other)
@@ -723,14 +723,14 @@ def _matmul_dispatch(self, other):
                         np.broadcast_shapes(batch_shape_A, batch_shape_B)
                     except ValueError:
                         raise ValueError("Batch dimensions are not broadcastable")
-            
+
                 return self._matmul_multivector(other)
             else:
                 raise ValueError(
                     f"{err_prefix} (n,..,k={N}),(k={other.shape[-2]},..,m)->(n,..,m)"
                 )
-        
-            
+
+
         if isscalarlike(other):
             # scalar value
             return self._mul_scalar(other)
@@ -751,7 +751,7 @@ def _matmul_dispatch(self, other):
                 raise ValueError(
                     f"{err_prefix} (n,..,k={N}),(k={other.shape[-2]},..,m)->(n,..,m)"
                 )
-            
+
             # If A or B has more than 2 dimensions, check for
             # batch dimensions compatibility
             if self.ndim > 2 or other.ndim > 2:
@@ -763,7 +763,7 @@ def _matmul_dispatch(self, other):
                         np.broadcast_shapes(batch_shape_A, batch_shape_B)
                     except ValueError:
                         raise ValueError("Batch dimensions are not broadcastable")
-            
+
             result = self._matmul_sparse(other)
 
             # reshape back if a or b were originally 1-D
@@ -783,11 +783,11 @@ def _matmul_multivector(self, other):
             if self.ndim == 1:
                 result = self.reshape(1, self.shape[0])._matmul_multivector(other)
                 return result.reshape(tuple(other.shape[:-2]) + tuple(other.shape[-1:]))
-            
+
             broadcast_shape = np.broadcast_shapes(self.shape[:-2], other.shape[:-2])
             self_shape = broadcast_shape + self.shape[-2:]
             other_shape = broadcast_shape + other.shape[-2:]
-            
+
             self = self._broadcast_to(self_shape)
             other = np.broadcast_to(other, other_shape)
             result_shape = broadcast_shape + self.shape[-2:-1] + other.shape[-1:]
@@ -797,7 +797,7 @@ def _matmul_multivector(self, other):
                                 np.concatenate(self.coords),
                                 self.data, other.ravel('C'), result)
             return result
-        
+
         if self.ndim == 2:
             result_shape = (self.shape[0], other.shape[1])
             col = self.col
@@ -813,6 +813,58 @@ def _matmul_multivector(self, other):
 
 
     def dot(self, other):
+        """Return the dot product of two arrays.
+
+        Strictly speaking a dot product involves two vectors.
+        But in the sense that an array with ndim >= 1 is a collection
+        of vectors, the function computes the collection of dot products
+        between each vector in the first array with each vector in the
+        second array. The axis upon which the sum of products is performed
+        is the last axis of the first array and the second to last axis of
+        the second array. If the second array is 1-D, the last axis is used.
+
+        Thus, if both arrays are 1-D, the inner product is returned.
+        If both are 2-D, we have matrix multiplication. If `other` is 1-D,
+        the sum product is taken along the last axis of each array. If
+        `other` is N-D for N>=2, the sum product is over the last axis of
+        the first array and the second-to-last axis of the second array.
+
+        Parameters
+        ----------
+        other : array_like (dense of sparse)
+            Second array
+
+        Returns
+        -------
+        output : array (sparse or dense)
+            The dot product of this array with `other`.
+            It will be dense/sparse if `other` is dense/sparse.
+
+        Examples
+        --------
+
+        >>> import numpy as np
+        >>> from scipy.sparse import coo_array
+        >>> A = coo_array([[1, 2, 0], [0, 0, 3], [4, 0, 5]])
+        >>> v = np.array([1, 0, -1])
+        >>> A.dot(v)
+        array([ 1, -3, -1], dtype=int64)
+
+        For 2-D arrays it is the matrix product:
+
+        >>> A = coo_array([[1, 0], [0, 1]])
+        >>> B = coo_array([[4, 1], [2, 2]])
+        >>> A.dot(B).toarray()
+        array([[4, 1],
+               [2, 2]])
+
+        For 3-D arrays the shape extends unused axes by other unused axes.
+
+        >>> A = coo_array(np.arange(3*4*5*6)).reshape((3,4,5,6))
+        >>> B = coo_array(np.arange(3*4*5*6)).reshape((5,4,6,3))
+        >>> A.dot(B).shape
+        (3, 4, 5, 5, 4, 3)
+        """
         if not (issparse(other) or isdense(other) or isscalarlike(other)):
             # If it's a list or whatever, treat it like an array
             o_array = np.asanyarray(other)
@@ -851,7 +903,7 @@ def dot(self, other):
         else:
             return self._sparse_dot(other)
 
-    
+
     def _sparse_dot(self, other):
         self_is_1d = self.ndim == 1
         other_is_1d = other.ndim == 1
@@ -865,7 +917,7 @@ def _sparse_dot(self, other):
         if self.shape[-1] != other.shape[-2]:
                 raise ValueError(f"shapes {self.shape} and {other.shape}"
                                  " are not aligned for n-D dot")
-        
+
         # Prepare the tensors for dot operation
         # Ravel non-reduced axes coordinates
         self_raveled_coords = _ravel_non_reduced_axes(self.coords,
@@ -876,17 +928,17 @@ def _sparse_dot(self, other):
         # Get the shape of the non-reduced axes
         self_nonreduced_shape = self.shape[:-1]
         other_nonreduced_shape = other.shape[:-2] + other.shape[-1:]
-        
+
         # Create 2D coords arrays
         ravel_coords_shape_self = (math.prod(self_nonreduced_shape), self.shape[-1])
         ravel_coords_shape_other = (other.shape[-2], math.prod(other_nonreduced_shape))
-        
+
         self_2d_coords = (self_raveled_coords, self.coords[-1])
         other_2d_coords = (other.coords[-2], other_raveled_coords)
 
         self_2d = coo_array((self.data, self_2d_coords), ravel_coords_shape_self)
         other_2d = coo_array((other.data, other_2d_coords), ravel_coords_shape_other)
-        
+
         prod = (self_2d @ other_2d).tocoo() # routes via 2-D CSR
 
         # Combine the shapes of the non-reduced axes
@@ -899,7 +951,7 @@ def _sparse_dot(self, other):
             prod_coords.extend(np.unravel_index(c, s))
 
         prod_arr = coo_array((prod.data, prod_coords), combined_shape)
-        
+
         # reshape back if a or b were originally 1-D
         # TODO: Move this logic before computation of prod_coords for efficiency
         if self_is_1d:
@@ -908,7 +960,7 @@ def _sparse_dot(self, other):
             prod_arr = prod_arr.reshape(combined_shape[:-1])
 
         return prod_arr
-    
+
     def _dense_dot(self, other):
         self_is_1d = self.ndim == 1
         other_is_1d = other.ndim == 1
@@ -941,13 +993,100 @@ def _dense_dot(self, other):
         return prod_arr
 
     def tensordot(self, other, axes=2):
+        """Return the tensordot product with another array along the given axes.
+
+        The tensordot differs from dot and matmul in that any axis can be
+        chosen for each of the first and second array and the sum of the
+        producs is computed just like for matrix multiplication, only not
+        just for the rows fo the first times the columns of the second. It
+        takes the dot product of the collection of vectors along the specified
+        axes.  Here we can even take the sum of the products along two or even
+        more axes if desired. So, tensordot is a dot product computation
+        applied to arrays of any dimension >= 1. It is like matmul but over
+        arbitrary axes for each matrix.
+
+        Given two tensors, `a` and `b`, and the desired axes specified as a
+        2-tuple/list/array containing two sequences of axis numbers,
+        ``(a_axes, b_axes)``, sum the products of `a`'s and `b`'s elements
+        (components) over the axes specified by``a_axes`` and ``b_axes``.
+        The `axes` input can be a single non-negative integer, ``N``;
+        if it is, then the last ``N`` dimensions of `a` and the first
+        ``N`` dimensions of `b` are summed over.
+
+        Parameters
+        ----------
+        a, b : array_like
+            Tensors to "dot".
+
+        axes : int or (2,) array_like
+            * integer_like
+              If an int N, sum over the last N axes of `a` and the first N axes
+              of `b` in order. The sizes of the corresponding axes must match.
+            * (2,) array_like
+              A 2-tuple of sequences of axes to be summed over, the first applying
+              to `a`, the second to `b`. The sequences must be the same length.
+              The shape of the corresponding axes must match between `a` and `b`.
+
+        Returns
+        -------
+        output : coo_array
+            The tensor dot product of this array with `other`.
+            It will be dense/sparse if `other` is dense/sparse.
+
+        See Also
+        --------
+        dot
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> import scipy.sparse
+        >>> A = scipy.sparse.coo_array([[[2, 3], [0, 0]], [[0, 1], [0, 5]]])
+        >>> A.shape
+        (2, 2, 2)
+
+        Integer axes N are shorthand for (range(-N, 0), range(0, N)):
+
+        >>> A.tensordot(A, axes=1).toarray()
+        array([[[[ 4,  9],
+                 [ 0, 15]],
+        <BLANKLINE>
+                [[ 0,  0],
+                 [ 0,  0]]],
+        <BLANKLINE>
+        <BLANKLINE>
+               [[[ 0,  1],
+                 [ 0,  5]],
+        <BLANKLINE>
+                [[ 0,  5],
+                 [ 0, 25]]]])
+        >>> A.tensordot(A, axes=2).toarray()
+        array([[ 4,  6],
+               [ 0, 25]])
+        >>> A.tensordot(A, axes=3)
+        array(39)
+
+        Using tuple for axes:
+
+        >>> a = scipy.sparse.coo_array(np.arange(60).reshape(3,4,5))
+        >>> b = np.arange(24).reshape(4,3,2)
+        >>> c = a.tensordot(b, axes=([1,0],[0,1]))
+        >>> c.shape
+        (5, 2)
+        >>> c
+        array([[4400, 4730],
+               [4532, 4874],
+               [4664, 5018],
+               [4796, 5162],
+               [4928, 5306]])
+
+        """
         if not isdense(other) and not issparse(other):
             # If it's a list or whatever, treat it like an array
             other_array = np.asanyarray(other)
 
             if other_array.ndim == 0 and other_array.dtype == np.object_:
-                # Not interpretable as an array; return NotImplemented so that
-                # other's __rmatmul__ can kick in if that's implemented.
+                # Not interpretable as an array; return NotImplemented
                 return NotImplemented
 
             try:
@@ -961,7 +1100,7 @@ def tensordot(self, other, axes=2):
         if any(self.shape[ax] != other.shape[bx]
                for ax, bx in zip(axes_self, axes_other)):
             raise ValueError("sizes of the corresponding axes must match")
-        
+
         if isdense(other):
             return self._dense_tensordot(other, axes_self, axes_other)
         else:
@@ -972,7 +1111,7 @@ def _sparse_tensordot(self, other, axes_self, axes_other):
         ndim_self = len(self.shape)
         ndim_other = len(other.shape)
 
-        # Prepare the tensors for tensordot operation       
+        # Prepare the tensors for tensordot operation
         # Ravel non-reduced axes coordinates
         self_non_red_coords = _ravel_non_reduced_axes(self.coords, self.shape,
                                                       axes_self)
@@ -988,7 +1127,7 @@ def _sparse_tensordot(self, other, axes_self, axes_other):
                               if ax not in axes_self)
         other_nonreduced_shape = tuple(other.shape[ax] for ax in range(ndim_other)
                                if ax not in axes_other)
-        
+
         # Create 2D coords arrays
         ravel_coords_shape_self = (math.prod(self_nonreduced_shape),
                                 math.prod([self.shape[ax] for ax in axes_self]))
@@ -1015,7 +1154,7 @@ def _sparse_tensordot(self, other, axes_self, axes_other):
 
         if coords == []:  # if result is scalar
             return sum(prod.data)
-            
+
         # Construct the resulting COO array with combined coordinates and shape
         return coo_array((prod.data, coords), shape=combined_shape)
 
@@ -1053,11 +1192,11 @@ def _matmul_sparse(self, other):
         The method converts input n-D arrays to 2-D block array format,
         uses csr_matmat to multiply them, and then converts the
         result back to n-D COO array.
-        
+
         Parameters:
         self (COO): The first n-D sparse array in COO format.
         other (COO): The second n-D sparse array in COO format.
-        
+
         Returns:
         prod (COO): The resulting n-D sparse array after multiplication.
         """
@@ -1067,7 +1206,7 @@ def _matmul_sparse(self, other):
         # Get the shapes of self and other
         self_shape = self.shape
         other_shape = other.shape
-        
+
         # Determine the new shape to broadcast self and other
         broadcast_shape = np.broadcast_shapes(self_shape[:-2], other_shape[:-2])
         self_new_shape = tuple(broadcast_shape) + self_shape[-2:]
@@ -1075,14 +1214,14 @@ def _matmul_sparse(self, other):
 
         self_broadcasted = self._broadcast_to(self_new_shape)
         other_broadcasted = other._broadcast_to(other_new_shape)
-        
+
         # Convert n-D COO arrays to 2-D block diagonal arrays
         self_block_diag = _block_diag(self_broadcasted)
         other_block_diag = _block_diag(other_broadcasted)
-        
+
         # Use csr_matmat to perform sparse matrix multiplication
         prod_block_diag = (self_block_diag @ other_block_diag).tocoo()
-        
+
         # Convert the 2-D block diagonal array back to n-D
         return _extract_block_diag(
             prod_block_diag,
@@ -1093,17 +1232,17 @@ def _matmul_sparse(self, other):
     def _broadcast_to(self, new_shape, copy=False):
         if self.shape == new_shape:
             return self.copy() if copy else self
-        
+
         old_shape = self.shape
 
         # Check if the new shape is compatible for broadcasting
         if len(new_shape) < len(old_shape):
             raise ValueError("New shape must have at least as many dimensions"
                              " as the current shape")
-        
-        # Add leading ones to shape to ensure same length as `new_shape` 
+
+        # Add leading ones to shape to ensure same length as `new_shape`
         shape = (1,) * (len(new_shape) - len(old_shape)) + tuple(old_shape)
-        
+
         # Ensure the old shape can be broadcast to the new shape
         if any((o != 1 and o != n) for o, n in zip(shape, new_shape)):
             raise ValueError(f"current shape {old_shape} cannot be "
@@ -1112,18 +1251,19 @@ def _broadcast_to(self, new_shape, copy=False):
         # Reshape the COO array to match the new dimensions
         self = self.reshape(shape)
 
+        idx_dtype = get_index_dtype(self.coords, maxval=max(new_shape))
         coords = self.coords
         new_data = self.data
         new_coords = coords[-1:]  # Copy last coordinate to start
         cum_repeat = 1 # Cumulative repeat factor for broadcasting
-        
+
         if shape[-1] != new_shape[-1]: # broadcasting the n-th (col) dimension
             repeat_count = new_shape[-1]
             cum_repeat *= repeat_count
             new_data = np.tile(new_data, repeat_count)
-            new_dim = np.repeat(np.arange(0, repeat_count), self.nnz)
+            new_dim = np.repeat(np.arange(0, repeat_count, dtype=idx_dtype), self.nnz)
             new_coords = (new_dim,)
-        
+
         for i in range(-2, -(len(shape)+1), -1):
             if shape[i] != new_shape[i]:
                 repeat_count = new_shape[i] # number of times to repeat data, coords
@@ -1135,15 +1275,15 @@ def _broadcast_to(self, new_shape, copy=False):
                 new_coords = tuple(np.tile(new_coords[i+1:], repeat_count))
 
                 # Create new dimensions and stack them
-                new_dim = np.repeat(np.arange(0, repeat_count), nnz)
+                new_dim = np.repeat(np.arange(0, repeat_count, dtype=idx_dtype), nnz)
                 new_coords = (new_dim,) + new_coords
             else:
                 # If no broadcasting needed, tile the coordinates
                 new_dim = np.tile(coords[i], cum_repeat)
                 new_coords = (new_dim,) + new_coords
-                
+
         return coo_array((new_data, new_coords), new_shape)
-        
+
 
 def _block_diag(self):
     """
@@ -1179,7 +1319,7 @@ def _extract_block_diag(self, shape):
 
     # Initialize new coordinates array
     new_coords = np.empty((len(shape), self.nnz), dtype=int)
-    
+
     # Calculate within-block indices
     new_coords[-2] = row % n_row
     new_coords[-1] = col % n_col
@@ -1212,7 +1352,7 @@ def _process_axes(ndim_a, ndim_b, axes):
             raise ValueError("axes indices are out of bounds for input arrays")
     else:
         raise TypeError("axes must be an integer or a tuple/list of integers")
-    
+
     axes_a = [axis + ndim_a if axis < 0 else axis for axis in axes_a]
     axes_b = [axis + ndim_b if axis < 0 else axis for axis in axes_b]
     return axes_a, axes_b
@@ -1225,13 +1365,13 @@ def _ravel_non_reduced_axes(coords, shape, axes):
     if not non_reduced_axes:
         # Return an array with one row
         return np.zeros_like(coords[0])
-    
+
     # Extract the shape of the non-reduced axes
     non_reduced_shape = [shape[ax] for ax in non_reduced_axes]
-    
+
     # Extract the coordinates of the non-reduced axes
     non_reduced_coords = tuple(coords[idx] for idx in non_reduced_axes)
-    
+
     # Ravel the coordinates into 1D
     return np.ravel_multi_index(non_reduced_coords, non_reduced_shape)
 
diff --git a/scipy/sparse/_dok.py b/scipy/sparse/_dok.py
index c2fc2f3ccc5a..f539589973c5 100644
--- a/scipy/sparse/_dok.py
+++ b/scipy/sparse/_dok.py
@@ -44,7 +44,7 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False, *, maxprint=None):
                 raise TypeError('Invalid input format.') from e
 
             if arg1.ndim > 2:
-                raise TypeError('Expected rank <=2 dense array or matrix.')
+                raise ValueError(f"DOK arrays don't yet support {arg1.ndim}D input.")
 
             if arg1.ndim == 1:
                 if dtype is not None:
diff --git a/scipy/sparse/tests/test_base.py b/scipy/sparse/tests/test_base.py
index 0f8fb0cd347a..3fa52e0d17f2 100644
--- a/scipy/sparse/tests/test_base.py
+++ b/scipy/sparse/tests/test_base.py
@@ -1057,11 +1057,11 @@ def test_sum_dtype(self):
         datsp = self.spcreator(dat)
 
         def check(dtype):
-            dat_mean = dat.mean(dtype=dtype)
-            datsp_mean = datsp.mean(dtype=dtype)
+            dat_sum = dat.sum(dtype=dtype)
+            datsp_sum = datsp.sum(dtype=dtype)
 
-            assert_array_almost_equal(dat_mean, datsp_mean)
-            assert_equal(dat_mean.dtype, datsp_mean.dtype)
+            assert_array_almost_equal(dat_sum, datsp_sum)
+            assert_equal(dat_sum.dtype, datsp_sum.dtype)
 
         for dtype in self.checked_dtypes:
             check(dtype)
@@ -1094,11 +1094,11 @@ def test_numpy_sum(self):
                      [-6, 7, 9]])
         datsp = self.spcreator(dat)
 
-        dat_mean = np.sum(dat)
-        datsp_mean = np.sum(datsp)
+        dat_sum = np.sum(dat)
+        datsp_sum = np.sum(datsp)
 
-        assert_array_almost_equal(dat_mean, datsp_mean)
-        assert_equal(dat_mean.dtype, datsp_mean.dtype)
+        assert_array_almost_equal(dat_sum, datsp_sum)
+        assert_equal(dat_sum.dtype, datsp_sum.dtype)
 
     def test_mean(self):
         keep = not self.is_array_test
diff --git a/scipy/sparse/tests/test_common1d.py b/scipy/sparse/tests/test_common1d.py
index 21f25b79d4cd..e1208d67cef6 100644
--- a/scipy/sparse/tests/test_common1d.py
+++ b/scipy/sparse/tests/test_common1d.py
@@ -43,6 +43,16 @@ def test_no_1d_support_in_init(spcreator):
         spcreator([0, 1, 2, 3])
 
 
+# Test init with nD dense input
+# sparrays which do not yet support nD
+@pytest.mark.parametrize(
+    "spcreator", [csr_array, dok_array, bsr_array, csc_array, dia_array, lil_array]
+)
+def test_no_nd_support_in_init(spcreator):
+    with pytest.raises(ValueError, match="arrays don't.*support 3D"):
+        spcreator(np.ones((3, 2, 4)))
+
+
 # Main tests class
 @pytest.mark.parametrize("spcreator", spcreators)
 class TestCommon1D:
diff --git a/scipy/sparse/tests/test_coo.py b/scipy/sparse/tests/test_coo.py
index 6e7f97bcd646..a163ce3d1460 100644
--- a/scipy/sparse/tests/test_coo.py
+++ b/scipy/sparse/tests/test_coo.py
@@ -418,7 +418,7 @@ def test_tuple_constructor_for_dim_size_zero():
 
     empty_arr = coo_array(([], ([], [], [], [])), shape=(4,0,2,3))
     assert_equal(empty_arr.toarray(), np.empty((4,0,2,3)))
-    
+
 
 @pytest.mark.parametrize(('shape', 'new_shape'), [((4,9,6,5), (3,6,15,4)),
                                                   ((4,9,6,5), (36,30)),
@@ -626,8 +626,8 @@ def test_nd_matmul_vector(mat_shape, vec_shape):
 
 
 mat_mat_shapes = [
-    ((2, 3, 4, 5), (2, 3, 5, 7)), 
-    ((0, 0), (0,)), 
+    ((2, 3, 4, 5), (2, 3, 5, 7)),
+    ((0, 0), (0,)),
     ((4, 4, 2, 0), (0,)),
     ((7, 8, 3), (3,)),
     ((7, 8, 3), (3, 1)),
@@ -670,7 +670,7 @@ def test_nd_matmul_sparse_with_inconsistent_arrays():
         sp_x @ sp_y
     with pytest.raises(ValueError, match="matmul: dimension mismatch with signature"):
         sp_x @ (sp_y.toarray())
-    
+
     sp_z = random_array((1,5,3,2), density=0.6, random_state=rng, dtype=int)
     with pytest.raises(ValueError, match="Batch dimensions are not broadcastable"):
         sp_x @ sp_z
@@ -688,7 +688,7 @@ def test_dot_1d_1d(): # 1-D inner product
     assert_equal(res, exp)
 
 
-def test_dot_sparse_scalar():    
+def test_dot_sparse_scalar():
     a = coo_array([[1, 2], [3, 4], [5, 6]])
     b = 3
     res = a.dot(b)
@@ -696,7 +696,7 @@ def test_dot_sparse_scalar():
     assert_equal(res.toarray(), exp)
 
 
-def test_dot_with_inconsistent_shapes(): 
+def test_dot_with_inconsistent_shapes():
     arr_a = coo_array([[[1, 2]], [[3, 4]]])
     arr_b = coo_array([4, 5, 6])
     with pytest.raises(ValueError, match="not aligned for n-D dot"):
@@ -708,12 +708,12 @@ def test_dot_with_inconsistent_shapes():
     ((3,2,4,7), (7,)), ((5,), (6,3,5,2)), # dot of n-D and 1-D arrays
     ((3,2,4,7), (7,1)), ((1,5,), (6,3,5,2)),
     ((4,6), (3,2,6,4)), ((2,8,7), (4,5,7,7,2)), # dot of n-D and m-D arrays
-    ((4,5,7,6), (3,2,6,4)), 
+    ((4,5,7,6), (3,2,6,4)),
 ]
 @pytest.mark.parametrize(('a_shape', 'b_shape'), dot_shapes)
-def test_dot_nd(a_shape, b_shape): 
+def test_dot_nd(a_shape, b_shape):
     rng = np.random.default_rng(23409823)
-    
+
     arr_a = random_array(a_shape, density=0.6, random_state=rng, dtype=int)
     arr_b = random_array(b_shape, density=0.6, random_state=rng, dtype=int)
 
@@ -739,9 +739,9 @@ def test_dot_nd(a_shape, b_shape):
     ((2,3,4), (2,3,4), ([0, 1, 2], [0, 1, 2])),
 ]
 @pytest.mark.parametrize(('a_shape', 'b_shape', 'axes'), tensordot_shapes_and_axes)
-def test_tensordot(a_shape, b_shape, axes): 
+def test_tensordot(a_shape, b_shape, axes):
     rng = np.random.default_rng(23409823)
-    
+
     arr_a = random_array(a_shape, density=0.6, random_state=rng, dtype=int)
     arr_b = random_array(b_shape, density=0.6, random_state=rng, dtype=int)
 
@@ -759,11 +759,11 @@ def test_tensordot(a_shape, b_shape, axes):
         assert_equal(res.toarray(), exp)
     else:
         assert_equal(res, exp)
-        
 
-def test_tensordot_with_invalid_args(): 
+
+def test_tensordot_with_invalid_args():
     rng = np.random.default_rng(23409823)
-    
+
     arr_a = random_array((3,4,5), density=0.6, random_state=rng, dtype=int)
     arr_b = random_array((3,4,6), density=0.6, random_state=rng, dtype=int)
 
@@ -787,7 +787,7 @@ def test_tensordot_with_invalid_args():
                           ((1,1,1), (4,5,6)), ((1,3,1,5,4), (8,2,3,9,5,4)),])
 def test_broadcast_to(actual_shape, broadcast_shape):
     rng = np.random.default_rng(23409823)
-    
+
     arr = random_array(actual_shape, density=0.6, random_state=rng, dtype=int)
     res = arr._broadcast_to(broadcast_shape)
     exp = np.broadcast_to(arr.toarray(), broadcast_shape)
diff --git a/scipy/sparse/tests/test_csr.py b/scipy/sparse/tests/test_csr.py
index a7ede4ab3c1b..6b011ad4fdce 100644
--- a/scipy/sparse/tests/test_csr.py
+++ b/scipy/sparse/tests/test_csr.py
@@ -203,7 +203,7 @@ def test_broadcast_to():
     assert_array_equal(res_c.toarray(), np.broadcast_to(b, (2,4)))
     assert_array_equal(res_d.toarray(), np.broadcast_to(b, (1,)))
     assert_array_equal(res_e.toarray(), np.broadcast_to(e, (4,0)))
-    
+
     with pytest.raises(ValueError, match="cannot be broadcast"):
         csr_matrix([[1, 2, 0], [3, 0, 1]])._broadcast_to(shape=(2, 1))
 


Code review comments:
diff_hunk: @@ -941,13 +993,100 @@ def _dense_dot(self, other):
         return prod_arr
 
     def tensordot(self, other, axes=2):
+        """Return the tensordot product with another array along the given axes.
+
+        The tensordot differs from dot and matmul in that any axis can be
+        chosen for each of the first and second array and the sum of the
+        producs is computed just like for matrix multiplication, only not
+        just for the rows fo the first times the columns of the second. It
path: scipy/sparse/_coo.py
text: @user1:
```suggestion
        just for the rows of the first times the columns of the second. It
```

diff_hunk: @@ -813,6 +813,58 @@ def _matmul_multivector(self, other):
 
 
     def dot(self, other):
+        """Return the dot product of two arrays.
+
+        Strictly speaking a dot product involves two vectors.
+        But in the sense that an array with ndim >= 1 is a collection
+        of vectors, the function computes the collection of dot products
+        between each vector in the first array with each vector in the
+        second array. The axis upon which the sum of products is performed
+        is the last axis of the first array and the second to last axis of
+        the second array. If the second array is 1-D, the last axis is used.
+
+        Thus, if both arrays are 1-D, the inner product is returned.
+        If both are 2-D, we have matrix multiplication. If `other` is 1-D,
+        the sum product is taken along the last axis of each array. If
+        `other` is N-D for N>=2, the sum product is over the last axis of
+        the first array and the second-to-last axis of the second array.
+
+        Parameters
+        ----------
+        other : array_like (dense of sparse)
path: scipy/sparse/_coo.py
text: @user1:
```suggestion
        other : array_like (dense or sparse)
```

diff_hunk: @@ -941,13 +993,100 @@ def _dense_dot(self, other):
         return prod_arr
 
     def tensordot(self, other, axes=2):
+        """Return the tensordot product with another array along the given axes.
+
+        The tensordot differs from dot and matmul in that any axis can be
+        chosen for each of the first and second array and the sum of the
+        producs is computed just like for matrix multiplication, only not
path: scipy/sparse/_coo.py
text: @user1:
```suggestion
        products is computed just like for matrix multiplication, only not
```