Particle Gibbs reproducibly gives different results on x64 vs x86

[penelopeysm]

Not for the first time, x86 CI gives us mysterious problems, in that some numerical tests fail on x86 despite all RNGs being thoroughly seeded (i.e. the results should be fully deterministic).

For example, using this runtests.jl on x64 and x86 GitHub runners will give different results at the end, see https://github.com/penelopeysm/Shaymin.jl/actions/runs/12412993388/job/34653933016

using Test
using StableRNGs
using Random
using Turing

Random.seed!(468)

@testset verbose = true "Shaymin.jl" begin
    @testset "using global seed" begin
        x1 = randn(3)
        @info x1
    end

    @testset "stablerng" begin
        x2 = randn(StableRNG(468), 3)
        @info x2
    end

    @testset "pg"  begin
        @model function f(y)
            a ~ Normal(0, 1)
            y ~ Normal(a, 1)
        end
        Random.seed!(468)
        alg = PG(15)
        chain = sample(StableRNG(468), f(1.5), alg, 50; progress=false)
        @show mean(chain[:a])
    end
end

# x64
[ Info: [0.07200886749732076, -0.0740437565595174, 0.6327762377562545]
[ Info: [1.2876157288026433, -0.2953479054222536, -1.205615981210787]
mean(chain[:a]) = 0.7475257036106626

# x86
[ Info: [0.07200886749732076, -0.0740437565595174, 0.6327762377562545]
[ Info: [1.2876157288026433, -0.2953479054222536, -1.205615981210787]
mean(chain[:a]) = 0.7973086809553678

Note that:

1. The results are deterministic if run repeatedly on the same architecture, so the problem isn't that the implementation doesn't use the provided rng;
2. The first two testsets with Random.randn(10) and rand(StableRNG(468), 10) are deterministic across architectures, so it's not a mistake in the implementation of the random number generator.

[penelopeysm]

Is it related to this?

    @testset "advancedps" begin
        x4 = randn(AdvancedPS.TracedRNG(), 3)
        @info x4
    end

x64:

[ Info: [0.9001334534074001, -0.21170514711276572, 0.04622435546537583]

x86:

[ Info: [-0.29490566974498955, 0.02019167249744647, 1.7979388207251714]

TracedRNG isn't deterministic even on the same architecture, though, so that doesn't match up with previous observations.

[yebai]

TracedRNG isn't deterministic even on the same architecture, so that doesn't match up with previous observations.

We will consider transferring TracedRNG to Libtask in TuringLang/Libtask.jl#177 and improve it so it is reproducible across architectures. cc @willtebbutt

[penelopeysm]

I think we should keep this one open; we haven't figured out the root cause yet (#2449 only really plasters over it 😄 ).

[penelopeysm]

(I've been spending a little bit of time on narrowing it down, but haven't quite figured it out yet.)

[penelopeysm]

Found the culprit, finally:

AdvancedPS.jl/src/rng.jl

Lines 33 to 41 in 379336b

	"""
	split(key::Integer, n::Integer=1)
	Split `key` into `n` new keys
	"""
	function split(key::Integer, n::Integer=1)
	T = typeof(key) # Make sure the type of `key` is consistent on W32 and W64 systems.
	return T[hash(key, i) for i in UInt(1):UInt(n)]
	end

Here, key is UInt64, but UInt(1) resolves to UInt64 or UInt32 depending on the architecture, which gives different results for the same input.