Parallel Simulation of Two-Dimensional Ising Models Using Probabilistic Cellular Automata

Abstract

We study the numerical simulation of the shaken dynamics, a parallel Markovian dynamics for spin systems with local interaction and transition probabilities depending on the two parameters q and J that “tune” the geometry of the underlying lattice. The analysis of the mixing time of the Markov chain and the evaluation of the spin-spin correlations as functions of q and J, make it possible to determine in the (q, J) plane a phase transition curve separating the disordered phase from the ordered one. The relation between the equilibrium measure of the shaken dynamics and the Gibbs measure for the Ising model is also investigated. Finally two different coding approaches are considered for the implementation of the dynamics: a multicore CPU approach, coded in Julia, and a GPU approach coded with CUDA.

1 Introduction

The Gibbs sampling of lattice spin models is a major task for statistical mechanics. The related numerical techniques are mainly based on Markov chain dynamics for single and cluster spin flip [5, 15, 17], and can be easily implemented by means of random mapping representation techniques [8].

The theory of parallel Markov chains as Probabilistic Cellular Automaton (PCA) dates back to 1989 [6]. These processes are characterized by a factorized transition matrix on the configuration space, and can be simulated numerically updating all the spins using the same random map [2]. A class of PCAs with the transition probabilities given in terms of a pair Hamiltonian with the spins simultaneously updated at each time step has been analyzed in [3, 10] and the corresponding PCA has been used to study the Ising model on planar graphs. We explore the computational possibilities of this model to generalize the random sampling algorithms for Ising spin systems on a set of geometries for the two-dimensional lattices.

Formally a PCA is a Markov Chain $({X_n})_{n\in {\mathbb {N}}}$ whose transition probabilities are such that given two generic configurations $\tau =(\tau _1,\ldots ,\tau _k)$ and $\sigma =(\sigma _1,\ldots ,\sigma _k)$

$$\begin{aligned} P\{X_n = \tau | X_{n-1} = \sigma \} = \prod _{i=1}^k P\{(X_n)_{i} = \tau _{i} | X_{n-1} = \sigma \} \end{aligned}$$

(1)

and the components of the “configuration” can be independently updated at each time n. The evolution of a Markov Chain of this type is well suited to be simulated on parallel processors and GPUs.

In this framework, a shaken dynamics for the PCA parameterized by J and q has recently been introduced [2], and a model has been proposed where the configurational variables are split into two groups $\tau =(\tau _1,\ldots ,\tau _k)$ and $\sigma =(\sigma _1,\ldots ,\sigma _k)$, where $\tau _i,\sigma _i\in \{-1,1\}$ for each i, and are arranged on a bipartite graph [1]. Different interactions among the $\tau$ and $\sigma$ variables give rise to the possibility of interpolation among different lattice geometries. The equilibrium measure of the shaken dynamics has been extensively investigated and the analytic form of the critical curve in the (q, J) parameters space has been found [1].

The PCA dynamics we consider is a parallel and irreversible version of the heath bath dynamics and is obtained by concatenating two different update rules [2] (shaken dynamics). By the Hamiltonian defined in [1], which depends on the (J, q) parameters, we identify numerically two regions of the space (J, q) characterized by different behaviors of the dynamics.

The elementary step of the shaken dynamics is naturally defined on the a finite subset $\Lambda$ of the square lattice ${\mathbb {Z}}^2$ and consists of a sequence of two inhomogeneous half steps. However, in both [1, 2] it has been pointed out that the shaken dynamics can be seen as an alternate dynamics on a subset of the honeycomb lattice. The proposed dynamics, although not faster than ad hoc dynamics for specific models, simulates a whole class of statistical mechanics models spanning from the one-dimensional Ising model to the square lattice and hexagonal one across all the intermediate models depending on the values of J and q. Some of the values of J and q are particularly interesting because they allow to use the shaken dynamics to simulate

• the Ising model on the isotropic hexagonal lattice for $J=q$

• (an approximation to) the Ising model on the square lattice for $q>> 1$

• the Ising model on a collection of weakly interacting unidimensional systems for small values of q.

The numerical investigation we put forward is aimed at:

• illustrating a simple heuristic method to numerically determine the critical curve

• evaluating the mixing time of the chain as a function of J and q

• studying the spin-spin correlations as a function of J and q.

Further, for $J = q$ we compare the mixing time of the shaken dynamics with that of a single spin flip dynamics for the Ising model on the hexagonal lattice and, for $q \gg 1$ we also compare the mixing time of the shaken dynamics with that of a single spin flip and an alternate parallel dynamics for the Ising model on the square lattice and evaluate the distance between the equilibrium measure of the shaken dynamics from the Gibbs measure for the Ising model on the square lattice.

