Many machine learning problems can be expressed as the optimization of some cost functional over a parametric family of probability distributions. It is often beneficial to solve such optimization problems using natural gradient methods. These methods are invariant to the parametrization of the family, and thus can yield more effective optimization. Unfortunately, computing the natural gradient is challenging as it requires inverting a high dimensional matrix at each iteration. We propose a general framework to approximate the natural gradient for the Wasserstein metric, by leveraging a dual formulation of the metric restricted to a Reproducing Kernel Hilbert Space. Our approach leads to an estimator for gradient direction that can trade-off accuracy and computational cost, with theoretical guarantees. We verify its accuracy on simple examples, and show the advantage of using such an estimator in classification tasks on Cifar10 and Cifar100 empirically.
Thank you for your revision and for the rebuttal. This a  strong submission with insightful angle on natural gradients and with provable guarantees. The authors improved a lot the manuscript and incorporated reviewers feedback. I am increasing my score to 8. 
 
###
Summary of the paper: 
The paper provides a way to estimate the natural Wasserstein gradient using Kernel estimators. The idea is neat and novel. Natural Wasserstein Gradient similar to the so called natural fisher gradients preconditions the gradient using a matrix that uses the local curvature of the manifold of the parametric distribution. 
Authors give variational forms of the Fisher information matrix of an explicit model , using  the variational form of the chi squared or the Fisher Rao divergence. Similarly authors give a variational form of the wasserstein natural gradient . Let theta be the parameter of the parametric implicit model, theta in R^q.  For a descent direction $u$, the variational form is obtained via finding an objective S,  $\sup_{f\in C^c_{\infty}} S(f, u ) = u^{\top}G_{W}u$, where $G_{W}$ is a form of  "Wasserstein information matrix".
Authors then propose to learn the function f in an RKHS and propose to find the descent direction by solving
 $\min_{u} <u, \nabla_{\theta} Loss (p_{\theta})>  + \sup_{f\in RKHS}  S(f, u ) + r(u)- \lambda ||f||^2_{rkhs}$
where r(u) is a quadratic regularizer on u. 
The sup problem has a closed form solution and can be approximated  using Nystrom approximation and randomization on dimensions. The problem in u has also a closed form solution , and one used u as the proxy to the natural  descent. 
Authors under some assumption  show that the  estimated natural W gradient in RKHS  is concentrated around the true one. 
Experiments on synthetic data and in classification on CIFAR 10 and CIFAR 100 shows that the preconditioning of the gradients that the method offers allows faster convergence in both well conditioned and ill conditioned initialization of the weights of the neural network.
 
Review : 
The paper is not easy to follow and the high level intuition how the method works is not well explained. 
It would be easier for the reader, to motivate the natural wasserstein descent from how one defines natural Fisher descent , where one seeks a first order approximation of $KL(p_{\theta},p_{\theta+ \epsilon u})$ as we perturb in the parameter space and this well known that this epsilon $u^{\top}F u$. 
Hence natural Gradient descent is :
$\min _{u} <u, \nabla_{\theta}Loss(p_{\theta})> + KL(p_{\theta},p_{\theta+  u}) \approx min _{u} <u, \nabla_{\theta}Loss(p_{\theta})> + u^{\top}F u$  
Now for the wasserstein distance one has also similarly:
$\min _{u} <u, \nabla_{\theta}Loss(p_{\theta})> + W^2_2(p_{\theta},p_{\theta+  u})$
and it is known that as epsilon goes to zero we have: 
 $W^2_2(p_{\theta},p_{\theta+ \epsilon u})/\epsilon = ||p_{\theta} - p_{\theta+ \epsilon u}||^2_{H^{-1}(p_{\theta})}+o(\epsilon) = \sup_{f} \int f (p_{\theta} - p_{\theta+\epsilon u}) - \frac{1}{2} \mathbb{E}_{p_{\theta}}||\nabla_x f(x)||^2+o(\epsilon) $
Now replacing with the implicit model as epsilon goes to zero we get the expression given in the paper using a simple taylor expansion:
$=   \sup_{f} \int <\nabla_{\theta} h_{\theta}^{\top}\nabla_x f(h_{\theta}),u > d\nu - \frac{1}{2}\mathbb{E} _{p_{\theta}}||\nabla_x f(x)||^2$
in a sense the paper is proposing to linearize $W^2_2$ around the perturbation in the parameter space of the implicit model and this can be done using  $|| .||_{H^{-1}(q)}$ , as pointed and used in many recent works.  then the paper proposes to approximate $|| .||_{H^{-1}(q)}$ in RKHS which was already proposed in Mroueh et al  in Sobolev Descent.  AISTATS 2019. 
We encourage the authors to layout in the beginning the derivations form this point of view which will make the paper easier to digest, the expression in Equation 7 seems mysterious and pulled out of a hat, but it is easier to understand by going to perturbation analysis usually done on KL for Fisher Natural gradient and to do it also here starting from the linearization of $W_2$ with $||.||_{H^{-1}(q)}$  , and how to approximate it in RKHS as it was already proposed in the literature in Mroueh et al Sobolev Descent. 
I read carefully the proofs of Proposition 1, 2, 3. I did not ready full the proofs of the concentration of the estimator , but they seem sensible as they follow usual bounding strategies in this context. 
Questions:
- There is nothing special about the wasserstein natural gradient flow variational form and implicit model, once can apply the same to the variational form of Fisher, that would be probably more efficient? It would be great to baseline this one ? 
-the constraint $\int f(x)p_{\theta}(x)=0$ is not imposed in the kernelized version?
- the method comes disappointing since it seems that the preconditioning that the Wasserstein gradient gives is not enough and $r(u)=u^{\top}D u$ is need where D is diagonal depends on T. Have you tried with $D=Identity$? it might be that the scaling of the gradients is coming only from that $D^{-1}$?
- Can you give timings for computing each gradient update and how it compares to regular SGD or diagonal approximation of Fisher natural gradient?
- Does one need preconditioned gradient if the network was self normalized (like batch norm or spectral norm etc)?   
Overall assessment: 
That is a good theoretical work with provable guarantees. The computational complexity of each gradient estimate is large which makes the method not quite appealing in practice.