diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index cc508437..ab7ed937 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -40,6 +40,7 @@ jobs: - YaoBlocks - YaoSym - YaoPlots + - YaoToEinsum steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 diff --git a/Makefile b/Makefile index 45e53f2d..e412c328 100644 --- a/Makefile +++ b/Makefile @@ -3,9 +3,9 @@ JL = julia --project default: init test init: - $(JL) -e 'using Pkg; Pkg.develop([Pkg.PackageSpec(path = joinpath("lib", pkg)) for pkg in ["YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "YaoAPI"]]); Pkg.precompile()' + $(JL) -e 'using Pkg; Pkg.develop([Pkg.PackageSpec(path = joinpath("lib", pkg)) for pkg in ["YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "YaoAPI", "YaoToEinsum"]]); Pkg.precompile()' init-docs: - $(JL) -e 'using Pkg; Pkg.activate("docs"); Pkg.develop([Pkg.PackageSpec(path = joinpath("lib", pkg)) for pkg in ["YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "YaoAPI"]]); Pkg.precompile()' + $(JL) -e 'using Pkg; Pkg.activate("docs"); Pkg.develop([Pkg.PackageSpec(path = joinpath("lib", pkg)) for pkg in ["YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "YaoAPI", "YaoToEinsum"]]); Pkg.precompile()' update: $(JL) -e 'using Pkg; Pkg.update(); Pkg.precompile()' @@ -22,10 +22,10 @@ test-%: $(JL) -e 'using Pkg; Pkg.test("$*")' test: - $(JL) -e 'using Pkg; Pkg.test(["YaoAPI", "YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "Yao"])' + $(JL) -e 'using Pkg; Pkg.test(["YaoAPI", "YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "Yao", "YaoToEinsum"])' coverage: - $(JL) -e 'using Pkg; Pkg.test(["YaoAPI", "YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "Yao"]; coverage=true)' + $(JL) -e 'using Pkg; Pkg.test(["YaoAPI", "YaoArrayRegister", "YaoBlocks", "YaoSym", "YaoPlots", "Yao", "YaoToEinsum"]; coverage=true)' servedocs: $(JL) -e 'using Pkg; Pkg.activate("docs"); using LiveServer; servedocs(;skip_dirs=["docs/src/assets", "docs/src/generated"])' diff --git a/Project.toml b/Project.toml index 4faa42d5..e1a7aa85 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,7 @@ YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" YaoPlots = "32cfe2d9-419e-45f2-8191-2267705d8dbc" YaoSym = "3b27209a-d3d6-11e9-3c0f-41eb92b2cb9d" +YaoToEinsum = "9b173c7b-dc24-4dc5-a0e1-adab2f7b6ba9" [weakdeps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" diff --git a/docs/make.jl b/docs/make.jl index 9a549392..7ca1b483 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -61,6 +61,7 @@ const PAGES = [ "man/cuda.md", "man/plot.md", "man/automatic_differentiation.md", + "man/yao2einsum.md", "man/simplification.md", "man/bitbasis.md", ], @@ -71,7 +72,7 @@ const PAGES = [ indigo = DocThemeIndigo.install(Yao) makedocs( - modules = [Yao, YaoAPI, YaoArrayRegister, YaoBlocks, BitBasis, YaoSym, YaoPlots, AD, Optimise], + modules = [Yao, YaoAPI, YaoArrayRegister, YaoBlocks, BitBasis, YaoSym, YaoPlots, YaoToEinsum, AD, Optimise], format = Documenter.HTML( prettyurls = ("deploy" in ARGS), canonical = ("deploy" in ARGS) ? "https://docs.yaoquantum.org/" : nothing, diff --git a/docs/src/man/yao2einsum.md b/docs/src/man/yao2einsum.md new file mode 100644 index 00000000..f8cce1f6 --- /dev/null +++ b/docs/src/man/yao2einsum.md @@ -0,0 +1,45 @@ +# Tensor network backend + +Simulating quantum circuits using tensor networks has been studied in the literature[^Markov2008][^Pan2022]. The `YaoToEinsum` package provides a convenient way to convert Yao circuits to tensor networks, which can be used for further analysis and optimization. + +## Tutorial +The main function is +```julia +yao2einsum(circuit; initial_state=Dict(), final_state=Dict(), optimizer=TreeSA()) +``` +which transforms a [`Yao`](https://github.com/QuantumBFS/Yao.jl) circuit to a tensor network that generalizes the hyper-graph (einsum notation). The return value is a `TensorNetwork` object. + +* `initial_state` and `final_state` are for specifying the initial state and final state. Left the qubits unspecified if you want to keep them as the open indices. +* `optimizer` is for specifying the contraction order optimizing algorithm of the tensor network. The default value is the `TreeSA()` algorithm that developed in [^Kalachev2021][^Liu2023]. Please check the README of [OMEinsumEinsumContractors.jl](https://github.com/TensorBFS/OMEinsumContractionOrders.jl) for more information. + +In the following example, we show how to convert a quantum Fourier transform circuit to a tensor network and contract it to +- Get the matrix representation of the circuit. +- Get the probability of measuring the zero state after applying the circuit on the zero state. + +```@repl +import Yao +using Yao.EasyBuild: qft_circuit +n = 10; +circuit = qft_circuit(n); # build a quantum Fourier transform circuit +network = Yao.yao2einsum(circuit) # convert this circuit to tensor network +reshape(Yao.contract(network), 1<0 for i=1:n]), final_state=Dict([i=>0 for i=1:n]), + optimizer=Yao.YaoToEinsum.TreeSA(; nslices=3)) # slicing technique +Yao.contract(network)[] ≈ Yao.zero_state(n)' * (Yao.zero_state(n) |> circuit) +``` + +## API +```@docs +yao2einsum +TensorNetwork +optimize_code +contraction_complexity +contract +``` + +## References +[^Pan2022]: Pan, Feng, and Pan Zhang. "Simulation of quantum circuits using the big-batch tensor network method." Physical Review Letters 128.3 (2022): 030501. +[^Kalachev2021]: Kalachev, Gleb, Pavel Panteleev, and Man-Hong Yung. "Recursive multi-tensor contraction for xeb verification of quantum circuits." arXiv preprint arXiv:2108.05665 (2021). +[^Markov2008]: Markov, Igor L., and Yaoyun Shi. "Simulating quantum computation by contracting tensor networks." SIAM Journal on Computing 38.3 (2008): 963-981. +[^Liu2023]: Liu, Jin-Guo, et al. "Computing solution space properties of combinatorial optimization problems via generic tensor networks." SIAM Journal on Scientific Computing 45.3 (2023): A1239-A1270. \ No newline at end of file diff --git a/lib/YaoAPI/src/blocks.jl b/lib/YaoAPI/src/blocks.jl index f407c273..24fd1744 100644 --- a/lib/YaoAPI/src/blocks.jl +++ b/lib/YaoAPI/src/blocks.jl @@ -305,21 +305,21 @@ Returns the most compact matrix form of given block, e.g ```jldoctest; setup=:(using Yao) julia> mat(X) -2×2 LuxurySparse.SDPermMatrix{ComplexF64, Int64, Vector{ComplexF64}, Vector{Int64}}: - 0.0+0.0im 1.0+0.0im - 1.0+0.0im 0.0+0.0im +LuxurySparse.SDPermMatrix{ComplexF64, Int64, Vector{ComplexF64}, Vector{Int64}} +(1, 2) = 1.0 + 0.0im +(2, 1) = 1.0 + 0.0im julia> mat(Float64, X) -2×2 LuxurySparse.SDPermMatrix{Float64, Int64, Vector{Float64}, Vector{Int64}}: - 0.0 1.0 - 1.0 0.0 +LuxurySparse.SDPermMatrix{Float64, Int64, Vector{Float64}, Vector{Int64}} +(1, 2) = 1.0 +(2, 1) = 1.0 julia> mat(kron(X, X)) -4×4 LuxurySparse.SDPermMatrix{ComplexF64, Int64, Vector{ComplexF64}, Vector{Int64}}: - 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im - 0.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im - 0.0+0.0im 1.0+0.0im 0.0+0.0im 0.0+0.0im - 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im +LuxurySparse.SDPermMatrix{ComplexF64, Int64, Vector{ComplexF64}, Vector{Int64}} +(1, 4) = 1.0 + 0.0im +(2, 3) = 1.0 + 0.0im +(3, 2) = 1.0 + 0.0im +(4, 1) = 1.0 + 0.0im julia> mat(kron(X, X) + put(2, 1=>X)) 4×4 SparseMatrixCSC{ComplexF64, Int64} with 8 stored entries: diff --git a/lib/YaoToEinsum/CITATION.bib b/lib/YaoToEinsum/CITATION.bib new file mode 100644 index 00000000..42462489 --- /dev/null +++ b/lib/YaoToEinsum/CITATION.bib @@ -0,0 +1,61 @@ +@article{Luo2020, + title = {Yao.jl: {Extensible}, {Efficient} {Framework} for {Quantum} {Algorithm} {Design}}, + volume = {4}, + shorttitle = {Yao.jl}, + url = {https://quantum-journal.org/papers/q-2020-10-11-341/}, + doi = {10.22331/q-2020-10-11-341}, + abstract = {Xiu-Zhe Luo, Jin-Guo Liu, Pan Zhang, and Lei Wang, +Quantum 4, 341 (2020). +We introduce \${\textbackslash}texttt\{Yao\}\$, an extensible, efficient open-source framework for quantum algorithm design. \${\textbackslash}texttt\{Yao\}\$ features generic and differentiable programming of quantum circuits. It a…}, + language = {en-GB}, + urldate = {2023-03-23}, + journal = {Quantum}, + author = {Luo, Xiu-Zhe and Liu, Jin-Guo and Zhang, Pan and Wang, Lei}, + month = oct, + year = {2020}, + note = {Publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften}, + pages = {341}, +} + +@article{Pan2022, + title = {Simulation of {Quantum} {Circuits} {Using} the {Big}-{Batch} {Tensor} {Network} {Method}}, + volume = {128}, + url = {https://link.aps.org/doi/10.1103/PhysRevLett.128.030501}, + doi = {10.1103/PhysRevLett.128.030501}, + abstract = {We propose a tensor network approach to compute amplitudes and probabilities for a large number of correlated bitstrings in the final state of a quantum circuit. As an application, we study Google’s Sycamore circuits, which are believed to be beyond the reach of classical supercomputers and have been used to demonstrate quantum supremacy. By employing a small computational cluster containing 60 graphical processing units (GPUs), we compute exact amplitudes and probabilities of 2×106 correlated bitstrings with some entries fixed (which span a subspace of the output probability distribution) for the Sycamore circuit with 53 qubits and 20 cycles. The obtained results verify the Porter-Thomas distribution of the large and deep quantum circuits of Google, provide datasets and benchmarks for developing approximate simulation methods, and can be used for spoofing the linear cross entropy benchmark of quantum supremacy. Then we extend the proposed big-batch method to a full-amplitude simulation approach that is more efficient than the existing Schrödinger method on shallow circuits and the Schrödinger-Feynman method in general, enabling us to obtain the state vector of Google’s simplifiable circuit with n=43 qubits and m=14 cycles using only one GPU. We also manage to obtain the state vector for Google’s simplifiable circuits with n=50 qubits and m=14 cycles using a small GPU cluster, breaking the previous record on the number of qubits in full-amplitude simulations. Our method is general in computing bitstring probabilities for a broad class of quantum circuits and can find applications in the verification of quantum computers. We anticipate that our method will pave the way for combining tensor network–based classical computations and near-term quantum computations for solving challenging problems in the real world.}, + number = {3}, + urldate = {2023-02-09}, + journal = {Physical Review Letters}, + author = {Pan, Feng and Zhang, Pan}, + month = jan, + year = {2022}, + note = {Publisher: American Physical Society}, + pages = {030501}, +} + +@article{Kalachev2021, + title = {Recursive {Multi}-{Tensor} {Contraction} for {XEB} {Verification} of {Quantum} {Circuits}}, + url = {http://arxiv.org/abs/2108.05665}, + abstract = {The computational advantage of noisy quantum computers have been demonstrated by sampling the bitstrings of quantum random circuits. An important issue is how the performance of quantum devices could be quantified in the so-called "supremacy regime". The standard approach is through the linear cross entropy (XEB), where the theoretical value of the probability is required for each bitstring. However, the computational cost of XEB grows exponentially. So far, random circuits of the 53-qubit Sycamore chip was verified up to 10 cycles of gates only; the XEB fidelities of deeper circuits were approximated with simplified circuits instead. Here we present a multi-tensor contraction algorithm for speeding up the calculations of XEB of quantum circuits, where the computational cost can be significantly reduced through a recursive manner with some form of memoization. As a demonstration, we analyzed the experimental data of the 53-qubit Sycamore chip and obtained the exact values of the corresponding XEB fidelities up to 16 cycles using only moderate computing resources (few GPUs). If the algorithm was implemented on the Summit supercomputer, we estimate that for the 20-cycles supremacy circuits, it would only cost 7.5 days, which is several orders of magnitudes lower than previously estimated in the literature.}, + author = {Kalachev, Gleb and Panteleev, Pavel and Yung, Man-Hong}, + year = {2021}, + note = {arXiv: 2108.05665}, + pages = {1--9}, +} + +@article{Markov2008, + title = {Simulating {Quantum} {Computation} by {Contracting} {Tensor} {Networks}}, + volume = {38}, + issn = {0097-5397}, + url = {https://epubs.siam.org/doi/abs/10.1137/050644756}, + doi = {10.1137/050644756}, + abstract = {Adiabatic quantum computation has recently attracted attention in the physics and computer science communities, but its computational power was unknown. We describe an efficient adiabatic simulation of any given quantum algorithm, which implies that the adiabatic computation model and the conventional quantum computation model are polynomially equivalent. Our result can be extended to the physically realistic setting of particles arranged on a two‐dimensional grid with nearest neighbor interactions. The equivalence between the models allows stating the main open problems in quantum computation using well‐studied mathematical objects such as eigenvectors and spectral gaps of sparse matrices.}, + number = {3}, + urldate = {2023-10-05}, + journal = {SIAM Journal on Computing}, + author = {Markov, Igor L. and Shi, Yaoyun}, + month = jan, + year = {2008}, + note = {Publisher: Society for Industrial and Applied Mathematics}, + pages = {963--981}, +} diff --git a/lib/YaoToEinsum/LICENSE b/lib/YaoToEinsum/LICENSE new file mode 100644 index 00000000..75029780 --- /dev/null +++ b/lib/YaoToEinsum/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 GiggleLiu and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/lib/YaoToEinsum/Project.toml b/lib/YaoToEinsum/Project.toml new file mode 100644 index 00000000..1f841942 --- /dev/null +++ b/lib/YaoToEinsum/Project.toml @@ -0,0 +1,28 @@ +name = "YaoToEinsum" +uuid = "9b173c7b-dc24-4dc5-a0e1-adab2f7b6ba9" +authors = ["GiggleLiu and contributors"] +version = "0.2.2" + +[deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +OMEinsum = "ebe7aa44-baf0-506c-a96f-8464559b3922" +YaoBlocks = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df" + +[weakdeps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + +[extensions] +YaoToEinsumCUDAExt = "CUDA" + +[compat] +CUDA = "4, 5" +OMEinsum = "0.8" +YaoBlocks = "0.13" +julia = "1.9" + +[extras] +SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "SymEngine"] diff --git a/lib/YaoToEinsum/README.md b/lib/YaoToEinsum/README.md new file mode 100644 index 00000000..f689653c --- /dev/null +++ b/lib/YaoToEinsum/README.md @@ -0,0 +1,14 @@ +# YaoToEinsum + +Convert [Yao](https://github.com/QuantumBFS/Yao.jl) circuit to tensor networks (einsum). + +## Installation + +`YaoToEinsum` is a [Julia language](https://julialang.org/) package. To install `YaoToEinsum`, please [open Julia's interactive session (known as REPL)](https://docs.julialang.org/en/v1/manual/getting-started/) and press ] key in the REPL to use the package mode, then type the following command + +```julia +pkg> add YaoToEinsum +``` + +## Using +Please refer to the [documentation](). \ No newline at end of file diff --git a/lib/YaoToEinsum/ext/YaoToEinsumCUDAExt.jl b/lib/YaoToEinsum/ext/YaoToEinsumCUDAExt.jl new file mode 100644 index 00000000..cf57007e --- /dev/null +++ b/lib/YaoToEinsum/ext/YaoToEinsumCUDAExt.jl @@ -0,0 +1,7 @@ +module YaoToEinsumCUDAExt +using CUDA, YaoToEinsum + +function CUDA.cu(tnet::TensorNetwork) + return TensorNetwork(tnet.code, tnet.tensors .|> CuArray) +end +end diff --git a/lib/YaoToEinsum/src/Core.jl b/lib/YaoToEinsum/src/Core.jl new file mode 100644 index 00000000..3d36cc22 --- /dev/null +++ b/lib/YaoToEinsum/src/Core.jl @@ -0,0 +1,62 @@ +""" + TensorNetwork + +A (generalized) tensor network representation of a quantum circuit. + +### Fields +* `code::AbstractEinsum`: The einsum code. +* `tensors::Vector`: The tensors in the network. +""" +struct TensorNetwork + code::AbstractEinsum + tensors::Vector +end +function Base.show(io::IO, c::TensorNetwork) + print(io, "TensorNetwork") + print(io, "\n") + print(io, contraction_complexity(c)) +end +function Base.show(io::IO, ::MIME"text/plain", c::TensorNetwork) + Base.show(io, c) +end +function Base.iterate(c::TensorNetwork, state=1) + if state > 2 + return nothing + elseif state == 1 + return (c.code, 2) + else + return (c.tensors, 3) + end +end + +""" + contract(c::TensorNetwork) + +Contract the tensor network, and return the result tensor. +""" +function contract(c::TensorNetwork) + return c.code(c.tensors...; size_info=uniformsize(c.code, 2)) +end + +""" + optimize_code(c::TensorNetwork, optimizer=TreeSA()) + +Optimize the code of the tensor network. + +### Arguments +* `c::TensorNetwork`: The tensor network. +* `optimizer::Optimizer`: The optimizer to use, default is `TreeSA()`. Please check [OMEinsumContractors.jl](https://github.com/TensorBFS/OMEinsumContractionOrders.jl) for more information. +""" +function OMEinsum.optimize_code(c::TensorNetwork, args...) + optcode = optimize_code(c.code, uniformsize(c.code, 2), args...) + return TensorNetwork(optcode, c.tensors) +end + +""" + contraction_complexity(c::TensorNetwork) + +Return the contraction complexity of the tensor network. +""" +function OMEinsum.contraction_complexity(c::TensorNetwork) + return contraction_complexity(c.code, uniformsize(c.code, 2)) +end \ No newline at end of file diff --git a/lib/YaoToEinsum/src/YaoToEinsum.jl b/lib/YaoToEinsum/src/YaoToEinsum.jl new file mode 100644 index 00000000..1b10d22b --- /dev/null +++ b/lib/YaoToEinsum/src/YaoToEinsum.jl @@ -0,0 +1,13 @@ +module YaoToEinsum + +using YaoBlocks, YaoBlocks.YaoArrayRegister, OMEinsum +using LinearAlgebra + +export yao2einsum +export TensorNetwork, optimize_code, contraction_complexity, contract +export TreeSA + +include("Core.jl") +include("circuitmap.jl") + +end diff --git a/lib/YaoToEinsum/src/circuitmap.jl b/lib/YaoToEinsum/src/circuitmap.jl new file mode 100644 index 00000000..a1c7de9e --- /dev/null +++ b/lib/YaoToEinsum/src/circuitmap.jl @@ -0,0 +1,209 @@ +struct EinBuilder{T} + slots::Vector{Int} + labels::Vector{Vector{Int}} + tensors::Vector{AbstractArray{T}} + maxlabel::Base.RefValue{Int} +end + +YaoBlocks.nqubits(eb::EinBuilder) = length(eb.slots) +function add_tensor!(eb::EinBuilder{T}, tensor::AbstractArray{T,N}, labels::Vector{Int}) where {N,T} + @assert N == length(labels) + push!(eb.tensors, tensor) + push!(eb.labels, labels) +end + +function EinBuilder(::Type{T}, n::Int) where T + EinBuilder(collect(1:n), Vector{Int}[], AbstractArray{T}[], Ref(n)) +end +newlabel!(eb::EinBuilder) = (eb.maxlabel[] += 1; eb.maxlabel[]) + +function add_gate!(eb::EinBuilder{T}, b::PutBlock{D,C}) where {T,D,C} + return add_matrix!(eb, C, mat(T, b.content), collect(b.locs)) +end +# general and diagonal gates +function add_matrix!(eb::EinBuilder{T}, k::Int, m::AbstractMatrix, locs::Vector) where T + if isdiag(m) + add_tensor!(eb, reshape(Vector{T}(diag(m)), fill(2, k)...), eb.slots[locs]) + elseif m isa YaoBlocks.OuterProduct # low rank + nlabels = [newlabel!(eb) for _=1:k] + K = rank(m) + if K == 1 # projector + add_tensor!(eb, reshape(Vector{T}(m.right), fill(2, k)...), [eb.slots[locs]...]) + add_tensor!(eb, reshape(Vector{T}(m.left), fill(2, k)...), [nlabels...]) + eb.slots[locs] .= nlabels + else + midlabel = newlabel!(eb) + add_tensor!(eb, reshape(Matrix{T}(m.right), fill(2, k)..., K), [eb.slots[locs]..., midlabel]) + add_tensor!(eb, reshape(Matrix{T}(m.left), fill(2, k)..., K), [nlabels..., midlabel]) + eb.slots[locs] .= nlabels + end + else + nlabels = [newlabel!(eb) for _=1:k] + add_tensor!(eb, reshape(Matrix{T}(m), fill(2, 2k)...), [nlabels..., eb.slots[locs]...]) + eb.slots[locs] .= nlabels + end + return eb +end +# swap gate +function add_gate!(eb::EinBuilder{T}, b::PutBlock{2,2,ConstGate.SWAPGate}) where {T} + lj = eb.slots[b.locs[2]] + eb.slots[b.locs[2]] = eb.slots[b.locs[1]] + eb.slots[b.locs[1]] = lj + return eb +end + +# projection gate, todo: generalize to arbitrary low rank gate +function add_gate!(eb::EinBuilder{T}, b::PutBlock{2,1,ConstGate.P0Gate}) where {T} + add_matrix!(eb, 1, YaoBlocks.OuterProduct(T[1, 0], T[1, 0]), collect(b.locs)) + return eb +end + +# projection gate, todo: generalize to arbitrary low rank gate +function add_gate!(eb::EinBuilder{T}, b::PutBlock{2,1,ConstGate.P1Gate}) where {T} + add_matrix!(eb, 1, YaoBlocks.OuterProduct(T[0, 1], T[0, 1]), collect(b.locs)) + return eb +end + + +# control gates +function add_gate!(eb::EinBuilder{T}, b::ControlBlock{BT,C,M}) where {T, BT,C,M} + return add_controlled_matrix!(eb, M, mat(T, b.content), collect(b.locs), collect(b.ctrl_locs), collect(b.ctrl_config)) +end +function add_controlled_matrix!(eb::EinBuilder{T}, k::Int, m::AbstractMatrix, locs::Vector, control_locs, control_vals) where T + if length(control_locs) == 0 + return add_matrix!(eb, k, m, locs) + end + sig = eb.slots[control_locs[1]] + val = control_vals[1] + for i=1:length(control_locs)-1 + newsig = newlabel!(eb) + add_tensor!(eb, and_gate(T, control_vals[i+1], val), [newsig,eb.slots[control_locs[i+1]],sig]) + sig = newsig + val = 1 + end + if !isdiag(m) + t1 = reshape(Matrix{T}(m), fill(2, 2k)...) + t2 = reshape(Matrix{T}(I, 1<1, 2=>1, 3=>0, 4=>1)` specifies a product state `|1⟩⊗|1⟩⊗|0⟩⊗|1⟩`. + - In the second interface, a state is specified as an `ArrayReg`, e.g. `Dict(1=>rand_state(1), 2=>rand_state(1))`. +If any qubit in initial state or final state is not specified, it will be treated as a free leg in the tensor network. +* `optimizer` is the optimizer used to optimize the tensor network. The default is `TreeSA()`. +Please check [OMEinsumContractors.jl](https://github.com/TensorBFS/OMEinsumContractionOrders.jl) for more information. + + +```jldoctest +julia> using Yao + +julia> c = chain(3, put(3, 2=>X), put(3, 1=>Y), control(3, 1, 3=>Y)) +nqubits: 3 +chain +├─ put on (2) +│ └─ X +├─ put on (1) +│ └─ Y +└─ control(1) + └─ (3,) Y + + +julia> yao2einsum(c; initial_state=Dict(1=>0, 2=>1), final_state=Dict(1=>ArrayReg([0.6, 0.8im]), 2=>1)) +TensorNetwork +Time complexity: 2^4.700439718141093 +Space complexity: 2^2.0 +Read-write complexity: 2^6.0 +``` +""" +function yao2einsum(circuit::AbstractBlock{D}; initial_state::Dict=Dict{Int,Int}(), final_state::Dict=Dict{Int,Int}(), optimizer=TreeSA()) where {D} + T = promote_type(ComplexF64, dict_regtype(initial_state), dict_regtype(final_state), YaoBlocks.parameters_eltype(circuit)) + vec_initial_state = Dict{Int,ArrayReg{D,T}}([k=>render_single_qubit_state(T, v) for (k, v) in initial_state]) + vec_final_state = Dict{Int,ArrayReg{D,T}}([k=>render_single_qubit_state(T, v) for (k, v) in final_state]) + yao2einsum(circuit, vec_initial_state, vec_final_state, optimizer) +end +dict_regtype(d::Dict) = promote_type(_regtype.(values(d))...) +_regtype(::ArrayReg{D,VT}) where {D,VT} = VT +_regtype(::Int) = ComplexF64 +render_single_qubit_state(::Type{T}, x::Int) where T = x == 0 ? zero_state(T, 1) : product_state(T, bit"1") +render_single_qubit_state(::Type{T}, x::ArrayReg) where T = ArrayReg(collect(T, statevec(x))) + +function yao2einsum(circuit::AbstractBlock{D}, initial_state::Dict{Int,<:ArrayReg{D,T}}, final_state::Dict{Int,<:ArrayReg{D,T}}, optimizer) where {D,T} + n = nqubits(circuit) + eb = EinBuilder(T, n) + openindices = Int[] + for k=1:n + if haskey(initial_state, k) + add_tensor!(eb, statevec(initial_state[k]), [eb.slots[k]]) + else + push!(openindices, eb.slots[k]) + end + end + add_gate!(eb, circuit) + openindices2 = Int[] + for k=1:n + if haskey(final_state, k) + add_tensor!(eb, statevec(final_state[k]), [eb.slots[k]]) + else + push!(openindices2, eb.slots[k]) + end + end + network = build_einsum(eb, vcat(openindices2, openindices)) + return optimizer === nothing ? network : optimize_code(network, optimizer, MergeVectors()) +end + +function build_einsum(eb::EinBuilder, openindices) + return TensorNetwork(EinCode(eb.labels, openindices), eb.tensors) +end + diff --git a/lib/YaoToEinsum/test/circuitmap.jl b/lib/YaoToEinsum/test/circuitmap.jl new file mode 100644 index 00000000..5e6f2c13 --- /dev/null +++ b/lib/YaoToEinsum/test/circuitmap.jl @@ -0,0 +1,339 @@ +module CircuitMapTest +using YaoToEinsum +using Test, OMEinsum +using YaoBlocks, YaoBlocks.YaoArrayRegister +using SymEngine + +const CPhaseGate{T} = ControlBlock{<:ShiftGate{T},<:Any} + +@const_gate ISWAP = PermMatrix([1,3,2,4], [1,1.0im,1.0im,1]) +@const_gate SqrtX = [0.5+0.5im 0.5-0.5im; 0.5-0.5im 0.5+0.5im] +@const_gate SqrtY = [0.5+0.5im -0.5-0.5im; 0.5+0.5im 0.5+0.5im] +# √W is a non-Clifford gate +@const_gate SqrtW = mat(rot((X+Y)/sqrt(2), π/2)) + +""" + singlet_block(θ::Real, ϕ::Real) + +The circuit block for initialzing a singlet state. +""" +singlet_block() = chain(put(2, 1=>chain(X, H)), control(2, -1, 2=>X)) + + +mutable struct FSimGate{T<:Number} <: PrimitiveBlock{2} + theta::T + phi::T +end +YaoBlocks.nqudits(fs::FSimGate) = 2 +YaoBlocks.print_block(io::IO, block::FSimGate) = print(io, "FSim(θ=$(block.theta), ϕ=$(block.phi))") + +function Base.:(==)(fs1::FSimGate, fs2::FSimGate) + return fs1.theta == fs2.theta && fs1.phi == fs2.phi +end + +function YaoBlocks.mat(::Type{T}, fs::FSimGate) where T + θ, ϕ = fs.theta, fs.phi + T[1 0 0 0; + 0 cos(θ) -im*sin(θ) 0; + 0 -im*sin(θ) cos(θ) 0; + 0 0 0 exp(-im*ϕ)] +end + +YaoBlocks.iparams_eltype(::FSimGate{T}) where T = T +YaoBlocks.getiparams(fs::FSimGate{T}) where T = (fs.theta, fs.phi) +function YaoBlocks.setiparams!(fs::FSimGate{T}, θ, ϕ) where T + fs.theta = θ + fs.phi = ϕ + return fs +end + +YaoBlocks.@dumpload_fallback FSimGate FSimGate +YaoBlocks.Optimise.to_basictypes(fs::FSimGate) = fsim_block(fs.theta, fs.phi) + +""" + fsim_block(θ::Real, ϕ::Real) + +The circuit representation of FSim gate. +""" +function fsim_block(θ::Real, ϕ::Real) + if θ ≈ π/2 + return cphase(2,2,1,-ϕ)*SWAP*rot(kron(Z,Z), -π/2)*put(2,1=>phase(-π/4)) + else + return cphase(2,2,1,-ϕ)*rot(SWAP,2*θ)*rot(kron(Z,Z), -θ)*put(2,1=>phase(θ/2)) + end +end + +cphase(nbits::Int, i::Int, j::Int, θ::T) where T = control(nbits, i, j=>shift(θ)) +qft_circuit(::Type{T}, n::Int) where T = chain(n, hcphases(T, n, i) for i = 1:n) +hcphases(::Type{T}, n, i) where T = chain(n, i==j ? put(i=>H) : cphase(n, j, i, 2T(π)/T(2^(j-i+1))) for j in i:n); +qft_circuit(n::Int) = qft_circuit(Float64, n) + +entangler_google53(::Type{T}, nbits::Int, i::Int, j::Int) where T = put(nbits, (i,j)=>FSimGate(T(π)/2, T(π)/6)) + +struct Lattice53 + labels::Matrix{Int} +end + +function Lattice53(;nbits::Int=53) + config = ones(Bool, 5, 12) + config[end,2:2:end] .= false + config[1, 7] = false + labels = zeros(Int, 5, 12) + k = 0 + for (i,c) in enumerate(config) + if c + k += 1 + labels[i] = k + k>=nbits && break + end + end + return Lattice53(labels) +end + +nbits(lattice::Lattice53) = maximum(lattice.labels) + +function Base.getindex(lattice::Lattice53, i, j) + 1<=i<=size(lattice.labels, 1) && 1<=j<=size(lattice.labels, 2) ? lattice.labels[i,j] : 0 +end +upperleft(lattice::Lattice53,i,j) = lattice[i-j%2,j-1] +lowerleft(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j-1] +upperright(lattice::Lattice53,i,j) = lattice[i-j%2,j+1] +lowerright(lattice::Lattice53,i,j) = lattice[i+(j-1)%2,j+1] + +function pattern53(lattice::Lattice53, chr::Char) + res = Tuple{Int,Int}[] + # i0, di, j0, dj and direction + di = 1 + (chr>'D') + dj = 2 - (chr>'D') + j0 = 1 + min(dj-1, mod(chr-'A',2)) + direction = 'C'<=chr<='F' ? lowerright : upperright + for j=j0:dj:12 + i0 = chr>'D' ? mod((chr-'D') + (j-(chr>='G'))÷2, 2) : 1 + for i = i0:di:5 + src = lattice[i, j] + dest = direction(lattice, i, j) + src!=0 && dest !=0 && push!(res, (src, dest)) + end + end + return res +end + +function print_lattice53(lattice, pattern) + for i_=1:10 + i = (i_+1)÷2 + for j=1:12 + if i_%2 == j%2 && lattice[i,j]!=0 + print(" ∘ ") + else + print(" ") + end + end + println() + for j=1:12 + if i_%2 == j%2 && lattice[i,j]!=0 + hasll = (lowerleft(lattice, i, j), lattice[i,j]) in pattern + haslr = (lattice[i,j], lowerright(lattice, i, j)) in pattern + print(hasll ? "/ " : " ") + print(haslr ? " \\" : " ") + else + print(" ") + end + end + println() + end +end + +""" + rand_google53([T=Float64], depth::Int; nbits=53) -> AbstactBlock + +Google supremacy circuit with 53 qubits, also know as the Sycamore quantum supremacy circuits. `T` is the parameter type. + +References +------------------------- +* Arute, Frank, et al. "Quantum supremacy using a programmable superconducting processor." Nature 574.7779 (2019): 505-510. +""" +rand_google53(depth::Int; nbits::Int=53) = rand_google53(Float64, depth; nbits) +function rand_google53(::Type{T}, depth::Int; nbits::Int=53) where T + c = chain(nbits) + lattice = Lattice53(nbits=nbits) + k = 0 + for pattern in Iterators.cycle(['A', 'B', 'C', 'D', 'C', 'D', 'A', 'B']) + push!(c, rand_google53_layer(T, lattice, pattern)) + k += 1 + k>=depth && break + end + return c +end + +function rand_google53_layer(::Type{T}, lattice, pattern) where T + nbit = nbits(lattice) + chain(nbit, chain(nbit, [put(nbit, i=>rand([SqrtW, SqrtX, SqrtY])) for i=1:nbit]), + chain(nbit, [entangler_google53(T, nbit,i,j) for (i,j) in pattern53(lattice, pattern)]) + ) +end + +""" + pair_ring(n::Int) -> Vector + +Pair ring entanglement layout. +""" +pair_ring(n::Int) = [i=>mod(i, n)+1 for i=1:n] + +""" + pair_square(m::Int, n::Int) -> Vector + +Pair square entanglement layout. +""" +function pair_square(m::Int, n::Int; periodic=false) + res = Vector{Pair{Int, Int}}(undef, (m-!periodic)*n+m*(n-!periodic)) + li = LinearIndices((m, n)) + k = 1 + for i = 1:2:m, j=1:n + if periodic || i li[i%m+1, j] + k+=1 + end + end + for i = 2:2:m, j=1:n + if periodic || i li[i%m+1, j] + k+=1 + end + end + for i = 1:m, j=1:2:n + if periodic || j li[i, j%n+1] + k+=1 + end + end + for i = 1:m, j=2:2:n + if periodic || j li[i, j%n+1] + k+=1 + end + end + res +end + +###################### rotor and rotorset ##################### +""" + merged_rotor(noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{1, ComplexF64} + +Single qubit arbitrary rotation unit, set parameters notrailing, noleading true to remove trailing and leading Z gates. + +!!! note + + Here, `merged` means `Rz(η)⋅Rx(θ)⋅Rz(ξ)` are multiplied first, this kind of operation if now allowed in differentiable + circuit with back-propagation (`:BP`) mode (just because we are lazy to implement it!). + But is a welcoming component in quantum differentiation. +""" +merged_rotor(::Type{T}, noleading::Bool=false, notrailing::Bool=false) where T = noleading ? (notrailing ? Rx(zero(T)) : chain(Rx(zero(T)), Rz(zero(T)))) : (notrailing ? chain(Rz(zero(T)), Rx(zero(T))) : chain(Rz(zero(T)), Rx(zero(T)), Rz(zero(T)))) + +""" + rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{nbit, ComplexF64} + +Arbitrary rotation unit (put in `nbit` space), set parameters notrailing, noleading true to remove trailing and leading Z gates. +""" +function rotor(::Type{T}, nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false) where T + rt = chain(nbit, [put(nbit, ibit=>Rz(zero(T))), put(nbit, ibit=>Rx(zero(T))), put(nbit, ibit=>Rz(zero(T)))]) + noleading && popfirst!(rt) + notrailing && pop!(rt) + rt +end + +rotorset(::Type{T}, ::Val{:Merged}, nbit::Int, noleading::Bool=false, notrailing::Bool=false) where T = chain(nbit, [put(nbit, j=>merged_rotor(T, noleading, notrailing)) for j=1:nbit]) +rotorset(::Type{T}, ::Val{:Split}, nbit::Int, noleading::Bool=false, notrailing::Bool=false) where T = chain(nbit, [rotor(T, nbit, j, noleading, notrailing) for j=1:nbit]) +rotorset(::Type{T}, mode::Symbol, nbit::Int, noleading::Bool=false, notrailing::Bool=false) where T = rotorset(T, Val(mode), nbit, noleading, notrailing) + +""" + variational_circuit([T=Float64], nbit[, nlayer][, pairs]; mode=:Split, do_cache=false, entangler=cnot) + +A kind of widely used differentiable quantum circuit, angles in the circuit is randomely initialized. + +### Arguments + +* `T` is the parameter type. +* `pairs` is list of `Pair`s for entanglers in a layer, default to `pair_ring` structure, +* `mode` can be :Split or :Merged, +* `do_cache` decides whether cache the entangler matrix or not, +* `entangler` is a constructor returns a two qubit gate, `f(n,i,j) -> gate`. + The default value is `cnot(n,i,j)`. + +### References + +1. Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Chow, J. M., & Gambetta, J. M. (2017). Hardware-efficient Quantum Optimizer for Small Molecules and Quantum Magnets. Nature Publishing Group, 549(7671), 242–246. https://doi.org/10.1038/nature23879. +""" +function variational_circuit(::Type{T}, nbit, nlayer, pairs; mode=:Split, do_cache=false, entangler=cnot) where T + circuit = chain(nbit) + + ent = chain(nbit, entangler(nbit, i, j) for (i, j) in pairs) + if do_cache + ent = ent |> cache + end + has_param = nparameters(ent) != 0 + for i = 1:(nlayer + 1) + i!=1 && push!(circuit, has_param ? deepcopy(ent) : ent) + push!(circuit, rotorset(T, mode, nbit, i==1, i==nlayer+1)) + end + circuit +end + +variational_circuit(::Type{T}, n::Int; kwargs...) where T = variational_circuit(T, n, 3, pair_ring(n); kwargs...) + +variational_circuit(::Type{T}, nbit::Int, nlayer::Int; kwargs...) where T = variational_circuit(T, nbit, nlayer, pair_ring(nbit); kwargs...) +variational_circuit(nbit::Int; kwargs...) = variational_circuit(Float64, nbit; kwargs...) +variational_circuit(nbit::Int, nlayer::Int; kwargs...) = variational_circuit(Float64, nbit, nlayer; kwargs...) + +@testset "YaoToEinsum.jl" begin + n = 5 + a = rand_unitary(4)[:, 1:2] + a1 = rand_unitary(4)[:, 2] + mb = matblock(OuterProduct(conj.(a), a)) + mb1 = matblock(OuterProduct(conj.(a1), a1)) + for c in [put(n, 2=>Y), put(n, 2=>ConstGate.P0), put(n, (3,2)=>mb1), put(n, (2, 1)=>mb), put(n, 2=>ConstGate.P1), put(n, (5,3)=>SWAP), put(n, (4,2)=>ConstGate.CNOT), put(n, (2,3,1)=>kron(ConstGate.CNOT, X)), + put(n, 2=>Z), control(n, -3, 2=>X), control(n, 3, 2=>X), control(n, (2, -1), 3=>Y), control(n, (4,1,-2), 5=>Z)] + @show c + C = chain([put(n, i=>Rx(rand()*2π)) for i=1:n]..., c) + code, xs = yao2einsum(C; optimizer=nothing) + optcode = optimize_code(code, uniformsize(code, 2), GreedyMethod()) + @test reshape(optcode(xs...; size_info=uniformsize(code, 2)), 1<rand_state(1) for i=1:n]) + reg = join([initial_state[i] for i=n:-1:1]...) + reg |> c + inner = (2,3) + focus!(reg, inner) + for final_state in [Dict([i=>rand_state(1) for i in inner]), Dict([i=>1 for i in inner])] + freg = join(YaoToEinsum.render_single_qubit_state(ComplexF64, final_state[3]), YaoToEinsum.render_single_qubit_state(ComplexF64, final_state[2])) + net = yao2einsum(c; initial_state=initial_state, final_state=final_state, optimizer=TreeSA(nslices=3)) + println(net) + @test vec(contract(net)) ≈ vec(transpose(statevec(freg)) * state(reg)) + end +end + +@testset "symbolic" begin + n = 5 + c = qft_circuit(n) + initial_state = Dict([i=>zero_state(Basic, 1) for i=1:n]) + code, xs = yao2einsum(c; initial_state=initial_state) + @test eltype(xs) == AbstractArray{Basic} +end + +@testset "fix to basic type" begin + c = chain(kron(X,X)) + @test (yao2einsum(c) |> first) isa OMEinsum.SlicedEinsum +end +end \ No newline at end of file diff --git a/lib/YaoToEinsum/test/runtests.jl b/lib/YaoToEinsum/test/runtests.jl new file mode 100644 index 00000000..fd313091 --- /dev/null +++ b/lib/YaoToEinsum/test/runtests.jl @@ -0,0 +1,5 @@ +using YaoToEinsum, Test + +@testset "circuitmap" begin + include("circuitmap.jl") +end \ No newline at end of file diff --git a/src/Yao.jl b/src/Yao.jl index 87d1ac4b..5b29208a 100644 --- a/src/Yao.jl +++ b/src/Yao.jl @@ -15,7 +15,7 @@ Extensible Framework for Quantum Algorithm Design for Humans. const 幺 = Yao using Reexport -@reexport using YaoArrayRegister, YaoBlocks, YaoSym, YaoPlots +@reexport using YaoArrayRegister, YaoBlocks, YaoSym, YaoPlots, YaoToEinsum using YaoArrayRegister.BitBasis, YaoAPI using YaoBlocks: diff --git a/test/runtests.jl b/test/runtests.jl index ba573ab1..5ec31bf9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Yao using YaoAPI using YaoArrayRegister using YaoBlocks +using YaoToEinsum using YaoSym using BitBasis using Documenter @@ -12,11 +13,12 @@ using Random include("easybuild/easybuild.jl") end -DocMeta.setdocmeta!(Yao, :DocTestSetup, :(using Yao, YaoAPI, YaoArrayRegister, YaoBlocks, YaoPlots, YaoSym, BitBasis); recursive=true) +DocMeta.setdocmeta!(Yao, :DocTestSetup, :(using Yao, YaoAPI, YaoArrayRegister, YaoBlocks, YaoPlots, YaoSym, YaoToEinsum, BitBasis); recursive=true) Documenter.doctest(YaoAPI; manual=false) Documenter.doctest(YaoArrayRegister; manual=false) Documenter.doctest(YaoBlocks; manual=false) Documenter.doctest(YaoSym; manual=false) Documenter.doctest(YaoPlots; manual=false) +Documenter.doctest(YaoToEinsum; manual=false) Documenter.doctest(Yao; manual=false)