PLEASE NOTE THIS IS PRE-RELEASE SOFTWARE

PastaQ.jl: design and benchmarking quantum hardware

PastaQ.jl is a Julia software toolbox providing a range of computational methods for quantum computing applications. Some examples are the simulation of quancum circuits, the design of quantum gates, noise characterization and performance benchmarking. PastaQ relies on tensor-network representations of quantum states and processes, and borrows well-refined techniques from the field of machine learning and data science, such as probabilistic modeling and automatic differentiation.

---

Install

The PastaQ package can be installed with the Julia package manager. From the Julia REPL, type ] to enter the Pkg REPL mode and run:

julia> ]

pkg> add PastaQ

PastaQ.jl relies on the following packages: ITensors.jl for low-level tensor-network algorithms, Optimisers.jl for stochastic optimization methods, Zygote.jl for automatic differentiation, and Observers.jl for tracking/recording metrics and observables. Please note that right now, PastaQ.jl requires that you use Julia v1.6 or later.

---

Documentation

• STABLE -- documentation of the most recently tagged version.
• DEVEL -- documentation of the in-development version.

Examples

We briefly showcase some of the functionalities provided by PastaQ.jl. For more in-depth discussion, please refer to the tutorials folder.

Simulating quantum circuits

The vast majority of tasks related to designing and benchmarking quantum computers relies on the capability of simulating quantum circuits, built out of a collection of elementary quantum gates. In PastaQ, a quantum gates is described by a data structure g = ("gatename", support, params), consisting of a gate identifier gatename (a String), a support (an Int for single-qubit gates or a Tuple for multi-qubit gates), and a set of gate parameters, such as rotations angles, whenever needed. A comprehensive set of elementary gates is provided, including Pauli operations, phase and T gates, single-qubit rotations, controlled gates, Toffoli gate and others. Additional user-specific gates can be easily added to this collection. Once a circuit is defined, it can be executed using the runcircuit function:

using PastaQ

# a quantum circuit in PastaQ
gates = [("X" , 1),                              # Pauli X on qubit 1
         ("CX", (1, 3)),                         # Controlled-X on qubits [1,3]
         ("Rx", 2, (θ = 0.5,)),                  # Rotation of θ around X
         ("Rn", 3, (θ = 0.5, ϕ = 0.2, λ = 1.2)), # Arbitrary rotation with angles (θ,ϕ,λ)
         ("√SWAP", (3, 4)),                      # Sqrt Swap on qubits [2,3]
         ("T" , 4)]                              # T gate on qubit 4
         
# run the circuit
ψ = runcircuit(gates)
# returns the MPS at the output of the quantum circuit: `|ψ⟩ = Û|0,0,…,0⟩`
# first the gate ("X" , 1) is applied, then ("CX", (1, 3)), etc.

# ------------------------------------------------------------------
# Output:
#  ITensors.MPS
#  [1] ((dim=2|id=456|"Qubit,Site,n=1"), (dim=1|id=760|"Link,n=1"))
#  [2] ((dim=1|id=760|"Link,n=1"), (dim=2|id=613|"Qubit,Site,n=2"), (dim=1|id=362|"Link,n=1"))
#  [3] ((dim=2|id=9|"Qubit,Site,n=3"), (dim=2|id=357|"Link,n=1"), (dim=1|id=362|"Link,n=1"))
#  [4] ((dim=2|id=980|"Qubit,Site,n=4"), (dim=2|id=357|"Link,n=1"))

In this next example, we create a circuit to prepare the GHZ state, and sample projective measurements in the computational basis. We then execture the circuit in the presence of noise, where a local noise channel is applied to each gate. A noise model is described as noisemodel = ("noisename", (noiseparams...)), in which case it is applied to each gate identically. To distinguish between one- and two-qubit gates, for example, the following syntax can be used: noisemodel = (1 => noise1, 2 => noise2). For more sophisticated noise models (such as gate-dependent noise), please refer to the documentation.

using PastaQ
using ITensors

# number of qubits
n = 20

# manually create a circuit to prepare GHZ state,
# or use built-in call `circuit = ghz(n)` 
circuit = Tuple[("H", 1)]
for j in 1:n-1
  push!(circuit, ("CX", (j, j+1)))
end

# run the circuit to obtain the output MPS
hilbert = qubits(n)
ψ = runcircuit(hilbert, circuit)


# sample projective measurements in the computational basis
@show getsamples(ψ, 5)

# define a noise model with different error rates for
# one- and two-qubit gates
noisemodel = (1 => ("depolarizing", (p = 0.01,)), 
              2 => ("depolarizing", (p = 0.05,)))

