# Fast quantum circuit cutting with randomized measurements

Angus Lowe  Xanadu, Toronto, ON, M5G 2C8, Canada  Matija Medvidović  Xanadu,
Toronto, ON, M5G 2C8, Canada  Center for Computational Quantum Physics,
Flatiron Institute, 162 5th Avenue, New York 10010, USA  Department of
Physics, Columbia University, New York 10027, USA  Anthony Hayes  Xanadu,
Toronto, ON, M5G 2C8, Canada  Lee J. O’Riordan  Xanadu, Toronto, ON, M5G 2C8,
Canada  Thomas R. Bromley  Xanadu, Toronto, ON, M5G 2C8, Canada  Juan Miguel
Arrazola  Xanadu, Toronto, ON, M5G 2C8, Canada  Nathan Killoran  Xanadu,
Toronto, ON, M5G 2C8, Canada

###### Abstract

We propose a new method to extend the size of a quantum computation beyond the
number of physical qubits available on a single device. This is accomplished
by randomly inserting measure-and-prepare channels to express the output state
of a large circuit as a separable state across distinct devices. Our method
employs randomized measurements, resulting in a sample overhead that is
$\widetilde{O}(4^{k}/\varepsilon^{2})$, where $\varepsilon$ is the accuracy of
the computation and $k$ the number of parallel wires that are “cut” to obtain
smaller sub-circuits. We also show an information-theoretic lower bound of
$\Omega(2^{k}/\varepsilon^{2})$ for any comparable procedure. We use our
techniques to show that circuits in the Quantum Approximate Optimization
Algorithm (QAOA) with $p$ entangling layers can be simulated by circuits on a
fraction of the original number of qubits with an overhead that is roughly
$2^{O(p\kappa)}$, where $\kappa$ is the size of a known balanced vertex
separator of the graph which encodes the optimization problem. We obtain
numerical evidence of practical speedups using our method applied to the QAOA,
compared to prior work. Finally, we investigate the practical feasibility of
applying the circuit cutting procedure to large-scale QAOA problems on
clustered graphs by using a $30$-qubit simulator to evaluate the variational
energy of a $129$-qubit problem as well as carry out a $62$-qubit
optimization.

##  I Introduction

In this work, we consider combining measurement outcomes from multiple quantum
circuits to estimate the expectation value at the output of a larger quantum
circuit. Computing such expectation values is a key step in many proposals for
quantum algorithms such as the QAOA [1, 2, 3], VQE [4], or in the estimation
of output probabilities, for example, for quantum classifiers [5, 6].

Given a circuit comprising $m$ gates, the standard procedure to accomplish
this task is to repeatedly run the desired circuit and measure in some basis,
for a total quantum runtime on the order of $m/\varepsilon^{2}$.
Alternatively, one may perform a full classical simulation of the circuit,
which is believed to require a runtime exponential in the number of qubits $n$
without strong assumptions about the structure of the circuit, for example, a
bounded treewidth [7], low stabilizer-rank [8, 9, 10], or low negativity [11].
Based on current hardware limitations both approaches become infeasible for
moderately large $n$: the time cost of the classical simulation is
prohibitive, while on the other hand current quantum devices are limited to
small numbers of qubits. This motivates considering a hybrid approach, where
results obtained on smaller quantum devices are used to solve the estimation
task.

For instance, Ref. [10] suggests constructing “virtual qubits” which enable
the simulation of an otherwise intractable computation. Similarly, Ref. [12]
presented a framework for simulating clustered quantum circuits. A third
framework, presented in Refs. [13, 14], focuses on removing entangling gates
between sub-circuits. For these methods, as well as in follow-up work [15, 16,
17, 18, 19], one takes advantage of large, disconnected components of the
quantum circuit which may be obtained by removing a subset of wires or gates.
We refer to such methods collectively as quantum circuit cutting. In related
work, Ref. [20] gives an algorithm for speeding up classical algorithms for
3SAT with small quantum devices.

Our main result is a fast quantum circuit cutting method based on separating
weakly entangled circuits through randomized measurements. Compared to the
approach in [12] in which single-qubit Pauli measurements are performed, our
method offers a quadratic improvement in the overhead in certain cases of
interest. As a consequence, the number of wires that can be cut under a fixed
budget on the runtime is approximately doubled using this approach.
Furthermore, our method likely outperforms all other proposed circuit cutting
methods for simulating quantum computation using circuits with appropriate
structure (see Table 1).

Our approach combines intuition which stems from quantum tomography as well as
the notion of a matrix-product state (MPS). Specifically, the algorithm can be
viewed as an attempt to classically simulate the entanglement between qudits
in an MPS with low bond dimension using repeated mid-circuit measurements.
This measurement procedure is described by measure-and-prepare channels
(equivalent to entanglement-breaking channels [21]), enabling the expression
of the post-measurement state as a separable state across distinct devices.
The measurements we choose for this task are based on unitary 2-designs, which
are related to the sample-optimal measurements in various quantum learning
settings [22, 23, 24, 25].

A point of departure from previous work is that our procedure uses
synchronized measurements and preparations. For circuits with an MPS-like
structure — as considered in [26] for example — this requirement may be
satisfied by repeatedly executing circuits on a single device. However, if the
so-called communication graph [12] of the sub-circuits contains a cycle, then
our method genuinely requires multiple circuits, separated in space rather
than time.

To understand the source of our speedup, we turn to the form of the equation
we employ in our proposal. In any dimension $d>0$, one may consider a
collection of linear maps
$\Phi_{i}:\mathsf{L}(\mathbb{C}^{d})\to\mathsf{L}(\mathbb{C}^{d})$ such that

$\mathrm{id}=\sum_{i}a_{i}\Phi_{i}\;,$ (1)

where $\mathrm{id}:\mathsf{L}(\mathbb{C}^{d})\to\mathsf{L}(\mathbb{C}^{d})$ is
the identity map. Throughout, we let $\mathsf{L}(\mathbb{C}^{d})$ denote the
set of linear operators acting on $\mathbb{C}^{d}$. See Appendix A for a
comprehensive set of definitions. As suggested in Ref. [12], we choose an
identity of this form such that the substitution of any individual map
$\Phi_{i}$ in a suitable location within the circuit enables the execution of
a quantum computation which would otherwise be intractable. We refer to this
as wire cutting. A similar identity is employed to decompose 2-qubit
operations into product operators in Refs. [10, 14, 13, 15], and we refer to
this approach as gate cutting. In either case, each application of the
identity results in a multiplicative overhead in the runtime scaling with the
1-norm of the coefficients, e.g., $\sum_{i}|a_{i}|$ in the case where the
identity takes the form in Eq. (1). This is similar to the effect of the
negativity quantity associated with fully classical simulation method of [11]
based on quasiprobability sampling.

As previously observed in [13, 17, 15], the overhead in circuit cutting
depends on the choice of linear maps used in the decomposition. In particular,
Ref. [15] shows a lower bound on the overhead for gate cutting in several
cases, while Ref. [17] proposed a simulated-annealing-based search for
solutions to Eq. (1) to attempt to minimize the above 1-norm quantity.
Inspired by recent work demonstrating the utility of randomized measurements
[25], we give an explicit decomposition of the form in Eq. (1) whose 1-norm
scales with the dimension of the subspace upon which the channel acts. We
argue that this bound results in our method outperforming the state-of-the-art
for a natural problem (see Table 1), and are able to prove that it is at most
(roughly) quadratically worse than the best possible wire cutting method by
deriving an information-theoretic lower bound.

Figure 1:  Example of our quantum circuit cutting scheme. One of two measure-
and-prepare channels ($z=0,1$) is applied at random. If $z=0$, a randomized
measurement based on a unitary 2-design (e.g., a random Clifford) $U$ is
performed. If $z=1$, the qubits are traced out, and a random basis state
$\left|y^{\prime}\right\rangle$ is prepared. The output is an unbiased
estimator of the expectation value at the output of the original circuit.

We then apply our observations to derive sufficient conditions for the
efficient hybrid simulation of circuits which arise in the Quantum Approximate
Optimization Algorithm (QAOA) [1, 27, 28]. This continues a line of research
where the goal is to identify restricted families of circuits which can be
simulated using smaller devices than might naively be expected. For instance,
Ref. [10] proposed a notion of sparse circuits for which a small number of
qubits could be efficiently traded for classical computation, while Refs. [12,
29] show that weakly interacting, clustered Hamiltonians can be simulated
efficiently for short times through wire cutting. We show that $p$-layer QAOA
circuits can be efficiently simulated in a hybrid manner up to depth $p=O(\log
n)$, provided the underlying graph encoding the optimization problem has a
small balanced vertex separator. This addresses an open problem from Ref. [12]
which asked for good ways to partition circuits based on the structure of the
problem at hand. Our observations should also carry over to Trotter-Suzuki-
based circuits for the Hamiltonian simulation of 2-local Hamiltonians whose
interaction graph has small vertex separators, as opposed to the edge
separators considered in [12]. We then show how to sample from the output
distribution of a QAOA circuit to recover bitstrings encoding optimal
solutions by combining measurement outcomes obtained from the smaller
circuits.

We remark that prior work [30] proposed utilizing balanced vertex separators
to split instances of the Max-Cut problems into smaller sub-problems for the
QAOA, without considering circuit cutting methods. In addition, Refs. [31, 32]
consider applying the circuit cutting method proposed in [12] to solve a
modified version of the QAOA. However, these algorithms require estimating the
measurement distribution of the sub-circuits as a subroutine, which limits the
practicality of such proposals.

We demonstrate practical speedups by classically simulating our quantum
circuit cutting method applied to small instances of the QAOA and comparing to
the procedure in [12]. These numerical experiments are enabled by leveraging
open-source circuit cutting functionality available in PennyLane [33]. We then
show that a circuit cutting procedure can be performed on large-scale QAOA
problems. Using the tensor-contraction based method introduced in [12] and
available in PennyLane, we break the full QAOA circuit into multiple circuit
fragments of at most $30$ qubits and execute them on a cluster of NVIDIA A100
40GB GPUs, acting as simulator substitutes for practical quantum hardware
devices. As a result, the variational energy of a 2-layer QAOA circuit of 129
qubits is evaluated, as well as a full 1- and 2-layer QAOA optimization of a
62-qubit circuit.

##  II Fast circuit cutting with randomized measurements

We first provide details for the computational model considered in this work.
A quantum circuit is a directed acyclic network of gates (vertices) connected
by wires (edges) which represent the qubits upon which the gates act. In
general, gates represent quantum channels, i.e., completely positive and
trace-preserving linear maps. The overall action of the circuit is some
quantum channel acting on the input qubits, which throughout this paper we
take to be $\rho_{0}^{\otimes n}:=(|0\rangle\langle 0|)^{\otimes n}$ for a
circuit with $n$ input qubits. The size of the circuit is the total number of
gates. As in Peng et al. [12], we adopt a model of quantum computation in
which the goal is to estimate the expectation value of some diagonal
observable $O_{f}=\sum_{x\in\\{0,1\\}^{n}}f(x)|x\rangle\langle x|$, where
$f:\\{0,1\\}^{n}\to[-1,1]$ is a classically efficiently computable function.

More specifically, given an input circuit on $n$ qubits representing the
action of some channel $\mathcal{N}$, the task is to estimate the value
$\textnormal{Tr}(O_{f}\mathcal{N}(\rho_{0}^{\otimes n}))$ to within additive
error $\varepsilon$. This is a more general version of the computational model
adopted in [10], and — up to single-qubit rotations prior to measurement —
encompasses the “quantum mean-value problem” defined in [34]. The model is
motivated by the fact that measurement in an alternative, efficiently
implementable basis, can be absorbed into the unitary operation of the
circuit, and an observable whose spectral norm is larger by a polynomial
factor only incurs a polynomial overhead in the runtime for the estimation
task. In the remainder of this paper we further assume that $f$ can be
evaluated in constant time, for succinctness in the statement of our results.

