new vqe and hadamard test

Author: GiggleLiu

examples/VQE.jl

@@ -0,0 +1,17 @@
+using Yao, Yao.Blocks
+using QuAlgorithmZoo
+using KrylovKit
+
+function ed_groundstate(h::MatrixBlock)
+    E, V = eigsolve(h |> mat, 1, :SR, ishermitian=true)
+    println("Ground State Energy is $(E[1])")
+    register(V[1])
+end
+
+N = 5
+c = random_diff_circuit(N, N, [i=>mod(i,N)+1 for i=1:N], mode=:Merged) |> autodiff(:QC)
+dispatch!(c, :random)
+hami = heisenberg(N)
+ed_groundstate(hami)
+
+vqe_solve(c, hami)

src/HadamardTest.jl

@@ -0,0 +1,46 @@
+export hadamard_test, hadamard_test_circuit, swap_test, swap_test_circuit, singlet_block, state_overlap_circuit
+
+"""
+see WiKi.
+"""
+function hadamard_test_circuit(U::MatrixBlock{N}) where N
+    chain(N+1, put(N+1, 1=>H),
+        control(N+1, 1, 2:N+1=>U),  # get matrix first, very inefficient
+        put(N+1, 1=>H)
+        )
+end
+
+function hadamard_test(U::MatrixBlock{N}, reg::AbstractRegister) where N
+    c = hadamard_test_circuit(U)
+    reg = join(reg, zero_state(1))
+    expect(put(N+1, 1=>Z), reg |> c)
+end
+
+swap_test_circuit() = hadamard_test_circuit(SWAP)
+swap_test(reg::AbstractRegister) = hadamard_test(SWAP, reg)
+
+function singlet_block(::Type{T}, nbit::Int, i::Int, j::Int) where T
+    unit = chain(nbit)
+    push!(unit, put(nbit, i=>chain(XGate{T}(), HGate{T}())))
+    push!(unit, control(nbit, -i, j=>XGate{T}()))
+end
+
+singlet_block(nbit::Int, i::Int, j::Int) = singlet_block(ComplexF64, nbit, i, j)
+singlet_block() = singlet_block(2,1,2)
+
+"""
+Estimation of overlap between multiple density matrices.
+
+PRL 88.217901
+"""
+function state_overlap_circuit(nbit::Int, nstate::Int, ϕ::Real)
+    N = nstate*nbit + 1
+
+    chain(N, put(N, 1=>H),
+        put(N, 1=>shift(ϕ)),
+        chain(N, [chain(N, [control(N, 1, (i+(k*nbit-nbit)+1, i+k*nbit+1)=>SWAP) for i=1:nbit]) for k=1:nstate-1]),  # get matrix first, very inefficient
+        put(N, 1=>H)
+        )
+end
+
+Yao.mat(ρ::DensityMatrix{1}) = dropdims(state(ρ), dims=3)

src/QCBM.jl

@@ -46,7 +46,8 @@ gradient of MMD two sample test loss, `db` must be contained in qcbm.
 `p0` is current probability distribution.
 """
 function mmdgrad(qcbm::QCBM, db::AbstractDiff; p0::Vector)
-    vstatdiff(()->psi(qcbm), db, Vstat(kmat(qcbm.kernel)), p0=p0) - 2*vstatdiff(()->psi(qcbm), db, Vstat(kmat(qcbm.kernel)*qcbm.ptrain))
+    statdiff(()->probs(qcbm) |> as_weights, db, StatFunctional(kmat(qcbm.kernel)), initial=p0 |> as_weights) -
+        2*statdiff(()->probs(qcbm) |> as_weights, db, StatFunctional(kmat(qcbm.kernel)*qcbm.ptrain))
 end
 
 """

src/QuAlgorithmZoo.jl

@@ -28,6 +28,7 @@ include("Grover.jl")
 include("PhaseEstimation.jl")
 include("HHL.jl")
 include("hamiltonian_solvers.jl")
+include("HadamardTest.jl")
 
 
 end # module

src/hamiltonian_solvers.jl

@@ -1,4 +1,4 @@
-export heisenberg, iter_groundstate!, itime_groundstate!
+export heisenberg, iter_groundstate!, itime_groundstate!, vqe_solve!
 
 """
     heisenberg(nbit::Int; periodic::Bool=true)
@@ -29,7 +29,7 @@ end
 """
     itime_groundstate!({reg::AbstractRegister}, h::MatrixBlock; τ::Int=20, tol=1e-4) -> AbstractRegister
 