# run a noisy circuit
ρ = runcircuit(hilbert, circuit; noise = noisemodel)
@show fidelity(ψ, ρ)
@show getsamples(ρ, 5)

# quantum processes can also be obtained.
# unitary MPO
U = runcircuit(circuit; process = true)
# Choi matrix
Λ = runcircuit(circuit; process = true, noise = noisemodel)

# ------------------------------------------------------------------
# Output:
#  5×20 Matrix{Int64}:
#  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
#  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
#  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
#  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
#  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
# 
#  fidelity(ψ, ρ) = 0.40840853095498975
# 
#  5×20 Matrix{Int64}:
#  1  1  1  1  1  1  0  0  0  0  0  0  1  0  0  0  1  1  1  1
#  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
#  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1
#  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
#  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1

There is a number of built-in circuits available, one examples being random circuits. In the following, we generate a one-dimensional random quantum circuits built out with a brick-layer geometry of alternative CX gates and layers of random single-qubit rotations:

n = 100 
depth = 20
circuit = randomcircuit(n; depth = depth, 
                           twoqubitgates = "CX", 
                           onequbitgates = "Rn")
@time ψ = runcircuit(circuit; cutoff = 1e-10)
@show maxlinkdim(ψ)
# ------------------------------------------------------------------
# Output:
#  89.375383 seconds (5.25 M allocations: 64.781 GiB, 9.98% gc time)
#  maxlinkdim(ψ) = 908

Variational quantum eingensolver

We show how to perform a ground state search of a many-body hamiltonian $H$ using the variational quantum eigensolver (VQE). The VQE algorithm, based on the variational principle, consists of an iterative optimization of an objective function $\langle \psi(\theta)|H|\psi(\theta)\rangle/\langle\psi(\theta)|\psi(\theta)\rangle$, where $|\psi(\theta)\rangle = U(\theta)|0\rangle$ is the output wavefunction of a parametrized quantum circuit $U(\theta)$.

In the following example, we consider a quantum Ising model with 10 spins, and perform the optimization by leveraging Automatic Differentiation techniques (AD), provided by the package Zygote.jl. Specifically, we build a variational circuit using built-in circuit-contruction functions, and optimize the expectation value of the Hamiltonian using a gradient-based approach and the LBFGS optimizer. The gradients are evaluated through AD, providing a flexible interface in defining custom variational circuit ansatze.

using PastaQ
using ITensors
using Random
using Printf
using OptimKit
using Zygote

N = 10   # number of qubits
J = 1.0  # Ising exchange interaction
h = 0.5  # transverse magnetic field

# Hilbert space
hilbert = qubits(N)

function ising_hamiltonian(N; J, h)
  os = OpSum()
  for j in 1:(N - 1)
    os += -J, "Z", j, "Z", j + 1
  end
  for j in 1:N
    os += -h, "X", j
  end
  return os
end

# define the Hamiltonian
os = ising_hamiltonian(N; J, h)

# build MPO "cost function"
H = MPO(os, hilbert)
# find ground state with DMRG

nsweeps = 10
maxdim = [10, 20, 30, 50, 100]
cutoff = 1e-10
Edmrg, Φ = dmrg(H, randomMPS(hilbert); outputlevel=0, nsweeps, maxdim, cutoff);
@printf("\nGround state energy from DMRG: %.10f\n\n", Edmrg)

# layer of single-qubit Ry gates
Rylayer(N, θ) = [("Ry", j, (θ=θ[j],)) for j in 1:N]

# brick-layer of CX gates
function CXlayer(N, Π)
  start = isodd(Π) ? 1 : 2
  return [("CX", (j, j + 1)) for j in start:2:(N - 1)]
end

# variational ansatz
function variationalcircuit(N, depth, θ⃗)
  circuit = Tuple[]
  for d in 1:depth
    circuit = vcat(circuit, CXlayer(N, d))
    circuit = vcat(circuit, Rylayer(N, θ⃗[d]))
  end
  return circuit
end

depth = 20
ψ = productstate(hilbert)

cutoff = 1e-8
maxdim = 50

# cost function
function loss(θ⃗)
  circuit = variationalcircuit(N, depth, θ⃗)
  Uψ = runcircuit(ψ, circuit; cutoff, maxdim)
  return inner(Uψ', H, Uψ; cutoff, maxdim)
end

Random.seed!(1234)

# initialize parameters
θ⃗₀ = [2π .* rand(N) for _ in 1:depth]

# run VQE using BFGS optimization
optimizer = LBFGS(; maxiter=50, verbosity=2)
function loss_and_grad(x)
  y, (∇,) = withgradient(loss, x)
  return y, ∇
end
θ⃗, fs, gs, niter, n