  An efficient sampling method, the pmmLang+RBM, is proposed to compute the
quantum thermal average in the interacting quantum particle system. Benefiting
from the random batch method (RBM), the pmmLang+RBM reduces the complexity due
to the interaction forces per timestep from $O(NP^2)$ to $O(NP)$, where $N$ is
the number of beads and $P$ is the number of particles. Although the RBM
introduces a random perturbation of the interaction forces at each timestep,
the long time effects of the random perturbations along the sampling process
only result in a small bias in the empirical measure of the pmmLang+RBM from
the target distribution, which also implies a small error in the thermal
average calculation. We numerically study the convergence of the pmmLang+RBM,
and quantitatively investigate the dependence of the error in computing the
thermal average on the parameters including the batch size, the timestep, etc.
We also propose an extension of the pmmLang+RBM, which is based on the
splitting Monte Carlo method and is applicable when the interacting potential
contains a singular part.
