Monte-Carlo

Monte-Carlo sampling is where the Newman-Ziff algorithm really displays its power. A random order of occupation is generated, and the onset of cluster wrapping is detected as more sites are added. Over S runs we record the number of times that cluster wrapping occurred for each number of occupied sites (1 -- N). This micro-canonical ensemble is then used to estimate the critical site occupation probability p_c as described in the method section.

The percolation simulations were run on lattices with up to 256x256 sites. The CPU times taken for each size are given in Table 2. The number of monte-carlo samples for each is constant (1000). The time is of course directly proportional to the number of samples.

The data in Table 2 is plotted on the following figure. It shows that the CPU time scales close to linearly with the number of lattice sites, as expected from the algorithm. It seems that the O(NlogN) behaviour expected from omitting path compression is not too significant at this scale.

We can take advantage of the exact results to evaluate methods of error estimation. The term 'error' here refers to statistical error in sampling p_c for finite L. One common method for monte-carlo simulations is to do several subsamples and then take their mean as the result; an error estimate is given by the variance between subsamples. An even simpler, intuitive method is just to look at results with increasing numbers of samples, estimating the error magnitude according to which decimal places are changing. Both these methods were applied to monte-carlo estimation of p_c on the 6x6 lattice and are shown in Table 3 along with the actual error.

There were 4 subsamples making up each estimate, and the error was taken as the average deviation from the final result (mean). In fact, taking the mean of S subsamples reduces the error by 1/sqrt(S) -- in this case 1/2 -- so these estimates are conservative, as can be seen by comparing to the actual error. The intuitive visual estimate is generally similar in this case.

The actual error with increasing samples on the 6x6 lattice are shown in the following plot, confirming a general scaling behaviour of 1/sqrt(S).

Ziff & Newman (2002) calculate the standard deviation of the p distribution for various estimators on open systems. With the estimator used here it is not clear how to do so efficiently. It was attempted by calculating p_c on the basis of each individual estimate of R_L(p) (each run of the algorithm on N matrices, having a single onset of percolation). The standard deviation could then be calculated as usual. Unexpectedly, these were found to change very little with the number of samples, only the linear dimension L -- these are given in Table 4. It was also found that the mean of all individual estimates was a poor approximation of the estimate derived from the final micro-canonical ensemble, differing in the 2nd or 3rd decimal place (presumably because of rounding errors).

By looking at error estimates and running times, tradeoffs were made to get reasonable accuracy across a range of lattice sizes. The best of each of these are given in Table 5, and plotted on the following figure. The subsample variance method was used to estimate errors, and as above this is conservative: the given values can be halved confidently.

As mentioned in the method section, the finite size scaling with this estimator is supposed to go as p - p_c ~ L^-11/4 where p is the estimate on an LxL lattice, and p_c is the infinite system value. Plotting the above estimates transformed by this power reveals that they don't fall on a straight line. In fact, by looking at a log-log plot it can be seen that they don't fit nicely to any constant exponent: L = 6--16 seems to scale differently to L = 32--256 (rough ranges). Perhaps at these small sizes it hasn't reached its asymptotic behaviour.

Considering only the larger systems then, the Matlab polyfit function was used to fit a line through the points L = 32, 64, 128, 256. This is shown below:

Most estimators, especially on open systems, have finite size scaling according to p - p_c ~ L^-7/4. This transformation was tried here, and the same line fitting was done, except that it was also possible to include L=16 in a reasonably good fit, as shown on the following plot:

The equation of the first fitted curve is
p = 0.5927491 + 0.3685080 L^(-11/4)
and that of the second is
p = 0.5927450 + 0.0150330 L^(-7/4).

The intercept of a line fitted to finite-size estimates of p gives an estimate of the infinite system value p_c. Since the statistical errors on fitted points were of order 1e-5 (conservatively), and the lattices involved were of reasonable size, we might expect similar accuracy for our estimate of p_c. Newman & Ziff (2001) give the critical site occupation probability to an accuracy of 8 digits: 0.59274621. From this we see that the first estimate above (0.5927491) has an error of 2.9e-6 and the second (0.5927450) has an error of 1.2e-6.

An obvious way to parallelise this algorithm is to partition the number of samples between processors, each generating a subsample of the micro-canonical ensemble. Aggregating them is a simple sum of arrays. An advantage of this is that each processor can generate a result from its subsample before aggregating, and the variance of these is an error estimate as described above. In fact, this was how the above error estimates were produced.

The parallel efficiency is good since the process is so simple. For example, the extrapolated time to run 1 million samples of a 256x256 lattice serially was 37373.5 seconds (10:22:53.5). Running it over 8 CPUs took 4640.8 seconds (1:17:20.8) which is 99.33% efficient.