# Hybrid discrete-continuous compilation of trapped-ion quantum circuits with
deep reinforcement learning

Francesco Preti  Forschungszentrum Jülich, Institute of Quantum Control
(PGI-8), D-52425 Jülich, Germany  Institute for Theoretical Physics,
University of Cologne, D-50937 Köln, Germany  Michael Schilling
Forschungszentrum Jülich, Institute of Quantum Control (PGI-8), D-52425
Jülich, Germany  Institute for Theoretical Physics, University of Cologne,
D-50937 Köln, Germany  Sofiene Jerbi  Institute for Theoretical Physics,
University of Innsbruck, 6020 Innsbruck, Austria  Dahlem Center for Complex
Quantum Systems, Freie Universität Berlin, Berlin, Germany.  Lea M.
Trenkwalder  Institute for Theoretical Physics, University of Innsbruck, 6020
Innsbruck, Austria  Hendrik Poulsen Nautrup  Institute for Theoretical
Physics, University of Innsbruck, 6020 Innsbruck, Austria  Felix Motzoi
Forschungszentrum Jülich, Institute of Quantum Control (PGI-8), D-52425
Jülich, Germany  Hans J. Briegel  Institute for Theoretical Physics,
University of Innsbruck, 6020 Innsbruck, Austria

(August 30, 2024)

###### Abstract

Shortening quantum circuits is crucial to reducing the destructive effect of
environmental decoherence and enabling useful algorithms. Here, we demonstrate
an improvement in such compilation tasks via a combination of using hybrid
discrete-continuous optimization across a continuous gate set, and
architecture-tailored implementation. The continuous parameters are discovered
with a gradient-based optimization algorithm, while in tandem the optimal gate
orderings are learned via a deep reinforcement learning algorithm, based on
projective simulation. To test this approach, we introduce a framework to
simulate collective gates in trapped-ion systems efficiently on a classical
device. The algorithm proves able to significantly reduce the size of relevant
quantum circuits for trapped-ion computing. Furthermore, we show that our
framework can also be applied to an experimental setup whose goal is to
reproduce an unknown unitary process.

##  I Introduction

The last decade has seen significant progress in the development of quantum
computing architectures [1]. While scalable fault-tolerant quantum computers
are still out of reach in the near future, noisy, intermediate-scale quantum
(NISQ) computers may already offer some benefits over classical ones for
specific computational tasks [2, 3]. In particular, variational algorithms [4,
5], where most of the operations depend on several continuous parameters, have
emerged as a suitable class of methods that could potentially achieve quantum
speed-up on NISQ devices.

Implementing a high-level quantum algorithm both on fault-tolerant and NISQ
devices requires adequate methods to compile it in the set of universal
quantum gates available to the hardware. While several frameworks for
compilation of quantum circuit-based algorithms on physical platforms are
being developed [6], common available approaches, such as heuristic and
automated search [7, 8], in many case cannot output an optimal circuit for a
specific target operation [9].

In digital quantum computers, compilers implement a general quantum circuit
through a discrete set of universal quantum gates. In real physical quantum
devices, however, and more generally in analog computation, an additional
layer of complexity is present, due to the necessity of optimizing continuous
gate parameters to reproduce target unitaries. These parameters may depend,
e.g., on the specific Hamiltonian employed in the quantum computing platform,
such as the XY-Hamiltonian or the Mølmer–Sørensen interaction for trapped ions
[10, 11], the cross resonance interaction for IBM quantum computers [12] or
the Fermi-Hubbard model for neutral atoms [13]. As a result, when taking into
account physical parameters, variational algorithms, and generally continuous
gate sets [14],  one must supplement the circuit compilation task with a
subsequent optimization of the parameters defining the individual constituent
gates. For the optimization of continuous parameters, we have several options
available, e.g., gradient-based algorithms [15], evolutionary algorithms [16]
and direct search [17]. For the compilation of the circuit gate structure, a
standard approach is given by methods based on the Solovay-Kitaev algorithm
[18], different circuit factorization strategies [19, 20], graph path
traversal algorithms, such as the $A^{*}$ algorithm [21], semi-definite
programming and and various machine learning methods, including reinforcement
learning [9]. More specifically, deep reinforcement learning has been recently
successfully implemented for the optimization of discrete quantum circuits
[22, 23, 24, 25, 26, 27, 28].

In reinforcement learning (RL) [29], an agent learns to maximize a properly
engineered reward signal by interacting with an environment, which encodes the
optimization task to solve. RL has already been applied to solve various
challenging tasks, including, e.g., surpassing human performance in certain
classes of games [30] or in complex, computationally expensive problems such
as protein folding [31]. Projective Simulation (PS) is a physics-inspired
framework for intelligent agents which has also been applied to solve RL tasks
in quantum physics [32, 33, 34, 35, 36, 37] and biology [38]. This model can
naturally be extended to a deep RL model [39, 40] and has found applications
in representation learning [41, 42].

In this work, we propose a unified approach to optimizing the placement and
parameter optimization of the gates in the circuit. We argue that such an
approach is both relevant to traditional compilation of a more expressive,
continuous gateset [14], as well as to variational circuits where the
experimental cost-function optimization may be simultaneously performed over
discrete and continuous degrees of freedom. We use a RL agent, based on the PS
framework, for the combinatorial optimization and a gradient-based optimizer
for the continuous optimization of the gate angles. In particular, we extend
the method proposed in [24] for the optimization of variational circuits in
quantum chemistry and [25] for state preparation to the case of unitary
compilation. We consider the framework where an agent, that has control over
an experimental platform, interacts with a black-box unitary process and
attempts to simulate it. The task of the RL agent is to optimize the position
of the gates on the circuit, whereas the gradient-based optimizer finds the
optimal set of continuous parameters that minimize a given cost function. The
reward function for the RL agent is constructed based on the results of the
continuous optimization.

We test the learning algorithm first on standard unitaries, such as Toffoli
gates, and then consider the task of quantum process simulation. The latter
can be conceived as an experimental black-box unitary approximation strategy
that allows the agent to reconstruct an unknown unitary process by compiling a
proper quantum circuit.

Independently of our circuit compilation results, we propose a method to
speed-up the simulation of quantum circuits based on an efficient
representation of trapped-ion gates replacing standard matrix exponentiation.
This method allows us to obtain analytic expressions for the ion-trap gates
and their gradients and also to simulate the given gate set on other quantum
computing platforms.

The paper is organized as follows: In Section II we introduce the quantum
circuit framework for trapped-ions. In Section III.1 we present our method to
compute fast analytic ion gates in simulation. In Section III.2 we discuss
continuous optimization methods and strategies to compute the gradients of the
cost function both in simulation and on a real quantum device. In Section
III.3 we introduce PS and its extensions in the context of RL methods. In
Section IV we introduce the problem of circuit synthesis and our hybrid RL-
continuous optimization method. In Section V we discuss the results of
applying the proposed method to the compilation of (black-box) unitaries.

Figure 1: General scheme of the proposed hybrid training loop, where the RL
agent (a) learns to optimize a (variational or traditional) quantum circuit.
The agent is equipped with a memory structure (b) that allows it to process
the information input by the environment (in the case of PS, the episodic and
compositional memory, ECM [32, 33]). Every perceptual input from the
environment triggers a sequence of transitions – showed in orange in (b) –
within the internal computational graph of the agent and modelled by
transition probabilities $p_{ij}$. The sequence of transitions connects a
percept clip $s_{t}$, which corresponds to the sequence of actions
$a_{1},...,a_{t-1}$ used to generate the circuit $V_{t}$, to an action clip
corresponding to $a_{t}$. By choosing an index $a_{t}$ for the action gates
$G_{1}=\text{MS},G_{2}=C_{xy},G_{3}=Z_{1},...,G_{n+2}=Z_{n}$, which
corresponds to the gate set introduced in Eqs. (1)-(3), the agent places a new
gate on the circuit (c). At each RL time step $t$, the circuit is used to
compute a cost function $C(\boldsymbol{\alpha})$ based on the fidelity – see
Eq. (4) – in simulation, or the Hilbert-Schmidt test – see Eq. (13) and (d) –
in experiments. The continuous optimizer (e) then outputs a guess for the
optimal parameters $\boldsymbol{\alpha}^{*}$ that correspond to the cost
function minimum for a specific circuit $V_{t}$. The minimal value of the cost
function is used to assign a reward (f) to the PS agent, thus closing the RL
training loop.

##  II Problem setting

In this work, we consider a specific set of gates, normally implemented in
trapped-ion quantum circuits, which is based on global Mølmer–Sørensen (MS)
gates, equatorial rotations acting on the entire register of qubits, as well
as local polar rotations. Trapped ions are among the most promising platforms
for quantum computing hardware [43]. They exhibit impressive coherence times
even in absence of dynamical decoupling and spin echo techniques [44] and have
been shown to allow for high-fidelity quantum gates [45, 46]. In a trapped-ion
quantum computer, ions – usually ^43Ca^+ – are confined in a Paul trap [47]
using a varying electromagnetic field. The ions are addressed individually by
a system of lasers aligned externally to the trap. The interaction of the
laser field with the ion motional and electronic degrees of freedom allows for
entangling operations. The laser pulses can be engineered to define the
following universal set of quantum gates [48, 49, 50]:

$\displaystyle\text{MS}(\theta,\phi)=\exp{-i\frac{\theta}{4}(S_{x}\cos\phi+S_{y}\sin\phi)^{2}}$
(1) $\displaystyle
C_{xy}(\theta,\phi)=\exp{-i\frac{\theta}{2}(S_{x}\cos\phi+S_{y}\sin\phi)}$ (2)
$\displaystyle Z_{j}(\theta)=\exp{-i\frac{\theta}{2}\sigma_{z}^{(j)}},$ (3)

where the first gate is a global MS gate [10, 11], the second gate is a
rotation of the entire register of qubits in the equatorial plane of the Bloch
sphere, and the third gate represents a single-qubit, and therefore local,
$\sigma_{z}$ rotation acting on qubit $j$. The operators
$S_{x}=\sum_{i=1}^{n}\sigma_{x}^{(i)}$,
$S_{y}=\sum_{i=1}^{n}\sigma_{y}^{(i)}$, $S_{z}=\sum_{i=1}^{n}\sigma_{z}^{(i)}$
are given by the Pauli operators acting on qubit $i$. These gates can be
easily generalized to the qudit case [51]. For the optimization of unitaries,
the gate overlap fidelity is a standard figure of merit [52]:

$\displaystyle F(U,V)=\frac{1}{d^{2}}|\tr{VU^{\dagger}}|^{2},$ (4)

where $U$ and $V$ are unitaries (one of them is the goal of the optimization)
and $d$ is the dimension of the Hilbert space – for $n$-qubit systems
$d=2^{n}$. Commonly, synthesising quantum circuits to reproduce an arbitrary
unitary $U$ requires having access to quantum computing hardware that
implements a universal gate set and running the optimization of the continuous
parameters.

If the circuit synthesis is performed on a classical computer, it is also
necessary to simulate the gate set efficiently. Simulating and optimizing
$n$-qubit collective gates such as those given in Eq. (1) and Eq. (2) is
generally considered more challenging than with just two-qubit entangling
gates and single qubit rotations [48]. A standard strategy is to progressively
increase the number of entangling MS gates on the circuit, accompanied by a
suitable number of single and multi-qubit rotations, and to progressively
optimize the gate parameters until an acceptable threshold of the figure of
merit is reached. This is a viable strategy to obtain one solution for a
quantum compilation problem on a trapped-ion device; it is however sub-optimal
with respect to the number of gates required. More efficient solutions exist,
but the optimization landscape may be difficult to navigate for various
algorithms [53]. In particular, as the optimization landscape with respect to
the gate angles is particularly vast and depends on the arrangement of the
gates on the circuit, it is often necessary to search through several
combinations of gate sequences. Random or automated search can be implemented
to reduce the number of gates present on the circuit by arbitrarily trying
several configurations [22, 48].

##  III Methods

In the following section, we discuss our approach to address the three
relevant aspects that characterize a circuit synthesis task: the efficient
computer-aided simulation of the relevant trapped-ion gates, the optimization
of the continuous gate parameters and the optimal arrangement of the gates on
the circuit. In addition, we address the possible implementation of our hybrid
RL-gradient based optimization method on real quantum devices.

###  III.1 Dynamics via fast exponentiation

In simulation, direct computation of the gates in Eqs. (1)-(2) is commonly
done via direct matrix exponentiation of a Hamiltonian, which is generally
slow for large matrices, since it requires several matrix multiplications,
each one with a complexity of at least $O(d^{2.8})$ [54], where $d$ is the
matrix dimension. Instead, here we show how to reduce the per-iteration cost
to quadratic by factorizing the exponential with respect to the rotation angle
$\theta$ and to the phase angle $\phi$. This is achieved via a spectral
decomposition where the phase-independent part can then be cached before the
optimization.

Mathematically, the factorization of the matrix exponential for a MS or
$C_{xy}$ gate can be written via spectral decomposition as

$\displaystyle U_{H}(\theta,\phi)=\exp{iH(\theta,\phi)}=$ (5)
$\displaystyle=V(\theta,\phi)\left(\sum_{l=0}^{d-1}e^{-i\lambda_{l}(\theta,\phi)}\outerproduct{l}{l}\right)V^{\dagger}(\theta,\phi),$

where $H(\theta,\phi)$ is a $\theta$\- and $\phi$-dependent Hamiltonian – see
exponents in Eqs. (1)-(2) –, $V$ is the matrix of eigenvectors with respective
eigenvalues $\lambda_{l}$.

Regrouping in terms of degenerate eigenvalues, and considering the case where
the set of eigenvalues are $\phi$-independent, i.e.,
$\lambda_{l}=\lambda_{l}(\theta)$, while a $\phi$-dependent unitary is applied
to the Hamiltonian, i.e., $H(\theta,\phi)=V(\phi)H(\theta)V^{\dagger}(\phi)$.
As a consequence, the eigenvectors matrix $V=V(\phi)$ is independent of
$\theta$, and we can rewrite

