Parallel Mersenne Twister

#1

Like many people, I’m pretty excited about the new multithreading capabilities in 1.3. But there’s one thing in the announcement that I’m not sure about, so it seems worth some discussion.

It says,

The approach we’ve taken with Julia’s default global random number generator ( rand() and friends) is to make it thread-specific. On first use, each thread will create an independent instance of the default RNG type (currently MersenneTwister ) seeded from system entropy.

I haven’t worked much with parallel random number generation, but ever since Hellekalek’s Don’t Trust Parallel Monte Carlo I’ve had a strong prior of suspicion when things sound too easy.

There are some ways to address this, for example O’Neill, 2014.

If there really is a problem with the current parallel RNG approach (is there?), what criteria should be most important when considering what approach to make the default? Which alternatives should be considered?

1 Like
Brew a Parallel RNG?
#2

MersenneTwister has an incredible period length of

julia> big(2)^19937-1
43154247973881626480552355163379198390539350432267115051652505414033306801376580911304513629318584665545269938257648835317902217334584413909528269154609168019007875343741396296801920114486480902661414318443276980300066728104984095451588176077132969843762134621790396391341285205627619600513106646376648615994236675486537480241964350295935168662363909047948347692313978301377820785712419054474332844529183172973242310888265081321626469451077707812282829444775022680488057820028764659399164766265200900561495800344054353690389862894061792872011120833614808447482913547328367277879565648307846909116945866230169702401260240187028746650033445774570315431292996025187780790119375902863171084149642473378986267503308961374905766340905289572290016038000571630875191373979555047468154333253474991046248132504516341796551470575481459200859472614836213875557116864445789750886277996487304308450484223420629266518556024339339190844368921018424844677042727664601852914925277280922697538426770257333928954401205465895610347658855386633902546289962132643282425748035786233580608154696546932563833327670769899439774888526687278527451002963059146963875715425735534475979734463100678367393327402149930968778296741391514599602374213629898720611431410402147238998090962818915890645693934483330994169632295877995848993366747014871763494805549996163051541225403465297007721146231355704081493098663065733677191172853987095748167816256084212823380168625334586431254034670806135273543270714478876861861983320777280644806691125713197262581763151313596429547763576367837019349835178462144294960757190918054625114143666384189433852576452289347652454631535740468786228945885654608562058042468987372436921445092315377698407168198376538237748614196207041548106379365123192817999006621766467167113471632715481795877005382694393400403061700457691135349187874888923429349340145170571716181125795888889277495426977149914549623916394014822985025331651511431278802009056808456506818877266609831636883884905621822262933986548645669080672191704740408891349835685662428063231198520436826329415290752972798343429446509992206368781367154091702655772727391329424277529349082600585884766523150957417077831910016168475685658673192860882070179760307269849987354836042371734660257694347235506301744118874141292438958141549100609752216882230887611431996472330842380137110927449483557815037586849644585749917772869926744218369621137675101083278543794081749094091043084096774144708436324279476892056200427227961638669149805489831121244676399931955371484012886360748706479568669048574782855217054740113945929622177502575565811067452201448981991968635965361551681273982740760138899638820318776303668762730157584640042798880691862640268612686180883874939573818125022279689930267446255773959542469831637863000171279227151406034129902181570659650532600775823677398182129087394449859182749999007223592423334567850671186568839186747704960016277540625331440619019129983789914712515365200336057993508601678807687568562377857095255541304902927192220184172502357124449911870210642694565061384919373474324503966267799038402386781686809962015879090586549423504699190743519551043722544515740967829084336025938225780730880273855261551972044075620326780624448803490998232161231687794715613405793249545509528052518010123087258778974115817048245588971438596754408081313438375502988726739523375296641615501406091607983229239827240614783252892479716519936989519187808681221191641747710902480633491091704827441228281186632445907145787138351234842261380074621914004818152386666043133344875067903582838283562688083236575482068479639546383819532174522502682372441363275765875609119783653298312066708217149316773564340379289724393986744139891855416612295739356668612658271234696438377122838998040199739078061443675415671078463404673702403777653478173367084844734702056866636158138003692253382209909466469591930161626097920508742175670306505139542860750806159835357541032147095084278461056701367739794932024202998707731017692582046210702212514120429322530431789616267047776115123597935404147084870985465426502772057300900333847905334250604119503030001704002887892941404603345869926367501355094942750552591581639980523190679610784993580896683299297681262442314008657033421868094551740506448829039207316711307695131892296593509018623094810557519560305240787163809219164433754514863301000915916985856242176563624771328981678548246297376249530251360363412768366456175077031977457534912806433176539995994343308118470147158712816149394421276614228262909950055746981053206610001560295784656616193252269412026831159508949671513845195883217147982748879261851417819979034417285598607727220866677680426090308754823803345446566305619241308374452754668143015487710877728011086004325892262259413968285283497045571062757701421761565262725153407407625405149931989494459106414660534305378576709862520049864880961144869258603473714363659194013962706366851389299692869491805172556818508298824954954815796063169517658741420159798754273428026723452481263569157307213153739781041627653715078598504154797287663122946711348158529418816432825044466692781137474494898385064375787507376496345148625306383391555145690087891955315994462944493235248817599907119135755933382121706191477185054936632211157222920331148502487563303118018805685073569841580518118710778653953571296014372940865270407021924383167290323231567912289419486240594039074452321678019381871219092155460768444573578559513613304242206151356457513937270939009707237827101245853837678338161023397586854894230696091540249987907453461311923963852950754758058205625956600817743007191746812655955021747670922460866747744520875607859062334750627098328593480067789456169602494392813763495657599847485773553990957557313200809040830036446492219409934096948730547494301216165686750735749555882340303989874672975455060957736921559195480815514035915707129930057027117286252843197413312307617886797506784260195436760305990340708481464607278955495487742140753570621217198252192978869786916734625618430175454903864111585429504569920905636741539030968041471

