Random batch methods (RBM) for interacting particle systems
Shi Jin Lei Li Jian-Guo Liu




1]Shi Jin shijin-m@sjtu.edu.cn
2]Lei Lileili2010@sjtu.edu.cn
3]Jian-Guo Liujliu@phy.duke.edu
[1,2]School of Mathematical Sciences, Institute of Natural Sciences, MOE-LSC, Shanghai Jiao Tong University, Shanghai, 200240, P. R. China.
[3]Department of Mathematics and Department of Physics, Duke University, Durham, NC 27708, USA.


Random Batch Methods (RBM) for interacting particle systems
    [
    
===========================================================





We develop Random Batch Methods for interacting particle systems with large number of particles. These methods use small but random batches for particle interactions, thus the computational cost is reduced from O(N^2) per time step to O(N), for a system with N particles with binary interactions.  On one hand, these methods are efficient Asymptotic-Preserving schemes for the underlying particle systems, allowing N-independent time steps and also capture, in the N →∞ limit, the solution of the mean field limit which are nonlinear Fokker-Planck equations; on the other hand, the stochastic processes generated by the algorithms can also be regarded as new models for the underlying problems. For one of the methods, we give a particle number independent error estimate under some special interactions. Then, we apply these methods to some representative problems in mathematics, physics, social and data sciences, including the Dyson Brownian motion from random matrix theory, Thomson's problem, distribution of wealth, opinion dynamics and clustering. Numerical results show that the methods can capture both the transient solutions and  the global equilibrium in these problems.





§ INTRODUCTION


In natural and social sciences, there are many collective behaviors resulted from a huge number of interrelated individuals. Examples include swarming or synchronization described by the Vicsek model <cit.>, flocking in school of fishes <cit.>, groups of birds <cit.>, chemotaxis of bacteria <cit.>, consensus clusters in opinion dynamics <cit.>, to name a few. 

We are interested to develop random algorithms for these interacting particle systems to reduce computational cost. 
While the ideas also work for second order systems like Langevin dynamics <cit.> and the model in <cit.>,  let us firstly focus on the following first order systems:

    d X^i=-∇ V(X^i) dt+∑_k: k≤ Jα_N,k∑_𝒞⊂{1,…, N}: i∈𝒞, |𝒞|=k K^i(𝒞) dt+σ dB^i, i=1,…, N,

where J is independent of N and α_N,k are some constants. In other words, the interaction acting on each particle is the superposition of the ones for which the number of particles involved is independent of N. Here, X^i∈ℝ^d are the labels for particles, -∇ V(·) is some given external field, K^i(·)'s are the interaction kernels and {B^i}'s are independent standard Brownian motions.  Such systems could be the overdamped limit of Langevin equations, where dX^i comes from the friction term so that X^i is the location of the ith particle and the terms on the right hand side are forces. Of course, X^i can also have other interpretations depending on applications, like the velocity <cit.>, the opinion or wealth (see section <ref>). We will loosely call X^i's the locations of particles, -∇ V the external force and K(·) the interacting forces for convenience in this paper, unless there are clear meanings. If σ=0, there is no diffusion effect. If σ>0, we have diffusion, and the equation is a stochastic differential equation (SDE) in Itô's sense <cit.>. We will mostly focus on J=2, so that the systems are of binary interactions (i.e. each particle interacts with others separately). System (<ref>) in the case of binary interaction can be written as

    d X^i=-∇ V(X^i) dt+1/N-1∑_j: j≠ i K(X^i-X^j) dt+σ dB^i, i=1,…, N.

We will present our algorithms and analysis for binary interactions (system (<ref>)) and consider that the initial values X_i(0)=:X_0^i are drawn independently from the same distribution:

    X_0^i ∼ν, i=1,…, N.

 Extensions to (<ref>) are straightforward, though more involved (see section <ref> and Remark <ref>).
Although there are examples of many body interactions like Kac interaction in spin glasses <cit.> for (<ref>), we remark that the binary interaction is much more important and ubiquitous in nature, like Coulomb's interaction between electrons and nuclei, due to the superposition principles for forces.  As another remark, we have assumed the additive noise for simplicity to illustrate our methodology. However, the discussion in this paper also applies for multiplicative noise (i.e. σ depends on x). (We in fact have one such example in Section <ref>.)

There may be two things we care about regarding (<ref>) and (<ref>), depending on applications. The first thing is the dynamics, for example, in the opinion dynamics, we would like to know how the consensus of opinions are developed.  Another thing is the equilibrium distribution of the N particles. If N is large, simulation of (<ref>) and (<ref>) is very expensive since for each time step during the evolution, the computational cost is of O(N^J).  One possible method for studying large and complex systems where individuals interact with each other is the mean field approximation <cit.>. In this approximation, the effect of surrounding particles is approximated by a consistent averaged force field so that one obtains a one body problem.
For example, regarding (<ref>) that we will focus on, in the N→∞ regime, the distribution μ of the particles formally satisfies the following nonlinear Fokker-Planck equation (see <cit.>)

    ∂_tμ=∇·(∇ V(x)μ)-∇·(μ K*μ)+1/2σ^2Δμ.

It is expected that μ is close to the empirical measure for (<ref>)

    μ_N(t):=1/N∑_i=1^N δ(x-X^i(t)),

which is a random measure. In fact, under certain assumptions on V and K, it can be shown that the uniform mean field limit holds (<cit.>).
Both the marginal distribution of X^1 in (<ref>), denoted as μ_N^(1), and (<ref>) are close to μ.

On the other hand, one often cares more about the mean field equations like (<ref>) and its invariant measure π, but the mean field equation and π are hard or expensive to solve or compute. Using the mean field approximation, one can generate some artificial particle systems of the form (<ref>) or (<ref>). In this sense, the interacting particle systems (<ref>) and (<ref>) are Monte-Carlo particle methods for the mean field equation and π. Clearly, when N (the number of numerical particles) is large, one still wants to reduce the computational cost.
Hence, no matter whether studying system (<ref>) (and (<ref>) ) is due to its own interest or due to numerical simulation of the mean field equation, it is highly desirable to design some efficient algorithms to solve the particle system (<ref>). 

Nowadays, in the era of big data, many stochastic algorithms have been developed to reduce the computational cost while keep certain accuracy.  Hence, one may borrow some ideas from these areas for the physical problems. The stochastic gradient descent (SGD) algorithm is developed to reduce the computational cost and for better exploring the high dimensional parameter space <cit.> for supervised learning <cit.>. In SGD, a small batch of samples are chosen each time to form the noisy gradient. A similar algorithm is the stochastic coordinate descent where only a few coordinates are updated each time <cit.>.  The idea of random batches also appears in the stochastic gradient Langevin dynamics (a Markov Chain Monte Carlo method (MCMC)), which can be applied for Bayesian inference <cit.>. Pretty much like SGD, the method uses random batch of data to update the parameters one needs to estimate. Besides, many random algorithms also appear for matrices, such as constructing approximate matrix decompositions by randomized Singular Value Decompositions <cit.> and the fast Monte Carlo algorithm for matrix-matrix product in low rank problems <cit.>. 

Motivated by the stochastic algorithms in machine learning and statistics, we develop Random Batch Methods (RBM) for interacting particle systems (<ref>) and (<ref>). The idea is quite simple: for a small duration of time, one randomly divides the particles into small batches and the interactions are only turned on inside each batch (see Section <ref> for details). Numerical experiments in later sections show that RBMs work quite well for a wide range of applications with binary interactions from physical, biological to data sciences.  These algorithms reduce the computational cost per time step from O(N^2) to O(N) for binary interactions (<ref>) (O(N^J) to O(N) for (<ref>)).  They not only recover the equilibrium well in variaous problems, but also approximate dynamics of measures with very singular interaction kernels (see for example section <ref>). Moreover in section <ref>, under suitable conditions on the external potential, interaction force and batch size, we prove for one of the RBMs that the numerical error has a bound that only depends on the time step,  independent of N and time. The key for the proof is Lemma <ref>, which guarantees that on average the random force is consistent with the full interaction.


The methodology of randomly choosing some objects has also been widely used in kinetic theory community. 
The multi-level Monte Carlo (MLMC) method <cit.> runs simulations with different number of particles for McKean-Vlasov equation using the idea of subsampling. The authors used smaller numbers of particles for repeated simulations with a subset of random variables (such as the Brownian motions) for the finest simulation to reduce the computational cost, while keeping good enough MLMC estimator. Different from <cit.>, RBMs use the idea of subsampling in computing interacting forces for one single simulation to approximate the distribution of particles.
The Direct Simulation Monte Carlo (DSMC) method proposed by Bird (<cit.>) uses randomly chosen simulation molecules to solve the Boltzmann equation with finite Knudsen number (see the work by Wagner <cit.> for the proof of convergence).  Moreover, Nanbu's simulation method directly derived from the Boltzmann equation has gained success and been proved to converge <cit.>.
In <cit.>, biological swarm models with random pair interactions have been studied and proven to converge to Kac-Boltzmann equation. In the reverse way to <cit.>,  the work by Albi and Pareschi (<cit.>) uses stochastic binary interaction algorithms for flocking and swarming dynamics. The algorithms in <cit.> are intrinsically doing similar things as our RBMs do; their Algorithm 4.3 is particularly alike RBM-1 (see section <ref>) with p=2. See a more recent application by Carrillo et al <cit.>. In this sense, our RBMs are generalizations of symmetric Nanbu (Algorithm 4.3) in <cit.>.
However, the numerical particle system in <cit.> was motivated by the mean-field limit, in a way similar to the DSMC methods of binary collisions for the Boltzmann equation.   It was aimed at solving the mean-field equation (corresponding to the right vertical line in Fig. <ref>).  RBMs are motivated by the idea of mini-batch methods from machine learning. They are Monte-Carlo type approximations directly to the (physical) particle systems (<ref>). They  correspond to the left horizontal line in Fig. <ref>. Moreover, one of the RBMs will be proved (under special conditions) to converge, when time steps go to zero, to the particle system (<ref>).  The construction of RBMs, as well as the proof of their convergence,  are obtained  without the knowledge of the mean-field limit.  Of course when N→∞ its density measure also converges to the mean-field limit equation, as will be proven using the Wasserstein distance. Therefore, the RBMs are a class of Asymptotic-Preserving schemes for particle systems (<ref>) or (<ref>), in the sense of 
<cit.> (see Fig. <ref>). Lastly, the idea of turning on interactions inside the batches can be easily extended. For example, one can change the batch size to adjust the noise levels. Moreover, RBMs can also be extended to (<ref>) involving more complex interactions.

Let us remark that RBMs can be straightforwardly applied to the second order particle systems. In Appendix <ref>, we establish similar error estimates for a Hamiltonian system, valid in finite time, and perform the numerical verification in section <ref>.

The rest of the paper is organized as follows. In section <ref>, we propose the Random Batch Method without replacement (RBM-1) and with replacement (RBM-r), and give some discussions. In section <ref>, we obtain an error estimate, which only depends on time step but not on N and time,  of the marginals for RBM-1 in the Wasserstein distance under some assumptions for the external and interacting forces. Though performed for binary interactions, the analysis should work for (<ref>) as well. In section <ref>, we do numerical tests to verify the theoretic results in section <ref> and use the Dyson Brownian motion to compare the RBMs.  In particular, in section <ref>, we compute the law of the Dyson Brownian motions from random matrix theory where the interaction kernel is singular. Our method can capture the evolution of distribution and the equilibrium semicircle law, and the two algorithms (RBM-1 and RBM-r) give comparable results. In section <ref>, we use the RBMs to solve some interesting application problems with binary interactions.  In section <ref>, we focus on Thomson's problems where we solve the dynamics on the sphere.  In section <ref>, we design stochastic dynamics for evolution of wealth and opinions. Lastly, we apply our randomized algorithms for efficient clusters in section <ref>. The paper is concluded in section <ref>.





§ THE RANDOM BATCH METHODS


In this section, we propose the RBMs by using random batches for the summation of the interacting force in (<ref>). The extensions to (<ref>) should be similar but are more involved, and we will give some discussion about this in section <ref>.
For the setup, we pick a short duration of time τ and consider the discrete time

    t_m=mτ.

Suppose we compute up to time T and the number of iteration for the stochastic algorithm is

    N_T=⌈T/τ⌉.

Clearly to simulate the ODE system (<ref>) directly, the complexity is 
O(N_T N^2). If N is large, this is expensive. Motivated by the stochastic algorithms in the machine learning community, we will use a randomized strategy.  




 §.§ The first algorithm (RBM-1)


A natural idea is that at each time t_m, we divide the N=np particles into n small batches with size p (p≪ N, often p=2) randomly, denoted by 𝒞_q, q=1,…, n, and then interact particles within each batch (For convenience, we have assumed that p divides N. In general, the last batch does not have to have size p.)
The algorithm is called RBM-1 (shown in <ref>). Each iteration contains two main steps: (1) Randomly shuffling and dividing the particles into n batches; (2) evolving with interactions only turned on inside batches.



Clearly, the update equation (<ref>) can be rewritten as

    dX^i=-∇ V(X^i) dt+1/N-1∑_j: j≠ iK(X^j-X^i)dt+σ dB^i
    +χ_m,i(X(t)) dt,

where

    χ_m,i(X(t))=1/p-1∑_j∈𝒞_q,j≠ iK(X^i-X^j)-1/N-1∑_j: j≠ iK(X^i-X^j).


For a given x=(x^1,…,x^N)∈ℝ^Nd that is independent of the random division, we have (see Lemma <ref> for the proof)

    𝔼χ_m,i(x)=0.

This is a key observation which eventually leads to the convergence of the algorithms in expectation. 
As a remark, the position X=(X^1,…, X^N) itself depends on the random division. Hence, in general,

    𝔼χ_m,i(X(t))≠ 0.



Regarding the complexity, note that random division into n batches of equal size can be implemented using random permutation. The latter can be realized in O(N) operations by Durstenfeld's modern revision of Fisher-Yates shuffle algorithm <cit.>. (In MATLAB, one can use 'randperm(N)' to generate a random permutation. Then, the first p elements are considered to be in the first batch, the second p elements are in the second batch, etc). 

If one is to simulate up to time T, the total number of time steps is N_T as in (<ref>). 
The cost for each batch per iteration is ≲ Cp^2, where C is the unit time independent of N,p,τ. Since there are n=N/p batches, so the total cost per iteration is like ≲ CpN. The total complexity is therefore

    O(pN_TN)=O(p/τNT).
 
The ration p/τ affects the total complexity. If p=2 which we will use for the numerical simulations due to the simplicity in implementations, the complexity is O(N_TN)=O(TN/τ). The cost is reduced significantly from the naive simulation.

In section <ref>, one will see that the RBM converges due to some averaging effect in time, and 
the rate follows in a similar way as Law of Large Number results. The error bound is like C√(τ/p). This means the ratio τ/p controls the error theoretically if the ODEs/SDEs can be solved accurately. Hence, if one wants the error to be smaller than some given accuracy ϵ, the computational cost is like O(NT/ϵ^2).
The dependence in N is linear! The power for ϵ is 2, which is typical for Monte Carlo methods. To improve this, one may consider some variance reduction techniques (see for example <cit.>).

RBM-1 is in spirit similar to the stochastic gradient descent in machine learning (<cit.>).
Recently, there are some analysis of SGD in the mathematical viewpoints and applications to physical problems <cit.>.
RBM-1 can be used both for simulating the evolution of the measure (<ref>) or (<ref>) and for sampling from the equilibrium state π.





 §.§ The Random Batch Method with replacement

RBM-1 requires the random division, and the elements in different batches are different.
This is in fact the sampling without replacement. If one allows replacement, we have RBM-r'.





RBM-r' can be reformulated as Random Batch Method with replacement (RBM-r) that has some flavor of  the stochastic coordinate descent method <cit.>. Here, the pseudo-time s is introduced for convenience and s_m=mτ. Roughly, t_m corresponds to s_mN/p.






For such type of methods, we expect to have 

    Ñ_T=1/pN N_T

iterations to get comparable behaviors. However, each step is very cheap: one only needs O(p^2) work. Hence, we still expect the complexity to be O(pNN_T).

Though the same as RBM-r', 
RBM-r does not have explicit concept of time, since the positions of the particles are not changing simultaneously. Intuitively,  RBM-r can only give sampling for the invariant measure for the nonlinear Fokker-Planck equation (<ref>) while unable to simulate the dynamics of the distribution.  However, as RBM-r' indicates,  N/p iterations correspond to time τ and thus a single step in RBM-1, so it might still approximate the evolution of distributions. In fact, the example later confirms this and m(N/p) iterations indeed give acceptable approximation for the distribution at time mτ.
As a last comment, we may sometimes choose random time step; for example τ∼ Exp(Δ t) (the exponential distribution with parameter Δ t) such that 𝔼τ=Δ t. Intuitively, this may help to increase the noise level and avoid being trapped in some local minimizers of V.




 §.§ Some discussions about the algorithms



In this subsection we make some discussion about our algorithms and complementary remarks. 


  * For system (<ref>), the RBMs can be similarly developed. The batch size should be larger than or equal to J (p≥ J). After the random batch is chosen, the interactions only happen inside the randomly chosen batch. One should also adjust the magnitude of the interaction such that an analogy of Lemma <ref> holds. In other words, the expectation of the random forces should equal to those in (<ref>). The complexity clearly is reduced from O(N^J) to O(N).


  * The sizes of batch do not have to be equal. One can even choose the sizes randomly.
One can also adjust the batch size to adjust the noise level. 


  * After the batches are obtained, one can do parallel computing to update the positions of X for RBM-1. However, one has to do regrouping after one time interval which stops the parallelism. There is no big efficiency improvement unless solving the SDEs is expensive. Regarding complexity,  there is no big difference in complexity. For RBM-r, to get p indices, one needs to generate p random numbers from [0, 1], and the total complexity is O(N) for N/p iterations. For RBM-1,  as in Durstenfld's implementation of random permutation, one needs to generate N random numbers from [0, 1] while do some swapping operations for each element. (This may be different from the intuition established from drawing balls from a bag. For randomly dividing N  balls, one can choose p from N balls first, then p from N-p, …. It seems that the complexity becomes smaller and smaller as the total number of balls is becoming smaller and smaller. This is indeed not the case because in computer one can use arrays to store the 'balls' and one only needs to generate the indices without touching these 'balls'.) 

In summary, there is no big difference in complexity for the two algorithms. The advantage of RBM-1 might be its ability for parallelism during one interval while the advantage of RBM-r is its simplicity so that it is likely more flexible for extensions. 


  * In Lemma <ref>, for x independent of the random division, we prove that

    Var(χ_m,i(x))=(1/p-1-1/N-1)Λ_i(x),

where Λ_i(x) is independent of p.
This means for larger p, the variance is smaller and the noise level is lower. This noise somehow reflects the fluctuation of the empirical measure μ̃_N^(1) around μ_N^(1).


  * If K is a singular forcing (like Coulomb), we can do splitting method and have

    dX^i= 1/p-1∑_j∈𝒞,j≠ i K(X^i-X^j) dt,
       dX^i=-∇ V(X^i) dt+σ dB^i.

If p=2, the singular forcing term can be often solved analytically. This is another advantage of the stochastic algorithm: for the N-particle system, if the forcing K is singular so that the problem is stiff, explicit scheme needs very small time step while implicit scheme is hard to invert. Using the stochastic algorithm plus time splitting, the evolution can be solved exactly, thus avoiding stability constraint.




  * Suppose one aims at the law of particles in (<ref>)
with very large N. By mean field limit, one can choose to
solve (<ref>) directly. A possible way is to use particle method for (<ref>)  where the number of particles is much smaller than N and the masses for the particles in the numerical method can even be different (the particle blob method <cit.>).  For the interacting particle system of this numerical purpose, our randomized algorithms also apply well. In this regard, we provided an efficient numerical particle method for the mean-field (or kinetic) equations. In fact this is exactly the starting point of the binary algorithms in <cit.>.



  * To better simulate the invariant measure in the case of σ=0, one may add the Brownian motion ϵ_N dB^i where ϵ_N decreases with N. For example, we may set 

    ϵ_N=1/N^γ, γ>0.



  * If the initial distribution is far from the equilibrium and we aim to get the global equilibrium, one may need many iterations for convergence to the equilibrium. A possible way is to use the Gibbs distribution corresponding to V 

    ν(dx)∝exp(-V(x)) dx

 for initialization to reduce the number of iterations. This distribution can be sampled using the Markov chain Monte Carlo (MCMC) methods <cit.>. Of course, for special V such as quadratic functions, the initial distribution can be sampled directly without MCMC.
 
 
  * Compared with traditional treecode and fast multipole method <cit.> where the interaction kernels are required to decay, the RBMs work also for kernels that do not necessarily decay to reduce the computational cost to O(N). They can also be applied for kernels of the general form K(x,y) which may not be symmetric and may not have translational invariance.





§ AN ERROR ANALYSIS FOR RBM-1


In this section, we study RBM-1 proposed in Section <ref>. In particular, we will check how close it is to the fully coupled system (<ref>). We leave the study of RBM-r to the future. However, as commented above, we expect RBM-r to also work when RBM-1 works. Recall that each iteraction consists of random division and evolution. The mechanism that makes RBMs work is Lemma <ref> and small step size. The philosophy is as following. Suppose that the system has certain chaotic property. When τ is small enough, the accumulative behavior along many time steps will be roughly comparable to the average behavior, which is (<ref>) thanks to Lemma <ref>. This is similar to Law of Large Numbers, but on time.   


We assume the following conditions on the confining and interacting potentials:

Suppose V is strongly convex on ℝ^d so that x↦ V(x)-r/2|x|^2 is convex, and ∇ V, ∇^2V have polynomial growth (i.e. |∇ V(x)|+|∇^2 V(x)|≤ C(1+|x|^q) for some q>0).
Assume K(·) is bounded, Lipschitz on ℝ^d with Lipschitz constant L and has bounded second order derivatives. Moreover,

    r>2L.



The condition r>2L is to ensure that the evolution group for the deterministic part of (<ref>) is a contraction. In particular, if 

    Ẋ^i=-∇ V(X^i)+1/N-1∑_j: j≠ iK(X^i-X^j),  Ẏ^i=-∇ V(Y^i)+1/N-1∑_j: j≠ iK(Y^i-Y^j),

 then

    d/dt∑_i=1^N|X^i-Y^i|^2≤ -(r-2L) ∑_i=1^N|X^i-Y^i|^2.

For such systems, one needs σ>0 for a nontrivial equilibrium; otherwise, all the particles will go to a single point.

Our goal is to prove that the distribution generated by RBM-1 is close to the 
distribution of the marginal distribution of the N particle system in the Wasserstein distance. We recall that the Wasserstein-2 distance is given by <cit.>

    W_2(μ, ν)=(inf_γ∈Π(μ,ν)∫_ℝ^d×ℝ^d|x-y|^2 dγ)^1/2,

where Π(μ,ν) means all the joint distributions whose marginal distributions are μ and ν respectively.

In order to achieve this goal, we consider the synchronous coupling between (<ref>) and (<ref>). 
We denote X^i the solution obtained by the N interacting particle system while X̃^i the solution obtained by RBM-1. Correspondingly, B̃^i will denote the Brownian motion used in RBM-1. Note that both (<ref>) and (<ref>) have exchangeability. This means, for example, the joint distribution of X^i's is symmetric. Consequently, X̃^i and X̃^j (for all i and j) have the same distribution.  We construct the coupling as follows. For i=1,2,…, N,

    X̃^i(0)=X^i(0)= X_0^i ∼ν,
        B^i(t)=B̃^i(t).

With this coupling, we define the error process

    Z^i(t):=X̃^i(t)-X^i(t).



Let ξ_m-1 denote the random division of batches at t_m-1, and  define the filtration {ℱ_m}_m≥ 0 by

    ℱ_m-1=σ(X_0^i, B^i(t), ξ_j; t≤ t_m-1, j≤ m-1).

In other words, ℱ_m-1 is the σ-algebra generated by the initial values X_0^i (i=1,…, N), B^i(t), t≤ t_m-1, and ξ_j, j≤ m-1. Hence, ℱ_m-1 contains the information of how batches are constructed for t∈ [t_m-1, t_m). We also introduce the filtration {𝒢_m}_m≥ 0 by

    𝒢_m-1=σ(X_0^i, B^i(t), ξ_j; t≤ t_m-1, j≤ m-2).

If we use σ(ξ_m-1) to mean the σ-algebra generated by ξ_m-1, the random division of batches at t_m-1, then ℱ_m-1=σ(𝒢_m-1∪σ(ξ_m-1)). Throughout this section, C will denote generic constants whose concrete value can change from line to line. We use · to represent the L^2(Ω) norm (Ω means the sample space):

    v=√(𝔼|v|^2).


We now state the theorem.

Suppose Assumption <ref> holds. With the coupling constructed above,

    sup_t≥ 0Z^1(t)≤ C√(τ/p-1+τ^2),

where C is independent of N, p and t. Consequently, let μ_N^(1)(t) be the first marginal for (<ref>) and μ̃^(1)_N be the first marginal for system (<ref>). Then

    sup_t≥ 0W_2(μ^(1)_N(t), μ̃^(1)_N(t))≤ C√(τ/p-1+τ^2)≤  C√(τ).




We need some preparation for the proof.  Lemma <ref> is a type of consistency lemma for RBMs while Lemma <ref> somehow shows the stability of RBM-1.

Consider p≥ 2 and a given fixed x=(x^1, …, x^N)∈ℝ^Nd. Then, for all i,

    𝔼χ_m,i(x)=0,

where the expectation is taken with respect to the random division of batches. Moreover, the variance is given by

    (χ_m,i(x))=𝔼|χ_m,i(x)|^2=(1/p-1-1/N-1)Λ_i(x),

where

    Λ_i(x) :=1/N-2∑_j: j≠ i| K(x^i-x^j)-1/N-1∑_k:k≠ iK(x^i-x^k)|^2.



We use I(i,j)=1 to indicate that i,j are in the same batch. 
We rewrite

    f_m,i:=1/p-1∑_j:∈𝒞_q,j≠ iK(x^i-x^j)
    =1/p-1∑_j: j≠ i K(x^i-x^j)1_I(i,j)=1.

We first show that

    𝔼1_I(i,j)=1=p-1/N-1.



As is well-known, there are 

    M(n):=(pn)!/(p!)^n n!

ways of dividing pn distinguishable objects into n batches of size p. 
For a given (i, j) pair, to compute 𝔼1_I(i,j)=1 (the probability that i, j are in the same batch), we compute the number of ways to group (i, j) together. We first choose p-2 objects from np-2 to group with i,j and then form a batch. Then, divide the remaining into n-1 batches.   Hence,
the number of ways to make i,j in the same batch is np-2  p-2M(n-1) and thus

    𝔼1_I(i,j)=1
    =np-2  p-2M(n-1)/M(n)=p-1/N-1.

This then proves that

    𝔼f_m,i=1/N-1∑_j:j≠ iK(x^i-x^j),

and thus (<ref>) follows.

Now, let us consider the variance. We first compute the second moment of ξ_m,i:

    𝔼|f_m,i|^2
    =   1/(p-1)^2∑_j:j≠ i|K(x^i-x^j)|^2 ℙ(I(i,j)=1)
       +1/(p-1)^2∑_j,k: j≠ k, j≠ i,k≠ iK(x^i-x^j)· K(x^i-x^k) ℙ(I(i,j)I(i,k)=1).

By similar argument,

    ℙ( I(i,j)I(i,k)=1)=np-3  p-3M(n-1)/M(n)
    =(p-1)(p-2)/(N-1)(N-2).

Hence,

    𝔼|χ_m,i(x)|^2
    =    𝔼|f_m,i|^2-(𝔼|f_m,i|)^2 
    
    =   (1/p-1-1/N-1)(1/N-1∑_j≠ i|K(x^i-x^j)|^2 
       -1/(N-1)(N-2)∑_i,j: j≠ k, j≠ i, k≠ iK(x^i-x^j)· K(x^i-x^k)) 
    
    =   (1/p-1-1/N-1)1/N-2∑_j: j≠ i|K(x^i-x^j)-1/N-1∑_k: k≠ iK(x^i-x^k)|^2.



The following simple fact of probability will be useful later.

Fix i∈{1,…, N}. Let 𝒞_θ be the random batch of size p that contains i in the random division. 
Let Y_i (1≤ i≤ N) be N random variables (or random vectors) that are independent of 𝒞_θ. Then,

    1/p-1∑_j∈𝒞_θ,j≠ iY_j≤max_j Y_j.



We again use I(i,j)=1 to indicate that i,j are in the same batch.  By the definition and independence,

    1/p-1∑_j∈𝒞_θ,j≠ iY_j^2
       =1/(p-1)^2𝔼|∑_j: j≠ iI(i,j)Y_j|^2 
       =1/(p-1)^2∑_j,k:j≠ i,k≠ i𝔼(I(i,j)I(i,k))𝔼(Y_j· Y_k)
       ≤ (max_j Y_j)^2 1/(p-1)^2∑_j,k:j≠ i,k≠ i𝔼(I(i,j)I(i,k))
       =(max_j Y_j)^2 1/p-1∑_j∈𝒞_θ,j≠ i 1 ^2.

The claim thus follows. Note that the independence is used in the second equality.


We now move onto the estimates of X and X̃.

Suppose the coupling constructed in (<ref>) and Assumption <ref> hold. Then, for any q≥ 2, there exists a constant C_q independent of N such that for any i

    sup_t≥ 0(𝔼|X^i(t)|^q+𝔼|X̃^i(t)|^q)≤ C_q.

Besides, for any m>0 and q≥ 2,

    sup_t∈ [t_m-1, t_m)|𝔼(|X̃^i(t)|^q |ℱ_m-1)|
    ≤ C_1|X̃^i(t_m-1)|^q+C_2.

holds almost surely.

Moreover, for t∈ [t_m-1, t_m) and all i, there exist a constant C>0 independent of ξ_m, m, N and an index q>0 such almost surely that 

    |𝔼(X̃^i(t)-X̃^i(t_m-1) | ℱ_m-1)|≤ C (1+|X̃^i(t_m-1)|^q)τ,

and in L^2(Ω) that

    𝔼(|X̃^i(t)-X̃^i(t_m-1)|^2 | ℱ_m-1) ≤ C τ.




Consider system (<ref>) first. By Itô's calculus,

    d/dt𝔼|X^i|^q=    q𝔼|X^i|^q-2(-X^i·∇ V(X^i)+1/N-1∑_j: j≠ iX^i· K(X^i-X^j)) 
       +1/2q(q+d-2)σ^2𝔼|X^i|^q-2.

Note that

    X^i·∇ V(X^i)=(X^i-0)· (∇ V(X^i)-∇ V(0))+X^i·∇ V(0) ≥ r|X^i|^2+X^i·∇ V(0).

Recalling also that K is bounded,

    d/dt𝔼|X^i|^q≤ -qr 𝔼|X^i|^q
    +q(K_∞+|∇ V(0)|)𝔼|X^i|^q-1
    +1/2q(q+d-2)σ^2𝔼|X^i|^q-2.

By Young's inequality,

    𝔼(|X^i|^q-1)
    ≤(q-1)ν/q𝔼(|X^i|^q)
    +1/qν^q-1,

the second term on the right hand side of (<ref>) is therefore controlled. 
If q=2, the last term on the right hand side of (<ref>) is controlled by a constant. Otherwise, one can apply Young's inequality similarly to control it with 𝔼(|X^i|^q). Clearly, when choosing ν fixed but small enough, 𝔼|X_i|^q can be uniformly bounded in time for any q≥ 2, and the bound is independent of N.

For X̃^i, we first consider a given random division so that the equation is given by

    dX̃^i=-∇ V(X̃^i)dt+1/p-1∑_j∈𝒞_θ,j≠ iK(X̃^i-X̃^j) dt+σ dB^i,

where 𝒞_θ is the random batch that contains i from the random division at t_m-1, or ξ_m-1.
Now, consider that t∈ [t_m-1, t_m). Conditioning on ℱ_m-1 and applying Itô's calculus on [t_m-1, t_m), one also has

    d/dt𝔼(|X̃^i|^q|ℱ_m-1)=q𝔼[|X̃^i|^q-2(-X̃^i·∇ V(X̃^i)+
    X̃^i·1/p-1∑_j∈𝒞_θ,j≠ iK(X̃^i-X̃^j)))| ℱ_m-1] 
    +1/2q(q+d-2)σ^2𝔼(|X̃^i|^q-2|ℱ_m-1).

