# From Massively Parallel Algorithms and Fluctuating Time Horizons to Non-equilibrium Surface Growth

## Abstract

We study the asymptotic scaling properties of a massively parallel algorithm for discrete-event simulations where the discrete events are Poisson arrivals. The evolution of the simulated time horizon is analogous to a non-equilibrium surface. Monte Carlo simulations and a coarse-grained approximation indicate that the macroscopic landscape in the steady state is governed by the Edwards-Wilkinson Hamiltonian. Since the efficiency of the algorithm corresponds to the density of local minima in the associated surface, our results imply that the algorithm is asymptotically scalable.

To efficiently utilize modern supercomputers requires massively parallel implementations of dynamic algorithms for various physical, chemical, and biological processes. For many of these there are well-known and routinely used serial Monte Carlo (MC) schemes which are based on the realistic assumption that attempts to update the state of the system form a Poisson process. The parallel implementation of these dynamic MC algorithms belongs to the class of parallel discrete-event simulations, which is one of the most challenging areas in parallel computing  and has numerous applications not only in the physical sciences, but also in computer science, queueing theory, and economics. For example, in lattice Ising models the discrete events are spin-flip attempts, while in queueing systems they are job arrivals. Since current special- or multi-purpose parallel computers can have $`10^410^5`$ processing elements (PE) , it is essential to understand and estimate the scaling properties of these algorithms.

In this Letter we introduce an approach to investigate the asymptotic scaling properties of an extremely robust parallel scheme . This parallel algorithm is applicable to a wide range of stochastic cellular automata with local dynamics, where the discrete events are Poisson arrivals. Although attempts have been made to estimate its efficiency under some restrictive assumptions , the mechanism which ensures the scalability of the algorithm in the “steady state” was never identified. Here we accomplish this by noting that the simulated time horizon is analogous to a growing and fluctuating surface. The local random time increments correspond to the deposition of random amounts of “material” at the local minima of the surface. This correspondence provides a natural ground for cross-disciplinary application of well-known concepts from non-equilibrium surface growth  and driven systems  to a certain class of massively parallel computational schemes. To estimate the efficiency of this algorithm one must understand the morphology of the surface associated with the simulated time horizon. In particular, the efficiency of this parallel implementation (the fraction of the non-idling processing elements) exactly corresponds to the density of local minima in the surface model. We show that the steady-state behavior of the macroscopic landscape is governed by the Edwards-Wilkinson (EW) Hamiltonian , implying that the density of the local minima does not vanish when the number of PEs goes to infinity. This ensures that the simulated time horizon propagates with a non-zero average velocity in the steady state. Thus the algorithm is asymptotically scalable! Further, based on the strong analogy between the evolution of the simulated time horizon and the single-step surface growth model , we describe the asymptotic scaling properties of the parallel scheme.

The difficulty of parallel discrete-event simulations is that the discrete events (update attempts) are not synchronized by a global clock. The traditional dynamic MC algorithms were long believed to be inherently serial, i.e., in spin language, the corresponding algorithm could attempt to update only one spin at a time. Lubachevsky nevertheless presented an approach for parallel simulation of these systems  without changing the dynamics of the underlying model . Applications of his scheme include modeling of cellular communication networks , ballistic particle deposition , and metastability and hysteresis in kinetic Ising models .

Here we consider the case of one-dimensional systems with only nearest-neighbor interactions (e.g., Glauber spin-flip dynamics) and periodic boundary conditions. We restrict ourselves to the case where each PE carries one site (e.g., one spin) of the underlying system. Non-zero efficiency for this algorithm implies non-zero efficiency for the case where a PE carries a block of sites . Also, the stochastic model for the simulated time horizon is an exact mapping only for the synchronous algorithm, where the main simulation cycles are executed in lock-step on each PE. Our goal is to show for this one-site-per-PE synchronous algorithm (which can be regarded as the worst-case scenario) that the efficiency does not go to zero as the number of PEs goes to infinity. The basic synchronous parallel scheme  is as follows.

