Simulating quantum circuits with Yao

Specific Domains Quantum
1

Hi there! I’m new in the Julia ecosystems and my first project is trying to simulate a particular random circuit using Yao. My issue at the moment is basically the efficiency of code which seems to be quite slow. I would really appreciate any suggestion/advice on how to improve the code.

function circuit_layer(n::Int, β::Array{Float64, 1}, epsilon::Float64)
    U = chain(n)
    append!(U, chain(n, put(n, loc=>chain(rand(rng,randG)((2*rand())*epsilon*π))) for loc = 1:n))
    for j = 1:n-1
        append!(U, chain(n, subroutine(n, rot(ZZ, β[j]), (j, j+1))))
    end
    return U
end

function get_entropy_evolution(pmax::Int, spacing::Float64, n::Int, nA::Int, epsilon::Float64)
    fp = trunc(log2(pmax)) 
    times = union(trunc.([2^i for i in 1:spacing:fp]))
    ψ = Yao.uniform_state(n)
    result = zeros(length(times))
    count = 0
    for t in 1:pmax
        β = (π*0.5)*(2.0*rand(n-1) .- 1.0)
        ψ |> circuit_layer(n, β, epsilon)
        if (t in times) == true
            count = count + 1
            result[count] = entropy(state(ψ), n, nA, 1)
        end
    end
    result
end

The code works as follows: 1.) Construct one layer of the circuit, that is captured by the function circuit_layer, 2.) Then “time evolve” the state and only compute entropy at certain points in time.

I know the function svdvals is costly but in principle, it shouldn’t be a problem for system sizes around N=20

1 Like
2

Hi @rmedinar do you mind provide a runnable script? your implementation gives

julia> circuit_layer(10, rand(9), 0.2)
ERROR: UndefVarError: rng not defined
Stacktrace:
 [1] (::var"#1#2"{Int64, Float64})(loc::Int64)
   @ Main ./none:0
 [2] iterate
   @ ./generator.jl:47 [inlined]
 [3] isempty
   @ ./essentials.jl:767 [inlined]
 [4] chain
   @ ~/.julia/packages/YaoBlocks/TIlDJ/src/composite/chain.jl:41 [inlined]
 [5] circuit_layer(n::Int64, β::Vector{Float64}, epsilon::Float64)
   @ Main ./REPL[2]:3
 [6] top-level scope
   @ REPL[4]:1

and I’m not sure what is randG

3

Thanks @Roger-luo for the quick reply. I’m really sorry, this is the first time I ask a question like this one and I forgot to add the whole runnable script. Here it goes:

using MKL, Yao, Random, YaoPlots,  LinearAlgebra, Base.Threads

const SINGLE_GATES = [Rx, Ry, Rz]
@const_gate ZZ = mat(kron(Z,Z))
rng = MersenneTwister()
randG = Random.Sampler(rng, SINGLE_GATES)

function entropy(psi::Array, n::Int, nA::Int, nrenyi::Int)
    state = svdvals(reshape(psi, (2^nA, 2^(n-nA)))) .^ 2
    if nrenyi==1
        return -sum( state.*log.(state) )
	else
		return 1.0/(1.0-nrenyi) * log( sum( state.^nrenyi))
    end
end

function circuit_layer(n::Int, β::Array{Float64, 1}, epsilon::Float64)
    U = chain(n)
    append!(U, chain(n, put(n, loc=>chain(rand(rng,randG)((2*rand())*epsilon*π))) for loc = 1:n))
    for j = 1:n-1
        append!(U, chain(n, subroutine(n, rot(ZZ, β[j]), (j, j+1))))
    end
    return U
end

function get_entropy_evolution(pmax::Int, spacing::Float64, n::Int, nA::Int, epsilon::Float64)
    fp = trunc(log2(pmax)) 
    times = union(trunc.([2^i for i in 1:spacing:fp]))
    ψ = Yao.uniform_state(n)
    result = zeros(length(times))
    count = 0
    for t in 1:pmax
        β = (π*0.5)*(2.0*rand(n-1) .- 1.0)
        ψ |> circuit_layer(n, β, epsilon)
        if (t in times) == true
            count = count + 1
            result[count] = entropy(state(ψ), n, nA, 1)
        end
    end
    result
end
4 
julia> @time get_entropy_evolution(10, 0.1, 20, 2, 0.1);
  5.816009 seconds (36.36 k allocations: 20.596 GiB, 8.29% gc time)

If you change subroutine to put, it will be much faster

julia> function circuit_layer(n::Int, β::Array{Float64, 1}, epsilon::Float64)
           U = chain(n)
           append!(U, chain(n, put(n, loc=>chain(rand(rng,randG)((2*rand())*epsilon*π))) for loc = 1:n))
           for j = 1:n-1
               append!(U, chain(n, put(n, (j,j+1)=>rot(ZZ, β[j]))))
           end
           return U
       end

julia> @time get_entropy_evolution(10, 0.1, 20, 2, 0.1);
  0.949597 seconds (21.34 k allocations: 129.222 MiB, 0.38% gc time)

Subroutine is for computing a chunk of subprogram, it is not designed for applying a two qubit gate. This is why it is so slow.

Also, I think the time for computing SVD is comparable to computing the circuit. If you check the profile

@profile get_entropy_evolution(10, 0.1, 20, 2, 0.1)

you will see computing the circuit is even more costly than computing svdvals. This is because the svd matrix has the same length as the state vector, while you applied a lot of gates.

3 Likes
5

Wow, thanks! I just checked and the improvement is quite noticeable. Sure, computing svd seems comparable to computing the circuit. In principle, I’m interested in reaching large circuit depth, i.e., pmax~ O(10^3) so I think that’s the bottle neck of the computation beyond the usual Hilbert space dimension.

1 Like