Using similar estimates,

    d/dt𝔼(|X̃^i|^q|ℱ_m-1)≤ -r_1  𝔼(|X̃^i|^q|ℱ_m-1)+C_1,

for some r_1>0 and  constant C_1 that are deterministic. This clearly yields (<ref>).
Taking expectation about the randomness in ℱ_m-1 on both sides, one then obtains the same inequality as for X. Hence, (<ref>) follows.


Equation (<ref>) gives

    𝔼(X̃^i(t)-X̃^i(t_m-1) | ℱ_m-1)
    =    -∫_t_m-1^t𝔼(∇ V(X̃^i)| ℱ_m-1) ds 
       +∫_t_m-1^t𝔼( 1/p-1∑_j∈𝒞_θ,j≠ iK(X̃^i-X̃^j) | ℱ_m-1) ds.

Note that K is bounded and |∇ V(x)|≤ C(1+|x|^q) for some q>0.   
Together with the first part of this lemma (moment control), this implies (<ref>).

For the last claim, Itô's formula implies that

    d/dt𝔼[(X̃^i(t)-X̃^i(t_m-1))^2|ℱ_m-1] 
    
    =2𝔼[ (X̃^i(t)-X̃^i(t_m-1))·(-∇ V(X̃^i)+1/p-1∑_j∈𝒞_θ,j≠ iK(X̃^i-X̃^j))
    |ℱ_m-1]
    +σ^2.

