MAINT: sparse: Adjust expm and expm_multiply to pydata/sparse input

[mtsokol]

Hi!

According to the"array type out == array type in" design rule (#18286) an array that is converted to SciPy sparse matrix, needs to be converted to the original type before returning.

As a result, all functions in sparse.csgraph (#19796) and sparse.linalg first convert input objects to SciPy sparse objects before running desired function.

This should also apply to additional arrays that are created along the way. For example, in sparse.linalg.inv, I that is created should be a SciPy object.

Right now for pydata/sparse input, in inv(A), A gets converted to SciPy object, but I is created as a pydata/sparse object and then converted to SciPy object. The first step for I should already create SciPy object.

P.S. I think that 3rd party array -> SciPy dispatching should also happen in expm_multiply? That's the only place where _ident_like is used.

[perimosocordiae]

Seems reasonable to me.

[lucascolley]

CI is not happy:

=================================== FAILURES ===================================
________________________________ test_inv[csc] _________________________________
scipy/sparse/tests/test_array_api.py:204: in test_inv
    C = spla.inv(B)
        B          = <2x2 sparse array of type '<class 'numpy.int64'>'
	with 2 stored elements in Compressed Sparse Column format>
scipy/sparse/linalg/_matfuncs.py:76: in inv
    Ainv = spsolve(A, I)
        A          = <2x2 sparse array of type '<class 'numpy.int64'>'
	with 2 stored elements in Compressed Sparse Column format>
        I          = <2x2 sparse array of type '<class 'numpy.int64'>'
	with 2 stored elements (1 diagonals) in DIAgonal format>
scipy/sparse/linalg/_dsolve/linsolve.py:305: in spsolve
    warn('spsolve is more efficient when sparse b '
E   scipy.sparse._base.SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
        A          = <2x2 sparse array of type '<class 'numpy.float64'>'
	with 2 stored elements in Compressed Sparse Column format>
        Afactsolve = <built-in method solve of SuperLU object at 0x167b25590>
        M          = 2
        N          = 2
        b          = <2x2 sparse array of type '<class 'numpy.float64'>'
	with 2 stored elements (1 diagonals) in DIAgonal format>
        b_is_sparse = True
        b_is_vector = False
        is_pydata_sparse = False
        permc_spec = None
        pydata_sparse_cls = None
        result_dtype = dtype('float64')
        use_umfpack = False

[mtsokol]

The PR is ready from my side!

[rgommers]

This doesn't look correct to me. It's possible I am missing something, but what it looks like is:

1. _ident_like is changed to always return a scipy.sparse matrix
2. in sparse.linalg.inv, that changes the type of I
3. spsolve is now called from within inv with one pydata-sparse array and one scipy.sparse matrix (unsupported in principle!)
4. the handling within spsolve is asymmetric in the two input arrays (only the second input is checked for being pydata-sparse-class), so this goes wrong

In principle the rationale for this change seems good: do the conversions to-from scipy.sparse instances only once. However, it's a little dangerous to do things like this when private functions within scipy.sparse call again public functions from the same module. I think the robust way to do this is to separate public from private better. And even then: mixed-type inputs should never happen.

[mtsokol]

@rgommers Thanks for feedback! After checking it once again I decided to revert _ident_like changes - as you pointed out, to avoid mixing input types. Instead I kept top-level public API changes in expm and expm_multiply to support pydata/sparse there. I changed the title of the PR. Right now:

• _ident_like creates identity matrix with the same type as input,
• expm and expm_multiply support pydata/sparse by dispatching to scipy sparse array and converting back before return (in expm),
• There's no mixing types in spsolve and inv.