Let $\mu(dx)=\sum_{i=1}^np_i\delta_{x_i}(dx)$ and $\nu(dy)=\rho(y)dy$ be two probability measures on $\mathbb R^d$. Consider the $2-$Wasserstein distance below:
$$W_2(\mu,\nu)^2 \quad := \quad \inf_{\pi\in\Pi(\mu,\nu)}~ \int_{\mathbb R^d\times\mathbb R^d}~ |x-y|^2\pi(dx,dy),$$
where $\Pi(\mu,\nu)$ denotes the collection of couplings $\pi$ of $\mu$ and $\nu$. Obviously $W_2(\mu,\nu)$ depends on the points $x_1,\ldots, x_n$, and we denote $\mu\equiv \mu(x_1,\ldots,x_n)$ and define the energy function
$$E(x_1,\ldots, x_n): = W_2^2\big(\mu(x_1,\ldots, x_n),\nu\big).$$
I am interested in the minimisation of $(x_1,\ldots, x_n)\mapsto E(x_1,\ldots, x_n)$. My first intuition is to compute the gradient of $E$ and to use Newton's gradient method (the domain of definition of $E$ is not convex). I would like to know whether there is any related reference on this problem, and all answers or comments are highly appreciated!