$\displaystyle U_{H}(\theta,\phi)=$
$\displaystyle\sum_{k=0}^{n_{\lambda}}V(\phi)\left(e^{-i\lambda_{k}(\theta)}\sum_{l_{k}}\outerproduct{l_{k}}{l_{k}}\right)V^{\dagger}(\phi)$
(6) $\displaystyle=$ $\displaystyle
P(\phi)\odot\sum_{k=0}^{n_{\lambda}}e^{-i\lambda_{k}(\theta)}D_{l},$ (7)

where $n_{\lambda}$ is the number of distinct eigenvalues, $\ket{l_{k}}$ are
the respective eigenvectors,
$D_{k}=V(0)\sum_{l_{k}}\outerproduct{l_{k}}{l_{k}}V^{\dagger}(0)$ and
$P(\phi)$ is a matrix of phase components. In our case, the columns of
$P(\phi)$ are all equal and are given by the vector
$\boldsymbol{p}(\phi)=\begin{pmatrix}1\\\ e^{i\phi}\end{pmatrix}^{\otimes n}$,
while $\odot$ represents element-wise (Hadamard) matrix multiplication. The
$C_{xy}$ and MS gates have few unique eigenvalues with
$\lambda_{k}=(2k-n)\theta,\ 0\leq k\leq n_{\lambda}=n$ and
$\lambda_{k}=(2k-n)^{2}\theta,\ 0\leq k\leq n_{\lambda}=\lceil
n+\frac{1}{2}\rceil$, respectively, making the computation with cached $D_{k}$
particularly efficient. A detailed derivation, together with a discussion of
the computational speedup of this representation, is available in Appendix A.

###  III.2 Continuous gradient-based optimization

We consider the optimization of an observable with respect to the continuous
parameters. Although single problems can be optimized effectively by
implementing an appropriate discretization method, in general this approach
can lead to a sub-optimal solution, since it introduces discontinuities in the
search space. Instead, we want to consider the dependency of the fidelity from
continuous parameters and calculate its gradient.

####  III.2.1 Cost function based on known target unitaries.

We consider a quantum circuit composed of a sequence of continuous gates with
angle parameters
$\boldsymbol{\alpha}=(\boldsymbol{\alpha}_{1},...,\boldsymbol{\alpha}_{L})$
and where each gate is composed of multiple rotation angles
$\boldsymbol{\alpha}_{1}\in\mathbb{R}^{M_{1}},\ldots,\boldsymbol{\alpha}_{L}\in\mathbb{R}^{M_{L}}$
and $M_{1},...,M_{L}\geq 1$ and $M=\sum_{m=1}^{L}M_{m}$. The circuit is then
given by

$\displaystyle
V(\boldsymbol{\alpha})=\prod_{l=1}^{L}V_{l}(\boldsymbol{\alpha}_{l}),$ (8)

where $V_{l}(\boldsymbol{\alpha}_{l})$ is an arbitrary parametric unitary with
parameters $\boldsymbol{\alpha}_{l}$.

The gradient with respect to a figure of merit, here the average gate fidelity
– see Eq. (4) –, can also be directly computed and are given by

$\displaystyle\nabla_{\boldsymbol{\alpha}_{l}}F=\frac{2}{d^{2}}\Re{\tr(\nabla_{\boldsymbol{\alpha}_{l}}V(\boldsymbol{\alpha})U^{\dagger})\tr(V(\boldsymbol{\alpha})^{\dagger}U)},$
(9)

with

$\displaystyle\nabla_{\boldsymbol{\alpha}_{l}}V(\boldsymbol{\alpha})=V_{1}(\boldsymbol{\alpha}_{l1})\cdots\nabla_{\boldsymbol{\alpha}_{l}}V_{l}(\boldsymbol{\alpha}_{l})\cdots
V_{L}(\boldsymbol{\alpha}_{L}),$ (10)

for $1\leq l\leq L$. The naive element-wise computation of the gradient
components uses $M\cdot L$ matrix multiplications, as the matrix in Eq. (9)
needs to be evaluated $M$ times and is given by the product of $L$ unitaries.
However, the gradient can be computed recursively, a method which is often
referred to as GRAPE [55], by storing the values of the intermediate unitaries
$W_{l}=W_{l-1}V_{l}^{\dagger}(\boldsymbol{\alpha}_{l}),W_{0}=I_{d}$ for
$l=1,...,L$ in the product

$\displaystyle\nabla_{\boldsymbol{\alpha}_{l}}F=\frac{2}{d^{2}}\Re{\tr(W_{l-1}\nabla_{\boldsymbol{\alpha}_{l}}V_{l}(\boldsymbol{\alpha}_{l})W_{l}^{\dagger}VU^{\dagger})\tr(V^{\dagger}U)}.$
(11)

The gradient computation method given in Eq. (11) computes first the
intermediate unitaries, for which we need $L$ matrix multiplications and uses
them to evaluate the gradient, which compared to Eq. (9) needs only $4M+L+1$
matrix multiplications and therefore scales linearly with the number of gates
and gate parameters in the circuit.

Eq. (11) allows us to compute the gradient of any cost function that resembles
the structure of Eq. (4) efficiently in numerical simulations. This can be
further sped up by similarly analytically computing the Hessian matrix of the
cost function [56] or by using GPU or TPU architectures [57]. Eq. (11) can
also be useful in an experimental setting when the target dynamics, such as a
desired state or unitary evolution, are known a priori.

####  III.2.2 Cost function based on black-box access to a target unitary.

In many NISQ-relevant algorithms, the structure of the quantum circuit need
not be prescriptive. Instead, we optimize a general circuit ansatz
variationally to minimize some performance cost function. In this setting, we
once again may optimize both the discrete and continuous degrees of freedom of
the circuit to minimize the cost function. One can conceive of situations
where the target circuit unitary $U$ may be given as a black-box process [58,
59], that is where any decomposition of the circuit is unknown and we do not
have a classical description of the unitary entries, or a circuit to compute
them. In this situation, the cost function computing the distance between the
black-box target $U$ and a given parametric quantum circuit
$V(\boldsymbol{\alpha})$ needs to be estimated.

A possible way to compute distances between a unitary implemented on a quantum
computer and a black-box unitary implemented by a different quantum system is
the Hilbert-Schmidt test (HST) [60], a generalization of the SWAP test for
state preparation, which evaluates the gate fidelity of Eq. (4) on a quantum
circuit. A related cost function for black-box quantum compilation is given by
the local Hilbert-Schmidt (LHS) test – see Fig. 1, panel (d) –, which has been
shown to be easier to optimize and less sensitive to Barren plateaus [61, 62].
Observe that the cost function $C_{\text{LHS}}(V,U)$ is bound from above and
below by the average gate fidelity given in Eq. (4),

$\displaystyle C_{\text{LHS}}(V,U)<C_{\text{HST}}(V,U)<nC_{\text{LHS}}(V,U),$
(12)

which implies that minimizing $C_{\text{LHS}}(V,U)$ also minimizes
$C_{\text{HST}}(V,U)$. Other cost functions have been proposed which are based
on so-called incoherent learning, i.e., where the quantum computer and the
black-box quantum system do not need to interact coherently [63].

For gradient-based optimization, we need to differentiate the cost functions
$C_{\text{HST}}$ or $C_{\text{LHS}}$ with respect to continuous gate
parameters. Unfortunately, the gradient method given in Eq. (11) requires
access to the intermediate unitaries, which are generally not available in
real-world scenarios. Moreover, we do not generally have direct access to the
gradient of the unitary operations. As a consequence, we need to use the cost
function itself to estimate its gradient with respect to the continuous
parameters, which proves slower. In fact, computing the gradient in Eq. (9)
with the method of finite differences requires $2M$ evaluations of the test
circuit, where $M$ is the total number of parameters. An alternative is
represented by parameter-dependent sampling strategies that estimate the
gradient via trigonometric interpolation [64].

Let us now consider the specific case of trapped ions with the gate set
introduced in Section II. There are three different types of gates ($Z$, MS
and $C_{xy}$), with four types of parameters shifts: the three rotation angles
of the $Z$ gate, the MS gate and the $C_{xy}$ gate, respectively, and the
phase term $\phi$ of the MS and $C_{xy}$ gates, which is the same for both
unitaries (see also Appendix B).

In the following, we consider a cost function on the Hilbert Schmidt test, but
the results can be applied to the local test as well. In the case of the $Z$
gate, only one parameter $\theta$ is present, whereas for the other two gates
we have two possible parameter-shifts. The experimental quantum cost function
$C_{\text{HST}}$ comparing the parameterized circuit $V(\boldsymbol{\alpha})$
in Eq. (8) with a target unitary $U$ is given by

$\displaystyle C_{HST}(\boldsymbol{\alpha})=1-\Tr{H(U\otimes
V(\boldsymbol{\alpha})^{*})\rho(U^{\dagger}\otimes
V(\boldsymbol{\alpha})^{T})},$ (13)

where $\rho$ and $H$ are projectors defined in Appendix B and Ref. [60]. We
consider here only the shift with respect to one gate at a time, which we name
$V_{1}(\boldsymbol{\alpha}_{1})$, while the remaining gates in the circuit,
i.e., $V_{2}(\boldsymbol{\alpha}_{2}),...,V_{L}(\boldsymbol{\alpha}_{L})$, are
fixed. This procedure can be repeated for each gate parameter. We first
consider the case where $V_{1}(\boldsymbol{\alpha}_{1}=\theta)=Z(\theta)$. The
exact derivative of $C_{\text{HST}}(\theta)$ is given by:

$\displaystyle\partialderivative{\theta}C_{\text{HST}}(\theta,Z)=C_{\text{HST}}(\theta+\frac{\pi}{2},Z)-C_{\text{HST}}(\theta-\frac{\pi}{2},Z)$
(14)

For $V_{1}(\boldsymbol{\alpha}_{1}=(\theta,\phi)^{T})=C_{xy}(\theta,\phi)$ and
following Ref. [64], the parameter derivative with respect to $\theta$ is
given by

$\displaystyle\left.\partialderivative{\theta}C_{\text{HST}}(\theta,\phi,C_{xy})\right|_{\theta=0}=$
(15)
$\displaystyle\sum_{l=1}^{2n}\frac{(-1)^{l-1}}{2\sin\left(\frac{2l-1}{2n}\pi\right)}C_{\text{HST}}\left(\theta+\frac{2l-1}{2n}\pi,\phi,C_{xy}\right),$

whereas for
$V_{1}(\boldsymbol{\alpha}_{1}=(\theta,\phi)^{T})=\text{MS}(\theta,\phi)$, the
parameter-shift rule with respect to $\theta$ becomes,

$\displaystyle\left.\partialderivative{\theta}C_{\text{HST}}(\theta,\phi,\text{MS})\right|_{\theta=0}=$
(16)
$\displaystyle\sum_{l=1}^{2\lceil\frac{n}{2}\rceil}\frac{(-1)^{l-1}}{2\sin\left(\frac{2l-1}{\lfloor\frac{n}{2}+1\rfloor}\pi\right)}C_{\text{HST}}\left(\theta+\frac{2l-1}{\lfloor\frac{n}{2}+1\rfloor}\pi,\phi,\text{MS}\right)\delta_{i}^{\text{floor}}.$

A more simplified rule can be derived for the parameter $\phi$, which is valid
for both the $C_{xy}$ and the MS gate

$\displaystyle\partialderivative{\phi}C_{\text{\tiny
HST}}(\theta,\phi)=C_{\text{\tiny
HST}}(\theta,\phi+\frac{\pi}{4})-C_{\text{\tiny
HST}}(\theta,\phi-\frac{\pi}{4})+$ (17)
$\displaystyle+\Big{(}\frac{1}{\sqrt{2}}-\frac{1}{2}\Big{)}C_{\text{\tiny
HST}}(\theta,\phi-\frac{\pi}{2})+\Big{(}\frac{1}{\sqrt{2}}+\frac{1}{2}\Big{)}C_{\text{\tiny
HST}}(\theta,\phi+\frac{\pi}{2}).$

Using the expressions given by Eq. (14) for the $Z$ gate, by Eq. (15) and Eq.
(17) for the $C_{xy}$ gate and Eq. (16) and Eq. (17) for the MS gate, we can
compute any gradient of cost functions estimated in experiments on quantum
devices that use cost functions such as Eq. (13). For more details about the
parameter-shift rules of the $C_{xy}$ and MS gates, and further possible
simplifications, see Appendix B.

###  III.3 Combinatorial Optimization

In this section, we consider the problem of determining an optimal arrangement
of the gates on the quantum circuit. This problem arises from the necessity of
minimizing the depth of a quantum circuit to reduce the total circuit
execution time, thereby reducing decoherence. And while the circuit can be
compiled in layers [65], it is difficult to determine the optimal size due to
the large number of possible optimal parameter configurations that produce the
same circuit. Therefore, it is necessary to search through various
combinations of gate arrangements to determine a minimal one. This task, which
partially falls in the realm of combinatorial optimization [66], is
particularly suitable for reinforcement learning algorithms [22].

####  III.3.1 Reinforcement Learning with Projective Simulation

Reinforcement learning describes a class of algorithms which use an agent-
environment interaction model to maximize a reward function. In particular,
the agent can be represented by a parametrized model, called policy, that
outputs action signals on the environment upon receiving observations as
inputs. For each action or sequence thereof, the environment returns
observations and reward signals. In gradient-based methods, the policy
parameters are updated in the direction that maximizes the discounted future
expected return [29]. Based upon the different tasks considered, a vast range
of algorithms and methods have been developed to tackle various environments
[67].

In the following section, we consider the PS architecture.

