1 Introduction↩︎

Classical strong simulation of quantum circuits remains indispensable: it forms a basis for circuit compilation and optimization pipelines, aids the study of noise in actual hardware, and helps triangulate truly quantum-hard regimes that are infeasible for classical machines. The difficulty lies in the fact that an \(n\)-qubit state is described by \(2^n\) amplitudes; therefore, a useful simulator must exploit structure rather than represent the state naively.

Tensor-network contraction gives a standard route to strong simulation of quantum circuits. If \(N_C\) is the tensor-network graph of a circuit \(C\), with tensor vertices and shared indices as edges, then a contraction order is naturally governed by the treewidth of the line graph \(L(N_C)\) of \(N_C\). We write \[\operatorname{cc}(N_C)=\operatorname{tw}(L(N_C))\] for this contraction complexity. The Markov–Shi simulation bound is single-exponential in the corresponding tensor-network width parameter [1] and linear in the circuit size.

The route we develop instead is the sum-of-powers (SOP): the expression of an amplitude as a sum over Feynman paths. This representation is classical [2], but it has recently been turned into competitive simulators by encoding it with knowledge-representation tools from artificial intelligence: e.g. binary decision diagrams, in FeynmanDD [3], [4], and weighted model counting [5]–[8]. Fixing the circuit’s input \(\boldsymbol{y}\) and output \(\boldsymbol{z}\) turns the circuit into a scalar \[Z=\frac{1}{R}\sum_{\boldsymbol{x}\in\{0,1\}^V}\omega_r^{f(\boldsymbol{x})},\] a sum of roots of unity over assignments to a set \(V\) of Boolean path variables; the exponent \(f\) is a degree-two polynomial: a constant, one linear term per variable, and a quadratic sign term for each interacting pair of variables (for \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\) the modulus is \(r=8\), see 3). The interacting pairs form the SOP variable graph \(G_C=(V,E)\), and the difficulty of evaluating \(Z\) is concentrated in the structure of this graph.

Our main result is that the right structural measure of \(G_C\) is its rank-width, a branching width parameter from structural graph theory. We give a dynamic program (4) that evaluates \(Z\), hence the amplitude, in time exponential only in the rank-width of \(G_C\) and polynomial in the circuit size.

1.0.0.1 Contributions.

We establish the following based on parameters in 1:

1. We give an explicit rank-decomposition dynamic program for quadratic SOPs, running in time \(O(nr4^k\operatorname{poly}(n)+r\log r)\) on an \(n\)-variable graph of rank-width \(k\) and modulus \(r\) (2, 3), when provided with a rank-width decomposition. Applied to circuits over Hadamard and diagonal gates, this is a rank-width fixed-parameter strong simulation algorithm (6). On the way, we also expose a treewidth fixed-parameter strong simulation algorithm by combining the insight that \(G_C\) is the “primal graph” of the SOP with classical graphical-model algorithms (3.2).

2. We relate the SOP variable graph to the parameters of the tensor-network and decision-diagram approaches: \(G_C\) is a minor of the tensor line graph, giving \(\operatorname{tw}(G_C)\le\operatorname{cc}(N_C)\) (1) and \(\operatorname{rw}(G_C)\le\operatorname{cc}(N_C) + 1\) [9]; and since rank-width is at most linear rank-width, we recover the decision-diagram guarantee of [4] (1) as a special case.

3. We construct circuit families with \(\operatorname{rw}(G_C)=O(1)\) while \(\operatorname{lrw}(G_C)\) and \(\operatorname{cc}(N_C)\) are unbounded (2), explicitly separating the rank-width, linear-rank-width, and treewidth approaches. On these families, our rankwidth-parameterized algorithm runs with fixed polynomial overhead, whereas decision-diagram and tensor-network guarantees both degrade.

Although we state the construction for the \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\) gate set, it requires only that the Hadamard gate is the sole non-diagonal single-qubit gate. In fact, arbitrary diagonal single-qubit gates contribute linear terms but leave \(G_C\) unchanged. The results therefore cover every circuit over Hadamard, \(\mathrm{CZ}\), and arbitrary phase rotation gates, and in particular all of Clifford\(+\mathrm{T}\) (6.2).

1.0.0.2 Other related work.

Weighted model counting has recently been proposed as a general framework for quantum physics computations, including partition functions [5]; the present work uses such encodings as one practical route and contributes a rank-width fixed-parameter algorithm for the underlying SOP. Recent and independent work has also identified rank-width as a useful parameter for quantum-circuit simulation, mainly through ZX-calculus and stabilizer-decomposition methods. Kuyanov and Kissinger [10] give a contraction algorithm for graph-like ZX diagrams whose claimed running time is \(\widetilde{O}(4^R)\) given a rank-decomposition of width \(R\). Codsi and Laakkonen [11] develop a framework connecting graph-width measures and stabilizer decompositions, obtaining slice-wise polynomial-time algorithms parameterized by treewidth or rank-width of an associated tensor network together with the number of non-Clifford operations. These results are complementary to ours: they show, from a ZX/tensor-network perspective, that rank-width is a natural structural parameter for simulation. Our contribution is different in that we work directly with the quadratic SOP representation underlying FeynmanDD-style simulation; give an explicit DP that runs on a rank-decomposition of the SOP variable graph \(G_C\); and relate this graph to both graphical-model treewidth and tensor-network contraction complexity. In particular, \(\mathrm{T}\) and other single-qubit diagonal gates enter the SOP only as unary phases and leave \(G_C\) unchanged, so the number of non-Clifford gates affects our running time only polynomially, through the circuit size, and never in the exponent.

2 Preliminaries↩︎

This section collects the two ingredients the rest of the paper relies on: the graph width parameters that will measure simulation cost and the quantum circuits we simulate together with their tensor-network representation. A reader already comfortable with these notions can skim and refer back as needed.

2.1 Graph notation and parameters↩︎

We first collect the graph notation and width parameters used throughout the paper. All graphs are finite, undirected, and simple unless explicitly stated otherwise. For a graph \(G=(V,E)\) and sets \(S,T\subseteq V\), \(\boldsymbol{A}_G[S,T]\) denotes the submatrix of the adjacency matrix of \(G\) with rows indexed by \(S\) and columns indexed by \(T\). Whenever ranks of adjacency submatrices are used, the field is \(\mathbb{F}_2\).

