I wonder if it might be possible to make use of NumPy's generalized linear algebra ufuncs in dask.

For example, here's numpy's syntax for solving M*x=b for 1000 different (M, b) pairs:

>>> M = np.random.rand(1000, 10, 10)
>>> b = np.random.rand(1000, 10)
>>> x = np.linalg.solve(M, b)
>>> x.shape
(1000, 10)

To confirm this is doing what we expect:

>>> np.allclose(b, [np.dot(Mi, xi)
                    for (Mi, xi) in zip(M, x)])
True

This seems like something that would naturally fit into dask for an appropriately-chunked array, in order to distribute such computations. Is there a way to accomplish this currently, or is it on the road-map?

Thanks!