Simulating Large Quantum Circuits on a Small Quantum Computer

TL;DR

This paper introduces a cluster-based simulation method enabling small quantum computers to simulate larger quantum systems or algorithms, with applications demonstrated on molecular energy estimation.

Contribution

The authors propose a novel cluster simulation scheme that decomposes quantum circuits into manageable clusters, allowing small quantum devices to simulate larger systems efficiently.

Findings

Simulation of large molecules using small quantum computers is feasible.

Cluster parameters $K$ and $d$ effectively characterize circuit complexity.

Experimental validation with BeH$_2$ molecule demonstrates practical applicability.

Abstract

Limited quantum memory is one of the most important constraints for near-term quantum devices. Understanding whether a small quantum computer can simulate a larger quantum system, or execute an algorithm requiring more qubits than available, is both of theoretical and practical importance. In this Letter, we introduce cluster parameters $K$ and $d$ of a quantum circuit. The tensor network of such a circuit can be decomposed into clusters of size at most $d$ with at most $K$ qubits of inter-cluster quantum communication. We propose a cluster simulation scheme that can simulate any $(K,d)$-clustered quantum circuit on a $d$-qubit machine in time roughly $2^{O(K)}$, with further speedups possible when taking more fine-grained circuit structure into account. We show how our scheme can be used to simulate clustered quantum systems -- such as large molecules -- that can be partitioned into multiple significantly smaller clusters with weak interactions among them. By using a suitable clustered ansatz, we also experimentally demonstrate that a quantum variational eigensolver can still achieve the desired performance for estimating the energy of the BeH$_2$ molecule while running on a physical quantum device with half the number of required qubits.

Tables

Notation	Definition
$C$	the quantum circuit
$f$	the classical post-processing function
$y$	the $n$-bit string obtained
	when measuring all qubits
$G$	the graph associated with
	the circuit (or the Hamiltonian)
$𝒜$	the set of tensors associated
	with the graph $G$
$K$	the number of edges
	between different clusters of $G$
$d$	the number of qubits
	sufficient for simulating each cluster
$g$	the graph obtained by contracting
	each cluster of $G$ to a single vertex
$\mathrm{cc}\left(g\right)$	the contraction complexity of $g$

Equations

$$T(G,\mathcal{A})={\mathbb{E}}_{y}f(y).$$

$$T(G,\mathcal{A})=\sum _{i=1}^8{c}_{i}T(G^{\prime },{\mathcal{A}}_i),$$

$$H=\sum _j{H}_j^{(1)}+\sum _j{H}_j^{(2)},\forall i,j:\Vert H_j^{(i)}\Vert \le 1,$$

$$U(\theta )=U_D({\theta }_D)U_{\text{ENT}}\cdots U_1({\theta }_1)U_{\text{ENT}}{U}_0({\theta }_0),$$

$$A(\rho {)}_{\alpha }$$

$${}_{\alpha }A(U{)}_{\beta }$$

$${}_{\beta }A(O)$$

$$\operatorname{Tr}\big{[}U\rho U^{\dagger}O\big{]}=\sum_{\alpha,\beta}A(\rho)_{\alpha}\cdot{}_{\alpha}A(U)_{\beta}\cdot{}_{\beta}A(O).$$

$$T(G,\mathcal{A})=\sum _{\alpha \in {\Pi }^E}\prod _{v\in V}A(v{)}_{\alpha (v)},$$

$$T(G,\mathcal{A})={\mathbb{E}}_{y}f(y),$$

$$T(G,\mathcal{A})=\sum _{i=1}^8{c}_{i}T(G^{\prime },{\mathcal{A}}_i).$$

$$I$$

$$A=\frac{\mathrm{Tr}(AI)I+\mathrm{Tr}(AX)X+\mathrm{Tr}(AY)Y+\mathrm{Tr}(AZ)Z}{2}.$$

$$O_1$$

$$O_2$$

$$O_3$$

$$O_4$$

$$O_5$$

$$O_6$$

$$O_7$$

$$O_8$$

$${\Phi }_i(A)=\mathrm{Tr}(A{O}_i){\rho }_i$$

$$A=\sum _{i=1}^8{c}_i{\Phi }_i(A).$$

$$\Phi =\sum _{i=1}^8{c}_i{\Phi }_i.$$

$${}_{\alpha }A(I{)}_{\beta }=\sum _{i=1}^8{c}_i\cdot {}_{\alpha }A(O_i)\cdot A({\rho }_i{)}_{\beta }.$$

$${}_{\alpha}A(I)_{\beta}=\operatorname{Tr}\big{[}M(\alpha)M(\beta)^{\dagger}\big{]}=\delta_{\alpha,\beta}.$$

$$\sum_{i=1}^{8}c_{i}\cdot\operatorname{Tr}\big{[}M(\alpha)O_{i}\big{]}\cdot\operatorname{Tr}\big{[}\rho_{i}M(\beta)^{\dagger}\big{]}.$$

$$\operatorname{Tr}\big{[}\big{\lparen}M(\alpha)\otimes M(\beta)^{\dagger}\big{\rparen}\cdot S\big{]}$$

$$S$$

$$=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right)$$

$$\displaystyle\operatorname{Tr}\big{[}\big{\lparen}M(\alpha)\otimes M(\beta)^{\dagger}\big{\rparen}\cdot S\big{]}$$

$$={\delta }_{\alpha ,\beta }$$

$$\mathrm{Pr}\left[\left|\bar{X}-\mathbb{E}[\bar{X}]\right|<\delta \right]\ge 1-2\exp \left(-\frac{N{\delta }^2}{2a^2}\right)>\frac{2}{3},$$

$$T(G,\mathcal{A})$$

$$c_s$$

$$T(G,\mathcal{A})=\mathbb{E}_{s}\big{[}8^{K}c_{s}T(G^{\prime},\mathcal{A}_{s})\big{]}.$$

Full text

*Introduction.—*Near-term quantum computing applications will focus on Noisy Intermediate-Scale Quantum (NISQ) devices Preskill (2018), where quantum memory is limited both in quantity and quality. To meet the memory need of such applications (e.g., quantum simulation Lloyd (1996); Cirac and Zoller (2012); O’Malley et al. (2016), quantum optimization Farhi et al. (2014); Moll et al. (2018), and quantum machine learning Rebentrost et al. (2014); Biamonte et al. (2017)), it is desirable to seek a way to perform computations that require more qubits than physically available, at the cost of additional affordable classical processing.

Trading classical computation for quantum computation is a well-motivated topic of long-standing interest. An extreme example of this is the (fully) classical simulation with no quantum computation at all, which is however limited to small dimensions, weak entanglement, or circuits with special gate sets Markov and Shi (2008); Vidal (2003); Chen et al. (2018a); Shi et al. (2006); Jozsa (2006). Recently, the possibility of trading classical computation for “virtual qubits” has been discussed in Bravyi et al. (2016). A systematic understanding of such trade-offs will be crucial for realizing near-term quantum applications.

In this Letter, we introduce a cluster simulation scheme, a general framework for simulating large quantum circuits on a quantum computer with a small amount of quantum memory. The performance of our simulation depends on the cluster parameters of the given circuits. In particular, we are inspired by the classical fragmentation methods and Quantum Mechanics/Molecular Mechanics (QM/MM) methods for simulating molecules Gordon et al. (2012); Li et al. (2007); Warshel and Levitt (1976); Li et al. (2008) that can be partitioned into multiple weakly interacting clusters of significantly smaller size (e.g., an oligosaccharide consisting of several monosaccharides). Following the spirit of Bravyi et al. (2016) and Gordon et al. (2012); Li et al. (2007); Warshel and Levitt (1976); Li et al. (2008), a natural definition of cluster parameters of a quantum circuit should capture the decomposability of the circuit into clusters of bounded size and limited inter-cluster interactions.

Our definition of the cluster parameters is guided by the above intuition, but with an important distinction. Instead of looking into the decomposability of any given circuit, we are concerned about the decomposability of the corresponding tensor network, which is inspired by the tensor-network-based classical simulation of quantum circuits Markov and Shi (2008); Arfken and Weber (2013); Orús (2014); Arad and Landau (2010). One significant advantage of our definition, as we will see below, is to use more flexible decompositions of tensor networks than are possible with simple partitioning of qubits (e.g., as in Bravyi et al. (2016)).

Given our definition of clustered circuits, our main contribution is a scheme to simulate the entire quantum circuit by simulating each cluster on a small quantum machine with classical post-processing. The key difference between our scheme and fully classical simulation schemes Bernstein and Vazirani (1997); Markov and Shi (2008); Vidal (2003); Chen et al. (2018a); Shi et al. (2006); Jozsa (2006); Markov et al. (2018); Boixo et al. (2017); Aaronson and Chen (2016); Chen et al. (2018b) is that we keep part of the computation quantum (i.e., unitary). In particular, we design a method to simulate inter-cluster quantum interactions by classical means. Comparing to “virtual qubits” in Bravyi et al. (2016), our technique can be deemed as trading classical computation for “virtual quantum communication”. The cluster simulation scheme applies to general quantum circuits, which distinguishes it from application-specific techniques for saving qubits Romero et al. (2017); Bravyi et al. (2017); Moll et al. (2016); Steudtner and Wehner (2018).