-project wave function to ground state by exp(-hτ). `tol` is for `expmv`.
+Imaginary time evolution method to get ground state, i.e. by projecting wave function to ground state by exp(-hτ). `tol` is for `expmv`.
 """
 itime_groundstate!(h::MatrixBlock; τ::Real=20, tol=1e-4) = reg -> itime_groundstate!(reg, h; τ=τ, tol=tol)
 function itime_groundstate!(reg::AbstractRegister, h::MatrixBlock; τ::Int=20, tol=1e-4)
@@ -44,3 +44,24 @@ function itime_groundstate!(reg::AbstractRegister, h::MatrixBlock; τ::Int=20, t
     reg
 end
 
+# a patch for Yao.expect, to make it faster
+function Yao.expect(op::AddBlock, reg::AbstractRegister{1})
+    sum(opi->expect(opi, reg), op)
+end
+
+"""
+    vqe_solve!(circuit::MatrixBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) -> circuit
+
+variational quantum eigensolver, faithful simulation with optimizer Adam(lr=0.01).
+"""
+function vqe_solve!(circuit::MatrixBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) where N
+    optimizer = Adam(lr=0.01)
+    dbs = collect(circuit, AbstractDiff)
+    params = parameters(circuit)
+    for i = 1:niter
+        grad = opdiff.(()->zero_state(N) |> circuit, dbs, Ref(hamiltonian))
+        dispatch!(circuit, update!(params, grad, optimizer))
+        println("Step $i, Energy = $(expect(hamiltonian, zero_state(N) |> circuit))")
+    end
+    circuit
+end

test/HadamardTest.jl

@@ -0,0 +1,41 @@
+using Yao, Yao.Blocks
+using Test
+using LinearAlgebra
+using QuAlgorithmZoo
+
+@testset "state overlap" begin
+    reg1 = rand_state(3) |> focus!(1,2)
+    rho1 = reg1 |> ρ
+    reg2 = rand_state(3) |> focus!(1,2)
+    rho2 = reg2 |> ρ
+    reg3 = rand_state(3) |> focus!(1,2)
+    rho3 = reg3 |> ρ
+    desired = tr(mat(rho1)*mat(rho2))
+    c = state_overlap_circuit(2, 2, 0)
+    res = expect(put(5, 1=>Z), join(join(reg2, reg1), zero_state(1)) |> c) |> tr
+    @test desired ≈ res
+    desired = tr(mat(rho1)*mat(rho2)*mat(rho3)) |> real
+    c = state_overlap_circuit(2, 3, 0)
+    res = expect(put(7, 1=>Z), reduce(⊗, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real
+    @test desired ≈ res
+    desired = tr(mat(rho1)*mat(rho2)*mat(rho3)) |> imag
+    c = state_overlap_circuit(2, 3, -π/2)
+    res = expect(put(7, 1=>Z), reduce(⊗, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real
+    @test desired ≈ res
+end
+
+@testset "hadamard test" begin
+    nbit = 4
+    U = put(nbit, 2=>Rx(0.2))
+    reg = rand_state(nbit)
+
+    @test hadamard_test(U, reg) ≈ real(expect(U, reg))
+
+    reg = zero_state(2) |> singlet_block()
+    @test swap_test(reg) ≈ -1
+
+    reg = zero_state(2)
+    @test swap_test(reg) ≈ 1
+    reg = product_state(2, 0b11)
+    @test swap_test(reg) ≈ 1
+end

test/PhaseEstimation.jl

@@ -40,7 +40,7 @@ end
     reg = apply!(join(reg2, reg1), circuit)
 
     # measure
-    res = breflect(M, measure(focus!(copy(reg), 1:M), 10)[1]) / (1<<M)
+    res = breflect(M, measure(focus!(copy(reg), 1:M); nshot=10)[1]) / (1<<M)
 
     @test res ≈ ϕ
     @test apply!(reg, circuit |> adjoint) ≈ join(reg2, reg1)

test/hamiltonian_solvers.jl

@@ -18,5 +18,14 @@ using QuAlgorithmZoo
     reg |> itime_groundstate!(h, τ=20)
     EG = expect(h, reg)/nbit/4
     @test isapprox(EG, -0.4564, atol=1e-4)
-end
 
+    # using VQE
+    N = 4
+    h = heisenberg(N)
+    E = eigen(h |> mat |> Matrix).values[1]
+    c = random_diff_circuit(N, 5, [i=>mod(i,N)+1 for i=1:N], mode=:Merged) |> autodiff(:QC)
+    dispatch!(c, :random)
+    vqe_solve!(c, h)
+    E2 = expect(h, zero_state(N) |> c)
+    @test isapprox(E, E2, atol=1e-1)
+end

test/runtests.jl

@@ -36,3 +36,7 @@ end
 @testset "hamiltonian solvers" begin
     include("hamiltonian_solvers.jl")
 end
+
+@testset "hadamard test" begin
+    include("HadamardTest.jl")
+end