Random Batch Method (RBM)

• Random Batch Method (RBM) is a Monte Carlo-type algorithm that replaces full interaction with random sub-batch interactions, ensuring unbiased force estimation.
• It reduces computational complexity from O(N²) to O(N) by computing interactions only among small, randomly chosen batches.
• RBM offers quantitative error estimates and has been applied in classical, quantum, and sampling settings to improve computational efficiency.

Random Batch Method (RBM) denotes a family of Monte-Carlo-type algorithms for large interacting systems in which the all-to-all interaction graph is replaced, on each short time interval, by a randomly generated sparse interaction graph on small batches. In the standard binary-interaction setting, one randomly partitions $N$ particles into batches of size $p\ll N$, computes only intra-batch interactions, and rescales the result so that the batch force is an unbiased estimator of the full mean-field force. This changes the per-step cost from $O(N^2)$ to $O(pN)$, hence to $O(N)$ when $p$ is fixed, while preserving the correct averaged interaction and admitting quantitative error estimates in several classical, quantum, and sampling settings (Jin et al., 2018, Jin et al., 2021).

1. Core formulation

A canonical first-order interacting particle system treated in the RBM literature is

$$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$

The computational bottleneck is the pairwise interaction term, whose direct evaluation costs $O(N^2)$ per time step because each particle interacts with all $N-1$ others (Jin et al., 2018).

RBM replaces the full interaction by a batch interaction. In the without-replacement formulation, assuming $N=np$, one randomly divides $p\ll N$0 into disjoint batches $p\ll N$1 of size $p\ll N$2, and for $p\ll N$3 evolves

$p\ll N$4

The full force

$p\ll N$5

is thus replaced by the random-batch force

$p\ll N$6

For a fixed configuration $p\ll N$7 independent of the random partition, the fluctuation

$p\ll N$8

satisfies

$p\ll N$9

and

$O(N^2)$0

with

$O(N^2)$1

Unbiasedness is the central consistency property, and the variance formula makes explicit how increasing $O(N^2)$2 reduces the random-batch fluctuation (Jin et al., 2018).

In applications to sampling and inference, the term “batch” has a specialized meaning: it is a batch of interacting particles, not a minibatch of data for estimating gradients of a loss. That distinction is explicit in the RBM-SVGD literature and is important for interpreting the source of stochasticity in RBM-based methods (Li et al., 2019).

2. Algorithmic realizations and computational structure

The standard implementation, often called RBM-1 or RBM without replacement, redraws a full random partition at every time step and updates every particle using only its current batchmates. Because each particle interacts with only $O(N^2)$3 others, the work becomes $O(N^2)$4, and the reshuffling itself can be done in $O(N^2)$5 time using a standard permutation algorithm (Jin et al., 2018, Golse et al., 2019).

A second realization is RBM-r, or RBM with replacement. In this variant, one repeatedly picks a single batch $O(N^2)$6 of size $O(N^2)$7, lets only the particles in that batch interact for a short time, freezes the others during that substep, and repeats with independently chosen batches. For first-order pairwise systems, this variant is identified with kinetic Monte Carlo. Its convergence is formulated after the natural pseudo-time rescaling

$O(N^2)$8

reflecting the fact that only a fraction $O(N^2)$9 of the particles is active at each substep (Cai et al., 2024).

A third realization appears in the quantum $O(pN)$0-body setting. There the interaction of each particle with all $O(pN)$1 others is replaced by interaction with $O(pN)$2 particles chosen at random at each time step, multiplied by

$O(pN)$3

For simplicity, the quantum analysis of (Golse et al., 2019) treats only the case $O(pN)$4: $O(pN)$5 is assumed even, particles are organized into $O(pN)$6 random pairs, and the pairs are reshuffled independently at the beginning of every step. This “random pairs” formulation is the batch analogue of a fully coupled quantum Hamiltonian and preserves the standard RBM philosophy of using a small random subset of interactions instead of all pairwise forces (Golse et al., 2019).

These variants share the same structural principle but lead to distinct analytical regimes. Without replacement, every particle is updated at each step and long-time bounds are typically sharper. With replacement, the analysis must handle random active times and often requires auxiliary time-change constructions. Pair-based quantum RBM adds a further constraint: the comparison is not performed in trace norm but through reduced density matrices and Wigner transforms (Cai et al., 2024, Golse et al., 2019).

3. Convergence theory

The first general strong-convergence theorem for RBM-1 was established under a contractive setting with strongly convex confinement and bounded Lipschitz interaction. If

$O(pN)$7

then

$O(pN)$8

and consequently

$O(pN)$9

The constant is independent of $O(N)$0, $O(N)$1, and $O(N)$2, which is the core particle-number-independent stability statement behind the claim that RBM is asymptotic-preserving for large interacting systems (Jin et al., 2018).