We apply our scheme to Hamiltonian simulation Vidal (2004); Berry et al. (2007, 2015); Low and Chuang (2017), particularly for clustered quantum systems, and Variational Quantum Eigensolvers (VQE), a popular candidate for showing quantum advantages on near-term quantum devices McClean et al. (2016); Barrett et al. (2013); Wecker et al. (2015). In both applications, we show that the particular quantum circuits have favorable cluster parameters that are amenable to our techniques. One can also interpret our technique as a hybrid variational ansatz in which a quantum computer is used for some parts of the circuit and a classical computer is used for other parts, which might be of independent interest.

Our scheme can easily be extended to allow limited inter-cluster quantum communication. It can hence be leveraged to perform general quantum circuits on modular quantum systems, a leading proposal of scalable quantum computers (e.g., Monroe et al. (2016); Jiang et al. (2007); Yao et al. (2012); Dai et al. (2020)).

Computational model.— We use the same computational model as in Bravyi et al. (2016); see Fig. 1 (Phase 1). An $m$-gate quantum circuit $C$ with 1- and 2-qubit gates is applied to $|0\rangle^{\otimes n}$ and all output qubits are measured in the computational basis. A classical post-processing function $f:\{0,1\}^{n}\to[-1,1]$ is then applied to the measurement outcomes. We assume that $f$ can be efficiently computed classically. We call the overall procedure a quantum-classical algorithm (QC algorithm) and denote it by $(C,f)$. Its expected output $\mathbb{E}_{y}f(y)$ is averaged over all measurement outcomes $y\in\{0,1\}^{n}$. The goal of our simulation is to approximate $\mathbb{E}_{y}f(y)$ within precision $\epsilon$ with high probability, say at least $2/3$.

*Clustered circuits.—*Any QC algorithm $(C,f)$ can be represented by a tensor network $(G,\mathcal{A})$ consisting of a directed graph $G(E,V)$ and a collection of tensors $\mathcal{A}=\{A(v):v\in V\}$. The vertices $V$ of $G$ represent individual gates (denoted by $\Box$), input qubits (denoted by $\lhd$), and observables (denoted by $\rhd$) as shown in Fig. 1 (Phase 2), whereas the flow of qubits is encoded by the directed edges $E$ of $G$. Note that each gate vertex $\Box$ has the same in- and out-degree (i.e., the same number of incoming and outgoing edges) whereas $\lhd$ vertices only have outgoing edges and $\rhd$ vertices only have incoming edges. For each $v\in V$, $A(v)$ is a tensor that encodes the matrix entries of the corresponding gate, state, or observable, and the value $T(G,\mathcal{A})$ of the tensor network $(G,\mathcal{A})$ coincides with the output expectation of the corresponding $(C,f)$ algorithm, i.e.,

[TABLE]

See the Supplemental Material for more details.

A QC algorithm $(C,f)$ is $(K,d)$-clustered if its tensor network $(G,\mathcal{A})$ has the following structure. Setting the final observable $O_{f}$ aside, we partition the remaining vertices of $G$ into clusters $S_{1},\dotsc,S_{r}$ and let $g$ be the $(r+1)$-vertex multi-graph obtained by contracting each cluster to a single vertex. Let $K$ be the number of edges in $g$ minus the in-degree of $O_{f}$ (intuitively, $K$ is the total number of qubits communicated between clusters) and let $d$ be the number of qubits sufficient for simulating each cluster.

While finding the minimal $d$ can be non-trivial (especially if qubits can be recycled after measurement), a good estimate of $d$ is $\max_{i}d(S_{i})$ where $d(S_{i})$ is the out-degree of cluster $S_{i}$. This is a valid upper bound on the minimal $d$ since $d(S_{i})$ is the number of $\lhd$ vertices in $S_{i}$ plus the number of incoming edges to $S_{i}$, which upper bounds the total number of qubits required to simulate $S_{i}$.

For example, in Fig. 1 (Phase 2), two parts of a partition $\{S_{1},S_{2}\}$ are indicated by blue and orange, respectively. Since only one qubit is sent from $S_{1}$ to $S_{2}$, $K=1$. We have $d(S_{1})=3$ due to two outgoing edges from $S_{1}$ to $O_{f}$ and one from $S_{1}$ to $S_{2}$. Similarly, $d(S_{2})=2$ and thus we can take $d=3$. This circuit is hence $(1,3)$-clustered.

Our framework generally allows for more flexibility in decomposing quantum circuits compared to Bravyi et al. (2016). Consider the $2n$-qubit example in Fig. 2. Assuming each block $B_{i}$ with depth $D$ is dense, i.e., contains two-qubit gates between all pairs of qubits, any partition of the initial $2n$ qubits induces at least $\Omega(n)$ gates across the parties and thus requires $\Omega(n)$ qubits of communication to implement the circuit. However, this is a $(2,n+1)$-clustered circuit with only two qubits of communication between the blue and orange clusters in Fig. 2. Furthermore, the depth is reduced from $3D$ to $2D$ when simulating each cluster separately.

*Cluster simulation scheme.—*To fulfill the above intuition, for any $(K,d)$-clustered circuit, we need to show how (i) each cluster can be simulated on a $d$-qubit quantum machine and (ii) how to simulate the interaction among clusters. We design an edge-cutting procedure to decompose tensor networks as shown in Fig. 1 (Phase 3). In particular, we replace each edge, modeled as a perfect channel for communicating a qubit, by a collection of tensor networks that reproduce the inter-cluster communication by operations within each cluster. This, however, comes at a cost of having to average over several runs. (In the spirit of Bravyi et al. (2016), this technique can be thought of as “virtual quantum communication.”)

[TABLE]

where $G^{\prime}$ differs from $G$ by removing $e$ and adding one $\lhd$ and one $\rhd$ vertex, each $c_{i}\in\{-\frac{1}{2},\frac{1}{2}\}$, and each $(G^{\prime},\mathcal{A}_{i})$ corresponds to a valid quantum circuit.

(All proofs in this Letter are deferred to the Supplemental Material 111See Supplemental Material for more details about (1) the precise definition of tensor network corresponding to clustered circuits; (2) the proofs of Lemma 1, Theorem 1, 2, 3; (3) the details of experiments about VQE, which includes Refs. Hoeffding (1963); Miller (1972); Mitarai and Fujii (2019); Eisert et al. (2010)..) By repeating this process and deleting more edges, the tensor network can eventually be partitioned into individual clusters. Each cluster will only have outgoing edges to $O_{f}$ and can hence be simulated by a $d$-qubit quantum computer plus classical processing of the measurement outcomes (Phase 4 in Fig. 1). We combine individual simulation results by a simple sampling procedure according to Eq. (2).

Our overall simulation scheme consists of several iterations of the following steps: (i) producing a classical description of a quantum circuit with $O(m)$ gates and $d$ qubits (potentially recycled during the circuit), (ii) running this circuit on $|0\rangle^{\otimes d}$, and (iii) classically post-processing the measurement outputs. The final step has to produce with probability at least 2/3 an $\epsilon$-approximation of $T(G,\mathcal{A})$.

The complexity of our scheme scales with the cluster parameters $(K,d)$ as well as the total number of qubits $n$ and gates $m$ in the original circuit $C$, and the desired additive simulation accuracy $\epsilon$. The total classical and quantum running time of our simulator is $O(Q\operatorname{poly}(n+m)/\epsilon^{2})$, for some exponentially scaling parameter $Q$. For simplicity, we ignore the polynomial part of the run-time and call this a $(Q,d)$-simulator. A fully classical simulator is thus a $(2^{O(n)},0)$-simulator, while a scalable quantum computer is a $(1,n)$-simulator with an exponentially improved total run-time. Our result can be deemed as a smooth trade-off between these two extreme cases.

In the special case when there are only two clusters, the number of qubits $K$ communicated among the clusters can be regarded as an upper bound on entanglement. Hence, the result relates the classical computation cost to the entanglement between the two clusters.

The efficiency of our simulation can be further improved for special classes of post-processing functions $f:\{0,1\}^{n}\to[-1,1]$. For example, consider decomposable $f$ satisfying $f(y)=\prod_{j=1}^{r}f_{j}(y_{j})$, where $y=y_{1}\dots y_{r}$ is a partition of the original $n$-bit string $y$ into substrings $y_{j}$ that correspond to outputs of different clusters 222We have assumed that the number of terms in the decomposition of $f$ agrees with the number of clusters $r$. If some cluster does not produce a qubit that feeds directly into the final observable, we can insert a fictitious function $f_{j}$ in the decomposition ($f_{j}$ has no arguments and is identically equal to $1$)., and $f_{j}(y_{j})\in[-1,1]$. Typical examples of such decomposable functions arise from Pauli observables in VQE Kandala et al. (2017) or estimating probabilities of specific output strings Pashayan et al. (2017). For such functions, we can replace $O_{f}$ by smaller tensors $O_{f_{j}}$ and include them in the corresponding clusters $S_{j}$. As a result, the induced graph $g$ no longer contains $O_{f}$. Nevertheless, we can still apply Lemma 1 to decompose each cluster and simulate it on a $d$-qubit quantum machine. However, inspired by Markov and Shi (2008), a more efficient scheme for combining individual simulations is now possible. Its complexity depends on $\operatorname{cc}(g)$—the contraction complexity of $g$—that is the minimum (over all possible contraction orders) of the maximum node degree during the procedure of contracting the graph to a single vertex.

Note that $\operatorname{cc}(g)\leq K$, where $K$ is the number of edges in $g$. Compared to $2^{O(K)}$ in Theorem 1, the factor $2^{O(\operatorname{cc}(g))}$ in Theorem 2 is a significant improvement for some families of graphs. For example, among constant-degree graphs with $n$ nodes, $\operatorname{cc}(g)=O(1)$ for trees and $\operatorname{cc}(g)=O(\sqrt{n})$ for planar graphs Bodlaender (1994), while $K$ can be as large as $O(n)$.