We now define the two different measure-and-prepare channels which are used in
our scheme. Let $\\{\left|j\right\rangle\\}_{j=1}^{d}$ denote the standard
basis for a $d$-dimensional complex Euclidean space. In any dimension $d>0$ we
may consider a random POVM $\\{U|j\rangle\langle j|U^{\dagger}\\}_{j=1}^{d}$
comprising rank-1 measurement operators where $U\in\mathsf{U}(\mathbb{C}^{d})$
is a random unitary operator (matrix-valued random variable) which forms a
unitary 2-design (see Appendix A). For a fixed choice of such a random POVM,
the first channel we consider
$\Psi_{0}:\mathsf{L}(\mathbb{C}^{d})\to\mathsf{L}(\mathbb{C}^{d})$ is then
defined as

$\displaystyle\Psi_{0}(X)=\operatorname{{\mathbb{E}}}_{U}\left[\sum_{j=1}^{d}\left\langle
j\right|U^{\dagger}XU\left|j\right\rangle U|j\rangle\langle
j|U^{\dagger}\right],$ (2)

for every $X\in\mathsf{L}(\mathbb{C}^{d})$. This mapping is a measure-and-
prepare channel, i.e., it describes a measurement performed on the
$d$-dimensional register (using the random POVM described above) followed by
preparing a new register in the state corresponding to the measurement outcome
that was observed. It would also suffice for our purposes to consider a fixed
rank-1 POVM whose measurement operators are based on state 2-designs, since
the resulting measure-and-prepare channel would be identical. However, we
focus on randomized basis measurements to make explicit how such a channel
could be implemented in practice.

| Type | Overhead | Communication
---|---|---|---
Thm. 1 in [10] | Gate cutting | $2^{O(nD)}$ | None
Thm. 1 in [12] | Wire cutting | $16^{k}$ | None
Thm. 5 in [14] | Gate cutting | $\min\\{9^{kD},16^{k}\\}$ | None
Eq. (5.1) in [15] | Gate cutting | $4^{kD}$ | LOCC (shared entanglement)
Thm. II.2 | Wire cutting | $4^{k}$ | Classical communication
Table 1: Upper bounds on the runtime of circuit cutting methods applied to
circuits of the form given in Fig. 1, omitting polynomial factors in $n$, $k$,
$D$, and $\varepsilon$. It is assumed that the task is to simulate a quantum
computation with this circuit using a device comprising at most $n/2+k$ qubits
(“doubling” the number of qubits). The “Type” column identifies whether the
procedure exploits decompositions of gates or of the identity operation
applied to wires. Here $D$ is a sparseness parameter, denoting the assumption
that the qubits being cut participate in at most $D$ multi-qubit gates. Note
that in this case, “Classical communication” may be achieved by recycling
qubits on a single quantum device.

The second measure-and-prepare channel we consider is the completely
depolarizing channel
$\Psi_{1}:\mathsf{L}(\mathbb{C}^{d})\to\mathsf{L}(\mathbb{C}^{d})$ defined as

$\displaystyle\Psi_{1}(X)=\textnormal{Tr}(X)\frac{\mathds{1}}{d},$ (3)

for every $X\in\mathsf{L}(\mathbb{C}^{d})$. The following lemma implies that
measure-and-prepare channels of the above form can be used to develop an
improved approach to circuit cutting.

######  Lemma II.1.

Let $d$ be a positive integer and $\Psi_{0}$, $\Psi_{1}$ be the channels
defined in Eqs. (2) and (3), respectively, acting on $d$-dimensional states.
Define the Bernoulli random variable $z\in\\{0,1\\}$ to be equal to $1$ with
probability $d/(2d+1)$. It holds that

$\displaystyle\mathrm{id}=(2d+1)\operatorname{{\mathbb{E}}}_{z}\left[\
(-1)^{z}\ \Psi_{z}\right].$ (4)

Eq. (4) suggests a sampling-based approach to wire cutting by randomly
inserting channels on the wires to be “cut”. This idea is depicted in Fig. 1,
and is realized by the procedure described in Algorithm 1 in Appendix B.2. The
following theorem then bounds the overhead incurred from using this method to
perform wire-cutting on a natural example of a clustered circuit. Once again,
we ignore the classical time to compute the post-processing function $f$ in
the statement of this theorem.

######  Theorem II.2 (Bipartitioning circuits).

Let $C$ be a size-$m$ quantum circuit acting on $n$ qubits which is a
composition of circuits $C_{A},C_{B}$ acting non-trivially on sets of qubits
$A,B\subseteq[n]$, respectively. If $|A\cap B|\leq k$, then quantum
computation using $C$ can be simulated to within accuracy $\varepsilon$ in
time $O(4^{k}(m+k^{2})/\varepsilon^{2})$ by a quantum circuit acting on at
most $\max\\{|A|,|B|\\}$ qubits using the procedure described in Algorithm 1.

The proof of the theorem is given in Appendix B.3. The theorem follows from
Lemma II.1 by using the fact that random $k$-qubit Clifford operators: (i)
form a unitary $3$-design [35, 36, 37], (ii) can be sampled efficiently [38],
and (iii) can be implemented using poly-depth circuits [39]. We remark that
while random Clifford operations on many qubits can be challenging to
implement in practice, we expect they may be feasible for small values of $k$,
which is the regime suitable for efficient circuit cutting.

When $|A|\approx|B|\approx n/2$ and $k$ is a small constant, this technique
effectively doubles the number of qubits which can be simulated given a
quantum device. In Table 1 we compare the overheads of several circuit cutting
methods assuming tightness of the upper bounds shown in each respective work
(see Sec. III.3 for numerical evidence that this is approximately the case for
the wire cutting method in Ref. [12]). In the table, the setting being
considered is a special case of that in Theorem II.2: we have a composition of
two circuits acting on wires belonging to overlapping sets $A$ and $B$ with
$|A|=|B|=n/2+k$, and the goal is to run the quantum computation using a device
limited to $n/2+k$ wires. The scenario is depicted in the left-hand side of
Fig. 1. To enable the comparison with gate-cutting methods, it is necessary
that we assume each circuit satisfies a certain sparseness property,
introduced in Ref. [8]. Namely, we assume that the qubits in $A\cap B$
participate in at most $D$ multi-qubit gates. Otherwise, gate-cutting methods
may incur an unbounded overhead to accomplish the simulation task.

A natural question is whether the method described in Algorithm 1 saturates
the best possible overhead for wire cutting. For example, one might wonder
whether a subexponential runtime, or runtime on the order of $2^{ck}$ for some
$0<c<1$ can be achieved. We rule out these possibilities using an information-
theoretic lower bound on the sample complexity of a task that reduces to wire-
cutting, which we describe in more detail in Appendix C. Roughly speaking, a
procedure for wire cutting enables one to succeed at a difficult quantum state
discrimination task by cutting $n$ wires, similar to the lower bounds on
classical shadows which appear in [23].

######  Theorem II.3 (Informal).

Any procedure for bipartitioning the circuit $C$ through wire cutting
necessarily incurs a worst-case overhead of $\Omega(2^{k}/\varepsilon^{2})$.

We now explain how the procedure generalizes to multiple applications of the
identity in Lemma II.1. Consider a depth-$L$, $n$-qubit quantum circuit
comprising rounds $1,\dots,L$ of commuting quantum gates which implements the
channel $\mathcal{N}$. We say that a subset of the wires in the circuit are
parallel if and only if there exists an $i\in\\{1,\dots,L-1\\}$ such that each
wire connects a gate applied in round $j\leq i$ to a gate in round $j\geq
i+1$. In a circuit diagram, this corresponds to a set of wires which may be
bisected by a single vertical line. Suppose that the circuit may be separated
into disconnected sub-circuits by removing $\ell$ subsets of parallel wires,
and let $k_{1},\dots,k_{\ell}$ denote the number of wires in each subset. Then
by repeated application of Lemma II.1, we have

$\begin{gathered}\textnormal{Tr}(O_{f}\mathcal{N}(\rho_{0}^{\otimes n}))=\\\
\underset{z_{1},\dots,z_{\ell}}{\mathbb{E}}\left[\textnormal{Tr}(O_{f}\mathcal{N}^{\prime}_{z_{1},\ldots,z_{\ell}}(\rho_{0}^{\otimes
n}))\prod_{j=1}^{\ell}(2d_{j}+1)(-1)^{z_{j}}\right].\end{gathered}$ (5)

Here, $z_{1},\dots,z_{\ell}$ are independent Bernoulli random variables equal
to $1$ with probability $d_{j}/(2d_{j}+1)$ where $d_{j}=2^{k_{j}}$, for each
$j\in[\ell]$, and $\mathcal{N}^{\prime}_{z_{1},\dots,z_{\ell}}$ is the
modified quantum channel which is implemented by substituting the channel
$\Psi_{z_{j}}$ for the identity map at the $j^{\text{th}}$ location in the
original circuit, for each $j\in[\ell]$. For fixed $z_{1},\dots,z_{\ell}$, an
unbiased estimator of the quantity
$\textnormal{Tr}(O_{f}\mathcal{N}^{\prime}_{z_{1},\dots,z_{\ell}}(\rho_{0}^{\otimes
n}))$ can be obtained by measuring the output wires of the modified circuit in
the computational basis. Furthermore, if removal of the wires in each of the
$\ell$ subsets separates the circuit into $r$ disconnected fragments, then the
circuit can be executed using $r$ quantum devices, each possibly with fewer
qubits than the original. See Appendix D for further details on this point.

This setting introduces the following drawbacks: (i) the overhead is
comparatively large when $k_{j}=1$ for every $j\in[\ell]$ (as opposed to when
$\ell$ is small and each $k_{j}$ is large) and (ii) classical communication
between multiple separate devices is required. The latter issue is avoided in
Algorithm 1 since the first circuit in the sequence of circuits has no
dependence on the results from the second circuit, and so its output wires
could be measured simultaneously with the wires to be cut. These measurement
results could then be stored and used at a later time to construct the
estimator, while qubits on that same device are recycled in order to execute
the second circuit. But consider the communication graph [12] representing the
flow of communication between different circuits. Then this recycling
procedure is not possible whenever the graph contains a cycle, which may be
the case in general. In this case, the circuit sends its partial measurement
results to another, and has to wait until receiving a measurement outcome
before proceeding with its remaining gates.

##  III Circuit cutting for QAOA

In this section, we combine the above results with observations related to the
structure of QAOA circuits, and perform numerical simulations of quantum
circuit cutting applied to the QAOA.

###  III.1 Structure of QAOA circuits for circuit cutting

Here, we establish sufficient conditions for the efficient application of
quantum circuit cutting to QAOA circuits for the Max-Cut problem. A key aspect
to consider is the wires which must be cut to separate a single layer of the
QAOA circuit for Max-Cut into multiple smaller fragments. For a given input
graph $G=(V,E)$ encoding an instance of the Max-Cut problem, the cost operator
takes the form

$H_{\mathcal{C}}=\sum_{(i,j)\in E}Z_{i}Z_{j}\;,$ (6)

where $Z_{i}$ is the single-qubit Pauli-Z operator acting on the $i$-th qubit.
The Max-Cut solution corresponds to the state with the minimal expectation
value of this cost operator. One may alternatively view this operator as a
2-local Hamiltonian whose interaction graph is precisely the input graph. The
QAOA consists of applying alternating layers of mixing unitaries based on
single-qubit rotations and cost unitaries based on the above cost operator.
Since the mixing unitary comprises single-qubit gates which do not influence
qubit connectivity, this part of the circuit can safely be ignored for the
purposes of this section. We may therefore focus on the cost unitary,

$U_{\mathcal{C}}(\gamma)=e^{-i\gamma H_{\mathcal{C}}},$ (7)

which is implemented by a product of commuting two-qubit $ZZ$ rotations, one
for each edge in $E$. The fact that these operations commute also implies that
the cost unitary can be implemented using multiple different circuits,
corresponding to permutations of the terms in the product. This observation
allows us to make the simplifying assumption that gates in the QAOA circuit
are applied in an order which respects a chosen partition of the gates for the
purpose of circuit cutting, i.e., all gates in a given partition of the edges
are applied contiguously in the circuit.

In the specific case of circuit bipartitions, we make use of the idea of a
balanced vertex separator of an undirected graph $G$, which is a subset of
vertices whose removal leaves two sets of disconnected vertices of roughly
equal size, for example, of size at most $2|V|/3$.

######  Theorem III.1 (Imprecise, see restatement in Appendix D).

