Matrix multiplication inconsistent for dask and sparse: `dask @ sparse` works, `sparse @ dask` fails

What is your issue?

The order of matrix multiplication matters for dask and sparse arrays, but probably shouldn't.

Here is an example:

import numpy as np
from sparse import COO
import xarray as xr

sparse_da = xr.DataArray(COO.from_numpy(np.arange(10)))
dask_da = xr.DataArray(np.arange(100).reshape(10, 10).chunk({"dim_1": 5})

# this works:
dask_da @ sparse_da

# this raises a TypeError for "unsupported types"
sparse_da @ dask_da

# this works as expected
sparse_da @ dask_da.as_numpy()

In the workflow where this is used, the dask array has no chunks along its common dimensions with the sparse array, so it seems like sparse @ chunk should be fine. Also, in this workflow, loading the dask array into memory or making the sparse array dense would use a very large amount of memory.

Start by running the reported example and comparing both matrix-multiplication orders. Trace how xarray DataArrays containing sparse and dask arrays handle each operand order; done means sparse_da @ dask_da works without loading the dask array into memory or densifying the sparse array.