Generative Ensemble Regression: Learning Particle Dynamics from Observations of Ensembles with Physics-Informed Deep Generative Models

For non-interacting particle systems, generating particle trajectories is relatively straightforward by directly applying the discretization of governing SODE or ODE. For example, if the particle trajectories are diffusion processes or Lévy processes, we can use the following forward Euler scheme:

where $\Delta t$ is the time step, $\bm{\xi}_{i}$ and $\bm{\zeta}_{\alpha,i}$ are i.i.d. standard Gaussian random variables and $\alpha$-stable random variables, respectively. We could represent $\bm{\mu}_{t}$ and $\bm{\sigma}_{t}$ with neural networks if they are unknown, or represent the unknown parameters with trainable variables if we know their parameterized form.

Our target is to tune the trainable variables in the generative model, including the parameters in $G$ and those for parameterizing $\bm{\mu}_{t}$ and $\bm{\sigma}_{t}$, so that the generated marginal distribution $\tilde{\rho}_{t_{i}}$ fits the data $\mathcal{D}_{i}$ for each $i$. We thus need to define a distance function $\mathsf{d}(\cdot,\cdot)$ to measure the difference between the two input distributions, which can be estimated from samples drawn from the two distributions. Consequently the loss function in non-interacting particle systems is defined as:

where $\hat{\rho}_{\mathcal{D}_{i}}$ is the empirical distribution induced from the sample set $\mathcal{D}_{i}$. We will refer to it as the distribution loss.

There could be many ways to define $\mathsf{d}$, including Wasserstein distances, maximum mean discrepancy, etc. In this paper we use two approaches to define $\mathsf{d}$.

First we choose the squared sliced Wasserstein-2 (SW) distance [26, 3] as the function $\mathsf{d}$:

where $W_{2}$ is the Wasserstein-2 distance, and ${\pi_{\bm{e}}}_{\#}\mu$ is the one dimensional distribution induced by projecting $\mu$ onto the direction $\bm{e}$, defined by

similarly for ${\pi_{\bm{e}}}_{\#}\nu$. $\mathcal{H}^{d-1}$ is the uniform Hausdorff measure on the sphere $\mathbb{S}^{d-1}$. In short, the squared sliced Wasserstein-2 distance is the expectation of the squared Wasserstein-2 distance between the two input measures projected onto uniformly random directions. The sliced Wasserstein-2 distance is exactly the Wasserstein-2 distance for one-dimensional distributions, but is easier to calculate for higher dimensional distributions. We present the details of the estimation in the Supplementary Information section S1.

We also use GANs to obtain $\mathsf{d}$. The generative model we introduced above can generate “fake” samples $\tilde{\bm{X}}_{t_{i}}$ at $\{t_{i}\}_{i=1}^{n}$, and for each $i$, we use a discriminator $D_{i}$ to discriminate generated samples $\tilde{\bm{X}}_{t_{i}}$ and real samples $\bm{X}_{t_{i}}$ from $\mathcal{D}_{i}$. The adversarial loss given by $D_{i}$ can act as a metric of the difference between $\tilde{\rho}_{t_{i}}$ and $\hat{\rho}_{\mathcal{D}_{i}}$. In particular, we use WGAN-GP [5] as our version of GANs in our paper, with

and the loss function for each discriminator $D_{i}$ is defined as

where $\rho_{\hat{\bm{x}}_{i}}$ is the distribution generated by uniform sampling on straight lines between pairs of points sampled from $\tilde{\rho}_{t_{i}}$ and $\hat{\rho}_{\mathcal{D}_{i}}$, and $\lambda=0.1$ is the gradient penalty coefficient. Here, $\mathsf{d}(\tilde{\rho}_{t_{i}},\hat{\rho}_{\mathcal{D}_{i}})$ can be mathematically interpreted as the Wasserstein-1 distance between the two input distributions. WGAN-GP is computationally more expensive than the sliced Wasserstein distance since we need to train the generative model and the discriminators iteratively, but it is more scalable to high dimensional problems, for which we made a comparison in the supplementary information section S1.

The straightforward discretization of the governing Equations 3 cannot be directly applied to generate “fake” particle trajectories in interacting particle systems, especially for those with strong nonlocal interactions, since the computational cost for one time step would be $O(N^{2})$, where $N$ is the number of particles required so that the system is close to the mean field limit, which makes the learning almost intractable. Instead, we propose to first employ a neural network $\tilde{\bm{\mu}}$ as a surrogate model of the velocity field $\bm{\mu_{t}}$ in the spatial-temporal domain to generate trajectories, and then apply an additional penalty in the loss function to enforce the velocity field $\tilde{\bm{\mu}}$ to be consistent with the Equations 3. In this paper we use the forward Euler scheme to generate trajectories:

where $\Delta t$ is the time step. By employing the surrogate model for velocity, the computational cost for generating the particle trajectories is linear, instead of quadratic, with the number of particles.

To enforce $\tilde{\bm{\mu}}$ to be consistent with our knowledge of the dynamics of Equations 3, we first randomly generate $K$ particle trajectories with Equation 12, denoted as $\{\{\tilde{\bm{X}}^{(k)}_{i\Delta t}\}_{i=0}^{I}\}_{k=1}^{K}$, and then calculate the “forces” applied on these $k$ particles by another $M$ randomly generated particles $\{\{\tilde{\bm{Y}}^{(m)}_{i\Delta t}\}_{i=0}^{I}\}_{m=1}^{M}$, using the formula of interactions in Equations 3, where we replace the unknown parameter $\alpha$ with a trainable variable. Meanwhile, from the velocity field $\tilde{\bm{\mu}}$ we can directly calculate the accelerations for these $k$ particles using a material derivative. The forces and accelerations should be consistent with each other. We can thus define an $L_{2}$ loss:

and we name it as the Newton loss. The time span $[0,I\Delta t]$ should cover the time for the latest snapshot $t_{n}$, and in this paper we simply set $I$ as the ceiling of $t_{n}/\Delta t$. To reduce computational cost, the average over $I+1$ time steps in $L_{\text{Newton}}$ can also be approximated by mini-batch, i.e., taking average over $B$ random time steps in each training step.

For one time step, the computational cost is $O(KM)$ for calculating the force terms and $O(K)$ for the acceleration terms. Note that $M$ should be larger than or equal to $N$, but $K$ can be much less than $N$, thus the total computational cost for $L_{\text{Newton}}$ can be much less than $O(N^{2})$, and this makes the learning tractable.

In the end, the loss for interacting particle systems will be a combination of the Newton loss and the distribution loss:

where the weight $\eta$ is set as 1 in this paper.

In order to provide flexibility and adapt the method to different problems, we can also modify the generated distribution and real data before feeding them to the distance function d and calculate the distribution loss. We present some examples here.

In some cases $\rho_{t}$ may have heavy tails, e.g., when the particle trajectories correspond to a Lévy process. The heavy tails could spoil the training, since the rare outlier samples could dominate the loss function. We can choose a suitable bounded map $h:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ to preprocess both the generated samples and the real data so that the heavy tails are removed. If the observations of the particle coordinates are noisy, we can also perturb the generated samples to add artificial noise. By scaling the artificial noise with trainable variables, the size of the noise in observations can also be learned. If we can only make observations of particles in a specific domain, i.e., the observations are truncated, we will also filter the generated samples using a corresponding mask so that the effective domain for the generated samples and observations are the same.

In this section, we test our method on the 1D problems using the sliced Wasserstein distance as $\mathsf{d}$. We first consider the SODE with Brownian motion and then with the $\alpha$-stable symmetric Lévy process as the stochastic term:

where $\alpha=1.5$, with $\rho_{0}=\mathcal{N}(0,0.04)$.

Note that we have no knowledge of $\rho_{0}$ in learning. Also, we assume we know that the stochastic term has a constant coefficient but we need to infer it; here the ground truth is $1.0$. To represent the unknown coefficient, we use a trainable variable rectified by a softplus function $\text{softplus}(x)=\ln(1+e^{x})$, which ensures positivity. As for the drift $\mu(x)=x-x^{3}$, we consider the following two cases for both SODE problems. In case 1, we know that the drift $\mu(x)$ is a cubic polynomial of $x$. In this case, we use a cubic polynomial $a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}$ to parameterize $\mu(x)$, and want to infer the four coefficients $a_{0},a_{1},a_{2}$ and $a_{3}$. In case 2, we only know that the drift $\mu(x)$ is a function of $x$ and hence we use a neural network to parameterize $\mu(x)$.

For the SODE problem with Brownian motion, we first prepare a pool of $10^{5}$ sample paths, then independently draw $10,000$ samples at $t=0.2,0.5,1.0$ from the pool as our training data. The results for both cases of drift parameterization are illustrated in Figure 3(a). In the figure, all the densities are estimated using Gaussian kernel density estimation with Scott’s rule, and the inferred densities and the reference densities come from $10^{6}$ samples produced by the generative model or simulation.

Both cases of drift parameterization provide a good inference of the diffusion coefficient, with an error less than $7\%$ after $2\times 10^{5}$ training steps, in all the runs. When using the cubic polynomial parameterization, the inferred drift fits well with the ground truth, with the relative $L_{2}$ error about $3\%$ in $[-3,3]$ averaged over three runs. The inferred drift using the neural network only fits the ground truth in the region between -1.5 and 1.5. This is reasonable since the particles mainly concentrate in this region, and we can hardly learn the drift outside this region, where the training data are sparse. However, we note that such an inference of drift is sufficiently good for an accurate time extrapolation of the particle distribution at $t=5$.

One interesting observation is that in case 1 the inferred densities are more accurate than those estimated from the training data. This is because our knowledge of the governing SODE bridges the limited samples in three snapshots, and as we infer the density at, e.g., $t=0.2$, we are not only utilizing the data at $t=0.2$ but also implicitly leveraging the data at $t=0.5$ and $1.0$. We can make an analogy in the context of linear regression: multiple noisy data are helpful for inferring the hidden ground truth, since they are bridged by the regressed linear function. We present a more detailed discussion on this topic in the supplementary information in section S2.

For the SODE problem with Lévy process, we prepare $1.5\times 10^{5}$ sample paths, then independently draw $10,000$ samples within the region $[-1000,1000]$ at $t=0.2,0.5,1.0$ from the pool as our training data. To prevent instability during the training, we apply double precision and clip the generated $\alpha$-stable random variable in Equation 2 between $-100$ and $100$.

We first test our method without pre-processing the samples as in the problem with Brownian motion. As we can see in the second row of Figure 3(b), the inferences are far away from the ground truth. This is due to the heavy tail of $\rho_{t}$ in the Lévy process: some samples far away from $0$, although rare, could dominate the loss function and spoil the training. To deal with this problem, we then apply a bounded map $h(x)=2\tanh(x/2)$ to all the generated and real samples before feeding them to $\mathsf{d}$. The results are shown in the first row of Figure 3(b) where we can see the inferences are much improved. In case 2 where the drift is parameterized by neural networks, the inferred drift outside $[-1.5,1.5]$ is better than that in the problem with Brownian motion. This is because the samples are more scattered in the Lévy case.

In this section, we test our method on 2D problems using the sliced Wasserstein distance as $\mathsf{d}$, with various scenarios of observations. We consider the following 2D SODE:

where

and

The parameters are set as $a_{0}=a_{1}=s_{0}=s_{1}=s_{2}=1.0$, $a_{2}=a_{3}=-0.5$. We set the initial distribution as $\rho_{0}=\mathcal{N}(0,\bm{I}_{2})$. We assume we know that the diffusion coefficient is a constant lower triangular matrix but we need to infer the three unknown parameters $s_{0},s_{1},s_{2}$. In particular, we use $(\text{softplus}(\tilde{s}_{0}),\tilde{s}_{1},\text{softplus}(\tilde{s}_{2}))$ to approximate $(s_{0},s_{1},s_{2})$, where $\tilde{s}_{0},\tilde{s}_{1},\tilde{s}_{2}$ are three trainable variables. Similar as in the 1D problems, we consider the following two cases for the drift. In case 1, we know the form of $\bm{\mu}$ and $\varphi$ in Equation 17 and 18 but we need to infer $a_{0},a_{1},a_{2}$ and $a_{3}$ ($(a_{0},a_{1})$ and $(a_{2},a_{3})$ are exchangeable). In case 2, we only know that $\bm{\mu}$ is a gradient of $\varphi$ in Equation 17, but we have no knowledge of $\varphi$. In this case, we use a neural network to parameterize $\varphi$.

For the training data, we prepare $10^{5}$ sample paths, and consider the following various scenarios of observations at $t=0,0.1,0.2,0.3,0.5,0.7,1.0$. The data are visualized in the supplementary information section S3.

• •

  Scenario 1: we assume that our observations are ideal, i.e., we make accurate observations of all the particle coordinates as the training data.

• •

  Scenario 2: we assume that our observations are noisy. Specifically, we make observations of all the particle coordinates, but each coordinate is perturbed by an i.i.d. random noise $\mathcal{N}(\bm{0},e^{2}\bm{I}_{2})$, where $e$ is set as $0.2$ but also need to be inferred during learning.

• •

  Scenario 3: we assume that our observations are truncated. Specifically, we make observations of the particle coordinates in $\Omega=(\infty,0.5)\times\mathbb{R}$, with the particles outside of $\Omega$ dropped.

For the first scenario, the inferred drift and diffusion are illustrated in Figure 4(a). In case 1, all the inferred drift and diffusion parameters approach the ground truth during the training, with relative error less than $0.3\%$ for each $a_{i}$ and less than $4\%$ for each $\sigma_{i}$ after $2\times 10^{5}$ training steps. In case 2, the inferred diffusion parameters are also close to the ground truth, with error comparable with that in case 1. In the region where the training data are dense (i.e., the density estimated from all the training data is no less than 0.05), the inferred drift field has a relative $L_{2}$ error about $13\%$. The inference is much worse in other regions, and the reason is the same as in the 1D problem: the neural network can hardly learn the drift field where the training data are sparse.

For the second and third scenarios, we apply the technique of perturbing and filtering the generated samples, respectively. In Figure 4(b) and Figure 4(c) we show the results of case 1 with the drift parameterized by four parameters. After $2\times 10^{5}$ training steps, all the parameters converge to ground truth. For scenario 2, the error is less than $0.5\%$ for each drift parameter $a_{i}$, less than $3\%$ for each diffusion parameter $s_{i}$, and about $5\%$ for the noise parameter $e$. For scenario 3, the error is less than $0.2\%$ for each drift parameter $a_{i}$, less than $1\%$ for each diffusion parameter $s_{i}$. We remark that we failed to learn the drift with neural network parameterization in the second and third scenario, suggesting a room for further improvements.

In this section, we test our method on higher dimensional non-interacting particle systems. Note that the sliced Wasserstein distance does not perform well in high dimensional problems, and we thus switch to the WGAN-GP to provide $\mathsf{d}$. The comparison between the sliced Wasserstein distance and WGAN-GP in high dimensional problems is presented in the supplementary information section S1.

We consider a $d$-dimensional SODE:

where

for $\mu^{(i)}$ as the $i$-th component of $\bm{\mu}\in\mathbb{R}^{d}$, and $X_{t}^{(i)}$ is the $i$-th component of $\bm{X}_{t}\in\mathbb{R}^{d}$. We set the diffusion coefficient matrix as

where the $2d-1$ nonzero entries $\{s_{i}\}_{i=1}^{d}$ and $\{s^{\prime}_{i}\}_{i=2}^{d}$ are set as 1. The initial distribution is $\rho_{0}=\mathcal{N}(\bm{0},0.04\bm{I}_{d})$. Due to the non-diagonal diffusion coefficient matrix, the motion in different dimensions are coupled.

We prepare $10^{5}$ sample paths and observe all the particle positions at $t=0.2,0.5,1.0$ as our training data. We use a cubic polynomial with four variables to parameterize $\mu^{(i)}$ as a function of $X_{t}^{(i)}$ for each $i$, while using $2d-1$ trainable variables to learn the nonzero entries in the diffusion coefficient matrix, with the diagonal entries rectified by a softplus function for positivity. In other words, we use $6d-1$ variables to parameterize the drift and diffusion for the $d$-dimensional SODE. The results are shown in Figure 5. Even for the 100-dimensional problem, after $6\times 10^{5}$ training steps, the average error of the diffusion coefficients is about $0.04$, and the average relative $L_{2}$ error of the drift in the interval $[-3,3]$ in each dimension is about $14\%$.

It may appear surprising that we can solve a 100-dimensional problem with only $10^{5}$ samples since usually an exponentially large number of samples is required to describe the distribution. However, we remark that the difficulty of the learning task is significantly reduced since our partial knowledge, i.e., the form of the SODE, is encoded into the generator.

In this section we consider the 1D and 2D cases of the interacting particle system with governing equations 3. To prepare the training data, we follow [23] and set the initial condition, including the density and velocity as

for the 1D case and

for the 2D case, where $\mathbbm{1}$ is the indicator function. Using the Velocity-Verlet method with time step 0.01, we perform simulations with 1024 and 9976 particles for the 1D and 2D cases, respectively, and generate data at $t=0.5,0.6,...,2.0$, i.e., 16 snapshots in total. The input radius $r$ of the influence function $\phi(r)$ is clipped to be at least $r_{\text{min}}=0.01$, both in simulation and learning, to avoid the singularity. While using the same number of particles and make observations at the same time instants as in [23], our method does not require knowledge of the initial condition or the velocity field in the data snapshots, as opposed to their method based on Bayesian optimization to infer $\alpha$.

We apply the sliced Wasserstein distance for the distribution loss, with the batch size equal to the number of particles in training data. When calculating the distance at each time instant, the generated and real distributions are normalized with the mean and standard deviation of the real distribution. For the Newton loss, we set $K=16$, $M=10000$ and $B=10$. We use $2\times\text{sigmoid}(\beta)$ to represent the inferred $\alpha$ so that the inference is bounded by 0 and 2, and the variable $\beta$ is initialized as 0.

The results are visualized in Figure 6. For the 1D case, the inferred $\alpha$ is 0.492 with ground truth 0.5 in the end of training. As a comparison, the inference is 0.480 for the 1D case in [23]. Our inferred velocity field also matches well with the reference ground truth at different time instants, even at $t=0$ and $0.25$ when we have no observations at all. This is because we incorporated the knowledge of dynamics in the learning system. For the 2D case, while the inferred velocity field also matches well with the ground truth, the inferred $\alpha$ is 0.469 with ground truth 0.5 in the end of training, i.e., about 6% error, much larger than in the 1D case. We attribute the error to the fact that the singularity problem for $\phi(r)$, where the order for $r$ is dimension-dependent, is more severe in the 2D case. In supplementary information section S4, we show that as we increase $r_{\min}$ from 0.01 to 0.1, the error is reduced from 6% to 2%. Alternatively, in the 2D case we also try to directly use the data at $t=0.5$ as the starting coordinates to generate trajectories for $t\geq 0.5$ ($M$ is also reduced to the data size $9976$), so that the initial generator $G$ is removed from the generative model. By doing so, the number of trainable variables is reduced and the learning becomes easier. The inferred $\alpha$ is 0.493 in the end of training. As a comparison, the inference is 0.513 for the 2D case in [23]. However, we remark that while this strategy helps to infer $\alpha$ without extra data, we cannot infer the velocity or density for $0\leq t<0.5$.

This work was supported by the PhILMS grant DE-SC0019453 and by the OSD/AFOSR MURI Grant FA9550-20-1-0358. We would like to thank Prof. Hui Wang and Ms. Tingwei Meng for carefully checking the proof of our theorem. We also want to thank Dr. Zhongqiang Zhang for helpful discussions.

The squared sliced Wasserstein-2 distance is the expectation of the squared Wasserstein-2 distance between the two input measures projected onto uniformly random directions. To estimate $SW_{2}(\mu,\nu)$ from samples of $\mu$ and $\nu$, we use the following process introduced in [3].

• •

  Draw samples independently from $\mu$ and $\nu$ with batch size $b$, denoted as $\mathcal{U}$ and $\mathcal{V}$.

• •

  Uniformly sample $m$ projection directions $\{\bm{e}_{j}\}_{j=1}^{m}$ in $\mathbb{R}^{d}$. In this paper we set $m=1000$.

• •

  For each random direction $\bm{e}_{j}$, project and sort the samples in $\mathcal{U}$ and $\mathcal{V}$ in the direction of $\bm{e}_{j}$, getting $\{u_{i,j}\}_{i=1}^{b}$ and $\{v_{i,j}\}_{i=1}^{b}$, where $u_{i,j}\leq u_{i+1,j}$ and $v_{i,j}\leq v_{i+1,j}$ for $i=1,2...b-1$. Calculate $L_{j}=\sum_{i=1}^{b}(u_{i,j}-v_{i,j})^{2}/b$.

• •

  Calculate $L=\sum_{j=1}^{m}L_{j}/m$ as the estimation of squared $SW_{2}(\mu,\nu)$.

Compared with GANs, the sliced Wasserstein distance does not need discriminators, and is more robust than WGAN-GP in low dimensional problems. Take the following 1D problem as an example. We consider the SODE:

with $\rho_{0}=\mathcal{N}(-0.5,0.5)$. We set $a=b=1$ so that the exact solution is

We have $10,000$ samples at $t=0.5$ and $1.5$, respectively, as the training data, and wish to infer the constant drift and diffusion coefficients. The input noise to the generator is uniform noise from -1 to 1. The results are shown in Figure 8.

We can clearly see oscillations of the inferred drift coefficient for WGAN-GP. This can be attributed to the two-player game between the generator and discriminator, which was also reported in [30].

We did not observe significant oscillations using WGAN-GP in higher dimensional problems. Indeed, WGAN-GP outperforms the sliced Wasserstein distance in high dimensional problems. Let us consider the following simple SODE problem as an example, where the motions are uncoupled between dimensions:

where $X_{t}^{(i)}$ is the $i$-th component of $\bm{X}_{t}\in\mathbb{R}^{d}$, with $\rho_{0}=\mathcal{N}(\bm{0},0.04\bm{I}_{d})$. We prepare $10^{5}$ sample paths and observe all the particle positions at $t=0.2,0.5,1.0$ as our training data. We use a cubic polynomial $a^{(i)}_{0}+a^{(i)}_{1}x+a^{(i)}_{2}x^{2}+a^{(i)}_{3}x^{3}$ to parameterize the $i$-th component of the drift. We use another trainable variable rectified by a softplus function to represent the diffusion coefficient $s^{(i)}=1$ in each direction. In total, $5d$ variables are used to parameterize the $d$-dimensional SODE. The results are shown in Figure 9. While the SW distance works for 2D problems, it does not scale well to high dimensional problems for which WGAN-GP gives much better inferences.

In this section, we wish to show that our method can make use of the SODE and multiple snapshots to reduce the error of density estimation from limited samples. We consider the one-dimensional SODE:

with $\rho_{0}=0.5\mathcal{N}(-0.5,0.3^{2})+0.5\mathcal{N}(0.5,0.3^{2})$.

We test the following two cases of training data sets. In both cases we have 1000 samples at $t=0.05,0.1,0.15,0.2,0.25$, but the particle trajectories are different.

• •

  Case 1: We make observations from different sets of trajectories at different time instants.

• •

  Case 2: We make observations from the same set of trajectories at different time instants.

We assume that we know Equation 28, but have no knowledge of $\rho_{0}$. We could make the analogy of performing linear regression, where we know the slope but do not know the intercept. In Figure 10 we show the comparison between the inferred densities and the densities estimated directly from data. All the densities are estimated using Gaussian kernel density estimation. The inferred densities and the ground truth come from $10^{5}$ samples produced by the generative model or simulation, with bandwidth decided by the Scott’s rule. To remove the effect of bandwidth selection in the density estimation directly from data, we perform a grid search of the optimal bandwidth factor via the $L_{2}$ error against ground truth, from 0.01 to 0.5 with grid size 0.01 (Scott’s rule suggests about 0.25).

The grid search strategy is actually infeasible in practice since we don’t know the ground truth, but it should perform better than any bandwidth selector. Despite that, in Figure 10(a) we can clearly see that in case 1 our inferred densities with naive Scott’s rule significantly outperform the ones estimated directly from training data with grid search. This cannot be attributed to the number of samples, since in case 2 the inferred densities with our method cannot perform better, as shown in Figure 10(b).

The different performances in two cases are reasonable. In case 1, our method can utilize the SODE and the observations at multiple time instants to reduce the error from limited samples. As we infer the density at $t=0.05$, we are not only utilizing the data at $t=0.05$ but also implicitly taking advantage of the data at later time instants. In case 2, since the observations come from the same sample paths, the observations at later time instants cannot provide additional information considering that the SODE is already known. Let us again make an analogy in the context of linear regression: if the observations have independent noise, multiple observations will be more helpful than a single observation. However, if the observations have the same noise, multiple observations cannot help us more than a single observation, if we already know the slope.

In Figure 11 we visualize the training data in 2D problems with various scenarios of observations.

In this section we study the effect of various clipping radius $r_{\text{min}}$ from 0.01 to 0.1, both in the simulation for data generation and the learning algorithm. We remark that changing $r_{\text{min}}$ in this range would not make a huge difference to the system behavior, since the dynamics are dominated by non-local interactions. In Figure 12 we show the inferred $\alpha$ during the training. As $r_{\text{min}}$ increase from 0.01 to 0.1, the inferred $\alpha$ in the end of training increases from 0.469 to 0.490, with the ground truth 0.5. This suggests that the error mainly comes from the singularity problem of $\phi(r)$ when $r$ is small.

The proof is based on the weak convergence: convergence of a sequence of probability measures $\{P_{n}\}_{n=1}^{\infty}$ to a probability measure $P$ in the Wasserstein-$p$ metric is equivalent to the weak convergence plus the convergence of $p$-th moments of $\{P_{n}\}_{n=1}^{\infty}$ to $P$ [31]. Here, the weak convergence mean $E_{P_{n}}[f(x)]$ converges to $E_{P}[f(x)]$ for each $f\in\mathcal{F}$, where $\mathcal{F}=\{\text{all bounded and continuous functions}\}$ or $\mathcal{F}=\{\text{all bounded and Lipschitz functions}\}$, which are equivalent[32]. Note that in Theorem 1 we assume that the domain for $X_{t}$ is a compact subset of Euclidean space for each $t$, thus $\|X_{1:T}\|^{p}$ is a bounded and continuous function of $X_{1:T}$. The convergence of $p$-th moment then directly comes from the weak convergence, i.e., we only need to show that $\{P_{n}\}_{n=1}^{\infty}$ converges to $P$ weakly. This is proved in Lemma 4 below.

The notations are inherited from Theorem 1. For simplicity, we will sometimes use $X,Y,Z$ to represent the nodes in sequence, and use $P^{XY}$ to represent the probability measure of $(X,Y)$, etc.