Figure 2: This figure represents the policy parametrization of the PS-LSTM
algorithm. In this implementation, the agent receives a sequential input at
time $t$ and constructs a policy by outputting weights $h_{t}(s,a,u_{t})$
corresponding to each action $a$, the current percept $s$ and the hidden state
of the LSTM network $u_{t}$. Afterwards, an action $a^{*}$ is sampled from the
policy constructed this way and the weight corresponding to this action is
propagated further in the hidden state of the LSTM network, as given in Eq.
(23). The weights are then reset when the episode terminates.

Projective Simulations (PS) is a framework for agency and decision making that
has also found applications as a RL agent [68, 34]. In that context, a PS
agent interacts with an environment by performing actions sampled from an
action space $A$, whereas the environment provides the agent with perceptual
inputs, which reside in a percept space $S$, and reward signals. Its central
feature is represented by a so-called episodic and compositional memory (ECM)
– see Fig. 1 (b) –, a graphical model consisting of a network, or weighted
graph, where the vertices are clips and represent, e.g., remembered percepts
or remembered actions, but also more general states of the agent’s memory. In
this framework, the agent creates a clip inside the ECM each time it receives
a previously unknown input from the environment, or each time it creates a new
action clip or a more abstract clip, e.g., through action composition. This
makes it is possible to create ECM networks with complex graph structures,
allowing for generalization capabilities [33, 42]. Each input triggers a
random walk between the nodes of the ECM that is governed by the edge weights
of the graph. We assume in the following that the ECMs created are two-layered
networks, with one layer representing percepts clips and the other action
clips and where a percept at time $t$ of the RL interaction is connected
directly to all action clips. The edge weights are initialized uniformly:

$\displaystyle\ \forall a\in A,\forall s\in S:h_{0}(s,a)=1.$ (18)

We define the probability distribution111In standard PS, the probability is
defined without the exponential factors. over percepts $s\in S$ and actions
$a\in A$ as

$\displaystyle p(a|s)=\frac{e^{\beta h_{t}(s,a)}}{\underset{a}{\sum}e^{\beta
h_{t}(s,a)}},$ (19)

for a edge weight $h_{t}(s,a)$ connecting $s$ to $a$. The PS algorithm
optimizes the policy, similar to other RL frameworks, through a reward signal,
which is provided by the environment to the agent. At each RL time-step $t$,
the update rule is given by

$\displaystyle h_{t+1}(s,a)=h_{t}(s,a)-\gamma(h_{t}(s,a)-1)+$ (20)
$\displaystyle+g_{t}(s,a)R_{t},$

where $R_{t}$ is the reward at time step $t$, $\gamma$ is a damping
coefficient that regularizes the result and reduces potential instances of
trapping in local minima, $g_{t}(s,a)$ are so-called glow values [69], which
are initially set to zero

$\displaystyle\ \forall a\in A,\forall s\in S:g_{0}(s,a)=0.$ (21)

They are set (or reset) to $g_{t}(s,a)=1$ at time $t$ if the corresponding
edge is traversed and decay in value for each time step where they are not
used according to the rule

$\displaystyle g_{t+1}(s,a)=(1-\eta)g_{t}(s,a),$ (22)

where $0\leq\eta\leq 1$. The glow mechanism helps distribute the rewards along
the entire chain of state-action transitions traversed by the agent in an
episode. We see that the edge weights that are rewarded positively grow,
thereby enhancing the probability that the same action is chosen again in the
future upon receiving a similar percept.

PS has been successfully extended to include more powerful computational
structures, such as deep feedforward energy-based networks [39, 41] and
recurrent networks [40]. The aforementioned architectures enable the PS agent
to update the policy using function approximators. This allows the agent, as
it is the case for deep $Q$-learning [70], to construct more powerful
representations of the percept- and action-spaces and achieve high
performances in environments with, e.g., continuous parametric percepts and
actions without the need of discretization strategies [29]. We focus here on
the Long-Short Memory Network (LSTM) architecture [71, 72], a recurrent neural
network that is equipped with an internal memory architecture that helps it
learn long-range correlations in a sequential input. In this case, the action
sampled at time $t$ also depends on the LSTM internal state $u_{t}$, i.e.,

$\displaystyle a^{*}\sim p(a|s,u_{t})=\frac{e^{\beta
h_{t}(s,a,u_{t})}}{\underset{a}{\sum}e^{\beta h_{t}(s,a,u_{t})}}$ (23)

with $s\in S$ and $a,a^{*}\in A$ and the $u$-value is updated w.r.t. the
sampled action as follows:

$\displaystyle u_{t+1}=h_{t}(s,a^{*},u_{t}).$ (24)

The term $u_{t}=u_{t}(s,a)$ also depends on percepts and actions, but as we
see in Eq. (24), only the $u$-value corresponding to the action sampled by the
agent at time $t$ are propagated as information to the next RL time step, as
shown in Fig. 2. This value is the one that carries relevant information about
the correlations in the sequential structure of the percepts. We would like to
stress that the presence of an additional non-linear term in the update rule
in Eq. (23) induces a different type of ECM structure, where the term
$u_{t}(s,a)$ represents an additional edge weight. The update of the reward
mechanism can be generalized starting from the standard PS reward update
mechanism to fit the training of neural network-based policies, in analogy
with the case of $Q$-learning [39, 40].

##  IV Circuit compilation

Quantum compilation is the general task of reproducing a general operation
$\mathcal{M}\in\mathbb{C}^{4^{n}\crossproduct 4^{n}}$ on a $n$-qubit quantum
processor. In particular, we consider the compilation of unitary operations
$U\in U(n)$.

####  IV.0.1 Layer-based compilation and heuristic search

A general approach for gate synthesis in trapped-ion circuits is discussed in
Ref. [48]. This approach is based on the universality of two-qubit entangling
gates for quantum computation [73]. As a consequence, it is possible to
compile a quantum algorithm in growing layers, i.e., by iteratively placing
entangling MS gates on the circuit followed in each case by a collective
rotation of the following type

$\displaystyle\mathcal{R}(\theta_{k},\phi_{k},...,\theta_{k+n+1},\phi_{k+1})=$
(25)
$\displaystyle=C_{xy}(\theta_{k},\phi_{k})\prod_{i=1}^{n}Z_{i}(\theta_{k+i})C_{xy}(\theta_{k+n+1},\phi_{k+1}),$

where the gates $C_{xy}$ and $Z$ are defined in Eqs. (2)-(3). For each layer
of gates present on the circuit, we have one MS gate and one general rotation
$\mathcal{R}$ – see Eq. (25) –, which consists of $n$ local $Z$ gates and two
$C_{xy}$ gates. The unitary added to the circuit at each layer $l$ is:

$\displaystyle
V_{l}(\boldsymbol{\alpha_{l}})=\text{MS}(\theta^{\text{MS}}_{l},\phi^{\text{MS}}_{l})\
\mathcal{R}(\theta_{l},\phi_{l},...,\theta_{l+n+1},\phi_{l+1}).$ (26)

The cost function

$\displaystyle
C(\bm{\alpha})=1-\frac{1}{d^{2}}\absolutevalue{\tr{\prod_{l=1}^{L}V(\bm{\alpha}_{l})U^{\dagger}}}^{2},$
(27)

based on the gate fidelity in Eq. (4), has to be minimized with respect to the
angle parameters $\boldsymbol{\alpha}$. After obtaining the optimal parameters
$\boldsymbol{\alpha}^{*}$, if the desired error threshold is reached, the
algorithm terminates, otherwise a new layer of gates is placed on the circuit,
up to a maximum number of layers $L_{\text{MS}}$. By running this procedure
iteratively with different angle parameter sets, we can search through the
optimization landscape to find different gate decompositions (see Algorithm 1
and Ref. [48]).

Algorithm 1 Exhaustive search with layer-based compilation

Input Target unitary $U$, threshold $\epsilon$, cost function $C$.
Output $V^{*},\bm{\alpha}^{*}$.

1:$\triangleright$ Construct layers:

2:for $L=1$ to $L_{\text{MS}}$ do

3: $\triangleright$ Add MS gate and rotations:

4: $V(\boldsymbol{\alpha})=\prod_{l=1}^{L}V_{l}(\boldsymbol{\alpha}_{l})$

5: $\triangleright$ All angles but the ones of MS gates:

6: $\bm{\alpha}=(\theta_{1},\phi_{1},...,\theta_{(n+2)L},\phi_{2L})^{T}$

7: $\tilde{K}=(n+2)(L+1)$

8: $\triangleright$ Loop over numbers of rotations:

9: for $k=1$ to $\tilde{K}$ do

10: $\triangleright$ Loop over combinations of indices:

11: for $j=1$ to $\binom{\tilde{K}}{k}$ do

12: $\triangleright$ Set angles with chosen indices to zero:

13: $\alpha_{\sigma(j,k)}=0$

14:
$\tilde{\bm{\alpha}}^{k}=(...,\alpha_{\sigma(j,k)}=0...,\alpha_{\sigma(j,k)}=0,...)^{T}$

15: $\triangleright$ Cost function according to Eq. (27):

16: $C(\tilde{\bm{\alpha}}^{k})=1-F(V,U)$

17:
$C^{*}=\min(C(\tilde{\bm{\alpha}}^{k}),\nabla_{\bm{\alpha}^{k}}C(\tilde{\bm{\alpha}}^{k}))$

18: if $C^{*}\leq\epsilon$ then

19: $\boldsymbol{\alpha}^{*}=(\tilde{\bm{\alpha}}^{k})^{*}$

20: $V^{*}=V$

21: break

22: else

23: continue

24: end if

25: end for

26: end for

27:end for

However, while this method can be applied to circuits with a small number of
layers and qubits, its use becomes impractical as these two parameters grow.
In fact, the number of total combinations to analyze for a given number of MS
layers $L_{\text{MS}}$ and a number of qubits $n$ is

$\displaystyle
N_{\text{combinations}}=\frac{2^{(n+2)(L_{\text{MS}}+1)}-1}{1-\frac{1}{4}^{(n+2)}}-L_{\text{MS}}^{2}.$
(28)

For a $3$-qubit circuit with $L$ layers of entangling gates, the number of
possible circuits is $N_{\text{combinations}}=1049591$. For a $4$-qubit
circuit with $5$ entangling gates the number of combinations is approximately
$N_{\text{combinations}}\sim 2^{36}$. We see that already for $4$-qubit gates,
a brute force approach is unfeasible. Moreover, due to the large search space,
even approaches based on random search are not a viable option, especially if
the number of entangling gates is large.
In the following section, we suggest a different approach to the optimization
scheme where the position of the discrete parameters is modified by a RL agent
equipped with a curriculum scheme, whereas the continuous parameters are
optimized at each iteration and using a pre-defined heuristic for angle
initialization. This can offer benefits in several situations where the number
of gates is large enough to make the use of automated search prohibitive but
not so large to make the problem completely intractable from the point of view
of continuous parameter optimization.

Algorithm 2 PS-based compilation

Input Target unitary $U$, Action set
$G_{1}=\text{MS},G_{2}=C_{xy},G_{3}=Z_{1},...,G_{n+2}=Z_{n}$, threshold
function $\epsilon_{t}$, cost function $C$
Output $V^{*},\bm{\alpha}^{*}$

1:for $e=1$ to $E_{\text{max}}$ do

2: $s_{1}=(0,...,0),V_{1}=I_{d}$

3: for $t=1$ to $L_{\text{max}}$ do

4: $\triangleright$ Sample action according to Eq. (19):

5: $a_{t}\sim\pi(a|s_{t})$

6: $V_{t+1}(\bm{\alpha}_{t})=G_{a_{t}}V_{t}(\bm{\alpha}_{t-1})$

7: $s_{t+1}=(a_{1},...,a_{t},0,...,0)$

8:
$\boldsymbol{\alpha}=(\boldsymbol{\alpha}_{1},...,\boldsymbol{\alpha}_{t})\sim
2\pi\cdot\mathcal{N}(\boldsymbol{0}_{\boldsymbol{\alpha}},\boldsymbol{I}_{\boldsymbol{\alpha}})$

9: $\triangleright$ According to Eq. (27):

10:
$\boldsymbol{\alpha}^{*}=\text{argmin}(C(\bm{\alpha}),\nabla_{\boldsymbol{\alpha}}C(\bm{\alpha}))$

11: $R_{t}=\begin{cases}2&\text{if $\epsilon_{\text{min}}\leq
C(\boldsymbol{\alpha}^{*})\leq\epsilon_{t}$}\\\ 10&\text{if
$C(\boldsymbol{\alpha}^{*})\leq\epsilon_{\text{min}}$}\\\
0&\text{otherwise}\end{cases}$

12: if $C^{*}\leq\epsilon_{t}$ then

13: $V_{*}=V_{t+1}$

14: break

15: end if

16: Update rule for all $h_{t+1}(s_{t},a_{t})$ (see Eq. (20))

17: end for

18: Update threshold $\epsilon_{t}$ according to Eq. (30)

19:end for

####  IV.0.2 Reinforcement learning-based compilation

In the following, we discuss the implementation of RL-based compilation. In
this system, the agent acts on the environment, which represents a quantum
circuit, by choosing a gate from, e.g., the gate set of Section II and placing
it on the circuit. As an observation, the agent receives information about the
environment internal state. This can vary depending on the task to be solved:
in Ref. [22], e.g., the input is a single-qubit unitary and the task is to
construct a pre-trained optimal compiler that constructs any unitary with
minimal average number of actions. In Refs. [40, 74], the agent receives an
encoded representation of the quantum circuit in terms of gates and qubits.
The circuit-based input has the advantage of scaling linearly with the number
of qubits for $n$-qubit entangling gates and its sequential structure makes it
suitable for recurrent or autoregressive policies [28, 25]. The perceptual
input of the PS-LSTM agent is given by the current gate on the circuit. The
action of the PS-LSTM agent is placing one further gate on the circuit. A
representation of the interaction between the agent and the quantum circuit
environment is given in Fig. 1. At each time step, the agent can place one or
more gates on the circuit. Then the circuit is mapped to the corresponding
unitary function, which depends on parameters $\boldsymbol{\alpha}$. The
gradient-based optimizer, i.e., the L-BFGS-B algorithm [75], is given the task
of finding the optimal set of parameters to maximize the fidelity between the
current circuit and a target unitary process. The reward function and the
percept for the next interaction step are constructed based on the result of
the optimization. The algorithm is shown in Algorithm 2 for the standard PS
method and can be easily generalized to a framework with deep networks [39,
40].