The paper is organized as follows. In the next section we define the lattice spin model and the PCA dynamics. In Sect. 3 we present the numerical findings concerning the model under investigation. Finally, in Sect. 4, we provide some details concerning the two implementations of the dynamics taken into account: the multicore CPU one and the GPU one.

2 The Model

Consider the Ising Hamiltonian on a graph G(V, E)

$$\begin{aligned} H_G(\sigma )=-\sum _{(x,y)\in E}J_{xy}\sigma _x\sigma _y \end{aligned}$$

(2)

where $\sigma _x\in \{-1,1\}$ for all the $x\in V$ and $J_{xy}\in \mathbb {R}^+$.

We assume that $V=\Lambda _1\cup \Lambda _2$, where $\Lambda _1$ and $\Lambda _2$ are finite squared subsets of the square lattice with $L^2$ sites and periodic boundary conditions

$$\begin{aligned} \Lambda =\Lambda _1=\Lambda _2=(\mathbb {Z}/L\mathbb {Z})^2 \end{aligned}$$

(3)

and all the edges in E have one endpoint in $\Lambda _1$ and the other in $\Lambda _2$. The $\sigma$ and $\tau$ variables denote the Ising configuration on the vertices of $\Lambda _1$ and $\Lambda _2$. Each $\sigma _u$, with $u\in \Lambda _1$, can be put in a one-to-one correspondence with $\tau _u$ with the same index $u\in \Lambda _2$.

Let $x=(i,j)$ be a vector of coordinate on the torus $({\mathbb Z}/L{{\mathbb {Z}}})^2$. Then

$$\begin{aligned} x^\uparrow =(i,j+1),\quad x^\rightarrow =(i+1,j),\quad x^\downarrow =(i,j-1),\quad x^\leftarrow =(i-1,j) \end{aligned}$$

(4)

are the coordinates of the four points at unit distance from x. Set $J_{xy}=J$ for all $(x,u)\in E$ with $x\ne u$ and $J_{xy}=q$ if $x=y$.

With this notation we write the Hamiltonian studied in [1, 2]

$$\begin{aligned} \begin{aligned} H(\sigma ,\tau )&= - \sum _{x \in \Lambda } \left[ J \sigma _x(\tau _{x^{\uparrow }} + \tau _{x^{\rightarrow }}) +q\sigma _x \tau _x\right] \\&= - \sum _{x \in \Lambda }\left[ J \tau _x(\sigma _{x^{\downarrow }} + \sigma _{x^{\leftarrow }}) +q\tau _x\sigma _x \right] \end{aligned} \end{aligned}$$

(5)

on the pairs of Ising configurations $\sigma$ on $\Lambda _1$ and $\tau$ on $\Lambda _2$. The interactions of this Hamiltonian can be visualized on the induced bipartite graph represented in Figs. 1 and 3. The parameter q is also referred to as the self interaction parameter.

Fig. 1

The lattices $\Lambda _1$, $\Lambda _2$ with the q (red) and J (black) interactions (Color figure online)

As pointed out in [1] a careful look to the Hamiltonian (5) and to the graph of Fig. 1 shows that the bipartite graph is isomorphic to the hexagonal lattice with edges J and q on whose vertices are arranged the variables $\sigma$ and $\tau$ as shown in Fig. 2.

Fig. 2

The hexagonal graph

Fig. 3

A representation of the hexagonal graph that highlights the relation with the two square lattices $\Lambda _{1}$ and $\Lambda _{2}$

The Gibbs measure at temperature $1/\beta$ for the Hamiltonian (5) is

$$\begin{aligned} \pi _2(\sigma ,\tau )={e^{-\beta H(\sigma ,\tau )} \over \sum _{(\sigma ,\tau )\in {{\mathcal {X}}}\times {{\mathcal {X}}}}e^{-\beta H(\sigma ,\tau )}} \end{aligned}$$

(6)

where ${{\mathcal {X}}}\times {\mathcal X}=\{-1,1\}^{|\Lambda |}\times \{-1,1\}^{|\Lambda |}$ is the configuration space of the variable $(\sigma ,\tau )$. The critical value of $\beta _c$ separates the ordered phase where all the spin have the same probability to take the values $+1$ or $-1$ from the ordered phase where the measure is polarized [4].

It has been proven in [1] that the function (plotted in Fig. 4) separating the ordered phase from the disordered one is

$$\begin{aligned} J_c(q)=\tanh ^{-1}\big (\sqrt{\tanh ^2q+1}-\tanh q\big ) \end{aligned}$$

(7)

Fig. 4

The critical curve $J_c(q)$

We observe that

$$\begin{aligned} \lim _{q\rightarrow \infty }J_c(q)=\tanh ^{-1}(\sqrt{2}-1)=0.4406867 \end{aligned}$$

