What are good random number generators for Monte Carlo calculations?

Question

Many problems in Computational Physics need the use of random number generators. When studying magnetic materials using the Heisenberg/Ising Hamiltonians (related questions about them can be seen here, here and here), one of the most used tools is the Metropolis method¹.

The Metropolis method can be summarized in the chart below (the image is linked directly from here):

In order to accept the energy due to the exchange in a spin state, a random number is generated and then compared with the probability of the change. So, this number must be as random as possible and with very low periodicity.

I just remember (back to the late 90's) that one of the used random number generators was the Kirkpatrick-Stoll² (also called R250).

My question: Are there any other modern random number generators better than R250?

(By "good", I mean: low correlation, long period and fast to code and run.)

1. Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 21 1087–1092 (1953) (DOI: 10.1063/1.1699114)
2. Kirkpatrick, S., and E. Stoll, A Very Fast Shift-Register Sequence Random Number Generator, J. Comp. Phys. 40 517-526 (1981) (DOI: 10.1016/0021-9991(81)90227-8)

Answer 1

Depending on how many random numbers you need in a short amount of time, it might be worth to consider using a cryptographic PRNG. In particular AES-CTR.

Now of course you might say that "but AES-CTR is soo slow". Actually it isn't. If you only re-key every GB or so of generated data you get a speed of

• about 12 CPU cycles needed for 16 bytes of random data per core on all Intel CPUs since Sandy Bridge, meaning you get $4/3\cdot f\cdot c$ bytes for $f$ being the CPU frequency and $c$ being the core count.¹
• about 6 CPU cycles needed for 16 bytes of random data per core on all AMD Ryzen-based CPUs, meaning you get $8/3\cdot f\cdot c$ bytes for $f$ being the CPU frequency and $c$ being the core count.¹
• about 12 CPU cycles needed for 64 bytes of random data per core on all Intel CPU with AVX-512 and the new VAES extension (currently only Ice Lake), meaning you get $16/3\cdot f\cdot c$ bytes for $f$ being the AVX-512 frequency and $c$ being the core count. Alternatively, Zen3 CPUs have 256-bit VAES as well and achieve the same efficiency.
• about 30GB/s on a high-end Nvidia Maxwell and newer GPUs. Newer research indicates even higher performance at almost 40GB/s on high-end Maxwell (970) and almost 110GB/s on high-end Turing (2070S).

Additionally the quality of the random bytes you get is top-notch with even a determined adversary having an expected chance of less than $2^{-32}$ for telling $2^{36}$ AES-produced bytes apart from truly random bytes for a secret key. Then to maintain that quality you simply re-key every second or so (of full-speed output) without losing too much of your security bound.

---

Footnote 1: These numbers are justified with the AESENC instruction throughput which is 1 for all supported Intel architectures and 0.5 for Ryzen according to this 380-page (PDF) sheet. The latency is typically overcome by computing up to 8 different AES blocks in parallel. The throughput is the sum of the reciprocal throughputs of: 1 XOR (initial key addition), 9 AESENC instructions, 1 AESENCLAST instruction, 1 ADD (AES input generation), the 11 round keys are assumed to be stored in XMM registers for all the computations and only computed + loaded once per re-keying.

Answer 2

It's been years since I've done Monte Carlo calculations (though it was more recent than the 90's!), so hopefully the information given below is still reasonably up-to-date. I've also had reason to look into pseudo-random number generation in the last few years, for other algorithms.

My question: Are there any other modern random number generators better than R250?

Yes. In fact, it's been known for a while that R250 has some bad properties and should be avoided. See e.g. Ziff, "Four-tap shift-register-sequence random-number generators," Computers in Physics 12, 385 (1998) (alternative arXiv link) and references therein. Among the issues with R250 (or in generalized notation R(103,250)) Ziff mentions:

More recently, Ferrenberg et al. [9] found that R(103, 250) leads to results being more than 100 standard deviations from the (known) true values, in simulations of the Ising model with the Wolff cluster-flipping Monte-Carlo algorithm.

The basic problem of two-tap generators R($a, b$) is that they have a built-in three-point correlation between $x_n$, $x_{n−a}$, and $x_{n−b}$, simply given by the generator itself, such that if any two of the $_n$ are known, the third follows directly from the recursion $x_n = x_{n−a} \otimes x_{n−b}$. While these correlations are spread over the size $p = \max(a, b, c, . . .)$ of the generator itself, they can evidently still lead to significant errors. These three-point correlations have been recently brought out clearly in a simulation by Schmid and Wilding

Other problems with this generator are also known. Compagner and Hoogland have shown how a pattern of all 1’s in the initialization string leads to complex (and beautiful) pattern of subsequent bits that persists for a surprisingly long time. Shchur et al. showed that, if an event occurs with a probability close to one (such as 31/32), it is not too unlikely for say 249 successive true outcomes to occur, which then leads to a very serious error at the 250-th step, when the R(103, 250) generator is used.

(by "good", I mean: low correlation, long period and fast to code and run)

The most popular multipurpose RNG these days is the Mersenne Twister, introduced in Matsumo, and Nishimura, "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator," ACM Transactions on Modeling and Computer Simulation 8, 3-30 (1998) (alternative non-paywall link). It has a long period, is reasonably fast, and is widely implemented in libraries. This works fine for sequential code, but it's easy to make mistakes when implementing it in a parallel environment, e.g. for GPU implementations. In such cases, make sure to use a suitable implementation. There's a good discussion of the requirements in L'Ecuyer, Munger, Oreshkin, and Simard, "Random numbers for parallel computers: Requirements and methods, with emphasis on GPUs," Mathematics and Computers in Simulation 135, 3-17 (2017).

Another more recent option is the so-called WELL generator (Well Equidistributed Long-period Linear) introduced in Panneton, L'Ecuyer, and Matsumoto, "Improved long-period generators based on linear recurrences modulo 2," ACM Transactions on Mathematical Software 32, 1-16 (2006) (alternative non-paywall link). In the Apache Commons Math user guide it is stated that,

The WELL generators are more modern then MersenneTwister (the paper describing than has been published in 2006 instead of 1998) and fix some of its (few) drawbacks. If initialization array contains many zero bits, MersenneTwister may take a very long time (several hundreds of thousands of iterations to reach a steady state with a balanced number of zero and one in its bits pool). So the WELL generators are better to escape zeroland as explained by the WELL generators creators.

For a general Monte Carlo simulation I would probably pick one of these two, since they've been tested and explored in practical applications for a while now. Pragmatically speaking, a reviewer probably won't challenge you on their use. Also (though I'm not sure how much to read into it), the Mersenne Twister has excitingly been found to compare well to a physical, quantum-based random number generator, see Ghersi, Parakh, and Mezei, "Comparison of a quantum random number generator with pseudorandom number generators for their use in molecular Monte Carlo simulations," Journal of Computational Chemistry 38, 2713-2720 (2017).

Anyway, probably neither is the best at any one of the qualities you're looking for. If you need really fast, something out of the Xorshift or Permuted congruential generator (PCG) families are probably the way to go. (While PCG has good statistical properties, some drawbacks are discussed in this Stats.SE question.) But if I were to go that route, I would consider checking the results using also a more conventional PRNG. That is probably a good idea anyway to make sure results don't depend on the specific choice of PRNG.

Final note: It should not be assumed that the mentioned pseudorandom number generators are cryptographically safe.