Similarly,

    𝔼[ (X̃^i(t)-X̃^i(t_m-1))·(-∇ V(X̃^i)+1/p-1∑_j∈𝒞_θ,j≠ iK(X̃^i-X̃^j))
    |ℱ_m-1] ≤ C

is uniformly bounded by Hölder's inequality for conditional expectations (though X̃^i(t)-X̃^i(t_m-1) is small, there is no need to get finer bound since we have extra σ^2). The last claim also follows.



The following gives the control on Z.

For t∈ [t_m-1, t_m), 

    Z^i(t)-Z^i(t_m-1)≤ C τ,
       |𝔼((Z^i(t)-Z^i(t_m-1))·χ_m,i(X(t)))|≤Λ_i_∞τ/p-1 
    +C(τ^2+Z^i(t)τ).



By the coupling in (<ref>), Z^i satisfies on t∈ [t_m-1, t_m)

    dZ^i=-(∇ V(X̃^i)-∇ V(X^i))dt
    +1/p-1∑_j∈𝒞_θ,j≠ iK(X̃^i-X̃^j) dt-1/N-1∑_j: j≠ iK(X^i-X^j) dt.

Since ∇ V has polynomial growth, the claim for Z^i(t)-Z^i(t_m-1) is then an easy consequence of the q-moment estimates in Lemma <ref>.