Suppose the graph $G=(V,E)$ with $|V|=n$ has some known balanced vertex
separator $S\subseteq V$ such that $|S|=\kappa$. Then quantum computation with
a $p$-layer QAOA circuit for Max-Cut on $G$ can be simulated to within
accuracy $\varepsilon$ using a pair of (non-entangled) quantum devices with
$n-\Omega(n)$ qubits in time $\widetilde{O}(2^{\alpha
p\kappa}/\varepsilon^{2})$ for some universal constant $\alpha$.

We once again ignore classical time to compute the objective function of a
given bit-string in the statement of this theorem. The proof is in Appendix D.
The theorem follows by combining observations about the structure of single-
layer QAOA circuits with Eq. (5) as well as an intermediate result concerning
how the number of wires to be cut scales with the number of layers $p$. We are
thus able to introduce $\Omega(n)$ “virtual qubits” in polynomial time on up
to $\log$-depth circuits, so long as the underlying graph for the QAOA has a
certain structure. In particular, this is a regime for which certain negative
results on the overall effectiveness of the QAOA do not apply. For example, it
is known that for $o(\log n)$-depth QAOA both the average- and worst-case
performance is limited compared to the best classical algorithms [40, 41, 42].

Based on the above observations, choosing suitable instances of the Max-Cut
problem for circuit cutting is related to finding graphs with small vertex
separators. Concretely, we should pick a class of graphs whose vertices are
clustered, such that the graph can be (vertex-)separated into roughly equal-
sized components. The number of vertices in the separator relates to the
number of cuts required in the overall circuit, and thus the overhead in
performing circuit cutting. We therefore suggest the class of graphs depicted
in Fig. 2 for benchmarking QAOA circuit cutting. These graphs are defined to
be those which can be written as $r$ sets of $n$ vertices connected through
$r-1$ groups of $k=O(1)$ vertices. These connecting groups of vertices form a
vertex separator of size $k(r-1)$, and the overhead for computing expectation
values scales exponentially with this quantity.

We may also consider the more general case in which there is no prior
information about the structure of the graphs available. While the problem of
finding the minimal balanced vertex separator is known to be NP-hard [43],
there is a polynomial-time classical algorithm for finding an
$O(\sqrt{\log\kappa})$ approximation to the optimal solution based on semi-
definite programming relaxation [44]. Here, $\kappa$ is the size of the
minimal vertex separator. Hence, even without prior knowledge of the minimal
vertex separator, one may employ such classical algorithms to search for
approximate solutions on the given input graph, and the overhead in Theorem
III.1 holds with $\kappa$ the minimal vertex separator, up to a factor of
$O(\sqrt{\log\kappa})$ in the exponent.

###  III.2 Sampling from cut circuits

In this section, we establish that it is possible to recover optimal solutions
from circuit fragments which have been used to perform the QAOA. We emphasize
that this sampling aspect of the algorithm has not been addressed in any
previous works on circuit cutting known to the authors. Our analysis leads to
a simple suggestion that helps bridge the gap between procedures for circuit
cutting and their application to problems of practical interest.

In the Max-Cut problem, the objective function maps bitstrings (partitions of
vertices) to their cost (number of edges cut), $f:\\{0,1\\}^{n}\to[0,M]$. Here
and throughout, we consider the graph $G=(V,E)$ with $n=|V|$ and $M=|E|\leq
n^{2}$. The goal of the QAOA is to return a bitstring $x$ such that $f(x)$ is
close to maximal, with high probability. The optimized QAOA circuit produces
the state $\left|\bm{\gamma},\bm{\beta}\right\rangle$ according to parameters
$\bm{\gamma},\bm{\beta}$. This results in the measurement distribution given
by $q(x)=|\langle x|\bm{\gamma},\bm{\beta}\rangle|^{2}$ such that
$\operatorname{{\mathbb{E}}}_{x\sim q}f(x)$ is ideally close to the maximum.
Let us denote this expectation by $\mu$ from now on, and note that
$\mu\in[0,M]$.

Given access to samples from $q$, it is straightforward to show that

$\Pr_{x\sim q}\left[x:f(x)\geq\mu\right]\geq\frac{1}{M}.$ (8)

In other words, after around $M$ trials one obtains a bit string $x$ which has
cost at least $\mu$. Now consider the case where the QAOA circuit is being
partitioned into disconnected fragments by removing $k$ wires in total. Using
Eq. (5) in the case where $\ell=k$ (individual wires are being cut), this
expression implies

$\displaystyle q(x)$ $\displaystyle=\textnormal{Tr}(|x\rangle\langle
x|\mathcal{N}(\rho_{0}^{\otimes n}))$ (9)
$\displaystyle=5^{k}\underset{z_{1},\dots,z_{k}}{\mathbb{E}}\left[\widetilde{q}(x|z_{1},\dots,z_{k})\prod_{j=1}^{k}(-1)^{z_{j}}\right]$
(10)

for every $x\in\\{0,1\\}^{n}$, where $z_{1},\dots,z_{k}$ are independent
Bernoulli random variables which determine the settings in a modified quantum
circuit, as described in Sec. II, and $\widetilde{q}(x|z_{1},\dots,z_{k})$ is
the probability of obtaining the outcome $x$ from the modified circuit
conditioned on fixed settings $z_{1},\dots,z_{k}$. Considering the absolute
value of each of the terms in the expectation, the expectation in the right-
hand side of Eq. (10) is at most

$\widetilde{q}(x):=\underset{z_{1},\dots,z_{k}}{\mathbb{E}}\left[\widetilde{q}(x|z_{1},\dots,z_{k})\right].$
(11)

This is the marginal probability of observing the outcome $x\in\\{0,1\\}^{n}$
in the modified circuit. Therefore, it holds that

$\widetilde{q}(S)\geq q(S)/5^{k}\quad\forall S\subseteq\\{0,1\\}^{n},$ (12)

where $q(S):=\sum_{x\in S}q(x)$ and likewise for the distribution
$\widetilde{q}$. Hence, using Eq. (8) we get

$\Pr_{x\sim\widetilde{q}}\left[x:f(x)\geq\mu\right]\geq\frac{1}{5^{k}M},$ (13)

which in turn implies one would have to wait at most $5^{k}$ times as long in
expectation to sample a bit string $x$ such that $f(x)\geq\mu$ using the
smaller circuit fragments. This is analogous to the overhead incurred for
expectation values: in both cases, the cost of circuit cutting is a need to
repeatedly execute the circuit by an amount that grows exponentially in the
number of wires cut.

We remark that repeating the above analysis using the wire cutting procedure
derived in the proof of Theorem 1 in Ref. [12] results in an improved overhead
of $4^{k}$ for the sampling task, and as mentioned previously does not require
synchronizing preparations to measurement outcomes in the modified circuits.

Figure 2:  An example problem graph used as input for QAOA simulations. This
graph contains $r=3$ clusters, $n=20$ nodes within each cluster and $k=1$
vertex separators between each cluster. Note that the cluster graphs are not
complete; an edge exists between any two nodes in a cluster with probability
of $0.7$. An edge exists between a node in a cluster and its neighbouring
vertex separator with probability $0.3$.

###  III.3 Numerical comparison for circuit bipartitions

Figure 3:  Convergence of the QAOA cost expectation as a function of the
number of shots. In both insets, yellow nodes indicate the wires which are
cut. Shaded intervals represent one standard deviation from the mean. Left:
$p=1$ layer, larger graph (13 qubits). Right: $p=2$ layers, medium-sized graph
(9 qubits).

Here we numerically benchmark our circuit cutting method against the methods
given in Ref. [12], and compare to the theoretical bounds. The supporting code
used to generate these numerical results can be found in [45]. We focus on a
family of graphs described in Sec. III.1, with specific examples depicted in
Fig. 3. These graphs can later be generalized, as seen on Fig. 2. The cost
variances for these graphs are examined at QAOA depth $p=1$ and $p=2$ as
functions of the number of required device shots. The results can be seen on
Fig. 3. All reference cost values have been calculated numerically exactly due
to relatively small test graphs (maximum 13 qubits).

For all graphs considered, QAOA parameters [1]
$\bm{\gamma}=\left(\gamma_{1},\ldots,\gamma_{p}\right)$ and
$\bm{\beta}=\left(\beta_{1},\ldots,\beta_{p}\right)$ have been fixed to their
optimal values

$\bm{\gamma}^{*},\bm{\beta}^{*}=\operatorname*{arg\,min}\left\langle\bm{\gamma},\bm{\beta}\right|H_{\mathcal{C}}\left|\bm{\gamma},\bm{\beta}\right\rangle,\;$
(14)

where $\left|\bm{\gamma},\bm{\beta}\right\rangle$ is the QAOA output state
given labeled parameters.

One of the main advantages of the circuit cutting method proposed in this work
is the faster convergence to the correct value of the target observable. To
meet the requirement of the cost function being in the interval $[-1,1]$, we
rescale the expression given in Eq. (6) by $M=|E|$, the number of edges in the
graph. Simulations were performed using PennyLane [33], an open-source Python
software framework for quantum differentiable programming, for both cutting
methods.

As discussed in Appendix D, a QAOA circuit with $p$ layers will generally
require at least $p$ separate cuts, one per QAOA layer. For $p\geq 2$, the
measure-and-prepare nature of both methods introduces mid-circuit measurements
that are not commonly supported on qubit simulators. To accommodate mixed
quantum states that would be generated by circuit cutting on a real quantum
device, one can either employ a mixed-state simulator or introduce auxiliary
qubits. Our implementation of the randomized channel method uses the former,
while the PennyLane implementation of the Pauli-based method of Peng et al.
[12] employs the latter.

The randomized method introduced in Sec. II outperforms the method based on
Pauli measurements [12] in terms of convergence speed and accuracy for all
considered graph sizes, cut sizes, and QAOA depths. Improvements range from
marginal to orders of magnitude. For larger cut sizes $k$, and especially at
QAOA depth $p=2$, at $10^{6}$ shots, the randomized method offers reliable
cost estimates, while the Pauli method still exhibits standard deviations
larger than the target interval $[-1,1]$ as shown in Fig. 3. We observe that
the randomized channel approach converges faster in all QAOA instances
checked, as shown on Fig. 4.

Figure 4:  Numerical experiment results for QAOA cost variance as a function
of the cut size at $10^{6}$ shots. Exponential scaling agrees with the bounds
presented in Sec. II. The parenthesized number indicates the cut size for a
given method at $p=1$.

###  III.4 Large-scale simulations

Figure 5:  Execution time vs. the number of qubits for an evaluation of a QAOA
circuit using circuit cutting. The simulation is distributed over $10$ GPU-
nodes (4 GPUs per node) and the input problem graph uses the fixed parameters
$p=2$, $n=25$, and $k=1$ while the number of clusters $r$ is varied to
increase the number of qubits. An analytic line is also included based upon
Eq. (67), where the number of unique sub-circuits $N(r)$ has been scaled by an
inferred constant execution time of $1.325$ seconds per circuit.

Figure 6:  Cost vs. steps for a full 62-qubit QAOA optimization using circuit
cutting. The input problem graph parameters are fixed at $n=20$, $r=3$, and
$k=1$. For $p=1$ QAOA layers, the optimization was run over 2 GPU-nodes and
took around 30 minutes, converging to the global minimum of $-33.8$. For
$p=2$, the optimization was run over 10 GPU-nodes and took around 12 hours. In
this case, a local minimum is obtained, although this may be improved upon by
careful choice of initial circuit parameters [46].

Having demonstrated the faster convergence of our circuit cutting method over
existing approaches for small-sized QAOA problems, we now switch focus toward
showing the applicability of a circuit cutting procedure to larger problem
sizes of above 50 qubits, a scale that is challenging for direct execution on
near-term quantum hardware devices. Our objective is to provide a proof-of-
principle evaluation and optimization of a large-scale QAOA problem when cut
into smaller fragments of at most 30 qubits. To overcome the limitations of
existing 30-qubit quantum hardware, we use a cluster of GPU-based simulators
as a substitute.

This section uses the circuit cutting method introduced in Theorem 2 of Peng
et al. [12, 16], which involves performing process tomography for each circuit
fragment and contracting the resulting tensors. The tensor-contraction-based
method is fully supported in PennyLane [33] and is compatible with both
simulator- and hardware-based devices. When performed using simulation, this
method is analogous to existing tensor network methods that aim to optimize
the order of contractions in the network [47].

