Exact algorithm for sampling the two-dimensional Ising spin glass

A sampling algorithm is presented that generates spin-glass configurations of the two-dimensional Edwards-Anderson Ising spin glass at finite temperature with probabilities proportional to their Boltzmann weights. Such an algorithm overcomes the slow dynamics of direct simulation and can be used to study long-range correlation functions and coarse-grained dynamics. The algorithm uses a correspondence between spin configurations on a regular lattice and dimer (edge) coverings of a related graph: Wilson’s algorithm [D. B. Wilson, Proceedings of the Eighth Symposium on Discrete Algorithms (SIAM, Philadelphia, 1997), p 258] for sampling dimer coverings on a planar lattice is adapted to generate samplings for the dimer problem corresponding to both planar and toroidal spin-glass samples. This algorithm is recursive: it computes probabilities for spins along a “separator” that divides the sample in half. Given the spins on the separator, sample configurations for the two separated halves are generated by further division and assignment. The algorithm is simplified by using Pfaffian elimination rather than Gaussian elimination for sampling dimer configurations. For $n$ spins and given floating point precision, the algorithm has an asymptotic run-time of $O\left(n^{3/2}\right)$; it is found that the required precision scales as inverse temperature and grows only slowly with system size. Sample applications and benchmarking results are presented for samples of size up to $n={128}^2$, with fixed and periodic boundary conditions.

Figure 1

Results of applying the sampling algorithm to an individual 2D Ising spin-glass sample, for temperatures $T=0.5,0.2,0.08,0.02$, for a single Gaussian spin-glass sample with fixed boundaries. The images show the variability of the spin assignments (top) and of the domain walls (bottom) over a range of temperatures, in a sample with $n={126}^2$ variable spins surrounded by a layer of fixed spins. At least 240 samples were generated at each temperature. The grayscale values indicate the probability of a given spin being fixed (upper row) or of neighboring spins being fixed relative to each other (lower row). For spin assignments, the darkest colors indicate that the spin is equally likely to be up or down, while light colors indicate that the spin occurs with a single alignment in nearly all sampled configurations. These alignments result from correlations with the fixed boundary spins. For the domain walls displayed in the lower row of images, the lines indicate the probability of relative domain walls between two configurations: the darkest lines indicate the bond dual to that domain wall has a 50% chance of opposite or equal relative orientations; where there is no line separating two spins (or only a very light one), the two spins have a very high probability of a single relative orientation, either aligned or opposite. Specifically, the bond satisfaction variance ${\mu }_{i,j}\left(1-{\mu }_{i,j}\right)$ is plotted along each dual edge, where ${\mu }_{i,j}$ is the frequency of the $J_{ij}{s}_i{s}_j$ being positive. Note that as $T$ decreases, the frequency of specific droplet excitations, outlined by domain walls, can either increase or decrease, reflecting the sensitivity of the configurations to temperature. This can be seen, for example, in two of the regions that are active at $T=0.02$, the approximately $20\times 20$ region in the far upper left and the approximately $30\times 60$ region at the center right: the spins in the former become more fixed as temperature decreases while the spins in the latter region become more variable (darker) when the temperature is decreased from $T=0.08$ to $T=0.02$.

Figure 2

(Color online) A depiction of the correspondence between domain-wall loops for an Ising spin system and dimer matchings on the decorated dual lattice $G$. (a) A spin system with fixed boundary conditions; an up arrow at location $i$ indicates $s_i=+1$ and a down arrow indicates $s_i=-1$. The dual lattice $G_0$ is indicated by the lines connecting the dual nodes. (b) A Fisher city replacement. Each dual lattice node is expanded to a Fisher city, a set of six nodes composed of two linked triangles, to generate the decorated lattice. For work on the square lattice, the bond strengths are set to be $w\left(e_{ij}\right)=0$ inside the city, and the bond strengths between the cities, indicated here by the notation $w_d$, $d=0,1,2,3$, are set according to Eq. (3). (c) An example dimer covering (i.e., perfect matching) $M$ on the decorated graph $G$. The thicker (also red) bonds with circular ends indicate edges in $M$. The domain walls, composed of dimers that connect distinct cities, are indicated by dashed lines. (d) When the cities are contracted out from $G$, the loops on $G_0$ remain. Given this choice of dimer covering $M$, the spins that are inside the domain walls are flipped to create the new sampled configuration.

