Supporting just numpy should be relatively easy. This will also work for `method="blockwise"` by default. We may want to rename `groupby_reduce` to `groupby_agg`? For dask proper, we'll need to use `dask.array.cumreduction` instead of `dask.array.blockwise` + `dask.array.reductions._tree_reduce`