For second-order interacting particle systems, analogous uniform-in-time strong convergence holds under hypocoercive contraction assumptions. In the mean-field scaling regime,

$O(N)$3

with $O(N)$4 independent of $O(N)$5 and $O(N)$6. The proof uses modified energies and a decomposition of the interaction discrepancy into Lipschitz terms and batch fluctuations (Jin et al., 2020).

A law-level improvement appears in the mean-field error analysis of RBM toward the nonlinear Fokker–Planck limit. For the joint law $O(N)$7 of the RBM particles and the tensorized mean-field law $O(N)$8, the rescaled relative entropy satisfies

$O(N)$9

By a Csiszár–Kullback–Pinsker type inequality this yields

$p$0

which improves the discretization-step rate from $p$1 to $p$2 in Wasserstein distance (Huang et al., 2024).

The mean-field analysis of the Cucker–Smale model further separates the limits $p$3 and $p$4. In that setting,

$p$5

for fixed $p$6, and

$p$7

The analysis replaces continuous-time propagation of chaos by discrete-time chaos propagation through the RBM flux map (Wang et al., 2024).

Quantum RBM requires a different comparison principle. For the full density operator $p$8 and the RBM density operator $p$9, the comparison is made at the level of the single-particle reduced density operators and their Wigner transforms. The main estimate in (Golse et al., 2019) shows that

$$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$0

is bounded by an $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$1 expression whose right-hand side does not depend on $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$2 or $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$3. This provides a weak phase-space convergence result that is uniform in both particle number and semiclassical parameter (Golse et al., 2019).

4. Variants, corrections, and generalized noise models

The basic RBM construction has generated several specialized variants designed to control variance, accommodate non-Gaussian forcing, or stabilize singular interactions.

Variant	Setting	Key statement
RBM-Lévy	Interacting particle systems with Lévy noise	Cost drops from $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$4 to $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$5; $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$6 under finite second moment, and $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$7 under finite first moment (Liu et al., 2024)
RBM-M	Strongly singular interaction kernels	Uses $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$8; the corrected interaction is asymptotically unbiased and has smaller variance up to lower-order terms (Zhao et al., 2024)
rvRBM	Nonlocal PDEs of Fokker–Planck type	Adds a control-variate correction from a cheap surrogate model with the same asymptotic state as the original interaction (Pareschi et al., 2023)
VR-RBL	Langevin dynamics	Estimates batch-force variance and subtracts it from the Brownian covariance to remove artificial heating and obtain a local truncation error of $$d X^i=-\nabla V(X^i)\,dt+\frac{1}{N-1}\sum_{j: j\neq i} K(X^i-X^j)\,dt+\sigma dB^i,\quad i=1,\ldots, N.$$9 in expectation (Xu et al., 2024)
Split-kernel RBM	Second-order systems and molecular dynamics	Decomposes $O(N^2)$0, computes the short-range singular part exactly, and applies RBM only to the smooth long-range part (Jin et al., 2020)

RBM-Lévy extends the shuffle-and-interact mechanism to jump-driven systems. The extension is nontrivial because the proof must handle discontinuous trajectories and, depending on the Lévy measure, possibly infinite second moments. The resulting Wasserstein estimates are uniform in time, and the authors explicitly connect the order of the Wasserstein metric to the maximal finite moment of the Lévy process (Liu et al., 2024).

RBM-M addresses a distinct failure mode: for strongly singular kernels, standard RBM can miss or overrepresent rare but strong interactions, so the batch estimator becomes unstable. The momentum-like exponential averaging

$O(N^2)$1

retains information from previous batch samplings and yields a smaller upper bound than standard RBM under the smooth assumptions used in the theorem. The paper emphasizes, however, that $O(N^2)$2 should typically be small because particle distributions change at every step (Zhao et al., 2024).

Variance-reduction strategies split into two main designs. In nonlocal kinetic equations, rvRBM subtracts and adds a cheap surrogate interaction $O(N^2)$3 and chooses the control-variate coefficient $O(N^2)$4 by variance minimization. In Langevin dynamics, VR-RBL estimates the covariance of the random batch force and modifies the Brownian increment so that fluctuation-dissipation balance is preserved to higher order. Both approaches keep the linear-scaling advantage of RBM while explicitly targeting the batching-induced fluctuation term (Pareschi et al., 2023, Xu et al., 2024).

5. Applications across computation, sampling, and PDEs

RBM has been deployed across a broad spectrum of interacting-system models. In Bayesian inference and particle-based variational sampling, RBM-SVGD applies batching to the SVGD particle ODE

$O(N^2)$5

reducing each iteration from $O(N^2)$6 to $O(N^2)$7 while retaining unbiasedness of the interaction term. For the continuous-time RBM-SVGD on the torus, the finite-time theorem gives

$O(N^2)$8

showing that the one-particle marginal converges to that of SVGD as the step size $O(N^2)$9 (Li et al., 2019).