2.1.0.1 Treewidth.

A tree decomposition of a graph \(G=(V,E)\) is a pair \((T,\mathcal{B})\) where \(T\) is a tree with vertex set \(V(T)\) and \(\mathcal{B}=(B_t)_{t\in V(T)}\) is a family of bags \(B_t\subseteq V\) satisfying the following conditions.

1. Every vertex appears in some bag: \(V=\bigcup_{t\in V(T)}B_t\).

2. Every edge is covered by some bag: for every \(uv\in E\), there is a node \(t\) with \(u,v\in B_t\).

3. Bags containing a fixed vertex form a connected subtree: for every \(v\in V\), the set \(\{t\in V(T)\mid v\in B_t\}\) is connected in \(T\).

The width of \((T,\mathcal{B})\) is \(\max_{t\in V(T)}|B_t|-1\). The treewidth of \(G\) is \[\operatorname{tw}(G)=\min_{(T,\mathcal{B})}\left(\max_{t\in V(T)}|B_t|-1\right),\] where the min ranges over all tree decompositions of \(G\). We use the standard facts that treewidth is monotone under taking minors and that \(\operatorname{tw}(K_t)=t-1\). Computing treewidth exactly is -hard in general. However, for every fixed integer \(k\), one can decide whether \(\operatorname{tw}(G)\le k\) and, if so, construct a tree decomposition of width at most \(k\) in polynomial time, indeed in linear time for fixed \(k\) [12].

2.1.0.2 Rank-width.

Let \(G=(V,E)\) be a graph. For \(X\subseteq V\), define the cut-rank of \(X\) in \(G\) by \[\rho_G(X)=\operatorname{rank}_{\mathbb{F}_2}\boldsymbol{A}_G[X,V\setminus X].\] A rank-decomposition of \(G\) is a pair \((T,\mu)\) where \(T\) is a subcubic tree and \(\mu\) is a bijection from the leaves of \(T\) to \(V\). Here, subcubic simply means that the maximal degree of a node in \(T\) is three.

Each edge \(e\) of \(T\) induces a bipartition, or cut, of \(V\). Indeed, deleting \(e\) splits \(T\) into two connected components. The leaves in these two components correspond, through \(\mu\), to two disjoint subsets \(X_e\) and \(\overline{X}_e\) of \(V\). Thus, \(V = X_e \uplus \overline{X}_e\), where \(\uplus\) denotes the disjoint union.

The width of an edge \(e\) in \(T\) is then \(\rho_G(X_e)\), the width of the decomposition is the maximum width of its edges, and finally, the rank-width of \(G\) is \[\operatorname{rw}(G)=\min_{(T,\mu)}\max_{e\in E(T)}\rho_G(X_e),\] where \(E(T)\) is edge set of \(T\).

Fixed-width rank decompositions can be found, if they exist, in polynomial time for each fixed width [13]–[15].

Linear rank-width is the linear-layout analogue of rank-width. A linear layout of \(G\) is an ordering \(\pi=(v_1,\ldots,v_n)\) of \(V\). Its width is \(\max_{1\le i<n}\rho_G(\{v_1,\ldots,v_i\}).\) The linear rank-width of \(G\) is \[\operatorname{lrw}(G)=\min_{\pi}\max_{1\le i<n}\rho_G(\{v_1,\ldots,v_i\}).\] Equivalently, it is rank-width restricted to caterpillar-shaped rank decompositions [15]. That is, when the decomposition tree \(T\) in the rank-decomposition \((T,\mu)\) is restricted to be a caterpillar: a tree whose nonleaf vertices form a path. Hence \[\operatorname{rw}(G)\le \operatorname{lrw}(G).\] The inequality can be strict on natural graph families. In particular, Adler and Kanté prove that the linear rank-width of a forest equals its pathwidth; complete binary trees therefore have bounded rank-width but linear rank-width growing with their height [16]. We will use this fact later to argue our new algorithm is provably better than previous approaches on some families of circuits.

2.2 Quantum circuits and gate sets↩︎

We recall quantum circuits and fix the elementary gate set on which we phrase our results. For an \(n\)-qubit system, the state space is \((\mathbb{C}^2)^{\otimes n} \cong \mathbb{C}^{2^n}\), with computational basis states \(\ket{x} = \ket{x_1\cdots x_n}\), \(x_i\in\{0,1\}\). A general \(n\)-qubit state is therefore \(\ket{\psi}=\sum_{x\in\{0,1\}^n} a_x \ket{x}\) with \(a_x\in\mathbb{C}\), and for normalized states \(\sum_{x\in\{0,1\}^n} |a_x|^2 = 1\). A quantum gate is a unitary linear operator acting on one or more qubits; when a gate acts on only a subset of the \(n\) qubits, it is understood to act as the identity on the remaining ones. A quantum circuit is a finite sequence of such gates acting on a fixed number of qubits.

2.2.0.1 Elementary gate set.

The circuits we consider are built over the set of elementary gates \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\), with standard matrices \[\mathrm{H}=\frac{1}{\sqrt2}\begin{bmatrix}1&1\\1&-1\end{bmatrix}, \qquad \mathrm{T}=\begin{bmatrix}1&0\\0&e^{i\pi/4}\end{bmatrix}, \qquad \mathrm{CZ}=\begin{bmatrix} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&-1 \end{bmatrix}.\] The gate set \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\) is universal for quantum computation and is also known as the Clifford + T gate set [17], where we note that the \(\mathrm{T}\) gate renders one Clifford gate, namely \(\text{S} = \mathrm{T}^2\), redundant.

2.2.0.2 Circuits.

We write an \(n\)-qubit circuit as \[C = U_m\cdots U_2U_1,\] where the gates are ordered as applied: \(U_1\) first, then \(U_2\), and so on up to \(U_m\). We refer to the index \(i\) as the position of the gate \(U_i\) in the circuit. The action of \(C\) on an input state \(\ket{\psi}\) is \(C\ket{\psi} = U_m\cdots U_2U_1\ket{\psi}\), and the matrix of the whole circuit is obtained by multiplying the gate matrices in reverse textual order, since the rightmost gate acts first.