So by using randjump to make sure the RNGs are spaced, they will not overlap.

The approach I use in VectorizedRNG is that I have a large number of different multipliers for the LCG component of the PCG. Each thread (or process if multiprocessing) will be given its own unique multiplier, and thus its own unique stream.

1 Like
#3

Here is an example that initializes the default RNGs using randjump (from my comment in Julia’s issue tracker):

using Future

Random.seed!(1)
for i in 2:Threads.nthreads()
    Random.THREAD_RNGs[i] = Future.randjump(Random.THREAD_RNGs[i - 1], big(10)^20)
end

(Note: this code snippet is not “thread safe”)

I think a related and more challenging question is how to make parallel algorithms composable and reproducible. Since Julia schedules tasks dynamically, executing above snippet in the beginning of your script does not guarantee reproducibility anymore if you use multithreaded algorithms touching the default RNG (or any other parallel RNG using thread-local state).

2 Likes
#4

Oh that’s interesting. Seems this could be easily avoided with thread-local RNG state. I’m sure there recent work in this direction (on my phone currently or I’d check)

#5

IIUC, Threads.@spwan can run the task in any thread. So, if your computation relies on a thread-local state, it relies on exactly how the tasks are scheduled (which is not deterministic).

1 Like
#6

It cannot be avoided with thread-local RNG state. It could be partially mitigated by using Task-local RNG state that uses seed!(rand(parent_task.RNG)) on each @spawn / @async. That way, the random streams are as predictable as the computational graph (i.e. are still non-deterministic with channels or with tasks that check whether others have computed; but code that only uses wait for synchronization has a chance).

Mersenne Twister is too large for that. If we used something like AES128-ctr or xoroshiro, then we could reasonably do that: < 64 bytes of state, <60 cycles to spin off a new instance from an old (overhead paid on every single task creation).

Regarding randjump: Is there any performance reason with Mersenne Twister for using randjump instead of just having entirely separate RNGs?

#7

