Adding some context to MRocklins answer:
A least-squares solution can be written in two ways [1], use the one for which the (to be) inverted matrix is smallest. Dask already has it implemented: dask.array.linalg.lstsq. Linear solvers have (slightly sub-)cubic complexity, so increases in the smallest dimension of A will cause you to hit a wall at some point. People then start using random linear algebra or gradient descend (see also the MIT opencourseware linear algebra series from Gilbert Strang, also available on youtube), for which you probably would have to write your own code (does not hurt to look around though, dask is a nice toolbox for this problem ;). Dask.array.reduction is probably a good function to use if you have to do it yourself.
[1]
x = (A^T A)^-1 A^T d
x = A^T (A A^T)^-1 d