RFC Supporting `scipy.sparse.sparray`

### Context

SciPy is now favoring sparse arrays (i.e. `scipy.sparse.sparray` and its subclasses) over sparse matrices (i.e. `scipy.sparse.spmatrix` and its subclasses) to enlarge the scope of matrices to $n$-dimensional data-structures since SciPy 1.8 (see https://github.com/scipy/scipy/pull/14822).

Sparse matrices now subclass their associated sparse arrays (see https://github.com/scipy/scipy/pull/18440).

scikit-learn has been supporting sparse matrices but now also needs to support SciPy sparse arrays.

### Proposed solutions

Ordered by preference:
 - Use `scipy.sparse.issparse` and the `format` attribute everywhere and not use `isinstance` at all
 - Use `isinstance` on private compatibility class defined to be `scipy.sparse.spmatrix` or `scipy.sparse.sparray` conditionally on SciPy's version (i.e. use `scipy.sparse.sparray` if available)
 - Rely on duck-typing or the class name to check for sparse arrays

cc @ivirshup

COMMENT:
Questions about how to handle migration from spmatrix to sparray when returning sparse:

- Many "feature" modules have classes that can hold either sparse matrix or sparse array. Often in `transform` or `fit` a sparse is created to be held or returned. If there is a sparse input for `X` and maybe even `y` we can match that sparse type. But what should we do when a dense is converted to sparse as part of `fit()` or `transform()`? I can get examples, or you can look in `text.py` `_hash.py`, `image.py`, etc. in the feature_extraction subpackage.
- About 12 functions/classes return sparse but have kwarg values like `sparse`, `sparse_input` `sparse_output`. Can we use those to determine what form of output? They are now bool. How to indicate spmatrix vs sparray?
- `linear_model.coordinate_descent` (and others) uses `sparse_coef` property or attribute. How to decide spmatrix/sparray? Make two properties?

Questions about specific functions:
- `datasets.load_svmlight_files` -> output is sparse. Does it matter which type?
- `dataset.fetch_rcv1` converts to csr_array for lexicographic order. Would a user need to be spmatrix?
- `barycenter_kneighbors_graph` -> returns newly created sparse object. how to handle?
- `classification._check_targets` -> converts ndarray to sparse. Existing comment lines question why we convert an ndarray to sparse just before returning. I agree. It is already using the memory. Should we switch to `ndarray`? How to handle spmatrix/sparray in the meantime?

I've looked through the other functions and think we can make the others work either with existing keyword args, or they just don't return sparse. So there might be ~30 functions affected between the above comments. 

I propose that:
- for `fit()` we should allow ourselves to keep whichever sparse format we like. User shouldn't care.
- for `transform` we look at type of `X`. If `X` is dense and converted to sparse, can we just shift to return sparray? Or is a kwarg needed?
- Any functions with existing bool inputs regarding "sparse" should use that kwarg to indicate sparray or spmatrix too. Maybe bool or string? Ideas?
- `sparse_coef` property can either always now be sparray. Or we make two properties like `sparse_array_coef`?

Then I'll need feedback on the specific functions with Questions above.
If you'd rather see this in PR form before committing one way or the other let me know and I'll put that together.

Thanks!

Labels: New Feature, API, Needs Decision, RFC