Figure 3

The node indexing and edge orientations within a Fisher city (left) and the corresponding elements of the $6\times 6$ submatrix of the Kasteleyn matrix $K$ (right). The numbering of nodes is shown for the first city listed in the dual lattice; subsequent cities have multiples of 6 added to their indices. In the case of the square spin lattice (indicated by the outer square bonds on the original spin lattice), all nonzero $K$ elements, including $d$, are set to unit magnitude. The labeling of the $0\to 2$ and $1\to 2$ edges indicate how the strengths can be modified in the case of the triangular lattice: in this case, one can set $d=\text{exp}\left[-\beta w\left(e_d\right)\right]$ to account for the diagonal bond $e_d$ perpendicular to the $2\to 3$ edge. The Kasteleyn matrix has row $a$ and column $b$ indices $a,b=0,\dots ,5$.

Figure 4

(Color online) An example of the nested dissection and the Kasteleyn matrix $K$ for a $6\times 6$ spin lattice sample surrounded by an outer layer of fixed spins. (a) The set of $8\times 8$ Ising spins sit on the sites of the light gray lattice of bonds of strength $J_{ij}$, where the diagonal bonds are indicated for the case of a triangular lattice. The graph $G$ on which the dimer sampling is computed is shown by the darker lines and circular nodes. The gray bands indicate the nested dissection used for these nodes: the lighter gray region contains the dimer separator $\Delta \subset G$ and is bordered by the middle two rows of spins, the spin separator $D$. The darker and medium bands indicate, in order, the subsequent subdivisions of the sample. (b) A display of the nonzero elements of the Kasteleyn matrix $K$, for a left-to-right and top-down ordering of the Kasteleyn cities. The nonzero elements of the $294\times 294$ Kasteleyn matrix for this are shown as black dots. The edges internal to the Kasteleyn cities are closest to the diagonal: further nonzero elements represent connections between the cities. (c) The permuted matrix $K$, where the cities are indexed according to a nested dissection. The gray regions in $K$ include connections contained within each of the separators of the nested dissection, with the same shades as in (a), and between the separators and other nodes. The procedure of Pfaffian elimination can at most affect elements within the gray regions and also the values near the diagonal, for the nodes not contained in the gray regions in (a). Spin values are assigned to $D$ by examining the part of $K^{-1}$ indexed by the nodes of $\Delta$, i.e., the lower right square submatrix contained within the light gray region.

Figure 5

Example of cross operations used for Pfaffian elimination. (i) A skew-symmetric matrix $K$. (ii) The result of a cross operation $K\to L\left(\alpha ,i,j\right)K{L}^T\left(\alpha ,q,r\right)$ of the first type, with $q=0$, $r=2$, $\alpha =-a$, applied to $K$. This is found by adding $\alpha =-K\left(q,r\right)/K\left(q,q+1\right)$ times column $q+1$ to column $r$ and then $\alpha$ times row $q+1$ to row $r$ to eliminate the element at location $\left(q,r\right)$. (iii) The result of the next cross operation, with $q=0$, $r=3$, and $\alpha =-b$. (iv) The result of two subsequent operations of the second type, where $\alpha \left(q,r\right)=-K\left(q,r\right)/K\left(q,q-1\right)$, for $q=1$, $r=2$ and $q=1$, $r=3$. These latter types of operation are not needed to compute $\text{Pf}\left(K\right)$ but are needed for finding ${\left[K^{-1}\right]}_{\Delta }$. The Pfaffian of $K$ is the product of the superdiagonal elements in even rows: here, $\text{Pf}\left(K\right)=\left(1\right)\left(1-ac+b\right)$.