(8)

is the critical value of $\beta$ for the Ising model on the square lattice, while on the point $J_c(q)=q$ the equation (7) gives the critical value for the Ising model on the hexagonal lattice $J=q=0.6585$. If $q\rightarrow 0$ the value of J satisfying (7) tends to $\infty$. This is consistent with the fact that the one-dimensional Ising model does not exhibit a phase transition. The physical interpretation of this limit can be given in the framework of the Random Cluster model, [7] which assigns the weight $1 - e^{-2q}$ to the edge q. With this interpretation, when $q \rightarrow 0$ (see Fig. 3) the honeycomb lattice becomes a collection of non-connected one-dimensional lattices.

Following [2] we let the system evolve as a Markov chain where the spins in $\Lambda _1$ and in $\Lambda _2$ are alternatively updated with a probability proportional to the exponential of the Hamiltonian of the target configuration in ${{\mathcal {X}}}\times {{\mathcal {X}}}$.

More precisely, using the notation

$$\begin{aligned} \begin{aligned} \overrightarrow{h_{x}}({\sigma })&= J({\sigma }_{{x}^{\uparrow }} + {\sigma }_{{x}^{\rightarrow }}) + q {\sigma }_{x}\\ \overleftarrow{h_{x}}({\sigma })&= J({\sigma }_{{x}^{\downarrow }} + {\sigma }_{{x}^{\leftarrow }}) + q {\sigma }_{x} \end{aligned} \end{aligned}$$

(9)

we consider a Markov chain on ${{\mathcal {X}}}\times {{\mathcal {X}}}$ with transition probabilities given by

$$\begin{aligned} P((\sigma , \tau ),(\sigma , \tau ^{\prime })) = \frac{e^{-H(\sigma ,\tau ^{\prime })}}{Z_{\sigma }} = \frac{e^{-\sum _{u\in \Lambda }\overleftarrow{h_{u}}({\sigma }) {\tau ^{\prime }}_{u}}}{Z_{\sigma }} = \prod _{u\in \Lambda }\frac{e^{\overleftarrow{h_{u}}({\sigma }) {\tau ^{\prime }}_{u}}}{2 \cosh {\overleftarrow{h_{u}}({\sigma })}} \end{aligned}$$

(10)

at odd times and

$$\begin{aligned} P((\sigma , \tau ),(\sigma ^{\prime }, \tau )) = \frac{e^{-H(\sigma ^{\prime },\tau )}}{Z_{\tau }} = \frac{e^{-\sum _{u\in \Lambda }\overrightarrow{h_{u}}({\tau }) {\sigma ^{\prime }}_{u}}}{Z_{\tau }} = \prod _{u\in \Lambda }\frac{e^{\overrightarrow{h_{u}}({\tau }) {\sigma ^{\prime }}_{u}}}{2 \cosh {\overrightarrow{h_{u}}({\tau })}} \end{aligned}$$

(11)

at even times where $Z_\sigma =\sum _{\eta \in {{\mathcal {X}}}} e^{-H(\sigma ,\eta )}$ and $Z_\tau =\sum _{\eta \in {{\mathcal {X}}}} e^{-H(\eta ,\tau )}$.

The factorization in Eqs. (10) and (11) and the mutual dependence of the variables $\sigma$ and $\tau$ allow for a straightforward parallel numerical implementation. In particular, to simulate the evolution of the chain it is possible to sample the value $\zeta \in \{-1, 1\}$ of the spin at site u with probability $P(\tau ^{\prime }_{u} = \zeta |\sigma ) = \frac{e^{\zeta \overleftarrow{h_{u}}({\sigma }) }}{2\cosh {\overleftarrow{h_{u}}({\sigma })}}$ at odd times and $P(\sigma ^{\prime }_{u} = \zeta |\tau ) = \frac{e^{\zeta \overrightarrow{h_{u}}({\tau }) }}{2\cosh {\overrightarrow{h_{u}}({\tau })}}$ at even times independently for all $u \in \Lambda$.

In this framework, the shaken dynamics introduced in [2] is obtained by looking at the evolution of the spin configuration in $\Lambda _1$. In other words, the shaken dynamics is the marginal of the alternate dynamics defined by Eqs. (10) and (11) and the shaken transition probabilties are

$$\begin{aligned} P^{\square }(\sigma , \sigma ^{\prime }) = \sum _{\tau } \frac{e^{-H(\sigma ,\tau )}}{Z_{\tau }} \frac{e^{-H(\sigma ^{\prime },\tau )}}{Z_{\tau }} \end{aligned}$$

In [2] it has been also shown that the equilibrium measure for this dynamics is

