Same random sequence on GPU and CPU?

is there a way to generate the same random sequence of numbers on a GPU as on the CPU? i was hoping i could just use rng = Random.default_rng() as the generator for both GPU and CPU vectors, but while it doesn’t emit an error, it doesn’t give the same sequence. what’s curious is it always works fine for short vectors, but never for long ones.

doesn’t work for long vectors:

julia> using Random, CUDA

julia> CUDA.allowscalar(true)

julia> rng = Random.default_rng()
TaskLocalRNG()

julia> sz=3000
3000

julia> v = Vector{Float64}(undef, sz);

julia> Random.seed!(rng, 1)     ### seed
TaskLocalRNG()

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [1.3124259911967384, -0.013263403885168697, -1.0268774653294301]

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [1.4277456563925504, -1.0891864209485216, -1.29325857771683]

julia> v = CuVector{Float64}(undef, sz);

julia> Random.seed!(rng, 1)        ### re-seed
TaskLocalRNG()

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [-0.2627388126694136, 0.6279014918085116, 1.0570439330176862]

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [0.47525705350476416, -0.8981703138311293, 0.5934767517052868]

but it does for short:

julia> sz=3
3

julia> v = Vector{Float64}(undef, sz);

julia> Random.seed!(rng, 1)         ### seed
TaskLocalRNG()

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [-0.2627388126694136, 0.6279014918085116, 1.0570439330176862]

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [0.05101573888172962, 1.674880612231772, 1.117523429442645]

julia> v = CuVector{Float64}(undef, sz);

julia> Random.seed!(rng, 1)         ### re-seed
TaskLocalRNG()

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [-0.2627388126694136, 0.6279014918085116, 1.0570439330176862]

julia> randn!(rng, v);

julia> @info v[1:3]
[ Info: [0.05101573888172962, 1.674880612231772, 1.117523429442645]

this is with 1.7-beta4 and 0-day master CUDA.jl

You need to re-seed the RNG. This works, but is slow as it just generates the random numbers on the CPU and uploads them to the GPU, the generators are not compatible. To use the GPU for random number generation, you need to use CUDA’s RNG.

what do you mean by “re-seed the RNG”? i called Random.seed!() twice in each example.

i’m fine with it being slow. just using this for tests to make sure the refactoring of CPU code to use the GPU yields the same results.

Sorry, misinterpreted your code snippet. But isn’t it the Array case that’s inconsistent here?

julia> sz = 1; a = Vector{Float64}(undef, sz); Random.seed!(rng, 1); Random.rand!(rng, a); a[1]
0.12187431157074147

julia> sz = 7; a = Vector{Float64}(undef, sz); Random.seed!(rng, 1); Random.rand!(rng, a); a[1]
0.12187431157074147

julia> sz = 8; a = Vector{Float64}(undef, sz); Random.seed!(rng, 1); Random.rand!(rng, a); a[1]
0.8701349878780872

vs

julia> sz = 1; a = CuVector{Float64}(undef, sz); Random.seed!(rng, 1); Random.rand!(rng, a); a[1]
0.12187431157074147

julia> sz = 7; a = CuVector{Float64}(undef, sz); Random.seed!(rng, 1); Random.rand!(rng, a); a[1]
0.12187431157074147

julia> sz = 8; a = CuVector{Float64}(undef, sz); Random.seed!(rng, 1); Random.rand!(rng, a); a[1]
0.12187431157074147

i too found it curious that the first elements in a randomly generated Array depended on the length. thought this might be a bug, but figured maybe i’m just assuming too much and that the algorithm generates the Array in it’s entirety or sth. should i submit an issue?

worth noting that random Arrays of different lengths are the same in julia 1.6.1 but not julia 1.7.0-beta4

$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.6.1 (2021-04-23)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using Random

julia> sz = 3; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(a)
3-element Vector{Float64}:
 0.23603334566204692
 0.34651701419196046
 0.3127069683360675

julia> sz = 8; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(a)
8-element Vector{Float64}:
 0.23603334566204692
 0.34651701419196046
 0.3127069683360675
 0.00790928339056074
 0.4886128300795012
 0.21096820215853596
 0.951916339835734
 0.9999046588986136

julia> 

$ julia17
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.0-beta4 (2021-08-24)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using Random

julia> sz = 3; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(a)
3-element Vector{Float64}:
 0.12187431157074147
 0.6644047635643547
 0.927097461598609

julia> sz = 8; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(a)
8-element Vector{Float64}:
 0.8701349878780872
 0.006457014103689707
 0.773618934964177
 0.28204319051637505
 0.5531597766343656
 0.08793093987132328
 0.9123350298255161
 0.4730983492512397

And you can replicate the CuArray behavior using SubArray, so this looks like a Julia issue:

julia> sz = 3; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(a)
3-element Vector{Float64}:
 0.12187431157074147
 0.6644047635643547
 0.927097461598609

julia> sz = 8; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(a)
8-element Vector{Float64}:
 0.8701349878780872
 0.006457014103689707
 0.773618934964177
 0.28204319051637505
 0.5531597766343656
 0.08793093987132328
 0.9123350298255161
 0.4730983492512397

julia> sz = 3; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(view(a, :))
3-element view(::Vector{Float64}, :) with eltype Float64:
 0.12187431157074147
 0.6644047635643547
 0.927097461598609

julia> sz = 8; a = Vector{Float64}(undef, sz); Random.seed!(1); Random.rand!(view(a, :))
8-element view(::Vector{Float64}, :) with eltype Float64:
 0.12187431157074147
 0.6644047635643547
 0.927097461598609
 0.02155985966702112
 0.823003041206761
 0.9864191953343232
 0.6825847490017343
 0.16307852197722794

dang. i’ve link this thread to https://github.com/JuliaLang/julia/issues/42150#issuecomment-915427303. maybe it deserves it’s own issue though…

i’ve submitted a new issue. https://github.com/JuliaLang/julia/issues/42165