bsc-quantic / Tenet.jl

Composable Tensor Network library in Julia
https://bsc-quantic.github.io/Tenet.jl/
Apache License 2.0
24 stars 1 forks source link

Unexpected results from `Circuit` computation #210

Closed amrFRE closed 1 month ago

amrFRE commented 1 month ago

Using Yao circuits and translating them to tensor networks in Tenet, in principle we should be able to obtain statevectors performing the contraction after merging the tensor networks related to the circuit and the state ansatz, written as a tensor too. However, for the GHZ and W circuits, involving entangling gates, results are very different:

using Yao
using Tenet

n_qubits = 3
# Create the W state circuit 
phase = 2*acos(1/sqrt(3))
circuit_W = chain(
    n_qubits, 
    put(1=>Yao.Ry(phase)),               
    Yao.control(1, 2=>Yao.H),        
    Yao.control(2, 3=>Yao.X),          
    Yao.control(1, 2=>Yao.X),      
    put(1=>Yao.X)               
)

# Create the GHZ state circuit
circuit_GHZ = chain(
    n_qubits, 
    put(1=>Yao.H),               
    Yao.control(1, 2=>Yao.X),        
    Yao.control(2, 3=>Yao.X)   
)

# Choose circuit
circuit = circuit_W     #circuit_GHZ

# Statevector from Yao circuit
SV_Yao = zero_state(n_qubits)
apply!(SV_Yao, circuit)

# Statevector from Tenet
psi000 = Tenet.Quantum(Tenet.Product(fill([1, 0], n_qubits))) #|000>
SV_Tenet = reshape(Tenet.contract(merge(Tenet.Quantum(circuit), psi000')), 2^n_qubits)

println(SV_Tenet) 
println(statevec(SV_Yao))

Statevectors from Yao are correct. The GHZ statevector from Tenet seems to have the proper amplitudes but on different basis states (the overlap with the |111> state is zero). The W statevector does not have proper amplitudes. This may affect expectation values too (when computing overlaps, for example).

mofeing commented 1 month ago

You're calling merge incorrectly: you're doing $\bra{000 \mid H}$ instead of $\ket{H \mid 000}$. Unlike the math/physics notation (i.e. from-right-to-left), merge goes from-left-to-right due to the semantics of Base.merge so you should do:

merge(psi000, Tenet.Quantum(circuit))

Check out that this way both Yao and Tenet agree with the result:

julia> statevec(SV_Yao)
8-element Vector{ComplexF64}:
               -0.0 + 0.0im
 0.5773502691896258 + 0.0im
 0.5773502691896257 + 0.0im
                0.0 + 0.0im
 0.5773502691896257 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

julia> SV_Tenet = reshape(Tenet.contract(merge(psi000, Quantum(circuit))), 2^n_qubits)
8-element reshape(::Tensor{ComplexF64, 3, Array{ComplexF64, 3}}, 8) with eltype ComplexF64:
                 0.0 + 0.0im
 -0.5773502691896257 + 0.0im
  0.5773502691896258 + 0.0im
                 0.0 + 0.0im
 -0.5773502691896257 + 0.0im
                 0.0 + 0.0im
                 0.0 + 0.0im
                 0.0 + 0.0im
jofrevalles commented 1 month ago

@mofeing In the Tenet result there are negative signs, and in the Yao one there aren't...

mofeing commented 1 month ago

@mofeing In the Tenet result there are negative signs, and in the Yao one there aren't...

Ah right.

jofrevalles commented 1 month ago

As a side note, this is not a problem of Yao, since the same problem happens with Quac:

julia> using Yao

julia> n_qubits = 3

julia> circuit_GHZ = chain(n_qubits, put(1=>Yao.H), Yao.control(1, 2=>Yao.X), Yao.control(2, 3=>Yao.X))

julia> SV_Yao = zero_state(n_qubits); apply!(SV_Yao, circuit_GHZ)

julia> statevec(SV_Yao)
8-element Vector{ComplexF64}:
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.7071067811865475 + 0.0im

 julia> julia> println(statevec(SV_Yao)' * statevec(ArrayReg(bit"111"))) # expected result = 1/√2
0.7071067811865475 + 0.0im

With Tenet:

julia> using Quac; using Tenet

julia> n_qubits = 3

julia> circuit_GHZ = Circuit(n_qubits); push!(circuit_GHZ, Quac.H(1))

julia> push!(circuit_GHZ, Quac.Control(1, Quac.X(2))); push!(circuit_GHZ, Quac.Control(2, Quac.X(3)))

julia> psi000 = Tenet.Quantum(Tenet.Product(fill([1, 0], n_qubits))) #|000>

julia> psi111 = Tenet.Quantum(Tenet.Product(fill([0, 1], n_qubits))) #|111>

julia> julia> Tenet.contract(merge(psi000, Tenet.Quantum(circuit_GHZ))) |> r -> reshape(r, 2^n_qubits) # GHZ |000>
8-element reshape(::Tensor{ComplexF64, 3, Array{ComplexF64, 3}}, 8) with eltype ComplexF64:
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

julia> SV_Tenet # expected result = 1/√2
0-dimensional Tensor{ComplexF64, 0, Array{ComplexF64, 0}}:
0.0 + 0.0im
amrFRE commented 1 month ago

Sorry, yes, I thought that was implicit. I started using Yao instead of Quac because of the deprecation you told me about :)

mofeing commented 1 month ago

mmm @jofrevalles in your example there seems to be an issue with the ordering of the values, unlike with @amrFRE's example whose issue is the signs.

jofrevalles commented 1 month ago

mmm @jofrevalles in your example there seems to be an issue with the ordering of the values, unlike with @amrFRE's example whose issue is the signs.

Yes, these are different circuits.

jofrevalles commented 1 month ago

Okay, the problem should be solved with PR #212. Maybe I will add some tests there to cover for this potential problem.

By the way, while trying to solve this problem I came across this weird thing, why did we define the Y gate like this? 😂

julia> mat(Yao.Y)
2×2 LuxurySparse.SDPermMatrix{ComplexF64, Int64, Vector{ComplexF64}, Vector{Int64}}:
 0.0+0.0im  0.0-1.0im
 0.0+1.0im  0.0+0.0im

julia> Matrix(Quac.Y())
2×2 Matrix{ComplexF64}:
 0.0+0.0im  -1.0+0.0im
 1.0+0.0im   0.0+0.0im

I am going to fix that in Quac. That is why we had some different results on my previous comment @mofeing

amrFRE commented 1 month ago

Yes, it seems that it works now. I have tested it for the GHZ and W circuits without issues.

Only a couple of things:

  1. inds_order = find_permutation(merged_circ.sites) should instead be inds_order = find_permutation(merged_circ)
  2. Changes are still on the branch fix/Yao-extension. I mean that something else had to change because running your example on master gives this minus sign weird result again.

Thank you!

jofrevalles commented 1 month ago
  • inds_order = find_permutation(merged_circ.sites) should instead be inds_order = find_permutation(merged_circ)

Oops, true! I messed up the argument here, thanks!

2. Changes are still on the branch fix/Yao-extension. I mean that something else had to change because running your example on master gives this minus sign weird result again.

Yes, I will wait until @mofeing wakes up and approves this PR. Maybe I will add a small test set but I think it is ready to merge.

Thanks, @amrFRE!