To specify which qubits a gate acts on, we use the standard subscript convention: \(\mathrm{H}_a\) and \(\mathrm{T}_a\) denote a single-qubit Hadamard or \(\mathrm{T}\) gate applied to qubit \(a\) (and the identity on every other wire), and \(\mathrm{CZ}_{a,b}\) denotes a \(\mathrm{CZ}\) gate applied to qubits \(a\) and \(b\) (and the identity elsewhere); since \(\mathrm{CZ}\) is symmetric in its two arguments, the order of \(a\) and \(b\) does not matter. For a qubit \(a\), the Hadamard depth \(k_a\) is the total number of \(\mathrm{H}\) gates that act on qubit \(a\) in \(C\). This count indexes the path variables on wire \(a\) in the sum-of-powers construction of 3.

2.2.0.3 Strong simulation.

Throughout the paper, by classical strong simulation of a circuit \(C\) we mean the task of computing, for given input and output computational basis states \(\ket{\boldsymbol{y}}\) and \(\ket{\boldsymbol{z}}\) with \(\boldsymbol{y},\boldsymbol{z}\in\{0,1\}^n\), the amplitude \[\braket{\boldsymbol{z}|C|\boldsymbol{y}}\in\mathbb{C},\] i.e.the entry of the unitary associated with \(C\) in row \(\boldsymbol{z}\) and column \(\boldsymbol{y}\). This is the quantity whose squared modulus \(|\braket{\boldsymbol{z}|C|\boldsymbol{y}}|^2\) is the probability of observing outcome \(\boldsymbol{z}\) when \(C\) is run on input \(\boldsymbol{y}\) and measured in the computational basis. Strong simulation is therefore strictly more demanding than weak simulation, which only requires producing samples from this output distribution.

2.3 Tensor networks↩︎

The standard route to strong simulation contracts the circuit’s tensor network, and we will measure our SOP-based algorithm against it (3.2); we therefore recall tensor networks and the contraction complexity that governs their cost.

A tensor of order \(k\) is a multilinear map \(A\colon [d_1]\times\cdots\times[d_k]\to\mathbb{C}\), where \([d]=\{0,\ldots,d-1\}\); each of the \(k\) argument positions is an index, and we write \(A[a_1,\ldots,a_k]\) for the entry of \(A\) at \((a_1,\ldots,a_k)\). A tensor network is a finite collection of tensors \(\mathcal{T}=\{A_1,\ldots,A_m\}\) together with a pairing of some of their index positions across different tensors. Each matched pair (whose two positions must share a dimension) is a single bond (or contracted) index, and the unmatched positions are the open (or free) indices; write \(b_1,\ldots,b_p\) for the bond indices and \(o_1,\ldots,o_\ell\) for the open ones.

A bond configuration \(\boldsymbol{\beta}=(\beta_1,\ldots,\beta_p)\) assigns a value \(\beta_i\in[d_{b_i}]\) to every bond index \(b_i\). Together with a choice \(\boldsymbol{j}=(j_1,\ldots,j_\ell)\) of values for the open indices, it fixes a value for every index position of every tensor, so each \(A_t\) has a well-defined entry. The value of the network is the order-\(\ell\) tensor obtained by summing, over all bond configurations, the product of these entries: \[V[j_1,\ldots,j_\ell] =\sum_{\beta_1,\ldots,\beta_p} \prod_{t=1}^m A_t[\boldsymbol{\beta},\boldsymbol{j}],\] where \(A_t[\boldsymbol{\beta},\boldsymbol{j}]\) denotes the entry of \(A_t\) read at the values that \((\boldsymbol{\beta},\boldsymbol{j})\) assigns to the indices of \(A_t\). For circuits every index is a qubit wire, so each dimension is \(d=2\).

2.3.0.1 Tensor-network graph.

Given a tensor network \(\mathcal{T}\), its tensor-network graph (a multigraph) \(N_{\mathcal{T}}\) has one vertex per tensor and one edge \(\{t,t'\}\) for each bond index shared between \(A_t\) and \(A_{t'}\); open indices are not represented by edges.

It is sometimes convenient to represent a quantum circuit \(C=U_m\cdots U_1\) over \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\) as a tensor network. In this context, each gate is a tensor: \(\mathrm{H}\) and \(\mathrm{T}\) are order-\(2\) tensors (one input, one output qubit index) and \(\mathrm{CZ}\) is order-\(4\) (two input, two output qubit indices). Adjacent gates sharing a wire are connected by a bond index; the circuit’s input and output wires are the open indices. The tensor-network graph \(N_C\) thus has one vertex per gate and one edge per wire segment connecting two gates.

2.3.0.2 Contraction complexity.

A contraction order is a sequence of pairwise tensor contractions that reduces the network to a single tensor. A contraction step merges two tensors and sums out their shared indices, producing a new tensor whose bond indices with the rest of the network are the union of those of the two; its cost is exponential in the number of indices of this intermediate tensor. This process is captured by the line graph \(L(N_C)\), whose vertices are the bond indices of \(N_C\) (its edges) and whose edges record co-occurrence in a common tensor: a contraction step corresponds exactly to a vertex-elimination step, in which the shared indices are removed and the remaining indices on the merged tensor form a clique. A full contraction order is therefore an elimination order on \(L(N_C)\), and Markov and Shi [1] showed that the minimum-cost order is governed by \(\operatorname{tw}(L(N_C))\). This motivates the definition \[\operatorname{cc}(N_C)\mathrel{\vcenter{:}}=\operatorname{tw}(L(N_C))\] of the contraction complexity of a circuit.

3 Quadratic sums-of-powers for circuits↩︎

In 2.2, we defined strong simulation as the task of evaluating the amplitude \(\braket{\boldsymbol{z}|C|\boldsymbol{y}}\). We now recall how, for circuits over \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\), this amplitude is a quadratic sum-of-powers (SOP): a scalar sum of roots of unity over Feynman paths [2]–[4]. The SOP, and the graph \(G_C\) it induces on its variables, are the central objects of the paper. After deriving them, we show in 3.2 that \(G_C\) is never more complex than the circuit’s tensor network, so already at this stage SOP evaluation is competitive with tensor-network contraction; the rank-width refinement that beats it comes in 4.

