Weird PRNG behavior

Specific Domains Statistics
1

It looks like randn!(rng, v) gives different values depending on the type not the eltype of v:

julia> using Random
julia> v1 = Vector{Float64}(undef, 200)
julia> m = Matrix{Float64}(undef, 200, 3)
julia> v2 = view(m, :, 3)
julia> randn!(MersenneTwister(42), v1)
julia> randn!(MersenneTwister(42), v2)
julia> (v1 - v2)'
1×200 adjoint(::Vector{Float64}) with eltype Float64:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  1.33325  0.0355064  1.24369

julia> versioninfo()
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) E-2288G CPU @ 3.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)

Is this expected behavior? The compatibility and reproducibility constraints of PRNGs in Julia are often surprising to me, but the sequence of 200 numbers produced simultaneously from identically and cleanly initialized PRNGs should be identical IMHO.

2

Yes, this is expected, so that Julia has some freedom to optimize certain cases. If you need a guarantee that the stream of random numbers is always the same sequence regardless of type, you should use GitHub - JuliaRandom/StableRNGs.jl: A Julia RNG with stable streams.

1 Like
3

This holds for the same number of elements of the same eltype in the same Julia session? That’s the surprising bit for me. Because that seems like it would be the same stream. I understand that there are many optimizations that can lead to e.g. randn(RNG, n) giving a different result than [rand(rng) for _ in 1:n] (because there is a difference between generating n numbers in sequence and generating n numbers at once), but this is generating the same number of the same type of elements at the same time.

This is also the case when I do a view on the original array:

julia> v3 = view(Vector{Float64}(undef, 200), :)
julia>  randn!(MersenneTwister(42), v3)
julia> (v1 - v3)'

I find that particularly surprising because it means that a view doesn’t act like indexing (up to mutability/allocation).

This is relevant for a use-case where the PRNG is user-specified, because it is now a breaking change to change internal code to use views.

4

Changing the rng stream is never considered a breaking change, at least in Julia Base.

5

I believe that what you observe is due to specializations on the specific Array type in the Random code (for optimisations). This is a conservative constrain, and this could certainly be extended to contiguous views of Array, which would then produce the same streams on arrays and views thereof. It’s just that no-one bothered so far with implemeting this optimisation for views.

1 Like