  The single-site dynamics are a canonical class of Markov chains for sampling
from high-dimensional probability distributions, e.g. the ones represented by
graphical models.
  We give a simple and generic parallel algorithm that can faithfully simulate
single-site dynamics. When the chain asymptotically satisfies the
$\ell_p$-Dobrushin's condition, specifically, when the Dobrushin's influence
matrix has constantly bounded $\ell_p$-induced operator norm for an arbitrary
$p\in[1,\infty]$, the parallel simulation of $N$ steps of single-site updates
succeeds within $O\left({N}/{n}+\log n\right)$ depth of parallel computing
using $\tilde{O}(m)$ processors, where $n$ is the number of sites and $m$ is
the size of graphical model. Since the Dobrushin's condition is almost always
satisfied asymptotically by mixing chains, this parallel simulation algorithm
essentially transforms single-site dynamics with optimal $O(n\log n)$ mixing
time to RNC algorithms for sampling. In particular we obtain RNC samplers, for
the Ising models on general graphs in the uniqueness regime, and for satisfying
solutions of CNF formulas in a local lemma regime. With non-adaptive simulated
annealing, these RNC samplers can be transformed routinely to RNC algorithms
for approximate counting.
  A key step in our parallel simulation algorithm, is a so-called "universal
coupling" procedure, which tries to simultaneously couple all distributions
over the same sample space. We construct such a universal coupling, that for
every pair of distributions the coupled probability is at least their Jaccard
similarity. We also prove this is optimal in the worst case. The universal
coupling and its applications are of independent interests.