3.1 Sum of powers form and associated graph↩︎

For a circuit \(C = U_m \cdots U_1\), inserting computational-basis resolutions of the identity between consecutive gates gives \[\braket{\boldsymbol{z}| C | \boldsymbol{y}} = \sum_{\boldsymbol{y}_1, \ldots, \boldsymbol{y}_{m-1}} \prod_{j=1}^{m} \braket{\boldsymbol{y}_j | U_j | \boldsymbol{y}_{j-1}}, \text{ where } \boldsymbol{y}_0 = \boldsymbol{y}\text{ and } \boldsymbol{y}_m = \boldsymbol{z}.\]

3.1.0.1 Local gate factors.

For the gate set \(\{\mathrm{H}, \mathrm{T}, \mathrm{CZ}\}\), every nonzero local matrix entry is a power of \(\omega_8 = e^{2\pi i / 8}\), up to the \(1/\sqrt{2}\) normalization carried by each Hadamard. For \(p, q, p_1, p_2, q_1, q_2 \in \{0,1\}\), \[\label{eqn:deltas} \begin{align} \braket{p | \mathrm{H}| q} &{}= 2^{-1/2} (-1)^{pq} = 2^{-1/2} \omega_8^{4pq},\\ \braket{p | \mathrm{T}| q} &{}= \delta_{p,q} \omega_8^{p}, \text{ and }\\ \braket{p_1 p_2 | \mathrm{CZ}| q_1 q_2} &{}= \delta_{p_1, q_1} \delta_{p_2, q_2} \omega_8^{4 p_1 p_2}. \end{align}\tag{1}\] Thus a Hadamard contributes a quadratic phase \(4pq\) together with a factor \(1/\sqrt{2}\); a \(\mathrm{T}\)-gate contributes a unary phase \(p\); and a \(\mathrm{CZ}\)-gate contributes a quadratic phase \(4 p_1 p_2\). All exponents are taken modulo \(8\).

3.1.0.2 Path variables.

The Kronecker deltas in the \(\mathrm{T}\) and \(\mathrm{CZ}\) entries (see 1 ) force each summand to vanish unless the value on every affected wire is unchanged across the gate, so the incoming and outgoing indices across such a diagonal gate may be identified. A Hadamard has no such delta, so the values immediately before and after it remain independent. Equivalently, cut each qubit wire at every Hadamard: on wire \(a\) this yields \(k_a+1\) maximal segments, where \(k_a\) is the number of Hadamards on wire \(a\), carrying Boolean variables \(x_{a,0},\ldots,x_{a,k_a}\). Writing \(S\) for the set of all segment variables and reading off the contributions above,

• a \(\mathrm{T}\)-gate on segment \(s\) contributes \(x_s\);

• a \(\mathrm{CZ}\)-gate touching segments \(s, t\) contributes \(4 x_s x_t\);

• a Hadamard from segment \(s\) to its successor \(s'\) contributes \(4 x_s x_{s'}\) and one factor \(1/\sqrt{2}\).

Summing these local exponents gives a raw quadratic phase polynomial \[f(\boldsymbol{x}) = \sum_{\mathrm{T}\text{ on } s} x_s + 4 \sum_{\mathrm{H}: s \to s'} x_s x_{s'} + 4 \sum_{\mathrm{CZ}\text{ on } s, t} x_s x_t \pmod{8},\] with all sums counted with multiplicity and all coefficients reduced modulo \(8\).

3.1.0.3 Pinning the boundary.

To compute \(\braket{\boldsymbol{z}| C | \boldsymbol{y}}\) we impose the boundary conditions on \(f\). Let \(s_{\mathrm{in}}(a)\) and \(s_{\mathrm{out}}(a)\) be the first and last segments on wire \(a\), and pin \(x_{s_{\mathrm{in}}(a)} = y_a\) and \(x_{s_{\mathrm{out}}(a)} = z_a\) for \(a = 1, \ldots, n\). If these substitutions assign two different values to the same segment variable, then no Feynman path is compatible with the boundary data and \(\braket{\boldsymbol{z}| C | \boldsymbol{y}} = 0\). Otherwise, let \(V \subseteq S\) be the set of unpinned segment variables and write \(f|_{\boldsymbol{y}, \boldsymbol{z}}\) for the polynomial obtained from \(f\) after substituting the pinned variables.

3.1.0.4 The pinned SOP.

Collecting terms in \(f|_{\boldsymbol{y}, \boldsymbol{z}}\), the constant terms combine into a single phase \(c \in \mathbb{Z}_8\), the surviving one-variable terms become coefficients \(b_v \in \mathbb{Z}_8\) (the \(\mathrm{T}\)-gate phases together with any sign interaction having exactly one pinned endpoint), and the surviving free–free quadratic interactions become the edges \(E\) of a graph on the unpinned variables \(V\). Writing \(f_C := f|_{\boldsymbol{y}, \boldsymbol{z}}\), \[\begin{gather} f_C(\boldsymbol{x}) = c + \sum_{v \in V} b_v\, x_v + 4 \sum_{\{u,v\} \in E} x_u x_v \pmod{8},\\ \braket{\boldsymbol{z}| C | \boldsymbol{y}} = \frac{1}{R} \sum_{\boldsymbol{x}\in \{0,1\}^V} \omega_8^{f_C(\boldsymbol{x})}, \end{gather}\] with \(R = (\sqrt{2})^{m_\mathrm{H}}\) and \(m_\mathrm{H}\) the number of Hadamards in \(C\). We call \(G_C = (V, E)\) the SOP variable graph of the pinned amplitude \((C, \boldsymbol{y}, \boldsymbol{z})\).

3.1.0.5 General quadratic SOP.

We now state the SOP and its associated counting problem in their general form, decoupled from any particular circuit. Let \(G=(V,E)\) be a graph, \(r\) an even positive integer, \(c\in\mathbb{Z}_r\), and \(b_v\in\mathbb{Z}_r\) for every \(v\in V\). Consider the quadratic SOP \[\label{eq:sop} \begin{align} Z ={}& \frac{1}{R}\sum_{\boldsymbol{x}\in\{0,1\}^V}\omega_r^{f(\boldsymbol{x})}, \text{ with}\\ f(\boldsymbol{x})={} &c+\sum_{v\in V}b_vx_v+\eta\sum_{\{u,v\}\in E}x_ux_v\pmod r, \end{align}\tag{2}\] where \(\omega_r\) is the \(r\)th root of unity. Throughout the paper the edge coefficient is fixed as \(\eta=r/2\); this is the value produced by the construction above for \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\) (where \(r=8\) and \(\eta=4\)), and is what makes each cross-term contribute only a sign \((-1)^{x_ux_v}\).