$$\begin{aligned} \pi _{s}(\sigma )=\frac{{Z}_\sigma }{Z} \end{aligned}$$

with $Z = \sum _{\sigma } Z_{\sigma }$.

In the remainder of the paper, we use the wording shaken dynamics when we are interested in the evolution on the sub–lattice $\Lambda _1$ whereas we call the dynamics on the hexagonal lattice subject to the transition probabilities (10) and (11) the alternate parallel dynamics (on the hexagonal lattice).

3 Simulation Results

3.1 Numerical estimation of the Critical Curve

As stated before, the critical curve (7) is the function that separates the ordered and the disordered phases. Above this line the values of the spins tend to be highly correlated whereas on the opposite side the value assumed by each spin is weakly dependent on the values of other spins. To determine whether the system is in the ordered or disordered phase we compute, over a large number of iterations, the average and the variance of the magnetization m on one of the two layer $\Lambda _i$, defined as

$$\begin{aligned} m = \frac{1}{|\Lambda _{i}|} \sum _{x \in \Lambda _i}\sigma _x \end{aligned}$$

(12)

We know by the Theorem 2.1 in [1] $\pi _{s}(m) = \pi _{2}(m)$, that the average magnetization (in $\Lambda _{1}$) of the shaken dynamics is the same as the average magnetization of the parallel alternate dynamics (on the hexagonal lattice $\Lambda _{1} \cup \Lambda _{2}$).

To estimate the critical curve, we use a method that is based on the following heuristic argument. Once the dynamics is started from the configuration consisting of “all minus”, it will reach in a relative short time a minimizer of the free energy. During its evolution, the dynamics overcomes the free energy barriers of magnitude $\Delta H$ with a probability that is exponentially small in $\Delta H$. As long the dynamics remains in the high temperature regime (disordered phase), the minimizers of the free energy have zero magnetization. As time evolves, the dynamics tends to return very rapidly towards there states with zero magnetization. Therefore, at high temperature, we can expect the average over time value of the magnetization to be around zero and its variance to stay rather small. In the low temperature regime (ordered phase), the free energy minimizers have either an average magnetization that is “close” to $-1$ or “close” to $+1$ (the colder the system the closer the magnetization of these minimizer is to $\pm 1$). These minimizers are separated by very high free energy barriers. Therefore if the dynamics starts from “all minus” we can expect it to reach a “close to $-1$” minimizer very rapidly and, afterwards, to return quickly to such minimizers. Being the probability that, at low temperature, the dynamics starting from “all minus” reaches a “close to $+1$” mininimizer, it is reasonable to expect that the time average of the magnetization stays close to $-1$ and that its variance stays small. Finally, when the system is around the critical temperature, it has both “close to $-1$” and “close to $+1$” minimizers. However, the energy barriers separating these minimizers are not very tall and, if the dynamics is let to evolve for a sufficiently long time, it is likely to move from areas of the state space “close to $-1$” minimizers to areas of the state space “close to $+1$” minimizers. If this happens, the time average of the magnetization is expected be very close to zero whereas its variance is likely to be big.

For the numerical simulation we take the lattice $\Lambda$ to be a $200 \times 200$ torus and compute the evolution of the shaken dynamics starting from the configuration $\sigma _0 = \{-1, -1, \ldots -1\}$ for $(J, q) \in \{(0, 2) \times (0,2)\}$ on a $80 \times 80$ grid. For each pair (J, q) we compute the average and the variance of the magnetization over 300,000 time steps.

Figure 5 shows the average and the variance of the magnetization as a function of J for $q = 0.6585$. It is evident that the average magnetization has a sharp transition around the point $J = 0.6585$ which is the critical value of J for the Ising model on the honeycomb lattice. Around the same point the variance of the magnetization has a spike whereas its value is negligible for values of J far from the critical point.

The results obtained on the whole grid (q, J) are summarized in Fig. 6.

Fig. 5

Average (a) and variance (b) of the magnetization as a function of J for $q=0.6585$

Fig. 6

Average (a) and variance (b) of the magnetization on the whole (q, J) grid

It is known that, at equilibrium, the average value of the magnetization fluctuates heavily only close to the critical value of the interactions (see [14] for a reference). The Fig. 7 shows that the variance of the magnetization is significantly different from zero only for points of the (q, J) plane in the vicinity of the pins on the curve (7). This show that, even for a small lattice, the magnetization fluctuates only close to the critical line and for the whole class of Ising models that can be described tuning the values of J and q.

Fig. 7

The bars are centered at those points in the (q, J) plane for which the variance of the magnetization is sufficiently large ($\ge 0.03$). The length of the bars is proportional to the variance of the magnetization

3.2 Coalescence Times and Perfect Sampling