Furthermore, one may design different types of RL environments based on the
way the agent places gates on the circuit. We develop here two different
approaches for two different types of RL environments: In the first
environment, we just mimic the structure of layer-based compilation, in which
we keep the entangling gates fixed and then let the agent vary the rotations
between them. This assumes that the RL interaction is defined by the number of
gates placed between each entangling layer and the total amount of entangling
layers considered. This architecture has the advantage of restricting the
search space of the RL agent, but it allows for less exploration of the cost
function landscape.

Figure 3: Quantum circuit optimization of different gates using the gate set
defined in Eqs. (1)-(3). (a), (e) show the average reward; (b), (f) the
average circuit size (thick line) and the size of the best circuit found so
far by all agents (dotted line); (c), (g) the number of optimal sequences per
episode; (d), (h) are histograms which show the distribution of the optimal
gate sequences over the circuit size. All lines refer to simulations of the
$3$-qubit Toffoli gate, first compiled on a $3$-qubit (green line, average
over 10 agents and 20 episodes) and then on a $4$-qubit circuit (orange line,
average over 5 agents and 20 episodes). Due to the $n$-qubit interaction of
the trapped-ion gates, the optimal sequence that generates the gate – which we
define as the shortest sequence whose infidelity falls below
$\epsilon_{\text{min}}=10^{-2}$ – increases in size from $10$ to $19$ gates.
Shaded regions in the plots represent the corresponding standard deviations.
While the average fidelity increases to reach the maximum for both circuits,
we see that the average circuit length appears to be higher than the shortest
circuit length. There are possible explanations for this: the shortest
sequences can be harder to optimize, so there is a higher chance that the
optimizer fails at outputting the cost function minimum for a given circuit
structure and the curriculum scheme, that modifies the problem online during
the training. Overall, we observe that the size of the optimal circuit
decreases as training progresses, thus validating the successful optimization
of the policy.

The second architecture does not restrict the action of the agent on the
quantum circuit, with one single exception: if the agent places the same type
of gate twice in a row on the circuit, these two gates are merged together in
one single gate. This is needed to avoid the agent getting stuck in a loop
where it keeps performing the same operation over and over again, without any
meaningful exploration of the optimization landscape from the perspective of
the continuous optimizer. We employ this architecture in our simulations.

In our implementation, the action $a_{t}$ chosen at time $t$ is the index of
certain gate in the gate set, whereas the corresponding percept $s_{t}$ is
given by the concatenation of previously chosen actions $(a_{1},...,a_{t-1})$.
Thus, the action space is given by the gate indices $A=\\{1,2,...,n+2\\}$ and
the percept space by the Cartesian product of $L_{\text{max}}$ action spaces:
$S=A^{\times L_{\text{max}}}$. The PS-LSTM algorithm, however, can also be
given just the action at time $t-1$ as percept, since the LSTM memory can
automatically capture the correlations between the elements in the sequence.
Here, we adopt a RL training procedure with a curriculum scheme as described
in Ref. [24]. This allows the agent to sufficiently explore the solution space
and gradually adapt the solution. More specifically, for each task of quantum
circuit optimization, we define a curriculum strategy where we reward the
agent at time step $t$ each time it finds a sequence with achieved minimal
infidelity $C(\boldsymbol{\alpha}^{*})$ lower than a chosen moving threshold
$\epsilon_{t}$ and a fixed threshold $\epsilon_{\text{min}}=10^{-2}$:

$\displaystyle R_{t}=\begin{cases}2&\text{if
$\epsilon_{\text{min}}<C(\boldsymbol{\alpha}^{*})<\epsilon_{t}$}\\\
10&\text{if $C(\boldsymbol{\alpha}^{*})<\epsilon_{\text{min}}$}\\\
0&\text{otherwise}\end{cases}.$ (29)

The episode terminates when the reward the cost function minimum in a given
time step falls below the threshold $\epsilon_{t}$ or when the maximal length
of the circuit per episode, $L_{\text{max}}$, is reached. The RL training
terminates upon reaching the maximum number of episodes $E_{\text{max}}$. The
reward scheme helps to progressively increase the fidelity throughout training
without allowing for too long circuits. The threshold is then lowered as
episodes progress based on previous rewards obtained by the agent. In our
implementation, we lower the threshold when it has been surpassed by the agent
at least 500 times using the following scheme:

$\displaystyle\epsilon_{t+1}=\begin{cases}\epsilon_{\text{min}}+\frac{1}{2}(\epsilon_{\text{min}}-\epsilon_{t})&\text{if
$\epsilon_{\text{min}}\leq\epsilon_{t}\leq 1$}\\\
\epsilon_{\text{min}}&\text{if $\epsilon_{t}\leq\epsilon_{\text{min}}$}\\\
1&\text{otherwise}.\end{cases}$ (30)

##  V Results

We test our algorithm on two relevant tasks of circuit optimization: Standard
gate compilation, which can be useful in particular for experimental
applications of frequently used gates (Toffoli, etc.) and the simulation of
black-box unitary processes with quantum circuits. The latter framework is
particularly interesting for quantum simulation and offers the possibility of
implementing our algorithm directly on an experimental setting where a black-
box unitary process has to be simulated by a quantum circuit with available
gates.

Figure 4: (a) Fidelity; (b) circuit size (thick lines) and size of the best
circuit found by all agents (dotted lines); (c) number of optimal sequences
per episode and (d) histogram of the optimal sequences – i.e., whose
infidelity falls below $\epsilon_{\text{min}}=10^{-2}$ – based on the circuit
size for the UCC operators defined in Eq. (31) for $\beta=\frac{\pi}{2}$ and
for 3, 4, 5, 6, 7, 8 qubits. The learning curves show the average over 20
agents, after which an average over a time window of 20 episodes is performed.
We observe how the average maximal fidelity reached by the agents decreases as
a function of the number of qubits, whereas the optimal circuit size
increases, which means that a longer circuit is necessary to synthesize the
desired operator. Moreover, the fraction of optimal sequences found by the
agent is also decreasing for higher numbers of qubits, as optimal policies
become more and more sparse and thus harder to discover for the PS-LSTM agent.
The hardest task for the PS-LSTM agent proves to be $7$-qubit UCC operator,
where the learning curve of the agents fails to converge, while, e.g., for the
$8$-qubit UCC unitary the agent can find sequences with very high fidelities
($1-F<10^{-5}$), even though the convergence is sub-optimal. We also observe
that for each learning curve there is a dip at some point in the training:
this is most likely due to the curriculum scheme, which modifies the reward
threshold during training and therefore influences the size of the optimal
circuits found by the agent. However, when considering the ensemble of
simulations, we see that the agents are capable of finding shorter and shorter
optimal circuits. Shaded regions in the plots represent the corresponding
standard deviations.

###  V.1 Example 1: gate compilation

As a first application of the method presented above, we consider the problem
of compiling a gate on an ion-trap quantum processor using the gate set given
in Eqs. (1)-(3). For this class of tasks, we choose two gate compilation
problems which can be of interest for typical applications, before passing to
a more general framework which considers black-box unitaries. First, we
consider a standard quantum computing gate, the $3$-qubit Toffoli gate [76,
77] – see Fig. 3 – for a $3$-qubit (green line, upper four plots) and a
$4$-qubit circuit (orange line, lower four plots). The blue dashed lines show
the optimal solutions found by the agent. In the $3$-qubit case, this solution
matches agrees with the results given in Ref. [48]. Fig. 3 (a), (e) show the
fidelity of the sequences produced and optimized by the agents as the learning
progresses, (b), (f) show the average circuit size and the size of the
shortest circuit found by the agent, (c), (g) show the average number of
optimal sequences, i.e., with fidelity higher than 0.99, found in each episode
(which in the best-case scenario should converge to one per episode) and (d),
(h) show histograms representing the number of sequences for different bins of
circuit size. We observe from the first and second plot in each row of Fig. 3
that the agents are capable of maximizing the fidelity and at the same time of
reducing the average circuit size in both cases. Moreover, while the average
circuit size is larger than the size of the best circuit, most likely due to
the presence of the curriculum and the influence of local minima on the angle
optimization, we also see that the agents find shorter and shorter optimal
circuits as the training progresses, hinting that the at least one agent in
the ensemble is indeed learning to optimize the circuit correctly. Moreover,
the third plot in each row shows us that the agents find an increasing number
of high-fidelity sequences during training. The minimal sequences found by the
agents are located in the leftmost tail of the distribution. The $3$-qubit
Toffoli gate implemented on a $4$-qubit circuit could have in principle a
relatively long generating quantum circuit in the chosen gate set, since the
4-qubit gate set also affects the qubits left unchanged by the Toffoli gate.
We also observe, in our simulations, that sequences generating this gate are
sparse in the action space, which could make it generally difficult to produce
this gate on a register of $n$ qubits without the help of MS and $C_{xy}$
gates acting only on subspaces of the register. We test whether our method can
discover a circuit of reasonable size. For large numbers of qubits and large
circuit sizes, we observe that although the algorithm can find shorter
circuits, it is still impaired by the computational overhead of simulating
such circuits exactly. In this case, the gate decomposition method introduced
in Section III.1 provides a useful speedup for RL simulations, especially if
compiled on GPU architectures.

###  V.2 Example 2: General quantum process simulation of a Hamiltonian model

As a second example, we compile parametric-type operators that play a relevant
role in quantum chemistry [78, 79], i.e., UCC-type operators. These are part
of a more general class of many-body operators [80]:

$\displaystyle
U(\beta)=\exp(i\beta\left(\prod_{i=1}^{n}(\sigma^{(i)}_{x}-i\sigma^{(i)}_{y})+\text{h.c.}\right)).$
(31)

Figure 5: These plots summarize the result of simulating the XXZ-model unitary
with our approach for different parameter configurations for $n=3,4,5$ qubits
(3 agents for each parameter sample). (a) shows the maximal value of the
fidelity found by the agents, (b) the mean fidelity, (c) the size of the
shortest circuit associated with the fidelity in (a) and (d) the average
circuit size, as a function of the coupling $J$ of the XXZ model for two
different configuration of the transverse coupling $\Delta$ = 0.5 with varying
$J$ and and $J$ = 0.5 with varying $\Delta$ (dark green, orange and blue and
pink, light green and yellow, respectively). The evolution time was fixed to
$\tau=0.25$ and the maximum size of the circuit to 50. The average fidelity is
calculated over the last 200 episodes. Vertical bars show the standard error
for the mean values of fidelity and circuit size. We see here that the circuit
size increases dramatically as the parameters $\Delta$ and $J$ increase. We
also see from the number of outliers in (c) and (g), that for certain values
of the parameters it is hard for the agents to exactly retrieve the optimal
circuit.

The class of operators defined by Eq. (31) resembles the collective rotations
used in trapped-ion gate sets given in Eqs. (1)-(3), they are however sparser.
Due to their importance, they represent our first candidate for process
simulation on the quantum circuit. The results for the simulation using PS-
LSTM of several such operators – $n=3,4,5,6,7,8$ – and for
$\beta=\frac{\pi}{2}$ is given in Fig. 4. We observe that the agent is able to
increase the fidelity up to the optimal value. We also see that the size of
the optimal circuit generating the corresponding operators increases with the
number of qubits. Furthermore, due to the growing number of local rotations
available to the agent and the presence of vanishing gradients, it becomes
increasingly difficult to find and reward optimal sequences. This also shows
the benefit of implementing curriculum strategies in such hybrid discrete-
continuous optimization problems, where the landscape both for the RL agent
and for the numeric optimizer become increasingly sparse. We also note that
the generation of $7$ and $8$-qubit UCC operators proves more challenging. In
particular, for the $7$-qubit UCC unitary, the agent is able to raise the
fidelity to values close to $F=0.99$, but it cannot significantly increase the
number of circuits with $F>0.99$ over the course of the training. In the case
of the $8$-qubit UCC operator, the agent instead proves capable of discovering
sequences with very high fidelity, i.e $1-F<10^{-5}$, although the average
fidelity worsens slightly towards the end of the training. We also see that an
ensemble of agents can successfully reduce the size of the best circuit, even
in those cases where the average performance worsens during training.

We would like now to consider a (black-box) unitary $U(\tau)=e^{-iH\tau}$,
where $\tau$ is the evolution time described by a Hamiltonian of the following
type:

$\displaystyle
H=-2h\sum_{i=1}^{n}\sigma^{(i)}_{z}-J\sum_{i=1}^{n-1}\Big{(}\sigma^{(i)}_{x}\otimes\sigma^{(i+1)}_{x}+\sigma^{(i)}_{y}\otimes\sigma^{(i+1)}_{y}$
(32) $\displaystyle+$
$\displaystyle\Delta\left(\sigma^{(i)}_{z}\otimes\sigma^{(i+1)}_{z}-\frac{1}{4}I\right)\Big{)},$

