Currently the dask array matmul implemenation is similar to the tensordot implemenation (understandably) but lacks a clever trick that tensordot uses to avoid contractions in blockwise (which can result in large pileups of intermediate results).

Tensordot uses only outer products in blockwise (there are no shared indices in the blockwise call) , and then uses a sum call to reduce these. Because sum uses tree reductions, this results in nicer memory usage. Matmul doesn't use this trick.

I briefly tried replacing the blockwise call in matmul with a call to our internal dot function, but various tests failed. It would be good to do one of two things:

1. Find a way to back matmul by tensordot, so that the reduction trick is shared with a common codebase
2. Replicate the reduction trick in matmul

This was originally reported by @eric-czech and @tomwhite in https://github.com/pystatgen/sgkit/issues/375#issuecomment-731060672