*Application to Hamiltonian simulation.—*One of the most promising potential applications of our result is the simulation of clustered quantum systems. Specifically, we consider quantum systems with geometric layouts where each qubit only interacts with $O(1)$ adjacent qubits. The corresponding interaction graph $G$ (i.e., qubits as vertices and interactions as edges) has constant degree. Assume further that qubits in $G$ can be grouped into $n$ parties and let $g$ be the induced graph obtained by contracting each party of $G$ to a single vertex. The Hamiltonian of such a system can be written as a sum of local terms

[TABLE]

where each term acts on at most two qubits and the superscripts $(1)$ and $(2)$ indicate that these qubits belong to a single party or two different parties, respectively. The interaction strength between all parties can be characterized by $h=\sum_{j}\lVert H_{j}^{(2)}\rVert$. We are interested in quantum systems with weak interaction strength (e.g., Fig. 3). Assume that the system is initialized in a product state $\rho=\rho_{1}\otimes\dotsb\otimes\rho_{n}$, where $\rho_{i}$ is an efficiently preparable state of the $i$-th party. Our goal is to approximate the following correlation function: \operatorname{Tr}\big{[}(O_{1}\otimes\dotsb\otimes O_{n})e^{-iHt}(\rho_{1}\otimes\dotsb\otimes\rho_{n})e^{iHt}\big{]}, where $t$ is the evolution time and $O_{i}$ is an efficiently measurable observable of the $i$-th part with eigenvalues in $[-1,1]$.

At a high level, the above result is obtained by applying Theorem 2 to Hamiltonian simulation circuits of $e^{-iHt}$ based on the Trotter-Suzuki approximation, but with the following important improvements. To obtain a better estimate of the cluster parameters $(K,d)$, we need to apply Lemma 1 to trim the tensor network beyond simulating inter-cluster communication, and to conduct a careful analysis of $d$ to allow recycling of qubits. Inspired by Haah et al. (2018), we also need to improve the naive error analysis and to obtain an error bound in terms of the interaction strength $h$.

The exponential dependence on $t$ 333This dependence on $t$ has been subsequently improved by using a $p$-th order product formula ($p>1$). See Childs et al. (2019) for details. seems necessary as suggested by hardness results of classical simulation of quantum circuits (e.g., Terhal and DiVincenzo (2004)). It was also previously known that a classical algorithm can estimate local observables in time exponential in the size of the light-cone Hastings (2004); Osborne (2006), i.e., the number of input qubits that could influence a particular output qubit, resulting in a similar run-time bound. (For Hamiltonian evolution we still have an effective light-cone due to Lieb-Robinson bounds, e.g., Haah et al. (2018).) Our approach is, however, strictly stronger in the sense that we could estimate correlations across the entire system, something that cannot be achieved by the light-cone argument.

*Application to VQE.—*Variational Quantum Eigensolver (VQE) is a variational method for finding the lowest eigenvalue of an $n$-qubit Hamiltonian $H$ by applying some parameterized circuit $U(\theta)$ to $|0\rangle^{\otimes n}$ and minimizing the expectation with $H$: $\min_{\theta}\;\langle 0|^{\otimes n}U(\theta)^{{\dagger}}HU(\theta)|0\rangle^{\otimes n}.$

This method has been proposed for solving optimization problems on quantum computers Wecker et al. (2015); McClean et al. (2016); Farhi et al. (2014) and, thanks to its short-depth circuits, has become a promising candidate to surpass the classical optimization methods and show quantum advantage on NISQ devices Preskill (2018); Peruzzo et al. (2014); O’Malley et al. (2016); Kandala et al. (2017); Moll et al. (2018).

In Kandala et al. (2017), Kandala et al. propose a class of hardware-friendly variational circuits $U(\theta)$ and experimentally demonstrate the effectiveness of their VQE method for addressing problems of small molecules and quantum magnetism, using up to $6$ qubits. Their ansatz $U(\theta)$ has the following form:

[TABLE]

where $U_{i}(\theta_{i})=\bigotimes_{j=1}^{n}U_{i}^{j}(\theta_{i}^{j})$ and each $U_{i}^{j}(\theta_{i}^{j})$ is a parameterized single-qubit gate applied on the $j$-th qubit out of $n$, $U_{\text{ENT}}$ is a fixed sequence of two-qubit gates meant for producing entanglement, and $D$ is the number of rounds.

In the context of current NISQ devices, we propose a way to reduce the number of qubits required for implementing $U(\theta)$ by using our cluster simulation scheme. This involves the following steps: (i) choosing a partition $\mathcal{P}=\{S_{1},\dotsc,S_{r}\}$ of $n$ qubits such that $\lvert S_{i}\rvert\leq d$ for each $i$; (ii) removing some entangling gates from $U_{\text{ENT}}$ that go across different parts of $\mathcal{P}$ to decrease $\operatorname{cc}(g)$, where $g$ is the graph induced by regarding each set $S_{i}$ as a node and each gate that acts across two sets as an edge; (iii) runing this $n$-qubit $U(\theta)$ using a $(2^{O(\operatorname{cc}(g))},d)$-simulator.

We report an experiment estimating the ground energy of the $\text{BeH}_{2}$ molecule; see Fig. 7. Using a 3-qubit physical device, we run the 6-qubit $U(\theta)$ from Kandala et al. (2017) and achieve the similar accuracy, thus demonstrating the potential of implementing VQE with limited quantum memory. Additional details about the experiment and the discussion of reducing $\operatorname{cc}(g)$ can be found in the Supplemental Material.

*Summary.—*In this Letter, we provide a systematic approach for simulating clustered quantum circuits with limited use of quantum memory. Our scheme is relevant to promising NISQ applications such as Hamiltonian simulation and VQE. By reducing the number of qubits and the depth of the circuit, it is particularly applicable to intermediate scale devices and potentially also improves the circuit’s robustness to correlated noise. We leave open the problem of determining the best $(K,d)$ or $(\operatorname{cc}(g),d)$ for a given quantum circuit (this may be related to the graph partitioning, graph clustering, and treewidth problems). Another direction is to develop more case-by-case optimization techniques for realistic applications under our scheme.

I Supplemental Material

I.1 I: Table of Notation

I.2 II: Tensor networks

Any quantum circuit $C$ can be represented by a directed graph $G=(V,E)$ with three types of vertices: input states (denoted by $\lhd$), quantum gates (denoted by $\Box$), and observables (denoted by $\rhd$), with directed edges indicating the flow of qubits between them. We can further associate to each vertex $v\in V$ a tensor $A(v)$ that encodes the matrix entries of the object represented by this vertex, and denote by $\mathcal{A}=\{A(v):v\in V\}$ the collection of all these tensors. This results in a tensor network $(G,\mathcal{A})$ that completely captures the algorithm described by the original quantum circuit $C$ Vidal (2003); Shi et al. (2006); Jozsa (2006); Markov and Shi (2008). In particular, any QC algorithm $(C,f)$ can also be represented by such a tensor network when we combine the final computational basis measurement and the classical post-processing function $f$ into a diagonal observable $O_{f}$ (see Figure 1, Phase 2 in the main paper).

Our index convention for $k$-qubit tensors is borrowed from Markov and Shi (2008). Consider $\alpha=(\alpha_{1},\dotsc,\alpha_{k})$ where each $\alpha_{j}=(\alpha_{j}^{1},\alpha_{j}^{2})$ is a pair of binary indices, i.e., $\alpha_{j}\in\Pi=\{0,1\}^{2}$ 444We need two indices per qubit since we are working with an operator rather than a vector basis.. Let $M(\alpha_{j})=|\alpha_{j}^{1}\rangle\langle\alpha_{j}^{2}|$ be the elementary matrix corresponding to $\alpha_{j}$ and let $M(\alpha)=\bigotimes_{j=1}^{k}M(\alpha_{j})=|\alpha_{1}^{1},\dotsc,\alpha_{k}^{1}\rangle\langle\alpha_{1}^{2},\dotsc,\alpha_{k}^{2}|$ be its $k$-qubit generalization. For a state $\rho$, gate $U$, and observable $O$, the entries of their tensors are computed as follows:

[TABLE]

where the location of subscripts indicates whether qubits are incoming or outgoing. For example,

[TABLE]

This corresponds to the following tensor network:

\rho$$U$$O$$\alpha$$\beta

where arrows indicate the flow of information.

For a general tensor network $(G,\mathcal{A})$, its value is

[TABLE]

where $\alpha=(\alpha_{e}:e\in E)$ and each $\alpha_{e}=(\alpha_{e}^{1},\alpha_{e}^{2})\in\Pi$ ranges over all possible index assignments 555Equivalently, one can think of each $\alpha\in\Pi^{E}$ as a function of the form $\alpha:E\to\Pi$. for edge $e$, and $\alpha(v)$ denotes the restriction of $\alpha$ to edges incident to $v$. In particular, $A(v)_{\alpha(v)}$ is the entry of $A(v)$ when all edges incident to $v$ are labeled according to $\alpha$ 666Our notation in Eq. 9 ignores the distinction between incoming and outgoing edges, but one has to take this into account when evaluating the tensor entries using Eqs. 5, 6 and 7.. Note that

[TABLE]

so a quantum circuit can be simulated by approximating the value of its tensor network Markov and Shi (2008).

I.3 III: Proof of Lemma 1

Lemma 1**.**