The size of the underlying system, and thus the number of PEs, is $`L`$. Update attempts at each site are independent Poisson processes with the same rate, independent of the state of the underlying system. Hence, the random time interval between two successive attempts is exponentially distributed. Without loss of generality we use time increments of mean one (in arbitrary units). The Poisson arrivals correspond to attempted instantaneous changes in the state of the site. In the parallel algorithm each PE generates its own local simulated time for the next update attempt, denoted by $`\tau _i(t)`$, $`i=1,2,\mathrm{},L`$. Here, $`t`$ is the discrete index of the parallel steps simultaneously performed by each PE. Initially $`\tau _i(0)=0`$ for each site, and an initial configuration for the underlying system is specified. The simulated time of the first update attempt is determined by $`\tau _i(1)=\tau _i(0)+\eta _i(0)`$, where $`\{\eta _i\}`$ are independent exponential variables. For parallel steps $`t1`$, each PE must compare its local simulated time to the local simulated times of its neighbors. If $`\tau _i(t)\mathrm{min}\{\tau _{i1}(t),\tau _{i+1}(t)\}`$, the change of state of the site is attempted (and decided by the rules of the underlying system), and its local simulated time is incremented by an exponentially distributed random amount, $`\tau _i(t+1)=\tau _i(t)+\eta _i(t)`$. Otherwise, the change of state is not attempted and the local simulated time remains the same, $`\tau _i(t+1)=\tau _i(t)`$, i.e., the PE waits (“idles”). The comparison of the nearest-neighbor simulated times, and waiting if necessary, ensures that information passed between PEs does not violate causality. The algorithm is obviously free from deadlock, since at worst the PE with the absolute minimum local time can make progress. After the initial step, the probability density of the simulated time horizon $`\{\tau _i(t)\}`$ is a continuous measure, so the probability that updates for two nearest-neighbor sites are attempted at the same simulated time, is of measure zero. When modeling the efficiency, we ignore communication times between PEs, since they typically contribute to a scalable overhead. Thus, the efficiency is simply the fraction of non-idling PEs (inherent utilization). This exactly corresponds to the density of local minima of the simulated time horizon.

The question naturally arises: Is it possible that the fraction of non-idling PEs goes to zero in the $`L`$$``$$`\mathrm{}`$ limit? This would obviously make the algorithm unscalable and the performance of the actual implementation poor if not disastrous. To study this problem, we focus on the evolution of the simulated time horizon $`\{\tau _i(t)\}`$, which is completely independent of the underlying model. The above algorithmic steps can be compactly summarized as:

$`\tau _i(t+1)`$ $`=`$ $`\tau _i(t)`$ (1)   
$`+`$ $`\mathrm{\Theta }\left(\tau _{i1}(t)\tau _i(t)\right)\mathrm{\Theta }\left(\tau _{i+1}(t)\tau _i(t)\right)\eta _i(t).`$ (2)   
The $`\eta _i(t)`$ are drawn from an exponential distribution independently at every time $`t`$ and site $`i`$, and independent of $`\{\tau _i(t)\}`$. Here $`\mathrm{\Theta }()`$ is the Heaviside step-function. This stochastic evolution model is very simple and easily simulated on a serial computer. Alternatively, one can consider the evolution of the local slopes, $`\varphi _i=\tau _i\tau _{i1}`$:

$`\varphi _i(t+1)\varphi _i(t)`$ $`=`$ $`\mathrm{\Theta }\left(\varphi _i(t)\right)\mathrm{\Theta }\left(\varphi _{i+1}(t)\right)\eta _i(t)`$ (3)   
$``$ $`\mathrm{\Theta }\left(\varphi _{i1}(t)\right)\mathrm{\Theta }\left(\varphi _i(t)\right)\eta _{i1}(t).`$ (4)   
The periodic boundary conditions for $`\{\tau _i\}`$ impose the constraint $`_{i=1}^L\varphi _i=0`$. In this representation the operator for the density of local minima is

$$u(t)=\frac{1}{L}\underset{i=1}{\overset{L}{}}\mathrm{\Theta }(\varphi _i(t))\mathrm{\Theta }(\varphi _{i+1}(t)).$$
(5)   
The process described by (3) is a microscopic realization of biased diffusion . The random amount of material “deposited” at a local minimum $`\tau _i`$ corresponds to the transfer of this amount from $`\varphi _{i+1}`$ to $`\varphi _i`$. Since the noise in (3) is independent of $`\{\varphi _i(t)\}`$, the average can be simply taken. This yields a transparent continuity equation, $`\varphi _i(t+1)\varphi _i(t)=(j_ij_{i1})`$, where the average current is $`j_i=\mathrm{\Theta }(\varphi _i)\mathrm{\Theta }(\varphi _{i+1})`$. Translational invariance implies that $`u=\mathrm{\Theta }(\varphi _i)\mathrm{\Theta }(\varphi _{i+1})`$, which is the same as the magnitude of the average current or the mean velocity of the surface.

