Version 2.1, March 2018
Thanks to the excellent paper of L'Ecuyer, Wambergue and Bourceret,
"Spectral Analysis of the MIXMAX Random Number Generators"
we now have the ability to determine the effective resolution of the generator in higher dimensions.
Until now, we have been giving strict mathematical guarantees about convergence of Monte-Carlo integrations only up to
a dimension D lesser or equal to the dimension of the generator itself, in D <= N.
The paper of L'Ecuyer contains a number of interesting results. Fortunately, of the existing generators,
the N=17 generator passes most of the criteria as formulated in this paper, while the other two, N=8 and N=240
can in theory give incorrect results, under some astronomically rare conditions. A closer examination, below, shows that
the N=17 and N=240 generators are competitive in terms of the spectral index in corresponding dimensions with the best available on the market. Armed with the new spectral analysis from L'Ecuyer etal, we find that we can do even better than before and present three new generators for the same set of N.
The spectral index concerns the lattice structure of the linear congruential system, and measures the spacing
between successive planes which cover the generated points. The order in which the points are visited is not of concern
to the spectral test, therefore it is best thought of as measuring the effective resolution and not randomness. It is the maximal
size of a hypercube which can be squeezed in between the planes covered by the generator, so in some sense it is the maximal gap in the coverage of the points. The smaller these gaps in multidimensional coverage, the better.
The complete theory together with all the proofs and advanced theorems is in Knuth's book.
The generators which L'Ecuyer himself recommends, known as MRG32k3 and MRG63k3, definitely have excellent
randomness and spectral properties. These are far superior to all the obsolete stuff which is still shipped in various software suites and even made its way into the ISO C++ standard.
We will make our comparisons with these MRGs. To begin with, the MRG32k3 has a native resolution of 32 bits, which is 10^-10. Thus, the MRG32k3 and all other 32-bit architecture generators, which is the majority of those in wide use, have a maximum resolution of 10^-10 in the best case. This caused failure in a recent large-scale particle-physics simulation, underlying reason being the inability to represent the random arrival times of particle hits on the picosecond scale with any of the 32-bit RNGs. MIXMAX and MRG63k3 have the resolution matching and exceeding that of the double-precision numbers, 10^-16, in dimensions up to the titular dimension N (the MRGs of L'Ecuyer have N=k=3).
Furthermore, the spectral index of all linear generators rapidly worsens with dimension, so that e.g. an LCG with a 31 or 32-bit modulus cannot possibly have the spectral index better than about 0.1 in any dimension greater than D=10. Typically searches of optimal multipliers have been used to find a multiplier which approaches this best-case value, but obviously an effective resolution of say 0.1 is not good enough regardless if its the best possible for some modulus. For MRG32k3 the situation is much better than for old-time LCGs, but even there the spectral index deteriorates above 0.1 in D=60.
Let us see the situation for MIXMAX generators. The Propositions 1-2 had been taken care of long time ago by
discarding the first element of the MIXMAX vector. Proposition 3 is closely related to the sparseness
of the inverse of the MIXMAX matrix, and this was the reason for the introduction of the three-parameter family in 2015
with substantially better properties. According to Proposition 4, the bound on the spectral index of the three-parameter MIXMAXes
such as the N=17, N=8 and N=240 generally goes as the inverse of the parameter m. We have the table
N=17 N=240 N=8
Prop4 1.4E-11 4.4E-16 1.1E-16
The resolution of the double-precision numbers being 10^-16, we need not worry if the effective resolution of an RNG
is better than that number. Thus, Proposition 4 is of no practical consequence whatsoever
for N=8 and N=240 and only marginally so for N=17. The minimum dimension in which this has an effect is D=N,
so that the bound on the spectral index due to Prop4 is better than that of the competing MRGs in all these cases, being approximately 10^-4 in D=17 for MRG32k3 and 10^-7 for MRG63k3.
Proposition 6 establishes a subtle deficiency in the currently shipping generators, once again in dimensions greater than N.
This is a number-theoretical effect, which exists if some small multiple of the parameter m is close to the modulus p.
The table of spectral indexes according to the formula goes as:
N=17 N=240 N=8
Prop6A 1.1E-8 3.6E-4 1.4E-3
The value for the N=17 generator is probably too small to affect any real simulation. There are just not many real-world problems where the quantity of interest varies substantially when the value of several (at least 5) of the input variables changes only by 10^-8. For N=240, the deficiency is more worrisome, but keep in mind that both MRGs have a resolution worse than 0.1 in such high dimensions. The N=8 generator is inferior to the MRGs of L'Ecuyer.
The challenge before us is clear: create MIXMAX generators with spectral indexes below the
threshold of resolution of double-precision numbers, which is 10^-16.
We have made a search of the parameter space of MIXMAX to add the new criterion of Proposition 6.
The new parameter sets, which allow to push the spectral index at or below the resolution of double-precision numbers
are as follows
N=19: m=2^53+2^45+1, s=0
N=240: m=3141592653589847, s=0
N=8: m=2^51+2^40+1, s=0
The corresponding new table of spectral indexes is:
N=19 N=240 N=8
Prop4 1.1E-16 3.2E-16 4.4E-16
Prop6A 1.1E-16 NA 8.8E-16
The relations given in the paper of L'Ecuyer do not constitute an exhaustive list of all those that could potentially produce appreciable spacing between hyperplanes, but we have checked that all the other possible relations, such as those analogous to Prop6B, remain subdominant for these new parameter sets. In other words, both Propositions 4 and 6A remain applicable to calculate the spectral index for these specific parameter sets (but not in general for arbitrary m and s). In dimensions exceeding 2*N the spectral index will gradually deteriorate, and other techniques are required to find the spectral indexes, this work will be published separately.
Because these remain in the family of the well-understood MIXMAX three-parameter matrix A(N,s,m) as described in our 2015 paper, the required changes to the existing code to realize these generators are minimal, and we are releasing these as MIXMAX version 2.1. The performance remains very good. We thank L'Ecuyer, Wambergue and Bourceret for making their results available to the community which allows us to push forward the boundaries of whats possible.