Let $(G(E,V),\mathcal{A})$ be a tensor network of a QC algorithm and let $uv$ be any edge in $G$. There is a set of eight coefficients $c_{i}\in\{-\frac{1}{2},\frac{1}{2}\}$, observables $O_{i}$, and states $\rho_{i}$, such that the following modification of the edge $uv$

u$$v$$\Longrightarrow$$\displaystyle\sum_{i=1}^{8}c_{i}$$u$$O_{i}$$\rho_{i}$$v

does not affect the value of the overall tensor network. Moreover, each term above yields some tensor network $(G^{\prime},\mathcal{A}_{i})$ that again corresponds to a valid quantum circuit, and their values are related as follows:

[TABLE]

Proof.

Let $A$ be any $2\times 2$ matrix and let

[TABLE]

be the Pauli matrices. Since the normalized Pauli matrices $\{I,X,Y,Z\}/\sqrt{2}$ form an orthonormal matrix basis,

[TABLE]

Let us expand each Pauli matrix in its eigenbasis. If we denote the Pauli marices by $O_{i}$, their eigenprojectors by $\rho_{i}$ and the corresponding eigenvalues by $2c_{i}$ so that $O_{1}=O_{2}=2(c_{1}\rho_{1}+c_{2}\rho_{2})$, etc., then

[TABLE]

where $|\pm\rangle=\frac{|0\rangle\pm|1\rangle}{\sqrt{2}}$ and $|{\pm}i\rangle=\frac{|0\rangle\pm i|1\rangle}{\sqrt{2}}$.

For each case, define a corresponding map

[TABLE]

that can be implemented by first measuring the Pauli observable $O_{i}$ and then preparing the corresponding eigenstate $\rho_{i}$. It follows immediately that

[TABLE]

If $\Phi$ denotes the trivial single-qubit channel that acts as $\Phi(A)=A$, for any $2\times 2$ matrix $A$, then Eq. 21 is equivalent to

[TABLE]

The same decomposition works even when $\Phi$ acts on a subsystem of a larger system.

To prove Lemma 1, note that adding the trivial map $\Phi$ on any edge of $G$ does not affect the value of the tensor network $(G,\mathcal{A})$. Thus our strategy for constructing the decomposition $\{(G^{\prime},\mathcal{A}_{i})\}_{i=1}^{8}$ is to explicitly include the trivial map $\Phi$ on the edge $uv$, and then decompose $\Phi$ according to Eq. 22. Intuitively, we are cutting the edge $uv$ and mending its both ends with additional tensors in several different ways. This gives us a collection of eight new tensor networks $(G^{\prime},\mathcal{A}_{i})$ that each correspond again to a quantum circuit: the $i$-th circuit measures the qubit produced at vertex $u$ with observable $O_{i}$ and supplies the qubit required at vertex $v$ in state $\rho_{i}$.

To guarantee that the modified tensor network has the same value as the original, we have to verify that

\alpha$$\beta$$I$$=$$\displaystyle\sum_{i=1}^{8}c_{i}$$\alpha$$O_{i}$$\rho_{i}$$\beta

where on the left-hand side we use conjugation by the identity matrix $I$ as a way of implementing the trivial channel $\Phi$. Following the notation of Eqs. 5, 6 and 7, this is equivalent to verifying

[TABLE]

According to Eq. 7, the left-hand side is

[TABLE]

From Eqs. 5 and 7, the right-hand side is

[TABLE]

We can also write this as

[TABLE]

where $S=\sum_{i=1}^{8}c_{i}(O_{i}\otimes\rho_{i})$. Note from Eq. 12 to Eq. 19 that

[TABLE]

is the two-qubit SWAP operator, so

[TABLE]

which proves Eq. 23. ∎

I.4 Hoeffding’s inequality

Our proofs of Theorems 1 and 2 rely on the convergence speed of averaged i.i.d. random variables. They use the following inequality from Hoeffding (1963), which we state here for the sake of reference.

Hoeffding’s inequality**.**

Suppose $X_{1},\dotsc,X_{N}$ are i.i.d. real random variables such that $|X_{i}|\leq a$ and let $\bar{X}=(X_{1}+\dotsb+X_{N})/N$. Then

[TABLE]

where the last inequality holds if $N\geq 4a^{2}/\delta^{2}$.

I.5 IV: Proof of Theorem 1

Theorem 1**.**

Any QC algorithm $(C,f)$ with a $(K,d)$-clustered circuit $C$ consisting of $r$ clusters can be simulated with precision $\epsilon$ by running a $d$-qubit quantum device $O\left(2^{4K}r/\epsilon^{2}\right)$ times. The total classical and quantum run-time of this simulation is $O\left(2^{4K}(n+m)/\epsilon^{2}\right)$ where $n$ and $m$ are the number of qubits and gates in $C$, respectively.

Proof.

From the discussion leading up to Eq. 10 we know that simulating a quantum circuit is equivalent to approximating $T(G,\mathcal{A})$, the value of its tensor network given in Eq. 9. Let ${E_{\text{cut}}}$ denote the set of edges that connect different clusters of $G$. Since $G$ is $(K,d)$-clustered, $\lvert{E_{\text{cut}}}\rvert=K$. Let us decompose each edge $e\in{E_{\text{cut}}}$ according to Lemma 1 and denote the modified graph by $G^{\prime}$. If we collect all indices appearing in this decomposition into a single string $s=(s_{e}:e\in{E_{\text{cut}}})$, where $s_{e}\in[8]=\{1,\dotsc,8\}$ represents the index $i$ of edge $e\in{E_{\text{cut}}}$, then

[TABLE]

where $c^{e}_{1},\dotsc,c^{e}_{8}\in\{-\frac{1}{2},\frac{1}{2}\}$ are the coefficients $c_{i}$ from Lemma 1. Regarding $s$ as a random variable sampled uniformly from $[8]^{{E_{\text{cut}}}}$,

[TABLE]

To estimate the value of $T(G^{\prime},\mathcal{A}_{s})$, for a given $s\in[8]^{{E_{\text{cut}}}}$, recall from Lemma 1 that each $(G^{\prime},\mathcal{A}_{s})$ is also a quantum circuit. Moreover, since $G$ is $(K,d)$-clustered, each $(G^{\prime},\mathcal{A}_{s})$ consists of $r$ pieces (not counting the final observable $O_{f}$), each of which can be executed on a $d$-qubit quantum computer. Our simulator runs these pieces sequentially. To obtain the correct result, we need to initialize and measure all qubits of each piece in the right basis. Our strategy is as follows. We initialize an input qubit of a given piece to $|0\rangle$ if it corresponds to an input qubit also in the original circuit. Similarly, we measure an output qubit in the standard basis if it is originally fed into the final observable $O_{f}$. The remaining input and output qubits are supposed to be communicated between different clusters. In the simulation, however, we respectively initialize them to Pauli eigenstates and measure using Pauli observables determined by $s$. Specifically, if $e\in{E_{\text{cut}}}$ is an edge between two clusters, its endpoints correspond to a pair of input / output qubits that are initialized / measured according to the value of $s_{e}\in[8]$ (see the proof of Lemma 1).

After executing all pieces as described above, we obtain a collection of Pauli measurement outcomes $\sigma^{e}_{s_{e}}\in\{-1,1\}$, one for each $e\in{E_{\text{cut}}}$, and a collection of strings $y_{1},\dots,y_{r}$ of total length $n$ 777Some of these strings might be empty. This happens when all outgoing qubits from some cluster are fed into other clusters and none are fed directly into the final measurement.. We can now compute the following quantity by evaluating the classical post-processing function $f$:

[TABLE]

Note that, for each $s\in[8]^{{E_{\text{cut}}}}$, $t_{s}\in[-1,1]$ is a random variable with expectation $T(G^{\prime},\mathcal{A}_{s})$. Hence, using Eq. 33 we can express the value of $T(G,\mathcal{A})$ as an expectation over a combined random variable $(s,t_{s})$, where $s$ is drawn uniformly from $[8]^{{E_{\text{cut}}}}$ and $t_{s}$ is obtained by running the corresponding quantum circuit as described above:

[TABLE]

Our simulation approximates $T(G,\mathcal{A})$ by producing $N$ samples $\{(s_{1},t_{s_{1}}),\dotsc,(s_{N},t_{s_{N}})\}$ of the random variable $(s,t_{s})$ and computing their average $\frac{1}{N}\sum_{i=1}^{N}8^{K}c_{s_{i}}t_{s_{i}}$. The magnitude of each term is bounded by

[TABLE]

where we used $|c_{s}|=\prod_{e\in{E_{\text{cut}}}}|c_{s_{e}}^{e}|=1/2^{K}$. By Hoeffding’s inequality (see above), we can achieve accuracy $\epsilon$ with success probability at least $\frac{2}{3}$ if the number of experiments $N$ satisfies

[TABLE]

Since each sample of $t_{s}$ requires running $r$ independent quantum circuits, where $r$ is the number of clusters, the total number of circuits run on the $d$-qubit quantum device is

[TABLE]

To bound the total running time, let us denote the number of qubits and gates in the $i$-th cluster by $n_{i}$ and $m_{i}$, respectively. The total processing time for simulating the $i$-th cluster once is $O(n_{i}+m_{i})$. Note that $\sum_{i=1}^{r}n_{i}+m_{i}=O(n+m)$, where $n$ is the total number of qubits and $m$ is the total number of gates of $C$. Hence the total complexity is

[TABLE]

as claimed. ∎

I.6 V: Proof of Theorem 2

Theorem 2**.**

