Simulating quantum computations with Tutte polynomials

Abstract

We establish a classical heuristic algorithm for exactly computing quantum probability amplitudes. Our algorithm is based on mapping output probability amplitudes of quantum circuits to evaluations of the Tutte polynomial of graphic matroids. The algorithm evaluates the Tutte polynomial recursively using the deletion–contraction property while attempting to exploit structural properties of the matroid. We consider several variations of our algorithm and present experimental results comparing their performance on two classes of random quantum circuits. Further, we obtain an explicit form for Clifford circuit amplitudes in terms of matroid invariants and an alternative efficient classical algorithm for computing the output probability amplitudes of Clifford circuits.

Introduction

There is a natural relationship between quantum computation and evaluations of Tutte polynomials^1,2. In particular, quantum probability amplitudes are proportional to evaluations of the Tutte polynomial of graphic matroids. In this paper, we use this relationship to establish a classical heuristic algorithm for exactly computing quantum probability amplitudes. While this problem is known to be #P-hard in general³, our algorithm focuses on exploiting structural properties of an instance to achieve an improved runtime over traditional methods. Previously it was known that this problem can be solved in time exponential in the treewidth of the underlying graph⁴.

The basis of our algorithm is a mapping between output probability amplitudes of quantum circuits and evaluations of the Tutte polynomial of graphic matroids^2,5,6. Our algorithm proceeds to evaluate the Tutte polynomial recursively using the deletion–contraction property. At each step in the recursion, our algorithm computes certain structural properties of the matroid in order to attempt to prune the computational tree. This approach to computing Tutte polynomials was first studied by Haggard et al.⁷. Our algorithm can be seen as an adaption of their work to special points of the Tutte plane where we can exploit additional structural properties.

The performance of algorithms for computing Tutte polynomials based on the deletion–contraction property depends on the heuristic used to decide the ordering of the recursion^7,8,9. We consider several heuristics introduced by Pearce et al.⁸ and an additional heuristic, which is specific to our algorithm. We present some experimental results comparing the performance of these heuristics on two classes of random quantum circuits corresponding to dense and sparse instances.

The correspondence between output probability amplitudes of quantum circuits and evaluations of Tutte polynomials also allows us to obtain an explicit form for Clifford circuit amplitudes in terms of matroid invariants by a theorem of Pendavingh¹⁰. This gives rise to an alternative efficient classical algorithm for computing output probability amplitudes of Clifford circuits.

This paper is structured as follows. We introduce matroid theory in “Matroid theory” and the Tutte polynomial in “The Tutte polynomial.“ In “The Potts model partition function,” “Instantaneous quantum polynomial time,” and “Quantum computation and the Tutte polynomial,” we establish a mapping between output probability amplitudes of quantum circuits and evaluations of the Tutte polynomial of graphic matroids. This is achieved by introducing the Potts model partition function in “The Potts model partition function,” instantaneous quantum polynomial (IQP) time circuits in “Instantaneous quantum polynomial time,” and a class of universal quantum circuits in “Quantum computation and the Tutte polynomial.” In “Efficient classical simulation of Clifford circuits,” we use this mapping to obtain an explicit form for Clifford circuit amplitudes in terms of matroid invariants. We also obtain an efficient classical algorithm for computing the output probability amplitudes of Clifford circuits. We describe our algorithm in “Algorithm overview” and present some experimental results in “Experimental results.” Finally, we conclude in “Discussion.”

Results

Matroid theory

We shall now briefly introduce the theory of matroids. The interested reader is referred to the classic textbooks of Welsh¹¹ and Oxley¹² for a detailed treatment. Matroids were introduced by Whitney¹³ as a structure that generalises the notion of linear dependence. There are many equivalent ways to define a matroid. We shall define a matroid by the independence axioms.

Definition 1 (matroid)

A matroid is a pair $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ consisting of a finite set ${{{\mathcal{S}}}}$, known as the ground set, and a collection ${{{\mathcal{I}}}}$ of subsets of ${{{\mathcal{S}}}}$, known as the independent sets, such that the following axioms are satisfied:

• (1) The empty set is a member of ${{{\mathcal{I}}}}$.

• (2) Every subset of a member of ${{{\mathcal{I}}}}$ is a member of ${{{\mathcal{I}}}}$.

• (3) If A and B are members of ${{{\mathcal{I}}}}$ and ∣A∣ > ∣B∣, then there exists an x ∈ A\B such that B ∪ {x} is a member of ${{{\mathcal{I}}}}$.

The rank of a subset A of ${{{\mathcal{S}}}}$ is given by the rank function $r:{2}^{{{{\mathcal{S}}}}}\to {\mathbb{N}}$ of the matroid defined by $r(A):=\max \left(| X| | X\subseteq A,X\in {{{\mathcal{I}}}}\right)$. The rank of a matroid M, denoted r(M), is the rank of the set S.

The archetypal class of matroids is vector matroids. A vector matroid $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ is a matroid whose ground set ${{{\mathcal{S}}}}$ is a subset of a vector space over a field ${\mathbb{F}}$ and whose independent sets ${{{\mathcal{I}}}}$ are the linearly independent subsets of ${{{\mathcal{S}}}}$. The rank of a subset of a vector matroid is the dimension of the subspace spanned by the corresponding vectors. We say that a matroid is ${\mathbb{F}}$-representable if it is isomorphic to a vector matroid over the field ${\mathbb{F}}$. A matroid is a binary matroid if it is ${{\mathbb{F}}}_{2}$-representable and is a ternary matroid if it is ${{\mathbb{F}}}_{3}$-representable. A matroid that is representable over every field is called a regular matroid.

