Skip to content

The PCG implementation provided by Numpy has significant, dangerous self-correlation #16313

@vigna

Description

@vigna

The PCG generator used by Numpy has a significant amount self-correlation. That is, for each sequence generated from a seed there is a large number of correlated, nonoverlapping sequences starting from other seeds. By "correlated" I mean that interleaving two such sequences and testing the result you obtain failures that did not appear in each sequence individually.

The probability that two generators out of large set of terminals get two of those sequences is nonnegligible. Why this happens from a mathematical viewpoint is well known but it is explained here in detail: http://prng.di.unimi.it/pcg.pgp (see "Subsequences within the same generator").

To show this problem directly, I wrote this simple C program reusing the Numpy code: http://prng.di.unimi.it/intpcgnumpy.c . The program takes two 128-bit states of two generators (with the same LCG constant or "stream") in the form of high and low bits, interleaves their output and writes it in binary form. Once we send it through PractRand, we should see no statistical failure, as the two streams should be independent. But if try to start from two states with the same 64 lower bits, you get:

./intpcgnumpy 0x596d84dfefec2fc7 0x6b79f81ab9f3e37b 0x8d7deae980a64ab0 0x6b79f81ab9f3e37b | stdbuf -oL ~/svn/c/xorshift/practrand/RNG_test stdin -tf 2 -te 1 -tlmaxonly -multithreaded
RNG_test using PractRand version 0.94
RNG = RNG_stdin, seed = unknown
test set = expanded, folding = extra

rng=RNG_stdin, seed=unknown
length= 128 megabytes (2^27 bytes), time= 2.2 seconds
  Test Name                         Raw       Processed     Evaluation
  BCFN(0+0,13-2,T)                  R= +27.6  p =  1.0e-13    FAIL
  BCFN(0+1,13-2,T)                  R= +68.0  p =  2.3e-34    FAIL !!!
  BCFN(0+2,13-3,T)                  R= +90.8  p =  8.8e-43    FAIL !!!
  BCFN(0+3,13-3,T)                  R=+120.6  p =  6.9e-57    FAIL !!!!
  DC6-6x2Bytes-1                    R=  +8.9  p =  4.0e-5   mildly suspicious
  DC6-5x4Bytes-1                    R= +15.7  p =  4.3e-9   very suspicious
  [Low1/8]BCFN(0+0,13-4,T)          R= +11.6  p =  4.9e-5   unusual
  ...and 1074 test result(s) without anomalies

You can even go lower—you just need the same 58 lower bits:

./intpcgnumpy 0x596d84dfefec2fc7 0x0579f81ab9f3e37b 0x8d7deae980a64ab0 0x6b79f81ab9f3e37b | stdbuf -oL ~/svn/c/xorshift/practrand/RNG_test stdin -tf 2 -te 1 -tlmaxonly -multithreaded

[...]
rng=RNG_stdin, seed=unknown
length= 32 gigabytes (2^35 bytes), time= 453 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low1/16]FPF-14+6/32:cross        R= +11.6  p =  4.0e-10   VERY SUSPICIOUS
  [Low1/32]FPF-14+6/32:cross        R= +16.5  p =  3.2e-14    FAIL
  [Low1/32]FPF-14+6/16:cross        R= +12.8  p =  3.8e-11   VERY SUSPICIOUS
  [Low1/64]FPF-14+6/64:cross        R=  +6.8  p =  4.8e-6   mildly suspicious
  [Low1/64]FPF-14+6/32:cross        R=  +6.0  p =  1.9e-5   unusual
  [Low1/64]FPF-14+6/16:cross        R=  +5.5  p =  5.8e-5   unusual
  [Low4/32]FPF-14+6/64:all          R=  +5.8  p =  5.9e-5   unusual
  [Low4/32]FPF-14+6/32:(0,14-0)     R=  +7.7  p =  1.0e-6   unusual
  [Low4/32]FPF-14+6/32:(1,14-0)     R=  +7.7  p =  9.1e-7   unusual
  [Low4/32]FPF-14+6/32:all          R=  +6.5  p =  1.3e-5   unusual
  [Low4/64]FPF-14+6/64:all          R=  +5.9  p =  5.1e-5   unusual
  [Low4/64]FPF-14+6/64:cross        R=  +8.2  p =  3.0e-7   suspicious
  [Low4/64]FPF-14+6/32:(0,14-0)     R=  +7.6  p =  1.0e-6   unusual
  [Low8/64]FPF-14+6/64:(0,14-0)     R= +17.0  p =  2.2e-15    FAIL
  [Low8/64]FPF-14+6/64:(1,14-0)     R=  +9.1  p =  5.1e-8   mildly suspicious
  [Low8/64]FPF-14+6/64:all          R= +12.7  p =  2.1e-11   VERY SUSPICIOUS
  [Low8/64]FPF-14+6/32:(0,14-0)     R= +12.8  p =  1.7e-11   VERY SUSPICIOUS
  [Low8/64]FPF-14+6/32:all          R= +11.0  p =  9.3e-10   VERY SUSPICIOUS
  ...and 1696 test result(s) without anomalies

Note that to get more the 50% probability that two generators start from two correlated seed (chosen at random) you need just about half a million generators starting at random (birthday paradox). And if you consider the probability that they do not exactly start from the same state, but have significant overlapping correlating sequences, you need much less.

Any sensible generator from the literature will not behave like that. You can choose adversarially any two starting states of MRG32k3a, SFC64, CMWC, xoshiro256++, etc., and as long as you generate nonoverlapping sequences you will not see the failures above. This is a major drawback that can pop up when a number of devices uses the generator and one assumes (as it should be) that pairwise those sequences should not show correlation. The correlation can induce unwanted behavior that is hard to detect.

Please at least document somewhere that the generator should not be used on multiple terminals or in a highly parallel environment.

The same can happen with different "streams", as the sequences generated by an LCG by changing the additive constant are all the same modulo a change of sign and an additive constant. You can see some discussion here: rust-random/rand#907 and a full mathematical discussion of the problem here: https://arxiv.org/abs/2001.05304 .

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions