In numpy's version of the least squares computation, passing ordinate values (b) in a two-dimensional array (M, K) results in a 2D array (N, K) for the coefficients and a 1D array (K,) for the residuals, one for each column.

Dask seems to be missing the second case. When passed a 2D array for b, it returns a (1, 1) array for the residuals, instead of the expected (K,). The docstring doesn't discuss it, but it states that the output is (1,). So it's more than a missing implementation.

I think the culprit is the following line:

AFAIU, it should read:

residuals = (residuals ** 2).sum(axis=0, keepdims=True)
# So the output would be either (1,) or (1, K)
# OR
residuals = (residuals ** 2).sum(axis=0, keepdims=(b.ndim == 1))
# So it copies numpy with (1,) or (K,)