An $n$-qubit and m-gate QC algorithm $(C,f)$ with a decomposable function $f$ and a $(K,d)$-clustered circuit $C$ can be simulated with precision $\epsilon$ by running a $d$-qubit quantum device O\big{\lparen}2^{O(d^{\prime}(g))}(r^{3}\log r)/\epsilon^{2}\big{\rparen} times. Here $g$ denotes the induced (multi-)graph obtained from $G$ by removing the vertex $O_{f}$ and contracting each cluster to a single vertex while keeping the edges between different clusters, $d^{\prime}(g)$ is the maximum degree of a vertex in $g$, and $r$ is the number of vertices in $g$ or clusters in $G$ (note that $r\leq n+m$). The total running time of this simulation is $T_{Q}+T_{C}$ where

[TABLE]

are the total quantum and classical running times, and $\operatorname{cc}(g)$ is the contraction complexity of $g$.

Proof.

Since $f$ is decomposable, $f(y)=\prod_{j=1}^{r}f_{j}(y_{j})$, so we can replace the final observable $O_{f}$ by $r$ independent smaller observables $O_{f_{1}},\dotsc,O_{f_{r}}$. Moreover, since each $O_{f_{j}}$ is connected to a different cluster, we can absorb it into the corresponding cluster. Let us denote this slightly modified tensor network by $(G,\mathcal{A})$ and note that it is also $(K,d)$-clustered just like the original one. If we cut each edge in ${E_{\text{cut}}}$ using Lemma 1 (just like we did in the proof of Theorem 1), the resulting decomposition is again given by Eqs. 31 and 32. However, this time the graph $G^{\prime}$ of the modified network fully breaks up into individual clusters 888Recall that ${E_{\text{cut}}}$ does not include the edges attached to $O_{f}$, so the original tensor network did not fully break up because the clusters connected to $O_{f}$ remained in the same connected component even when the edges in ${E_{\text{cut}}}$ were removed.:

S_{2}$$S_{3}$$S_{1}$$T(G,\mathcal{A})\;=\;\displaystyle\sum_{s\in[8]^{{E_{\text{cut}}}}}c_{s}$$T(G^{\prime},\mathcal{A}_{s})

If we denote these clusters by $S_{1},\dotsc,S_{r}$, we can further decompose each term as a product over them:

[TABLE]

where $G^{\prime j}$ is the $j$-th cluster (or the $j$-th component) of $G^{\prime}$ and $s(S_{j})$ denotes the restriction of $s$ to the edges incident to cluster $S_{j}$ (just like in Eq. 9). In particular, note that $\bigcup_{j=1}^{r}\mathcal{A}^{j}_{s(S_{j})}=\mathcal{A}_{s}$ and that each \big{\lparen}G^{\prime j},\mathcal{A}^{j}_{s(S_{j})}\big{\rparen} is a piece of the larger network $(G^{\prime},\mathcal{A}_{s})$ as well as a tensor network on its own.

Operationally, the value T\big{\lparen}G^{\prime j},\mathcal{A}^{j}_{s(S_{j})}\big{\rparen} in Eq. 42 is the expectation of the circuit of the $j$-th component of $(G^{\prime},\mathcal{A}_{s})$ whose all input qubits are initialized to specific eigenvectors of Paulis (determined by $s$) and all output qubits are measured by specific Pauli observables (also determined by $s$) or the observables $O_{f_{j}}$ whose eigenvalues lay in $[-1,1]$. As a consequence,

[TABLE]

We define a new (multi)-graph $g$ obtained by contracting each $S_{j}$ of graph $G$ to a single node. The vertices and edges of the graph $g$ can be identified with $\{S_{1},\dotsc,S_{r}\}$ and ${E_{\text{cut}}}$, respectively. For simplicity, let us abbreviate the vertex $S_{j}$ as $j$. Then, we define a new tensor network $(g,\mathscr{a})$ associated with the graph $g$. For each vertex $j$ of $g$, we assign a tensor $a(j)\in\mathscr{a}$ whose indices range over $[8]=\{1,\dotsc,8\}$ and whose entries are given by

[TABLE]

where $s(j)$ is the restriction of $s\in[8]^{{E_{\text{cut}}}}$ to $S_{j}$ and $c_{s(j)}$ is the product of $c^{e}_{s_{e}}$ (see Eq. 32), with $e$ ranging over all edges in ${E_{\text{cut}}}$ connecting $S_{j}$ and any $S_{j^{\prime}}$ with $j^{\prime}>j$. Note that the indices of edges in $(g,\mathscr{a})$ range over $[8]=\{1,\dotsc,8\}$, while the indices of edges in $(G,\mathcal{A})$ range over $\Pi=\{0,1\}^{2}$.

Intuitively, $a(j)_{s(j)}$ is obtained by contracting the $j$-th cluster of $\mathcal{A}_{s}$ and absorbing the coefficients $c^{e}_{s_{e}}$ into the cluster with the smaller index (each edge $e\in{E_{\text{cut}}}$ connects two different clusters and we always absorb into the one with the smallest index). Following Eq. 9, the value of $(g,\mathscr{a})$ is

[TABLE]

where we used equations (44), (42), and (31), and the fact that each $c^{e}_{s_{e}}$ occurs exactly once (from the cluster with the smallest index out of the two connected by the edge $e$). Because of this identity, we can instead focus on computing the value of the smaller tensor network $(g,\mathscr{a})$.

Our method of computing $T(g,\mathscr{a})$ starts by estimating all entries $a(j)_{s(j)}$ of tensors in $\mathscr{a}=\{a(1),\dotsc,a(r)\}$, which can be done by running a $d$-qubit quantum computer multiple times to estimate $T(G^{\prime j},\mathcal{A}^{j}_{s(j)})$. This yields an approximation $\tilde{\mathscr{a}}$ of $\mathscr{a}$. Then we contract $(g,\tilde{\mathscr{a}})$ using a classical algorithm. To obtain $T(g,\mathscr{a})$ with precision $\epsilon$, we need to compute the entries of each $a(j)$ to precision roughly $2^{-O(K)}$ where $K$ is the number of edges of $g$. We can improve this trivial estimate by using the fact that these tensors come from quantum circuits.

Recall from Eq. 43 that \big{\lvert}T(G^{\prime j},\mathcal{A}^{j}_{s(j)})\big{\rvert}\leq 1, for any $j\in\{1,\dotsc,r\}$ and $s\in[8]^{{E_{\text{cut}}}}$. Hence,

[TABLE]

from the definition in Eq. 44. More generally, consider a new tensor $a(L)$ obtained by contracting any subset $L\subseteq\{1,\dotsc,r\}=V(g)$ of vertices of $g$ and their associated tensors $\{a(v):v\in L\}$ to a single vertex. Just like we argued above Eq. 43, the subset $L$ is also associated with the expectation of some quantum circuit with bounded observables. Therefore, Eq. 49 extends to any subset of vertices:

[TABLE]

Using this, we can complete the proof by applying Lemma 2 (below)—it bounds by how much a tensor network’s value can deviate when each tensor is perturbed entry-wise by some small amount. Indeed, this is precisely what we need to bound $|T(g,\tilde{\mathscr{a}})-T(g,\mathscr{a})|$ where $\tilde{\mathscr{a}}$ is our estimate of $\mathscr{a}$.

Let us flesh out this strategy in more detail. Note that the total number of entries in all tensors is upper bounded by $D=r8^{d^{\prime}}$, where $r$ is the number of tensors and $d^{\prime}$ is the maximum degree of the graph $g$. To estimate a single entry of a given tensor, let us use the following parameters:

[TABLE]

where $\delta$ is the desired accuracy of our estimate, $N$ is the number of times we run the corresponding quantum circuit, and $\epsilon>0$ can be chosen arbitrarily (it determines the accuracy of our estimate of the tensor network’s value).

By Hoeffding’s inequality (see above), the probability that the estimate is not within $\delta$ of the actual value is at most

[TABLE]

By union bound, the probability that at least one out of $D$ estimates if off by more than $\delta$ is at most

[TABLE]

In other words, all estimates are correct with probability at least $2/3$, yielding a tensor network $(g,\tilde{\mathscr{a}})$ that is entry-wise $\delta$-close to the original network $(g,\mathscr{a})$. By Lemma 2 (see below), the values of these tensor networks differ by

[TABLE]

where $\epsilon$ can be chosen arbitrarily in $(0,1]$.

To analyze the total cost, note that our simulator executes at most $DN$ quantum circuits, where $D$ is an upper bound on the total number of entries of all tensors in the network and $N$ is the number of circuits run to estimate each entry. Expressing this in terms of the relevant parameters,

[TABLE]

where $r=\lvert V(g)\rvert$ and $d^{\prime}$ is the maximum degree of a vertex in $g$.

To compute the total quantum running time, note that $O(n+m)$ quantum operations are needed for running all $r$ clusters once (by a similar argument as in Theorem 1). Hence, we can replace a factor of $r$ by $O(n+m)$ and the total quantum running time is bounded by

[TABLE]

After all entries are estimated, we can classically compute the value of the tensor network $(g,\tilde{\mathscr{a}})$ in classical time $T_{C}=\operatorname{poly}(r)2^{O(\operatorname{cc}(g))}$ by finding a near-optimal tree decomposition of the line graph of $g$ Markov and Shi (2008). Noting that $d^{\prime}(g)\leq\operatorname{cc}(g)$, i.e., the maximum degree is at most the contraction complexity for any graph $g$, leads to the simplified formula in the main part of the paper. ∎

Lemma 2**.**

