Inconsistent results between Yao and Qiskit

Hi everyone. I have been trying to simulate the growth of entanglement entropy in random unitary circuits using Yao for a while now. When comparing with old results from my Qiskit code I found out about some inconsistencies in the results, even though I basically translated my Qiskit code to Yao. Moreover, when sampling over unitaries, I get NaNs for the entropy. I think this could be associated with some precision issues but I don’t know how to tackle them. Let me explain with a small example code:

using Yao, YaoExtensions
using Random
using QuantumInformation, SharedArrays, Statistics, Distributed

function entropy(psi, n, nA)
    state = svdvals(reshape(psi, (2^nA, 2^(n-nA)))) .^ 2
    return -sum( state.*log.(state) )
end

# Circuit definition: 1.) Random single-qubit rotation gates and 2.) Controlled-Z gates as entanglers
# I used periodic boundary conditions in this particular case

function ruc(n::Int, θ::Array{Float64, 1})
    U = chain(n)
    append!(U, chain(n, put(n, loc=>chain(rand(rng,randG)( θ[loc] ))) for loc = 1:n))
    for j = 1:n-1
         append!(U, chain(n, cz(n, j, j+1)))
    end
    append!(U, chain(n, cz(n, 1)))
    return U
end

# pmax = circuit depth, n = num of qubits, nA=2, and eps controls the angles of
# the single-qubit rotation gates

function entropyEvolution(pmax::Int, n::Int, nA::Int, eps::Float64 )
    ψ = Yao.uniform_state(n)
 
    ent = Float64[]
    θ = (2*rand(n) .- 1)*(π*eps)
    
    for t in 1:pmax
        
        U = ruc(n, θ)
        ψ |> U
        
        push!(ent, entropy(state(ψ), n, nA))
    end
   ent
end

test_cz = SharedArray{Float64, 2}(100, p);
test_cz_ϵ = SharedArray{Float64, 2}(100, p);

@distributed for i in 1:100
    test_cz[i,:] = entropyEvolution(p, 10, 2, 1.0)
end

@distributed for i in 1:100
    test_cz_ϵ[i,:] = entropyEvolution(p, 10, 2, 0.05)
end

After running this code I get something like this for the entanglement entropy
Screenshot 2021-06-22 at 15.12.08
However, running the same code in Qiskit doesn’t show such oscillations in the entropy as well as no NaN values. Am I missing something here? Any help in immensely appreciated.

from a first glance, your randG seems to be something missing, I assume it’s a random gate set. not sure which gate set it is, but if you are not using a constant gate set, then this could be the problem since you aren’t actually sampling random distributed arbitrary unitaries but just a list of gates.

But again (like your previous questions…), your code doesn’t run, I’m not sure where does the NaNs actually come from.


update: by setting const randG = [Rx, Ry, Rz] and p = 10 I can run your code, but this doesn’t give me NaNs, the results do not contain NaN, I’m not sure about the oscillations due to unfamiliarity of what you are doing however

julia> any(isnan.(test_cz_ϵ))
false

julia> any(isnan.(test_cz))
false

Sorry again, I seem to always forget to add some variable declaration so that the code runs without issue. At the moment, I don’t have any explanation for the oscillations in the entanglement. Regarding the NaNs, I definitely continued getting them in my simulations. Now, since their appearance is always at a small circuit depth, I thought this could be associated with numerical precision. For example, if I get two pure build-in states in Yao:

using Yao, LinearAlgebra

function entropy(psi, n, nA)
    state = svdvals(reshape(psi, (2^nA, 2^(n-nA)))) .^ 2
    return -sum( state.*log.(state) )
end

ψ1= Yao.zero_state(12)
ψ2 = Yao.uniform_state(12)

@show entropy(ψ1, 12, 2)
@show entropy(ψ2, 12, 2)

I get a NaN as result in the first case and a very small negative value for the second state. As a simple test, I just redefined the entropy function such that only singular values above a certain threshold are considered, which luckily seems to fix the NaN issue.