To gain some insight into the evolution of the surface, we perform a naive coarse-graining by taking an ensemble average on (3) and replacing $`\mathrm{\Theta }(\varphi )`$ with a smooth representation. The procedure is independent of the actual form of the representation. We use $`\mathrm{\Theta }_\kappa (\varphi )=(1/2)[\mathrm{tanh}(\varphi /\kappa )+1]`$ so $`lim_{\kappa 0}\mathrm{\Theta }_\kappa (\varphi )=\mathrm{\Theta }(\varphi )`$. To leading order in $`\varphi /\kappa `$,

$`\varphi _i(t+1)\varphi _i(t)`$ $`=`$ $`{\displaystyle \frac{1}{4\kappa }}\varphi _{i+1}2\varphi _i+\varphi _{i1}`$ (6)   
$``$ $`{\displaystyle \frac{1}{4\kappa ^2}}\varphi _i(\varphi _{i+1}\varphi _{i1}).`$ (7)   
Taking the naive continuum limit one obtains

$$_t\widehat{\varphi }=\frac{^2\widehat{\varphi }}{x^2}\lambda \frac{}{x}\widehat{\varphi }^2$$
(8)   
for the coarse-grained field, $`\widehat{\varphi }`$, where roughly speaking $`\lambda `$, the coefficient of the nonlinear term, carries the details of the coarse-graining procedure. Equation (8) is the nonlinear biased diffusion or Burgers’ equation . Through $`\widehat{\varphi }=\widehat{\tau }/x`$ it is simply related to the deterministic part of the KPZ equation for the coarse-grained surface height fluctuations ,

$$_t\widehat{\tau }=\frac{^2\widehat{\tau }}{x^2}\lambda \left(\frac{\widehat{\tau }}{x}\right)^2.$$
(9)   
To capture the fluctuations one typically extends the above equations with appropriate noise, i.e., conserved for (8), and non-conserved for (9). This implies that the evolution of the simulated time horizon is KPZ-like. In one dimension the steady state of such systems (on coarse-grained length scales) is governed by the EW Hamiltonian , $`_{\mathrm{EW}}`$$``$$`𝑑x(\widehat{\tau }/x)^2`$. This corresponds to a simple random-walk surface, where the coarse-grained slopes are independent in the thermodynamic limit, yielding $`1/4`$ for the density of local minima. Obviously this value will be different for our specific microscopic model. However it cannot vanish: a zero density of local minima in the $`L`$$``$$`\mathrm{}`$ limit would imply that it is zero at all length scales. This would contradict our finding that the steady state at the coarse-grained level is governed by the EW Hamiltonian. The non-zero density of local minima is an important characteristic of this (steady-state) universality class. It ensures that our specific microscopic surface propagates with a non-zero average velocity in the steady state. Models belonging to other universality classes do not necessarily have non-zero extremal-point densities. For example, we can show  that the density of local minima vanishes for a one-dimensional curvature-driven random Gaussian surface.

We now present our MC results to test the coarse-graining approach. First we follow the time evolution of the width $`w^2(t)`$$`=`$$`(1/L)_{i=1}^L\left[\tau _i(t)\overline{\tau }(t)\right]^2`$, where $`\overline{\tau }(t)`$$`=`$$`(1/L)_{i=1}^L\tau _i(t)`$ and the average $``$ is taken over many independent runs. After the early-time regime, which is strongly affected by the intrinsic width, and before saturation, we find $`w^2(t)`$$``$$`t^{2\beta }`$ \[Fig. 1(a)\]. Although the system exhibits very strong corrections to scaling, for our largest system, $`L`$$`=`$$`10^5`$, we find $`\beta `$$`=`$$`0.326\pm 0.005`$, which includes within two standard errors the exact KPZ exponent, $`1/3`$. In the steady state the width is stationary and $`w^2`$$``$$`L^{2\alpha }`$ for large $`L`$. Here the corrections to scaling are somewhat smaller than in the earlier regime, and the above scaling is obeyed for $`L`$$``$$`10^3`$ with $`\alpha `$$`=`$$`0.49\pm 0.01`$. This agrees with the prediction that the long-distance behavior is governed by the EW Hamiltonian, for which $`\alpha `$$`=`$$`1/2`$. Further, plotting rescaled width $`w^2/L^{2\alpha }`$ vs rescaled time $`t/L^z`$, with $`z`$$`=`$$`\alpha /\beta `$, confirms dynamic scaling for the intermediate-to-late time crossover  \[Fig. 1(a) inset\]. We also measured the average steady-state structure factor of $`\{\tau _i\}`$, finding the expected $`1/k^2`$ behavior for small wave vectors $`k`$. The spatial two-point correlations decay linearly for $`\{\tau _i\}`$ and are short ranged for the slopes, $`\{\varphi _i\}`$.