Every finite graph G = (V,E) induces a matroid $M(G)=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ as follows. Let the ground set ${{{\mathcal{S}}}}$ be the set of edges E and let the independent sets ${{{\mathcal{I}}}}$ be the subsets of E that are a forest, i.e., they do not contain a simple cycle. It is easy to check that M(G) satisfies the independence axioms. The rank of a subset A of a cycle matroid is ∣V∣ − κ(A), where κ(A) denotes the number of connected components of the subgraph with edge set A. The rank of the cycle matroid M(G), denoted r(M(G)) or simply r(G), is the rank of the set E. The matroid M(G) is called the cycle matroid of G. We say that a matroid is graphic if it is isomorphic to the cycle matroid of a graph.

Graphic matroids are regular. To see this, consider assigning to the graph G an arbitrary orientation D(G), that is, for each edge e = {u,v} in G, we choose one of u and v to be the positive end and the other one to be the negative end. Then, construct the oriented incidence matrix of G with respect to the orientation D(G).

Definition 2 (oriented incidence matrix)

Let G = (V,E) be a graph and let D(G) be an orientation of G. Then, the oriented incidence matrix of G with respect to D(G) is the ∣V∣ × ∣E∣ matrix ${A}_{D(G)}={({a}_{ve})}_{| V| \times | E| }$ whose entries are