Let $(g,\mathscr{a})$ be a tensor network with $g=(V,E)$ and $|V|=r$. Assume that, for any subset $L\subseteq V$, all entries of the tensor $a(L)$ obtained by contracting $L$ to a single vertex have norm at most one. Let $\tilde{\mathscr{a}}$ be another tensor network that is entry-wise close to $\mathscr{a}$, i.e., for all $v\in V$ and $s(v)\in[8]^{d^{\prime}(v)}$,

[TABLE]

for some $0\leq\delta\leq 1/(r8^{d^{\prime}})$, where $d^{\prime}(v)$ denotes the degree of a vertex $v$ and $d^{\prime}=\max_{v}d^{\prime}(v)$ denotes the maximum degree of $g$. Then

[TABLE]

Proof.

We denote the entry-wise perturbation by $\Delta(v)_{s(v)}=\tilde{a}(v)_{s(v)}-a(v)_{s(v)}$, so $|\Delta(v)_{s(v)}|\leq\delta$. Then, following Eq. 9,

[TABLE]

We expand the first sum, cancel its first term, and parametrize the remaining terms by subsets $L\subsetneq V$ that correspond to those $v\in V$ for which we take the entries of $a(v)$ as opposed to $\Delta(v)$:

[TABLE]

Let $E(L)=\{(u,v)\in E:u,v\in L\}$. Recall that the value of $a(v)_{s(v)}$ depends only on the edges incident to $v$. Since $v\in V\setminus L$ in the second product, the value of $\Delta(v)_{s(v)}$ is not influenced by the edges in $E(L)$. We can thus decompose the string $s$ into two parts, $s_{1}$ and $s_{2}$, where $s_{2}$ is indexed by the edges in $E(L)$ and $s_{1}$ is indexed by the remaining edges:

[TABLE]

The bracketed expression is equal to $a(L)_{s(L)}$ which are the entries of the tensor obtained by contracting $L$ to a single vertex. Since $\lvert a(L)_{s(L)}\rvert\leq 1$ for all $s$, by triangle inequality we get that

[TABLE]

Since $|E-E(L)|\leq d^{\prime}\cdot|V-L|$, where $d^{\prime}$ is the maximum degree,

[TABLE]

Note that if $x\leq 1/r$ then

[TABLE]

Since $\delta\leq 1/(r8^{d^{\prime}})$,

[TABLE]

which completes the proof. ∎

I.7 Appendix VI: Proof of Theorem 3

Consider an $N$-qubit system with qubits grouped into $n$ parties. Let

[TABLE]

be a Hamiltonian on this system, where each term $H^{(i)}_{j}$ acts on at most two qubits and $i$ is the number of affected parties (the $i=1$ terms act within a single party while the $i=2$ terms act between two different parties). Denote by $g$ the $n$-vertex graph whose vertices represent the parties and whose edges correspond to interaction terms $H_{j}^{(2)}$ among different parties.

Starting with a product state $\rho_{1}\otimes\dotsb\otimes\rho_{n}$ across all parties and letting the system evolve for time $t$ with the Hamiltonian $H$, our goal is to approximate the expectation

[TABLE]

of local observables $O_{i}$ on each party.

Theorem 3**.**

Let $H$ be the Hamiltonian in Eq. 81 where each qubit is affected by at most a constant number of terms $d^{\prime}$. Assume that the states $\rho_{i}$ and observables $O_{i}$ in Eq. 82 can be implemented efficiently, and that the eigenvalues of each $O_{i}$ lie within $[-1,1]$. Then the correlation function in Eq. 82 can be approximated by an

[TABLE]

where $h_{B}=\sum_{j}\lVert H_{j}^{(2)}\rVert$ is the total inter-party interaction strength, $t\geq 0$ is the desired evolution time, $\operatorname{cc}(g)$ is the contraction complexity of the interaction graph $g$ between the parties, $\epsilon$ is the desired accuracy, and $d$ is the number of qubits of the largest party.

Proof.

The overall structure of the proof is as follows. We first implement the evolution of $H$ in the circuit model using a nested Trotter approximation based on Lemmas 6 and 7 (below). Then we use Lemma 8 to decompose the circuit and remove the interactions among parties so that each party can be simulated separately on a small quantum device.

Let us write the Hamiltonian as $H=A+B$ where $A$ and $B$ consist of all $H^{(1)}_{j}$ and $H^{(2)}_{j}$ terms, respectively. Using a nested Trotter approximation, we can decouple $H$ into these two groups and then into individual terms $H^{(i)}_{j}$.

First, let us decouple $A$ and $B$. Letting $m_{1}=\lceil 4d^{\prime}h_{B}t^{2}/\epsilon\rceil+\lceil 2h_{B}^{2}t^{2}/\epsilon\rceil$, Lemma 6 implies that

[TABLE]

Next, we approximate $e^{-iAt/m_{1}}$ and $e^{-iBt/m_{1}}$ by breaking them up into individual terms. For $e^{-iBt/m_{1}}$, we directly use Lemma 7:

[TABLE]

For $e^{-iAt/m_{1}}$, we first set

[TABLE]

where $h_{A}=\sum_{j}\lVert H_{j}^{(1)}\rVert$ and then use Lemma 7 and Eq. 120:

[TABLE]

Combining everything together, the overall Trotter approximation of $e^{-iHt}$ is given by

[TABLE]

and its incurred error is

[TABLE]

To construct a QC algorithm $(C,f)$ for approximating the correlation function in Eq. 82, we need to prepare the states $\rho_{i}$, approximate the unitary evolution by $H$, and measure the observables $O_{i}$. We implement these three steps as follows.

Recall that $N$ is the total number of qubits of all $n$ parties. Starting from $|0\rangle^{\otimes N}$, we first prepare the initial state $\rho=\rho_{1}\otimes\dotsb\otimes\rho_{n}$ by a sequence of single- and two-qubit gates (to account for mixed $\rho_{i}$, a random sequence of gates can be used.) By assumption, this can be done efficiently, and we ignore $\operatorname{poly}(N)$ factors in the simulator complexity.

Next, we can approximately implement $e^{-iHt}$ by applying the Trotter unitary $\tilde{U}$ from Eq. 88. Since each factor $e^{-iH_{j}^{(2)}t/m_{1}}$ and $e^{-iH_{j}^{(1)}t/(m_{1}m_{2})}$ can be implemented by a single- or two-qubit gate, $\tilde{U}$ consists of $O(Jm_{1}m_{2})=O(Jh_{A}^{2}t^{2}/\epsilon)=\operatorname{poly}(N)t^{2}/\epsilon$ gates, where $J$ is the total number of terms $H_{j}^{(i)}$ in the system (assume $J=\operatorname{poly}(N)$).

Finally, for each observable $O_{i}$ with spectral decomposition $O_{i}=\sum_{k}\alpha^{(i)}_{k}|\phi_{k}\rangle\langle\phi_{k}|$, we can find (again, by assumption) an efficiently implementable unitary $V_{i}$ such that

[TABLE]

where $|k\rangle$ is the standard basis. Let $f_{i}(k)=\alpha^{(i)}_{k}$. To implement $O=O_{1}\otimes\dotsb\otimes O_{n}$, we can first implement the basis change $V_{i}$ for each party, measure each party in the standard basis to obtain outcomes $y_{1},\ldots,y_{n}$, and then output $f(y)=\prod_{i=1}^{n}f_{i}(y_{i})$.

This procedure amounts to a QC algorithm $(C,f)$ whose output expectation is

[TABLE]

where $\tilde{U}$ is a Trotter approximation of $e^{-iHt}$. To bound the error from replacing $e^{-iHt}$ with $\tilde{U}$, observe that

[TABLE]

where the second-to-last inequality uses the assumption that $\|O\|\leq 1$.

The final step is to evoke Lemma 8 (below). Partitioning the qubits in the QC algorithm $(C,f)$ according to the parties in the Hamiltonian simulation, each party has at most $d$ qubits. We denote by $g^{\prime}$ the $n$-vertex graph whose vertices represent the parties and whose edges correspond to the two-qubit gates between different parties in the circuit $C$. Denote by $d_{g^{\prime}}(i,j)$ the number of edges between vertices $i$ and $j$ in $g^{\prime}$. Since each $H_{j}^{(2)}$ contributes $O(m_{1})=O(h_{B}^{2}t^{2}/\epsilon)$ gates to the Trotter approximation $\tilde{U}$,

[TABLE]

where $g$ is the interaction graph of the Hamiltonian $H$ and $d_{g}(i,j)$ is the number of edges between vertices $i$ and $j$ in $g$. Hence, the contraction complexities of $g^{\prime}$ and $g$ are related as follows:

[TABLE]

By Lemma 8, we conclude that a

[TABLE]

exists for the Hamiltonian $H$. ∎

Lemma 3** (Lemma 4 in Haah et al. (2018)).**

Let $A_{t}$ and $B_{t}$ be continuous time-dependent Hermitian operators, and let $U_{t}^{A}$ and $U_{t}^{B}$ with $U_{0}^{A}=U_{0}^{B}=I$ be the corresponding time evolution unitaries. Then the following hold:

• (a)

$W_{t}={U_{t}^{B}}^{\dagger}U_{t}^{A}$* is the unique solution of*

i\frac{dW_{t}}{dt}=\big{\lparen}{U_{t}^{B}}^{\dagger}(A_{t}-B_{t})U_{t}^{B}\big{\rparen}W_{t}* and $W_{0}=I$.*

• (b)

If $\lVert A_{s}-B_{s}\rVert\leq\delta$ for all $s\in[0,t]$, then

$\lVert U_{t}^{A}-U_{t}^{B}\rVert\leq t\delta$.

