Do you know please how to write this in a matrix form in order to be able to calculate this amount efficiently in numpy? D0 is a (n0,n0) matrix, D is a (n,n) matrix, mu0 is (n0,1) and mu is (n,1). Maybe it is not necessary to write this in a matrix form and there is an easy command with numpy to write 4 loops, I do know know.
$$ \overline{\mathrm{GW}}_2^2(\mathbf D_0, \mathbf D) := \min_{\Gamma \in \overline{\mathcal M}} \sum_{ijk\ell}(\mathbf D_{0ij}- \mathbf D_{k\ell})^2 \mathbf \Gamma_{ik} \mathbf\Gamma_{j\ell} \boldsymbol\mu_{0i}\boldsymbol\mu_{0j}\boldsymbol\mu_k\boldsymbol\mu_\ell. $$