To assess whether the number of steps for which a Markov chain is run is large enough for its distribution to be close to the equilibrium distribution, the mixing time should be considered. For a Markov chain $({X_n})_{n\in {\mathbb {N}}}$ with state space ${\mathcal {X}}$ and stationary distribution $\pi$, the mixing time is defined as

$$\begin{aligned} T_{\text {mix}}= T_{\text {mix}}(\varepsilon ) = \min \{n > 0 : \Vert \mu _\sigma ^n - \pi \Vert _{\text {TV}} < \varepsilon \, \forall \, \sigma \in {\mathcal {X}}\} \end{aligned}$$

where $\mu _\sigma ^n$ is the distribution of $X_n$ conditioned on $X_0 = \sigma$, $\Vert \mu - \nu \Vert _{\text {TV}}$ denotes the total variation distance between the probability measures $\mu$ and $\nu$ and $\varepsilon$ is some “small” number (for instance $e^{-1}$). For a reference on mixing times see, for instance [11]. Determining useful bounds for the mixing time of a Markov chain is, in general, a quite challenging task. However, indication on the mixing time of a Markov chain can be gathered looking at the coalescence times (see [8] for a reference).

Consider two Markov chains $({X_n})_{n\in {\mathbb {N}}}$ and $({Y_n})_{n\in {\mathbb {N}}}$ living on the same state space ${\mathcal {X}}$ and consider the coupling $({Z_n})_{n\in {\mathbb {N}}} = (X_n, Y_n)$ obtained by letting $X_n$ and $Y_n$ to evolve with the same update function and the same sequence of random numbers. More formally, set $X_{n+1} = f(X_{n}, W_{n})$ and $Y_{n+1} = f(Y_{n}, W_{n})$ with $(W_{n})_{n \in {\mathbb {N}}}$ a sequence of i.i.d. random variables. Note that the update function f and the sequence of random variables $(W_{n})_{n \in {\mathbb {N}}}$ is the same for both processes (for an introduction on the coupling method see [12]). Further assume that the update function is chosen is a way such that $P_Z(X_n = Y_n) \rightarrow 1$ as $n \rightarrow \infty$.

We define the coalescence time T between $X_n$ and $Y_n$ as $T = \min \{n\in {\mathbb {N}}: X_n = Y_n\}$. Note that, since $X_n$ and $Y_n$ evolve with the same update function and the same sequence of random numbers $X_n = Y_n$ for all $n > T$. This definition extends naturally to a collection of K chains $X_n^k$ with $k \in 1\ldots K$.

The total variation distance between the distribution at time t between two coupled Markov chains can be bounded using the following lemma.

Lemma 1

Let $X^{t}$ and $Y^{t}$ be two Markov chains with state space S with initial distributions $\mu ^{0}$ and $\nu ^{0}$, respectively, and $Z^{t}=\left( X^{t}, Y^{t}\right)$ be a coupling on them. Then

$$\begin{aligned} \left\| \mu ^{t}-\nu ^{t}\right\| _{T V} \le P\left( X^{t} \ne Y^{t}\right) \end{aligned}$$

(13)

The proof is straightforward and we presented here for completeness.

Proof

Given any $A \subset S$, we have

$$\begin{aligned} \begin{aligned} \mu ^{t}(A)-\nu ^{t}(A)&=P\left( X^{t} \in A\right) -P\left( Y^{t} \in A\right) \\&=P\left( X^{t} \in A, Y^{t} \notin A\right) -P\left( Y^{t} \in A, X^{t} \notin A\right) \\&\le P\left( X^{t} \in A, Y^{t} \notin A\right) \\&\le P\left( X^{t} \ne Y^{t}\right) \end{aligned} \end{aligned}$$

(14)

The bound for $\nu ^{t}(A)-\mu ^{t}(A)$ is determined likewise. $\square$

Thanks to the previous lemma the mixing time of the chain $({X_n})_{n\in {\mathbb {N}}}$ can be estimated by the coalescence time of the chains $({X_n^k})_{n\in {\mathbb {N}}}$ for $k = 1\ldots |{\mathcal {X}}|$, all defined on the state space ${\mathcal {X}}$, where chain $X_n^k$ has initial distribution concentrated on state k.

To effectively determine the coalescence time of the shaken dynamics, however, it is not necessary to run $2^{|\Lambda |}$ copies of the Markov chains, but it is possible to use the so called sandwiching technique since the shaken dynamics preserves the partial ordering between configurations¹ . In other words, it can be directly checked that if $X_0^k \le X_0^l$ than $X_n^k \le X_n^l$ for all $n > 0$ (see, again, [8], for a reference). To determine the coalescence time it is therefore sufficient to look at the coalescence times of two chains starting, respectively, from $\sigma ^{\text {top}} = \{1, 1, \ldots , 1\}$ (the largest possible configuration) and $\sigma ^{\text {bot}} = \{-1, -1, \ldots , -1\}$ (the smallest possible one).