The SOP value is then recovered via \[Z=\frac{1}{R}\sum_{j=0}^{r-1}N_j\omega_r^j.\]

For a circuit over \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\) this is the instance with \(r=8\) and \(G=G_C\), the coefficients \(b_v\) and constant \(c\) being those identified above.

3.2 SOP evaluation matches contraction complexity↩︎

Recall from 2.3 the tensor-network graph \(N_C\) of a circuit, its line graph \(L(N_C)\), and the contraction complexity \(\operatorname{cc}(N_C)=\operatorname{tw}(L(N_C))\) governing the Markov–Shi tensor-contraction bound [1]. We relate the SOP variable graph \(G_C\) of 3 to this quantity. The resulting bound \(\operatorname{tw}(G_C)\le\operatorname{cc}(N_C)\) lets us compare the rank-width FPT dynamic program of 4 with tensor contraction (5), and it yields a treewidth FPT dynamic program for SOP evaluation directly, as explained below.

Since we only have fixed-arity tensor networks (that arise from quantum circuits) one may also compare to the treewidth of \(N_C\) itself [1]. If every tensor has arity at most \(\Delta\), then a tree decomposition of \(N_C\) of width \(t\) can be converted into one of \(L(N_C)\) of width at most \(\Delta(t+1)-1\) by replacing each bag with all incident tensor indices. Without bounded arity, this implication is false: a star with one high-arity tensor has treewidth one as a tensor graph, but a line graph that is a clique.

Thus \(\operatorname{tw}(G_C)\le\operatorname{cc}(N_C)\): the SOP variable graph is never harder than the tensor-network contraction. Quadratic SOP evaluation is moreover a special case of the standard sum–product problem on a graphical model: the variables are Boolean and the unary and binary factors are exactly the coefficients of 2 , so the SOP variable graph \(G_C\) is the primal graph of the instance. By bucket elimination (equivalently, tree-decomposition dynamic programming) such instances can be evaluated in time \(2^{O(\operatorname{tw}(G_C))}\operatorname{poly}(n)\) [18], [19]. The link between path-sum formulas and graphical-model primal graphs is also noted in [4]; here 1, together with the observation that \(G_C\) is the primal graph, lets us invoke this framework directly to obtain a treewidth FPT dynamic program for SOP evaluation that already matches the Markov–Shi contraction-complexity guarantee. 4 refines the parameter from treewidth to the rank-width of \(G_C\).

4 Rank-based dynamic program for SOPs↩︎

This section gives the paper’s main algorithm: an exact evaluation of a quadratic SOP that is fixed-parameter tractable in the rank-width of its variable graph \(G_C\). Since rank-width can be much smaller than treewidth or linear rank-width, this is what lets the SOP route outperform both tensor contraction and decision-diagram simulators on the circuit families of 5. We give a dynamic-programming algorithm that processes a rooted rank-decomposition \((\mathcal{T},\mu)\) of \(G_C\) from the leaves to the root in 4.1. Then, in 4.2 we propose an optimization based on Fourier-modes.

4.1 SOP evaluation by rank-decomposition↩︎

The algorithm processes a rank-decomposition of \(G_C\) from the leaves upward, keeping at each node a compact table that records how the partial assignments below it interact with the rest of the graph. Because each cut of a width-\(k\) decomposition has rank at most \(k\), only \(2^k\) interaction patterns can occur, so the tables stay small—this is the source of the fixed-parameter guarantee.

4.1.0.1 Setup.

We work with a rooted rank-decomposition \((\mathcal{T},\mu)\) of \(G_C\): \(\mathcal{T}\) is a tree of maximum degree \(3\) whose leaves are in bijection with \(V\) via \(\mu\), with one designated internal node \(r_{\mathcal{T}}\) as the root (obtained, if necessary, by subdividing an arbitrary edge of an unrooted decomposition and rooting at the new vertex). For a node \(u\) of \(\mathcal{T}\), let \(X_u\subseteq V\) be the set of vertices corresponding to leaves of the subtree rooted at \(u\), and let \(\overline{X}_u=V\setminus X_u\). At the root, \(X_{r_{\mathcal{T}}}=V\) and \(\overline{X}_{r_{\mathcal{T}}}=\varnothing\); at a leaf \(u\) with \(\mu(u)=v\), \(X_u=\{v\}\). At each node \(u\), the DP table records how assignments inside \(X_u\) look to the outside.

For \(\boldsymbol{z}\in\{0,1\}^{X_u}\), define the boundary signature \[\sigma_u(\boldsymbol{z})=\boldsymbol{z}^\top\boldsymbol{A}[X_u,\overline{X}_u]\in\mathbb{F}_2^{\overline{X}_u}.\] This vector records, for each outside vertex, the parity of its selected neighbors inside \(X_u\). Define the internal residue \[\phi_u(\boldsymbol{z})= \sum_{v\in X_u} b_vz_v + \eta\sum_{\{v,w\}\in E,\;v,w\in X_u} z_vz_w \pmod r.\] The constant \(c\) is not included in \(\phi_u\) and is added back at the root.