Figure 6

(Color online) An example of the dimer assignment procedure for a separator $\Delta$ that is three cities wide. Initially, no edges are matched (left part of $k=0$). The first choice, $k=0$, is between the two edges inside $\Delta$ that are incident upon the far left node. In the example shown, the lower bond (connecting node 0 to node 2; see Fig. 3) is chosen, as shown on the left of the $k=1$ section of the figure. At this stage, one has computed matrices $A_1$, $M_1$, and $V_1$. The comparison of the next two possible matchings, shown on the right part of the $k=1$ subfigure, compares the inclusion of the $\left(q,r\right)=\left(3,5\right)$ and (3,4) edges. In some cases, as in the first $k=3$ panel, a choice is forced and $A_k$, $M_k$, and $V_k$ need not be updated. The $k=4$ choices are forced as an even number of domain walls must cross the separator so that an even number of the top nodes and an even number of the bottom nodes are unmatched.

Figure 7

(Color online) Application of the result of the example dimer assignment from Fig. 6 to the spin assignment. (a) The initial spin configuration with fixed spins on the boundary. The portion of $G$ used to compute $K$, ${\left[K^{-1}\right]}_{\Delta }$ and the dimer assignment is indicated in gray. The middle row of Fisher cities composes $\Delta$, the dimer separator. (b) The sample dimer assignment (partial matching) for $\Delta$ from Fig. 6 superimposed on $G$. (c) Extra choices in the matching are forced by the matching internal to $\Delta$. These additional dimers cut across the bonds separating spins in $D$, the two spin rows parallel to $\Delta$. (d) In the last step, the modifiable spins are updated. The update is based upon the portion of domain walls forced by the partial matching in (c). Moving from left to right, for example, from the two fixed spins on the middle of the left side, a spin is reversed if an odd number of dimers extending from $\Delta$ are crossed.

Figure 8

(Color online) The plot shows the sample averages of the number of bits of precision $B_{avg}$ required to obtain the correct sampling, for system sizes $n=L^2=6^2$ through ${126}^2$ with fixed boundaries, as a function of inverse temperature $\beta$, for Gaussian disorder. The lines indicate linear fits of the form $B_{avg}=c\beta$. The number of bits needed for periodic boundary conditions (not shown) are very close to these same lines, at each system size. To find an accurate result with high confidence, one can use twice the average needed value: this was sufficient for all samples $\left(>{10}^4\right)$ that we examined.

Figure 9

(Color online) Run time, measured in seconds, to generate a single configuration, as a function of system size $L$, using a 2.4 GHz Intel Core 2 Duo processor (MacBook Pro). Double precision (64 bit floating point) data are indicated with triangles, while multiprecision results for $B=512$, 2048, and 8192 bits are indicated by squares, diamonds, and circles, respectively. Samples are generated for $L\le 128$ with fixed boundaries (closed symbols) and for $L\le 64$ with periodic boundaries (filled symbols). The sample-to-sample fluctuation of disorder realization is less than 0.1% of the run time so error bars are not shown. The solid line indicates the form of the expected dependence of run time on system size for a given fixed precision, that is, $\sim L^3$.

Figure 10

Relative domain walls found in an individual 2D Ising spin-glass sample with ${64}^2$ spins, periodic boundary conditions, and unit variance Gaussian disorder, for temperature $T=0.16$. The bond satisfaction probabilities were estimated by averaging over 660 samples. As in Fig. 1, the lines indicate the probability of relative domain walls between two configurations: the darkest lines indicate where the bond dual to that domain wall has a nearly 50% chance of opposite or equal relative orientations; where there is no line or a light line separating two spins, the two spins have a very high probability of a single relative orientation, either aligned or opposite.