Lemma 4**.**

For Hermitian matrices $A$ and $B$,

[TABLE]

where $\mathrm{ad}_{B}$ is the adjoint map $\mathrm{ad}_{B}(A)=[B,A]$ and $\mathrm{ad}_{B}^{k}$ is the $k$-th power of $\mathrm{ad}_{B}$:

[TABLE]

Proof.

Note that $e^{-i(A+B)t}$ is generated by the Hamiltonian $A+B$ and $e^{-iBt}e^{-iAt}$ is generated by the Hamiltonian $e^{-iBt}(A+B)e^{iBt}$ by part (a) of Lemma 3. According to Lemma 5.3 in Miller (1972),

[TABLE]

Using this identity we get, for any $s\in[0,t]$,

[TABLE]

By part (b) of Lemma 3,

[TABLE]

which is the desired result. ∎

Lemma 5**.**

Let $H=A+B=\sum_{j}H^{(1)}_{j}+\sum_{j}H^{(2)}_{j}$ where each $H_{j}$ is a one- or two-qubit Hamiltonian with $\lVert H_{j}\rVert\leq 1$. If each qubit is involved in at most $d^{\prime}$ terms and $h_{B}=\sum_{j}\lVert H_{j}^{(2)}\rVert$ denotes the total interaction strength, then

[TABLE]

Proof.

Since $\lVert[X,Y]\rVert\leq 2\lVert X\rVert\lVert Y\rVert$ for any matrices $X$ and $Y$,

[TABLE]

To bound $\lVert[A,B]\rVert$, note that

[TABLE]

For a fixed $j$, assume $H_{j}^{(2)}$ acts on qubits $a$ and $b$. Then $[H_{i}^{(1)},H_{j}^{(2)}]$ can be non-zero only when $H_{i}^{(1)}$ also acts on $a$ or $b$. The number of such $H_{i}^{(1)}$ that do not commute with $H_{j}^{(2)}$ is bounded by $2d^{\prime}$. Hence

[TABLE]

Combining everything together we get $\lVert\mathrm{ad}_{B}^{k}(A)\rVert\leq(2h_{B})^{k-1}\cdot 4h_{B}d^{\prime}=2^{k+1}d^{\prime}h_{B}^{k}$. ∎

Lemma 6**.**

Assume the same setting as in Lemma 5. Then, for any positive $\epsilon\leq d^{\prime}t$,

[TABLE]

whenever $m\geq\lceil{4d^{\prime}h_{B}t^{2}/\epsilon\rceil}$.

Proof.

First, by combining Lemmas 4 and 5,

[TABLE]

where the last inequality follows from $\epsilon\leq d^{\prime}t$ and $1+\frac{1}{2}+\frac{1}{4}+\dotsc=2$. Note that $U^{m}-V^{m}=U^{m-1}(U-V)+U^{m-2}(U-V)V+\dotsb+(U-V)V^{m-1}$ so, for any unitary operations $U$ and $V$,

[TABLE]

where the last inequality follows from the triangle inequality and $\lVert U\rVert=\lVert V\rVert=1$. Using this,

[TABLE]

which is the desired bound. ∎

Lemma 7**.**

Let $H=\sum_{j}H_{j}$ and $h=\sum_{j}\lVert H_{j}\rVert$. For any positive $\epsilon\leq ht$,

[TABLE]

whenever $m\geq\lceil 2h^{2}t^{2}/\epsilon\rceil$.

Proof.

Let $L$ denote the number of terms in $H$. By direction expansion,

[TABLE]

If $m\geq\lceil 2h^{2}t^{2}/\epsilon\rceil$,

[TABLE]

as desired. ∎

Lemma 8**.**

Let $(C,f)$ be a QC algorithm and $\mathcal{Q}=\{Q_{1},\dotsc,Q_{r}\}$ be a partition of its qubits. Let $g^{\prime}$ be the $r$-vertex (multi-)graph obtained by representing each subset $Q_{i}$ as a vertex $i$ and each two-qubit gate in $C$ acting across sets $Q_{i}$ and $Q_{j}$ as an edge $(i,j)$. Let $K$ be the number of edges in $g^{\prime}$ and let $d=\max_{i}|Q_{i}|$ be the size of the largest set of qubits. Then $(C,f)$ has a $(2^{O(K)},d)$-simulator. If in addition $f$ is decomposable, i.e., $f(y)=\prod_{i=1}^{r}f_{i}(y_{i})$ where $y_{i}$ is the outcome of the standard basis measurement on qubits $Q_{i}$, and $f_{i}(y_{i})\in[-1,1]$ then $(C,f)$ has a $(2^{O(\operatorname{cc}(g^{\prime}))},d)$-simulator.

Proof.

Let $(G,\mathcal{A})$ be the tensor network of $(C,f)$. Our strategy is to construct a clustering of $G$ whose induced graph $g$ is very similar to the graph $g^{\prime}$ given in the Lemma. The two simulators are then obtained by applying Theorems 1 and 2, respectively.

We obtain the desired clustering of $G$ from the graph $g^{\prime}$. Since each cluster corresponds to either a vertex or an edge in $g^{\prime}$, we get two types of clusters. Let us describe them in more detail.

For each vertex $i$ of $g^{\prime}$, we introduce a cluster $S_{i}$ that contains tensors that act only on qubits $Q_{i}$. That is, $S_{i}$ consists of all vertices of $G$ that correspond to either an initial qubit state or a single- or a two-qubit gate in the set $Q_{i}$.

Similarly, we introduce a separate cluster also for each edge of $g^{\prime}$ (recall that the edges of $g^{\prime}$ correspond to two-qubit gates between different parts of $\mathcal{Q}$). We denote these clusters by $T_{j}$.

Since the number of edges of $g^{\prime}$ is $K$, the overall clustering of $G$ is given by

[TABLE]

where (as usual) the final observable $O_{f}$ has been left out of the clustering. This partition induces a graph $g$ 999Although we reuse the symbol $g$, here it has a slightly different meaning than in Theorem 3 where this Lemma is evoked. Here the graph $g$ is induced from a circuit $C$ while in Theorem 3 it is induced by interactions of a Hamiltonian $H$. The former is obtained from the latter via a Trotter approximation, see Eq. 95. that is very similar to $g^{\prime}$. In particular, the number of edges in $g$ (or the number of qubits communicated between clusters of $\mathcal{P}$) is $O(K)$.

Let us verify that the number of qubits needed for simulating each cluster of $\mathcal{P}$ is also $d$. Following the strategy of Theorem 1, we can repeatedly apply Lemma 1 to cut the tensor network into separate clusters. Since all interactions between qubits in different parts of $\mathcal{Q}$ are mediated by two-qubit gates $U_{j}$ represented by the clusters $T_{j}$, splitting the tensor network into clusters $\mathcal{P}$ effectively amounts to cutting these two-qubit gates $U_{j}$ out of the circuit. For each $U_{j}$, we apply Lemma 1 four times:

U_{j}$$=$$\displaystyle\sum_{k_{1},\dotsc,k_{4}=1}^{8}\!\!\!a_{k}$$O_{k_{1}}^{1}$$O_{k_{2}}^{2}$$\rho_{k_{1}}^{1}$$\rho_{k_{2}}^{2}$$O_{k_{3}}^{3}$$O_{k_{4}}^{4}$$\rho_{k_{3}}^{3}$$\rho_{k_{4}}^{4}$$U_{j}

where $a_{k}=c_{k_{1}}c_{k_{2}}c_{k_{3}}c_{k_{4}}\in\{-\frac{1}{2^{4}},\frac{1}{2^{4}}\}$ is the product of the four coefficients from Lemma 1 and $O_{k_{1}}^{1},\dotsc,O_{k_{4}}^{4}$ and $\rho_{k_{1}}^{1},\dotsc,\rho_{k_{4}}^{4}$ are the four associated observables and states, respectively. Note that we can simulate each $U_{j}$ classically since only two qubits are involved. This takes care of all clusters $T_{j}$.

The remaining clusters $S_{i}$ are associated to subsets of qubits $Q_{i}$. Let us argue that $|Q_{i}|$ qubits are enough for simulating $S_{i}$ if qubit recycling is permitted. Indeed, for each input qubit of $S_{i}$ that results from cutting out some two-qubit gate $U_{j}$ (e.g., $\rho_{k_{3}}^{3}$ or $\rho_{k_{4}}^{4}$ in the above figure), there always exists an output qubit of $S_{i}$ (e.g., the qubit measured by $O_{k_{1}}^{1}$ or $O_{k_{2}}^{2}$ in the above figure) that can be recycled to be the new input qubit after the measurement. Hence, no additional qubits – other than the original $|Q_{i}|$ ones – are needed for simulating the cluster $S_{i}$.

To summarize, the original parameters $K$ and $d$ of the graph $g^{\prime}$ defined in the statement of the Lemma translate into $O(K)$ and $d$ for the graph $g$ induced by the partition $\mathcal{P}$ of $G$. Hence, the circuit $C$ is $(O(K),d)$-clustered and thus by Theorem 1 has a $(2^{O(K)},d)$-simulator.

If $f$ is decomposable, we can get a $(2^{O(\operatorname{cc}(g))},d)$-simulator by using Theorem 2. We can find a good contraction order for $g$ based on $g^{\prime}$. Indeed, we can first absorb each cluster $T_{j}$ that represents some two-qubit gate $U_{j}$ into one of its adjacent clusters $S_{i}$ or $S_{i^{\prime}}$ – this effectively introduces a double edge between $S_{i}$ and $S_{i^{\prime}}$. Repeating this for each $T_{j}$ reduces the graph $g$ to $g^{\prime}$. Hence $\operatorname{cc}(g)=O(\operatorname{cc}(g^{\prime}))$, yielding a $(2^{O(\operatorname{cc}(g^{\prime}))},d)$-simulator for $(C,f)$. ∎

