# Wide Quantum Circuit Optimization with Topology Aware Synthesis

Mathias Weiden  Department of Electrical Engineering and Computer Science,
University of California, Berkeley  Justin Kalloor  Department of Electrical
Engineering and Computer Science, University of California, Berkeley  John
Kubiatowicz  Department of Electrical Engineering and Computer Science,
University of California, Berkeley  Ed Younis  Computational Research
Division, Lawrence Berkeley National Laboratory  Costin Iancu  Computational
Research Division, Lawrence Berkeley National Laboratory

# Wide Quantum Circuit Optimization with Topology Aware Synthesis

Mathias Weiden  Department of Electrical Engineering and Computer Science,
University of California, Berkeley  Justin Kalloor  Department of Electrical
Engineering and Computer Science, University of California, Berkeley  John
Kubiatowicz  Department of Electrical Engineering and Computer Science,
University of California, Berkeley  Ed Younis  Computational Research
Division, Lawrence Berkeley National Laboratory  Costin Iancu  Computational
Research Division, Lawrence Berkeley National Laboratory

###### Abstract

Unitary synthesis is an optimization technique that can achieve optimal multi-
qubit gate counts while mapping quantum circuits to restrictive qubit
topologies. Because synthesis algorithms are limited in scalability by their
exponentially growing run time and memory requirements, application to
circuits wider than 5 qubits requires divide-and-conquer partitioning of
circuits into smaller components. In this work, we will explore methods to
reduce the depth (program run time) and multi-qubit gate instruction count of
wide (16-100 qubit) mapped quantum circuits optimized with synthesis. Reducing
circuit depth and gate count directly impacts program performance and the
likelihood of successful execution for quantum circuits on parallel quantum
machines.

We present TopAS, a topology aware synthesis tool built with the _BQSKit_
framework that preconditions quantum circuits before mapping. Partitioned
subcircuits are optimized and fitted to sparse qubit subtopologies in a way
that balances the often opposing demands of synthesis and mapping algorithms.
This technique can be used to reduce the depth and gate count of wide quantum
circuits mapped to the sparse qubit topologies of Google and IBM. Compared to
large scale synthesis algorithms which focus on optimizing quantum circuits
after mapping, TopAS is able to reduce depth by an average of 35.2% and CNOT
gate count an average of 11.5% when targeting a 2D mesh topology. When
compared with traditional quantum compilers using peephole optimization and
mapping algorithms from the Qiskit or $t|ket\rangle$ toolkits, our approach is
able to provide significant improvements in performance, reducing CNOT counts
by 30.3% and depth by 38.2% on average.

##  I Introduction

Modern quantum machines are subject to high levels of environmental noise, are
difficult to control, and consist of only tens of qubits. Quantum error
correction and mitigation are widely regarded as vital components in the
development of practical quantum computers [16]. However, because quantum
error correction is not yet realistically implementable, it is essential that
quantum programs be optimized so that run times and operation counts are
minimized, as each operation adds additional error into the system and longer
run times increase the probability that qubits will decohere. Operations that
involve multiple qubits, such as the _CNOT_ , _iSWAP_ , _Toffoli_ , and _CZ_
gates, are typically far more error prone and expensive to implement compared
to single qubit rotation gates [6]. For this reason, multi-qubit gate counts
can be used as a metric for the error introduced into quantum programs.

Optimization techniques that aim to decrease the multi-qubit gate count of
quantum circuits are becoming increasingly used and studied [12, 25, 18].
Unitary synthesis is a method of decomposing a quantum circuit into a simpler
quantum circuit that fits some specified gate set. Synthesis seeks to produce
shorter circuits with reduced depth and gate count in order to improve the
likelihood that circuits are executed correctly.

The runtime of unitary synthesis algorithms scales exponentially with the
number of qubits $n$. For current state of the art unitary synthesis
algorithms such as the QSearch/LEAP algorithm, only circuits up to 5 qubits
are realistically decomposable. Larger circuits must thus be partitioned into
smaller subcircuits. After synthesizing each subcircuit independently,
subcircuits in the original circuit are replaced with their synthesized
versions.

Other than sensitivity to environmental noise, machines in the Noisy
Intermediate Scale Quantum (NISQ) era are defined by their limited sizes and
connectivities (which physical qubits are allowed to interact). The
connectivity of qubits in machines using superconducting qubits is especially
restrictive. Examples of realistic and popular superconducting qubit physical
topologies are illustrated in Figure 1.

(a) Mesh

(b) Linear

(c) Falcon

Figure 1: Example Physical Qubit Topologies. The 2D mesh topology is based off
quantum machines provided by Google, while the falcon topology follows IBM’s
heavy hexagonal connectivity scheme.

Each vertex in a physical qubit topology graph represents a hardware qubit.
The set of edges between these hardware qubits is the set of supported multi-
qubit interactions for that quantum machine. Quantum algorithms are typically
designed assuming that each qubit is able to interact directly with all other
qubits: logical circuits assume all to all connectivity. In order for a
quantum processor to execute these densely connected logical circuits, a qubit
mapping algorithm must insert SWAP gates routing all interactions along edges
in the physical topology.

As quantum circuits become wider and physical topologies remain sparse, the
number of routing gates inserted by these mapping algorithms grows quickly. It
is therefore desirable to approximate and transform wide quantum circuits
using unitary synthesis so that the total number of multi-qubit gates and
depth of quantum circuits is reduced.

In this paper, we present TopAS, a qubit topology aware synthesis tool that
applies unitary synthesis to wide quantum circuits to reduce multi-qubit gate
counts and circuit depth. TopAS first partitions a logical quantum circuit,
synthesizes each partitioned subcircuit independently, reassembles the
optimized logical circuit, then uses a mapping algorithm to ensure that all
operations are legal according to some specified qubit topology. When
targeting the mesh physical topology, TopAS is able to produce circuits with
an average of 30.1% fewer CNOT gates than _Qiskit_ [4] and $t|ket\rangle$
[20], and 11.5% fewer CNOT gates than other large scale synthesis techniques
such as QGo [23]. TopAS improves upon the QGo algorithm by optimizing before
mapping and by preconditioning circuit partitions so that the fully mapped
results are less deep and require fewer multi-qubit gates.

The remainder of this paper is structured as follows: Section II discusses the
general processes of unitary optimization, quantum circuit partitioning, and
compares tools that apply synthesis before and after mapping. Section III
presents the design choices made for the TopAS tool. Section IV compares the
TopAS tool to other optimization and mapping techniques. Finally, Section V
provides commentary and discussion about advancements that could further
improve the performance of wide quantum circuit synthesis tools.

##  II Background

###  II-A Quantum Computing Basics

The fundamental unit of information in a quantum computer is the qubit, which
can be represented as a vector of the form

$|\psi\rangle=\alpha\begin{bmatrix}1\\\
0\end{bmatrix}+\beta\begin{bmatrix}0\\\ 1\end{bmatrix}$

where $\alpha$ and $\beta$ are complex numbers such that
$|\alpha|^{2}+|\beta|^{2}=1$ and $[1\;0]^{T}$ and $[0\;1]^{T}$ are orthonormal
basis vectors representing two distinct quantum states. The state of a quantum
system with $n$ qubits in states
$|\psi_{1}\rangle,|\psi_{2}\rangle,\dots|\psi_{n}\rangle$ lies in a
$2^{n}\times 2^{n}$ Hilbert space, and can be represented by the tensor
product
$|\psi_{1}\rangle\otimes|\psi_{2}\rangle\otimes\dots\otimes|\psi_{n}\rangle$.
This quantum state $|\psi\rangle$ can be evolved by use of $2^{n}\times 2^{n}$
unitary operators[14].

The quantum circuit model represents this unitary as a series of quantum gates
[3]. Single qubit and multi-qubit gates act on qubits which are drawn as
horizontal wires (see Figure 3). The number of wires or qubits $n$ is called
the circuit width, while the critical path length or depth of the circuit is
$T$. The depth is typically drawn as the x-axis of the circuit. For the scope
of this paper, a universal gate set of {U3, CNOT} is assumed. Note that SWAP
gates can trivially be decomposed into 3 CNOTs as shown in Figure 2.

Figure 2: How to implement a SWAP operation using a SWAP gate, 3 CNOTs, and 3
CNOTs and 4 Hadamards.

###  II-B Quantum Circuit Synthesis

Given a target $2^{n}\times 2^{n}$ unitary $U$ and an error threshold
$\epsilon$, a unitary synthesis algorithm builds a new circuit whose unitary
$U_{S}$ satisfies the inequality $\|U-U_{S}\|\leq\epsilon$ [25]. Most recent
synthesis tools base their distance metric off the Hilbert-Schmidt inner
product as it is computationally inexpensive [2, 11, 9]:

$\|U-U_{S}\|_{HS}=Tr(U^{{\dagger}}U_{S})$

There are broadly two types of synthesis, top-down and bottom-up. Top-down
synthesis techniques are rule-based, and aim to break down large unitaries
into smaller ones. Usually, these algorithms are quick, but the resultant
circuit depth grows exponentially, which limits their effectiveness. Bottom-up
synthesis starts with an empty circuit and gradually adds gate until a
solution is found. Techniques such as QSearch/LEAP use an A* heuristic search
to find an optimal depth approximation for the overall unitary [21].

Furthermore, synthesis algorithms such as QSearch also accept a coupling graph
as input so that the connectivity between qubits may be specified. The
coupling graph provided restricts the synthesis algorithm by requiring that it
only place multi-qubit gates that obey this connectivity. Doing so means that
the resulting circuit is fully mapped to the specified coupling graph.
Previous work [21] has demonstrated that for small circuits, synthesis
algorithms are able to map circuits to restrictive qubit topologies using
fewer gates than other compiler tools such as _Qiskit_ [4] and $t|ket\rangle$
[20]. Thus when provided a coupling graph and a small width circuit, synthesis
is able to completely remove the need for a mapping algorithm. Synthesis
algorithms are ultimately limited by the size of the solution search trees,
and the size of the unitary matrices for which the Hilbert-Schmidt distance
must be calculated.

###  II-C Quantum Circuit Fidelity and Performance

In the NISQ era, the probability that a quantum circuit is successfully
executed depends heavily on the number of multi-qubit gates and the execution
time of circuits. This is because multi-qubit gates have a high probability of
introducing noise and because quantum states tend to decohere after a certain
amount of time [8]. Producing quantum circuit implementations with both fewer
multi-qubit gates and lower depth is therefore desirable, and signals that an
implementation has a higher chance of executing properly. For $N$ partitions
and a synthesis threshold $\epsilon=10^{-10}$ per partition, the total circuit
approximation error $N\epsilon$ is typically between $10^{-8}$-$10^{-7}$ in
Hilbert-Schmidt distance for our selection of benchmarks. This means although
synthesis introduces approximation errors into circuit implementations, the
error threshold is sufficiently low that gate and decoherence errors will
dominate overall error. For this reason, we present both CNOT gate counts and
circuit depth to evaluate circuit implementations (see Section IV).

Reducing depth is analogous to reducing the run time of a classical program,
and thus also acts as a direct measure of circuit performance assuming a
quantum machine is capable of executing gates in parallel. Mitigating other
sources of error, such as superconducting crosstalk [17], remains an important
area of research but is not considered in our evaluation.

###  II-D Mapping Quantum Circuits

Quantum circuits are typically designed with the assumption that the
underlying hardware that may run the circuit can support interactions between
any pair of qubits. However, due to technological restrictions modern quantum
machines typically have very sparse inter-physical qubit connectivity. Several
examples of realistic sparse physical qubit topologies are illustrated in
Figure 1.

In order to ensure that a quantum circuit can be run on a quantum computer, a
mapping algorithm is used to transform the circuit so that it conforms to some
restricted qubit connectivity. The process of mapping quantum circuits to
physical topologies happens in two stages: placement or layout and routing.
During the placement phase, an assignment of logical qubits in the quantum
circuit to physical qubits in the qubit topology is created. In the next
phase, the routing phase, quantum SWAP operations are inserted to ensure that
all multi-qubit interactions specified in the quantum circuit may take place
along edges in the physical topology [12, 18]. Throughout both phases of
circuit mapping, the goal is to minimize the number of SWAP operations
inserted into the final circuit. This is because for machines using the {CNOT,
U3} gate set, a SWAP is implemented using 3 CNOT gates (Figure 2). As quantum
circuits become wider and more densely connected, the overhead in CNOT count
due to routing SWAP gates quickly rises. However, minimizing the number of
SWAP operations is known to be NP-Hard [19], so most algorithms use heuristic-
based approaches that are sub-optimal, especially for wider circuits.

###  II-E Quantum Circuit Partitioning

Given a quantum circuit and a positive integer $k$ (called the partition width
or block size), a partitioning algorithm divides the circuit into subcircuits
(also called blocks or partitions) of width at most $k$. Valid partitions
consist only of gates that act on those specific qubits in the partition. As
soon as a gate crosses a partition boundary, the partition must be cut
vertically along the violating wire. A new partition can then be created, or
more gates can be added so long as they act only on qubits still within the
partition. An example circuit with width 3 partitions is shown in Figure 3.

Figure 3:  Example quantum circuit with 5 qubits. Each qubit is represented by
a horizontal line, and the timing of the circuit is described by the x-axis.
The circuit has been partitioned into 3 subcircuits.

For the purposes of unitary synthesis, the goal of a partitioning algorithm is
to form as few partitions as possible, and thus ensure that each of these
partitions is as large as possible. Our experiments show that these large
partitions tend to see a larger reduction in CNOTs as compared to smaller
partitions. This effect is illustrated in Figure 5.

###  II-F Post-Mapping vs. Logical Circuit Synthesis

Previous tools to apply unitary synthesis optimization to wide quantum
circuits, such as QGo [23], follow a _post-mapping_ synthesis flow as
illustrated in Figure 4(a). After preliminary circuit optimizations, the
logical circuit is mapped to a specified physical qubit topology. This mapped
quantum circuit is then partitioned into subcircuits, each of which are
independently synthesized. The synthesized subcircuits are then reassembled
into the complete optimized and mapped circuit.

Although post-mapping synthesis schemes are able to reduce depth and multi-
qubit gate count, there are several pitfalls that limit their effectiveness.
First, because mapped circuits typically contain more gates than their
unmapped counterparts, partitioning algorithms tend to form far more
partitions on mapped circuits. Because each of these partitions must be
individually synthesized, this approach is typically more time consuming than
partitioning and synthesizing the logical, unmapped, circuit. Second, these
methods are very sensitive to the quality of the mapping algorithm used, as
they can only reduce inefficiencies introduced into mapped circuits. Doing the
same on the logical circuit instead allows for mapping to be done on an (often
much) shorter depth and gate count approximate circuit. Logically synthesized
circuits typically have fewer gates, which means fewer SWAP gates are
typically needed during routing.

These points motivate the reasoning behind adopting a pre-mapping, or logical
circuit synthesis flow. The process of logical circuit synthesis for wide
quantum circuits is illustrated in Figure 4(b). After some initial quick
circuit optimizations, the logical circuit is partitioned into synthesizable
subcircuits. Each of these subcircuits is paired with a synthesis subtopology
to which the subcircuits are mapped. After each subcircuit is synthesized, the
circuit is reconstructed and finally fully mapped to the specified physical
qubit topology.

(a) Post-mapping synthesis optimization.

(b) Pre-mapping or logical circuit synthesis optimization.

Figure 4: Program flow diagrams for wide circuit optimization tools using
unitary synthesis post- and pre-mapping.

###  II-G Synthesis Subtopology Selection

Each partitioned subcircuit must be assigned a graph $G_{S}=(V,E_{S})$ that
specifies the connectivity between qubits in the partition. The synthesis
algorithm chooses multi-qubit gates that correspond to edges in the set
$E_{S}$, mapping the subcircuit to this _synthesis subtopology_ $G_{S}$. The
synthesis subtopology selection process in logical circuit synthesis tools
diverges from that in post-mapping synthesis tools such as QGo. Since post-
mapping tools assume circuits are already routed, the interactions between the
qubits in a partition are guaranteed to fit the physical topology. Therefore,
we can choose $G_{S}$ to be the subgraph induced by the physical qubits in a
partition.

Logical synthesis affords more choice here as the partitioned circuit is not
yet routed to obey the restrictions of the physical hardware. Each synthesis
subtopology is a graph of order equal to the number of qubits in that
partition. The only other requirement for synthesis subtopologies is that they
are connected. This is because if two qubits interact in the partitioned
circuit, but there is no path between them in the synthesis subtopology, the
synthesis algorithm will be unable to find a circuit whose unitary
approximates the original subcircuit’s unitary. As synthesis algorithms are
mostly limited to operating in the 3-5 qubit partition width ranges, the
number of possible synthesis subtopologies is relatively limited.

Our experiments show that the choice of synthesis subtopology can have drastic
effects on the number of multi-qubit gates in the optimized subcircuit, as
well as on the performance of the routing algorithm during the final mapping
process. The goal of the synthesis subtopology selection process is thus to
choose a graph $G_{S}$ for each partition such that the synthesis algorithm is
able to produce low gate count subcircuits and the mapping algorithm is able
to more efficiently place and route the subcircuit. The structure of the
subtopology $G_{S}$ also heavily impacts the depth of synthesized subcircuits.
Depth is heavily dependent on the degrees of vertices in the subtopology
$G_{S}$. Subtopology selection thus also must carefully consider the
connectivity within a partition to ensure that the depth of the optimized
circuits is minimized. Section III-C further details the choices made for the
subtopology selection algorithm for the TopAS tool.

##  III Topology Aware Synthesis

We present TopAS, a physical qubit topology aware synthesis tool. By
partitioning and synthesizing subcircuits from the logical circuit, the number
of multi-qubit operations in a quantum circuit is reduced before mapping
begins. The choice to synthesize logical quantum circuits enables TopAS to
produce mapped circuits with fewer total CNOT gates than many other
optimization tools. TopAS selects synthesis subtopologies for partitioned
subcircuits that are sparse and easily embedded within the underlying physical
qubit topology on which the circuit will be run. Synthesis subtopology
selection is done in such a way that reduces the number of operations needed
for both computation and mapping.

The TopAS tool uses a _scan_ partitioning strategy introduced by the authors
of QGo [23]. The synthesis algorithm used is the QSearch/LEAP algorithm [21].
Working versions of the partitioning and synthesis algorithms are provided by
the _BQSKit_ tool [24]. A flowchart of the TopAS tool’s execution is displayed
in Figure 4(b).

###  III-A Logical Quantum Circuit Partitioning

The partitioning algorithm views the input circuit as a two dimensional grid,
where qubits are represented by rows, and time steps in the circuit are
represented by columns. Each element of the grid either contains a quantum
operation or is empty. The width of this grid is $n$, the number of qubits in
the circuit. The length is $T$, the depth or critical path length of the
circuit. The _scan_ partitioner algorithm sequentially _scans_ through all as
yet unpartitioned gates in a circuit, examines each possible grouping of $k$
or fewer qubits, and picks the one that lends the largest partition. The
_scan_ partitioner tends to form few partitions, and often forms partitions
with a high average number of multi-qubit gates. As mentioned in Section II-E
and illustrated in Figure 5, subcircuits with a large number of multi-qubit
gates tend to yield the most reduction by synthesis. The distribution of
partitions produced by the _scan_ partitioning algorithm for a subset of
benchmarks is also shown in Figure 5.

Note that partitions consist only of subcircuits with interacting qubits. A
circuit logical connectivity graph is used to enforce this policy. This graph
has a vertex for each qubit in the circuit. The presence of an edge $(u,v)$
indicates that there is a multi-qubit gate between qubits $u$ and $v$ at some
point in the circuit. Given a partition width of $k$, candidate partitions
consist of all connected subgraphs in this circuit logical connectivity graph
of order at most $k$. In the worst case, a circuit where all qubits interact
with each other, there are $O{n\choose k}$ candidate partitions.

Figure 5:  Partition frequency and mean normalized partition size after
optimization as a function of CNOT count. Normalized partition size is
calculated using the number of CNOTs in the synthesized subcircuit divided by
CNOTs in the original partition. A value less than 1 indicates that synthesis
reduces those partitions. The partitioner predominately forms subcircuits with
6-11 CNOTs, but is also able to form many subcircuits with 16-17 CNOTs.

###  III-B Synthesizing Logical Quantum Circuits

As logical circuits do not contain gates needed to conform to some restrictive
physical topology, they tend to be shorter than their mapped counterparts.
Shorter circuits tend to yield fewer partitioned subcircuits, as there are
fewer total gates to partition. Partitioning circuits before mapping thus
conserves the number of subcircuits produced by the partitioning algorithm. As
fewer subcircuits must be independently synthesized in this case, less time is
often required to optimize circuits that are partitioned before mapping.
Circuits with fewer partitions also tend to have more accurate output
circuits. Each partitioned circuit is synthesized independently, so each
synthesis procedure will produce a mapped circuit whose unitary representation
will be within some Hilbert-Schmidt norm distance $\epsilon$ away from the
target unitary. The authors of [15] showed that the total error in a circuit
that consists of synthesized subcircuits is bounded by the sum of each
subcircuit’s distance. If a circuit is partitioned into $N$ subcircuits,
because each partition is synthesized to within a distance $\epsilon$, the
total distance of the partitioned and synthesized circuit is bounded by
$N\epsilon$. Decreasing $N$, as is done in the logical partitioning case, thus
improves the accuracy of the synthesized circuit. Results comparing the upper
bound on total circuit errors for several synthesized benchmarks are shown in
Table I.

###  III-C Topology Aware Subtopology Selection

The main contribution of this work is a subtopology selection strategy and
compilation workflow built around it. It allows for partitioned subcircuits to
be preconditioned in such a way that balances the opposing demands of
synthesis and mapping. As discussed in Section II-G, the goal of the
subtopology selection process is to assign a graph $G_{S}$ to each partitioned
subcircuit. The vertices in $G_{S}$ represent qubits in the subcircuit, while
edges describe the allowed interactions between qubits in the output
synthesized subcircuit.

Figure 6:  The number of SWAP operations within partitions synthesized to
different synthesis subtopologies for the mesh, falcon, and linear physical
topologies. Simpler subtopologies result in fewer SWAP, but may produce
circuits with more CNOTs.

(a) Line

(b) Star

(c) Ring

(d) Kite

(e) Theta

(f) Complete

Figure 7:  Possible order 4 synthesis subtopologies. The graphs a, b, and c
are embedded within the mesh physical topology. Graphs d, e, and f are not.
Only a and b are embedded in the falcon topology, while a is the only embedded
subgraph of the linear topology.  Figure 8: Neighbor aware subtopology
selection encourages the reuse of edges between sequentially executing
partitions.

Ideally, a quantum circuit is synthesized so that it uses only edges present
in the physical topology for multi-qubit interactions. Before mapping, the
placement and thus allowed interactions between qubits is not known. However,
synthesizing densely connected subcircuits to sparse graph structures can
decrease the need for SWAP operations during the mapping phase. Limiting the
set of synthesis subtopologies and carefully selecting sparse subtopologies
for each subcircuit allows us to balance the opposing demands of the routing
and synthesis algorithms, leading to fewer total CNOT gates. Figure 6 shows
the average number of internal SWAP gates per CNOT for a variety of synthesis
subtopologies targeting three physical topologies. Internal SWAP gates are
routing operations that are inserted between the first and last gate of a
partition. If each subcircuit is synthesized to an embedded subgraph of the
physical topology and is executed atomically, there is a solution to the
mapping problem that requires zero SWAP operations on the subcircuit’s qubits
during its execution. However, because routing often disrupts the execution of
subcircuits, it is unlikely for there to be exactly zero internal SWAP
operations even for subtopologies that are embedded within the physical
topology. The number of internal SWAP gates observed for a partition is
normalized by the number of CNOTs in the partition. This allows us to more
effectively compare subcircuits of various sizes. Typically, sparser
subtopologies require more gates for computation than denser subtopologies,
but they require fewer SWAP gates during mapping.

The TopAS tool restricts the set of graphs from which to choose synthesis
subtopologies to the set of connected graphs that are embedded in the target
physical topology. TopAS uses a partition width of 4 qubits. The order 4
candidate subtopologies are illustrated in Figure 7.

QSearch, the synthesis algorithm at the core of the TopAS optimization flow,
aims to synthesize circuits using as few multi-qubit gates as possible. Often
times, circuit depth is minimized as a consequence of this goal but is not the
primary metric of success. Subtopology selection is the only variable that
TopAS uses to reduce circuit depth directly. For certain subtopologies
choices, circuits produced tend to contain many gates that must be executed
sequentially. The most significant contributor to this property is the
possible graph matchings present in a given subtopology. In topologies such as
the star (Figure 7b), any choice of a single edge is a maximum matching. This
means that for the star subtopology, no two CNOTs can be executed in parallel.
Other subtopologies like the line and ring (Figure 7a and 7c) have larger
maximum matchings, and thus allow more parallelism within synthesized
subcircuits. Although they typically require more routing gates, subtopologies
containing more edges tend to allow for lower depth circuits.

The qubit interactions within a partitioned subcircuit can be described using
a weighted undirected graph $G_{L}=(V,E_{L})$. For each multi-qubit
interaction that occurs between qubits $k_{1},k_{2}\in[K]$ in the subcircuit,
there is an edge $(k_{1},k_{2},w)\in E_{L}$, where the weight $w$ corresponds
to the number of times that interaction occurs. Synthesis algorithms such as
QSearch allow for the connectivity between qubits to be specified by an
unweighted undirected graph $G_{S}=(V,E_{S})$. For a maximum partition size of
$k$, $V=[k]$ and $E_{S}\subseteq V\times V$.

Our experiments indicate that synthesis is best able to reduce the multi-qubit
gate count of subcircuits when the logical connectivities and synthesis
subtopologies match. A kernel function [7] is used to quantify the similarity
between the graphs. The scoring function $K:G_{L}\times G_{S}\rightarrow[0,1]$
examines each edge in the synthesis subtopology $G_{S}$ and logical
connectivity $G_{L}$. From the edges, two vectors $v_{L},v_{P}\in
R^{k(k-1)/2}$ are constructed. The vector $v_{P}$ is simply an indicator
vector, with a $1$ at each element that corresponds to a present edge in
$G_{S}$. Element $i$ of $v_{L}$ contains the weight $w$ associated with edge
$i$ in $G_{L}$. The normalized inner product

$\text{similarity}(v_{P},v_{L})=\frac{v_{P}^{T}v_{L}}{\sum\limits_{i}v_{L}(i)}$

is then returned to quantify the similarity between the graphs. For each
partition, all $k!$ permutations of qubit labels are evaluated for each of the
candidate synthesis subtopologies. The permuted subtopology with the highest
kernel function score is considered the best candidate subtopology for the
purposes of producing the smallest output subcircuit. The effectiveness of the
similarity kernel function was evaluated by synthesizing partitions of width 4
to all order 4 subgraphs embedded within the mesh physical topology (Figures
7a, 7b, and 7c). In this scenario, the similarity function was able to
identify the subtopology that would produce the synthesized subcircuits with
the fewest CNOT gates 89% of the time.

If multiple subtopologies have the same similarity score, that with the fewest
edges is preferred. In order favor sparse subtopologies, each similarity is
multiplied by a bias factor. The biases that produced the best results for the
mesh topology were $1.0$ (line), $1.0$ (star), and $0.8$ (ring).

When selecting subtopologies, TopAS also considers the impact of subtopology
choice for the immediately preceding and succeeding partitions. Although
partitions are synthesized separately, they are not executed in isolation. If
a partition assumes a physical edge exists between two qubits, and the
partition immediately following also assumes this physical edge exists, fewer
total SWAPs may be needed between the execution of these two partitions. The
TopAS tool therefore checks each subcircuit’s neighboring partitions for
qubits that are shared. If interactions from the neighboring partitions occur
between shared qubits, they are added as edges to the current partition’s
logical connectivity graph $G_{L}$. Interactions from neighboring partitions
carry the same weight as interactions within a partition, even if the qubits
in the neighbor partitions do not interact directly within the partition
itself. This policy was chosen as larger amounts of edge sharing between
partitions greatly reduced the final CNOT count of optimized circuits. TopAS
uses this _neighbor aware_ version of the logical connectivity graph in the
computation of the similarity function. An example of two sequentially
executing partitions with subtopologies that have been selected to reuse an
edge is shown in Figure 8.

Algorithm 1 Subtopology Selection

Input: graph $G_{L}$, set of graphs $G_{neighbors}$, $bias$ function
Output: graph $G_{S}$

for all graphs $G_{n}$ in $G_{neighbors}$ do

if $shared\\_edge$ in $G_{n}$ and in $G_{L}$ then

$G_{overlap}\leftarrow shared\\_edge$

end if

end for

$v_{L}\leftarrow$ edge weight vector for $G_{overlap}\cup G_{L}$

$candidates\leftarrow$ order $\leq 4$ subgraphs in physical topology

for all graphs $G_{P}$ in $candidates$ do

$v_{P}\leftarrow$ indicator vector for $G_{P}$

$score=$similarity$(v_{P},v_{L})$

if $bias(G_{P})\times score>high\\_score$ then

$high\\_score\leftarrow bias(G_{P})\times score$

$G_{S}\leftarrow G_{P}$

end if

end for

return $G_{S}$

###  III-D Partition Replacement and Mapping

As indicated in Figure 5, it is sometimes the case that synthesized
subcircuits grow instead of reduce in size. When this happens, TopAS considers
replacing the synthesized subcircuit with the original subcircuit. Always
choosing the logical subcircuit when it has fewer CNOTs can be problematic in
the final routing pass, as the original subcircuit may be more difficult to
route than the synthesized version. Thus, TopAS considers differences in both
the number of CNOTs and logical connectivities of the two subcircuits. TopAS
adopts a policy of replacing the synthesized subcircuit with the original if
it shows more than 30% (empirically determined best threshold) fewer CNOTs or
has a more routable logical connectivity. The routability of a subcircuit’s
logical connectivity is determined by the average number of of internal SWAP
operations per CNOT for a given physical topology (see Figure 6).

After subtopology selection and synthesis take place, TopAS maps the
synthesized quantum circuit to a physical topology. For placement and routing,
the SABRE Layout and SABRE Swap [12] mapping algorithms are used.

##  IV Results

|  |  | Total | CNOTs per
---|---|---|---|---
| Mapping | Topology | Error | Partition
| QGo | Falcon | 1.38E-08 | 15.3
mult 16 | QGo | Mesh | 1.51E-08 | 12.9
| TopAS | - | 7.20E-09 | 15.4
| QGo | Falcon | 1.17E-07 | 12.5
mult_32 | QGo | Mesh | 1.15E-07 | 11.7
| TopAS | - | 4.96E-08 | 15.1
| QGo | Falcon | 5.83E-08 | 7.3
qft_64 | QGo | Mesh | 5.18E-08 | 7.8
| TopAS | - | 2.32E-08 | 8.1
| QGo | Falcon | 8.09E-08 | 8.1
qft_100 | QGo | Mesh | 8.72E-08 | 7.4
| TopAS | - | 3.85E-08 | 8.0
| QGo | Falcon | 9.09E-08 | 6.5
add_65 | QGo | Mesh | 7.02E-08 | 7.3
| TopAS |  | 3.53E-08 | 7.0
| QGo | Falcon | 1.26E-07 | 11.8
add_101 | QGo | Mesh | 1.23E-07 | 9.4
| TopAS | - | 6.18E-08 | 7.0
| QGo | Falcon | 1.42E-07 | 7.61
tfim_40 | QGo | Mesh | 1.57E-07 | 6.3
| TopAS | - | 1.05E-07 | 8.1
| QGo | Falcon | 5.17E-07 | 8.4
tfim_100 | QGo | Mesh | 4.71E-07 | 7.7
| TopAS | - | 2.70E-07 | 8.0
| QGo | Falcon | 1.83E-07 | 8.0
hubbard_18 | QGo | Mesh | 1.60E-07 | 6.5
| TopAS |  | 9.56E-08 | 3.7
| QGo | Falcon | 3.86E-07 | 10.3
shor_26 | QGo | Mesh | 3.32E-07 | 11.6
| TopAS | - | 1.52E-07 | 13.9

TABLE I:  Upper bound of total circuit error (sum of Hilbert-Schmidt distances
per partition) and average CNOTs per partition for each benchmark. A synthesis
threshold of $\epsilon=10^{-10}$ was used for both QGo and TopAS.

The performance of TopAS is compared to other optimization and mapping tools
by measuring the depth and total CNOT operation counts of optimized and mapped
circuits. Here, we compare TopAS (as described in Section III) to the _Qiskit_
, $t|ket\rangle$, and QGo tools. In all cases, circuits are optimized using
$t|ket\rangle$’s full peephole optimization. Afterwards, either _Qiskit_ or
$t|ket\rangle$ are used to map the quantum circuit to the specified physical
topology. The QGo tool was used to optimize either the _Qiskit_ or
$t|ket\rangle$ mapped circuits, whichever contained fewer CNOT operations. The
synthesis tool used is the QSearch/LEAP algorithm [21]. The original QGo
algorithm was re-implemented using the partitioning and synthesis algorithms
provided by the _BQSKit_ toolkit [24].

Figure 9:  Comparison of relative CNOT count (top) and depth (bottom) for
circuits mapped to the Google style mesh physical topology. CNOT count and
depth is shown relative to optimized circuits mapped using _Qiskit_ ’s SABRE
Swap algorithm. Lower depth corresponds directly to better program runtimes on
machines with parallel gate execution, lower depth and CNOT count correspond
to improved execution probability on noisy machines.

Figure 10:  Comparison of relative CNOT count (top) and depth (bottom) for
circuits mapped to the IBM style falcon physical topology. CNOT count and
depth is shown relative to optimized circuits mapped using _Qiskit_ ’s SABRE
Swap algorithm. Lower depth corresponds directly to better program runtimes on
machines with parallel gate execution, lower depth and CNOT count correspond
to improved execution probability on noisy machines.

Performance was evaluated using a variety of wide quantum circuit benchmarks.
The SupermarQ benchmark suite [22] proposed a _volume_ metric to describe a
benchmark suite’s diversity. Our selection of benchmarks has a SupermarQ
volume of $7.46\times 10^{-6}$. This small selection of circuits outperforms
many other benchmark suites in SupermarQ volume due to its focus on real
algorithms. Our suite’s SupermarQ volume score is limited by the fact that
none of our benchmarks include intermediate measurements. This decision was
made as the primary focus in the evaluation of synthesis based optimization
tools is purely depth and multi-qubit gate count reduction. Each benchmark
assumes qubits are measured after the final program gate is executed.

The _qft_ and _shor_ circuits were generated using _Qiskit_. The _add_ and
_mult_ circuits were generated with [10]. The _hubbard_ benchmark was
generated with help from [13]. TFIM circuits were generated using the ArQtic
tool [1]. Benchmarks are labeled so that the number following the benchmark
name indicates the number of qubits in the circuit. These benchmarks were
selected as they were found to be relatively easily generated, and represent a
variety of wide circuits that may soon be executable on quantum machines.
Results for the QGo and TopAS tools were collected using NERSC’s _Perlmutter_
supercomputer.

Figures 9 and 10 show the relative CNOT gate count and depth of circuits
optimized and mapped to the mesh and falcon physical topologies illustrated in
Figure 1. The CNOT gate count and depth are normalized by the results obtained
for circuits optimized then mapped using _Qiskit_ ’s SABRE mapping algorithms.
A synthesis threshold distance of $\epsilon=1\times 10^{-10}$ and a partition
width of $k=4$ was used for both QGo and TopAS.

On average, the TopAS tool is able to reduce CNOT count compared to _Qiskit_
and $t|ket\rangle$ by 35.9% and 24.3% respectively when targeting a mesh
physical topology. For the falcon physical topology, TopAS outperforms
_Qiskit_ and $t|ket\rangle$ by an average of 33.4% and 19.1% respectively. For
the _tfim_100_ circuit, TopAS outperforms _Qiskit_ by 62.1% and $t|ket\rangle$
by 53.1%. For the _shor_26_ benchmark, TopAS outperforms _Qiskit_ by 38.7% and
$t|ket\rangle$ by 30.8%. Reductions in circuit depth are even greater,
averaging 38.6% and 37.6% lower compared to _Qiskit_ and $t|ket\rangle$ across
all benchmarks.

In most cases, the TopAS tool is able to produce circuits with fewer CNOTs and
much lower depth than the QGo tool, especially for the mesh physical topology.
On average, TopAS outperforms QGo by 11.5% in CNOT count and 34.3% in depth
when targeting the mesh physical topology. For the falcon physical topology,
TopAS only outperforms QGo by an average of 0.8% in CNOT count, but maintains
a 37.3% average depth reduction.

QGo optimized circuits always reduce CNOT count compared to unsynthesized
input circuits, but sometimes lead to deeper circuits. This is because when
synthesizing pure SWAP gates, the QSearch algorithm tends to find the more
expensive depth 5 implementation illustrated in Figure 2. TopAS is able to
avoid this pitfall by synthesizing partitions that do not contain SWAP gates,
thus resulting in much reduced circuit depths.

As discussed in Section III-B, working on the logical, unmapped, quantum
circuit also allows TopAS to form fewer partitions compared to QGo. This
results lower upper bounds on the amount of synthesis error for each
benchmark. Synthesis error upper bounds are given in Table I. This table also
lists the average number of CNOTs per partition for all benchmarks.

TopAS is able to outperform QGo in CNOT count for all benchmarks when
targeting the mesh physical topology except for the _hubbard_18_ circuit
benchmark. The discrepancy in performance for this benchmark is primarily
explained by the average size of partitions formed. Table I shows that for
this circuit, TopAS is only able to form partitions with an average of 3.7
CNOT gates. Figure 5 illustrates that typically, partitions with this number
of CNOTs tend not to reduce. Partitions formed by QGo tend to include SWAP
gates, which greatly increases the number of CNOTs per partition (6.5 CNOT
gates on average for _hubbard_18_). For the _hubbard_18_ circuit, TopAS’s
logical partitioning was simply not able to form few enough partitions to
compensate for the small average number of CNOT gates per partition. Using a
larger partition width, or a partitioning algorithm that is able to produce
larger partitions on average would likely improve the performance of TopAS for
the _hubbard_ benchmark. This point also partially explains the drop in TopAS’
performance compared to QGo for benchmarks targeting the falcon physical
topology.

Figure 11:  Number of CNOTs in QFT (top) and TFIM (bottom) circuits mapped to
the 2D mesh physical topology as a function of circuit width.

The scalability of the TopAS tool compared to _Qiskit_ , $t|ket\rangle$, and
QGo is illustrated in Figure 11. Each tool was used to optimize and map QFT
and TFIM circuits with widths of 10-100 qubits. Despite performing poorly in
comparison at small circuit widths, TopAS maintains an advantage compared to
other optimization tools as the circuit width increases. In the case of
certain circuits such as the TFIM circuits, TopAS’ advantage increases with
circuit width.

##  V Discussion

TopAS is able to produce circuits with fewer multi-qubit gates and lower depth
by partitioning logical quantum circuits and matching subcircuits with sparse
qubit subtopologies in such a way that balances the demands of synthesis and
mapping algorithms. Partitioning and using unitary synthesis to optimize the
logical quantum circuit allows for smaller circuits to be passed to mapping
algorithms. In the case of the _tfim_40_ and _tfim_100_ benchmarks, the
logical connectivity graphs are simply linear path graphs of order 40 and 100
respectively. These circuits have implementations that do not require any SWAP
operations, as they are directly embeddable within the mesh physical topology.
However, mapping algorithms implemented in _Qiskit_ and $t|ket\rangle$ are
unable to find these placements, likely due to the large widths and depths of
these circuits. By optimizing the logical quantum circuits, TopAS is able
shorten circuits in such a way that allows for these mapping algorithms to
find better placements. This advantage in mapping performance is maintained in
the falcon topology case, where the TFIM circuits are no longer directly
embeddable in the physical topology.

Using synthesis to map partitions to restrictive subtopologies also improves
the performance of the mapping algorithm. Mapping to subtopologies that are
easily embedded within the physical topology means that fewer SWAPs are needed
within the execution time of partitions. This effect is illustrated in Figure
6. The _neighbor aware_ subtopology selection mechanism illustrated in Figure
8 helps to maximize the reuse of edges in the physical topology between the
execution of partitions. The combination of these design choices allows for
TopAS to outperform other optimization and mapping tools.

Section IV demonstrates how TopAS is able to reduce both CNOT count and
circuit depth significantly compared to other tools. These metrics both play a
major role in the likelihood that a circuit will be executed correctly. The
amount of error due to synthesis (shown in Table I) is far lower than that
introduced by gate noise and decoherence. In total, optimizing wide circuits
with TopAS therefore greatly increases circuit fidelity compared to other
tools.

Because mapping is done after partitioning and synthesis in the TopAS program
flow, it is not limited to using a single mapping algorithm. We observed that
although TopAS synthesizes partitions to subtopologies embedded within the
physical topology, it is often the case that the mapping algorithm disrupts
the execution of synthesized subcircuits. When partitions are executed
atomically, there is a solution to the routing problem such that no SWAP gates
are needed during the partition’s execution. A mapping algorithm that is aware
of higher level structures than the primitive gate set may therefore further
improve the performance of TopAS. Such a tool was put forth by the authors of
[5], but it is only able to consider Toffoli gates instead of arbitrary $k$
qubit unitary operations. A full _partition aware_ mapping algorithm is
therefore likely necessary to reap the full benefits of this strategy.

Although TopAS’ set of valid synthesis subtopologies only includes graphs that
are embedded within the physical topology, it does not consider whether a set
of subtopologies can be packed into the physical topology simultaneously. This
fact explains the discrepancy in TopAS’ performance between the mesh and the
falcon physical topologies. For example, although the star subtopology is
embedded within the falcon physical topology, far fewer star graphs can be
packed into the falcon topology than the mesh topology. A subtopology
selection process that more precisely weighs the frequency with which
subgraphs appear in a physical topology would thus likely further improve the
performance of TopAS.

The two primary factors that limit the scalability of TopAS are the runtime of
the partitioning algorithm used and the runtime of the synthesis algorithm
used. In the case where qubits in the logical circuit interact with all other
qubits, the runtime of the _scan_ partitioner is $O{n\choose k}$, where $n$ is
the circuit width and $k$ is the partition width. This poor scaling limited
the width of benchmark circuits tested to 100 qubits. Other partitioning
schemes capable of forming large partitions on average are therefore needed to
optimize larger quantum circuits with TopAS. Algorithms such as QFAST [25]
effectively increase the width of circuits that can be optimized, but in our
experience tend not to perform as well as QSearch/LEAP in the range of 1-25
CNOT gates. Increasing the partition width tends to increase the number of
CNOTs per partition, which may provide more opportunity for improvement.
Increasing the partition width also has the effect of adding more graphs to
the set of valid synthesis subtopologies. Supposing that the optimization
algorithm scales, this may further improve the performance of TopAS by
allowing for a greater amount of mapping to be handled by the optimization
algorithm.

As shown in Figure 11, the number of CNOTs in QFT circuits optimized with the
TopAS tool grows roughly linearly with the number of qubits in the range
shown. As the QFT is a common component of more complex quantum algorithms,
the scaling of this circuit is of particular importance. With a width of 100
qubits, the TopAS optimized QFT circuit contains 23.9% and 11.7% fewer CNOTs
than the _Qiskit_ and $t|ket\rangle$ optimized and mapped circuits. Notably,
the rate of increase in CNOTs per qubit in the QFT circuit seems to grow more
slowly for TopAS than for _Qiskit_ and $t|ket\rangle$. Because QGo only
optimizes circuits that have already been mapped by these tools, we expect
that TopAS’ performance advantage over QGo can improve at larger (100-1000
qubit) circuit widths. This effect can clearly be seen in the TFIM circuit
scaling results.

Although out of the scope for the current paper, an obvious next question
would be to see what happens when TopAS is combined with QGo. We have clearly
shown that a single application of TopAS’ logical circuit synthesis almost
always outperforms a single round of QGo optimization. When applying both
TopAS and QGo to the _hubbard_18_ benchmark mapped to the 2D mesh physical
topology (a benchmark which TopAS fails to reduce CNOT count more than QGo for
the mesh topology), CNOT count was reduced by approximately 5% compared to QGo
alone. Thus, combining multiple levels of optimization has the potential to
further improve the performance of quantum circuits.

##  VI Conclusion

In this work, we have presented TopAS, a topology aware synthesis tool that
optimizes wide quantum circuits. By optimizing quantum circuits using unitary
synthesis before they are mapped to restrictive qubit topologies, TopAS
preconditions circuits so that they are minimized by synthesis and made easier
to route by mapping algorithms. TopAS is able to reduce CNOT count and circuit
depth, and thus improves circuit performance, compared to other state of the
art synthesis based optimization tools targeting wide circuit optimization.
TopAS also outperforms the optimization and mapping frameworks provided by
_Qiskit_ and $t|ket\rangle$. Because CNOT count and circuit depth are
significantly reduced, the likelihood that circuits optimized with TopAS
execute successfully on noisy, realistic, near term qubit topologies is
greatly increased. Further reduction is possible on TopAS optimized circuits
by applying successive rounds of post-mapping circuit optimization.

## Acknowledgements

This work was supported by the DOE under contract DE-5AC02-05CH11231 through
the Office of Advanced Scientific Computing Research (ASCR) Quantum Algorithms
Team and Accelerated Research in Quantum Computing programs, and by the NSF
Challenge Institute for Quantum Computation (CIQC) program under award
OMA-2016245.

## References

  * [1] L. Bassman, C. Powers, and W. A. de Jong, “ArQTiC: A full-stack software package for simulating materials on quantum computers,” 2021. [Online]. Available: http://arxiv.org/abs/2106.04749
  * [2] 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)_. IEEE, 2020, pp. 223–234.
  * [3] D. E. Deutsch, “Quantum Computational Networks,” _Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences_ , vol. 425, no. 1868, pp. 73–90, 1989.
  * [4] Q. Developers, “Qiskit: An Open-source Framework for Quantum Computing,” Jan. 2019. [Online]. Available: https://doi.org/10.5281/zenodo.2562111
  * [5] C. Duckering, J. M. Baker, A. Litteken, and F. T. Chong, “Orchestrated trios: compiling for efficient communication in quantum programs with 3-qubit gates,” in _Proceedings of the 26th ACM International Conference on Architectural Support for Programming Languages and Operating Systems_. ACM, apr 2021. [Online]. Available: https://doi.org/10.1145%2F3445814.3446718
  * [6] J. Gambetta and S. Sheldon, “Cramming more power into a quantum device,” _IBM Research Blog_ , Mar. 2019. [Online]. Available: https://www.ibm.com/blogs/research/2019/03/power-quantum-device/
  * [7] T. Hofmann, B. Schölkopf, and A. J. Smola, “Kernel methods in machine learning,” vol. 36, no. 3, 2008. [Online]. Available: http://arxiv.org/abs/math/0701907
  * [8] Kandala _et al._ , “Error Mitigation Extends the Computational Reach of a Noisy Quantum Processor,” _Nature_ , vol. 567, no. 7749, 2019.
  * [9] S. Khatri, R. LaRose, A. Poremba, L. Cincio, A. T. Sornborger, and P. J. Coles, “Quantum-Assisted Quantum Compiling,” _Quantum_ , vol. 3, p. 140, 2019\.
  * [10] H. Khetawat, M. Hassan, A. Neri, A. Rodrigues, and T. Wong, “QArithmetic.” [Online]. Available: https://github.com/hkhetawat/QArithmetic
  * [11] V. Kliuchnikov, A. Bocharov, and K. M. Svore, “Asymptotically Optimal Topological Quantum Compiling,” _Physical review letters_ , vol. 112, no. 14, p. 140504, 2014.
  * [12] G. Li, Y. Ding, and Y. Xie, “Tackling the qubit mapping problem for nisq-era quantum devices,” in _Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems_ , 2019, pp. 1001–1014.
  * [13] J. R. McClean, K. J. Sung, I. D. Kivlichan, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner, T. Hardikar, V. Havlíček, O. Higgott, C. Huang, J. Izaac, Z. Jiang, X. Liu, S. McArdle, M. Neeley, T. O’Brien, B. O’Gorman, I. Ozfidan, M. D. Radin, J. Romero, N. Rubin, N. P. D. Sawaya, K. Setia, S. Sim, D. S. Steiger, M. Steudtner, Q. Sun, W. Sun, D. Wang, F. Zhang, and R. Babbush, “OpenFermion: The electronic structure package for quantum computers,” 2019. [Online]. Available: http://arxiv.org/abs/1710.07629
  * [14] M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” 2002\.
  * [15] T. Patel, E. Younis, C. Iancu, W. de Jong, and D. Tiwari, _QUEST: Systematically Approximating Quantum Circuits for Higher Output Fidelity_. New York, NY, USA: Association for Computing Machinery, 2022, p. 514–528. [Online]. Available: https://doi.org/10.1145/3503222.3507739
  * [16] J. Preskill, “Quantum computing in the nisq era and beyond,” _Quantum_ , vol. 2, p. 79, Aug 2018. [Online]. Available: http://dx.doi.org/10.22331/q-2018-08-06-79
  * [17] M. Sarovar, T. Proctor, K. Rudinger, K. Young, E. Nielsen, and R. Blume-Kohout, “Detecting crosstalk errors in quantum information processors,” _Quantum_ , vol. 4, p. 321, sep 2020. [Online]. Available: https://doi.org/10.22331%2Fq-2020-09-11-321
  * [18] A. Shafaei, M. Saeedi, and M. Pedram, “Optimization of quantum circuits for interaction distance in linear nearest neighbor architectures,” in _Proceedings of the 50th Annual Design Automation Conference on - DAC ’13_. ACM Press, 2013, p. 1. [Online]. Available: http://dl.acm.org/citation.cfm?doid=2463209.2488785
  * [19] M. Y. Siraichi, V. F. d. Santos, S. Collange, and F. M. Q. Pereira, “Qubit allocation,” in _Proceedings of the 2018 International Symposium on Code Generation and Optimization_ , 2018, pp. 113–125.
  * [20] S. Sivarajah, S. Dilkes, A. Cowtan, W. Simmons, A. Edgington, and R. Duncan, “t— ket¿: A retargetable compiler for nisq devices,” _Quantum Science and Technology_ , 2020.
  * [21] E. Smith, M. G. Davis, J. M. Larson, E. Younis, C. Iancu, and W. Lavrijsen, “LEAP: Scaling numerical optimization based synthesis using an incremental approach,” 2021. [Online]. Available: http://arxiv.org/abs/2106.11246
  * [22] T. Tomesh, P. Gokhale, V. Omole, G. S. Ravi, K. N. Smith, J. Viszlai, X.-C. Wu, N. Hardavellas, M. R. Martonosi, and F. T. Chong, “Supermarq: A scalable quantum benchmark suite,” 2022. [Online]. Available: https://arxiv.org/abs/2202.11045
  * [23] X.-C. Wu, M. G. Davis, F. T. Chong, and C. Iancu, “Reoptimization of quantum circuits via hierarchical synthesis,” in _2021 International Conference on Rebooting Computing (ICRC)_ , 2021, pp. 35–46.
  * [24] E. Younis, C. C. Iancu, W. Lavrijsen, M. Davis, and E. Smith, “Berkeley quantum synthesis toolkit (bqskit) v1,” Tech. Rep., 2021. [Online]. Available: https://bqskit.lbl.gov
  * [25] E. Younis, K. Sen, K. Yelick, and C. Iancu, “QFAST: Quantum synthesis using a hierarchical continuous circuit space,” _2021 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , 2020. [Online]. Available: http://arxiv.org/abs/2003.04462


# Data-Efficient Control Barrier Function Refinement

Bolun Dai1, Heming Huang1, Prashanth Krishnamurthy1, Farshad Khorrami1
1Control/Robotics Research Laboratory, Electrical & Computer Engineering
Department, Tandon School of Engineering, New York University, Brooklyn, NY,
11201<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Control barrier functions (CBFs) have been widely used for synthesizing
controllers in safety-critical applications. When used as a safety filter, it
provides a simple and computationally efficient way to obtain safe controls
from a possibly unsafe performance controller. Despite its conceptual
simplicity, constructing a valid CBF is well known to be challenging,
especially for high-relative degree systems under nonconvex constraints.
Recently, work has been done to learn a valid CBF from data based on a
handcrafted CBF (HCBF). Even though the HCBF gives a good initialization
point, it still requires a large amount of data to train the CBF network. In
this work, we propose a new method to learn more efficiently from the
collected data through a novel prioritized data sampling strategy. A priority
score is computed from the loss value of each data point. Then, a probability
distribution based on the priority score of the data points is used to sample
data and update the learned CBF. Using our proposed approach, we can learn a
valid CBF that recovers a larger portion of the true safe set using a smaller
amount of data. The effectiveness of our method is demonstrated in simulation
on a unicycle and a two-link arm.

##  I Introduction

With the rapid integration of automated systems in our daily lives, the
ability to ensure operation safety for automatic systems has gained more
attention [1] [2]. In optimal control, the safety requirements can be
formulated as optimization constraints. Following this line of thought, one of
the most popular approaches is model predictive control (MPC). However, the
main obstacle to directly applying MPC is its solution time. MPC is more
commonly used in a multi-layered control architecture, where MPC generates a
reference trajectory and a reactive control tracks the reference trajectory at
a higher frequency. One issue with this approach is that the constraint is
only active when possible unsafe scenarios might occur within the preview
horizon. When the time horizon is small, this leads to abrupt control actions,
which is not ideal for most safety-critical applications. Increasing the time
horizon is also not ideal since the solution time would also increase.

Recently, CBF-based methods [3] have become increasingly popular in the field
of safety-critical control. CBFs act as a certificate of the safety of a state
with respect to a safe set. For zeroing CBFs, a state is safe when the CBF at
that state is positive and unsafe when it is negative. Additionally, CBFs can
be used to alter unsafe control actions into safe control actions in a
minimally invasive way. Since CBFs are commonly used within a quadratic
program (QP) [4], the solution time is much faster than MPC. Additionally, it
is natural for CBF-based controllers to consider constraint boundaries even
when the system is far from them, thus, generating smooth control actions.
Although CBFs have many advantages, constructing a valid CBF that recovers the
true safe set is still challenging [5][6]. However, finding an HCBF that
recovers only a portion of the true safe set is relatively easy. Taking
advantage of this insight, work has been done in learning-based CBF refinement
starting from an HCBF [7]. This approach makes constructing CBFs much easier,
even for nonconvex safe sets. However, this approach requires a large dataset
collected online, which can be both time-consuming and expensive.

Similar to learning CBFs, other learning-based methods also suffer from high
sample complexity [8], e.g., deep reinforcement learning (DRL). Experience
replay (ER) is a commonly used method in off-policy DRL algorithms [9][10] to
reduce the sample complexity by recycling data collected from previous
episodes. When using ER, the data is stored in a replay buffer and sampled
when updating the control policy. In the vanilla version of ER, the data is
sampled uniformly. Prioritized experience replay (PER) [11] is then proposed
as a low sample complexity version of ER. For PER, the stored data is sampled
based on a priority score. Data with higher priorities are more likely to be
sampled, and vice versa. The use of PER greatly reduces the sample complexity
of off-policy DRL algorithms [12] [13]. In this work, we propose integrating a
PER strategy into learning-based CBF refinement.

In this paper, we propose an algorithmic approach to enable lower sample
complexity for learning a CBF starting from an HCBF. The main contribution of
this paper is threefold: (1) we developed a method that integrates PER into
the learning-based CBF refinement pipeline; (2) we provide a theoretical
analysis of how PER affects the learning process; (3) we demonstrate the
effectiveness of our approach using simulation studies. The remainder of this
paper is structured as follows. In Section II, a background summary regarding
CBFs and PER is provided. In Section III, the data-efficient CBF refinement
problem is formulated. In Section IV, we describe the CBF learning algorithm
along with the prioritized data sampling pipeline. In Section V, The
effectiveness of our method is demonstrated in simulation on a unicycle and a
two-link arm.. Finally, in Section VI, the paper is concluded.

##  II Preliminaries

In this section, we present a brief introduction to CBF-based control and PER.

###  II-A Control Barrier Function

Consider a set $\mathcal{C}$ that is the 0-superlevel set of a continuously
differentiable function
$h:\mathcal{C}\subset\mathcal{D}\rightarrow\mathbb{R}$, yielding

$\displaystyle\mathcal{C}$
$\displaystyle=\\{x\in\mathcal{D}\subset\mathbb{R}^{n}\mid h(x)\geq 0\\}$ (1a)
$\displaystyle\partial\mathcal{C}$
$\displaystyle=\\{x\in\mathcal{D}\subset\mathbb{R}^{n}\mid h(x)=0\\}$ (1b)
$\displaystyle\mathrm{Int}(\mathcal{C})$
$\displaystyle=\\{x\in\mathcal{D}\subset\mathbb{R}^{n}\mid h(x)>0\\}$ (1c)

where $\partial\mathcal{C}$ represents the boundary of $\mathcal{C}$ and
$\mathrm{Int}(\mathcal{C})$ represents the interior of $\mathcal{C}$.
Additionally, we assume that $\mathrm{Int}(\mathcal{C})$ is not an empty set,
i.e., $\mathrm{Int}(\mathcal{C})\neq\emptyset$, and $\mathcal{C}$ does not
contain any isolated points. We refer to this set $\mathcal{C}$ as the safe
set. Then, consider a system in control affine form

$\dot{x}=f(x)+g(x)u$ (2)

where the state is $x\in\mathbb{R}^{n}$, the control is $u\in\mathbb{R}^{m}$,
the drift is $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, and the control
influence matrix is $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times m}$. Both
$f$ and $g$ are also locally Lipschitz continuous. The closed-loop dynamics of
the system in (2) is

$\dot{x}=f_{\mathrm{cl}}(x)=f(x)+g(x)\pi(x)$ (3)

where the feedback controller is

$u=\pi(x)$ (4)

and $\pi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is locally Lipschitz
continuous. For any initial state $x_{0}\in\mathbb{R}^{n}$, there exists an
interval of existence

$I(x_{0})=[t_{0},t_{\mathrm{max}})$ (5)

such that $x(t)$ is the unique solution to (3) on $I(x_{0})$. When
$t_{\mathrm{max}}=\infty$, the system in (3) is forward complete. Using the
concepts mentioned above, the notions of forward invariance, safety with
respect to a set, and CBFs are defined as follows.

######  Definition 1 (Forward Invariance & Safety).

The system defined in (3) is forward invariant with respect to a set
$\mathcal{C}\subset\mathbb{R}^{n}$ if for every $x_{0}\in\mathcal{C}$, we have
$x(t)\in\mathcal{C}$ for all $t\in I(x_{0})$. A system that is forward
invariant with respect to $\mathcal{C}$ is said to be safe with respect to
$\mathcal{C}$. A controller that makes a closed-loop system safe with respect
to $\mathcal{C}$ is said to be safe with respect to $\mathcal{C}$.

######  Definition 2 (Control Barrier Function [3]).

Let $\mathcal{C}$ be the 0-superlevel set of a continuously differentiable
function $h:\mathcal{D}\rightarrow\mathbb{R}$, with $\partial h/\partial x\neq
0$ for all $x\in\partial\mathcal{C}$. Then, $h$ is a control barrier function
(CBF) on $\mathcal{C}$ if there exists an extended class
$\mathcal{K}_{\infty}$ function $\alpha:\mathbb{R}\rightarrow\mathbb{R}$ such
that for all $x\in\mathcal{D}$ the system defined in (2) satisfies

$\sup_{u\in\mathcal{U}}\Big{[}\frac{\partial h(x)}{\partial
x}\Big{(}f(x)+g(x)u\Big{)}\Big{]}\geq-\alpha(h(x))$ (6)

with $\mathcal{U}$ being the admissible set of controls.

The CBF constraint in (6) is used in CBF-based quadratic programs (CBFQPs)
[14] that find the closest safe control action to a possibly unsafe control
action

$\displaystyle\min_{u\in\mathcal{U}}\ $ $\displaystyle\
\|u-\pi_{\mathrm{perf}}(x)\|^{2}$ (7) $\displaystyle\mathrm{subject\ to}\ $
$\displaystyle\ \Big{[}\frac{\partial h(x)}{\partial
x}\Big{(}f(x)+g(x)u\Big{)}\Big{]}\geq-\alpha(h(x))$

with the possibly unsafe control action being generated by a performance
controller $\pi_{\mathrm{perf}}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$.
Assuming no control constraints (a common assumption in the CBF literature)
and a valid CBF, the CBFQP will always be feasible, and the corresponding
controller is Lipschitz continuous [15].

###  II-B Prioritized Experience Replay

Unlike vanilla ER, where the data is sampled uniformly, for PER, the data are
sampled based on a priority score $p$. In off-policy DRL algorithms, the
priority score is commonly chosen as the temporal difference (TD) error

$\mathrm{TD\ Error}=Q-(r+Q^{\prime})$ (8)

with $Q\in\mathbb{R}$ being the action-value function at this state,
$r\in\mathbb{R}$ being the reward, and $Q^{\prime}\in\mathbb{R}$ the action-
value function at the next state. Then, the sampling probability for each data
point can be computed as

$P(i)=\frac{p_{i}^{\alpha_{p}}}{\sum_{i=0}^{N}p_{i}^{\alpha_{p}}}$ (9)

where $P(i)$ represents the probability of sampling the $i$-th data point,
$p_{i}$ is the priority score for the $i$-th data point, $N$ is the number of
data points within the replay buffer, and $\alpha_{p}$ regulates how PER
behaves, with $\alpha_{p}=0$ corresponding to uniform sampling and
$\alpha_{p}=1$ corresponding to full prioritized sampling. The implementation
details of PER will be presented in Section IV-B.

##  III Problem Formulation

In this section, we present the problem formulation for data-efficient CBF
refinement. Consider the constraints

$c_{i}(x)\geq 0$ (10)

with $c_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $i=1,\cdots,n_{c}$.
Define $\mathcal{S}_{i}$ as

$\mathcal{S}_{i}=\\{x\mid c_{i}(x)>0\\}$ (11)

and the intersections of all $\mathcal{S}_{i}$’s as $\mathcal{S}$

$\mathcal{S}=\bigcap_{i=1}^{n_{c}}{\mathcal{S}_{i}}.$ (12)

Then, the safe set $\mathcal{C}$ is defined as the largest forward invariant
set in $\mathcal{S}$, which leads to the relationship

$\mathcal{C}\subseteq\mathcal{S}.$ (13)

Given the true safe set of our problem $\mathcal{C}$, we assume access [7] to
an HCBF $\widehat{h}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ that is continuously
differentiable, locally Lipschitz continuous and satisfies

$\widehat{\mathcal{C}}=\\{x\mid\widehat{h}(x)\geq 0\\}\subseteq\mathcal{C}.$
(14)

Without loss of generality, with $h$ defined in (1), we can assume the
relationship

$h(x)=\widehat{h}(x)+\Delta{h}(x)$ (15)

with $\Delta h:\mathbb{R}^{n}\rightarrow\mathbb{R}$ being a continuously
differentiable and locally Lipschitz continuous function. Additionally, we
assume that $\widehat{h}$ and $\Delta{h}$ have the same relative degree as
$h$. This is a mild assumption given that the relative degree of systems can
be inferred using first principles. In this paper, we aim to find an
algorithmic approach that learns $\Delta{h}$ in a data-efficient manner.

##  IV Method

In this section, we present our proposed approach for the data-efficient
refinement of CBFs. First, we describe the process of learning-based CBF
refinement. Then, we show how to improve the CBF learning procedure using PER.
Finally, we discuss how to evaluate the effectiveness of our proposed
approach.

###  IV-A Learning-Based CBF Refinement

Since the only unknown in (15) is $\Delta{h}$, we can use a neural network
parameterized by $\theta$ to estimate $\Delta{h}(x)$, which is denoted as
$\Delta\widehat{h}(x\mid\theta)$. Then, the learned CBF is defined as

$\widetilde{h}(x\mid\theta)=\widehat{h}(x)+\Delta\widehat{h}(x\mid\theta).$
(16)

We use a deep differential network (DDN) [16] to represent
$\Delta\widehat{h}(x\mid\theta)$, which given input $x$, outputs both
$\Delta\widehat{h}(x\mid\theta)$ and the analytical
$\partial\Delta\widehat{h}(x\mid\theta)/\partial x$. If smooth activation
functions are used, DDNs will be continuously differentiable. Another
challenge when learning the CBF is having no ground truth CBF values. However,
it is easier to recognize whether a state is safe using the constraints in
(10). We then write loss functions to enforce the conditions in (1), yielding,

$\displaystyle\mathcal{L}_{+}(\theta)$
$\displaystyle=\frac{1}{B}\sum_{x_{i}\in\mathcal{X}_{+}}\max\Big{(}0,-\widetilde{h}(x_{i}\mid\theta)\Big{)}$
(17a) $\displaystyle\mathcal{L}_{-}(\theta)$
$\displaystyle=\frac{1}{B}\sum_{x_{i}\in\mathcal{X}_{-}}\max\Big{(}0,\widetilde{h}(x_{i}\mid\theta)\Big{)}$
(17b)

with $B$ representing the batch size, $\mathcal{L}_{+}$ representing the loss
for safe states, $\mathcal{L}_{-}$ representing the loss for unsafe states,
$\mathcal{X}_{+}$ being the dataset containing safe interactions, and
$\mathcal{X}_{-}$ being the dataset containing unsafe interactions. It is much
easier to determine whether a state is safe using (10). One major issue with
the loss functions in (17) is that a sample $x$ that satisfies

$\widetilde{h}(x\mid\theta)\equiv 0$ (18)

gives zero loss. To better guide the learning, a distance function $d(x)$ is
introduced [7] that increases while going inside the interior of the safe set
and decreases while going inside the interior of the unsafe set. Using the
distance function, we have the new loss function as

$\displaystyle\mathcal{L}_{+}(\theta)$
$\displaystyle=\frac{1}{B}\sum_{x_{i}\in\mathcal{X}_{+}}\max\Big{(}0,-\widetilde{h}(x_{i}\mid\theta)+d(x_{i})\Big{)}$
(19a) $\displaystyle\mathcal{L}_{-}(\theta)$
$\displaystyle=\frac{1}{B}\sum_{x_{i}\in\mathcal{X}_{-}}\max\Big{(}0,\widetilde{h}(x_{i}\mid\theta)-d(x_{i})\Big{)}.$
(19b)

Additionally, to ensure that control actions exists at each point within the
safe set that satisfies (6), we use the loss function

$\displaystyle\mathcal{L}_{\nabla h}(\theta)=$ $\displaystyle\
\frac{1}{B}\sum_{x_{i}\in\mathcal{X}_{\nabla}}\max\Big{(}0,\frac{\partial\widetilde{h}(x_{i}\mid\theta)}{\partial
x}\dot{x}_{i}$ $\displaystyle-\alpha(\widetilde{h}(x_{i}\mid\theta))\Big{)}$
(20)

with $\mathcal{X}_{\nabla}$ being the replay buffer for the CBF constraint
loss. Then, we have the final loss function as

$\mathcal{L}(\theta)=\mathcal{L}_{+}(\theta)+\lambda\mathcal{L}_{-}(\theta)+\mathcal{L}_{\nabla
h}(\theta)$ (21)

where $\lambda\in\mathbb{R}_{+}$ weights the importance of the unsafe loss.
Since it is much more dangerous to recognize an unsafe state as a safe one
than the other way around, usually $\lambda>1$.

The training process is separated into multiple epochs, where each epoch
consists of one episode and multiple neural network updates. During each
episode, a CBFQP controller using the learned CBF is used to collect state-
action tuples, which are then stored in the replay buffers. After each episode
is finished, data are sampled from the replay buffers in batches to compute
the loss function in (21) and its gradient. Then, using a stochastic gradient
descent (SGD) algorithm, e.g., ADAM [17], the weight vector of the neural
network is updated. For each update, new batches of data will be sampled. Note
that, during training, the system would encounter unsafe states. We can
utilize simulation environments (e.g., PyBullet or MuJoCo) or create safe
versions of the real scenario (e.g., replace obstacles with virtual
projections) to avoid endangering the real system.

###  IV-B Prioritized Sampling

In this section, we will discuss how to incorporate a PER buffer into the CBF
learning pipeline. To utilize PER, we first need to construct a priority
score. For the unsafe buffer, we choose

$p_{i}^{-}=\max\Big{(}0,\widetilde{h}(x_{i}\mid\theta)-d(x_{i})\Big{)}.$ (22)

For the safe buffer, we choose

$p_{i}^{+}=\max\Big{(}0,-\widetilde{h}(x_{i}\mid\theta)+d(x_{i})\Big{)}.$ (23)

For the CBF constraint buffer, we choose

$p_{i}^{\nabla}=\max\Big{(}0,\frac{\partial\widetilde{h}(x_{i}\mid\theta)}{\partial
x}\dot{x}_{i}-\alpha(\widetilde{h}(x_{i}\mid\theta))\Big{)}.$ (24)

To make the prioritized buffer computationally more efficient, we can use a
sum tree data structure [11] as shown in Fig. 1, which is popular among PER
[11] implementations. In the sum tree, each leaf node represents a data point
with a value corresponding to the data point priority. Each parent node has
two children. The parent node value is the sum of the children node’s values.
The value of the root node is the sum of all leaf node values, which
corresponds to $\sum_{i=1}^{N}{p_{i}}$. When adding a data point, it occupies
the next available leaf node, and the values of all of its parent nodes are
updated. This update has a time complexity of $O(\log{n})$. When sampling a
data point, a random number $\xi$ is sampled from $[0,\sum_{N}p_{i}]$. Then,
starting from the root node, $\xi$ is first compared with the value of the
left child node; if smaller, the left child node is considered; if larger,
$\xi\leftarrow\xi-p_{l}$, with $p_{l}$ being the value of the left children
node, and the right child node is considered. This operation is recursively
applied until one of the leaf nodes is reached. This process is equivalent to
drawing a line from zero to $\sum_{N}p_{i}$ and letting each leaf node occupy
a space with length $p_{i}$. Then the probability of $\xi$ landing inside the
portion belonging to node $i$ has the probability of $p_{i}/\sum_{N}p_{i}$,
which is illustrated in Fig. 1.

Figure 1: This figure illustrates the structure of a sum tree with four leaf
nodes and its corresponding probability distribution.

Following the formulation in [11], a fake priority is assigned when data is
first added to the buffer; the priority is then updated when the loss function
is computed. The fake priority is chosen to be
$p_{\mathrm{fake}}=p_{\mathrm{max}}+\epsilon$ with $p_{\mathrm{max}}$ being
the largest priority currently present in the PER buffer and
$\epsilon\in\mathbb{R}_{+}$ being a small value. This choice of
$p_{\mathrm{fake}}$ ensures new data has the largest priority. Additionally,
the relationship

$BN_{U}\geq T$ (25)

must hold, with $N_{U}$ being the number of updates per epoch and $T$ being
the number of time steps within each episode. These conditions ensure each new
data point will be sampled at least once and an actual priority value is
assigned. Next, we discuss how the prioritized sampling scheme can expedite
learning in the CBF refinement setting.

######  Remark 1.

When using (19) and (IV-A), the replay buffer can be separated into two parts,
one that has zero losses and the other that has nonzero losses. For the zero
loss portion, the sampling priority is the lowest, and the corresponding
gradient value is zero.

######  Remark 2.

For the nonzero loss portion, (19) and (IV-A) have the same form as

$\mathcal{L}_{L1}=\frac{1}{\bar{N}}\sum_{i=1}^{\bar{N}}|F_{\theta}(x_{i})-y_{i}|$
(26)

with $\bar{N}\in\mathbb{Z}_{+}$ being the size of the replay buffer,
$F_{\theta}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n_{o}}$ the neural network,
$y_{i}\in\mathbb{R}^{n_{o}}$ the training labels, and $n_{o}$ the dimension of
the output ($n_{o}=1$ in our case).

######  Theorem 1.

For a constant $c$ that depends on $\theta$, we have

$\displaystyle\
c\cdot\mathbb{E}_{(x,y)\sim\mathrm{Prioritized}}\Big{[}\nabla_{\theta}|F_{\theta}(x)-y|\Big{]}$
$\displaystyle=\mathbb{E}_{(x,y)\sim\mathrm{Uniform}}\Big{[}\frac{1}{2}\nabla_{\theta}(F_{\theta}(x)-y)^{2}\Big{]}$
(27)

###### Proof.

$\displaystyle\
\mathbb{E}_{(x,y)\sim\mathrm{Prioritized}}\Bigg{[}\frac{\partial|F_{\theta}(x)-y|}{\partial\theta}\Bigg{]}$
$\displaystyle=\ $ $\displaystyle\
\frac{1}{\sum_{i=1}^{\bar{N}}|F_{\theta}(x_{i})-y_{i}|}\sum_{i=1}^{\bar{N}}|F_{\theta}(x_{i})-y_{i}|\frac{\partial|F_{\theta}(x_{i})-y_{i}|}{\partial\theta}$
$\displaystyle=\ $ $\displaystyle\
\frac{1}{\sum_{i=1}^{\bar{N}}|F_{\theta}(x_{i})-y_{i}|}\sum_{i=1}^{\bar{N}}|F_{\theta}(x_{i})-y_{i}|$
$\displaystyle\
\cdot\frac{\partial\Big{(}(F_{\theta}(x_{i})-y_{i})^{2}\Big{)}^{\frac{1}{2}}}{\partial(F_{\theta}(x_{i})-y_{i})^{2}}\frac{\partial(F_{\theta}(x_{i})-y_{i})^{2}}{\partial\theta}$
$\displaystyle=\ $ $\displaystyle\
\frac{1}{2\sum_{i=1}^{\bar{N}}|F_{\theta}(x_{i})-y_{i}|}\sum_{i=1}^{\bar{N}}\frac{\partial(F_{\theta}(x_{i})-y_{i})^{2}}{\partial\theta}.$
(28)

Then, setting

$c=\frac{1}{\bar{N}}\sum_{i=1}^{\bar{N}}|F_{\theta}(x_{i})-y_{i}|$ (29)

completes the proof. ∎

Theorem 1 states that by using prioritized sampling, the loss function is
equivalent to an L2 loss for the non-zero loss portion, while the loss
function is an L1 loss when using uniform sampling, as in [18].

######  Theorem 2.

Define the absolute prediction error (APE) for the $i$-th data point in the
replay buffer under L1 loss with learning rate $\eta\in\mathbb{R}_{+}$ at time
$t$ as

$\delta_{t}(i)=|F_{\theta_{t}}(x_{i})-y_{i}|,$ (30)

and under L2 loss with learning rate $\eta$ at time $t$ as

$\widetilde{\delta}_{t}(i)=|\widetilde{F}_{\theta_{t}}(x_{i})-y_{i}|.$ (31)

Let $F_{\theta_{0}}(x)=\widetilde{F}_{\theta_{0}}(x)$. If
$\delta_{0}(i)=\widetilde{\delta}_{0}(i)\geq 1$, $\forall i=1,\cdots,\bar{N}$,
then there exists $\epsilon(i)\in(0,1]$, such that for an initial APE
$\delta_{0}$ that satisfies $\epsilon(i)\leq\delta_{0}$, the time it takes to
be reduced to $\epsilon(i)$ under L1 loss, defined as

$t_{\epsilon}=\min_{t}\\{t\geq 0:\delta_{t}(i)\leq\epsilon(i)\\},$ (32)

is longer than the time it takes under L2 loss, defined as

$\widetilde{t}_{\epsilon}=\min_{t}\\{t\geq
0:\widetilde{\delta}_{t}(i)\leq\epsilon(i)\\},$ (33)

i.e., $\widetilde{t}_{\epsilon}\leq t_{\epsilon}$.

###### Proof.

Given

$\frac{\partial\delta_{t}(i)}{\partial
F_{\theta_{t}}(x_{i})}=\frac{\partial\widetilde{\delta}_{t}(i)}{\partial
F_{\theta_{t}}(x_{i})}=\frac{\partial\mathcal{L}_{L1}(F_{\theta_{t}},y)}{\partial
F_{\theta_{t}}(x_{i})}=\mathrm{sgn}\\{F_{\theta_{t}}(x_{i})-y_{i}\\}$ (34)

and

$\frac{\partial\mathcal{L}_{L2}(F_{\theta_{t}},y)}{\partial
F_{\theta_{t}}(x_{i})}=F_{\theta_{t}}(x_{i})-y_{i},$ (35)

we have

$\displaystyle\frac{d\delta_{t}(i)}{dt}$
$\displaystyle={\frac{\partial\delta_{t}(i)}{\partial
F_{\theta_{t}}(x_{i})}\frac{dF_{\theta_{t}}(x_{i})}{dt}}$
$\displaystyle={\frac{\partial\delta_{t}(i)}{\partial
F_{\theta_{t}}(x_{i})}\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}(-\eta)\Big{(}\frac{\partial\mathcal{L}_{L1}(F_{\theta_{t}},y)}{\partial
F_{\theta_{t}}(x_{i})}\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}\Big{)}^{T}}$
$\displaystyle=-\eta\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}\Big{(}\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}\Big{)}^{T}$ (36)

and

$\displaystyle\frac{d\widetilde{\delta}_{t}(i)}{dt}$
$\displaystyle=\frac{\partial\widetilde{\delta}_{t}(i)}{\partial
F_{\theta_{t}}(x_{i})}\frac{dF_{\theta_{t}}(x_{i})}{dt}$
$\displaystyle={\frac{\partial\widetilde{\delta}_{t}(i)}{\partial
F_{\theta_{t}}(x_{i})}\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}(-\eta)\Big{(}\frac{\partial\mathcal{L}_{L2}(F_{\theta_{t}},y)}{\partial
F_{\theta_{t}}(x_{i})}\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}\Big{)}^{T}}$
$\displaystyle=-\eta\widetilde{\delta}_{t}(i)\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}\Big{(}\frac{\partial
F_{\theta_{t}}(x_{i})}{\partial\theta_{t}}\Big{)}^{T}$ (37)

with $\partial F_{\theta_{t}}(x_{i})/\partial\theta_{t}$ being a row vector.
This shows that when $\widetilde{\delta}_{t}(i)>1$, the APE will decrease
faster under L2 loss than under L1 loss. Then, based on the comparison lemma
[19], if $\delta_{0}(i)$ is greater than one, the APE under L2 loss will reach
APE equal to one faster than under L1 loss. The two APEs will intersect at
some value lower than one. Thus, for some value $\epsilon(i)\in(0,1]$, we have
the relationship $\widetilde{t}_{\epsilon}\leq t_{\epsilon}$. ∎

######  Remark 3.

For the CBF losses in (19) and (IV-A), when $\delta_{0}(i)<1$, L1 loss
converges faster to zero than L2 loss.

######  Remark 4.

For the CBF losses in (19) and (IV-A), the goal is not to converge to zero but
go past zero. For example, in the case of $\mathcal{L}_{-}$, the goal is not
to make $\widetilde{h}(x_{i}\mid\theta)-d(x_{i})=0$ but rather
$\widetilde{h}(x_{i}\mid\theta)-d(x_{i})\leq 0$.

Theorem 2 shows when APE is larger than one, instantaneously, L2 loss
converges faster, and when the APE is less than one, L1 loss converges faster.
In the learning-based CBF refinement setting, using the proposed prioritized
sampling and priority assignment scheme, the effect is a combination of both
L1 and L2 loss. To understand this effect, we show that the data sampling
process within one epoch can be seen as having two phases. For the first
phase, the majority of the data sampled are newly collected data points, and
according to our empirical studies, the majority have an APE of less than one.
In the second phase, after most of the new data points are sampled, if (25) is
satisfied, data points with the largest APEs will more likely get sampled. If
$B$ and $N_{U}$ are well chosen, the amount of data sampled in the second
phase with an APE lower than one will be small, and the majority of the data
sampled with have an APE larger than one. Therefore, this data collection and
sampling scheme can be approximately seen as applying L1 loss with uniform
sampling to low APE data points while applying L1 loss with prioritized
sampling (i.e., effectively L2 loss) to high APE data points. This makes our
proposed approach faster than using either L1 or L2 loss with uniform
sampling.

Figure 2: This figure compares the state trajectories generated by the learned
controller using uniform sampling and prioritized sampling for the unicycle
system. The darker blue square represents the square obstacle. The light blue
circle represents the unsafe region corresponding to the HCBF. Figure 3: This
figure compares the evaluation metric for using uniform and prioritized
sampling over the training process on a unicycle system. The vertical dashed
lines denote the epoch of convergence.

###  IV-C Evaluation Metric

To quantify the effectiveness of our proposed approach, we would need to find
a metric that represents the size of the learned safe set. We use an approach
inspired by Hamilton-Jacobi (HJ) reachability-based methods [20]. First, we
define a distance between the current state and the safe set boundary

$\ell(x)=\min\\{c_{1}(x),\cdots,c_{r}(x)\\}.$ (38)

Then, the evaluation metric is defined as

$m_{e}=\min_{t\in[0,T]}\Big{\\{}\ell(x(t))\mid
x(0)\in\mathcal{X}_{0},\dot{x}(t)=f_{cl}(x\mid\theta)\Big{\\}}$ (39)

with $\mathcal{X}_{0}$ being the set of initial states and
$f_{cl}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ the closed-loop system with
control actions generated by the learned CBF controller.

##  V Simulation

In this section, we show the effectiveness of our approach using simulation
studies.

###  V-A Unicycle

The task is to reach a target state while avoiding collision with a square
obstacle. The unicycle system has the dynamics

$\begin{bmatrix}\dot{x}\\\ \dot{y}\\\
\dot{\phi}\end{bmatrix}=\begin{bmatrix}\cos\phi&0\\\ \sin\phi&0\\\
0&1\end{bmatrix}\begin{bmatrix}v\\\ \omega\end{bmatrix}$ (40)

with $x$ denoting the position of the unicycle along the $x$ axis, $y$
denoting the position along the $y$ axis, and $\phi$ denoting the heading of
the unicycle. The control inputs are the linear velocity $v$ and the angular
velocity $\omega$. Given that it is difficult to construct an HCBF for a
square obstacle, we write an HCBF for a circular obstacle [21] that contains
the square obstacle

$\widehat{h}(\mathbf{x})=x^{2}+y^{2}+2xl\cos\phi+2yl\sin\phi+l^{2}-r^{2}$ (41)

where $r$ is the radius of the constructed circular obstacle and $l$ is a
predefined lookahead distance. This HCBF corresponds to the safe set under the
constraint

$x^{2}+y^{2}\geq r^{2}.$ (42)

During training, we use a proportional controller as the performance
controller

$\pi(x,y,\phi)=\begin{bmatrix}K_{v}e\\\ K_{\omega}(\beta-\phi)\end{bmatrix}$
(43)

with $K_{v}$ and $K_{\omega}$ being the control gains, and

$\displaystyle e$
$\displaystyle=\sqrt{(x-x_{\mathrm{des}})^{2}+(y-y_{\mathrm{des}})^{2}}$ (44a)
$\displaystyle\beta$
$\displaystyle=\mathrm{atan}2(y_{\mathrm{des}}-y,x_{\mathrm{des}}-x).$ (44b)

In our experiments, we set $K_{v}=0.75$ and $K_{\omega}=3.0$. The distance
function is chosen to be

$d(x)=\max(|x|,|y|)-s/2$ (45)

with $s$ being the side length of the square obstacle. For this example,
$\ell(x)$ is the same as the distance function. We trained two learned CBF
controllers, one using vanilla ER and the other one using PER. To ensure a
fair comparison, we fixed all random seeds. During training, we set
$\lambda=100$, and use a learning rate of $10^{-4}$. We trained the learned
controllers for 200 epochs, where at each epoch, we collected one episode of
data (each episode has 1000 time steps) and updated the learned controllers 10
times using a batch size of 128. The state trajectory generated by the two
learned controllers at different epochs is shown in Figure 2, and the
evaluation metric computed at each epoch is shown in Fig. 3. It can be easily
seen from both figures that the controller using prioritized sampling is able
to enlarge the estimated safe set faster. Additionally, we can see from Fig.
3, that after convergence, the evaluation metric for both of the controllers
is the same. Thus, we can conclude that using prioritized sampling expedited
the training while not making the performance of the learned controller
deteriorate.

###  V-B Two-Link Arm

Figure 4: This figure compares the end-effector trajectories generated by the
learned controller using uniform sampling and prioritized sampling. Figure 5:
Comparison of the evaluation metric between using uniform and prioritized
sampling over the training process on a two-link arm system. The vertical
dashed lines denote the epoch of convergence. Figure 6: Screenshot of the two-
link arm simulation environment. The task is to move from the initial position
to the target position while not colliding with the wall. The first link is
light purple and the second link is dark green. Each tile on the ground is
$1m\times 1m$.

The two-link arm system has the dynamics

$\begin{bmatrix}\dot{q}\\\ \ddot{q}\end{bmatrix}=\begin{bmatrix}\dot{q}\\\
-M^{-1}(q)[C(q,\dot{q})+G(q)]\end{bmatrix}+\begin{bmatrix}0\\\
M^{-1}(q)\end{bmatrix}\tau$ (46)

with $q\in\mathbb{R}^{2}$ representing the joint angles,
$\tau\in\mathbb{R}^{2}$ representing the joint torques,
$M(q)\in\mathbb{R}^{2\times 2}$ representing the inertia matrix (which is
positive definite), $C(q,\dot{q})\in\mathbb{R}^{2\times 1}$ representing the
Coriolis and centrifugal terms, and $G(q)\in\mathbb{R}^{2\times 1}$
representing the gravitational terms. The task is to move the end-effector
from the initial position to the target position while avoiding collision with
a wall, which is illustrated in Fig. 6. The link length is $1$m for both of
the links. The wall is positioned at $y=1.75$m. Both the two-link arm and the
wall can be seen as convex polytopes. Thus, the task is a polytopic obstacle
avoidance task. To deal with this task, work has been done using discrete-
time-CBF-based MPC [22] and finding a conservative kinematic controller [23].
In this example, we learn a CBF-based dynamic controller to solve this task.

Given the difficulty in constructing an HCBF for polytopic constraints, we
construct an HCBF that constraints the position of the end-effector to be less
than $\beta$ along $\bar{n}\in\mathbb{R}^{2}$

$\widehat{h}(x)=-\bar{n}^{T}J(q)\dot{q}+\gamma(\beta-\bar{n}^{T}\mathrm{FK}(q))$
(47)

with $J(q)\in\mathbb{R}_{3\times 2}$ being the Jacobian matrix for the end-
effector, $\mathrm{FK}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ the forward
kinematics function for the end-effector, $\bar{n}\in\mathbb{R}^{2}$ the
normal vector of the edge of the polytopic obstacle,
$d_{\mathrm{max}}\in\mathbb{R}$ the distance between the base of the two-link
arm and the obstacle edge along $\bar{n}$, and $\gamma\in\mathbb{R}_{+}$. In
this example $n=2$. The HCBF in (47) corresponds to the constraint

$p_{ee}^{T}\bar{n}\leq\beta$ (48)

with $p_{ee}\in\mathbb{R}^{2}$ representing the position of the end-effector.
For this example, we choose $\beta=1.5$ and $\bar{n}=[0,1]^{T}$. During
training, we set

$d(x)=\mathrm{Distance}(x)-\mathrm{Penetration}(x)$ (49)

where $\mathrm{Distance}(x)$ is computed using the Gilbert Johnson Keerthi
(GJK) algorithm [24] and $\mathrm{Penetration}(x)$ is computed using the
expanding polytope algorithm (EPA) [25]. We choose $\ell(x)$ to be the same as
the distance function. To ensure safety during training, we artificially
enlarge the unsafe set and treat all interactions within the new unsafe set as
unsafe. During training, we set $\lambda=100$, and we use a learning rate of
$10^{-4}$. We trained the learned controllers for 80 epochs, where at each
epoch, we collected one episode of data (each episode has 3000 time steps) and
updated the learned controllers 30 times using a batch size of 128. The value
of $d(x)$ along the generated trajectory at different epochs is shown in Fig.
4, and the evaluation metric over the learning process is shown in Fig. 5. We
can see that the learned controller using prioritized sampling converges
faster while not affecting performance.

##  VI Conclusion

In this paper, we proposed an algorithmic approach to learning a CBF in a
data-efficient fashion, starting from an HCBF. By utilizing prioritized
sampling, our proposed method is able to find an estimation of the true safe
set while using a smaller amount of data compared to uniform sampling. We
demonstrated the effectiveness of our approach on a unicycle and a two-link
arm system. In future works, we plan to study the integration of a learned
performance controller into the proposed framework as well as perform
experimental studies on robotic systems testbeds.

## References

  * [1] B. Dai, P. Krishnamurthy, A. Papanicolaou, and F. Khorrami, “State constrained stochastic optimal control using LSTMs,” in _Proceedings of American Control Conference, New Orleans, LA_ , May 2021, pp. 1294–1299.
  * [2] B. Dai, V. R. Surabhi, P. Krishnamurthy, and F. Khorrami, “Learning locomotion controllers for walking using deep FBSDE,” _CoRR_ , vol. abs/2107.07931, 2021.
  * [3] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in _Proceedings of European Control Conference, Naples, Italy_ , June 2019, pp. 3420–3431.
  * [4] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs with application to adaptive cruise control,” in _Proceedings of IEEE Conference on Decision and Control, Los Angeles, CA_ , December 2014, pp. 6271–6278.
  * [5] J. J. Choi, D. Lee, K. Sreenath, C. J. Tomlin, and S. L. Herbert, “Robust control barrier-value functions for safety-critical control,” in _Proceedings of IEEE Conference on Decision and Control, Austin, TX_ , December 2021, pp. 6814–6821.
  * [6] S. Chen, Z. Wu, and P. D. Christofides, “Statistical machine-learning-based predictive control using barrier functions for process operational safety,” _Computers and Chemical Engineering_ , vol. 163, p. 107860, 2022.
  * [7] B. Dai, P. Krishnamurthy, and F. Khorrami, “Learning a better control barrier function,” in _Proceedings of IEEE Conference on Decision and Control, Cancún, Mexico_ , December 2022, pp. 945–950.
  * [8] S. Wei, P. Krishnamurthy, and F. Khorrami, “Neural Lyapunov control for nonlinear systems with unstructured uncertainties,” in _Proceedings of American Control Conference, San Diego, CA_ , May 2023, (To Appear).
  * [9] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, pp. 529–533, 2015.
  * [10] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” in _Proceedings of International Conference on Learning Representations, San Juan, PR_ , May 2016.
  * [11] T. Schaul, J. Quan, I. Antonoglou, and D. Silver, “Prioritized experience replay,” in _Proceedings of International Conference on Learning Representations, San Juan, PR_ , May 2016.
  * [12] Y. Pan, J. Mei, A.-M. Farahmand, M. White, H. Yao, M. Rohani, and J. Luo, “Understanding and mitigating the limitations of prioritized experience replay,” in _Proceedings of Uncertainty in Artificial Intelligence, Eindhoven, Netherlands_ , August 2022, pp. 1561–1571.
  * [13] S. Fujimoto, D. Meger, and D. Precup, “An equivalence between loss functions and non-uniform sampling in experience replay,” in _Proceedings of Advances in Neural Information Processing Systems, Virtual_ , vol. 33, December 2020, pp. 14 219–14 230.
  * [14] T. Gurriet, A. Singletary, J. Reher, L. Ciarletta, E. Feron, and A. D. Ames, “Towards a framework for realizable safety critical control through active set invariance,” in _Proceedings of the ACM/IEEE International Conference on Cyber-Physical Systems, Porto, Portugal_ , April 2018, pp. 98–106.
  * [15] B. Morris, M. J. Powell, and A. D. Ames, “Sufficient conditions for the Lipschitz continuity of qp-based multi-objective control of humanoid robots,” in _Proceedings of the IEEE Conference on Decision and Control, Florence, Italy_ , December 2013, pp. 2920–2926.
  * [16] M. Lutter, C. Ritter, and J. Peters, “Deep Lagrangian networks: Using physics as model prior for deep learning,” in _Proceedings of International Conference on Learning Representations, New Orleans, LA_ , May 2019\.
  * [17] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proceedings of International Conference on Learning Representations, San Diego, CA_ , May 2015.
  * [18] A. Robey, H. Hu, L. Lindemann, H. Zhang, D. V. Dimarogonas, S. Tu, and N. Matni, “Learning control barrier functions from expert demonstrations,” in _IEEE Conference on Decision and Control, Jeju Island, South Korea_ , December 2020, pp. 3717–3724.
  * [19] H. K. Khalil, _Nonlinear Control_. Pearson New York, 2015.
  * [20] S. Bansal, M. Chen, S. L. Herbert, and C. J. Tomlin, “Hamilton-Jacobi reachability: A brief overview and recent advances,” in _Proceedings of IEEE Conference on Decision and Control, Melbourne, Australia_ , December 2017, pp. 2242–2253.
  * [21] Y. Emam, P. Glotfelter, Z. Kira, and M. Egerstedt, “Safe model-based reinforcement learning using robust control barrier functions,” _CoRR_ , vol. abs/2110.05415, 2021.
  * [22] A. Thirugnanam, J. Zeng, and K. Sreenath, “Duality-based convex optimization for real-time obstacle avoidance between polytopes with control barrier functions,” in _Proceedings of American Control Conference, Atlanta, GA_ , June 2022, pp. 2239–2246.
  * [23] A. Singletary, W. Guffey, T. G. Molnar, R. Sinnet, and A. D. Ames, “Safety-critical manipulation for collision-free food preparation,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 10 954–10 961, 2022.
  * [24] E. G. Gilbert, D. W. Johnson, and S. S. Keerthi, “A fast procedure for computing the distance between complex objects in three-dimensional space,” _IEEE Journal on Robotics and Automation_ , vol. 4, no. 2, pp. 193–203, 1988.
  * [25] G. Van Den Bergen, “Proximity queries and penetration depth computation on 3d game objects,” in _Proceedings of Game Developers Conference, San Jose, CA_ , vol. 170, 2001.