Moreover, 

    |∇ V(X̃^i)-∇ V(X^i)|≤∫_0^1|(X̃^i-X^i)·∇^2V((1-z)X̃^i+zX^i)| dz,

and ∇^2V((1-z)X̃^i+zX^i) is controlled by C(|X̃^i|^q_1+|X^i|^q_1) for some q_1>0. Hence,

    𝔼|∇ V(X̃^i)-∇ V(X^i)||χ_m,i(X(t'))|
       ≤χ_m,i(X(t'))_∞X̃^i-X^i (𝔼(|X̃^i|^q_1+|X^i|^q_1)^2)^1/2
       ≤ CZ^i(t),

by Lemma <ref>, where t'∈ [t_m-1, t_m] is arbitraray.

Integrating (<ref>) in time over [t_m-1, t], then dotting with χ_m,i(X(t)), and taking the expectation,  one gets  

    | 𝔼((Z^i(t)-Z^i(t_m-1))·χ_m,i(X(t))) |
    ≤ C∫_t_m-1^tZ^i(s) ds+
    ∫_t_m-1^t𝔼[(1/p-1∑_j∈𝒞_θ,j≠ iδ K_ij(s))·χ_m,i(X(t))] ds
    +𝔼∫_t_m-1^tχ_m,i(X(s))·χ_m,i(X(t)) ds,

where

    δ K_ij(s):=K(X̃^i(s)-X̃^j(s))-K(X^i(s)-X^j(s)).


Since Z^i(s)≤Z^1(t)+Cτ, the first term is controlled by CZ^i(t)τ+Cτ^2. Since K is Lipschitz continuous,

    |δ K_ij(s)|≤ L(|Z^i(s)|+|Z^j(s)|).

Hence,

    (1/p-1∑_j∈𝒞_θ,j≠ iδ K_ij(s))≤ L(Z^i(s)+1/p-1∑_j∈𝒞_θ,j≠ i|Z^j(s)|).

Now, Z^j(s) depends on 𝒞_θ so it cannot be estimated directly. 
Dotting (<ref>) with Z^i, one has 1/2d/dt|Z^i|^2≤ C|Z^i| almost surely,
which gives 

    |Z^i(s)|≤ |Z^i(t_m-1)|+Cτ

 almost surely. Since Z^i(t_m-1) is independent of 𝒞_θ, Lemma <ref> then gives us that
 
    1/p-1∑_j∈𝒞_θ,j≠ i|Z^j(s)|≤Z^1(t_m-1)+Cτ≤Z^1(t)+Z^1(t)-Z^1(t_m-1)+Cτ.

 Since Z^1(t)-Z^1(t_m-1)≤ Cτ by what has been just proved, this term is bounded by Z^1(t)+Cτ.
 Hence,
 
    (1/p-1∑_j∈𝒞_θ,j≠ iδ K_ij(s))≤ 2L(Z^1(t)+Cτ).



Since X is independent of 𝒞_θ, applying Lemma <ref>,

    𝔼(χ_m,i(X(s))·χ_m,i(X(t)))
    ≤1/p-1Λ_i_∞,
 
and the claim follows.


Now, we are ready to prove Theorem <ref>.

With the coupling (<ref>), the continuous process Z^i satisfies for t∈ [t_m-1, t_m) that

    dZ^i =    -(∇ V(X̃^i)-∇ V(X^i))dt 
       +1/N-1∑_j: j≠ i(K(X̃^i-X̃^j)-K(X^i-X^j)) dt +χ_m,i(X̃(t))dt.

Using the strong convexity of V and Lipscthitz continuity of K, 

    1/2d𝔼∑_i=1^N(Z^i)^2
    ≤ -(r-2L)𝔼∑_i=1^N (Z^i)^2 dt
    +∑_i=1^N 𝔼Z^i(t)·χ_m,i(X̃(t)) dt,

where X̃=(X̃^1, …, X̃^N)∈ℝ^Nd.  Due to the exchangeability, 

    R(t):=1/N∑_i=1^N 𝔼Z^i(t)·χ_m,i(X̃(t))=𝔼Z^1(t)·χ_m,1(X̃(t)).

For notational convenience, we introduce

    u(t):=1/N𝔼∑_i=1^N(Z^i)^2=Z^1(t)^2.

Then,

    d/dtu(t)≤ -2(r-2L)u(t)+2R(t).


We now estimate R(t). We rewrite it as

    R(t)=    𝔼Z^1(t_m-1)·χ_m,1(X̃(t_m-1))
    +𝔼Z^1(t_m-1)· (χ_m,1(X̃(t))-χ_m,1(X̃(t_m-1)))
          + 𝔼(Z^1(t)-Z^1(t_m-1))·χ_m,1(X(t)) 
        +𝔼(Z^1(t)-Z^1(t_m-1))· (χ_m,1(X̃(t))-χ_m,1(X(t)))
    
              =:   I_1+I_2+I_3+I_4.


By Lemma <ref>, we have 

    I_1     = 𝔼(𝔼(Z^1(t_m-1)·χ_m,1(X̃(t_m-1)) )|𝒢_m-1) 
       =𝔼(Z^1(t_m-1)·𝔼 (χ_m,1(X̃(t_m-1))
    |𝒢_m-1) )=0.

The second equality holds because Z^1(t_m-1) is adapted to 𝒢_m-1 (note that Z^1 is continuous in time, so Z^1(t_m-1^-)=Z^1(t_m-1)). Note that 𝒢_m-1 is independent of the random division at t_m-1 and this is why we can get I_1=0 by regarding X̃(t_m-1) as a given fixed point.  Equation (<ref>) is the consistency that ensures convergence. 


For term I_2, 

    I_2    =𝔼( Z^1(t_m-1)·𝔼(χ_m,1(X̃(t))-χ_m,1(X̃(t_m-1)) |ℱ_m-1) ) 
       ≤ CZ^1(t_m-1)𝔼(χ_m,i(X̃(s))-χ_m,i(X̃(t_m-1)) 
    |ℱ_m-1).

Clearly,

    𝔼(χ_m,i(X̃(s))-χ_m,i(X̃(t_m-1)) |ℱ_m-1)
    =1/p-1∑_j∈𝒞_θ,j≠ i𝔼(δK̃^ij|ℱ_m-1)
    -1/N-1∑_j:j≠ i𝔼(δK̃^ij|ℱ_m-1),

where 

    δK̃^ij=K(X̃^i(s)-X̃^j(s))-K(X̃^i(t_m-1)-X̃^j(t_m-1)).

Denote δX̃^j:=X̃^j(s)-X̃(t_m-1). Performing Taylor expansion around t_m-1, one has

    δK̃^ij=∇ K(X̃^i(t_m-1)-X̃^j(t_m-1))·(δX̃^i-δX̃^j)
    +1/2M:(δX̃^i-δX̃^j)⊗(δX̃^i-δX̃^j),

with M being a random variable (matrix) bounded by ∇^2K_∞.
By (<ref>), one finds that
 
    |𝔼(δK̃^ij|ℱ_m-1)|
     ≤ CL(1+|X̃^i(t_m-1)|^q+|X̃^j(t_m-1)|^q)τ+C𝔼(|δX̃^i|^2+|δX̃^j|^2  |ℱ_m-1).

By (<ref>) and the fact that X̃^i(t_m-1) is independent of 𝒞_θ, Lemma <ref> then can control 

    I_2≤ CZ^1(t_m-1)τ≤ CZ^1(t)τ+Cτ^2.


By Lemma <ref>, the third term on the right hand side of (<ref>) is bounded by

    I_3  ≤ Cτ/p-1+CZ^1(t)τ+Cτ^2.


The I_4 term can be controlled by

    I_4≤Z^1(t)-Z^1(t_m-1)( 1/p-1∑_j∈𝒞_θδ K_ij(t)
    +1/N-1∑_j:j≠ iδ K_ij(t) ),

where δ K_ij(t) is defined in (<ref>).
The first term inside the parenthesis on the right hand side is estimated by (<ref>). The second term is 
straightforward. Hence, noticing Lemma <ref>, one has 

    I_4≤Z^1(t)-Z^1(t_m-1)(4LZ^1(t)+Cτ)≤ CZ^1(t)τ+Cτ^2.



Hence,

    R(t) ≤ CZ^1(t)τ+Cτ/p-1+Cτ^2

and thus

    d/dtu(t)
    ≤ -2(r-2L)u(t)+C√(u(t))τ+Cτ/p-1+Cτ^2.

This inequality therefore implies that

    sup_t≥ 0𝔼|X̃^1-X^1|^2=sup_t≥ 0u(t) ≤ Cτ/p-1+Cτ^2.

The closeness between μ_N^(1) and μ̃_N^(1) in W_2 distance is a simple application of definition (<ref>).



We point out that the constant C in Theorem <ref> stays bounded as σ→ 0. The proof in Theorem <ref> is valid for σ=0, when there is no Brownian motion.  In I_1 and I_2, we used
χ_m,i(X̃(t_m-1)) to split the terms. If we use χ_m,i(X(t)) to split the terms and use similar techniques
as in I_3, I_4,  then there is no requirement on ∇^2 K but we need r>6L.

The first part in Theorem <ref>  in fact says that we can approximate the trajectories of the particles, which is the convergence in strong sense.





It is straightforward to conclude that RBM-1 can approximate the mean field measure (solution to (<ref>)) since the mean field result is well known (see, for example, <cit.>). In fact, as shown in Fig. <ref>, the triangle inequality yields the following.

Suppose Assumption <ref> holds, then

    W_2(μ̃_N^(1)(t), μ(t))≤ C(√(τ)+N^-1/2+ϵ)

for any ϵ>0.

The error bound between μ̃_N^(1) and μ_N^(1) is C√(τ) uniformly in time with C independent of N. The error bound between μ_N^(1) and μ is given by CN^-1/2+ϵ (<cit.>). 
As shown in the figure, if one takes N→∞ first, one may get some measure μ̃(t), which is the mean field limit of the RBM-1 method. Analyzing this measure will be an interesting problem for the future.

As shown in Fig. <ref>, the RBMs are a class of Asymptotic-Preserving schemes <cit.> for particle system (<ref>) under the scalings of mean field limit, in that they approximate the particle system (<ref>) with an error of O(τ) independent of N, and as N →∞, they become good approximations to the limiting mean-field equation (<ref>).




If the deterministic flow does not have the contraction property, we may only be able to prove the convergence of our methods on finite time interval [0, T]. Moreover, the mean field limit by Dobrushin's estimate holds on finite time interval <cit.>. Consequently, it is possible to show that 

    W_2(μ̃^(1)_N(t), μ(t))≤ C(T)(√(τ)+N^-1/2), ∀ t∈ [0, T],

for general interaction potentials, where the coefficient depends on T now.





For N≫ 1 when the theory of mean field limit takes effect, the chaos could be created or propagated (in terms of the distribution of exchangeable particles). In other words, the j-marginal will be the tensor product of j copies of one particle distribution. In this case, the joint distribution could be of low rank, which helps RBMs: the number of time steps N_T can be much smaller than N/p to get reasonably good approximation ability. This is definitely something that can be explored in the future.





The fast Monte Carlo method in <cit.> uses subsampling in columns or rows for matrix-matrix multiplication. The idea of subsampling as unbiased estimator is the same as the one for RBMs. 
However, besides subsampling, another key idea of RBM is re-subsampling at later time intervals so that the dynamics will be averaged out in time, which differs from the random algorithms in matrix computations. This Law of Large Numbers type feature used in time guarantees the convergence of RBMs (as indicated by the Law of Large Number type error bound in Theorem <ref>). 





The analysis here should also work for (<ref>), provided some analogies of Lemma <ref> hold (as discussed in section <ref>). For techanical needs, one will assume r large enough so that the semigroup given by the deterministic flow is again a contraction. Besides, the analysis can be carried over directly to the O(1) time simulation of second order systems such as Hamiltonian dynamics, and we provide such an example in Appendix <ref>.




If one considers particles associated with charges or masses, they are not exchangeable and the current error analysis does not work directly. The consistency lemma (Lemma <ref>) in fact has nothing to do with exchangeability and can be generalized to interactions that dependent on the specific chosen particles.  The current proof of convergence seems to rely on exchangeability, but one can consider the charge or mass as an extended coordinate <cit.>. With this new insight, the new 'particles' become exchangeable. We believe RBMs also work for such systems but we feel it better to leave this as a subsequent project. 







§ NUMERICAL VERIFICATION


In this section, we run some numerical tests to evaluate the RBMs and verify our theory in section <ref>. 
The first example is a simple artificial example to test the dependence of the errors on N and τ, which satisfies the conditions in Theorem <ref>. The second example is the Dyson Brownian motion, and it does not satisfy the conditions in Theorem <ref> as the kernel is singular. RBMs work well for both examples, which implies that the algorithms can be effective for systems far beyond the systems mentioned in section <ref>.



 §.§ A simple test example




We now use the following simple test example to check how the error in RBM-1 depends on N, τ and T (T is the time point where we compute the numerical solutions). Here, the spatial dimension is 1 (d=1)

    Ẋ^i=-β X^i+1/N-1∑_j: j≠ iX^i-X^j/1+|X^i-X^j|^2.

The interaction is clearly smooth, bounded and with bounded derivatives. Moreover, it has a long-range interaction. 

In principle, to evaluate E(T)=√(𝔼|X̃^1(T)-X^1(T)|^2), we need to run many independent experiments and use empirical mean for the approximation. Doing this is clearly very expensive. Alternatively, we only run one experiment and use

    Ê(T):=√(1/N∑_i=1^N|X̃^i(T)-X^i(T)|^2)

to approximate E(T).


In Fig. <ref>, we show the numerical results for T=1. The initial distribution is taken from

    ρ_0(x)=√(4-x^2)/2π,

by the Metropolis-Hastings MCMC algorithm <cit.>.
The reference solution X^i(T) is obtained by solving the fully coupled system using the forward Euler scheme with τ=2^-15. The solution X̃^i(T) is generated by RBM-1 with p=2. Each step is solved by the forward Euler method with τ from 2^-7 to 2^-4. We considered N=50, 500, 2000 respectively.  We plot the error Ê(T) versus τ for these three N values. The first picture in Fig. <ref> is for β=0 while the second picture in Fig. <ref> is for β=1. Clearly, the error is insensitive to the change of N. When N is small, like N=50, the fluctuation in the error is kind of clear. When N is large, in the log-log scale, the curve is already close to straight lines with slope approximately 0.5, meaning that the error indeed decays like √(τ). 




In Fig. <ref>, we take N=500. The reference solution is again computed by solving the fully coupled system using the forward Euler scheme with τ=2^-15. The algorithm is performed by taking τ=2^-7. If there is confining potential, the error stays bounded as T increases. However, if there is no confining potential, the error clearly grows, consistent with Remark <ref>. This is indeed natural even for usual ODE discretization for the fully coupled system (<ref>). In fact, if there is no confining potential, the numerical error grows with T for the forward Euler method.




To test the practical performance of RBM-1, we show the error versus true CPU time in Fig. <ref> for N=100, 10^3, 10^4 respectively. The experiment was done using MATLAB R2015b on macOS Mojave system with 1.6 GHz Intel Core i5 processor and 4GB DDR3 memory. The data were collected by computing the times and errors using forward Euler again for steps τ=2^-6-2^-1, and then plotting errors versus times.  The lines with circles are for RBM-1 while lines with diamonds are for direct computation using forward Euler discretization. One can see clearly that
for the error tolerance ϵ≳ 10^-2 when random algorithms are acceptable, RBM-1 saves a lot of time
compared with the direct ODE solvers. For example, when N=10^3, 
RBM-1 is faster to achieve the given error tolerance ϵ≥ 10^-2; for N=10^4, RBM is much faster in the regime shown in the figure.
However, since RBMs converge with only 1/2 order, the time needed for RBM can be larger if one desires very accurate result and N is not big, and even for this case, one can consider some variance reduction method to improve the accuracy of RBMs. For N≳ 10^6, direct method takes very long time (it scales like N^2), and uses much memory resource, so the direct method is already unacceptable. RBMs can be run in a reasonable time amount, and this is exactly one of the advantages of RBMs.

We now modify the above example to a second order system as studied in Appendix <ref>, which is a Hamiltonian system, to show that RBM-1 also works for second order systems on finite time. How to develop random algorithms for Hamiltonian systems for long time is an interesting question, and we think it better to leave it for future research. In particular, we consider the following system of equations for i=1,…, N:

    {   Ẋ^i=V^i,
       V̇^i=1/N-1∑_j: j≠ iX^i-X^j/1+|X^i-X^j|^2.
    .


For time integration, we use the Verlet scheme

    {   X_1^i=X_0^i+V_0^i τ+1/2F_0^iτ^2,
       X_n+1^i=2*X_n^i-X_n-1^i+F_n^iτ^2,  n≥ 1.
    .

For the reference solution, F^i is the full force using all the particles, while for the random algorithms, F^i is computed using RBM-1 with p=2 as before. We sample the initial positions again from ρ_0(x)=√(4-x^2)/2π while sample the initial velocities  from the normal distribution  𝒩(0, 1).
The error as in (<ref>) is again computed at T=1. We consider N=50, 500, 2000 respectively and the reference solution is obtained using τ=2^-15.  The results are shown in Fig. <ref>. The curves are close to straight lines in loglog scale with slope approximately 0.5, meaning that the errors indeed decay like √(τ) for second order systems as well. 








 §.§ The Dyson Brownian motion

Now, we consider a typical example in random matrix theory <cit.> to test the difference between RBM-1 and RBM-r. 
The random matrix we consider is a Hermitian matrix valued Ornstein-Uhlenbeck process 

    dA=-β A dt+1/√(N)dB,

where the matrix B is a Hermitian matrix consisting of some Brownian motions. In particular, the diagonal elements are independent standard Brownian motions. The off-diagonal elements in the upper triangular half are of the form 1/√(2)(B_R+i B_I) where B_R and B_I are independent standard Brownian motions. The lower triangular half elements are determined using the Hermitian property. By Itô's calculus <cit.>, it can be shown that the eigenvalues of A satisfy the following  system of SDEs (1≤ j≤ N), called the Dyson Brownian motion:

    dλ_j(t) =-βλ_j(t) dt+1/N∑_k: k≠ j1/λ_j-λ_kdt
    +1/√(N) dB_j,

where {B_j}'s are independent standard Brownian motions. The Brownian motion effect is small when N is large. This system therefore should have similar effects as system (<ref>) with σ=0. The limiting equation for N→∞ is given by <cit.>

    ∂_tρ(x,t)+∂_x(ρ(u-β x))=0,   u(x, t)=π(Hρ)(x, t),

where ρ is the density for λ as N→∞, H(·) is the Hilbert transform on ℝ, and π=3.14… is the circumference ratio.

Below we consider

    β=1.

It is shown that the corresponding limiting equation (<ref>) has an invariant measure, given by the semicircle law:

    ρ(x)=1/π√(2-x^2)






To numerically test the behavior of our methods, we note an analytic solution to the limiting equation (<ref>)

    ρ(x, t)=√(2σ(t)-x^2)/σ(t)π,  σ(t)=1+e^-2t.

The prefactor for the interaction in the equation is 1/N. For convenience, we change it to 1/(N-1) without introducing significant difference. For each iteration, since the force is singular, we adopted the splitting strategy mentioned in section <ref>. 
In particular, define

    X^ij:=X^i-X^j.

The SDE solving step in RBM-1 is given by


  * 
    Y_m^i=1/2(X_m-1^i+X_m-1^j)
                          +sgn(X_m-1^ij)√(|X_m-1^ij|^2+4τ),
    
            Y_m^j=1/2(X_m-1^i+X_m-1^j)
                          -sgn(X_m-1^ij)√(|X_m-1^ij|^2+4τ).


  * 
    X_m^i=Y_m^i-τ  Y_m^i+√(τ/N)z^i,  
           X_m^j=Y^j(t_m)-τ  Y_m^j+√(τ/N)z^j.


Here, z^i,z^j∼𝒩(0, 1).

In Fig. <ref>, we show the numerical results using RBM-1 in section <ref>. The initial data (setting t=0 in (<ref>)) are sampled using the Metropolis Hastings algorithm <cit.>. We plot the results at t=0.5 and t=5. The number of particles is N=10^5 while we use τ=10^-3 for time step. As can be seen, RBM-1 can successfully recover the evolution of distribution and the equilibrium semicircle law (<ref>), as desired.
In Fig. <ref>, the results of RBM-r are shown. Again, we take N=10^5 and τ=10^-3. 
Within one iteration,  the same splitting scheme above is used. 
We find that RBM-r indeed has comparable results with RBM-1. 
Though RBM-r seemingly cannot simulate the dynamics of the distributions, N/2 iterations in fact has comparable behavior for time τ. This interesting observation confirms that  RBM-r can capture the dynamics for some examples.

Since both stochastic algorithms give similar behaviors, in later examples, we only use one of them to implement for each example. If we care more about the dynamical behavior, we use RBM-1 (see the two examples in Section <ref>). Otherwise, we use RBM-r. (In fact, the two algorithms do not show significant difference, even for evolutional problems.)






§ APPLICATIONS


In this section, we apply RBMs to some examples from physics, social and data sciences. 
Most of these examples do not fit into the theoretic framework of the analysis in section <ref>. For example, in section <ref>, the kernel is singular;  in <ref>, the equilibrium still exists though the confining potential does not exist. Even though the analysis in section <ref> does not apply for some applications, RBMs turn out to work quite well. For the wealth problem in <ref>, we can recover the equilibrium extremely well, though the interaction kernel does not decay and there is no confining potential.
On one hand, the positive results give more supports to the algorithms; on the other hand, for some applications, the stochastic algorithms can be regarded as new models for the underlying problems. 



 §.§ Charged particles on the sphere


The traditional Thomson problem is to determine the stable configuration of N electrons on a sphere.
When N becomes large, this could lead to the so-called spherical crystals (<cit.>). The configuration may have some meta-states (local minimizers of the energy surface). When the number of particles is large, the spherical crystals have defects due to the topology of the sphere <cit.>. 

In the N→∞ limit, hopefully, we will have a continuous distribution of charges on the sphere ρ(·). The problem then corresponds to determining ρ such that the energy 

    E(ρ)=1/2∬_S× S1/|x-y|ρ(x)ρ(y) dS_xdS_y

is minimized. It is unclear how the energies corresponding to local minimizers are distributed (if there are any).

Regarding charges with surface densities, let us make a mathematical remark. Suppose S is a surface that divides the whole space ℝ^d into two halves. Assume there is a continuous distribution of charges on S with density ρ. Let φ^±(x) be the limits of the potential on the two sides of S, and φ(x):=φ^+(x)+φ^-(x). Then, one has

    (-Δ)_S^1/2φ+s(φ)=ρ,

where (-Δ)_S^1/2 is the 1/2 fractional Laplacian on S and s(φ) is some pseudo-differential operators with a symbol of degree lower than 1. In other words, to the leading order, the 1/2 fractional Laplacian of φ equals ρ. In the case that S is a plane or a circle in 2D plane, s(φ)=0. In general, s(φ)≠ 0. In fact, by the jump condition of electric fields, 

    ρ=E^+· n^++E^-· n^-=-∂φ^+/∂ n^+-∂φ^-/∂ n^-.

It is well-known that the Dirichlet to Neumann operator ℒ is related to the 1/2-fractional Laplacian
by

    ℒφ=-(-Δ)_S^1/2f+r(f),

where n is the normal vector pointing into the side where the harmonic extension is performed and r(f) includes lower terms. In the case that S is a plane r(f)=0 <cit.>. (In fact, if S is the base of cylinders, ℒ=-(-Δ)_S^1/2 as well. See <cit.>. )
 In the case of the unit circle in 2D plane, one can refer to <cit.>. For spheres in higher dimensions, s≠ 0. With (<ref>) and (<ref>), (<ref>) follows. 
 

  For the unit circle in 2D plane, one has a simple way to see (<ref>) with r(φ)=0. Let f(z) be analytic inside Ω such that φ=Re(f(z)).  On the boundary, z=e^iθ
f^-=∑_k≥ 0f̂_k e^ikθ,
so that f(z)=∑_k≥ 0f̂_k z^n for z inside the disk. Clearly, 

    ∂φ^-/∂ n^-=Re(∂ f/∂ n^-)
    =-Re(zf'(z))|_z∈∂Ω=-Re(∑_k≥ 0|k|f̂_k e^ikθ).

The claim is true for the interior domain. For the exterior domain, one has ∑_k≤ 0 and
∂φ^+/∂ n^+=Re(∂ f/∂ n^+)
=Re(zf'(z))|_z∈∂Ω. The claim again holds.



     Interacting particle systems on the sphere can be realized experimentally by beads on water droplets immersed in oil  <cit.>. By adjusting the environmental solution, the interacting kernel K(·) can also be changed, which does not have to be the Coulomb interaction. For such systems, the particles clearly have heat exchange with the enviroment so that the interacting particle system may be described by certain Langevin equations on the sphere:

    dX^i=V^i dt,
        m dV^i=-γ V^i dt+ℙ_S(1/N-1∑_j≠ iF(X^i-X^j) ) dt
    +√(2D) dB_S^i

Here, ℙ_S is the projection onto the tangent space of the sphere while B_S^i is the spherical Brownian motion to guarantee that the particle stays on the sphere. For theories of SDEs on manifolds, one may refer
to <cit.>. 
Here, D and γ must be related as in the classical fluctuation-dissipation theorem <cit.>. 

In the overdamped limit and with suitable scaling, we then have interacting particle system on the sphere (D_1=D/γ^2)

    dX^i=ℙ_S(1/N-1∑_j≠ iF(X^i-X^j)) dt+√(2D_1) dB_S^i

Numerically discretizing such SDEs on the sphere is an interesting topic which we leave for the future. In this work, 
we consider the Coulomb interaction with σ :=√(2D_1)=0, and use RBM-r as the stochastic strategy.  The following simple scheme for the SDE on the sphere is then adopted.


  * Randomly picking two indices. Then, solve the following for time t∈ [t_m-1, t_m)

    dX^i=∑_j: I(i,j)=1X^i-X^j/|X^i-X^j|^3 dt,

where I(i,j)=1 means that i,j are in the same batch. This can be solved analytically. In particular, defining 

    v̂_m=(X_m-1^i-X_m-1^j)/|X_m-1^i-X_m-1^j|,

one then has:

    X_m^i=(X_m-1^i+X_m-1^j)+v̂_m(|X_m-1^i-X_m-1^j|^3+6τ)^1/3,
       X_m^j=(X_m-1^i+X_m-1^j)-v̂_m(|X_m-1^i-X_m-1^j|^3+6τ)^1/3



  * Project the obtained points back to the sphere by dividing its magnitude.

The reason for setting σ=√(2D_1)=0 is that we would like to explore energy stable configurations. 
We desire low temperature regime for the ground state. Besides, the stochastic algorithm also introduces randomness so that we still have chance to get out of the local minimizers.







To check whether the method can give the desired ground state approximately, we randomly choose initial points on the sphere and run the above stochastic algorithms for enough iterations.
We do many experiments and check whether we always obtain the same final energy level.



 In Fig. <ref>, we show the numerical results in two experiments. 
 The number of particles are chosen as N=60 and N=800 respectively. The initial points are chosen randomly. The time step is chosen as τ=10^-4. As before, we regard the 'time' to be τ after N/2 iterations.
 For N=60, we see that in the eventual near stable configuration, each particle has 5 or 6 neighbors, and this agrees with the known results by physicists <cit.>. This configuration is quite different from the fullerene C_60 structure which is induced by the special properties of Carbon atoms.
For the N=800 case, the particles are roughly distributed uniformly. For both figures, there is only one stable energy level during the whole process. This means the system was only trapped in the final stable configuration.

 
 
 To check whether there are other possible stable energy configurations, we collect in Fig. <ref> the energies for N=800 after 3*(N/2)/τ iterations (T=3) in 25 experiments.  The simulation shows that one can find the ground state of the configuration almost surely using the stochastic algorithm and it is not easy to be trapped in local minimizers, if there are any.
  As studied by the physicists, there are many energy levels for the Thomson sphere. However, the numerical results here seem to suggest that the stochastic algorithms can obtain the ground state with high probability and the local minimizers of the energy landscape probably has small energy barriers. 
 Maybe, some interesting phenomena happen for large N's which needs further investigation.
 
 




 §.§ Two examples from economics and social science



In this section, we apply RBM-1 for two important models in social sciences, namely the evolution of wealth <cit.> and opinion dynamics <cit.>. The obtained stochastic processes not only are algorithms for the original models, but also can be viewed as new models which consider the fact that only a few individuals commute during a short time.



  §.§.§ Stochastic dynamics of wealth


We consider the model proposed by Degond et al. <cit.>, which tries to understand
the evolution of N market agents with two attributes: the economic configuration X^i and its wealth Y^j.

    Ẋ^i=V(X^i, Y^i),
        dY^i=-1/N-1∑_k:k≠ iξ_ikΨ(|X^i-X^k|)∂_yϕ(Y^i-Y^k) dt+√(2D) Y^i dB^i.

The first equation describes the evolution of the economic configuration, which is driven by the local Nash equilibrium and it is related to mean-field games <cit.>. The second equation describes the evolution of the wealth, which contains two mechanisms: the trading model proposed by Bouchaud and Mezart <cit.>, and the geometric Brownian motion in finance proposed by Bachelier in 1900 <cit.>. The quantity √(2D) is the volatility. The function ϕ is the trading interaction potential, while ξ_ikΨ(|X^i-X^k|) is the trading frequency. Often one assumes that ξ_ik depends on the number of trading agents in the economic neighborhoods of i and k:

    ξ_ik=ξ(ρ^i,Ψ+ρ^k,Ψ/2),  ρ^i,Ψ=1/N-1∑_ℓ≠ jΨ(|X^ℓ-X^i|).

The mean field Fokker-Planck equation is given by

    ∂_t f+∂_x(V(x, y)f)+∂_y(F f)=D∂_yy(y^2 f)

where 

    F(x, y, t)=-∫_x'≥ 0, y'≥ 0ξ(1/2(ρ^Ψ(x,t)+ρ^Ψ(x', t)))Ψ(|x-x'|) 
    ×∂_yϕ(y-y')f(x',y',t) dx' dy',

and 

    ρ^Ψ(x, t)=∫_x'>0, y'>0Ψ(|x-x'|)f(x', y', t)dx' dy'.


Now, if one considers the homogeneous case where the wealth dynamics is independent of the position in the economic configuration space, then Ψ is a constant. In this case, the dynamics is reduced to the interacting particle system, except that one has multiplicative noise

    dY^i=-κ/N-1∑_k:k≠ i∂_yϕ(Y^i-Y^k) dt+√(2D) Y^i dB^i,

where κ:=Ψξ(1/2(ρ^iΨ+ρ^k,Ψ)) is now a constant. The mean field equation is now given by

    ∂_tρ(y)=-∂_y(F(y)ρ(y))+D∂_yy(y^2ρ(y)),

where

    F(y)=-κ∫_y≥ 0∂_yϕ(y-y')ρ(y')dy'.

The equilibrium distribution is given by

    ρ_∞(y)∝exp(-α(y)/D),

where α satisfies

    ∂_yα(y)=-1/y^2F(y)+2D/y.



We now apply RBM-1 with p=2 to (<ref>) and have
for t∈ [t_m-1, t_m)

    dY^i=-κ∂_yϕ(Y^i-Y^θ) dt+√(2D) Y^i dB^i, i=1,…, N,

where θ is a random index that is grouped with i in the random division.
In some sense, the stochastic dynamics described by this algorithm can model what is happening in the real world: each agent only trades with a small number of random agents at a time. Then, after some time interval, the agents trade with others. Hence, (<ref>) is not just an algorithm but also it can be viewed as a new model.

For numerical test, choose the quadratic trading interaction as in <cit.>

    ϕ(y)=1/2y^2.

This interaction function may not be practical as it increases with y (intuitively, as y→∞, it should go to zero). The good thing is that with this interaction function, one can find the equilibrium distribution of wealth for (<ref>) as

    ρ_∞(y)=(κη/D)^κ/D+1/Γ(κ/D+1)y^-(2+κ/D)exp(-κη/Dy)1_y>0 .

This distribution is the inverse Gamma distribution and agrees with the Pareto power law for large y.
Here, η is the mean wealth. 

Now, we take κ=D=1 and consider the random dynamics (<ref>):

    dY^i=-  (Y^i-Y^θ) dt+√(2) Y^i dB^i, i=1,…, N .





In Fig. <ref>, we plot the empirical distribution of the wealth with N=10^5 agents. We choose τ=10^-3 and do the simulation to T=3. The SDE again is solved by the splitting strategy. Note that the splitting scheme preserves the mean wealth. For the test, we choose initial data as X^i=|Y^i| with Y^i∼𝒩(0, 1).
The reference curve is (<ref>) with 

    η=√(2)/√(π).

Clearly, the numerical results agree perfectly with the expected wealth distribution.



  §.§.§ Stochastic opinion dynamics

In this section, we consider some stochastic revisions of the opinion dynamics in <cit.>, where the following two models are mentioned


    d/dtX^i=α1/N-1∑_j≠ iϕ(|X^j-X^i|) (X^j-X^i)

and

    d/dtX^i=α∑_j≠ iϕ(|X^j-X^i|)/∑_kϕ(|X^k-X^i|) (X^j-X^i).


Here, ϕ is called the influence function.
These models are introduced for the emergence of consensus of opinions. Note that in the original
model, the prefactor is 1/N and we use 1/(N-1) here. The effect of chis change is minor.

Here, (<ref>) is not convenient for stochastic algorithms because of the denominator. Instead, we consider RBM-1 applied to (<ref>), and have the following:

    dX^i=αϕ(|X^θ-X^i|) (X^θ-X^i) dt
    +ϵ_N dB^i,

where θ, again as before, is a random index that is fixed for t∈ [t_m-1, t_m).  If ϵ_N=0, this is the stochastic algorithm of (<ref>) directly.
The parameter ϵ_N is to represent the random fluctuation on its opinions. If the number of agents is large, we believe this should be small.
In fact, this stochastic algorithm seems closer to what is happening in the world: each person only talks to one or few people at a time. Combining its previous opinion, it forms a new opinion.

Let us consider the influence function

    ϕ(r)=χ_[0, 1]

as in <cit.>.


In Fig. <ref>,  we show the numerical results. The figure shows the results of three experiments with the same initial data. We choose N=10^3 and τ=10^-4, α=40. 
Clearly, the final stationary consensus is sensitive with respect to initial distribution. With the same initial data, though the dynamics is stochastic, the main behavior is the same for the three experiments.  There are four main clusters of consensus. However, interestingly, in some experiments (like the third picture), there are may be some individuals that do not belong to any cluster, which seems to be the case in real world: some individuals are isolated at the early stage, and after the main clusters of consensus form, they are not affected by these groups since they are so different. The randomness introduced by the algorithms does not quite affect the main clusters of consensus, and only a few individuals might behave differently due to the randomness. After certain time, when the clusters of consensus are formed, the randomness does not play any roles any more: the individuals only talk to members in their own clusters.


In Fig. <ref>,  we show the numerical results for the stochastic opinion dynamics with Brownian motion
ϵ_N=1/N^1/3. The N, τ, α parameters are the same.
 The evolution of clusters of consensus is roughly the same with or without the Brownian motion. However, Brownian motions indeed introduce fluctuation of opinions within the clusters. This means that the fluctuation is not very important when the main clusters of opinions are formed. 










 §.§ Clustering through interacting particle system


In this section, we consider using the interacting particle system (<ref>) for clustering, as discussed in <cit.>. 
The idea is like this: consider N particles with a given adjacency matrix
A=(a_ij), a_ij≥ 0. Then, we construct the interacting particle system as

    d/dtX^i= α/N-1∑_j: j≠ i(a_ij-β)(X^j-X^i)

where α and β are some parameters. This is designed such that the particles with positive a_ij-β attract with each other so that they tend to gather together, while those with negative a_ij-β repel each other so that they separate. The hoping is that the intrinsic clusters will emerge automatically.

With the RBM-r, the computational cost is significantly reduced and this then becomes a practical method.
Each time, we pick a random set 𝒞 that contains p=2 elements. The dynamics we consider is then

    d/dtX^i=α (a_iθ-β)(X^θ-X^i)

where θ is again the random index that is constant on [t_m-1, t_m). Clearly, the random batch can be picked from those with nonzeros of a_ij only to improve the efficiency.


The update formula is 

    X_m^i=[(X_m-1^i+X_m-1^j)+(X_m-1^i-X_m-1^j)exp(-2α (a_ij-β)τ)]/2;
    
    X_m^j=[(X_m-1^i+X_m-1^j)-(X_m-1^i-X_m-1^j)exp(-2α  (a_ij-β)τ)]/2 .







  §.§.§ Cluster for stochastic block model


Let us consider the stochastic block model (<cit.>). The model is like this: suppose there are k clusters. For i, j in the same cluster, ℙ(a_ij=1)=p and ℙ(a_ij=0)=1-p. Otherwise, ℙ(a_ij=1)=q and ℙ(a_ij=0)=1-q. The entries a_ij are assumed to be independent. 
We assume that we only know the adjacent matrix in one experiment (if we know the matrices for several experiments, we can then combine them to get more accurate clusters). Clearly, the adjacent matrix is noisy. We are going to test whether or not we can still recover the clusters using the noisy adjacent matrix. 


In Fig. <ref>, we show the adjacent matrix from one experiment. In this example, we have N=1200 particles. The first 200 particles are designed to be in the first cluster, the next 400 are in the second cluster and the last 600 are in the third cluster. The probabilities are chosen as 

    p=0.7,  q=0.3.
 
In the matrix, yellow dots mean 1 while blue dots mean 0.
As expected,  most of the entries in the off-diagonal blocks are 0 with some yellow 'dust' scattered inside them.  
Most of the entries in the diagonal blocks are 1 with blue dots inside them.


In the experiment, we set β=1/2 and α=40, τ=10^-3.
We initialize their positions randomly on [0, 50]. 
The numerical results in an experiment are shown in Fig. <ref>. 
 Again the "time" is regarded as τ after N/2 iterations.
 From the figure, it is clear that the clusters can be recovered correctly though the adjacent matrix is noisy. 



  §.§.§ Reordering for sparse matrix


As another example, let us consider reordering sparse matrices as a byproduct of the clustering.
The point is that large a_ij entry tends to group the two indices together. 
If we use the terminal X_i's to sort, the reordered matrix can have large entries near the diagonal. If there are several distinct clusters, we will then have diagonal block matrix.

Consider the matrix given by

    A_1=BB^T+I,

where B is the 'west0479' matrix, which is a sparse matrix in the standard database of  MATLAB. Consequently, A_1 is a sparse matrix. Since matrix A_1 can have negative entries, we define A=(a_ij) with a_ij=|A_1(i,j)| to get a suitable adjacent matrix. Since A is sparse, we do sampling over the nonzero entries only.


Fig. <ref> shows the matrices min(P, 10) ("min" here means the entry-wise minimum), where the meaning of P is as follows: in the first figure, P is A; in the second figure, P is the reordered matrix using our strategy at T=5; P in the third figure is the matrix with approximated minimal degree ordering.
Clearly, our strategy gathers all the big entries near the diagonal, which means particles with strong interactions indeed form clusters. 


The ordering resulted from our clustering strategy works not that good as the minimal degree ordering using the criterion from sparse matrix theory. For example, in the Cholesky decomposition, there are 37896 nonzero entries. Meanwhile, by the approximated minimal degree ordering,  there are 14493 nonzero entries in the Cholesky decomposition.  However, the ordering using clustering gathers big entries near diagonals, which can yield better stability and may be advantageous for some applications.





§ CONCLUSIONS



We have developed Random Batch Methods for interacting particle systems with large number of particles and they reduce the computational cost significantly for N (N≫ 1) particles from O(N^J) to O(N) per time step. For RBM-1, the method without replacement, we have given a particle number independent error estimate under some special interactions. We have applied these methods to some representative problems of binary interactions in math, physics, social science and data science, and numerical results have supported our theory and expectations.
The random algorithms are powerful for systems with large number of individuals and high dimensions. 
 
As well accepted, in stochastic gradient descent, adding momentum could help to find flatter minimizers and improve results. In other words, the Langevin dynamics seems better for optimization and sampling <cit.>. Hence, considering the interacting particles with mass might be better for sampling the invariant measure of the nonlinear Fokker-Planck equation. This is left for future research.
Besides, there are many interesting projections ahead, for example,  proof of convergences for more general external and interacting potentials, and for RBMs with replacements.  It is also interesting to develop similar particle methods for the mean field equations, whenever they are available, as was done in <cit.> but for more general mean field equations.




§ ACKNOWLEDGEMENT

S. Jin was partially supported by the NSFC grant No. 31571071.  The work of L. Li was partially sponsored by NSFC 11901389, Shanghai Sailing Program 19YF1421300 and NSFC 11971314. The work of J.-G. Liu was partially supported by KI-Net NSF RNMS11-07444 and NSF DMS-1812573.

 





§ AN ERROR ANALYSIS FOR A HAMILTONIAN SYSTEM



In this section, we give a convergence proof of the RBM-1 for the following second order (Hamiltonian) system

    Ẋ^i=V^i,
       V̇^i= b(X^i)+1/N-1∑_j: j≠ iK(X^i-X^j),

where the dot means time derivative. Note that the variable 'V' here is reserved for the velocity so we use b(·) to represent the external force field corresponding to -∇ V in the main text.


Similarly, the processes generated by RBM-1 are determined by the following ODEs for t∈ [t_m-1, t_m)

    Ẋ̃̇^i=U^i,
       U̇^i= b(X̃^i)+1/p-1∑_j∈𝒞_q: j≠ iK(X̃^i-X̃^j),

where 𝒞_q is the batch in which i lives.


Assume that the initial data are drawn randomly and independently:

Systems (<ref>) and (<ref>) share the same initial data X^i(0)=X̃^i(0)=X_0^i and 
V^i(0)=U^i(0)=V_0^i, and the initial data (X_0^i, V_0^i) are i.i.d sampled from some common distribution.



With the setup, one can show the convergence of RBM-1 similarly.

Let Assumption <ref> hold. Suppose b(·)∈ C^1(ℝ^d) with b and ∇ b being bounded, and K is bounded and Lipschitz continuous. Then
the RBM-1 converges on [0, T] in the sense that 

    sup_0≤ t≤ T(𝔼1/N∑_j=1^N(|X^j(t)-X̃^j(t)|^2+|V^j(t)-U^j(t)|^2))^1/2   = 
    sup_0≤ t≤ T(𝔼(|X^i(t)-X̃^i(t)|^2+|V^i(t)-U^i(t)|^2))^1/2   ≤ C(T)√(τ/p-1),

for all i∈{1,…, N} and some C(T)>0 independent of N.



To prove this, one similarly introduces the error processes:

    Z^i=X̃^i-X^i,
        W^i=U^i-V^i



Clearly, the error processes satisfy for t∈ [t_m-1, t_m) that

    Ż^i =W^i,
       Ẇ^i = ( b(X̃^i)-b(X^i) )
    +1/N-1∑_j: j≠ iδ K_ij(t)
    +χ_m,i(X̃),

where again

    δ K_ij=K(X̃^i-X̃^j)-K(X^i-X^j),

and the error of the interaction force is given by

    χ_m,i(x):=1/p-1∑_j∈𝒞_θ,j≠ iK(x^i-x^j)
    -1/N-1∑_j: j≠ i K(x^i-x^j).

The notation 𝒞_θ is used as we want to emphasize that it is a random set over all experiments instead of the realization in one experiment (for which we use 𝒞_q).


The following lemmas are some preparation for the proof of the convergence, which are analogies of the lemmas in Section <ref> but are much easier due to the absence of Brownian motions.

For any q>1, the q-moments are bounded almost surely. More precisely, there exists C_q(T)>0 such that it holds almost surely that

    sup_0≤ t≤ T(|X̃^i(t)|^q+|U^i(t)|^q) ≤ C_q(T),

and 

    sup_0≤ t≤ T(|X^i(t)|^q+|V^i(t)|^q)≤ C_q(T).


The proof is very straightforward for which we omit the details. For example, regarding system (<ref>), it holds
almost surely that for t∈ [t_m-1, t_m]:

    d/dt(|X̃^i|^q+|U^i|^q)    =q|X̃^i|^q-2X̃^i· U^i+ q|U^i|^q-2 U^i·(b(X̃^i)+1/p-1∑_j∈𝒞_θ: j≠ iK(X̃^i-X̃^j))
       ≤ C|X̃^i|^q-1|U^i|+C|U^i|^q-1
       ≤ C_1(|X̃^i|^q+|U^i|^q)+C_2.

This will give the uniform bound of |X̃^i|^q+|U^i|^q on [0, T].

With the moments bound, it is easy to conclude the following for which we omit the proof.

For all t∈ [t_m-1, t_m),

    X̃^i(t)-X̃^i(t_m-1)+U^i(t)-U^i(t_m-1)≤ C(T)τ,

where · again means the L^2(Ω) norm.
Also, almost surely, it holds that

    |U^i(t)-U^i(t_m-1)|+|V^i(t)-V^i(t_m-1)|≤ C(T)τ.





Define again

    u(t)=𝔼1/N∑_j=1^N (|Z^j|^2+|W^j|^2).

Due to exchangeability, 

    u(t)=𝔼 (|Z^i|^2+|W^i|^2), ∀ i∈{1, …, N}.


For t∈ [t_m-1, t_m), taking time derivative, one has

    d/dtu=2𝔼(Z^i· W^i+W^i·(b(X̃^i)-b(X^i))
    +W^i·1/N-1∑_j: j≠ iδ K_ij)
    
    +2𝔼 W^i(t)·χ_m,i(Y(t)).

Since K is Lipschitz continuous, 

    𝔼(Z^i· W^i++W^i·(b(X̃^i)-b(X^i))
    +W^i·1/N-1∑_j: j≠ iδ K_ij)
       ≤ C𝔼(|Z^i| |W^i|)
    +C1/N-1∑_j: j≠ i𝔼|W^i||Z^j| 
       ≤ Cu(t).

The second inequality above is due to exchangeability.

Now, one can do the same trick as in Section <ref>:

    𝔼 W^i(t)·χ_m,i(X̃(t)) =𝔼W^i(t_m-1)·χ_m,i(X̃(t_m-1))
    +𝔼W^i(t_m-1)· (χ_m,i(X̃(t))-χ_m,i(X̃(t_m-1)))
    
    +𝔼(W^i(t)-W^i(t_m-1))·χ_m,i(X(t))
    
    +𝔼(W^i(t)-W^i(t_m-1))· (χ_m,i(X̃(t))-χ_m,i(X(t)))
    =: I_1+I_2+I_3+I_4.


Similarly, I_1=0, due to Lemma <ref>. Also,

    I_2≤W^i(t_m-1)1/p-1∑_j∈𝒞_θ, j≠ i |δ K^ij|
    +1/N-1∑_j: j≠ i |δ K^ij| ,

where

    δ K^ij:=K(X̃^i(t)-X̃^j(t))-K(X̃^i(t_m-1)-X̃^j(t_m-1)).

Clearly, due to the Lipschitz continuity of K and the equation of X̃^i, it holds

    |δ K^ij|≤ C(∫_t_m-1^t |U^i(s)| ds
    +∫_t_m-1^t |U^j(s)| ds )
    ≤ C(|U^i(t_m-1)|τ+|U^j(t_m-1)|τ+τ^2 ),

where the almost sure bound of U^i(t)-U^i(t_m-1) in Lemma <ref> has been used.
Applying Lemma <ref> regarding the norm of random sum, one therefore obtains

    I_2≤ CW^i(t_m-1)τ≤ C(W^i(t)τ+τ^2).


Similar as in Section <ref>, one has

    I_3≤ C√(u)τ+Λ_i_∞(1/p-1-1/N-1)τ.


Similarly as I_2, one can control I_4 by

    I_4≤W^i(t)-W^i(t_m-1)1/p-1∑_j∈𝒞_θ, j≠ i |δ K_ij|
    +1/N-1∑_j: j≠ i |δ K_ij| ,

where

    δ K_ij(t)=K(X̃^i(t)-X̃^j(t))-K(X^i(t)-X^j(t)).

Using the Lipschitz continuity of K again and the equation of X^i, X̃^i, one has

    |δ K_ij|≤  |X̃^i-X^i|(t_m-1)+|U^i(t_m-1)-V^i(t_m-1)|τ
    
    +|X̃^j-X^j|(t_m-1)+|U^j(t_m-1)-V^j(t_m-1)|τ+Cτ^2.

Since W^i(t)-W^i(t_m-1)≤ Cτ and applying again Lemma <ref>, one has

    I_4≤ C√(u)τ+Cτ^2.

The claim then follows by Grönwall inequality as in Section <ref>.



unsrt