As before, we focus on instances from the class of graphs introduced in Sec.
III.1 for the Max-Cut problem, consisting of $r$ clusters of $n$ nodes
connected in a chain by smaller clusters of $k$ nodes, as exemplified in Fig.
2. The choice of problem parameters $r$, $n$ and $k$ as well as the number of
QAOA layers $p$ determines the total number of unique sub-circuits $N$
required for circuit cutting as well as the maximum number of qubits $m$ among
all sub-circuits, as discussed in Appendix E.

The simulations detailed here were performed on the NERSC Perlmutter
supercomputer using a cluster of NVIDIA A100 40GB GPUs running PennyLane’s
cuQuantum-enabled [48] Lightning-GPU simulator device. Each GPU node in
Perlmutter is equipped with 4 NVIDIA A100 40GB GPUs, an AMD EPYC 7763 CPU, and
networked with the HPE Slingshot-10 interconnect. Each GPU supports executing
circuits of up to $m=30$ qubits. The $N$ circuit fragments generated due to
cutting were distributed within the cluster using the Ray parallel execution
library [49]. Supporting utility functions, data and methods used for this
section are available in [45].

We first detail the execution of QAOA circuits with parameters $p=2$, $n=25$
and $k=1$ for a range of cluster numbers $r\in\\{3,4,5\\}$, resulting in
circuits with $77$, $103$, and $129$ qubits, respectively. Such circuits could
not be executed directly on a $30$-qubit device, but by using circuit cutting
we can break the circuit up into multiple smaller circuit fragments of number
$N(r)$, depending on the choice of $r$. In this case, $N(3)=1812$, $N(4)=3540$
and $N(5)=5268$, growing linearly with $r$ as can be seen in Appendix E. To
confirm the linear scaling numerically, we used simulation-based circuit
cutting to evaluate the expectation value of the QAOA circuit with a randomly
selected set of circuit parameters $\bm{\gamma}$ and $\bm{\beta}$ for all
choices of $r$ and recorded the execution time, resulting in the plot shown in
Fig. 6. An analysis of the classical efficiency of our distributed GPU-based
simulation is provided in Appendix F.

We also investigate the ability to optimize a large-scale QAOA problem using
simulated circuit cutting. Given the ability to evaluate a circuit, it is
straightforward to extend access to gradients using methods like parameter
shift [50] or finite difference, which require repeated circuit evaluation. We
optimize a 62-qubit QAOA problem using the gradient descent algorithm,
providing gradients by combining circuit-cutting-inspired parallelization with
the multi-parameter finite-difference method. Using the class of clustered
graphs with parameters $n=20$, $r=3$ and $k=1$, we consider QAOA circuits with
$p=1$ and $p=2$ layers.

In the case of $p=1$, an exact expression for the QAOA cost is available [51,
52], from which we are able to obtain the analytic optimal value to use as a
baseline comparison. It is shown in Fig. 6 that an optimization using
simulated circuit cutting is able to achieve this optimal cost value. For
$p=2$ there is no analytic expression known for the cost function, making the
global minimum harder to find. Nevertheless, the cutting-based optimization
results in a decreasing cost function and tends towards a value of $-0.51$.
Although the result does not achieve a global minimum, since the $p=1$ case is
able to reach a lower value of $-33.8$, this provides a demonstration of the
ability to perform optimization using circuit cutting methods at larger depths
$p$.

##  IV Conclusion and outlook

In this work, we presented a circuit cutting method based on randomized
measurements that provides a quadratic runtime improvement over the current
state-of-the-art for circuits where multiple neighbouring wires are cut
simultaneously. Our method requires classical communication between circuit
fragments to coordinate measurement outcomes and state preparation. For
circuits with matrix-product-state structure, the synchronization between
circuit fragments can be carried out by repeatedly executing smaller circuits,
even on a single device. For general circuits, our algorithm requires multiple
circuit executions across many devices, separated in space rather than time.

Our results raise the question of whether additional improvements may be
possible in this setting. With this in mind, we derived an information-
theoretic lower bound that is quadratically lower than our method. It is an
open question whether such a lower bound is in fact achievable or whether it
can be further tightened.

We chose the QAOA algorithm as a testbed for the new method. An emphasis was
put on the practical and operational aspects for real devices – we give an
algorithm based on randomized measure-and-prepare channels as well as a way to
sample the cut circuit. To the best of our knowledge, the latter has not been
addressed in any previous works, which instead focus on estimating expectation
values. Numerical simulations indicate a large speedup over previous state-of-
the-art for circuit cutting.

Finally, as a more general exploration of circuit cutting methods, we give
results on large-scale simulations of QAOA circuits over many GPUs acting as
proxies for individual $\sim$ 30-qubit devices. We performed forward passes of
up to 129 qubits and full optimization procedures for 62-qubit circuits. These
results demonstrate that our software implementation of circuit cutting
methods, built on top of PennyLane, indeed enables small-scale quantum devices
to successfully emulate the results of a large circuit. Although additional
difficulties may arise when employing quantum hardware instead of simulators,
these numerical experiments are a testament to the practicality of large-scale
circuit cutting workflows.

## Acknowledgements

We thank Zeyue Niu, Zain H. Saleem, Michael A. Perlin, and Yuri Alexeev for
useful discussions.

This work was supported under DARPA project HR0011-21-9-0073, Quantum
Compilation of Unitaries and Tensors and used resources of the National Energy
Research Scientific Computing Center (NERSC), a U.S. Department of Energy
Office of Science User Facility located at Lawrence Berkeley National
Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award
DDR-ERCAP0022246 under the QIS@Perlmutter program. MM acknowledges support
from the CCQ graduate fellowship in computational quantum physics. The
Flatiron Institute is a division of the Simons Foundation.

## References

  * Farhi  _et al._ [2014a] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014a), arXiv:1411.4028 [quant-ph] .
  * Farhi and Harrow [2019] E. Farhi and A. W. Harrow, Quantum supremacy through the quantum approximate optimization algorithm (2019), arXiv:1602.07674 [quant-ph] .
  * Hadfield  _et al._ [2019] S. Hadfield, Z. Wang, B. O’Gorman, E. G. Rieffel, D. Venturelli, and R. Biswas, From the quantum approximate optimization algorithm to a quantum alternating operator ansatz, Algorithms 12, 10.3390/a12020034 (2019).
  * Peruzzo  _et al._ [2013] A. Peruzzo _et al._ , A variational eigenvalue solver on a quantum processor. eprint, arXiv preprint arXiv:1304.3061 (2013).
  * Schuld  _et al._ [2020] M. Schuld, A. Bocharov, K. M. Svore, and N. Wiebe, Circuit-centric quantum classifiers, Physical Review A 101, 032308 (2020).
  * Lloyd  _et al._ [2020] S. Lloyd, M. Schuld, A. Ijaz, J. Izaac, and N. Killoran, Quantum embeddings for machine learning, arXiv preprint arXiv:2001.03622 (2020).
  * Markov and Shi [2008] I. L. Markov and Y. Shi, Simulating quantum computation by contracting tensor networks, SIAM Journal on Computing 38, 963 (2008).
  * Bravyi and Gosset [2016] S. Bravyi and D. Gosset, Improved classical simulation of quantum circuits dominated by clifford gates, Phys. Rev. Lett. 116, 250501 (2016).
  * Bravyi  _et al._ [2019] S. Bravyi, D. Browne, P. Calpin, E. Campbell, D. Gosset, and M. Howard, Simulation of quantum circuits by low-rank stabilizer decompositions, Quantum 3, 181 (2019).
  * Bravyi  _et al._ [2016] S. Bravyi, G. Smith, and J. A. Smolin, Trading classical and quantum computational resources, Phys. Rev. X 6, 021043 (2016).
  * Pashayan  _et al._ [2015] H. Pashayan, J. J. Wallman, and S. D. Bartlett, Estimating outcome probabilities of quantum circuits using quasiprobabilities, Phys. Rev. Lett. 115, 070501 (2015).
  * Peng  _et al._ [2020] T. Peng, A. W. Harrow, M. Ozols, and X. Wu, Simulating large quantum circuits on a small quantum computer, Physical Review Letters 125, 10.1103/physrevlett.125.150504 (2020).
  * Mitarai and Fujii [2021a] K. Mitarai and K. Fujii, Overhead for simulating a non-local channel with local channels by quasiprobability sampling, Quantum 5, 388 (2021a).
  * Mitarai and Fujii [2021b] K. Mitarai and K. Fujii, Constructing a virtual two-qubit gate by sampling single-qubit operations, New Journal of Physics 23, 023021 (2021b).
  * Piveteau and Sutter [2022] C. Piveteau and D. Sutter, Circuit knitting with classical communication (2022).
  * Perlin  _et al._ [2021] M. A. Perlin, Z. H. Saleem, M. Suchara, and J. C. Osborn, Quantum circuit cutting with maximum-likelihood tomography, npj Quantum Information 7, 10.1038/s41534-021-00390-6 (2021).
  * Wiersema  _et al._ [2022] R. Wiersema, L. Guerini, J. F. Carrasquilla, and L. Aolita, Circuit connectivity boosts by quantum-classical-quantum interfaces (2022).
  * Avron  _et al._ [2021] J. Avron, O. Casper, and I. Rozen, Quantum advantage and noise reduction in distributed quantum computing, Phys. Rev. A 104, 052404 (2021).
  * Ayral  _et al._ [2020] T. Ayral, F.-M. Le Régent, Z. Saleem, Y. Alexeev, and M. Suchara, Quantum divide and compute: Hardware demonstrations and noisy simulations, in  _2020 IEEE Computer Society Annual Symposium on VLSI (ISVLSI)_ (2020) pp. 138–140.
  * Dunjko  _et al._ [2018] V. Dunjko, Y. Ge, and J. I. Cirac, Computational speedups using small quantum devices, Phys. Rev. Lett. 121, 250501 (2018).
  * Horodecki  _et al._ [2003] M. Horodecki, P. W. Shor, and M. B. Ruskai, Entanglement breaking channels, Reviews in Mathematical Physics 15, 629 (2003), https://doi.org/10.1142/S0129055X03001709 .
  * Guţă  _et al._ [2020] M. Guţă, J. Kahn, R. Kueng, and J. A. Tropp, Fast state tomography with optimal error bounds, Journal of Physics A: Mathematical and Theoretical 53, 204001 (2020).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020).
  * Chen  _et al._ [2021] S. Chen, W. Yu, P. Zeng, and S. T. Flammia, Robust shadow estimation, PRX Quantum 2, 030348 (2021).
  * Elben  _et al._ [2022] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox (2022).
  * Huggins  _et al._ [2019] W. Huggins, P. Patil, B. Mitchell, K. B. Whaley, and E. M. Stoudenmire, Towards quantum machine learning with tensor networks, Quantum Science and Technology 4, 024001 (2019).
  * Farhi  _et al._ [2014b] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm applied to a bounded occurrence constraint problem, arXiv (2014b).
  * Farhi and Harrow [2016] E. Farhi and A. W. Harrow, Quantum supremacy through the quantum approximate optimization algorithm, arXiv (2016).
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Theory of trotter error with commutator scaling, Physical Review X 11, 10.1103/physrevx.11.011020 (2021).
  * Li  _et al._ [2021] J. Li, M. Alam, and S. Ghosh, Large-scale quantum approximate optimization via divide-and-conquer (2021).
  * Saleem  _et al._ [2021] Z. H. Saleem, T. Tomesh, M. A. Perlin, P. Gokhale, and M. Suchara, Quantum divide and conquer for combinatorial optimization and distributed computing (2021), arXiv:2107.07532 [quant-ph] .
  * Tang and Martonosi [2022] W. Tang and M. Martonosi, Scaleqc: A scalable framework for hybrid computation on quantum and classical processors (2022).
  * Bergholm  _et al._ [2018] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, M. S. Alam, S. Ahmed, J. M. Arrazola, C. Blank, A. Delgado, S. Jahangiri, _et al._ , Pennylane: Automatic differentiation of hybrid quantum-classical computations, arXiv:1811.04968 (2018).
  * Bravyi  _et al._ [2021] S. Bravyi, D. Gosset, and R. Movassagh, Classical algorithms for quantum mean values, Nature Physics 17, 337 (2021).
  * Kueng and Gross [2015] R. Kueng and D. Gross, Qubit stabilizer states are complex projective 3-designs (2015).
  * Webb [2016] Z. Webb, The clifford group forms a unitary 3-design, Quantum Info. Comput. 16, 1379–1400 (2016).
  * Zhu [2017] H. Zhu, Multiqubit clifford groups are unitary 3-designs, Phys. Rev. A 96, 062336 (2017).
  * Van Den Berg [2021] E. Van Den Berg, A simple method for sampling random clifford operators, in  _2021 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (2021) pp. 54–59.
  * Aaronson and Gottesman [2004] S. Aaronson and D. Gottesman, Improved simulation of stabilizer circuits, Phys. Rev. A 70, 052328 (2004).
  * Bravyi  _et al._ [2020] S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, Obstacles to variational quantum optimization from symmetry protection, Phys. Rev. Lett. 125, 260505 (2020).
  * Farhi  _et al._ [2020a] E. Farhi, D. Gamarnik, and S. Gutmann, The quantum approximate optimization algorithm needs to see the whole graph: Worst case examples (2020a), arXiv:2005.08747 [quant-ph] .
  * Farhi  _et al._ [2020b] E. Farhi, D. Gamarnik, and S. Gutmann, The quantum approximate optimization algorithm needs to see the whole graph: A typical case (2020b), arXiv:2004.09002 [quant-ph] .
  * Bui and Jones [1992] T. N. Bui and C. Jones, Finding good approximate vertex and edge partitions is np-hard, Information Processing Letters 42, 153 (1992).
  * Feige  _et al._ [2008] U. Feige, M. Hajiaghayi, and J. R. Lee, Improved approximation algorithms for minimum weight vertex separators, SIAM Journal on Computing 38, 629 (2008), https://doi.org/10.1137/05064299X .
  * [45] https://github.com/XanaduAI/randomized-measurements-circuit-cutting, (2022).
  * Lykov  _et al._ [2022] D. Lykov, J. Wurtz, C. Poole, M. Saffman, T. Noel, and Y. Alexeev, Sampling frequency thresholds for quantum advantage of quantum approximate optimization algorithm, arXiv preprint arXiv:2206.03579 (2022).
  * Gray and Kourtis [2021] J. Gray and S. Kourtis, Hyper-optimized tensor network contraction, Quantum 5, 410 (2021).
  * Fang  _et al._ [2022] L. Fang, A. Hehn, H. Bayraktar, and S. Stanwyck, Nvidia/cuquantum: cuquantum v22.05.0 (2022).
  * Moritz  _et al._ [2018] P. Moritz, R. Nishihara, S. Wang, A. Tumanov, R. Liaw, E. Liang, M. Elibol, Z. Yang, W. Paul, M. I. Jordan, _et al._ , Ray: A distributed framework for emerging $\\{$AI$\\}$ applications, in  _13th USENIX Symposium on Operating Systems Design and Implementation (OSDI 18)_ (2018) pp. 561–577.
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Evaluating analytic gradients on quantum hardware, Physical Review A 99, 032331 (2019).
  * Wang  _et al._ [2018] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, Quantum approximate optimization algorithm for MaxCut: A fermionic view, Phys. Rev. A 97, 022304 (2018).
  * Medvidović and Carleo [2021] M. Medvidović and G. Carleo, Classical variational simulation of the quantum approximate optimization algorithm, npj Quantum Information 7, 1 (2021).
  * Watrous [2018] J. Watrous,  _The Theory of Quantum Information_ (Cambridge University Press, 2018).
  * Cover and Thomas [2005] T. M. Cover and J. A. Thomas,  _Elements of Information Theory_ (Wiley, 2005).
  * Buscemi and Datta [2010] F. Buscemi and N. Datta, The quantum capacity of channels with arbitrarily correlated noise, IEEE Transactions on Information Theory 56, 1447 (2010).
  * Sason and Verdú [2016] I. Sason and S. Verdú, $f$ -divergence inequalities, IEEE Transactions on Information Theory 62, 5973 (2016).
  * Fano [1966] R. M. Fano,  _Transmission of Information: a Statistical Theory of Communications_ (MIT Press, 1966).
  * Haah  _et al._ [2017] J. Haah, A. W. Harrow, Z. Ji, X. Wu, and N. Yu, Sample-optimal tomography of quantum states, IEEE Transactions on Information Theory 63, 5628 (2017).
  * Lowe [2021] A. Lowe,  _Learning Quantum States Without Entangled Measurements_ , Master’s thesis (2021).
  * Shoukourian  _et al._ [2014] H. Shoukourian, T. Wilde, A. Auweter, and A. Bode, Predicting the energy and power consumption of strong and weak scaling HPC applications, Supercomputing Frontiers and Innovations 1, 20 (2014).

