From network reliability to the Ising model: A parallel scheme for estimating the joint density of states

Network reliability is the probability that a dynamical system composed of discrete elements interacting on a network will be found in a configuration that satisfies a particular property. We introduce a reliability property, Ising feasibility, for which the network reliability is the Ising model's partition function. As shown by Moore and Shannon, the network reliability can be separated into two factors: structural, solely determined by the network topology, and dynamical, determined by the underlying dynamics. In this case, the structural factor is known as the joint density of states. Using methods developed to approximate the structural factor for other reliability properties, we simulate the joint density of states, yielding an approximation for the partition function. Based on a detailed examination of why naïve Monte Carlo sampling gives a poor approximation, we introduce a parallel scheme for estimating the joint density of states using a Markov-chain Monte Carlo method with a spin-exchange random walk. This parallel scheme makes simulating the Ising model in the presence of an external field practical on small computer clusters for networks with arbitrary topology with $\sim {10}^6$ energy levels and more than ${10}^{308}$ microstates.

Figure 1

(a) An Ising-feasible configuration with three spin-up vertices (red dots) and eight discordant spin pairs or edges (solid line segments). (b) A configuration that is not Ising feasible because of the inconsistent edges (red line segments). Independently choosing edges and spin-up vertices will rarely produce an Ising-feasible configuration, but any set of randomly chosen spin-up vertices uniquely determines a set of edges.

Figure 2

(a) Conditional state distribution of $p\left(k|v\right)$ sampled by a naïve Monte Carlo simulation on a $16\times 16$ square lattice. The conditional probability is normalized separately for each value of $v$. The color reflects the value of $p\left(k|v\right)$ in logarithmic scale. As $p\left(k|v\right)$ is symmetric along $v=N/2$, the simulation is only done for $v\le N/2$. (b) The peaks of $p\left(k|v\right)$ for various lattices have the same functional form: $y\left(x\right)=2x\left(1-x\right)$, where $x\equiv v/N$ and $y\equiv k/M$.

Figure 3

The exact joint density of states $g\left(k,v\right)$ computed via exhaustive enumeration on a $5\times 5$ square lattice. Because a naïve Monte Carlo only samples points near the peak of each curve, and the tails of many curves are as important as the peaks of others, it severely underestimates the univariate density of states $g\left(k\right)\equiv {\sum }_{v}g\left(k,v\right)$ at $k=16$.

Figure 4

The running time for estimating the joint density of states $g\left(k,v\right)$ on 2D lattices of different sizes, from $N=8\times 8$ to $N=24\times 24$. The sequential Wang-Landau algorithm (blue) needs to cover $O\left(N^2\right)$ energy states, and become impractical for large systems. The spin-exchange WL algorithm (red) divides the energy states into $N/2$ energy slices, and the running time is bounded by the energy slice with the most number of states at $v=N/2$. This energy slice only contains $O\left(N\right)$ states. By dividing the energy slice into six equal-sized, 75% overlapping energy windows, the running time is reduced even further (yellow).

Figure 5

(a) The joint density of states $g\left(k,v\right)$ of a $32\times 32$ square lattice. (b) The univariate density of states $g\left(k\right)\equiv {\sum }_{v}g\left(k,v\right)$ estimated using a spin-exchange MCMC algorithm, compared with the exact analytical result. (c) The internal energy $U=-\partial lnZ/\partial \beta$, and (d) the heat capacity $C=\partial U/\partial T$ with coupling strength $J=0.5$.

Figure 6

(a) A Cayley tree with degree $d=3$ and number of shells $r=6$. (b) The joint density of states $g\left(k,v\right)$ of the Cayley tree with $d=3$ and $r=8(N=765)$. As $d=3$, there are many inaccessible states. (c) The heat capacity of the Ising model on the Cayley tree. The red curve is from the exact solution at $B=0$.

Figure 7

Compare the heat capacity with the known results at $B=0$ (without external field). (a) On a $32\times 322\mathrm{D}$ square lattice. The exact solution is calculated from the exact density of states $g\left(k\right)$ from Beale [30]. (b) On a $d=3$ and $r=6$ Cayley tree. The exact solution is computed as $C/k=-{\beta }^2{d}^{2}lnZ/d{\beta }^2$, where $lnZ$ is given by Ref. [32].