which is usually referred to as the XXZ model [81]. We want to analyze how the
optimal quantum circuit found by the RL agent changes when we vary the
Hamiltonian parameters. For the XXZ model, this parameter variation represents
for instance the transition from a XX model when $\Delta=0$ to a XXX model for
$\Delta=1$, or from $ZZ$ interaction for $J=0$ to a XXZ model for $J>0$. In
fact, we observe in different parameter regions different behaviours of the
circuit structure. Results of several RL runs for two different
configurations, one with fix coupling $J=0.5$ and varying $\Delta$ and another
one with fix $\Delta=0.5$ and varying $J$ of the XXZ model unitary are shown
in Fig. 5 for 3, 4 and 5 qubits. Plots (a), (e) show the maximal fidelity
found by the agents for each run, plots (b), (f) the mean fidelity, plots (c),
(g) the number of gates in the shortest circuit among those with the highest
fidelity and plots (d), (h) the average circuit size, all represented in their
functional dependence from the coupling $J$ and transverse coupling $\Delta$
of the XXZ model for two different configuration of the transverse coupling
$\Delta$ = 0.5 with varying $J$ and and $J$ = 0.5 with varying $\Delta$ (dark
green, orange and blue and pink, light green and yellow, respectively). We
also see how rapidly the circuit size grows in the presence of entanglement,
for example from $J=0$ to $J=1$, whereas the increase in average circuit size
seems less pronounced as we vary $\Delta$. We observe that the optimal circuit
size can also decrease, for example when $\Delta=1$ and $J=0.5$. We also see
some instabilities and high variance in both the average circuit size and the
size of the best circuit discovered, though it is hard to determine whether
the agent fails at finding a shorter circuit for a specific parameter or the
problem becomes suddenly harder to represent with the given gate set due to
the parameter variation.

##  VI Conclusion

In this work we construct a framework to learn both classically simulated and
black-box unitaries on a quantum circuit using RL and unconstrained
optimization. Our simulation is specifically tailored to an ion-trap
architecture based on collective gates and local rotations. We demonstrate the
synthesis of optimal circuits for Toffoli gates and UCC operators for varying
numbers of qubits. As instances of black-box unitaries, we also consider
Hamiltonian simulations of the XXZ model. More specifically, we study the
convergence and the quality of the solutions found by agent and optimizer as
we modify relevant parameters of the underlying black-box unitary, such as the
coupling in the XXZ Hamiltonian. After testing the algorithm on different
unitary process simulation scenarios, we observe that the optimization of
circuits is generally possible even for large numbers of gates, it is however
difficult to foresee how sparse the cost function minima can be as we increase
the number of qubits and reduce the sparsity in the corresponding Hamiltonian.
Possible improvements include combining the RL search with a graph traversal
algorithm to have a more efficient exploration of the discrete action space
[36, 27]. We expect that this approach can be applied to different
architectures beyond trapped ions, and more carefully engineered reward
functions can be developed to further enhance the discovery of optimal
circuits.  Our unified optimization, combining circuit compilation and unitary
synthesis, may find use both in classical pre-optimization and in experimental
circuit learning, be it for in-situ (variational) algorithms or module-
compilation tasks.

##  VII Data and Code availability

The code and the data are available at the following link
https://github.com/franz3105/RL_Ion_gates.

## Acknowledgement

We would like to thank Markus Schmitt, Robert Zeier and Marco Canteri for the
useful discussions. FP, MS and FM acknowledge support from the German Federal
Ministry of Education and Research (BMBF) within the framework programme
”Quantum technologies – from basic research to market” (Project QSolid, Grant
No. 13N161), the European Union’s Horizon Programme (HORIZON-
CL4-2021-DIGITALEMERGING-02-10) Grant Agreement 101080085 QCFD, by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – Cluster of Excellence Matter and Light for
Quantum Computing (ML4Q) EXC 2004/1 – 390534, and from the Jülich
Supercomputing Center through the JUWELS and JURECA cluster. SJ, HPN, LMT and
HJB acknowledge support from the Austrian Science Fund (FWF) through the
projects DK-ALM:W1259-N27 and SFB BeyondC F7102. SJ also acknowledges the
Austrian Academy of Sciences as a recipient of the DOC Fellowship. Simulations
were realized using the libraries NUMBA [82], JAX [83], PYTORCH [84] and SCIPY
[85]. The robot picture in Fig. 1 was generated by the authors using
Midjourney.

## References

  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Bravo-Prieto  _et al._ [2019] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, Variational quantum linear solver, arXiv preprint at arXiv:1909.05820 (2019).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 4213 (2014).
  * Cerezo  _et al._ [2021a] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021a).
  * Kreppel  _et al._ [2022] F. Kreppel, C. Melzer, J. Wagner, J. Hilder, U. Poschinger, F. Schmidt-Kaler, and A. Brinkmann, Quantum circuit compiler for a shuttling-based trapped-ion quantum computer (2022).
  * Venturelli  _et al._ [2019] D. Venturelli, M. Do, B. O’Gorman, J. Frank, E. Rieffel, K. E. C. Booth, T. Nguyen, P. Narayan, and S. Nanda, Quantum circuit compilation: An emerging application for automated reasoning, in  _Scheduling and Planning Applications Workshop_ (2019).
  * Wille  _et al._ [2020] R. Wille, S. Hillmich, and L. Burgholzer, Efficient and correct compilation of quantum circuits, in  _2020 IEEE International Symposium on Circuits and Systems (ISCAS)_ (2020) pp. 1–5.
  * Maronese  _et al._ [2021] M. Maronese, L. Moro, L. Rocutto, and E. Prati, Quantum compiling, arXiv preprint arXiv:2112.00187 (2021).
  * Sørensen and Mølmer [1999] A. Sørensen and K. Mølmer, Quantum computation with ions in thermal motion, Physical Review Letters 82, 1971–1974 (1999).
  * Mølmer and Sørensen [1999] K. Mølmer and A. Sørensen, Multiparticle entanglement of hot trapped ions, Phys. Rev. Lett. 82, 1835 (1999).
  * Malekakhlagh [2021] M. Malekakhlagh, How the Cross Resonance Gate Works, https://thequantumaviary.blogspot.com/2021/07/how-cross-resonance-gate-works.html (2021), [Online; accessed 06-January-2022].
  * Dallaire-Demers and Wilhelm [2016] P.-L. Dallaire-Demers and F. K. Wilhelm, Quantum gates and architecture for the quantum simulation of the fermi-hubbard model, Phys. Rev. A 94, 062304 (2016).
  * Preti  _et al._ [2022] F. Preti, T. Calarco, and F. Motzoi, Continuous quantum gate sets and pulse-class meta-optimization, PRX Quantum 3, 040311 (2022).
  * Daoud  _et al._ [2022] M. S. Daoud, M. Shehab, H. M. Al-Mimi, L. Abualigah, R. A. Zitar, and M. K. Y. Shambour, Gradient-based optimizer (gbo): A review, theory, variants, and applications, Archives of Computational Methods in Engineering (2022).
  * Simon [2013] D. Simon,  _Evolutionary optimization algorithms_ (Wiley-Blackwell, Hoboken, NJ, 2013).
  * Nelder and Mead [1965] J. A. Nelder and R. Mead, A Simplex Method for Function Minimization, The Computer Journal 7, 308 (1965).
  * Dawson and Nielsen [2006] C. M. Dawson and M. A. Nielsen, The solovay-kitaev algorithm, Quantum Info. Comput. 6, 81–95 (2006).
  * Vatan and Williams [2004] F. Vatan and C. Williams, Optimal quantum circuits for general two-qubit gates, Phys. Rev. A 69, 032315 (2004).
  * Wille  _et al._ [2016] R. Wille, O. Keszocze, M. Walter, P. Rohrs, A. Chattopadhyay, and R. Drechsler, Look-ahead schemes for nearest neighbor optimization of 1d and 2d quantum circuits, in  _2016 21st Asia and South Pacific Design Automation Conference (ASP-DAC)_ (2016) pp. 292–297.
  * Davis  _et al._ [2020] M. G. Davis, E. Smith, A. Tudor, K. Sen, I. Siddiqi, and C. Iancu, Towards optimal topology aware quantum circuit synthesis, in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (2020) pp. 223–234.
  * Moro  _et al._ [2021] L. Moro, M. G. A. Paris, M. Restelli, and E. Prati, Quantum compiling by deep reinforcement learning, Communications Physics 4, 178 (2021).
  * Fösel  _et al._ [2021] T. Fösel, M. Y. Niu, F. Marquardt, and L. Li, Quantum circuit optimization with deep reinforcement learning, arXiv preprint arXiv:2103.07585 (2021).
  * Ostaszewski  _et al._ [2021a] M. Ostaszewski, L. M. Trenkwalder, W. Masarczyk, E. Scerri, and V. Dunjko, Reinforcement learning for optimization of variational quantum circuit architectures, in  _Advances in Neural Information Processing Systems_, Vol. 34, edited by M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan (Curran Associates, Inc., 2021) pp. 18182–18194.
  * Sivak  _et al._ [2022] V. Sivak, A. Eickbusch, H. Liu, B. Royer, I. Tsioutsios, and M. Devoret, Model-free quantum control with reinforcement learning, Physical Review X 12 (2022).
  * Borah  _et al._ [2021] S. Borah, B. Sarma, M. Kewming, G. J. Milburn, and J. Twamley, Measurement-based feedback quantum control with deep reinforcement learning for a double-well nonlinear potential, Phys. Rev. Lett. 127, 190403 (2021).
  * Zhang  _et al._ [2020] Y.-H. Zhang, P.-L. Zheng, Y. Zhang, and D.-L. Deng, Topological quantum compiling with reinforcement learning, Phys. Rev. Lett. 125, 170501 (2020).
  * Yao  _et al._ [2021] J. Yao, L. Lin, and M. Bukov, Reinforcement learning for many-body ground-state preparation inspired by counterdiabatic driving, Phys. Rev. X 11, 031070 (2021).
  * Sutton and Barto [2018] R. S. Sutton and A. G. Barto,  _Reinforcement Learning: An Introduction_, 2nd ed. (The MIT Press, 2018).
  * Silver  _et al._ [2017] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, Y. Chen, T. Lillicrap, F. Hui, L. Sifre, G. van den Driessche, T. Graepel, and D. Hassabis, Mastering the game of go without human knowledge, Nature 550, 354 (2017).
  * Jumper  _et al._ [2021] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, A. Bridgland, C. Meyer, S. A. A. Kohl, A. J. Ballard, A. Cowie, B. Romera-Paredes, S. Nikolov, R. Jain, J. Adler, T. Back, S. Petersen, D. Reiman, E. Clancy, M. Zielinski, M. Steinegger, M. Pacholska, T. Berghammer, S. Bodenstein, D. Silver, O. Vinyals, A. W. Senior, K. Kavukcuoglu, P. Kohli, and D. Hassabis, Highly accurate protein structure prediction with alphafold, Nature 596, 583 (2021).
  * Briegel and De las Cuevas [2012] H. J. Briegel and G. De las Cuevas, Projective simulation for artificial intelligence, Scientific Reports 2, 400 (2012).
  * Melnikov  _et al._ [2017] A. A. Melnikov, A. Makmal, V. Dunjko, and H. J. Briegel, Projective simulation with generalization, Scientific Reports 7, 14430 (2017).
  * Melnikov  _et al._ [2018] A. A. Melnikov, H. Poulsen Nautrup, M. Krenn, V. Dunjko, M. Tiersch, A. Zeilinger, and H. J. Briegel, Active learning machine learns to create new quantum experiments, Proceedings of the National Academy of Sciences 115, 1221 (2018).
  * Wallnöfer  _et al._ [2020] J. Wallnöfer, A. A. Melnikov, W. Dür, and H. J. Briegel, Machine learning for long-distance quantum communication, PRX Quantum 1, 010301 (2020).
  * Dalgaard  _et al._ [2020a] M. Dalgaard, F. Motzoi, J. J. Sørensen, and J. Sherson, Global optimization of quantum dynamics with alphazero deep exploration, npj Quantum Information 6, 6 (2020a).
  * Nautrup  _et al._ [2019] H. P. Nautrup, N. Delfosse, V. Dunjko, H. J. Briegel, and N. Friis, Optimizing Quantum Error Correction Codes with Reinforcement Learning, Quantum 3, 215 (2019).
  * Ried  _et al._ [2019] K. Ried, T. Müller, and H. J. Briegel, Modelling collective motion based on the principle of agency: General framework and the case of marching locusts, PLOS ONE 14, 1 (2019).
  * Jerbi  _et al._ [2021] S. Jerbi, L. M. Trenkwalder, H. P. Nautrup, H. J. Briegel, and V. Dunjko, Quantum enhancements for deep reinforcement learning in large spaces, PRX Quantum 2 (2021).
  * Preti [2020] F. Preti, Deep projective simulation and state preparation, master thesis (2020).
  * Nautrup  _et al._ [2022] H. P. Nautrup, T. Metger, R. Iten, S. Jerbi, L. M. Trenkwalder, H. Wilming, H. J. Briegel, and R. Renner, Operationally meaningful representations of physical systems in neural networks, Machine Learning: Science and Technology 3, 045025 (2022).
  * Eva  _et al._ [2022] B. Eva, K. Ried, T. Müller, and H. J. Briegel, How a minimal learning agent can infer the existence of unobserved variables in a complex environment, Minds and Machines (2022).
  * Bruzewicz  _et al._ [2019] C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage, Trapped-ion quantum computing: Progress and challenges, Applied Physics Reviews 6, 021314 (2019).
  * Harty  _et al._ [2014] T. P. Harty, D. T. C. Allcock, C. J. Ballance, L. Guidoni, H. A. Janacek, N. M. Linke, D. N. Stacey, and D. M. Lucas, High-fidelity preparation, gates, memory, and readout of a trapped-ion quantum bit, Phys. Rev. Lett. 113, 220501 (2014).
  * Ballance  _et al._ [2016] C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, and D. M. Lucas, High-fidelity quantum logic gates using trapped-ion hyperfine qubits, Phys. Rev. Lett. 117, 060504 (2016).
  * Erhard  _et al._ [2019] A. Erhard, J. J. Wallman, L. Postler, M. Meth, R. Stricker, E. A. Martinez, P. Schindler, T. Monz, J. Emerson, and R. Blatt, Characterizing large-scale quantum computers via cycle benchmarking, Nature Communications 10 (2019).
  * Paul [1990] W. Paul, Electromagnetic traps for charged and neutral particles, Rev. Mod. Phys. 62, 531 (1990).
  * Martinez  _et al._ [2016] E. A. Martinez, T. Monz, D. Nigg, P. Schindler, and R. Blatt, Compiling quantum algorithms for architectures with multi-qubit gates, New Journal of Physics 18, 063029 (2016).
  * Gaebler  _et al._ [2012] J. P. Gaebler, A. M. Meier, T. R. Tan, R. Bowler, Y. Lin, D. Hanneke, J. D. Jost, J. P. Home, E. Knill, D. Leibfried, and D. J. Wineland, Randomized benchmarking of multiqubit gates, Phys. Rev. Lett. 108, 260503 (2012).
  * Monroe  _et al._ [2014] C. Monroe, R. Raussendorf, A. Ruthven, K. R. Brown, P. Maunz, L.-M. Duan, and J. Kim, Large-scale modular quantum-computer architecture with atomic memory and photonic interconnects, Phys. Rev. A 89, 022317 (2014).
  * Ringbauer  _et al._ [2022] M. Ringbauer, M. Meth, L. Postler, R. Stricker, R. Blatt, P. Schindler, and T. Monz, A universal qudit quantum processor with trapped ions, Nature Physics 18, 1053 (2022).
  * Nielsen [2002] M. A. Nielsen, A simple formula for the average gate fidelity of a quantum dynamical operation, Physics Letters A 303, 249 (2002).
  * Dalgaard  _et al._ [2022] M. Dalgaard, F. Motzoi, and J. Sherson, Predicting quantum dynamical cost landscapes with deep learning, Phys. Rev. A 105, 012402 (2022).
  * Strassen [1969] V. Strassen, Gaussian elimination is not optimal, Numerische Mathematik 13, 354 (1969).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: design of nmr pulse sequences by gradient ascent algorithms, Journal of Magnetic Resonance 172, 296 (2005).
  * Dalgaard  _et al._ [2020b] M. Dalgaard, F. Motzoi, J. H. M. Jensen, and J. Sherson, Hessian-based optimization of constrained quantum control, Physical Review A 102, 042612 (2020b).
  * Morningstar  _et al._ [2022] A. Morningstar, M. Hauru, J. Beall, M. Ganahl, A. G. Lewis, V. Khemani, and G. Vidal, Simulation of quantum many-body dynamics with tensor processing units: Floquet prethermalization, PRX Quantum 3 (2022).
  * Chuang and Nielsen [1997] I. L. Chuang and M. A. Nielsen, Prescription for experimental determination of the dynamics of a quantum black box, Journal of Modern Optics 44, 2455 (1997).
  * Bužek and Hillery [1996] V. Bužek and M. Hillery, Quantum copying: Beyond the no-cloning theorem, Phys. Rev. A 54, 1844 (1996).
  * Khatri  _et al._ [2019] S. Khatri, R. LaRose, A. Poremba, L. Cincio, A. T. Sornborger, and P. J. Coles, Quantum-assisted quantum compiling, Quantum 3, 140 (2019).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018).
  * Cerezo  _et al._ [2021b] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature Communications 12, 1791 (2021b).
  * Jerbi  _et al._ [2023] S. Jerbi, J. Gibbs, M. S. Rudolph, M. C. Caro, P. J. Coles, H.-Y. Huang, and Z. Holmes, The power and limitations of learning quantum dynamics incoherently (2023), arXiv preprint:2303.12834 .
  * Wierichs  _et al._ [2022] D. Wierichs, J. Izaac, C. Wang, and C. Y.-Y. Lin, General parameter-shift rules for quantum gradients, Quantum 6, 677 (2022).
  * Ostaszewski  _et al._ [2021b] M. Ostaszewski, E. Grant, and M. Benedetti, Structure optimization for parameterized quantum circuits, Quantum 5, 391 (2021b).
  * Mazyavkina  _et al._ [2021] N. Mazyavkina, S. Sviridov, S. Ivanov, and E. Burnaev, Reinforcement learning for combinatorial optimization: A survey, Computers & Operations Research 134, 105400 (2021).
  * Li [2017] Y. Li, Deep reinforcement learning: An overview (2017).
  * Boyajian  _et al._ [2020] W. L. Boyajian, J. Clausen, L. M. Trenkwalder, V. Dunjko, and H. J. Briegel, On the convergence of projective-simulation–based reinforcement learning in markov decision processes, Quantum Machine Intelligence 2, 13 (2020).
  * Mautner  _et al._ [2015] J. Mautner, A. Makmal, D. Manzano, M. Tiersch, and H. J. Briegel, Projective simulation for classical learning agents: A comprehensive investigation, New Generation Computing 33, 69 (2015).
  * Mnih  _et al._ [2015] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, S. Legg, and D. Hassabis, Human-level control through deep reinforcement learning, Nature 518, 529 (2015).
  * Hochreiter [1991] S. Hochreiter, Untersuchungen zu dynamischen neuronalen Netzen. Diploma thesis (1991).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, Long short-term memory, Neural Computation 9, 1735 (1997).
  * DiVincenzo [1995] D. P. DiVincenzo, Two-bit gates are universal for quantum computation, Phys. Rev. A 51, 1015 (1995).
  * Bukov  _et al._ [2018] M. Bukov, A. G. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Reinforcement learning in different phases of quantum control, Physical Review X 8 (2018).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, On the limited memory bfgs method for large scale optimization, Math. Program. 45, 503 (1989).
  * Toffoli [1980] T. Toffoli, Reversible computing, in  _Automata, Languages and Programming_, edited by J. de Bakker and J. van Leeuwen (Springer Berlin Heidelberg, Berlin, Heidelberg, 1980) pp. 632–644.
  * Monz  _et al._ [2009] T. Monz, K. Kim, W. Hänsel, M. Riebe, A. S. Villar, P. Schindler, M. Chwalla, M. Hennrich, and R. Blatt, Realization of the quantum toffoli gate with trapped ions, Phys. Rev. Lett. 102, 040501 (2009).
  * Taube and Bartlett [2006] A. G. Taube and R. J. Bartlett, New perspectives on unitary coupled-cluster theory, International Journal of Quantum Chemistry 106, 3393 (2006).
  * Bartlett and Musiał [2007] R. J. Bartlett and M. Musiał, Coupled-cluster theory in quantum chemistry, Rev. Mod. Phys. 79, 291 (2007).
  * Motzoi  _et al._ [2017] F. Motzoi, M. P. Kaicher, and F. K. Wilhelm, Linear and logarithmic time compositions of quantum many-body operators, Phys. Rev. Lett. 119, 160503 (2017).
  * Takahashi [1999] M. Takahashi,  _Thermodynamics of One-Dimensional Solvable Models_ (Cambridge University Press, 1999).
  * Lam  _et al._ [2015] S. K. Lam, A. Pitrou, and S. Seibert, Numba: A llvm-based python jit compiler, in  _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_ (2015) pp. 1–6.
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Li  _et al._ [2017] J. Li, X. Yang, X. Peng, and C.-P. Sun, Hybrid quantum-classical approach to quantum optimal control, Phys. Rev. Lett. 118, 150503 (2017).
  * Barhate [2021] N. Barhate, Minimal pytorch implementation of proximal policy optimization, https://github.com/nikhilbarhate99/PPO-PyTorch (2021).

##  Appendix A Fast analytic ion gates

In this section, we derive the representations for the MS and $C_{xy}$ gates
based on their spectral decomposition and we show how these can be further
simplified. We also provide descriptions of the gate gradients based on the
optimized representations. The $n$-qubit XY-Rotation and MS Hamiltonians are
given by

$\displaystyle\hat{H}_{\text{XY}}(n,\phi)=$ $\displaystyle
S_{x}\cos(\phi)+S_{y}\sin(\phi),$ (33)
$\displaystyle\hat{H}_{\text{MS}}(n,\phi)=$
$\displaystyle\hat{H}_{\text{XY}}(n,\phi)^{2}=\left(S_{x}\cos(\phi)+S_{y}\sin(\phi)\right)^{2},$
(34)

We will first focus on solving the XY-rotation gate and then generalize our
solution to the MS gate, using the fact that the MS-Hamiltonian is the square
of the XY-Hamiltonian. A representation of the $XY$ unitary in terms of its
real and imaginary entries for $n=8$, $\theta=\frac{3}{4}\pi$ and $\phi=2\pi$
is given in Fig. 6.

###  A.1 General Approach - Constructing the XY-rotation Gate

Due to the lack of multi-qubit interaction terms in the XY-rotation
Hamiltonian, the eigenvectors and eigenvalues of the $n$-qubit case can be
constructed from the single qubit closed-form solution, allowing us to
factorize the evolution into $\phi$ dependent terms and $\theta$ dependent
terms ($\theta$ can be understood as the time evolution parameter, using a
standard notation of the literature of trapped-ion quantum computing). This
allows us to calculate the associated unitaries from element-wise operations.
We decompose the Hamiltonian into

$\displaystyle\hat{H}_{\text{XY}}(n,\phi)=\hat{V}_{n}(\phi)\Lambda_{\text{XY}}\hat{V}_{n}^{\dagger}(\phi),$
(35)

with the $\phi$ dependent $\hat{V}_{n}$ and the diagonal matrix of the
eigenvalues $\Lambda_{\text{XY}}$. For the unitaries the eigenvalues then
capture the $\theta$ dependence via
$\exp(-i\Lambda_{\text{XY}}\frac{\theta}{2})$.

##### Eigenvalues

The eigenvalues are constructed from the single qubit case using the following
Kronecker sum notation $\hat{A}^{\oplus n}=\sum_{i=1}^{n}\mathds{I}^{\otimes
i-1}\otimes\hat{A}\otimes\mathds{I}^{\otimes n-i}$,

$\displaystyle\Lambda_{\text{XY}}=\text{diag}\left(2\mathbf{b}_{n}-n\right),\quad\text{with}\quad\mathbf{b}_{n}=\begin{pmatrix}0\\\
1\end{pmatrix}^{\oplus n}.$ (36)

We find that the $i$-th eigenvalue can be constructed from the binary Hamming
weight vector $\mathbf{b}_{n}$, because the $i$-th component of the binary
Hamming weight $b_{n}(i)$ corresponds to the number of qubits with $\ket{1}$
at index $i$ of the Hilbert space.

##### Eigenvectors

The eigenvector matrices $\hat{V}_{n}(\phi)$ are also constructed from the
single qubit case, but using the tensorproduct,

$\displaystyle\hat{V}_{n}(\phi)=\frac{1}{\sqrt{2^{n}}}\begin{pmatrix}1&1\\\
-e^{i\phi}&e^{i\phi}\end{pmatrix}^{\otimes
n}=\hat{P}_{n}(\phi)\odot\frac{1}{\sqrt{2^{n}}}\underbrace{(\hat{\mathds{I}}_{2}+i\sigma_{y})^{\otimes
n}}_{=\hat{S}_{n}}.$ (37)

We furthermore decompose them into an elementwise product ($\odot$) of a sign
$\hat{S}_{n}$ and phase $\hat{P}_{n}(\phi)$ component. The phase components
$\hat{P}_{n}(\phi)$ are column-independent
$\hat{P}_{n}(\phi)=[\mathbf{p}_{n}(\phi),\mathbf{p}_{n}(\phi),\cdots,\mathbf{p}_{n}(\phi)]$,
where the vectors themselves can be regarded as the element wise complex
exponential $v_{\phi}$ of the binary Hamming weight vector $\mathbf{b}_{n}$

$\displaystyle\mathbf{p}_{n}(\phi)=\begin{pmatrix}1\\\
e^{i\phi}\end{pmatrix}^{\otimes
n}=\exp(i\phi\mathbf{b}_{n})=v_{\phi}(\mathbf{b}_{n}).$ (38)

##### Calculating the Unitary

We can now use the column independence of the phase components
$\hat{P}_{n}(\phi)$, to commute it with the eigenvalues. This allows us to
separate the phase component further. For the unitary of the XY-gate we then
find

$\displaystyle C_{xy}(\phi,\theta)=$
$\displaystyle\exp\left(-i\hat{H}_{\text{XY}}(\phi)\theta\right)=\hat{V}_{n}(\phi)\exp(-i\Lambda_{\text{XY}})\theta)\hat{V}_{n}^{\dagger}(\phi)$
$\displaystyle=$
$\displaystyle\frac{1}{2^{n}}\underbrace{\hat{P}_{n}(\phi)\hat{P}_{n}^{\dagger}(\phi)}_{=v_{\phi}(\hat{C})}\hat{S}_{n}\exp(-\nicefrac{{i\theta}}{{2}}\Lambda_{\text{XY}})\hat{S}_{n}^{\dagger}.$
(39)

We substitute the phase term product with the previously defined element wise
complex exponential of $\hat{C}=\hat{P}_{n}-\hat{P}_{n}^{T}$. Furthermore we
split the sign components into degeneracy subspaces one for each different
eigenvalue in $\lambda_{\text{XY}}$. We decompose
$\hat{S}_{n}=\sum_{i=0}^{n}\hat{S}_{\lambda_{i}}$, where
$\hat{S}_{\lambda_{i}}$ retains the values of $\hat{S}_{n}$ for columns with
binary Hamming weight $i$, but is zero for all other columns. From this we can
then construct cached matrices
$\hat{D}_{\lambda_{i}}=\hat{S}_{\lambda_{i}}\hat{S}_{\lambda_{i}}^{\dagger}$,
so that we can rewrite

$\displaystyle\hat{S}_{n}\exp(-\nicefrac{{i\theta}}{{2}}\Lambda_{\text{XY}})\hat{S}_{n}^{\dagger}=\sum_{i=0}^{n}e^{-i\lambda_{\text{XY},i}\frac{\theta}{2}}\hat{D}_{\lambda_{i}}.$
(40)

From this we conclude

$\displaystyle C_{xy}(\phi,\theta)=$
$\displaystyle\frac{1}{2^{n}}v_{\phi}(\hat{C})\odot\sum_{i=0}^{n}\exp(-i\lambda_{\text{XY,i}}\theta)\phantom{,}\hat{D}_{\lambda_{i}}.$
(41)

The element wise exponential $v_{\phi}(\hat{C})$ is constructed efficiently by
reusing $2n+1$ phase values, because $C_{ij}\in[0,1,\cdots,2n]$. Each of the
required computations require $\mathcal{O}\left((2^{n})^{2}\right)$
operations, resulting in a total complexity of
$\mathcal{O}\left((n+1)(2^{n})^{2}\right)$.

Figure 6:  Real and imaginary components of the unitary matrix of the
$n$-qubit XY-rotation gate with $\phi=2\pi$ and
$\theta=\nicefrac{{3}}{{4}}\pi$. The colormap describes the value range of the
matrix elements.

###  A.2 Generalization to the MS gate

The approach used for $C_{xy}$ gate can be generalized to the MS gate. Due to
its Hamiltonian being the square of the aforementioned Hamiltonian, the two
gates share their eigenvectors, while the eigenvalues of the $C_{xy}$ gate are
squared $\lambda_{\text{XY,i}}^{2}=\lambda_{\text{MS, i}}$. This reduces the
number of different eigenvalues from $n+1$ to
$\lceil\frac{n}{2}\rceil+\mod(2+1,2)$. The expansion in Eq. (41) is
transformed into

$\displaystyle\text{MS}(\phi,\theta)=$
$\displaystyle\frac{1}{2^{n}}v_{\phi}(\hat{C})\odot\bigg{(}\sum_{i=0}^{\lceil\frac{n}{2}\rceil}\exp(-i\lambda_{\text{MS,i}}\theta)\phantom{,}\overbrace{\left(\hat{D}_{\lambda_{i}}+\hat{D}_{\lambda_{n-i}}\right)}^{=\tilde{D}_{i}}\delta^{\text{floor}}_{i}\bigg{)}.$
(42)

where

$\displaystyle\delta^{\text{floor}}_{i}=\begin{cases}1&\text{if $0\leq
i\leq\frac{n}{2}$}\\\ 0&\text{otherwise}.\end{cases}$ (43)

As both expansions, Eq. (41) and Eq. (42), rely solely on operations local to
the matrix elements, these operations are ideal for parallelisation. In Fig. 7
we show the improvement in terms of computation time (a) and speedup (b) of
our method compared to standard matrix exponentiation as a function of the
number of qubits : in (a) the orange line shows the computation time of the
matrix exponential, whereas the green line shows our method implemented
without caching the parameter-independent matrices, and the blue line our
method again, but with cached matrices; (b) shows instead the ratio between
the computation time of the cached method and standard matrix exponentiation.
We see that we can reach a speedup between 50 and 250 for $n\leq 12$, which
proves particularly useful in our quantum-circuit simulations: In fact, these
require large numbers of cost function evaluations, due to the presence of
both the reinforcement learning agent and the gradient-based optimization.

Figure 7: (a) Walltime for MS gate construction with standard matrix
exponentiation and with the analytical approach of this paper as a function of
the number of qubits. The orange line corresponds to the function scipy.expm,
the blue and green lines represent the analytical approach during caching and
after caching, respectively; (b) speedup of the new approach, once the
parameter-indedependent matrices have been cached, as a function of the number
of qubits. Here we see that the speedup is not constant as the number of
qubits changes, but we have a maximum speedup around 5 qubits, which then
decreases until we reach 8 qubits, and then increases again. This effect is
most likely due to an underlying slowdown in the underlying fundamental
operations, e.g. due to the processor cache being filled up quickly, so that
access to the RAM becomes necessary. Moreover, our method is implemented only
using NUMBA, i.e. compiled PYTHON code, whereas the exponential matrix
function of SCIPY profits from underlying time-critical routines written in C,
C++ and Fortran. It seems, however, that the speedup grows steadily for values
larger than 8 qubits.

###  A.3 Gradients and Hessians

The gradient of the $C_{xy}$ gate can then be computed analytically via

$\displaystyle\partialderivative{\phi}C_{xy}(\phi,\theta)=$ $\displaystyle
i\hat{C}\frac{1}{2^{n}}v_{\phi}(\hat{C})\odot\sum_{i=0}^{n}\exp(-i\lambda_{\text{XY,i}}\theta)\phantom{,}\hat{D}_{\lambda_{i}},$
(44) $\displaystyle\partialderivative{\theta}C_{xy}(\phi,\theta)=$
$\displaystyle\frac{1}{2^{n}}v_{\phi}(\hat{C})\odot\sum_{i=0}^{n}-i\lambda_{\text{XY,i}}\exp(-i\lambda_{\text{XY,i}}\theta)\phantom{,}\hat{D}_{\lambda_{i}}.$
(45)

For the derivative by $\phi$ the sum can be reused from the previous gate
construction in Eq. (41) further reducing computational demands. For the MS
gate we find

$\displaystyle\partialderivative{\phi}\text{MS}(\phi,\theta)=$ $\displaystyle
i\hat{C}\frac{1}{2^{n}}v_{\phi}(\hat{C})\odot\bigg{(}\sum_{i=0}^{\lceil\frac{n}{2}\rceil}\exp(-i\lambda_{\text{MS,i}}\theta)\phantom{,}\overbrace{\left(\hat{D}_{\lambda_{i}}+\hat{D}_{\lambda_{n-i}}\right)}^{=\tilde{D}_{i}}\delta^{\text{floor}}_{i}\bigg{)}$
(46) $\displaystyle\partialderivative{\theta}\text{MS}(\phi,\theta)=$
$\displaystyle\frac{1}{2^{n}}v_{\phi}(\hat{C})\odot\bigg{(}\sum_{i=0}^{\lceil\frac{n}{2}\rceil}-i\lambda_{\text{MS,i}}\exp(-i\lambda_{\text{MS,i}}\theta)\phantom{,}\overbrace{\left(\hat{D}_{\lambda_{i}}+\hat{D}_{\lambda_{n-i}}\right)}^{=\tilde{D}_{i}}\delta^{\text{floor}}_{i}\bigg{)}.$
(47)

Hessians can be derived direcly by applying the chain rule a second time in
Eqs. (44)-(47).

Figure 8: Average speedup (30 shots) required to compute the gradient for: (a)
the same circuit (50 gates) with increasing numbers of qubits. (b) a 6-qubit
(blue) and 7-qubit (green) gate set for increasing number of random gates on
the circuit. Shaded regions show the standard deviation for each point
(vertical bars are connected together to form a poligon). The orange line
represents the gradient compiled with NUMBA on CPU (Intel(R) Xeon(R) Gold 6240
CPU @ 2.60GHz), the second one the gradient compiled with JAX on GPU (NVIDIA
Tesla V100S-PCIE-32GB). We see a slowdown in the speedup as we move past 8
qubits: this is probably due to computational overheads that do not profit
from the GPU parallel calculations. In fact, we see that the slowdown is
present only as we increase the Hilbert space dimension and not as we increase
the number of gates.

##  Appendix B Parameter-shift rules for ion gates

The method for gate computation introduced above allows to express the
derivative of a quantum cost function in terms of so-called parameter-shifts
rules. These have been studied in the context of variational quantum circuit,
quantum control and quantum machine learning. The gates contained in the gate
set $Z_{1}(\theta),...,Z_{n}(\theta)$, $C_{xy}(\theta,\phi)$,
$\text{MS}(\theta,\phi)$. In general, expectation values with respect to
parametrized quantum circuits will have the form

$\displaystyle C(\boldsymbol{\alpha})=\bra{\psi}V(\boldsymbol{\alpha})\rho
V^{\dagger}(\boldsymbol{\alpha})\ket{\psi},$ (48)

for a given quantum state $\ket{\psi}$, where $V(\boldsymbol{\alpha})$ is a
gate in the gateset defined in Section (II). The cost function depends on a
specific target operator $\rho$.

We first show that the cost function in Eq. (48) is equivalent to the cost
function of Eq. (13), as well as to the expectation value of the local and
non-local Hilbert Schmidt circuit, which is given in Ref. [60]:

$\displaystyle C_{HST}(\boldsymbol{\alpha})=1-\Tr{H(U\otimes
V(\boldsymbol{\alpha})^{*})\rho(U^{\dagger}\otimes
V(\boldsymbol{\alpha})^{T})},$ (49)

with $\rho=H=\ket{\phi_{+}}_{A,B}\bra{\phi_{+}}_{A,B}$ for the Hilbert Schmidt
test, where $A$ and $B$ are the subystems for unitary $U$ and $V$ and

$\displaystyle\ket{\phi_{+}}_{A,B}=\frac{1}{\sqrt{d}}\sum_{\boldsymbol{j}}\ket{\boldsymbol{j}_{A}}\otimes\ket{\boldsymbol{j}_{B}}.$
(50)

###### Proof.

We show the equivalence for the Hilbert Schmidt test. The proof for the local
Hilbert Schmidt test is analogous.
We see that

$\displaystyle\Tr{H(UV^{\dagger}\otimes I)\rho(V\otimes
U^{*})}=\frac{1}{d}\sum_{\boldsymbol{j}}\Tr{\ket{\boldsymbol{j}_{A}}\bra{\boldsymbol{j}_{A}}UV^{\dagger}\ket{\boldsymbol{j}_{A}}\bra{\boldsymbol{j}_{A}}VU^{\dagger}}\Tr{\ket{\boldsymbol{j}_{B}}\otimes\bra{\boldsymbol{j}_{B}}},$

which can be represented as a sum of squared amplitudes of with the same
structure of (48). ∎

As a result, the cost functions considered here have to a common description,
similar to Eq. (48). We derive and comment the corresponding parameter-shift
rule for each one of these gates. We consider here the shift of one single
parameter at a time, i.e., a scalar rotation angle $\theta$. The gradient is
constructed by shifting each parameter according to its own specific
parameter-shift rule. The derivatives of the function can be written as:

$\displaystyle\partialderivative{\theta}C(\theta)=\bra{\psi}\partialderivative{\theta}V(\theta)UV^{\dagger}(\theta)+V(\theta)U\partialderivative{\theta}V^{\dagger}(\theta)\ket{\psi}.$
(51)

Assuming the gate has a generator $V(\theta)=e^{iG\theta}$, where $G$ is a
parameter-independent hermitian matrix, then we have:

$\displaystyle\partialderivative{\theta}C(\theta)=\bra{\psi}V(\theta)i[G,U]V^{\dagger}(\theta)\ket{\psi}.$
(52)

In the simplest case, we consider, e.g., $G=\sigma_{z}^{(i)}$ and obtain [86]:

$\displaystyle[\sigma_{z}^{(i)},U]=e^{i\sigma_{z}^{(i)}\pi/2}Ue^{-i\sigma_{z}^{(i)}\pi/2}-e^{-i\sigma_{z}^{(i)}\pi/2}Ue^{i\sigma_{z}^{(i)}\pi/2},$
(53)

which leads to

$\displaystyle\partialderivative{\theta}C(\theta)=C\left(\theta+\frac{\pi}{2}\right)-C\left(\theta-\frac{\pi}{2}\right),$
(54)

a parameter-shift rule which is valid for the local $Z$-rotations. For the
other two gates, we have to differentiate with respect to both the parameters
$\theta$ and $\phi$. For $\theta$, applying Eq. (15) in Ref. [64], we have:

$\displaystyle\left.\partialderivative{\theta}C(\theta,\phi)\right|_{\theta=0}=\sum_{l=1}^{2n}\frac{(-1)^{l-1}}{2\sin\left(\frac{2l-1}{2n}\pi\right)}C\left(\frac{2l-1}{2n}\pi,\phi\right).$
(55)

By using the general decomposition of the $C_{xy}$ gates, parameter-shift
rules can be obtained directly by studying how the eigenvalues and
eigenvectors of the unitaries vary as a function of the continuous parameters.
Since the $C_{xy}$ gate has $n+1$ distinct eigenvalues we will need at most
$2(n+1)$ samples to calculate the derivative – see [64]. We observe that the
matrices $\hat{V}_{n}(\phi)$ in Eq. (37) can be decomposed further in
elementary gates, such as:

$\displaystyle\hat{V}_{n}(\phi)=(\mathcal{P}(\phi)HX)^{\otimes
n}=\mathcal{P}(\phi)^{\otimes n}H^{\otimes n}X^{\otimes n},$ (56)

where

$\displaystyle\mathcal{P}(\phi)=\begin{pmatrix}1&0\\\
0&e^{i\phi}\end{pmatrix},H=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix},X=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}.$ (57)

Let us first consider the $C_{xy}$ gate. The corresponding Hamiltonian has
$n+1$ different eigenvalues with energy $\lambda_{i}=2i-n$ and
$c_{k}(\phi)=v_{k}(\phi)^{\dagger}Uv_{k}(\phi)$ are the overlaps between the
target operator $U$ and the gate eigenvectors $v_{k}(\phi)$ of the gate. For
the set of difference pairs we have $\lambda_{i}-\epsilon_{j}=2(i-j):=2l$.
Then Eq. (55) becomes

$\displaystyle\left.\partialderivative{\theta}C(\theta,\phi)\right|_{\theta=0}=\sum_{k=1}^{n}\sum_{l=1}^{2n}\frac{(-1)^{l-1}}{2\sin\left(\frac{2l-1}{2n}\pi\right)}c_{k}(\phi)e^{2ik(\frac{2l-1}{2n})}.$
(58)

.

This expression can be simplified by evaluating the sum at $\theta=0$ with
respect to the index $l$.