##  Appendix A Preliminaries

We explicitly define random variables throughout, including matrix-valued
random variables where relevant. We let $x\in A$ denote a random variable $x$
which takes values in the set $A$. For any positive integer $N$ we let $[N]$
denote the set $\\{1,\dots,N\\}$. We use sans-serif font to denote subspaces
of complex finite-dimensional Euclidean vector spaces, or operators acting on
these spaces. For any positive integer $d$ we let $\mathsf{L}(\mathbb{C}^{d})$
denote the set of square linear operators acting on $\mathbb{C}^{d}$,
$\mathsf{H}(\mathbb{C}^{d})$ the subset of operators in
$\mathsf{L}(\mathbb{C}^{d})$ which are Hermitian,
$\mathsf{Psd}(\mathbb{C}^{d})$ the subset of operators in
$\mathsf{H}(\mathbb{C}^{d})$ which are positive semidefinite, and
$\mathsf{D}(\mathbb{C}^{d})$ the subset of operators in
$\mathsf{Psd}(\mathbb{C}^{d})$ which have unit trace (quantum states). We also
let $\mathsf{U}(\mathbb{C}^{d})$ denote the set of unitary operators acting on
$\mathbb{C}^{d}$ and $\mathsf{S}(\mathbb{C}^{d})$ the set of unit-norm vectors
in $\mathbb{C}^{d}$ (vector-representation of pure states). We distinguish
between identity operators acting on $\mathbb{C}^{d}$ and on
$\mathsf{L}(\mathbb{C}^{d})$ by adopting the notation $\mathds{1}$ for the
former and $\mathrm{id}$ for the latter, using subscripts to denote the space
on which the operators act if this is not clear from context. Our main theorem
will make use of structured POVMs based on unitary $t$-designs (in the case
where $t=2$), which we define below.

######  Definition A.1 (Unitary $t$-design).

For positive integers $t,d>0$, a _unitary $t$-design_ is a random unitary
operator $U\in\mathsf{U}(\mathbb{C}^{d})$ which for any
$X\in\mathsf{L}(\mathbb{C}^{d})^{\otimes t}$ satisfies

$\displaystyle\operatorname{{\mathbb{E}}}_{U}\left[U^{\otimes
t}X(U^{\dagger})^{\otimes
t}\right]=\int_{\mathsf{U}(\mathbb{C}^{d})}V^{\otimes
t}X(V^{\dagger})^{\otimes t}\ \mathrm{d}\mu(V)$ (15)

where $\mu$ is the Haar measure on $\mathsf{U}(\mathbb{C}^{d})$.

In the case where $U$ in the above definition is a discrete random variable
taking values in
$\\{U_{1},U_{2}\dots,U_{m}\\}\subset\mathsf{U}(\mathbb{C}^{d})$ with
probabilities $p_{1},p_{2},\dots,p_{m}$, respectively, one may associate with
this unitary $t$-design the ensemble of unitaries
$\\{(p_{i},U_{i})\\}_{i=1}^{m}$. We denote by $\textnormal{Tr}_{j}(X)$ the
partial trace over the $j^{\text{th}}$ space, where it should be made clear
from context which subspace the index $j$ corresponds to. We also use the
notation $\textnormal{Tr}_{S}(\cdot)$ for a set $S\subset\mathbb{Z}_{+}$ to
denote tracing out a subspace corresponding to multiple qubits on a quantum
circuit.

##  Appendix B Proof of main result

###  B.1 Proof of Lemma II.1

The lemma is restated for convenience. Recall that we have

$\displaystyle\Psi_{0}(X)=\operatorname{{\mathbb{E}}}_{U}\left[\sum_{j=1}^{d}\left\langle
j\right|U^{\dagger}XU\left|j\right\rangle U|j\rangle\langle
j|U^{\dagger}\right]$ (2)

for every $X\in\mathsf{L}(\mathbb{C}^{d})$, where
$U\in\mathsf{U}(\mathbb{C}^{d})$ is a unitary 2-design, and

$\displaystyle\Psi_{1}(X)=\textnormal{Tr}(X)\mathds{1}/d\quad$ (3)

for every $X\in\mathsf{L}(\mathbb{C}^{d})$.

######  Lemma II.1.

Let $d$ be a positive integer and $\Psi_{0}$, $\Psi_{1}$ be the channels
defined in Eqs. (2) and (3), respectively, acting on $d$-dimensional states.
Define the Bernoulli random variable $z\in\\{0,1\\}$ to be equal to $1$ with
probability $d/(2d+1)$. It holds that

$\displaystyle\mathrm{id}=(2d+1)\operatorname{{\mathbb{E}}}_{z}\left[\
(-1)^{z}\ \Psi_{z}\right].$ (16)

###### Proof.

Since we assume that $U\in\mathsf{U}(\mathbb{C}^{d})$ is a unitary 2-design,
for any $\left|v\right\rangle\in\mathsf{S}(\mathbb{C}^{d})$ we have

$\displaystyle\operatorname{{\mathbb{E}}}_{U}\left[(U|v\rangle\langle
v|U^{\dagger})^{\otimes 2}\right]$
$\displaystyle=\int_{\mathsf{U}(\mathbb{C}^{d})}(V|v\rangle\langle
v|V^{\dagger})^{\otimes 2}\ \mathrm{d}\mu(V)$ (17)
$\displaystyle=\frac{1}{d(d+1)}\left(\mathds{1}\otimes\mathds{1}+W\right)$
(18)

where $W\in\mathsf{U}((\mathbb{C}^{d})^{\otimes 2})$ is the swap operator, and
the second line is a well-known identity that follows by exploiting the
permutation-invariance of Haar integrals of this form (cf. Eq. (7.179) in
[53]). Using Eq. (18) as well as the linearity of trace, we may write the
action of the first channel $\Psi_{0}$ on a linear operator
$X\in\mathsf{L}(\mathbb{C}^{d})$ as

$\displaystyle\Psi_{0}(X)$
$\displaystyle=\operatorname{{\mathbb{E}}}_{U}\left[\sum_{j=1}^{d}\left\langle
j\right|U^{\dagger}XU\left|j\right\rangle U|j\rangle\langle
j|U^{\dagger}\right]$ (19)
$\displaystyle=\textnormal{Tr}_{1}\left((X\otimes\mathds{1})\sum_{j=1}^{d}\operatorname{{\mathbb{E}}}_{U}\left[(U|j\rangle\langle
j|U^{\dagger})^{\otimes 2}\right]\right)$ (20)
$\displaystyle=\frac{1}{d+1}\left[\textnormal{Tr}_{1}\left((X\otimes\mathds{1})\right)+\textnormal{Tr}_{1}\left((X\otimes\mathds{1})W\right)\right]$
(21) $\displaystyle=\frac{1}{d+1}\left(\textnormal{Tr}(X)\mathds{1}+X\right),$
(22)

where the final line follows from the identity
$\textnormal{Tr}_{1}((A\otimes\mathds{1})W)=A$ for any square linear operator
$A$. Also, $\Psi_{1}(X)=\textnormal{Tr}(X)\mathds{1}/d$ for every
$X\in\mathsf{L}(\mathbb{C}^{d})$ by definition. Substituting into Eq. (22), we
have

$\displaystyle\Psi_{0}(X)=\frac{1}{d+1}\left[d\Psi_{1}(X)+X\right]\quad\forall
X\in\mathsf{L}(\mathbb{C}^{d})$ (23)

which, upon rearranging, gives

$\displaystyle X$ $\displaystyle=(d+1)\Psi_{0}(X)-d\Psi_{1}(X)$ (24)
$\displaystyle=(2d+1)\left(\frac{d+1}{2d+1}\Psi_{0}(X)-\frac{d}{2d+1}\Psi_{1}(X)\right)$
(25)
$\displaystyle=(2d+1)\operatorname{{\mathbb{E}}}_{z}\left[(-1)^{z}\Psi_{z}(X)\right]$
(26)