We observe that while leveraging on the coupling between Markov chains, it is possible to perform an unbiased sampling from the equilibrium distribution of a Markov chain using the Propp–Wilson algorithm [13], which requires two copies of the Markov chain to be run with the same update function and the same sequence of random numbers.

In the original version of the Propp–Wilson algorithm the chains must be run from the past. The starting time must be chosen so that the chains have coalesced by time zero.

Note that, in general, running the chains into the future until coalescence takes place leads to a biased sampling. However, it is possible to modify the algorithm to run it with read once randomness (see [16] and see also [8, Chapter 12] for a clear reference). The details of this modification of the algorithm are quite delicate. However, a key ingredient of this version of the algorithm requires to run pairs of coupled processes $(X_n, Y_n)$ and $(X'_n, Y'_n)$ until they both have coalesced. The number of times this procedure is required is distributed as a geometric random variable.

To study the coalescence times of the shaken dynamics, the simulations were run taking $\Lambda$ to be a $32 \times 32$ square lattice. This means that the induced hexagonal lattice $\Lambda _1 \cup \Lambda _2$ has $32 \times 32 \times 2$ points.

We computed the average coalescence time for values of J and q close to the critical line $J_c(q)$. The results obtained are summarized in Fig. 8.

Fig. 8

Logarithm of the average coalescence time for values of J and q close to the critical curve

For $J=q$, the shaken dynamics is the marginal of the alternate dynamics on the isotropic hexagonal lattice. More properly, pairs of configurations $(\sigma , \tau )$ with $\tau$ the configuration obtained from $\sigma$ by performing the first half step of the shaken dynamics can be regarded as spin configurations on the honeycomb lattice. The equilibrium distribution of these pairs is the Gibbs measure of the Ising model on the isotropic hexagonal lattice (see Theorem 2.1 in [2]). Therefore it makes sense to compare the mixing time of the shaken dynamics with the mixing time of a single spin flip dynamics defined on the hexagonal lattice and whose stationary distribution is the Gibbs measure. As a reference we take the heat bath dynamics defined as follows:

$$\begin{aligned} P(\sigma , \sigma ^{\prime }) = {\left\{ \begin{array}{ll} \frac{1}{|\Lambda |} \frac{e^{h_{x}(\sigma )\sigma _i}}{2 \cosh (e^{h_{x}(\sigma )})}&{} \text {if } \sigma ^{\prime }= \sigma ^{x} \\ 1 -\sum _{x\in \Lambda } P(\sigma , \sigma ^{\prime }) &{} \text {if } \sigma = \sigma ^{\prime }\\ 0 &{} \text {otherwise} \end{array}\right. } \end{aligned}$$

where $\sigma ^{x}$ is the configuration obtained from $\sigma$ by flipping the spin at site u and $h_{u}(\sigma ) = \sum _{y\sim x} J \sigma _y$. Also the heat bath dynamics preserves the partial ordering between configurations and, hence, also in this case it is sufficient to simulate the evolution of two chains one starting from all spins set to $+1$ and one starting from all spins set to $-1$.

Note that it is possible to argue that the parallel alternate dynamics studied here is a parallel version of the single spin flip heat bath described above.

The results obtained for several values of J (and, consequently, q) are shown in Fig. 9. Note that for the single spin flip dynamics the value shown in the chart is the number of steps divided by $2|\Lambda |$ so that, for both algorithms, we are comparing the total number of “attempted spin flips”.

The parallel alternate dynamics is faster than the single spin flip one even if the single spin flip one is “renormalized” with the volume of the box as described above.

Fig. 9

Sample average of the coalescence time for $J=q$ (hexagonal lattice)

In [2], Theorem 2.3 it has been shown that, for large values of q, the equilibrium distribution of the shaken dynamics approaches the Gibbs measure for the Ising model on the square lattice. More precisely it has been proven that, if

$$\begin{aligned} \lim _{|\Lambda | \rightarrow \infty } e^{-2q}|\Lambda | = 0, \end{aligned}$$

then, for J sufficiently large,

$$\begin{aligned} \lim _{|\Lambda | \rightarrow \infty } \Vert \pi _s - \pi _G\Vert _{\text {TV}} = 0, \end{aligned}$$

where $\pi _G$ is the Gibbs measure for the Ising model on the square lattice. Therefore it makes sense to evaluate numerically the goodness of this approximation as q increases. To this purpose we consider two observable: the magnetization m and the energy $H(\sigma )$. For both observable we compare their sample mean and sample standard deviation over samples drawn from the equilibrium distribution of the shaken dynamics with the sample mean and the sample standard deviation of two other reference dynamics having the Gibbs measure as stationary distribution. One of the two reference dynamics taken into account is, again, the heat bath dynamics. The other dynamics is a parallel version of the heat bath dynamics that updates, alternatively, the spins on the odd and the even sites of the lattice. The latter is the equivalent for the square lattice of the alternate parallel dynamics on the hexagonal lattice defined by equations 10 and 11. Theorem 2.2 in [2] states that the equilibrium measure of this dynamics is, indeed, the Gibbs measure on the square lattice. For all these dynamics, samples are drawn using the Propp-Wilson algorithm introduced above. Several values of J close to the critical value for the Ising model on the square lattice and the results obtained are summarized in Figs. 10, 11, 12 and 13.

The data suggests that, for $q \ge 2.5$ the approximation provided by the shaken dynamics is quite good.

Fig. 10

Sample average of the magnetization for several values of J

Fig. 11

Sample standard deviation of the magnetization for several values of J

Fig. 12

Sample average of the energy $H(\sigma )$ for several values of J

Fig. 13

Sample standard deviation of the energy $H(\sigma )$ for several values of J

We also estimated the time required to approach the equilibrium distribution by comparing the coalesce time of the shaken dynamics with those of the two other reference dynamics. Also in this case the number of steps required by the single spin flip dynamics is renormalized with the volume of the box $\Lambda$. The results obtained are summarized in Fig. 14. It is apparent that, though more flexible, the shaken dynamics becomes slower than “specialized” algorithms as the accuracy of the approximation increases.

Fig. 14

Sample average of the coalescence time (number of steps) for several values of J

Fig. 15

$J = 0.44$, $q = 3.0$

Parts (b) of Figs. 15, 16 and 17 show configurations drawn from the equilibrium distribution of the alternate parallel on the hexagonal whereas parts (a) show the corresponding sub-configurations on the sublattice $\Lambda _1$. These sub-configurations are, therefore, drawn from the equilibrium distribution of the shaken dynamics. In Fig. 15 it is possible to observe that the spins linked by a q-edge have almost always the same value. This is in good accordance with the fact that stationary measure of the shaken dynamics is close to the Gibbs measure for the Ising model on the square lattice. On the other side, Fig. 17 is consistent with the fact that for q very small the equilibrium measure of the shaken dynamics tends to that of a colection of weakly dependent unidimensional Ising models.

Fig. 16

$J = 0.6585$, $q = 0.6585$

Fig. 17

$J = 2.0$, $q = 0.03$

3.3 Correlations

Theorem 2.4 in [1] establishes that, if q is sufficiently small, $\pi (\sigma _{0,0}, \sigma _{\ell ,\ell }) < \pi (\sigma _{0,\ell }, \sigma _{\ell ,0})$ where $\pi$ is the equilibrium measure of the shaken dynamics and $\sigma$ is, therefore, a spin configuration living on $\Lambda _1$. In other words, the theorem states that the SW-NE correlations are weaker than the NW-SE ones if the self interaction is weak. We also expect that the SW-NE and the NW-SE correlations tend to be similar for large values of q, that is for those values of the pair (q, J) for which the equilibrium distribution of the shaken dynamics approaches the Gibbs measure of the Ising model on the square lattice. The results of the analysis of the SW-NE and the NW-SE correlations as $\ell$ varies with $\Lambda$ a $32 \times 32$ square box are shown in Table 1.

Table 1 Spin-spin correlations

All pairs (q, J) taken into account correspond to points of the q, J plane close to the critical curve $J_c(q)$. It is possible to observe that, as q decreases, the SW-NE correlations become, indeed, smaller than the NW-SE ones, whereas, for q large the two are quite similar. Further, if the pair (q, J) is below the critical curve the correlations decay quite rapidly. On the other hand, if (q, J) is above $J_c$ the correlations are significant also for larger values of $\ell$.

4 Implementation Details

To approximate numerically the critical curve $J_c(q)$, we take samples for different values of J and q. The code used for the simulation are written in Julia [9] and the simulations are performed through 80 thread processors running, in parallel, the simulation on 80 couples of values (q, J) in the range of $(q,J) \in (0,2)\times (0,2)$. The Hamiltonian is defined on a square $200\times 200$ lattice. Statistics are collected over 300,000 iterations. Figure 6 shows that the chosen simulation parameter is good enough to approximate the critical curve.

The elementary step of the shaken dynamics described in the previous section has been simulated by the Algorithm 3. A spin configuration is updated via a sequence of two similar half steps. The computation of the vector of local fields h that drives the transition probabilities of each spin is alternatively carried out using the functions collectUR and collectDL which determine the up-right and down-left contribution as in Eq. (9).

The algorithm 3 is the complete update in two steps of the shaken dynamics, which is more general than the one used in [10].

The choice of collecting the statistics over 300,000 time steps (after a warm up time of 300,000 additional stime stesps) turned out to be good enough, and the results show, unmistakably,the separation of the two phases (ordered and disordered).

We implemented the algorithm 3 in two parallel ways. A CUDA² implementation of a parallel heat bath for large dimension Lattice spin, and a Julia [9] implementation on a single CPU to be used on a multiprocessor systems (trivial parallel on a multi data input). Both have been optimized to handle our problem, and used to simulate the shaken dynamics of the PCA, a quasi-similar behavior was observed during our experiment.

4.1 Parallel Single-GPU Code

The general heat bath procedure has been implemented on SIMD (Single Instruction Multiple Data) system. To optimize the code exploiting the CUDA memory architecture we implemented three kernels for the functions collectUR, collectDL and for updating the configuration. We used the default random generator from curand library.

The collect function computes the transition probabilities in the given direction. Each thread handles one spin on the lattice field. The principal use of the global memory is the four square spin lattice, the two configuration sigma ($\sigma$) and tau ($\tau$), the fields which handles the Hamiltonian computation and the random-unit contains random uniform variables.

In our implementation all the operation are performed on register cache memory. We did not use shared memory for the random-unit. The code was written to run on the Nvidia-GPU Tesla P100 using by 16GB video memory, using 4 matrices of dimension $L\times L$, two for the lattice spin field (single byte), and two for the collected fields and for the random uniform (four byte). All the matrices are allocated on the global memory.

For the management of the memory, before allocating the memory of the 4 matrices, the code used approximately 303 MB, leaving 15973.250000 MB. We used $2*4*L*L + 2*L*L$ bytes, but we can not go beyond $10^5$ for this GPU.

The purpose of the CUDA implementation is to work on large dimension which allows us to observe the statistical behavior of the shaken dynamic, also for real time simulation.

Figure 18 shows a state of configuration captured at 60-th iteration on a simulation of the shaken dynamics for PCA lattice square with dimension $512\times 512$, under the temperature $J = 0.99$ and the external field $q=0.5$.

Fig. 18

GPU sample: $L=512$, $J=0.99$, $q=0.5$, iteration = 60th

4.1.1 Benchmarking

To measure the performance of our GPU code it is not fair to compare it with the single-CPU implementation from Julia. We have implemented a serial version of the shaken dynamics in a lower level language, a captured sample for a square lattice spins of size $512\times 512$ with $J=0.99$ and $q=0.5$ is given in Fig. 19.

Fig. 19

CPU sample: $L=512$, $J=0.99$, $q=0.5$, iteration = 60th

For our simple measurement, we set the parameters $J=0.44$ and $q=0.66$ and to have more significant value we measure it in milliseconds. We compare the two implementation for different dimension L, hence the size of the square lattice which is $L^{2}$.

Fig. 20

Running-time in function of the size L

For this benchmark (Fig. 20), we used an Nvidia graphic card Tesla P-100 vs single core of the CPU Intel(R) Xeon(R) CPU E5-2698 v4 @ 2.20GHz. We measure the time for one update execution. As we have mentioned the GPU memory is limited so that we did the experiment under this condition for the size of the square lattice spin.

We observe that the GPU is much faster than the CPU with a factor of 500 as far as the lattice size grows. We can see that the CPU time looks linear while the case for GPU is when the size is more than $2048\times 2048$.

5 Conclusions

5.1 Summary

The present work is a numerical experiment on the 2D Ising model. In particular the tasks are mainly focused on the shaken dynamics for approximating numerically the critical curve which relate the parameters J and q, it separates the two phases on the region of parameters (J, q). Also we are able to evaluate numerically the fact that the equilibrium distribution of this dynamics are close to Gibbs measure on the region when q is large [2]. We provide two different parallel implementation perspective of the algorithm in which we give a benchmark to identify a speed-up that we can gain on using GPU.

As a MCMC algorithm we also give a comparison on the convergence to the equilibrium state of the algorithm by means of coalescence time. In this purpose we compare the coalescence time between the alternate and shaken dynamics on the critical line (red line in Fig. 4), this is followed by further discussion on the numerical aspect of the algorithm.

5.2 Work in Progress

Most of the implementations we have used in this project are in Julia [9]. As pointed out above, the paper is a numerical supports for the two papers [2] and [1]. The code is quite complete for an academic use on the simulation of 2D Ising model, especially it contains the class of all the dynamics we have studied in this project and their methods. Our attempt is to provide a Julia library that can be used as framework to study the dynamics of the planar Ising model.