I.8 Appendix VII: VQE experiments

We consider VQE algorithms using the following parameterized circuit with $n$ qubits and $D$ layers proposed in Kandala et al. (2017):

[TABLE]

where $U_{i}(\theta_{i})=\bigotimes_{j=1}^{n}U_{i}^{j}(\theta_{i}^{j})$ and each $U_{i}^{j}(\theta_{i}^{j})$ is a parameterized single-qubit gate applied on the $j$-th qubit:

[TABLE]

where $Z_{j}(\beta)=\exp(-i\beta Z)$ and $X_{j}(\beta)=\exp(-i\beta X)$ are single-qubit $X$ and $Z$ rotations applied to the $j$-th qubit. Since $U_{0}(\theta_{0})$ is applied directly to the input $|0\rangle^{n}$, the first $Z$ rotation gates can be removed:

[TABLE]

so in total $\theta$ has $(3D+2)n$ parameters. Each layer of single-qubit gates in Eq. 133 is followed by the following layer of entangling gates:

[TABLE]

where $\mathrm{CZ}(i,j)$ denotes the controlled $Z$ gate between qubits $i$ and $j$, defined as $\mathrm{CZ}|a,b\rangle=(-1)^{ab}|a,b\rangle$, for all $a,b\in\{0,1\}$.

We run the VQE circuit from Kandala et al. (2017) with $n=6$ and $D=1$ on the 5-qubit “ibmq_ourense” cloud quantum computer provided by the IBM Quantum Experience IBM (2020). Our experiment uses up to three qubits. Our goal is to solve the optimization problem $\min_{\theta}F(\theta)$ with

[TABLE]

where $\theta$ has $(3+2)\cdot 6=30$ parameters and the minimal eigenvalue of $H$ approximates the ground energy of $\text{BeH}_{2}$ with an interatomic distance of 1.7 Å 101010We choose $H$ as a diagonal matrix, which both provides low measurement complexity and desired approximation precision of the energy..

We optimize the parameters of $\theta$ using the simultaneous perturbation stochastic approximation (SPSA) method. Let $\theta(k)$ denote the parameters at the $k$-th iteration of SPSA. In each iteration, a random vector $r(k)$ is generated, where each element of $r(k)$ is drawn uniformly from $\{1,-1\}.$ Let $\theta^{\pm}(k)=\theta(k)\pm c_{k}r(k)$, where $c_{k}$ are hyper-parameters. The gradient of $\theta(k)$ is estimated by

[TABLE]

where $F$ is evaluated by the quantum device. Then, $\theta$ is updated as follows:

[TABLE]

where $a_{k}$ controls the learning rate. We set $c_{k}=0.3/\sqrt{k}$ and $a_{k}=k^{-0.3}$, and perform $200$ iterations.

In Kandala et al. (2017), $F(\theta)$ is evaluated by running $U(\theta)$ on a $6$-qubit quantum computer. Using our cluster simulation scheme, we evaluate $F$ using only $3$ qubits. To do so, we partition the six qubits of $U(\theta)$ into two sets: $\{1,2,3\}$ and $\{4,5,6\}.$ We then decompose the only two-qubit gate CZ$(3,4)$ that acts between these two sets using the gadget proposed in Mitarai and Fujii (2019) (see Fig. 6 111111We have corrected typos that appear in Fig. 2 of Mitarai and Fujii (2019).). Note that $Z$, $RZ(\pi/2)$, and $(I+aZ)/2$ are all single-qubit operations 121212To implement $(I+aZ)/2$ on the IBM Quantum Experience, we directly measure the qubit in the $Z$-basis and do classical post-processing., so the circuits on $\{1,2,3\}$ and $\{4,5,6\}$ can be run separately using only $3$ qubits. For simulating each run of $U(\theta)$, we need to execute $12$ different circuits on a $3$-qubit computer. For each of these circuits, we produce $8000$ samples to estimate their output. The value of $F(\theta)$ is then obtained by aggregating and post-processing these outputs.

The resulting accuracy of our estimate is similar to the one obtained using the full $6$-qubit circuit Kandala et al. (2017), including the evaluation of $\theta^{-}(k),\theta^{+}(k)$ and $\theta(k)$. See Fig. 5. This demonstrates the potential of our cluster simulation scheme for decreasing the number of qubits required in near-term VQE applications. Our source code can be found at: https://github.com/TianyiPeng/Partiton_VQE.

I.9 Appendix VIII: Reducing $\operatorname{cc}(g)$ in VQE

To reduce the number of circuits required to run, we propose the following steps when applying our scheme to the VQE circuit in Eq. 133

choose a partition $\mathcal{P}=\{S_{1},\dotsc,S_{r}\}$ of $n$ qubits such that $\lvert S_{i}\rvert\leq d$ for each $i$;

• 2.

remove some entangling gates from $U_{\text{ENT}}$ that go across different parts of $\mathcal{P}$ and denote the modified sequence by $U_{\text{ENT}}^{\prime}$;

• 3.

replace some occurrences of $U_{\text{ENT}}$ in the original circuit $U(\theta)$ by $U_{\text{ENT}}^{\prime}$ and denote the modified circuit by $U^{\prime}(\theta)$.

Given the circuit $U^{\prime}(\theta)$ with fixed parameters $\theta$, we estimate $\langle 0|^{\otimes n}U^{\prime{\dagger}}(\theta)HU^{\prime}(\theta)|0\rangle^{\otimes n}$ by performing Pauli measurements and summing according to the Pauli decomposition of $H$:

[TABLE]

where $\sigma_{i,j}\in\{X,Y,Z,I\}$ is a Pauli matrix. Since the measurement is Pauli observable, we can use a QC circuit $(C,f)$ to compute the expectation value, where $f$ is decomposable. Denote $g$ by regarding each set $S_{i}$ as a node and each gate that acts cross two sets as an edge. By Lemma 8, the modified circuit has a $(2^{O(\operatorname{cc}(g))},d)$-simulator where $\operatorname{cc}(g)$ is the contraction complexity of the graph $g$.

This leads to a trade-off in time versus performance: replacing $U_{\text{ENT}}$ by $U_{\text{ENT}}^{\prime}$ decreases $\operatorname{cc}(g)$ and hence reduces the simulation time; however, this also decreases the expressive power of the circuit and thus degrades the quality of the solution. This may be justified in physically meaningful cases, e.g.,when the true ground state obeys an area law Eisert et al. (2010).

Our numerical observations suggest that for a small scale quantum circuit, replacing some $U_{\text{ENT}}$ by $U_{\text{ENT}}^{\prime}$ does not influence the performance significantly. Specifically, we let $U_{\text{ENT}}$ be a fixed sequence of two-qubits gates for producing entanglement:

[TABLE]

where $\mathrm{CNOT}(i,j)$ is the CNOT gate with the $i$-th qubit as control and the $j$-th qubit as target.

We partition our $n$ qubits into two sets: $S_{1}=\{1,\dotsc,\frac{n}{2}\}$ and $S_{2}=\{\frac{n}{2}+1,\dotsc,n\}$ and construct $U_{\text{ENT}}^{\prime}$ form $U_{\text{ENT}}$ by removing CNOT gates that cross the two sets:

[TABLE]

Replacing $U_{\text{ENT}}$ by $U_{\text{ENT}}^{\prime}$ decreases the number of gates across the two sets.

Let $R\subseteq\{1,\dotsc,D\}$ denote the set of locations where we replace $U_{\text{ENT}}$ by $U_{\text{ENT}}^{\prime}$ and let $D_{1}$ denote the number of remaining $U_{\text{ENT}}$ gates. In our numerical experiments, $n=6$ and we compare three different regimes:

• •

$D=9,D_{1}=9,R=\{\}$;

• •

$D=9,D_{1}=3,R=\{1,2,4,6,8,9\}$;

• •

$D=3,D_{1}=3,R=\{\}$.

For each run of the experiment, we generate a Hamiltonian

[TABLE]

where $\alpha_{i}$ is drawn uniformly at random from $[-1,1]$ and $\sigma_{i,j}$ is a random Pauli matrix drawn uniformly from the set $\{I,X,Y,Z\}$, with all random variables drawn independently. We set the number of terms to be $L=50$. After generating the Hamiltonian $H$, we execute $10^{5}$ steps of the simultaneous perturbation stochastic approximation (SPSA) method to determine the parameters $\theta$. At step $t$, we record the relative error of the current expectation value and the optimal eigenvalue:

[TABLE]

where $\theta(t)$ are the current parameter values and $v_{\mathrm{opt}}=\min_{|\phi\rangle}\langle\phi|H|\phi\rangle$ is the smallest eigenvalue of $H$ obtained by explicitly diagonalizing the Hamiltonian. We run each experiment $100$ times and average the relative error over these runs. The code can be found at https://github.com/TianyiPeng/Partiton_VQE.

We report our numerical observations in Fig. 7. Our results show that, for the problem at hand, replacing some $U_{\text{ENT}}$ by $U_{\text{ENT}}^{\prime}$ does not greatly influence the performance of the VQE. This suggests that our scheme is worth pursuing when implementing VQE with limited quantum memory. Because of its heuristic nature, theoretical analysis of our scheme is difficult and awaits future exploration.