Going backward in random number generation

Hi,

I have a function sample() that generate trajectories from random numbers. Before calling this function, I would do

Random.seed!(1)

My problem is that I call rand() one more time than necessary in sample() given the way my algorithm is implemented. This is part of my package PiecewiseDeterministicMarkovProcesses.jl.

Hence, if I call rand() after my function sample(), I missed a random number in the sequence of calls to rand().

My question is then: how can I go back once in the random generator?

I found empirically that this works, but it is a bit hacky:

using Random
Random.seed!(1)
rand() #0.23603334566204692
# this following call is erroneous, I want to remove it from the stack of calls
rand() #0.34651701419196046
Random.GLOBAL_RNG.idxF-=1
rand() # gives 0.34651701419196046 as expected

Thank you for your help,

Best regards

I am not sure if I understand, why would be a problem to miss a random number of the sequence? And why the solution is not to remove the extra and unnecessary call?

Why does this matter?

Mostly for debugging purposes

@rveltz I hope you’re not getting the feeling people are dismissive. As a reader, here’s what I got from this thread:

Naively, a hacky solution should be good enough for debugging purposes. After all, at some point you found the bug, you fixed it, and you remove the debugging hack. So maybe you can be more specific?

Another reason for asking that is that being able to go backward in a random number generator arguably defeats its purpose, depending on what other constraints you have. But we have no way of weighing that without a description of those other constraints!

Can you give us something to go on?

The PCG family of RNGs has an advance! function which can be used with a negative argument to backtrack. See package RandomNumbers.

Julia-1.1.0> using RandomNumbers.PCG

Julia-1.1.0> seed = 1234567
1234567

Julia-1.1.0> rng = PCGStateOneseq(UInt64, seed)
PCGStateOneseq{UInt128,Val{:XSH_RS},UInt64}(0xa10d40ffc2b1e573e589b22b2450d1fd)

Julia-1.1.0> rand(rng)
0.5716257379273757

Julia-1.1.0> rand(rng)
0.9945058856417783

Julia-1.1.0> advance!(rng, -1)
PCGStateOneseq{UInt128,Val{:XSH_RS},UInt64}(0xa84bc492b096c6e1ef3582ace3953880)

Julia-1.1.0> rand(rng)
0.9945058856417783
4 Likes

I was expected something like the answer of @greg_plowman. Indeed, I did not dive into the rand implementation, so I am not sure my hacky code does the thing I want

PRNGS don’t generally have a way to go backward in the sequence. Some are designed with an O(1) method for skipping to a point in the sequence, but that’s a rarity as far as I know.

What you’ve found with GLOBAL_RNG.idxF is an interesting implementation detail of the way that the MersenneTwister is optimized for speed: it generates random numbers in batches and stores the batch in a cache. When you get to the batch size, idxF will wrap around to 1, and your trick above will result in an index out of bounds error. It would be possible to design the cache for any PRNG to always allow going back in the sequence by up to a given number of slots but it’s a pretty unusual requirement.

4 Likes

Thank you a lot for this clarification!

I think a clean way for generic RNGs may be a wrapper type

struct BackwardRandom{T,R}
    buffer::Vector{T}
    rng::R
end

that implements rand by calling it on rng and keeping the last N draws in buffer, which you could then query, or use in place of calling rand on rng the desired number of times.

4 Likes

I don’t know if this applies for your needs, but you can copy the global RNG:

r = copy(Random.GLOBAL_RNG)
rand() # 0.29854174909400233
copy!(Random.GLOBAL_RNG, r)
rand() # 0.29854174909400233
1 Like

That’s interesting!

I think it is better to not use the global random generator. That way you don’t have to worry about another process calling rand and messing up your sequence.

julia> rng = MersenneTwister(1234);

julia> rand!(rng, zeros(5))
5-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592
 0.5662374165061859
 0.4600853424625171
 0.7940257103317943
3 Likes

I’m probably late to the party here but I did do a very similar thing as proposed by Tamas here by implementing a wrapperRNG for the RandomNumers.jl package. Instead of saving the sequence itself I saved the current state in “snapshots” so that you can then “jump” back in the sequence by finding the closest snapshot pre the point you want to move to and then advance it. If you use one of the RNG’s in the package that have small state this can be both mem and cpu efficient IMHO.

2 Likes

This seems like it could be useful if you have to abandon a computation and want to redo part of it with the same sequence. With a small state it seems quite practical.

Another approach, for MersenneTwister, would be to “jump forward” (using Future.randjump) by a number of steps a bit less than the period. But given that MersenneTwister uses caches, this is not totally straightforward to implement (see e.g. #25129). Moreover, I didn’t succeed after trying few minutes, so I may be overlooking something (I’m not sure what the exact period is, I thought it was 2^19937-1, but in the “dSFMT.h” file, it says a multiple of this number…).

I really recommend just using a different PRNG if you want to go backwards or make snapshots. RandomNumbers.jl offers various alternatives.

2 Likes

Do you also have a trick out there if one wants to backtrack random numbers, more than 1 samples ago (possibly a large number of samples!)?

Thanks!