Is it OK to do this? My understanding is that randjump is the only way to get “independent” RNG stream. Or at least that’s the one recommended by RNG researchers. But I have no deep understanding of this issue.

1 Like
#8

Statistical RNG research is all about using terribly bad fast PRNGs for monte carlo simulations in a way that hopefully does not distort results too much. As such, it is a question that depends on the specifics of how the RNG is implemented and of how the random numbers are used.

I come from a crypto background. There, our standard assumption is that a CSPRNG is indistinguishable from true random, within reasonable [1] limits on the amount of computation and without knowing the secret key / initial internal state. Under that assumption I can give you a simple proof that this procedure is ok: If you could run a computation where this is not OK, i.e. where the results differ from true random, then you would have demonstrated a distinguishing attack. This would lead to the world retiring the specific CSPRNG.

This is a very powerful and simple assumption, and is a useful heuristic even if it is known to be wrong. If such a procedure turns out to be not OK for a real statistical PRNG and real computation, then this would teach us something deep about the specific statistical PRNG and computation.

I don’t know whether that is OK for Mersenne Twister, nor can I say whether that is OK for xoroshiro, under reasonable assumptions on what kind of computations are done on the random numbers. Both generators are known to not be cryptographically secure, i.e. we know extremely fast computations that regenerate the internal state from a few observed outputs or the PRNG.

In other words, unqualified “use randjump” is bullshit. “use randjump instead of new states for Mersenne Twister and monte carlo simulations” or “use randjump instead of new states for xoroshiro and monte carlo simulations” might be good advice, though.

[1] “reasonable limits” vary from “not possible on realistic hardware in realistic timeframes” to “not possible between big bang and heat death of the universe, assuming current scientific understanding of crypto / complexity theory, fundamental physics and cosmology”. The claimed period of 2^20000 for Mersenne Twister is ridiculous; the visible universe between big bang and heat death is probably too small for 2^1000 computations, so the claimed period is “effectively never”. We can get “effectively never” for much fewer bits of internal state, though.

1 Like
#9

I appreciate the detailed explanation. So, using task-local CSPRNG and combining seed!(rand(parent_task.RNG)) with @spwan sounds like a very powerful pattern to have composable parallel algorithms (although I suppose it requires that the ratio of single-thread speed of non-CS PRNG and CSPRNG is much smaller than number of CPU cores; or that RNG is not bottleneck).

#10

@foobar_lv2 Thanks for the detail, I agree this is very helpful context.

Crypto applications are important, but do you think it should be a priority for general use? I expect there’s some high cost to pay, or every application would use crypto-secure PRNGs. Crypto-oriented users are likely to check details of the PRNG, while most users “just want some random numbers”. I’d think the best approach could be to have secure PRNGs available, but for the default to be faster, non-secure, but with good statistical properties. But then, that’s my application area so I do have some bias :slight_smile:

I’d be interested to hear anyone’s experience or opinions on PCG. On paper it seems the best of everything (save crypto). Are there reasons to prefer others over it?

1 Like
#11

Maybe relevant discussion.

There are the newer RNGs, like PCG and xoroshiro, that look strictly superior to mersenne twiser, but I am not an expert on that. If one has hardware support (for aes), then there are very fast cryptographic RNGs as well. Unfortunately, some architectures lack the hardware acceleration (ancient x86, some arm, maybe some atom). Vagaries of the cross-compilation build process might force generic binaries to use slow software implementations. But all modern CPU where you would realistically like to run non-trivial computations have hardware acceleration.

A fork-on-task model of RNG states would incur an overhead on every single task creation. Hence, this can be done at most for the default RNG, and only when using an RNG with fast seeding and small state, i.e. not mersenne twister (can work for PCG, xoroshiro, etc, and most cryptographic RNG).

It is arguable whether the price of default CSPRNG is worth the gain of mitigating some security flaws. Main price would be small slow-down of rand-heavy code that uses the default RNG on sane CPUs, and big slow-down on a few specific CPUs. Main gain would be to mitigate an entire class of security bugs (secondary gain is that reasoning about randomness becomes easier).