for every $X\in\mathsf{L}(\mathbb{C}^{d})$. This proves the claim. ∎

###  B.2 Bipartitioning pseudocode

We describe in full detail the circuit cutting procedure used in Theorem II.2
using the pseudocode below, in Algorithm 1. As in the statement of the
theorem, $C$ is an $n$-qubit circuit composed of sub-circuits $C_{A}$, $C_{B}$
acting on $A,B\subseteq[n]$ respectively, and $|A\cap B|=k$. Recalling our
notation for the model of quantum computation introduced in Sec. II, for any
function $f:\\{0,1\\}^{n}\to[-1,1]$ we have
$O_{f}=\sum_{x\in\\{0,1\\}^{n}}f(x)|x\rangle\langle x|$. The overall action of
the quantum circuit $C$ is to implement a quantum channel which we denote by
$\mathcal{N}$, and the goal is to estimate
$\textnormal{Tr}(O_{f}\mathcal{N}(\rho_{0}^{\otimes n}))$. Note also that we
assume that $A\cup B=[n]$ for simplicity.

Algorithm 1 Circuit cutting procedure in Theorem II.2.

Input: $f:\\{0,1\\}^{n}\to[-1,1]$, quantum circuit $C$ as in Theorem II.2,
accuracy parameter $\varepsilon$

Output: Random variable $Y$ s.t.
$\operatorname{{\mathbb{E}}}[Y]=\textnormal{Tr}(O_{f}\mathcal{N}(\rho_{0}^{\otimes
n}))$ and $|Y|\leq O(2^{k})$ with certainty

1:$z\leftarrow 1$ with probability $2^{k}/(2^{k+1}+1)$, $0$ otherwise

2:Initialize wires $A$ to $\left|0\right\rangle_{A}$

3:Apply circuit $C_{A}$ to wires $A$

4:if $z=0$ then

5: $V\leftarrow$ random Clifford operator on $k$ qubits

6: Apply circuit $V^{\dagger}$ to wires $A\cap B$

7: $y\leftarrow$ measurement of wires $A\cap B$

8:else

9: $y\leftarrow\text{Unif}(\\{0,1\\}^{k})$

10:end if

11:$x_{A}\leftarrow$ measurement of wires $[n]\backslash B$

12:Initialize wires $B$ to $\left|y\right\rangle_{A\cap
B}\left|0\right\rangle_{[n]\backslash A}$

13:Apply circuit $C_{B}$ to wires $B$

14:$x_{B}\leftarrow$ measurement of wires $B$

15:$x\leftarrow x_{A}x_{B}$

16:$Y\leftarrow(2^{k+1}+1)(-1)^{z}f(x)$

17:return $Y$

###  B.3 Proof of Theorem II.2

We first restate the result for convenience.

######  Theorem II.2.

Let $C$ be a size-$m$ quantum circuit acting on $n$ qubits which is a
composition of circuits $C_{A},C_{B}$ acting non-trivially on sets of qubits
$A,B\subseteq[n]$, respectively. If $|A\cap B|\leq k$, then quantum
computation using $C$ can be simulated to within accuracy $\varepsilon$ in
time $O(4^{k}(m+k^{2})/\varepsilon^{2})$ by a quantum circuit acting on at
most $\max\\{|A|,|B|\\}$ qubits using the procedure described in Algorithm 1.

###### Proof.

Define $d:=2^{k}$. To prove the claim regarding time and space resources
required, first note that the random variable $Y=(2^{k+1}+1)(-1)^{z}f(x)$ in
Algorithm 1 is bounded in magnitude by $|Y|\leq 2d+1$ with certainty. Hence by
Hoeffding’s Inequality the sample mean of
$N=O(d^{2}/\varepsilon^{2})=O(4^{k}/\varepsilon^{2})$ iterations of Algorithm
1 suffices to estimate the expectation value of $Y$ to within accuracy
$\varepsilon$ with high probability. Now, suppose we have a device comprising
$\max\\{|A|,|B|\\}$ qubits. If $z=0$, the procedure up to and including Line
11 can be performed on this device in time $O(m+k^{2})$. This follows from the
fact that there exists an efficient procedure to sample a depth-$O(k\log(k))$
circuit which implements a random Clifford operator, running in time
$O(k^{2})$ [38]. If $z=1$ then no sampling of random Cliffords is required,
and the procedure yields $x_{A}$ and $y$ in time $O(m)$. In Line 12, the same
device may then be re-initialized to the state
$\left|y\right\rangle\left|0\right\rangle$. Finally, the application of the
circuit $C_{B}$ takes time $O(m)$ for a total runtime of $O(m+k^{2})$ to
produce a single sample.

It remains to show that $Y$ is an unbiased estimator of
$\textnormal{Tr}(O_{f}\mathcal{N}(\rho_{0}^{\otimes n}))$, where $\mathcal{N}$
is the channel implemented by the circuit $C$. Consider the channels
$\Psi_{0}$ and $\Psi_{1}$ defined in Eqs. (2) and (3), respectively. Lines 5-7
and 12 perform the action of the channel $\Psi_{0}$ on qubits in $A\cap B$.
This follows from the fact that the random Clifford operators comprise a
3-design [35, 37, 36]. Similarly, lines 9 and 12 perform the action of the
channel $\Psi_{1}$ on qubits in $A\cap B$. Deferring the measurement in Line
11, the procedure therefore has the following action on the initial state for
a fixed $z\in\\{0,1\\}$:

$\displaystyle\rho_{0}^{\otimes n}$
$\displaystyle\to(\mathcal{U}_{A}\otimes\mathrm{id}_{[n]\backslash
A})(\rho_{0}^{\otimes n})$ (Line 3) (27)
$\displaystyle\to(\mathrm{id}_{[n]\backslash
B}\otimes\Psi_{z}\otimes\mathrm{id}_{[n]\backslash
A})(\mathcal{U}_{A}\otimes\mathrm{id}_{[n]\backslash A})(\rho_{0}^{\otimes
n})$ (Lines 4-10) (28)
$\displaystyle\to\tau^{(z)}:=(\mathrm{id}_{[n]\backslash
B}\otimes\mathcal{U}_{B})(\mathrm{id}_{[n]\backslash
B}\otimes\Psi_{z}\otimes\mathrm{id}_{[n]\backslash
A})(\mathcal{U}_{A}\otimes\mathrm{id}_{[n]\backslash A})(\rho_{0}^{\otimes
n}).$ (Line 13) (29)

Here, $\mathcal{U}_{A}$ and $\mathcal{U}_{B}$ are unitary channels acting on
wires (qubits) $A$ and $B$ respectively, such that the overall action of the
circuit is given by

$\displaystyle\mathcal{N}(\rho_{0}^{\otimes n})=(\mathrm{id}_{[n]\backslash
B}\otimes\mathcal{U}_{B})(\mathcal{U}_{A}\otimes\mathrm{id}_{[n]\backslash
A})(\rho_{0}^{\otimes n}).$ (30)

Abusing notation and letting $z$ denote the random variable defined in Line 1
of Algorithm 1, by Lemma II.1 and the linearity of channels we have

$\displaystyle\operatorname{{\mathbb{E}}}[(2d+1)(-1)^{z}\tau^{(z)}]=\mathcal{N}(\rho_{0}^{\otimes
n}).$ (31)

Therefore, by the linearity of trace it holds that

$\displaystyle\operatorname{{\mathbb{E}}}[Y]$
$\displaystyle=\operatorname{{\mathbb{E}}}\left[\operatorname{{\mathbb{E}}}\left[(2d+1)(-1)^{z}f(x)\
|\ z\right]\right]$ (32)
$\displaystyle=\textnormal{Tr}\left(O_{f}\operatorname{{\mathbb{E}}}\left[(2d+1)(-1)^{z}\tau^{(z)}\right]\right)$
(33) $\displaystyle=\textnormal{Tr}\left(O_{f}\mathcal{N}(\rho_{0}^{\otimes
n})\right)$ (34)

as required. ∎

##  Appendix C An information-theoretic lower bound

Let us begin by formally defining what constitutes a procedure for wire
cutting. Let $\mathcal{H}_{A}$, $\mathcal{H}_{B}$, and $\mathcal{H}_{C}$ be
complex Euclidean spaces of dimension $d_{A}$, $d_{B}$, $d_{C}$ corresponding
to registers $A$, $B$, $C$, and let
$\mathcal{H}:=\mathcal{H}_{A}\otimes\mathcal{H}_{B}\otimes\mathcal{H}_{C}$. We
require that any procedure for wire cutting succeeds at the following task
using a fixed measurement strategy on the register $B$. Given (i) a unitary
channel $\mathcal{U}_{2}$ acting on states of the register $(B,C)$, (ii) a
diagonal observable $O_{f}\in\mathsf{H}(\mathcal{H})$ such that $\left\lVert
O_{f}\right\rVert\leq 1$, and (iii) an unknown state of the register $(A,B,C)$
of the form $\rho\otimes\sigma$ where
$\rho\in\mathsf{D}(\mathcal{H}_{A}\otimes\mathcal{H}_{B})$ and
$\sigma\in\mathsf{D}(\mathcal{H}_{C})$, a procedure for wire cutting first
performs a measurement on the register $B$, leaving the register $A$ in some
reduced state $\rho_{A}$. Then, the register $B$ is prepared in the state
$\tau_{y}$ conditioned on the outcome $y$ having been obtained. Next, the
unitary channel $\mathcal{U}_{2}$ is applied to the state on the composite
register $(B,C)$. Finally all registers are measured in the standard basis,
resulting in an outcome $z\in[d_{A}d_{B}d_{C}]$. A successful procedure for
wire cutting outputs an $\varepsilon$-accurate estimate of the expectation
value
$\textnormal{Tr}(O_{f}(\mathds{1}_{A}\otimes\mathcal{U}_{2})(\rho\otimes\sigma))$
using $N$ i.i.d. copies of $z$.

Note that a successful procedure for this task using $N\leq
O(\kappa^{2}/\varepsilon^{2})$ iterations is implied by the existence of a
decomposition of the identity of the form

$\displaystyle\mathrm{id}=\sum_{i}a_{i}\Phi_{i}$ (35)

for some collection of measure-and-prepare channels
$\Phi_{i}:\mathsf{L}(\mathcal{H}_{B})\to\mathsf{L}(\mathcal{H}_{B})$ and real
$a_{i}$ satisfying

$\displaystyle\sum_{i}|a_{i}|=\kappa.$ (36)

This follows from the fact that such a decomposition could be used to perform
the wire cutting task described above, resulting in an unbiased estimator of
the expectation value which is bounded in magnitude by $\kappa$. Hence, our
lower bound on $N$ will also imply a lower bound on the quantity $\kappa$.

To prove our lower bound, let us consider the case where $d_{A}=d_{C}=1$. Then
we must be able to determine $\textnormal{Tr}(O_{f}\mathcal{U}_{2}(\rho))$ to
within error $\varepsilon$ using a procedure of the form described above.
Namely, we measure $\rho\in\mathsf{D}(\mathcal{H}_{B})$ using some POVM
$\\{M_{y}\\}$ and prepare the state $\tau_{y}\in\mathsf{D}(\mathcal{H}_{B})$
conditioned on receiving the outcome $y$. We then apply $\mathcal{U}_{2}$ on
this state and measure it in the computational basis receiving outcome
$z\in[d_{B}]$, and this procedure is repeated $N$ times using $N$ identical
copies of the state $\rho$.

Since there is only one non-trivial register we will drop the subscript $B$
from now on. Define

$\displaystyle\rho_{j}:=\frac{2\varepsilon}{d}\Pi_{j}+(1-\varepsilon)\frac{\mathds{1}}{d},\quad
j=0,1$ (37)

where

$\displaystyle\Pi_{0}:=\sum_{i=1}^{d/2}|i\rangle\langle
i|,\qquad\Pi_{1}:=\mathds{1}-\Pi_{0}.$ (38)

These states will allow us to construct a difficult state discrimination
problem which reduces to a successful procedure for wire cutting. The lower
bound on the number of iterations required will then follow by bounding the
information gained from the measurement statistics in a single iteration.