$\displaystyle\left.\partialderivative{\theta}C(\theta,\phi)\right|_{\theta=0}=\sum_{k=1}^{n}c_{k}(\phi)\frac{\csc\left(\frac{\pi}{2n}\right)e^{\frac{\pi
i(k(4n+2)-n)}{n}}\left(e^{\frac{2\pi
i(n+1)(n-2k)}{n}}-(-1)^{n+\frac{1}{n}}\right)}{2n\left((-1)^{n}e^{4\pi
ik}+e^{2\pi in}\right)}.$ (59)

For the MS gate we obtain instead the following expression:

$\displaystyle\left.\partialderivative{\theta}C(\theta,\phi)\right|_{\theta=0}=\sum_{k=1}^{n}c_{k}(\phi)\left\\{-\frac{\csc\left(\frac{\pi}{2n}\right)e^{8\pi
ikn-\frac{\pi i(n-2k)^{2}}{n}}\left((-1)^{n+\frac{1}{n}}e^{\frac{2\pi
i(n+1)(n-2k)^{2}}{n}}-1\right)}{2n\left((-1)^{n}e^{2\pi
i\left(4k^{2}+n^{2}\right)}+e^{8\pi ikn}\right)}\right\\}.$ (60)

For the derivative of Eq. (48) with respect to $\phi$, we just have to
consider the gate $\mathcal{P}(\phi)$ in Eq. (56), whose single-qubit gate
generator is given in Eq. (61). We observe that due to the simple structure of
the gate, the derivative of the $C_{x,y}$ gate with respect to $\phi$ is
simply given by:

$\displaystyle\mathcal{P}(\phi)=\exp{\frac{i\phi}{2}(I_{n}-S_{z})}=\exp{\frac{i\phi}{2}}\exp{-\frac{i\phi}{2}S_{z}},$
(61)

where $S_{z}=\sum_{i=1}^{n}\sigma_{z}^{(i)}$ is the collective $Z$-operator
acting on the entire qubit register. After inserting Eq. (61) into Eq. (48)
for a $C_{xy}$ gate, and using the representation in Eq. (35), we obtain

$\displaystyle
C(\theta,\phi)=\bra{\psi}\left(\exp{-\frac{i\phi}{2}S_{z}}B_{\text{XY}}(\theta)\exp{\frac{i\phi}{2}S_{z}}U\exp{-\frac{i\phi}{2}S_{z}}B_{\text{XY}}(\theta)^{\dagger}\exp{\frac{i\phi}{2}S_{z}}\right)\ket{\psi},$
(62)

where $B_{\text{XY}}(\theta)=H^{\otimes n}X^{\otimes
n}\left(\sum_{l=0}^{n}\exp{-i\lambda_{l}\theta}\ket{l}\bra{l}\right)X^{\otimes
n}H^{\otimes n}$ – the same equation holds for the MS gate, one needs just to
replace $\lambda_{\text{XY}}$ with the corresponding eigenvalues
$\lambda_{MS}$. Considering that $\mathcal{P}(\phi)$ has two eigenvalues, we
can write

$\displaystyle\exp{-\frac{i\phi}{2}S_{z}}=e^{i\phi}\Pi^{(1)}_{z}+e^{-i\phi}\Pi^{(2)}_{z},$
(63)

where $\Pi^{(1)}_{z}$ and $\Pi^{(2)}_{z}$ are the projectors on the respective
degenerate eigenspaces. By using Eq. (63) in Eq. (62), we see that Eq. (62)
has the following form

$\displaystyle
C(\theta,\phi)=b_{0}(\theta)+b_{1}(\theta)e^{-i\phi}+b_{2}(\theta)e^{i\phi}+b_{3}(\theta)e^{2i\phi}+b_{4}(\theta)e^{-2i\phi},$
(64)

where $b_{i}(\theta),i=0,1,2,3,4$ are $\phi$-independent coefficients. Thus,
we can write the derivative of $C$ with respect to $\phi$ as

$\displaystyle\partialderivative{\phi}C(\theta,\phi)=-ib_{1}(\theta)e^{-i\phi}+ib_{2}(\theta)e^{i\phi}+b_{3}(\theta)2ie^{2i\phi}-2ib_{4}(\theta)e^{-2i\phi},$
(65)

which can be written as a linear combination of
$C(\theta,\phi\pm\frac{\pi}{2})$ and $C(\theta,\phi\pm\frac{\pi}{4})$ cost
functions:

$\displaystyle\partialderivative{\phi}C(\theta,\phi)=\tilde{C}(\theta,\phi+\frac{\pi}{4})-\tilde{C}(\theta,\phi-\frac{\pi}{4})+(\frac{1}{\sqrt{2}}-\frac{1}{2})\tilde{C}(\theta,\phi-\frac{\pi}{2})+(\frac{1}{\sqrt{2}}+\frac{1}{2})\tilde{C}(\theta,\phi+\frac{\pi}{2}).$
(66)

##  Appendix C Algorithmic implementation details

Our framework consists of two main parts: The agent (we consider a PS-LSTM
agent with LSTM cells and a linear layer, but any RL agent is a viable
option), which should learn to construct a proper representation of the
quantum circuit and the optimizer, which has to be equipped with a proper
gradient function. The gradient function is constructed according to the
standard GRAPE procedure and using the gate representations for the $C_{xy}$
and the MS gates given in Appendix A. We use and test two different versions
of the gradient: The first one is compiled using NUMBA [82], a library for
fast python code, the second one employs JAX to allow execution on GPU. A
comparison of the two gradient functions is given in Fig. 8. We observe that
the GPU-based function allows for a certain speed-up, in particular as the
number of gates on the circuit increases. The main advantage of NUMBA lies in
a faster and more straightforward implementation of the parallel optimization
runs with different seeds.

At time $t$ of the agent-environment interaction, the agent receives an input
(percept) from the environment and outputs an action which is executed onto
the environment. The agent-environment interaction has the following
structure.

Action: The action of the agent corresponds to placing one of the available
gates in the gate set (entangling or non-entangling) onto the quantum circuit.
The gate is represented by an integer $a\in{1,...,|A|}$. The array of all the
integers chosen by the agent up to time $t$ forms the quantum circuit
structure at time $t$.

Percept: As input to the RL agent, we generally use a one-hot encoding of the
circuit. For a circuit of length $L$ and $|A|$ different gates, the percept
$s_{t}\in\\{0,1\\}^{L}\crossproduct\\{0,1\\}^{|A|}$ has entries equal to:

$\displaystyle(s_{t})_{i,j}=\begin{cases}1&\text{if $i,j=a,t$}\\\
0&\text{otherwise.}\end{cases}$ (67)

However, for the PS-LSTM agent (or other agents that are modified
analogously), a more suitable input can also be used. Since the update of the
internal state of the PS-LSTM agent follows the RL steps, the agent can just
accept the input at time $t$ as a percept, since previous information about
past agent-environment interactions is still processed by the the internal
state of the LSTM network. Therefore, the percept becomes the one-hot encoded
vector

$\displaystyle(s_{t})_{j}=\begin{cases}1&\text{if $j=a$}\\\
0&\text{otherwise.}\end{cases}$ (68)

This percept only gives the agent partial information about the internal state
of the environment, which turns the problem from an MDP (Markov Decision
Process) into a so-called POMDP problem (Partially Observable Markov Decision
Process) [40]. Other agents than those derived by PS may use different inputs,
based on their specific convergence properties in dealing with different
environments. In general, if the agent can accept a recurrent or
autoregressive network as a policy, the percept given in Eq. (68) may be used.
It could be, however, that the algorithm needs to be modified appropriately to
function with this type of percept. The first type of input has been used for
the simulation of UCC operators, whereas the second one has been implemented
in all the other simulations. In general, both percepts lead to similar
results, but the second one should be preferred for the PS-LSTM architecture
and most importantly does not scale with the size of the circuit, leading
therefore to faster forward passages in the policy network.

Reward and curriculum: Throughout the development of this manuscript, several
reward systems were tested.

Figure 9: Quantum circuit representation of two possible RL environments for
quantum circuit optimization. The first environment (a) resembles the
structure of the layer-based compilation, that is, the entangling MS gates are
fixed and the agent can place rotations on the circuit between the entangling
layers. The second one (b) has no pre-defined circuit structure but rather
leaves the agent completely free to place any gate withing the gate set on the
quantum circuit, with only one additional simplification: two gates of the
same type placed immediately one after the other are automatically merged to
form one single gate. This is done to prevent the agent from getting stuck in
loops, i.e., local minima, where it keeps choosing the same gate over and over
again. In general, the first circuit reduces the size of the action space and
therefore the possible shapes of the corresponding cost function, hence
reducing exploration in favour of a more standardized search.

In general, we used a two-step reward that assigns a smaller reward value when
the cost function minimum falls below the actual curriculum threshold
$\epsilon_{t}$ and a larger reward value when the cost function minimum falls
below the global target threshold $\epsilon_{\text{min}}$. Both curriculum
thresholds can be adjusted depending on the gate synthesis problem to be
tackled. The episode terminates when the cost function minimum in a given time
step falls below the threshold $\epsilon_{t}$ or when the maximal length of
the circuit per episode, $L_{\text{max}}$, is reached. The RL training
terminates upon reaching the maximum number of episodes $E_{\text{max}}$. The
reward scheme helps to progressively increase the fidelity throughout training
without allowing for too long circuits. The threshold is then lowered as
episodes progress based on previous rewards obtained by the agent. In our
implementation, we lower the threshold when it has been surpassed by the agent
at least 500 times using the following scheme:

$\displaystyle\epsilon_{t+1}=\begin{cases}\epsilon_{\text{min}}+\frac{1}{2}(\epsilon_{\text{min}}-\epsilon_{t})&\text{if
$\epsilon_{\text{min}}\leq\epsilon_{t}\leq 1$}\\\
\epsilon_{\text{min}}&\text{if $\epsilon_{t}\leq\epsilon_{\text{min}}$}\\\
1&\text{otherwise}\end{cases}$ (69)

Environments: We propose two different types of agent-circuit interaction. In
the first type – see Fig. 9 (a) –, which is more similar to layer-based
compilation, the agent only places rotations on the circuit, while the
entangling gates are fixed in position. While the episode progresses, more
entangling gates are added to the environment, until a maximum
$L_{\text{max}}$ of gates on the circuit is reached. The circuit is also
simplified automatically by the environment, in such a way that if the agent
places the same two gates on the circuit one after the other, they are reduced
to one single gate. This is implemented in order to prevent the agent from
getting trapped in local minima where it places the same gate over and over
again on the circuit, therefore reducing the necessary training time. In the
second type – see Fig. 9 (b) –, entangling gates are also given as a possible
action on the environment. This second type of environment leaves more room
for exploration, but it is also more challenging for the agent. The
simulations presented in this work were realized by using only the first type
of environment.

Agents: In order to test the effectiveness of our algorithm, we employ
different agents for testing. We consider the standard state-of-the-art PPO
algorithm [87] and the two versions of PS with deep energy-based neural
networks PS-DEBN and PS-LSTM. We observe that PPO performs slightly worse than
PS-LSTM, but better than PS-FNN, but is probably due to the non-recurrent
version of the PPO algorithm implemented. Due to the large action and percept
space, we did not include standard PS in the comparisons, since this would
require to store large $h$-matrices for each percept-action transition,
leading to slow computation and eventually memory overflow. This algorithm can
however be used in circuit synthesis problems with a small number of gates and
qubits.

Optimizers: As an unconstrained optimizer, in this work we only consider the
L-BFGS-B algorithm, as it is implemented in SCIPY [85]. The number of
iterations of the optimizers for each RL interaction is set to 100. Both
optimizers can run a given number of optimization attempts in parallel with
different seeds (a technique usually referred to as random restart), which
should prevent the optimizer from getting trapped in local minima.

Simulations In this work we consider three different groups of simulations:
the synthesis of a $3$-qubit Toffoli gate on a $3$-qubit and $4$-qubit
circuit, the synthesis of UCC (Unitary Coupled Cluster) operators and the
synthesis of the unitary of the XXZ Hamiltonian with varying Hamiltonian
parameters. In the case of the XXZ Hamiltonian and the Toffoli gate, the
simulation is realized on CPUs (72 Intel(R) Xeon(R) Gold 6240 CPU @ 2.60GHz),
whereas the simulation of the UCC operators is performed on GPUs (4 NVIDIA
A100 Tensor Core GPU with 40 GB). This means, the library used to compute cost
functions and gradients in order to optimize them is NUMBA for the CPU-based
simulations and JAX for GPU-based ones. The reason is that NUMBA allows for
faster parallelization with the JOBLIB library that proves difficult to
achieve with JAX, thereby enabling us a better exploration of the cost
function landscape for small numbers of qubits, whereas JAX is significantly
faster for larger numbers of qubits. In particular, we set the number of
optimization runs per RL iteration to 10 when we employ the NUMBA version,
whereas we use only 1 when we employ the code using JAX. The curriculum
threshold was kept to $\epsilon_{\text{min}}=10^{-2}$. We used two-layered
LSTM networks implementing the policy given in Eq. (23) and 128 neurons, a
training batch size of 64, a learning rate of 0.01, a curriculum update window
of 500. The agent is trained with a replay-memory upon finding a gate sequence
with infidelity lower than the curriculum threshold, and also every 100 agent-
environment interactions if the replay memory is large enough. The target
network is updated every 50 agent-environment interactions. The number of
iterations of the optimizer in the hybrid RL-continuous simulations is fixed
to 100. The temperature parameter $\beta$ of the softmax distibution that
parametrized the PS-LSTM policy – see Eq. (19) – is annealed from
$\beta=10^{-3}$ to $\beta=1$ according to a linear schedule based on the
number of episodes, in order to force the agent, towards the end of the
training, to consider shorter and shorter circuits, thereby reducing the
exploration. Other parameters vary based on the simulation considered. The
data and hyperparameters can be found in
https://github.com/franz3105/RL_Ion_gates.