$${a}_{ve}=\left\{\begin{array}{ll}+1,&{{{\rm{if}}}}\ v\ {{{\rm{is}}}}\ {{{\rm{the}}}}\ {{{\rm{positive}}}}\ {{{\rm{end}}}}\ {{{\rm{of}}}}\ e;\\ -1,&{{{\rm{if}}}}\ v\ {{{\rm{is}}}}\ {{{\rm{the}}}}\ {{{\rm{negative}}}}\ {{{\rm{end}}}}\ {{{\rm{of}}}}\ {{{\rm{e}}}};\\ 0,&{{{\rm{otherwise}}}}.\end{array}\right.$$

(1)

The rows of the oriented incidence matrix A_D(G) correspond to the vertices of G and the columns correspond to the edges of G. Each column contains exactly one +1 and exactly one −1 representing the positive and negative ends of the corresponding edge. If the column space of A_D(G) is the ground set of a vector matroid, then it is easy to see that a subset is independent if and only if it is a forest in G. Hence, the oriented incidence matrix provides a representation of a graphic matroid over every field.

A minor of a matroid M is a matroid that is obtained from M by a sequence of deletion and contraction operations.

Definition 3 (deletion)

Let $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ be a matroid and let e be an element of the ground set. Then, the deletion of M with respect to e is the matroid $M\setminus \{e\}=({{{\mathcal{S}}}}^{\prime} ,{{{\mathcal{I}}}}^{\prime} )$ whose ground set is ${{{\mathcal{S}}}}^{\prime} ={{{\mathcal{S}}}}\setminus \{e\}$ and whose independent sets are ${{{\mathcal{I}}}}^{\prime} =\{I\subseteq {{{\mathcal{S}}}}\setminus \{e\}| I\in {{{\mathcal{I}}}}\}$.

The deletion of an element from the cycle matroid of a graph corresponds to removing an edge from the graph.

Definition 4 (contraction)

Let $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ be a matroid and let e be an element of the ground set. Then, the contraction of M with respect to e is the matroid $M/\{e\}=({{{\mathcal{S}}}}^{\prime} ,{{{\mathcal{I}}}}^{\prime} )$ whose ground set is ${{{\mathcal{S}}}}^{\prime} ={{{\mathcal{S}}}}\setminus \{e\}$ and whose independent sets are ${{{\mathcal{I}}}}^{\prime} =\{I\subseteq {{{\mathcal{S}}}}\setminus \{e\}| I\cup \{e\}\in {{{\mathcal{I}}}}\}$.

The contraction of an element from the cycle matroid of a graph corresponds to removing an edge from the graph and merging its two endpoints.

An element e of a matroid is said to be a loop if {e} is not an independent set and said to be a coloop if e is contained in every maximally independent set. If an element e of a matroid is either a loop or a coloop then the deletion and contraction of e are equivalent.

The Tutte polynomial

We shall now briefly introduce the Tutte polynomial, which is a well-known invariant in matroid and graph theory.

Definition 5 (Tutte polynomial of a matroid)

Let $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ be a matroid with rank function $r:{2}^{{{{\mathcal{S}}}}}\to {\mathbb{N}}$. Then, the Tutte polynomial of M is the bivariate polynomial defined by

$${{{\rm{T}}}}(M;x,y):=\mathop{\sum}\limits_{A\subseteq {{{\mathcal{S}}}}}{(x-1)}^{r(M)-r(A)}{(y-1)}^{| A| -r(A)}.$$

(2)

The Tutte polynomial may also be defined recursively by the deletion–contraction property.

Definition 6 (deletion–contraction property)

Let $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ be a matroid. If M is the empty matroid, i.e., ${{{\mathcal{S}}}}=\varnothing$, then

$${{{\rm{T}}}}(M;x,y)=1.$$

(3)

Otherwise, let e be an element of the ground set. If e is a loop, then

$${{{\rm{T}}}}(M;x,y)=y{{{\rm{T}}}}(M\setminus \{e\};x,y).$$

(4)

If e is a coloop, then

$${{{\rm{T}}}}(M;x,y)=x{{{\rm{T}}}}(M/\{e\};x,y).$$

(5)

Finally, if e is neither a loop nor a coloop, then

$${{{\rm{T}}}}(M;x,y)={{{\rm{T}}}}(M\setminus \{e\};x,y)+{{{\rm{T}}}}(M/\{e\};x,y).$$

(6)

The deletion–contraction property immediately gives an algorithm for recursively computing the Tutte polynomial. This algorithm is in general inefficient, but the performance may be improved by using isomorphism testing to reduce the number of recursive calls¹⁴. The performance of this algorithm depends on the heuristic used to choose elements of the ground set^7,8,9. Björklund et al.¹⁵ showed that the Tutte polynomial can be computed in time exponential in the number of vertices.

The Tutte polynomial of a graph may be recovered by considering the Tutte polynomial of the cycle matroid of a graph and using the fact that the rank of a subset A of a cycle matroid is ∣V∣ − κ(A), where κ(A) denotes the number of connected components of the subgraph with edge set A.

Definition 7 (Tutte polynomial of a graph)

Let G = (V,E) be a graph and let κ(A) denote the number of connected components of the subgraph with edge set A. Then, the Tutte polynomial of G is a polynomial in x and y, defined by

$${{{\rm{T}}}}(G;x,y):=\mathop{\sum}\limits_{A\subseteq E}{(x-1)}^{\kappa (A)-\kappa (E)}{(y-1)}^{\kappa (A)+| A| -| V| }.$$

(7)

The Tutte polynomial is trivial to evaluate along the hyperbola (x − 1)(y − 1) = 1 for any matroid. In the case of graphic matroids, Jaeger et al.¹⁶ showed that the Tutte polynomial is #P-hard to evaluate, except along this hyperbola and when (x,y) equals one of nine special points.

Theorem 1 (Jaeger et al.¹⁶)

The problem of evaluating the Tutte polynomial of a graphic matroid at an algebraic point in the (x,y)-plane is #P-hard except when (x − 1)(y − 1) = 1 or when (x,y) equals one of (1,1), (−1,−1), (0,−1), (−1,0), (i,−i), (−i,i), (j,j²), (j,j²), or (j²,j), where $j=\exp (2\pi {{{\rm{i}}}}/3)$. In each of these exceptional cases the evaluation can be done in polynomial time.

Vertigan¹⁷ extended this result to vector matroids.

Theorem 2 (Vertigan¹⁷)

The problem of evaluating the Tutte polynomial of a vector matroid over a field ${\mathbb{F}}$ at an algebraic point in the (x,y)-plane is #P-hard except when (x − 1)(y − 1) = 1, (x,y) equals (1,1), or when

• (1) $| {\mathbb{F}}| =2$ and (x,y) equals one of (−1,−1), (0,−1), (−1,0), (i,−i), or (−i,i);

• (2) $| {\mathbb{F}}| =3$ and (x,y) equals one of (j,j²) or (j²,j), where $j=\exp (2\pi {{{\rm{i}}}}/3)$; or

• (3) $| {\mathbb{F}}| =4$ and (x,y) equals (−1,−1).

In each of these exceptional cases, except when (x,y) equals (1,1), the evaluation can be done in polynomial time.

Snook¹⁸ showed that when (x,y) equals (1,1) and ${\mathbb{F}}$ is either a finite field of fixed characteristic or a fixed infinite field, then evaluating the Tutte polynomial is #P-hard. It is an open problem to understand the complexity of evaluating the Tutte polynomial at (1,1) over any fixed field.

The Potts model partition function

The Potts model is a statistical physical model described by an integer $q\in {{\mathbb{Z}}}^{+}$ and a graph G = (V,E), with the vertices representing spins and the edges representing interactions between them. A set of edge weights ${\{{\omega }_{e}\}}_{e\in E}$ characterises the interactions and a set of vertex weights ${\{{\upsilon }_{v}\}}_{v\in V}$ characterises the external fields at each spin. A configuration of the model is an assignment σ of each spin to one of q possible states. The Potts model partition function is defined as follows.

Definition 8 (Potts model partition function)

Let $q\in {{\mathbb{Z}}}^{+}$ be an integer and let G = (V,E) be a graph with the weights ${{\Omega }}={\{{\omega }_{e}\}}_{e\in E}$ assigned to its edges and the weights ${{\Upsilon }}={\{{\upsilon }_{v}\}}_{v\in V}$ assigned to its vertices. Then, the q-state Potts model partition function is defined by

$${{{{\rm{Z}}}}}_{{{{\rm{Potts}}}}}(G;q,{{\Omega }},{{\Upsilon }}):=\mathop{\sum}\limits_{\sigma \in {{\mathbb{Z}}}_{q}^{V}}{w}_{G}(\sigma ),$$

(8)

where

$${w}_{G}(\sigma )=\exp \left(\mathop{\sum}\limits_{\{u,v\}\in E}{\omega }_{\{u,v\}}\delta ({\sigma }_{u},{\sigma }_{v})+\mathop{\sum}\limits_{v\in V}{\upsilon }_{v}\delta ({\sigma }_{v})\right).$$

(9)

The Potts model partition function with an external field is equivalent to the zero-field case on an augmented graph $G^{\prime} =(V^{\prime} ,E^{\prime} )$. To construct $G^{\prime}$ from G, for each of the connected components ${\{{C}_{i}\}}_{i = 1}^{\kappa (E)}$ of G add a new vertex u_i and for every vertex v ∈ V(C_i) add an edge e_v = {u_i, v} with the weight υ_v assigned to it. Then, we have the following proposition.

Proposition 3

$${{{{\rm{Z}}}}}_{{{{\rm{Potts}}}}}(G;q,{{\Omega }},{{\Upsilon }})={q}^{-\kappa (E)}{{{{\rm{Z}}}}}_{{{{\rm{Potts}}}}}(G^{\prime} ;q,{{\Omega }}\cup {{\Upsilon }},0).$$

(10)

A similar proposition appears in Welsh’s monograph¹⁹; we prove Proposition 3 in Supplementary Note 1.

It will be convenient to consider the Potts model with weights that are all positive integer multiples of a complex number θ. We shall implement this model on the augmented graph $G^{\prime}$ with all weights equal to θ by replacing each edge with the appropriate number of parallel edges. Let us denote the partition function of this model by ${{{{\rm{Z}}}}}_{{{{\rm{Potts}}}}}(G^{\prime} ;q,\theta )$. Then, we have the following proposition relating the partition function of this model to the Tutte polynomial of the augmented graph $G^{\prime}$.

Proposition 4

$${{{{\rm{Z}}}}}_{{{{\rm{Potts}}}}}(G^{\prime} ;q,\theta )={q}^{\kappa (E^{\prime} )}{({{{{\rm{e}}}}}^{\theta }-1)}^{r(G^{\prime} )}{{{\rm{T}}}}\left(G^{\prime} ;x,y\right),$$

(11)

where $x=\frac{{{{{\rm{e}}}}}^{\theta }+q-1}{{{{{\rm{e}}}}}^{\theta }-1}$ and y = e^θ.

In particular, the q-state Potts model partition function is related to the Tutte polynomial along the hyperbola (x − 1)(y − 1) = q. For a proof of Proposition 4, we refer the reader to Welsh’s monograph¹⁹ [Section 4.4].

The two-state Potts model partition function specialises to the Ising model partition function.

Definition 9 (Ising model partition function)

Let G = (V,E) be a graph with the weights ${{\Omega }}={\{{\omega }_{e}\}}_{e\in E}$ assigned to its edges and the weights ${{\Upsilon }}={\{{\upsilon }_{v}\}}_{v\in V}$ assigned to its vertices. Then, the Ising model partition function is defined by

$${{{{\rm{Z}}}}}_{{{{\rm{Ising}}}}}(G;{{\Omega }},{{\Upsilon }}):=\mathop{\sum}\limits_{\sigma \in {\{-1,+1\}}^{V}}{w}_{G}(\sigma ),$$

(12)

where

$${w}_{G}(\sigma )=\exp \left(\mathop{\sum}\limits_{\{u,v\}\in E}{\omega }_{\{u,v\}}{\sigma }_{u}{\sigma }_{v}+\mathop{\sum}\limits_{v\in V}{\upsilon }_{v}{\sigma }_{v}\right).$$

(13)

Proposition 5

$${{{{\rm{Z}}}}}_{{{{\rm{Potts}}}}}(G;2,{{\Omega }},{{\Upsilon }})={w}_{G}{{{{\rm{Z}}}}}_{{{{\rm{Ising}}}}}\left(G;\frac{{{\Omega }}}{2},\frac{{{\Upsilon }}}{2}\right),$$

(14)

where ${w}_{G}=\exp \left(\frac{1}{2}{\sum }_{e\in E}{\omega }_{e}+\frac{1}{2}{\sum }_{v\in V}{\upsilon }_{v}\right)$.

Proof. The proof follows from some simple algebra.

Instantaneous quantum polynomial time

We shall now briefly introduce the class of commuting quantum circuits, known as IQP time circuits⁵. These circuits exhibit many interesting mathematical properties. In particular, the output probability amplitudes of IQP circuits are proportional to evaluations of the Tutte polynomial of binary matroids². IQP circuits comprise only gates that are diagonal in the Pauli-X basis and are described by an X-program.

Definition 10 (X-program)

An X-program is a pair (P,θ), where $P={({p}_{ij})}_{m\times n}$ is a binary matrix and θ ∈ [−π,π] is a real angle. The matrix P is used to construct a Hamiltonian of m commuting terms acting on n qubits, where each term in the Hamiltonian is a product of Pauli-X operators

$${{{{\rm{H}}}}}_{(P,\theta )}:=-\theta \mathop{\sum }\limits_{i=1}^{m}\mathop{\bigotimes}\limits_{j = 1}^{n}{X}_{j}^{{p}_{ij}}.$$

(15)

Thus, the columns of P correspond to qubits and the rows of P correspond to interactions in the Hamiltonian.

An X-program induces a probability distribution ${{{{\mathcal{P}}}}}_{(P,\theta )}$ known as an IQP distribution.

Definition 11 (${{{{\mathcal{P}}}}}_{(P,\theta )}$)

For an X-program (P,θ) with $P={({p}_{ij})}_{m\times n}$, we define ${{{{\mathcal{P}}}}}_{(P,\theta )}$ to be the probability distribution over binary strings x ∈ {0, 1}ⁿ, given by

$${{{\mathrm{{Pr}}}}}[x]:=| {\psi }_{(P,\theta )}(x){| }^{2},$$

(16)

where

$${\psi }_{(P,\theta )}(x)=\left\langle x\right|\exp \left(-{{{\rm{i}}}}{{{{\rm{H}}}}}_{(P,\theta )}\right)\left|{0}^{n}\right\rangle .$$

(17)

The principal probability amplitude ψ_(P,θ)(0ⁿ) of an IQP distribution is directly related to an evaluation of the Tutte polynomial of the binary matroid whose ground set is the row space of P.

Proposition 6

Let (P,θ) be an X-program with $P={({p}_{ij})}_{m\times n}$. Let $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ be the binary matroid whose ground set ${{{\mathcal{S}}}}$ is the row space of P, then

$${\psi }_{(P,\theta )}(0)={{{{\rm{e}}}}}^{{{{\rm{i}}}}\theta (r(M)-m)}{({{{\rm{i}}}}\sin (\theta ))}^{r(M)}{{{\rm{T}}}}\left(M;x,y\right),$$

(18)

where $x=-{{{\rm{i}}}}\cot (\theta )$ and y = e^2iθ.

A similar result may be obtained for the other probability amplitudes. This can easily be seen when $\theta =\frac{\pi }{2k}$ for $k\in {{\mathbb{Z}}}^{+}$, by firstly letting P∥^kx be the matrix obtained from P by appending k rows identical to x, and then observing that ${\psi }_{(P,\theta )}(x)=-{{{\rm{i}}}}{\psi }_{(P{\parallel }^{k}x,\theta )}({0}^{n})$. For a proof of Proposition 6 and a treatment of the general θ case, we refer the reader to ref. ² [“Discussion”].

We shall consider X-programs that are induced by a weighted graph.

Definition 12 (graph-induced X-program)

For a graph G = (V,E) with the weights ${\left\{{\omega }_{e}\in [-\pi ,\pi ]\right\}}_{e\in E}$ assigned to its edges and the weights ${\left\{{\upsilon }_{v}\in [-\pi ,\pi ]\right\}}_{v\in V}$ assigned to its vertices, we define the X-program induced by G to be an X-program ${{{{\mathcal{X}}}}}_{G}$ such that

$${{{{\rm{H}}}}}_{{{{{\mathcal{X}}}}}_{G}}=-\mathop{\sum}\limits_{\{u,v\}\in E}{\omega }_{\{u,v\}}{X}_{u}{X}_{v}-\mathop{\sum}\limits_{v\in V}{\upsilon }_{v}{X}_{v}.$$

(19)

Any X-program can be efficiently represented by a graph-induced X-program⁵. The principal probability amplitude ${\psi }_{{{{{\mathcal{X}}}}}_{G}}({0}^{n})$ of the IQP distribution generated by a graph-induced X-program is directly related to the Ising model partition function of the graph with imaginary weights.

Proposition 7

Let G = (V,E) be a graph with the weights ${{\Omega}}={\left\{{\omega}_{e}\in [-\pi ,\pi ]\right\}}_{e\in E}$ assigned to its edges and the weights ${{\Upsilon }}={\left\{{\upsilon }_{v}\in [-\pi ,\pi ]\right\}}_{v\in V}$ assigned to its vertices, then

$${\psi }_{{{{{\mathcal{X}}}}}_{G}}\left({0}^{| V| }\right)=\frac{1}{{2}^{| V| }}{{{{\rm{Z}}}}}_{{{{\rm{Ising}}}}}(G;{{{\rm{i}}}}{{\Omega }},{{{\rm{i}}}}{{\Upsilon }}).$$

(20)

Proposition 7 is well known^20,21; we provide a proof in Supplementary Note 2. It will be convenient to consider graph-induced X-programs ${{{{\mathcal{X}}}}}_{G(\theta )}$ with weights that are all positive integer multiples of a real angle θ. As in “Results,” this model can be implemented on the augmented graph $G^{\prime} =(V^{\prime} ,E^{\prime} )$ with all weights equal to θ by replacing each edge with the appropriate number of parallel edges. Let us denote the graph-induced X-program of this model by ${{{{\mathcal{X}}}}}_{G^{\prime} (\theta )}$. Then, we have the following proposition.

Proposition 8

$${\psi }_{{{{{\mathcal{X}}}}}_{G(\theta )}}\left({0}^{| V| }\right)={\psi }_{{{{{\mathcal{X}}}}}_{G^{\prime} (\theta )}}\left({0}^{| V^{\prime} | }\right).$$

(21)

We prove Proposition 8 in Supplementary Note 3. We also have the following proposition relating the principal probability amplitude to the Tutte polynomial of the augmented graph.

Proposition 9

$${\psi }_{{{{{\mathcal{X}}}}}_{G^{\prime} }(\theta )}\left({0}^{| V^{\prime} | }\right)={{{{\rm{e}}}}}^{{{{\rm{i}}}}\theta \left(r(G^{\prime} )-| E^{\prime} | \right)}{({{{\rm{i}}}}\sin (\theta ))}^{r(G^{\prime} )}{{{\rm{T}}}}\left(G^{\prime} ;x,y\right),$$

(22)

where $x=-{{{\rm{i}}}}\cot (\theta )$ and y = e^2iθ.

We prove Proposition 9 in Supplementary Note 4. Note that if we let $M=({{{\mathcal{S}}}},{{{\mathcal{I}}}})$ be the binary matroid whose ground set ${{{\mathcal{S}}}}$ is the column space of the orientated incidence matrix ${A}_{D(G^{\prime} )}$ of $G^{\prime}$ with an arbitrary orientation $D(G^{\prime} )$ assigned to it, then we can use Proposition 6 to obtain Proposition 9.

Quantum computation and the Tutte polynomial

In this section, we show that quantum probability amplitudes may be expressed in terms of the evaluation of a Tutte polynomial. We achieve this by showing that output probability amplitudes of a class of universal quantum circuits are proportional to the principal probability amplitude of some IQP circuit.

It will be convenient to define the following gate set.

Definition 13 (${{{{\mathcal{G}}}}}_{\theta }$)

For a real angle θ ∈ [−π,π], we define ${{{\mathcal{G}}}}(\theta )$ to be the gate set

$${{{{\mathcal{G}}}}}_{\theta }:=\{H,{{{{\rm{e}}}}}^{{{{\rm{i}}}}\theta X},{{{{\rm{e}}}}}^{{{{\rm{i}}}}\theta XX}\},$$

(23)

where H denotes the Hadamard gate.

It is easy to see that the gate set ${{{{\mathcal{G}}}}}_{\frac{\pi }{4}}$ generates the Clifford group and the gate set ${{{{\mathcal{G}}}}}_{\frac{\pi }{8}}$ is universal for quantum computation.

In the IQP model it is easy to implement the gates e^iθX and e^iθXX. So in order to implement the entire gate set ${{{{\mathcal{G}}}}}_{\theta }$, it remains to show that we can implement the Hadamard gate. This can be achieved by the use of postselection when $\theta =\frac{\pi }{4k}$ for $k\in {{\mathbb{Z}}}^{+}$⁶. To apply a Hadamard gate to the target state ${\left|\alpha \right\rangle }_{t}$ consider the following Hadamard gadget. Firstly, introduce an ancilla qubit in the state ${\left|0\right\rangle }_{a}$ and apply the gate ${{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{4}{({\mathbb{I}}-X)}_{t}{({\mathbb{I}}-X)}_{a}}$ to ${\left|\alpha \right\rangle }_{t}{\left|0\right\rangle }_{a}$. Then, measure qubit t in the computational basis and postselect on an outcome of 0. The output state of this gadget is then $H{\left|\alpha \right\rangle }_{a}$.

We shall consider quantum circuits that comprise gates from the set ${{{{\mathcal{G}}}}}_{\frac{\pi }{4k}}$ for an integer $k\in {{\mathbb{Z}}}^{+}$. Let C_k,n,m denote such a circuit that acts on n qubits and comprises m Hadamard gates. Further let ${{{{\mathcal{X}}}}}_{G}({C}_{k,n,m})$ denote the graph-induced X-program that implements the circuit C_k,n,m by replacing each of the m Hadamard gates with the Hadamard gadget. Then, we have the following proposition.

Proposition 10

$$\left\langle {0}^{n}\right|{C}_{k,n,m}\left|{0}^{n}\right\rangle ={\sqrt{2}}^{m}{\psi }_{{{{{\mathcal{X}}}}}_{G}({C}_{k,n,m})}\left({0}^{n+m}\right).$$

(24)

Proof. The proof follows immediately from the application of the Hadamard gadgets.

Any quantum amplitude may therefore be expressed as the evaluation of a Tutte polynomial by Propositions 8–10.

Efficient classical simulation of Clifford circuits

In this section, we show how the correspondence between quantum computation and evaluations of the Tutte polynomial provides an explicit form for Clifford circuit amplitudes in terms of matroid invariants, namely, the bicycle dimension and Brown’s invariant. This gives rise to an efficient classical algorithm for computing the output probability amplitudes of Clifford circuits. We note that it was first observed by Shepherd² that to compute the probability amplitude of a Clifford circuit, it is sufficient to evaluate the Tutte polynomial of a binary matroid at the point (x,y) equals (−i,i), which can be efficiently computed by Vertigan’s algorithm¹⁷. We proceed with some definitions.

Let V be a linear subspace of ${{\mathbb{F}}}_{2}^{n}$. The bicycle dimension and Brown’s invariant are defined as follows.

Definition 14 (bicycle dimension)

The bicycle dimension of V is defined by

$$d(V):=\dim (V\cap {V}^{\perp }).$$

(25)

Definition 15 (Brown’s invariant)

If $| {{{\rm{supp}}}}(x)| \equiv 0\,({{\mathrm{mod}}}\,\,4)$ for all x ∈ V ∩ V^⊥, then Brown’s invariant σ(V) is defined to be the smallest integer such that

$$\mathop{\sum}\limits_{x\in V}{i}^{| {{{\rm{supp}}}}(x)| }={\sqrt{2}}^{d(V)+\dim (V)}{{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{4}\sigma (V)}.$$

(26)

The following theorem of Pendavingh¹⁰ provides an explicit form for the Tutte polynomial of a binary matroid at (−i,i) in terms of the bicycle dimension and Brown’s invariant.

Theorem 11 (Pendavingh¹⁰)

Let V be a linear subspace of ${{\mathbb{F}}}_{2}^{{{{\mathcal{S}}}}}$ and let M(V) be the corresponding binary matroid with ground set ${{{\mathcal{S}}}}$. If $| {{{\rm{supp}}}}(x)| \equiv 0\,({{\mathrm{mod}}}\,\,4)$ for all x ∈ V ∩ V^⊥, then

$${{{\rm{T}}}}(M(V);-{{{\rm{i}}}},{{{\rm{i}}}})={\sqrt{2}}^{d(V)}{{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{4}(2| S| -3r(M)-\sigma (V))}.$$

(27)

Otherwise, T(M(V); −i,i) = 0. Further, T(M(V); −i,i) can be evaluated in polynomial time.

As an immediate consequence of Theorem 11, we obtain an explicit form for Clifford circuit amplitudes in terms the bicycle dimension and Brown’s invariant of the corresponding matroid. Furthermore, we obtain an efficient classical algorithm for computing the output probability amplitudes of Clifford circuits. For similar results of this flavour see refs. ^22,23.

Algorithm overview

We shall now use the correspondence between quantum computation and evaluations of the Tutte polynomial to establish a heuristic algorithm for computing quantum probability amplitudes. To compute a probability amplitude, it is sufficient to compute the Tutte polynomial of a graphic matroid at $x=-{{{\rm{i}}}}\cot \left(\frac{\pi }{4k}\right)$ and $y={{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{2k}}$ for an integer k ≥ 2^2,5,6. Our algorithm will use the deletion–contraction property to recursively compute the Tutte polynomial. At each step in the recursion, the algorithm will compute certain structural properties of the graph in order to attempt to prune the computational tree. Our algorithm can be seen an adaption of the work of Haggard et al.⁷ to special points of the Tutte plane.

We note that our approach differs from tensor network-based methods, which involve the contraction of a graph with tensors assigned to its vertices. These methods have been used to simulate quantum computations while exploiting structural properties of the graph^{4,24,25,26,27}. However, our approach allows us to exploit an alternative class of structural properties. We proceed by describing the key aspects of our algorithm.

To improve the performance of our algorithm, we shall use the following deletion–contraction formula for multigraphs.

Proposition 12

Let G = (V,E) be a multigraph and let e be a multiedge of G with multiplicity ∣e∣. If e is a loop, then

$${{{\rm{T}}}}(G;x,y)={y}^{| e| }{{{\rm{T}}}}(G\setminus \{e\};x,y).$$

(28)

If e is a coloop, then

$${{{\rm{T}}}}(G;x,y)=\left(x+\mathop{\sum }\limits_{i=1}^{| e| -1}{y}^{i}\right){{{\rm{T}}}}(G/\{e\};x,y).$$

(29)

Finally, if e is neither a loop nor a coloop, then

$${{{\rm{T}}}}(G;x,y)={{{\rm{T}}}}(G\setminus \{e\};x,y)+\left(\mathop{\sum }\limits_{i=0}^{| e| -1}{y}^{i}\right){{{\rm{T}}}}(G/\{e\};x,y).$$

(30)

Proposition 12 can easily be proven from the deletion–contraction formula by induction; we omit the proof.

If U is the underlying graph of G, then the number of recursive calls may be bounded by $O\left({2}^{| E(U)| }\right)$. Alternatively, we may bound the number of recursive calls in terms of the number of vertices plus the number of edges s = ∣V(U)∣ + ∣E(U)∣ in the underlying graph. The number of recursive calls R_s is then bounded by R_s ≤ R_s−1 + R_s−2, which is precisely the Fibonacci recurrence. Hence, the number of recursive calls is bounded by $O\left({\phi }^{| V(U)| +| E(U)| }\right)$, where $\phi =\frac{1+\sqrt{5}}{2}$ is the golden ratio²⁸. A careful analysis shows that the number of recursive steps is bounded by $O\left(\tau (U)\cdot | E(U)| \right)$, where τ(U) denotes the number of spanning trees in U¹⁴.

At each step in the recursion, we use the multigraph deletion–contraction formula to remove all multiedges that correspond to either a loop or a coloop in the underlying graph. This process contributes a multiplicative factor to the proceeding evaluation. Note that when G is a graph whose underlying graph is a looped forest, then every edge in the underlying graph is either a loop or a coloop. Hence, we obtain the following formula for the Tutte polynomial of G.

Corollary 13

Let G = (V,E) be a multigraph whose underlying graph U is a looped forest. Further, for each edge e in U, let ∣e∣ denote its multiplicity in G. Then

$${\rm{T}}(G;x,y)=\prod _{\begin{array}{c}e\in E(U)\\ {\rm{loop}}\end{array}}{y}^{|e|}\prod _{\begin{array}{c}e\in E(U)\\ {\rm{coloop}}\end{array}}(x+\mathop{\sum }\limits_{i=1}^{|e|-1}{y}^{i}).$$

(31)

Proof. The proof follows immediately from Proposition 12. ◻

There are a number of techniques that we can use to simplify the graph at each step in the recursion. Firstly, we may remove any isolated vertices, since they do not contribute to the evaluation.

Secondly, when $x=-{{{\rm{i}}}}\cot \left(\frac{\pi }{4k}\right)$ and $y={{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{2k}}$ for an integer $k\in {{\mathbb{Z}}}^{+}$, we may replace each multiedge with a multiedge of equal multiplicity modulo 4k. To account for this, we multiply the proceeding evaluation by a efficiently computable factor. Specifically, we invoke the following proposition.

Proposition 14

Fix $k\in {{\mathbb{Z}}}^{+}$. Let G = (V,E) be a multigraph and let $G^{\prime} =(V^{\prime} ,E^{\prime} )$ be the graph formed from G by taking the multiplicity of each multiedge in G modulo 4k. Then

$${{{\rm{T}}}}(G;x,y)={\left({{{\rm{i}}}}{{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{4k}}\sin \left(\frac{\pi }{4k}\right)\right)}^{\kappa (E)-\kappa (E^{\prime} )}{{{\rm{T}}}}(G^{\prime} ;x,y),$$

(32)

where $x=-{{{\rm{i}}}}\cot \left(\frac{\pi }{4k}\right)$ and $y={{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{2k}}$.

We prove Proposition 14 in Supplementary Note 5.

The Tutte polynomial of a multigraph whose edge multiplicities are all integer multiples of an integer $k\in {{\mathbb{Z}}}^{+}$ may be evaluated at the point $x=-{{{\rm{i}}}}\cot \left(\frac{\pi }{4k}\right)$ and $y={{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{2k}}$ in polynomial time. This can be seen by the following proposition.

Proposition 15

Fix $k\in {{\mathbb{Z}}}^{+}$. Let G = (V,E) be a multigraph whose edge multiplicities are all integer multiples of k. Further let $G^{\prime} =(V^{\prime} ,E^{\prime} )$ be the graph formed from G by taking the multiplicity of each multiedge in G divided by k. Then

$${{{\rm{T}}}}(G;x,y)={\left(\sqrt{2}{{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi (1-k)}{4k}}\sin \left(\frac{\pi }{4k}\right)\right)}^{-r(G)}{{{\rm{T}}}}(G^{\prime} ;-{{{\rm{i}}}},{{{\rm{i}}}}),$$

(33)

where $x=-{{{\rm{i}}}}\cot \left(\frac{\pi }{4k}\right)$ and $y={{{{\rm{e}}}}}^{\frac{{{{\rm{i}}}}\pi }{2k}}$.

We prove Proposition 15 in Supplementary Note 6 and note that this is a special consequence of the k-thickening approach of Jaeger et al.¹⁶. The Tutte polynomial may then be efficiently computed by Vertigan’s algorithm¹⁷; we call such a multigraph a Vertigan graph. We may therefore prune the computational tree whenever the graph is a Vertigan graph with respect to k. Note that this corresponds to quantum circuits comprising gates from the Clifford group.

The Tutte polynomial factorises over components.

Proposition 16

Let G = (V,E) be a graph with connected components $C={\{{C}_{i}\}}_{i = 1}^{k}$, then

$${{{\rm{T}}}}(G;x,y)=\mathop{\prod }\limits_{i=1}^{k}{{{\rm{T}}}}({C}_{i};x,y).$$

(34)

Proposition 16 can easily be proven from the deletion–contraction formula; we omit the proof. At each step in the deletion–contraction recursion, if the graph is disconnected, then we may use this property to prune the computational tree and hence improve performance.

An identical result holds for biconnected components.

Proposition 17 (Tutte²⁹)

Let G = (V,E) be a graph with biconnected components $B={\{{B}_{i}\}}_{i = 1}^{k}$, then

$${{{\rm{T}}}}(G;x,y)=\mathop{\prod }\limits_{i=1}^{k}{{{\rm{T}}}}({B}_{i};x,y).$$

(35)

Proposition 17 can easily be proven from the deletion–contraction formula. For a proof, we refer the reader to ref. ²⁹ [“Discussion”]. Similarly to the connected component case, we may use this property to prune the computational tree and improve performance. Note that the biconnected components of a graph may be listed in time linear in the number of edges via depth-first search³⁰.

The Tutte polynomial of a multigraph whose underlying graph is a cycle may be computed in polynomial time by invoking the following proposition.

Proposition 18 (Haggard et al.⁷)

Let G = (V,E) be a multigraph whose underlying graph U is an n-cycle with edges indexed by the positive integers. Further, for each edge e in U, let ∣e∣ denote its multiplicity in G. Then

$$\begin{array}{lll}{{{\rm{T}}}}(G;x,y)&=&\mathop{\sum }\limits_{k=1}^{n-2}\left(\mathop{\prod }\limits_{j=k+1}^{n}{y}_{x}\left(| {e}_{j}| \right)\mathop{\prod }\limits_{j=1}^{k-1}{y}_{1}\left(| {e}_{j}| \right)\right)\\ &&+\,{y}_{x}\left(| {e}_{n}| +| {e}_{n-1}| \right)\mathop{\prod }\limits_{j=1}^{n-2}{y}_{1}\left(| {e}_{j}| \right),\end{array}$$

(36)

where ${y}_{x}(j):=x+\mathop{\sum }\nolimits_{i = 1}^{j-1}{y}^{i}$.

Proposition 18 can easily be proven from the deletion–contraction formula. For a proof, we refer the reader to ref. ⁷ [Theorem 4]. We may use this proposition to prune the computational tree whenever the underlying graph is a cycle.

The Tutte polynomial of a planar graph along the hyperbola (x − 1)(y − 1) = 2 may be evaluated in polynomial time via the Fisher–Kasteleyn–Temperley algorithm^31,32,33. We may therefore use this algorithm to prune the computational tree whenever the underlying graph is planar. Note that we may test whether a graph is planar in time linear in the number of vertices³⁴.

The performance of our algorithm depends on the heuristic used to select edges. We shall consider six edge-selection heuristics: vertex order, minimum degree, maximum degree, minimum degree sum, maximum degree sum, and non-Vertigan. These edge-selection heuristics were first studied by Pearce et al.⁸, with the exception of non-Vertigan, which is specific to our algorithm.

Vertex order: The vertices of the graph are assigned an ordering. A multiedge is selected from those incident to the lowest vertex in the ordering and whose other endpoint is also the lowest vertex of any incident in the ordering. For contractions, the vertex inherits the lowest of the positions in the ordering.

Minimum degree: A multiedge is selected from those incident to a vertex with minimal degree in the underlying graph.

Maximum degree: A multiedge is selected from those incident to a vertex with maximal degree in the underlying graph.

Minimum degree sum: A multiedge is selected from those whose sum of degrees of its endpoints is minimal in the underlying graph.

Maximum degree sum: A multiedge is selected from those whose sum of degrees of its endpoints is maximal in the underlying graph.

Non-Vertigan: A multiedge is selected from those whose multiplicity is not an integer multiple of k; we call such a multiedge non-Vertigan. Using this edge-selection heuristic, the number of recursive calls may be bounded by $O\left({2}^{\nu (G)}\right)$, where ν(G) denotes the number of non-Vertigan multiedges in G. This is due to the fact that both the deletion and contraction operation reduce the number of non-Vertigan multiedges by at least one. We note that this is similar to the sum-over-Cliffords approach studied in refs. ^35,36,37.

There are many other methods that may improve the performance of our algorithm, which we do not study. We shall proceed by discussing some of these.

Isomorphism testing: During the computation the graphs encountered and the evaluation of their Tutte polynomial is stored. At each recursive step, we test whether the graph is isomorphic to one already encountered, and if so, we use the evaluation of the isomorphic graph instead. Haggard et al.⁷ showed that isomorphism testing can lead to an improvement in the performance of computing Tutte polynomials. Note that this may not be as effective when the input is a multigraph.

Almost planar: At each step in the recursion, we may test whether the graph is close to being planar, and if so, select edges in such a way that the deletion and contraction operations give rise to a planar graph. For example, if the graph is apex, that is, it can be made planar by the removal of a single vertex, then we may select a multiedge incident to such a vertex. Similarly, if the underlying graph is edge apex or contraction apex³⁸, then we may select a multiedge such that the deletion or the contraction operation gives rise to a planar graph.

k-connected components: Similarly to the connected and biconnected component case, we may compute the Tutte polynomial in terms of its k-connected components^39,40.

Experimental results

In this section, we present some experimental results comparing the performance of the edge-selection heuristics on two classes of random quantum circuits. Our experiments were performed using SageMath 9.0⁴¹. The source code and experimental data are available at ref. ⁴².

The first class we consider corresponds to random instances of IQP circuits induced by dense graphs. Specifically, an instance is an IQP circuit induced by a complete graph with edge weights chosen uniformly at random from the set $\{\frac{m\pi }{8}| m\in {\mathbb{Z}}/8{\mathbb{Z}}\}$. This class of IQP circuits is precisely that was studied in ref. ⁴³, where it is conjectured that approximating the corresponding amplitudes up to a multiplicative error is #P-hard on average.

The second class we consider corresponds to random instances of IQP circuits induced by sparse graphs. Specifically, an instance is an IQP circuit induced by a random graph where each of the possible edges is included independently with probability 1/2 and with edge weights chosen uniformly at random from the set $\{\frac{m\pi }{8}| m\in {\mathbb{Z}}/8{\mathbb{Z}}\}$.

We run our algorithm using each of the edge-selection heuristics to compute the principal probability amplitude of 64 random instances of both the dense and sparse class on 12 vertices. The performance of each edge-selection heuristic is measured by counting the number of leaves in the computational tree. Our experimental data are presented in Supplementary Note 7. We find that the non-Vertigan edge-selection heuristic performs particularly well for the dense class and the maximum degree sum edge-selection heuristic performs particularly well for the sparse class.

Discussion

We established a classical heuristic algorithm for exactly computing quantum probability amplitudes. Our algorithm is based on mapping output probability amplitudes of quantum circuits to evaluations of the Tutte polynomial of graphic matroids. The algorithm evaluates the Tutte polynomial recursively using the deletion–contraction property while attempting to exploit structural properties of the matroid. We considered several edge-selection heuristics and presented experimental results comparing their performance on two classes of random quantum circuits. Further, we obtained an explicit form for Clifford circuit amplitudes in terms of matroid invariants and an alternative efficient classical algorithm for computing the output probability amplitudes of Clifford circuits.