The table at node \(u\) is stored sparsely. Its keys are pairs \((\boldsymbol{\sigma},s)\), where \(s\in\mathbb{Z}_r\) and \(\boldsymbol{\sigma}\) is an occurring boundary signature, and it stores \[D_u[\boldsymbol{\sigma},s] = \#\{\boldsymbol{z}\in\{0,1\}^{X_u}\mid \sigma_u(\boldsymbol{z})=\boldsymbol{\sigma},\;\phi_u(\boldsymbol{z})=s\}.\] Although \(\sigma_u(\boldsymbol{z})\) is written as a vector in \(\mathbb{F}_2^{\overline{X}_u}\), it always lies in the image of the linear map \[\{0,1\}^{X_u}\rightarrow \mathbb{F}_2^{\overline{X}_u}, \qquad \boldsymbol{z}\mapsto \boldsymbol{z}^\top\boldsymbol{A}[X_u,\overline{X}_u],\] namely the row space of \(\boldsymbol{A}[X_u,\overline{X}_u]\) over \(\mathbb{F}_2\). For every nonroot node \(u\), the edge from \(u\) to its parent in the rooted rank-decomposition induces the cut \((X_u,\overline{X}_u)\). Since the decomposition has width \(k\), \(\operatorname{rank}_{\mathbb{F}_2}\boldsymbol{A}[X_u,\overline{X}_u]\le k\) and therefore at most \(2^k\) boundary signatures can occur at \(u\). At the root, \(\overline{X}_u=\varnothing\), so there is only the empty signature. Thus each table has at most \(r2^k\) nonzero states.

4.1.0.2 Leaves.

If \(u\) is a leaf corresponding to vertex \(v\), then \(X_u=\{v\}\). The assignment \(z_v=0\) contributes \[D_u[0,0]\mathrel{+}=1,\] and the assignment \(z_v=1\) contributes \[D_u[\boldsymbol{A}[v,V\setminus\{v\}],b_v]\mathrel{+}=1.\] The use of \(\mathrel{+}=\) matters: the two states may coincide.

4.1.0.3 Joins.

Let \(t\) be an internal node with children \(L\) and \(R\). Write \[X_t=X_L\uplus X_R, \qquad Y=\overline{X}_t=V\setminus X_t.\] Then \(\overline{X}_L=X_R\cup Y\) and \(\overline{X}_R=X_L\cup Y\). Take a left state \((\boldsymbol{\alpha},p)\) and a right state \((\boldsymbol{\beta},q)\), with \[\boldsymbol{\alpha}\in\mathbb{F}_2^{X_R\cup Y}, \qquad \boldsymbol{\beta}\in\mathbb{F}_2^{X_L\cup Y}.\] Let \(\boldsymbol{a}\in\{0,1\}^{X_L}\) and \(\boldsymbol{b}\in\{0,1\}^{X_R}\) realize these states. The new quadratic terms are exactly the edges between \(X_L\) and \(X_R\). Their parity is \[\chi(\boldsymbol{\alpha},\boldsymbol{\beta})=\boldsymbol{a}^\top\boldsymbol{A}[X_L,X_R]\boldsymbol{b}\in\mathbb{F}_2.\] This value is well-defined.

The parent signature and residue are \[\boldsymbol{\gamma}=\boldsymbol{\alpha}|_Y+\boldsymbol{\beta}|_Y\in\mathbb{F}_2^Y, \qquad s=p+q+\eta\chi(\boldsymbol{\alpha},\boldsymbol{\beta})\pmod r.\] The join does not create an unbounded set of parent signatures. If \(\boldsymbol{\alpha}\) and \(\boldsymbol{\beta}\) are realized by assignments \(\boldsymbol{a}\in\{0,1\}^{X_L}\) and \(\boldsymbol{b}\in\{0,1\}^{X_R}\), then \[\boldsymbol{\gamma} = \boldsymbol{a}^\top\boldsymbol{A}[X_L,Y]+\boldsymbol{b}^\top\boldsymbol{A}[X_R,Y] = (\boldsymbol{a},\boldsymbol{b})^\top\boldsymbol{A}[X_t,Y] = \sigma_t(\boldsymbol{a},\boldsymbol{b}).\] Thus, every parent signature produced by the join lies in the row space of \(\boldsymbol{A}[X_t,Y]\). Since \(Y=\overline{X}_t\), the cut \((X_t,Y)\) is again a cut induced by an edge of the rooted rank-decomposition unless \(t\) is the root. Therefore this row space has size at most \(2^k\) at every nonroot internal node, and size \(1\) at the root.

The join update is \[D_t[\boldsymbol{\gamma},s]\mathrel{+}=D_L[\boldsymbol{\alpha},p]D_R[\boldsymbol{\beta},q].\]

4.2 Fourier-mode implementation↩︎

The residue table of the previous subsection carries an avoidable factor \(r\) in the join. We remove it by working with the discrete Fourier transform of the residue counts on \(\mathbb{Z}_r\), bringing the running time to the \(O(nr4^k\operatorname{poly}(n)+r\log r)\) bound quoted in 6.

4.2.0.1 Residue counts and their Fourier transform.

The vector of residue counts \(N=(N_0,N_1,\ldots,N_{r-1})\) is a function on the cyclic group \(\mathbb{Z}_r\). Using the characters \(\chi_a(j)=\omega_r^{aj}\) for \(a,j\in\mathbb{Z}_r\), define the discrete Fourier transform \[\label{eq:residue-dft} \widehat N(a)=\sum_{j=0}^{r-1}N_j\chi_a(j) =\sum_{j=0}^{r-1}N_j\omega_r^{aj}.\tag{3}\] This is the usual finite Fourier transform; the same transform and its inverse underlie the fast Fourier transform [20]. Here the group is only \(\mathbb{Z}_r\), but the same principle is useful in what follows: instead of storing coefficients indexed by residues, we store evaluations against the \(\chi_a\).

3 is also exactly an amplitude evaluation. Grouping the assignments \(\boldsymbol{x}\) by their residue \(f(\boldsymbol{x})\pmod r\) gives \[\widehat N(a) =\sum_{j=0}^{r-1}\#\{\boldsymbol{x}\in\{0,1\}^V\mid f(\boldsymbol{x})=j\pmod r\}\omega_r^{aj} =\sum_{\boldsymbol{x}\in\{0,1\}^V}\omega_r^{a f(\boldsymbol{x})},\] since every assignment \(\boldsymbol{x}\) contributes to exactly one residue class. In particular, the original SOP value satisfies \(\widehat N(1)=\sum_{\boldsymbol{x}}\omega_r^{f(\boldsymbol{x})}=RZ\), so amplitude evaluation gives one Fourier coefficient and evaluating the scaled polynomials \(af\) for all \(a\in\mathbb{Z}_r\) gives all of them. By orthogonality of characters of \(\mathbb{Z}_r\), the transform inverts as \[\label{eq:inverse-residue-dft} N_t=\frac{1}{r}\sum_{a=0}^{r-1}\omega_r^{-at}\widehat N(a), \qquad t\in\mathbb{Z}_r.\tag{4}\] Computing all residue counts and computing all scaled amplitude evaluations are therefore equivalent up to an \(r\)-point Fourier transform, performed in \(O(r\log r)\) when \(r\) has a suitable factorization (and in \(O(r^2)\) by direct inversion otherwise — lower-order in the fixed-\(r\) regime).

4.2.0.2 Per-mode dynamic program.

This suggests running the rank-decomposition DP one Fourier mode at a time: for each \(a\in\mathbb{Z}_r\), track only one complex number per boundary signature, computing \(\widehat N(a)\) in a single pass, and recover all \(N_j\) at the root via 4 . Concretely, for each \(a\in\mathbb{Z}_r\) define \[A_u^{(a)}[\boldsymbol{\sigma}] = \sum_{\substack{\boldsymbol{z}\in\{0,1\}^{X_u}\\ \sigma_u(\boldsymbol{z})=\boldsymbol{\sigma}}} \omega_r^{a\phi_u(\boldsymbol{z})}.\] At the root, \[\widehat N(a)=\sum_{\boldsymbol{x}\in\{0,1\}^V}\omega_r^{a f(\boldsymbol{x})} =\omega_r^{ac}A_{r_{\mathcal{T}}}^{(a)}[\varnothing],\] and the counts are recovered by \[\label{eq:inverse-dft} N_j=\frac{1}{r}\sum_{a=0}^{r-1}\omega_r^{-aj}\widehat N(a) =\frac{1}{r}\sum_{a=0}^{r-1}\omega_r^{a(c-j)}A_{r_{\mathcal{T}}}^{(a)}[\varnothing].\tag{5}\]

For a fixed mode \(a\), the leaf initialization becomes \[A_u^{(a)}[0]\mathrel{+}=1, \qquad A_u^{(a)}[\boldsymbol{A}[v,V\setminus\{v\}]]\mathrel{+}=\omega_r^{ab_v}.\] At a join, \[A_t^{(a)}[\boldsymbol{\gamma}] \mathrel{+}= A_L^{(a)}[\boldsymbol{\alpha}] A_R^{(a)}[\boldsymbol{\beta}] \omega_r^{a\eta\chi(\boldsymbol{\alpha},\boldsymbol{\beta})}.\] Since \(\eta=r/2\), \[\omega_r^{a\eta\chi}=\omega_r^{a(r/2)\chi}=(-1)^{a\chi}.\] Thus the interaction kernel depends only on the parity of \(a\): even modes have no signed cross term, and odd modes have the same \((-1)^{\chi}\) kernel.

4.2.0.3 Implementation details.

All of the Fourier modes share the same rank-decomposition, signature spaces, restriction maps, parent-signature maps, and cross-term table \(\chi(\boldsymbol{\alpha},\boldsymbol{\beta})\). These maps can be precomputed once for each join node. The only mode-dependent leaf data are the unary phases \(\omega_r^{ab_v}\), and the only mode-dependent join data are the two kernels given by the parity of \(a\). In practice one therefore batches even modes and odd modes separately. The asymptotic bound remains \(O(nr4^k\operatorname{poly}(n))\), but the constant factors are smaller than in a literal execution of \(r\) independent programs.

5 Breaking through the linear rank-width and contraction complexity barriers↩︎

The dynamic program above is parameterized by the rank-width of the SOP variable graph. This is a branching parameter: the decomposition tree may split the variables recursively. Linear decision-diagram approaches [4] use a variable ordering and are therefore governed by a linear layout parameter, such as linear rank-width. Tensor-network contraction is governed by the contraction complexity \(\operatorname{cc}(N_C)=\operatorname{tw}(L(N_C))\). The following construction separates these parameters at the level of the SOP variable graph, so the advantage of 4 is not an artifact of a particular benchmark but provable on explicit circuit families.

We first establish that for an arbitrary simple graph we can construct a quantum circuit such that its SOP variable graph is the graph input graph. This avoids relying on a special feature of a particular benchmark family (cf. the separation from [4]).

We next define a family of graphs with bounded rank-width, unbounded linear rank-width, and arbitrarily large treewidth. Let \(B_h\) be the complete binary tree of height \(h\). For an integer \(t\ge 1\), let \[\Gamma_{h,t}=B_h[K_t]\] be the graph obtained from \(B_h\) by replacing every vertex \(v\) by a clique \(K_t^v\) of \(t\) twins, and by replacing every tree edge \(\{u,v\}\) by the complete bipartite graph between \(K_t^u\) and \(K_t^v\). Equivalently, \(\Gamma_{h,t}\) is the lexicographic blow-up of \(B_h\) by a \(t\)-clique.

The case \(t=1\) already separates rank-width from linear rank-width: the SOP variable graph is a complete binary tree, so the rank-width is bounded while the linear rank-width grows with the height. Increasing \(t\) makes the tensor-network contraction complexity grow as well, without changing the rank-width bound. For instance, choosing \(t=2^h\) gives circuits of size \(\Theta(2^h t^2)\) whose rank-width parameter is constant and whose contraction complexity is \(\Omega(t)\). This is a separation between the structural parameters used by the three algorithmic routes; it is not a lower bound against simulators using additional structure outside these parameters.

6 Consequences for quantum-circuit simulation↩︎

We now turn from the algorithm to what it means for simulation. We state the rank-width guarantee for circuit amplitudes (6.1), show that it covers far more than \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\), and place it alongside the decision-diagram and model-counting simulators (6.3).

6.1 Rank-width FPT simulation↩︎

Applying the rank-decomposition dynamic program of 4 to the SOP variable graph \(G_C\) of a circuit \(C\) over \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\) evaluates the amplitude \(\langle\boldsymbol{z}|C|\boldsymbol{y}\rangle\). By [3], [4] the amplitude equals the SOP value \(Z\) on \(G_C\), so 3 computes it in time \(O(|C|\,r4^k\operatorname{poly}(|C|)+r\log r)\) given a rank-decomposition of \(G_C\) of width \(k\); for a fixed gate set and modulus this is \(2^{O(k)}\operatorname{poly}(|C|)\), fixed-parameter tractable in the rank-width of \(G_C\). The algorithm uses only the graph, its coefficients, and the rank-decomposition, so the same bound holds for any quadratic SOP, not only circuit amplitudes.