To accomplish this, we make use of elementary facts from information theory
(see, for example, Ref. [54].) We adopt standard notation for the mutual
information and conditional mutual information between two random variables,
as well as the entropy and conditional entropy of a random variable. We let
$D(P\;\|\;Q)$ denote the relative entropy between two distributions $P$, $Q$
such that $\mathrm{supp}(P)\subseteq\mathrm{supp}(Q)$. We also make use of the
$\chi^{2}$-divergence between two discrete distributions $P$, $Q$ over the
sample space $\mathcal{X}$ defined through

$\displaystyle
D_{\chi^{2}}(P\;\|\;Q):=\sum_{x\in\mathcal{X}}Q(x)\left(\frac{P(x)}{Q(x)}-1\right)^{2}=\left(\sum_{x\in\mathcal{X}}\frac{P(x)^{2}}{Q(x)}\right)-1.$
(39)

In addition to more standard facts from information theory, we require the
following two results, which we state without proof. The first is a special
case of Lemma 6 in Ref. [55], while the second is a looser version of Eq. (5)
in Ref. [56].

######  Lemma C.1.

Let $P_{X,Y}$ be a discrete distribution over the sample space
$\mathcal{X}\times\mathcal{Y}$, and let $P_{X}$, $P_{Y}$ be the marginal
distributions over $\mathcal{X}$, $\mathcal{Y}$, respectively. For any pair of
discrete distributions $Q_{1}$ over $\mathcal{X}$ and $Q_{2}$ over
$\mathcal{Y}$, it holds that

$\displaystyle D(P_{X,Y}\;\|\;P_{X}\otimes P_{Y})\leq
D(P_{X,Y}\;\|\;Q_{1}\otimes Q_{2}).$ (40)

######  Lemma C.2.

Let $P$, $Q$ be discrete distributions over the sample space $\mathcal{X}$
such that $\mathrm{supp}(P)\subseteq\mathrm{supp}(Q)$. It holds that

$\displaystyle D(P\;\|\;Q)\leq\frac{1}{\ln(2)}D_{\chi^{2}}(P\;\|\;Q).$ (41)

Let $U\in\mathsf{U}(\mathcal{H})$ be a Haar-random unitary operator and
$x\in\\{0,1\\}$ be uniformly random. Suppose we are given as input to our
problem the unitary channel $\mathcal{U}:X\mapsto U^{\dagger}XU$, observable
$\Pi_{0}\in\mathsf{H}(\mathcal{H})$, and the state
$U\rho_{x}U^{\dagger}\in\mathsf{D}(\mathcal{H})$. Then the wire cutting
procedure produces standard basis outcomes $z_{1},\dots,z_{N}\in[d]$ which are
independent given $x$ and $U$, as well as intermediate measurement results
$y_{1},\dots,y_{N}$ which are independent given $x$ and $U$. The final output
of the procedure after the $N$ iterations enables one to correctly determine
the value of $x\in\\{0,1\\}$ with high probability. This follows since without
loss of generality we can assume the output will be an
$\varepsilon/3$-accurate estimate of the quantity
$\textnormal{Tr}(\Pi_{0}U^{\dagger}U\rho_{x}U^{\dagger}U)=\textnormal{Tr}(\Pi_{0}\rho_{x})$,
with high probability. Now, because

$\displaystyle\textnormal{Tr}(\Pi_{0}\rho_{0})=\frac{1}{2}+\frac{\varepsilon}{2}\qquad\textnormal{and}\qquad\textnormal{Tr}(\Pi_{0}\rho_{1})=\frac{1}{2}-\frac{\varepsilon}{2}$
(42)

such an $\varepsilon/3$-accurate estimate of
$\textnormal{Tr}(\Pi_{0}\rho_{x})$ allows one to infer the value of $x$ with
certainty. By Fano’s Inequality [57], the mutual information between the
random variables obtained from identical iterations of this procedure
$z:=(z_{1},\dots,z_{N})$ and the choice of state $x$ satisfies

$\displaystyle I(x:z)\geq\Omega(1).$ (43)

On the other hand, we can upper bound the left-hand side of the above as
follows. We have

$\displaystyle I(x:z)$ $\displaystyle\leq I(x:z,U)$ (Data-Processing
Inequality) (44) $\displaystyle=I(x:U)+I(x:z|U)$ (chain rule for mut. inf.)
(45) $\displaystyle=I(x:z_{1},\dots,z_{N}|U)$ ($x$ and $U$ indep.) (46)
$\displaystyle=\sum_{i=1}^{N}I(x:z_{i}|U,z_{i-1},\dots,z_{1})$ (chain rule for
mut. inf.) (47)
$\displaystyle=\sum_{i=1}^{N}H(z_{i}|U,z_{i-1},\dots,z_{1})-H(z_{i}|x,U,z_{i-1},\dots,z_{1})$
(defn. of mut. inf.) (48)
$\displaystyle\leq\sum_{i=1}^{N}H(z_{i}|U)-H(z_{i}|x,U)$ (49)
$\displaystyle=\sum_{i=1}^{N}I(x:z_{i}|U)$ (50)

where the second to last line follows because conditioning reduces entropy,
and $z_{i}$ is independent of $z_{i-1},\dots,z_{1}$ given $x$ and $U$ by the
assumptions made in our criteria for what constitutes a wire cutting
procedure. Hence, it remains to bound the individual mutual information terms.

Conditioning on $U$ we find $x\to y_{i}\to z_{i}$ forms a Markov Chain, using
the fact that for fixed $y_{i}$ and $U$ we defined $z_{i}$ to be the outcome
obtained from measuring the state $U^{\dagger}\tau_{y_{i}}U$ in the
computational basis. By Data-Processing Inequality we have

$\displaystyle I(x:z_{i}|U)\leq I(x:y_{i}|U).$ (51)

Define $\tilde{y}$ to be the random variable $y_{i}$ conditioned on $U$, which
has conditional distribution given by
$P_{\tilde{y}|x=x^{\prime}}(y)=\textnormal{Tr}(M_{y}U\rho_{x^{\prime}}U^{\dagger})$
for each $x^{\prime}\in\\{0,1\\}$ and marginal distribution
$P_{\tilde{y}}(y)=\sum_{x^{\prime}=0}^{1}P_{\tilde{y}|x=x^{\prime}}(y)$. Also,
let $P_{x}=(\frac{1}{2},\frac{1}{2})$, let $P_{x,\tilde{y}}$ be the joint
distribution, and define a distribution $Q$ given by
$Q(y)=\textnormal{Tr}(M_{y})/d$. Then the right-hand side of the above is

$\displaystyle\operatorname{{\mathbb{E}}}_{U}\left[D(P_{x,\tilde{y}}\;\|\;P_{x}\otimes
P_{\tilde{y}})\right]$
$\displaystyle\leq\operatorname{{\mathbb{E}}}_{U}\left[D(P_{x,\tilde{y}}\;\|\;P_{x}\otimes
Q)\right]$ (by Lemma C.1) (52)
$\displaystyle=\frac{1}{2}\operatorname{{\mathbb{E}}}_{U}\left[D(P_{\tilde{y}|x=0}\;\|\;Q)+D(P_{\tilde{y}|x=1}\;\|\;Q)\right]$
(53)
$\displaystyle=\operatorname{{\mathbb{E}}}_{U}\left[D(P_{\tilde{y}|x=0}\;\|\;Q)\right].$
(54)

In the second line we made use of the unitary invariance of the Haar measure,
which implies that the two terms in the sum in Eq. (53) are equal by the fact
that there exits a unitary operator $V\in\mathsf{U}(\mathcal{H})$ such that
$\rho_{1}=V\rho_{0}V^{\dagger}$. We then have

$\displaystyle\operatorname{{\mathbb{E}}}_{U}\left[D(P_{\tilde{y}|x=0}\;\|\;Q)\right]$
$\displaystyle\leq\frac{1}{\ln(2)}\operatorname{{\mathbb{E}}}_{U}\left[D_{\chi^{2}}(P_{\tilde{y}|x=0}\;\|\;Q)\right]$
(by Lemma C.2) (55)
$\displaystyle=\frac{1}{\ln(2)}\left[\left(\sum_{y}\frac{\operatorname{{\mathbb{E}}}_{U}\left[P_{\tilde{y}|x=0}(y)^{2}\right]}{Q(y)}\right)-1\right].$
(56)

The Haar integral in the numerator of each of the terms above is evaluated
implicitly in previous arguments for lower bounds of this type [58, 23], and
explicitly as Lemma 3.2.6 in [59], by which it holds that

$\displaystyle\operatorname{{\mathbb{E}}}_{U}\left[P_{\tilde{y}|x=0}(y)^{2}\right]\leq
Q(y)^{2}\left(1+\frac{\varepsilon^{2}}{d+1}\right)$ (57)

for every outcome $y$, and hence

$\displaystyle\frac{1}{\ln(2)}\left[\left(\sum_{y}\frac{\operatorname{{\mathbb{E}}}_{U}\left[P_{\tilde{y}|x=0}(y)^{2}\right]}{Q(y)}\right)-1\right]\leq\frac{\varepsilon^{2}}{\ln(2)(d+1)}.$
(58)

In summary, we have

$\displaystyle I(x:z_{i}|U)\leq\frac{\varepsilon^{2}}{\ln(2)(d+1)}$ (59)

for every $i\in[N]$ so that, by Eq. (50) and the lower bound in Eq. (43) we
have

$\displaystyle\Omega(1)\leq I(x:z)\leq\frac{N\varepsilon^{2}}{\ln(2)(d+1)}$
(60)

which holds if and only if $N=\Omega(d/\varepsilon^{2})$. For a system
comprising $k$ qubits, this is $\Omega(2^{k}/\varepsilon^{2})$. This proves
the bound in Theorem II.3.

As discussed previously, this implies a lower bound on an optimization problem
defined as follows. Let $\mathcal{M}\mathcal{P}(\mathbb{C}^{d})$ denote the
set of measure-and-prepare channels acting on $d$-dimensional quantum states.
Define $p^{\star}$ through:

$\displaystyle p^{\star}:=\inf_{a_{i},\Phi_{i}}$
$\displaystyle\sum_{i}|a_{i}|$ (61) s.t.
$\displaystyle\mathrm{id}-\sum_{i}a_{i}\Phi_{i}=0$
$\displaystyle\Phi_{i}\in\mathcal{M}\mathcal{P}(\mathbb{C}^{d})\quad\forall i$
$\displaystyle a_{i}\in\mathbb{R}\quad\forall i.$

It holds that $p^{\star}\geq\Omega(\sqrt{d})$.

##  Appendix D Proof of Theorem III.1

In this appendix we offer proofs for Theorem III.1 as well as some related
results in the form of lemmas. We state a notational convenience at the
outset: in order not to confuse undirected input graphs for the Max-Cut
problem with the (multi-) graphs representing quantum circuits, we reserve the
term fragments to refer to subgraphs of the latter, obtained by removing a
subset of wires. Recall that we are interested in these fragments since
applying the identity in Lemma II.1 (or in the proof of Theorem 1 of Ref.
[12]) to each of the removed wires allows us to simulate the circuit using
smaller devices corresponding to the fragments.

The QAOA ansatz with $p$ layers is a circuit whose structure is repeated $p$
times in succession. To aid our discussion on the relationship between problem
instances and the complexity of circuit cutting, it will be helpful to reduce
the scope of our analysis to the case of a single layer. The following
observation, which is not specific to the QAOA ansatz, motivates this choice.
For an $n$-qubit circuit, we define the support of a set of its gates $A$ to
be the subset of qubits upon which they act non-trivially,
$\text{supp}(A)\subseteq[n]$. Suppose this circuit is separated into fragments
$F_{1},\dots,F_{r}$ upon removing a subset of wires. Then we refer to the
support of the subset of gates belonging to the fragment $F_{j}$ as the
support of the fragment itself, denoted $\text{supp}(F_{j})\subseteq[n]$, for
each $j\in[r]$. This quantity represents the size of each of the sub-circuits
formed in the circuit cutting procedure.

######  Lemma D.1.