In quantum Monte Carlo, RBM was introduced to alleviate the cost of two-body interactions in both variational and diffusion Monte Carlo. For the overdamped Langevin formulation of VMC, updating the position of each particle requires only $N-1$0 operations, and the per-step cost for $N-1$1 particles is reduced from $N-1$2 to $N-1$3. In path-integral quantum thermal sampling, pmmLang+RBM reduces the interaction-force complexity per timestep from $N-1$4 to $N-1$5, where $N-1$6 is the number of beads and $N-1$7 is the number of particles (Jin et al., 2020, Ye et al., 2021).

In collective behavior and control, RBM has been used for flocking, consensus, and guiding problems. For the guiding model of evaders and repulsive drivers, the all-to-all evader interaction is replaced by a sequence of random sparse interaction graphs. The resulting surrogate is embedded in gradient-based optimal control and then in a model predictive control strategy, producing a semi-feedback controller for the full system while reducing the forward cost from $N-1$8 to $N-1$9 (Ko et al., 2020).

In stochastic-statistical closure for multiscale turbulence, RBM is used to batch the interaction of fluctuation modes in high-order moment equations. The direct stochastic forecast cost $N=np$0 and covariance-evolution cost $N=np$1 are reduced to $N=np$2 and $N=np$3, respectively, with $N=np$4. The reduced-order version links many small-scale fluctuation modes to ensemble samples of dominant leading modes and is used in one-layer and two-layer Lorenz ’96 systems (Qi et al., 2023).

RBM has also moved beyond particle systems into PDEs on graphs. One line of work applies RBM directly at the PDE level, before any space or time discretization, by randomly activating only one subgraph batch on each short time window. For parabolic equations on graphs, the mean-square error is first order in the RBM step size and uniform in time (Hernández, 29 Aug 2025). A second line follows a discretize+RBM strategy: first discretize the heat equation on a graph by finite elements, then decompose the discrete operator into graph-based blocks and randomly activate them. For the state approximation,

$N=np$5

and for the optimal-control version the RBM switching parameter must be chosen more restrictively, with $N=np$6 in the convergence regime described in the paper (Hernández et al., 13 Jun 2025).

6. Induced stochasticity, limitations, and interpretation

RBM is not a neutral sparsification of the interaction graph. Because the batch force is only an estimator of the full force, the method injects additional stochasticity into the dynamics. In two toy models with mean-field phase transition—the Curie–Weiss model and a double-well McKean–Vlasov system—the effective dynamics show that RBM adds extra noise, does not destroy the phase transition, but shifts the critical temperature downward, equivalently increasing the critical inverse temperature (Guillin et al., 2023).

That induced fluctuation can be benign, useful, or problematic depending on the regime. The foundational RBM paper explicitly notes that the stochastic processes generated by the algorithms can also be regarded as new models for the underlying problems, especially in social or economic systems where agents interact only with a few random peers at a time (Jin et al., 2018). In contrast, Langevin simulations expose a more delicate effect: the variance of the random batch force acts like an additional noise source and can produce artificial heating unless corrected by a fluctuation-dissipation-consistent variance estimate (Xu et al., 2024).

Several analytical limitations remain model-dependent. In the quantum $N=np$7-body theory of (Golse et al., 2019), the theorem is for the expected RBM dynamics rather than for a fixed random realization, and it is proved only for the simplest pair-batch variant with $N=np$8. In the Lévy setting, the analysis does not yet cover space-inhomogeneous noise, singular interaction kernels, or multiplicative jump noise (Liu et al., 2024). For SVGD, the finite-time convergence proof is available, but stronger long-time bounds would require contraction properties of the SVGD flow (Li et al., 2019). For the guiding problem, the RBM-MPC strategy is supported heuristically and numerically, while a rigorous error theory in that nonlinear control setting remains open (Ko et al., 2020).

The with-replacement variant illustrates a further tradeoff. For first-order systems it admits an explicit Wasserstein-2 convergence rate after pseudo-time rescaling, and for the deterministic Cucker–Smale model it admits a uniform-in-time, uniform-in-$N=np$9 error bound. But the Cucker–Smale analysis also shows that RBM-r has a weaker exponential decay factor than RBM-1, so replacement is analytically and empirically less favorable in that setting (Cai et al., 2024, Wang et al., 25 Jan 2025).

Taken together, these results support a precise interpretation: RBM is a unifying stochastic approximation paradigm for high-dimensional interaction operators. Its central benefit is linear-scaling interaction evaluation; its central cost is the artificial fluctuation produced by batching. Theoretical work on RBM therefore revolves around three recurring themes—unbiasedness of the batch estimator, quantitative control of the batch variance, and structural corrections when the added randomness interacts unfavorably with singular kernels, jump noise, phase transitions, or fluctuation-dissipation balance (Jin et al., 2021, Xu et al., 2024).