To further probe the universal properties of the surface in the steady state we construct the full width distribution, $`P(w^2)`$. The EW class is characterized by a universal scaling function, $`\mathrm{\Phi }(x)`$, such that $`P(w^2)=w^2^1\mathrm{\Phi }(w^2/w^2)`$ , where

$$\mathrm{\Phi }(x)=\frac{\pi ^2}{3}\underset{n=1}{\overset{\mathrm{}}{}}(1)^{n1}n^2e^{\frac{\pi ^2}{6}n^2x}.$$
(10)   
Systems with $`L`$$``$$`10^3`$ show convincing data collapse \[Fig. 1(b)\] onto this exact scaling function.

Next we estimate the efficiency of the algorithm. When simulating the system described by (1) and measuring the average local minimum density $`u(t)`$, we observe that for every system size it monotonically decreases as a function of time and approaches a constant slightly smaller than $`1/4`$ for large systems. Using the close similarity between our model and the single-step solid-on-solid surface-growth model , we can understand the scaling behavior for the steady-state average $`u`$ and the fluctuations $`\sigma ^2=u^2u^2`$. In the single-step model the height differences (i.e., the local slopes) are restricted to $`\pm 1`$, and the evolution consists of particles of height $`2`$ being deposited at the local minima. The advantage of this model is that it can be mapped onto a hard-core lattice gas for which the steady-state probability distribution of the configurations is known exactly . This enables one to find arbitrary moments of the local minimum density operator, analogous to (5). For the single-step model we find that $`u_L=(1/4)(11/L)^1=1/4+1/(4L)+𝒪(1/L^2)`$, and $`\sigma _L^2=1/(16L)+𝒪(1/L^2)`$. We propose that the scaling of the local minimum density for large $`L`$ in our model follows the same form, i.e.,

$$u_Lu_{\mathrm{}}L^1,\sigma _LL^{1/2}.$$
(11)   
Our reasons are: (i) both models in their steady states belong to the EW universality class (short-range correlated local slopes), which guarantees that $`u_{\mathrm{}}`$ is non-zero; (ii) the constraint $`_{i=1}^L\varphi _i=0`$ in our model and the conservation of the number of particles in the lattice gas will produce similar finite-size effects. Our simulation results show very good agreement with (11) \[Fig. 2(a)\]. Further, extrapolating to $`L`$$``$$`\mathrm{}`$, yields $`u_{\mathrm{}}`$$`=`$$`0.24641\pm (7\times 10^6)`$. The scaling relation (11) also implies that $`u`$ is a self-averaging macroscopic quantity: its full distribution $`P_L(u)`$, for large $`L`$, is Gausian with the parameters in (11) \[Fig. 2(b)\]. Thus, in the $`L`$$``$$`\mathrm{}`$ limit it approaches a delta-function centered about $`u_{\mathrm{}}`$.

Finally, we point out the lack of up-down symmetry in our model in the steady-state. This is most easily noticeable at short distances, either by looking at a snapshot \[Fig. 3(a)\], or through the high degree of asymmetry in the nearest-neighbor two-slope distribution \[Fig. 3(b)\]: the hilltops are sharp and the valley-bottoms are flattened. Such stationary-state skewness is generally observable in one-dimensional KPZ growth, but has only recently received serious attention .

In summary, we studied the asymptotic scaling properties of a general parallel algorithm by regarding the simulated time horizon as a non-equilibrium surface. We conclude that the basic algorithm (one site per PE) is scalable for one-dimensional arrays. The same correspondence can be applied to model the performance of the algorithm for higher-dimensional logical PE topologies. While this will involve the typical difficulties of surface-growth modeling, such as an absence of exact results and very long simulation times, it establishes potentially fruitful connections between two traditionally separate research areas.

We thank S. Das Sarma, S. J. Mitchell, and G. Brown for stimulating discussion. We acknowledge the support of DOE through SCRI-FSU, NSF-MRSEC at UMD, and NSF through Grant No. DMR-9871455.