Speed is the more common concern, but security incidents have way bigger impact than a handful of wasted cycles. Really good code would be unaffected either way (if optimized for speed and random generation is a bottleneck, then the code should already use a faster RNG than mersenne twister; if security is relevant, then the code should already use a CSPRNG).

I personally think that adopting a fast CSPRNG as default, and exporting an even faster non-cryptographic RNG is a no-brainer. Most security people would agree, and some other modern languages agree as well, e.g. swift. But reasonable people can disagree on this point; and the julia user-base is more “technical computing” than swift “general purpose”, so the decision is much harder for us than for swift.

2 Likes
#12

I would go for a fast, non-crypto RNG as the default, given that Julia is mostly targeted at scientific/numerical computing.

Regardless of the defaults, my experience is that it really pays off to keep the RNG explicit as an argument to all functions I write. It may sound a bit tedious, but I have always regretted not doing this later on, as it makes parallelization and unit testing much easier later on.

3 Likes
#13

This is excellent best-practice :+1:

#14

This is a great point, I need to update Soss to do this.

#15

Counter-based RNGs (e.g. from https://github.com/sunoru/Random123.jl) also help a lot, in my experience, since they make it very easy to partition the RNG while still using a common seed.

3 Likes
RFC: Make an abstraction for custom task parallel schedulers (GSoC proposal)
#16

I think I agree with most of what foobar_lv2 has said above. But from my amateur level reading of the literature (e.g. https://en.wikipedia.org/wiki/Mersenne_Twister#cite_note-40), seed!(rand(parent_task.RNG)) is specifically what you must avoid ever doing. IIUC, this will give you N separate but correlated random number sequences. Whereas randjump will give you uncorrelated (but limited) streams. Or initialization with a CSPRNG will give you unlimited uncorrelated streams (per foobar_lv2’s “powerful and simple assumption” above). One way to achieve the latter may be to write this call as seed!(aes_hash(rand(parent_task.RNG))) as mentioned at ttps://en.wikipedia.org/wiki/Cryptographically_secure_pseudorandom_number_generator#Designs_based_on_cryptographic_primitives, though I’ve done relatively little research here so I won’t promise that’s sufficient either.

#17

I am not sure that is a concern in practice for non-crypto applications.

The periodicity of MT is approximately 10^{6000}. Assume a single draw takes 1 ns, so you get no more then 10^{8 + 9} a year. With judicious use of jumps, you should be fine for even massive parallel workloads.

#18

That’s fine, and randjump is cool, but note that our current method requires less bookeeping, is conceptually simpler, and doesn’t have that little asterisk about “judicious use” (what is that?). I actually fail to see any reason to actually use randjump unless you’re in some alternative universe where cryptographic-quality hashes don’t exist.

1 Like
#19

And the fact that there exist papers / discussions on this kind of thing makes my point: Ensure that the PRNG is indistinguishable from true random, and bury the entire academic discourse on this kind of stuff. If one chooses a weak PRNG, i.e. a PRNG that fails for certain types of computations, then all uses / schemes will forever need to ask themselves and document “Am I doing a forbidden computation? How is the set of forbidden computations propagated to downstream users of my output?”. If one instead adopts the random oracle model, then nothing is forbidden, and the amount of scientific literature that people need to read fits on a napkin (“you may mentally model the RNG as true random; if this model turns out incorrect, then you will be compensated by having discovered the cryptographic sensation of the decade”).

In order to reap these benefits (mental clarity), we don’t even need a strong CSPRNG with adequate security margins; we only need one that passes the laugh test. E.g. reduced round number / simplified key schedule AES-CTR variants are fine for that. If one forks global RNG-state on task creation, one pays <15 throughput cycles + 16 bytes on semi-modern hardware (one could e.g. have a single key shared by all threads/tasks, and only fork the counter, since 2^128 \approx \infty).

1 Like