6.2 Diagonal gates and the Clifford\(+\mathrm{T}\) gate set↩︎

The path-sum construction of 3 uses only two features of the gate set, and both extend well beyond \(\{\mathrm{H},\mathrm{T},\mathrm{CZ}\}\); consequently the results of the paper apply to a substantially larger class of circuits.

The first feature is that \(\mathrm{H}\) is the only nondiagonal gate. As recalled in 3, writing \(\mathrm{H}_{p,q}=\tfrac{1}{\sqrt2}(-1)^{pq}\) shows that each \(\mathrm{H}\) pairs the path value \(x\) before it with a fresh value \(x'\) after it, contributing the quadratic sign \((-1)^{xx'}=\omega_r^{(r/2)xx'}\) and a scalar \(1/\sqrt2\), while \(\mathrm{CZ}\) contributes the same quadratic sign between two wires. We keep the standard normalized \(\mathrm{H}\): the \(1/\sqrt2\) factors are constants, collected once per Hadamard into the normalization \(R\). These quadratic signs are exactly the edges of the SOP variable graph \(G_C\).

The second feature is that \(\mathrm{T}\) contributes only a diagonal phase. Replacing it by an arbitrary diagonal single-qubit gate \[\mathrm{D}_{\alpha,\beta}=\begin{bmatrix}\alpha&0\\0&\beta\end{bmatrix}\] changes nothing structurally: being diagonal, \(\mathrm{D}_{\alpha,\beta}\) leaves the value on its wire unchanged and contributes the unary factor \(\alpha\) at \(x=0\) and \(\beta\) at \(x=1\). It thus couples no wires and adds no edge to \(G_C\), only a unary weight. Folding the scalar \(\alpha\) into \(R\) leaves weights \(1\) and \(\beta/\alpha\); when these phases are commensurate—powers of a common root of unity \(\omega_r\)—they are exactly the unary coefficients \(\omega_r^{b_v}\) of 2 . Diagonal gates of this form include the phase gates \(\mathrm S\), \(\mathrm{T}\), and \(\mathrm Z\), global-phase gates, and \(\mathrm Z\)-rotations through rational multiples of \(\pi\).

Hence every circuit over \(\{\mathrm{H},\mathrm{D}_{\alpha,\beta},\mathrm{CZ}\}\) produces, after pinning the input and output basis states, a quadratic SOP of the form 2 : the \(\mathrm{H}\) and \(\mathrm{CZ}\) gates give the cross-terms with coefficient \(r/2\), the diagonal gates give the unary coefficients \(b_v\), and any interaction whose endpoints are both pinned folds into the constant \(c\). Since \(\mathrm{CZ}\) and \(\mathrm{CX}\) differ only by local Hadamards (2.2), this gate set realizes \(\{\mathrm{H},\mathrm{T},\mathrm{CX}\}\); and because the Clifford group is generated by \(\mathrm{H}\), the phase gate \(\mathrm S=\mathrm{D}_{1,i}\), and \(\mathrm{CX}\) [17], it realizes all of Clifford\(+\mathrm{T}\), whose phases are \(8\)th roots of unity so that again \(r=8\). Every result of the paper then holds verbatim for such circuits—the rank-width dynamic program (2, 3), the treewidth and contraction-complexity bounds of 3.2, and the simulation guarantee of 6.1—since these depend only on the quadratic SOP and its variable graph \(G_C\), not on the specific gates that produced it.

6.3 Knowledge-representation simulators↩︎

The SOP is exactly the kind of object that knowledge-representation methods evaluate, and our results pin down what each costs. Decision-diagram simulators [3], [4] process the variables in a fixed linear order, so on the SOP core they are governed by linear rank-width; since \(\operatorname{rw}(G)\le\operatorname{lrw}(G)\) with an unbounded gap on natural families (2), the rank-width dynamic program of 4 strictly strengthens this dependence while recovering it as a special case (1).

The SOP can equally be handed to a weighted model counter [5]–[8]. The standard encoding introduces one Boolean sign variable \(s_{uv}\) per edge \(\{u,v\}\in E\), constrained by \(s_{uv}\leftrightarrow(x_u\wedge x_v)\) and weighted \(-1\) when true, so that summing it out reproduces the cross-term sign \((-1)^{x_ux_v}\) of 2 ; the amplitude is then a weighted model count of the resulting formula \(F_C\). Each sign variable is local to one edge, so the treewidth of the formula’s incidence graph, i.e. the bipartite variable-clause graph \(I(F_C)\) whose edges represent variable inclusion in clauses, stays within a constant of that of \(G_C\).

Combined with \(\operatorname{tw}(G_C)\le\operatorname{cc}(N_C)\) (1), a standard model-counting dynamic program [22], [23] therefore simulates \(C\) in time \(2^{O(\operatorname{cc}(N_C))}\operatorname{poly}(|C|)\), matching Markov–Shi. Rank-width, by contrast, must be measured on \(G_C\) and not on \(F_C\). This is because the sign variables are harmless for treewidth but not for rank-width. Concretely, \(K_n\) has rank-width \(1\), while the incidence graph of its encoding forces rank-width linear in \(n\). Model counting can thus be seen as a route that exposes the SOP, but the rank-width algorithm must run on the variable graph itself.

6.3.0.1 Take-home message for practical matters:

General-purpose decision-diagram and model-counting engines can serve as the simulation workhorse, and the SOP structure lets one dispatch a specialized fixed-parameter algorithm (i.e. the new rank-width dynamic program of 4, or treewidth bucket elimination) whenever the corresponding width of \(G_C\) is small. Established weighted model counters, like Ganak [24], do implement heuristics based on treewidth decompositions. Developing improved heuristics based on treewidth and rank-width decompositions for simulation via model counting is an interesting direction for future work.

References↩︎

---

1. Contact: guillermo.perez@uantwerpen.be↩︎

2. This work was supported by the FWO/NWO KR2iQS project (G076326N), a collaboration between the University of Antwerp and Leiden University.↩︎