If a circuit can be separated into $r$ fragments $F_{1},\dots,F_{r}$ by
removing $\ell$ sets of at most $\kappa$ parallel wires (see Sec. II), then a
composition of $p$ layers of that circuit can be separated into $r$ fragments
$F_{1}^{\prime},\dots,F_{r}^{\prime}$ by removing at most $(2p-1)\ell$ sets of
at most $\kappa$ parallel wires. Furthermore,
$\textnormal{supp}(F_{j})=\textnormal{supp}(F_{j}^{\prime})$ for each
$j\in[r]$.

###### Proof.

First consider the (intact) single-layer circuit, where each gate is labelled
by the fragment to which it belongs when it is separated. That is, we assign a
label $j$ to each gate if it is in the fragment $F_{j}$ for some $j\in[r]$. We
may construct a new directed multi-graph by contracting all the edges (wires)
incident on gates with the same label. This is the communication graph in Ref.
[12], so we adopt this terminology here as well. The edges in the
communication graph represent those edges which separate the original circuit
into $r$ fragments when removed. We then choose an identical partition of the
gates in each of the $p$ repeated layers, resulting in a sequence of $p$
identical circuits whose gates are labelled using the labels $1,\dots,r$. The
partition of the gates in the $p$-layer circuit is then defined by these
labels, i.e., the fragment $F_{j}^{\prime}$ contains all gates labelled by
$j$, for each $j\in[r]$. It is clear that
$\textnormal{supp}(F_{j})=\textnormal{supp}(F_{j}^{\prime})$ for all
$j\in[r]$.

Let us now turn to the bound on the number and size of sets of parallel wires
in the $p$-layer communication graph. There is a contribution of $p\ell$ sets
of $\kappa$ parallel wires in the communication graph since each layer adds
$\ell$ sets of $\kappa$ inter-fragment wires. It remains to bound the number
of wires in the communication graph connecting gates in two different layers.
Suppose the layers are labelled $1,\dots,p$. For any $i\in[p-1]$, the wires
corresponding to the $a^{\text{th}}$ qubit connect gates with different labels
between layers $i$ and $i+1$ if only if the final gate acting upon the
$a^{\text{th}}$ qubit in the first layer has a different label from the first
gate acting upon the $a^{\text{th}}$ qubit. This happens only if there is a
wire in the single-layer communication graph which corresponds to the
$a^{\text{th}}$ qubit. Repeating this reasoning for each qubit leads to a
contribution to the total number of multi-graph edges that is bounded from
above by $\ell\kappa$ for each $i\in[p-1]$. Furthermore, the wires connecting
gates in layer $i$ to layer $i+1$ are parallel by definition, and so these can
be arbitrarily partitioned into at most $\ell$ sets of at most $\kappa$ wires
for each $i\in[p-1]$. In total, we have at most $p\ell+(p-1)\ell=2p-1$ sets of
at most $\kappa$ parallel wires which, upon removal, yield the desired
fragments $F_{1},\dots,F_{r}$. ∎

######  Lemma D.2.

Suppose there exists a partition of the edge set $E(G)$ into $r$ disjoint
subsets $E_{1},\dots,E_{r}$ and let $g_{1},\dots,g_{r}$ denote the subgraphs
of $G$ corresponding to these subsets of edges, respectively. Also, let
$1\leq\gamma\leq{r\choose 2}$ denote the number of pairs of indices
$i,j\in[r]$, $i<j$ such that $V(g_{i})\cap V(g_{j})\neq\emptyset$. Then there
exists a single-layer QAOA circuit for Max-Cut on $G$ that can be separated
into fragments $F_{1},\dots,F_{r}$ by removing $\gamma$ sets of at most
$\kappa$ parallel wires, where $\kappa:=\max_{i\neq j}|V(g_{i})\cap V(g_{j})|$
and $\mathrm{supp}(F_{j})=V(g_{j})\ \forall j\in[r]$.

###### Proof.

Note that in the single-layer QAOA circuit for Max-Cut, the vertex set $V$
represents the qubits in the circuit. The edge set $E$ represents the
commuting 2-qubit gates acting between qubits. Consider the circuit in which
all gates in $E_{1}$ are applied first, followed by $E_{2}$ and so on. Any
pair of wires connecting gates in $E_{i}$ to gates in $E_{j}$ for some $i<j$
are then parallel by construction. Furthermore, the set of qubits acted upon
by gates in both $E_{i}$ and $E_{j}$ is given by $V(g_{i})\cap V(g_{j})$, and
for each element in this set there is at most one wire connecting a gate in
$E_{i}$ to a gate in $E_{j}$ due to the ordering of gates we have chosen.
Hence, there are at most $\gamma$ nonempty sets of parallel wires which
connect gates in $E_{i}$ to $E_{j}$ for some $i<j$, each of which has at most
$\kappa$ elements. Removing these wires produces the fragments
$F_{1},\dots,F_{r}$ whose gates are given by $E_{1},\dots,E_{r}$,
respectively, such that $\text{supp}(F_{j})=V(g_{j})$ for every $j\in[r]$ as
claimed. ∎

We now prove the theorem using the above lemmas. For any $\delta\in[1/2,1]$ a
$\delta$-balanced vertex separator of a graph $G=(V,E)$ is a subset of
vertices $S\subseteq V$ whose removal leaves two sets of disconnected vertices
each of size at most $\delta|V|$.

######  Theorem III.1.

Suppose the graph $G=(V,E)$ with $|V|=n$ has some known $2/3$-balanced vertex
separator $S\subseteq V$ such that $|S|=\kappa$, and assume $\kappa\leq n/6$.
Then quantum computation with a $p$-layer QAOA circuit for Max-Cut on $G$ can
be simulated using a pair of (non-entangled) quantum devices each with at most
$5n/6$ qubits in time

$\displaystyle O\left(\frac{4^{(2p-1)(\kappa+2)}\textnormal{poly}(p\cdot
n)}{\varepsilon^{2}}\right).$ (62)

###### Proof.

We partition the edge set $E$ in the following way. Let $\widetilde{g}_{1}$,
$\widetilde{g}_{2}$ be the disconnected subgraphs which remain upon removing
the vertices $S$, and assume without loss of generality that

$\displaystyle|V(\widetilde{g}_{2})|\leq\lfloor
n/2\rfloor\leq|V(\widetilde{g}_{1})|\leq\lfloor 2n/3\rfloor.$ (63)

Let $g_{2}$ be the subgraph of $G$ corresponding to the vertices
$V(\widetilde{g}_{2})\cup S$. We then define a partition of $E$ into disjoint
sets $E_{1}$, $E_{2}\subseteq E$ by letting $E_{2}:=E(g_{2})$ and
$E_{1}:=E\backslash E_{1}$. Let $g_{1}$ be the subgraph of $G$ corresponding
to the edge set $E_{1}$. We have

$\displaystyle|V(g_{j})|\leq\lfloor 2n/3\rfloor+\kappa\leq\lfloor
5n/6\rfloor,\quad j=1,2$ (64) and $\displaystyle V(g_{1})\cap V(g_{2})=S.$
(65)

By Lemma D.2 a single layer of the QAOA circuit ansatz can be separated into
fragments $F_{1}$, $F_{2}$ by removing a single set of at most $\kappa$
parallel wires, and $|\text{supp}(F_{j})|=|V(g_{j})|\leq\lfloor 5n/6\rfloor$
is a bound on the number of qubits in the support of fragment $j$. Also, by
Lemma D.1, $p$ layers of this circuit ansatz can be separated111That the wires
to be removed in order to achieve this separation can be efficiently
determined follows from the explicit procedure for doing so outlined in the
proofs of Lemmas D.1 and D.2. into fragments $F_{1}^{\prime}$,
$F_{2}^{\prime}$ by removing at most $2p-1$ sets of $\kappa$ parallel wires,
with $|\text{supp}(F_{j}^{\prime})|\leq\lfloor 5n/6\rfloor$ for $j=1,2$. Then,
since the cost operator $\mathcal{C}$ for the Max-Cut problem is a diagonal
observable the expected cost $\langle\mathcal{C}\rangle$ can be written as

$\displaystyle\langle\mathcal{C}\rangle=\underset{z_{1},\dots,z_{2p-1}}{\mathbb{E}}\left[\langle\mathcal{C}^{\prime}_{z_{1},\dots,z_{2p-1}}\rangle\prod_{j=1}^{2p-1}(2^{\kappa+1}+1)(-1)^{z_{j}}\right]$
(66)

by Eq. (5). Here, $\langle\mathcal{C}^{\prime}_{z_{1},\dots,z_{2p-1}}\rangle$
is the expected value of a suitably modified circuit in which the channel
$\Psi_{z_{i}}$ is applied to the $i^{\text{th}}$ group of parallel wires for
each $i\in[2p-1]$, and $z_{1},\dots,z_{2p-1}$ are Bernoulli random variables
as described in Lemma II.1. The modified circuits can be executed on two
quantum devices with at most $\max_{j}|\text{supp}(F_{j}^{\prime})|\leq\lfloor
5n/6\rfloor$ qubits and classical communication. Using the random Clifford
strategy described in the proof of Theorem II.2, each of these circuits has at
most $\text{poly}(p\kappa|E|)$ gates. Measuring in the computational basis
allows one to construct an unbiased estimator of the expectation value
$\langle\mathcal{C}\rangle$ in a manner analogous to Line 7 in Algorithm 1.
This estimator is bounded in magnitude by $(2^{\kappa+1}+1)^{(2p-1)}$ with
certainty, so by Hoeffding’s Inequality
$(2^{\kappa+1}+1)^{2(2p-1)}/\varepsilon^{2}\leq
4^{(2p-1)(\kappa+2)}/\varepsilon^{2}$ repetitions of this procedure suffices
to estimate $\langle\mathcal{C}\rangle$ to within additive error
$\varepsilon$, with high probability. ∎

Figure 7: Schematic diagram of the idea behind Theorem II.2. The graph on the
left is the input to the Max-Cut problem, which has a balanced vertex
separator given by the vertices in yellow. The two-qubit gates are applied in
an order which respects the partition of the edges in the graph. We apply the
resolution of the identity for circuit cutting on the two wires indicated by
the dashed line.

##  Appendix E Circuit cutting complexity for QAOA with clustered graphs

We consider the execution of a multi-layered circuit resulting from encoding a
problem graph of the type introduced in Sec. III.1 for Max-Cut using QAOA.
There are two key components of the simulation to consider: the number of
unique sub-circuits $N$ required for circuit cutting, and the number of qubits
in each fragment circuit $m$. For a circuit of $p$ layers and a problem graph
of $r$ clusters, $n$ nodes within each cluster, and $k$ vertex separators, we
find:

$\displaystyle N=N(p,r,k)$
$\displaystyle=3^{pk}4^{(p-1)k}+(r-2)12^{(2p-1)k}+3^{(p-1)k}4^{pk},$ (67)

The number of qubits required to simulate any fragment is upper bounded by
$m=n+(3p-1)k$ (including additional auxiliary qubits due to mid-circuit
measurements). In Sec. III.4, by treating a GPU as a simulator analog of a
quantum hardware device, we aim for one GPU per fragment execution, i.e.,
$m=n+(3p-1)k\leq 30$.

##  Appendix F Analyzing the efficiency of simulated circuit cutting

This appendix considers how the simulation run time scales for performing
circuit cutting on QAOA circuits, for the problem graphs introduced in Sec.
III.1, as the number of classical resources is varied.

Increasing the number of GPUs used for the parallelized execution of the $N$
fragment sub-circuits demonstrates reasonable strong-scaling behaviour [60],
as it decreases the overall runtime for a QAOA circuit execution. However,
given that the number of circuit executions per GPU reduces, we reach
approximately 15-times speed-up using 32 GPU-nodes as depicted in Fig. 8,
yielding an efficiency of approximately 47% at our largest scale evaluation.
These inefficiencies are to be expected since increasing the number of GPU
nodes also increases the time needed for scheduling, transmission, device
setup overheads, and serial execution components which all contribute to the
overall runtime.

Figure 8:  Execution time and relative speed-up vs. the number of GPU nodes
for an execution of a 79-qubit QAOA circuit. The input problem graph
parameters are fixed at $p=1$, $r=3$, $n=25$, and $k=2$, while the number of
GPU nodes used in the distributed execution of the fragment circuits is
varied. The ideal linear scaling relationship is also provided with reference
to the single-node performance, and illustrates the divergence from ideal
strong-scaling behaviour at large node numbers.