The complexity of approximating complex-valued Ising and Tutte partition functions

We study the complexity of approximately evaluating the Ising and Tutte partition functions with complex parameters. Our results are partly motivated by the study of the quantum complexity classes BQP and IQP. Recent results show how to encode quantum computations as evaluations of classical partition functions. These results rely on interesting and deep results about quantum computation in order to obtain hardness results about the difficulty of (classically) evaluating the partition functions for certain fixed parameters. The motivation for this paper is to study more comprehensively the complexity of (classically) approximating the Ising and Tutte partition functions with complex parameters. Partition functions are combinatorial in nature and quantifying their approximation complexity does not require a detailed understanding of quantum computation. Using combinatorial arguments, we give the first full classification of the complexity of multiplicatively approximating the norm and additively approximating the argument of the Ising partition function for complex edge interactions (as well as of approximating the partition function according to a natural complex metric). We also study the norm approximation problem in the presence of external fields, for which we give a complete dichotomy when the parameters are roots of unity. Previous results were known just for a few such points, and we strengthen these results from BQP-hardness to #P-hardness. Moreover, we show that computing the sign of the Tutte polynomial is #P-hard at certain points related to the simulation of BQP. Using our classifications, we then revisit the connections to quantum computation, drawing conclusions that are a little different from (and incomparable to) ones in the quantum literature, but along similar lines.

Abstract. We study the complexity of approximately evaluating the Ising and Tutte partition functions with complex parameters. Our results are partly motivated by the study of the quantum complexity classes BQP and IQP. Recent results show how to encode quantum computations as evaluations of classical partition functions. These results rely on interesting and deep results about quantum computation in order to obtain hardness results about the difficulty of (classically) evaluating the partition functions for certain fixed parameters. The motivation for this paper is to study more comprehensively the complexity of (classically) approximating the Ising and Tutte partition functions with complex parameters. Partition functions are combinatorial in nature and quantifying their approximation complexity does not require a detailed understanding of quantum computation. Using combinatorial arguments, we give the first full classification of the complexity of multiplicatively approximating the norm and additively approximating the argument of the Ising partition function for complex edge interactions (as well as of approximating the partition function according to a natural complex metric). We also study the norm approximation problem in the presence of external fields, for which we give a complete dichotomy when the parameters are roots of unity. Previous results were known just for a few such points, and we strengthen these results from BQP-hardness to #P-hardness. Moreover, we show that computing the sign of the Tutte polynomial is #P-hard at certain points related to the simulation of BQP. Using our classifications, we then revisit the connections to quantum computation, drawing conclusions that are a little different from (and incomparable to) ones in the

Introduction
We study the Ising and Tutte partition functions, which are wellknown partition functions arising in combinatorics and statistical physics (see, for example, Sokal (2005)). Early works which studied the complexity of (exactly) evaluating these partition functions Jaeger et al. (1990) considered both real and complex parameters. Applications in statistical mechanics actually require consideration of complex numbers because the possible points of physical phase transitions occur exactly at real limit points of complex zeroes of these partition functions (see Sokal's explanation in Section 5 "Complex Zeros of Z G : Why should we care?" Sokal (2005)). However, given the difficulty of completely resolving the complexity of the approximation problem, most works which comprehensively studied the complexity of approximately evaluating these partition functions Goldberg & Jerrum (2008, 2014; Jerrum & Sinclair (1993) restricted attention to real parameters. A notable counter-example is the paper of Bordewich et al. (2005) which studied normalised additive approximations for #P functions including these partition functions. Bordewich et al. were motivated by a result of Freedman et al. (2003) showing that an approximate evaluation of the Jones polynomial associated with a particular complex parameter (a 5th root of unity) can be used to simulate the quantum part of any algorithm in the quantum complexity class BQP, which is the class of decision problems solvable by a quantum computer in polynomial time with bounded error. The relevance of this result to the partition functions that we study follows from a result of Thistlethwaite (1987), showing that the Jones polynomial is essentially a specialisation of the Tutte partition function.
Recently, there have been several papers showing how to encode quantum computations as evaluations of partition functions. These results rely on interesting and deep results about quantum computation to obtain hardness results about the difficulty of (classically) evaluating Ising and Tutte partition functions. For example, Kuperberg (2015) used three results in quantum computation (a density theorem from Freedman et al. (2002), the Solovay-Kitaev theorem (see Nielsen & Chuang (2004)), and PostBQP=PP Aaronson (2005)) to demonstrate the #P-hardness of a certain kind of approximation of the Jones polynomial. His theorem is repeated later as Theorem 5.1, where it is discussed in more detail. He also derived related results about multiplicative approximations of the Tutte polynomial for certain real parameters.
IQP stands for "Instantaneous Quantum Polynomial time". It is characterised by a class of quantum circuits introduced by Shepherd & Bremner (2009). Fujii & Morimae (2013 showed how to encode IQP circuits as instances of the Ising model. Thus, they were able to use a quantum complexity result of (Bremner et al. 2011, Corollary 3.3) (showing that weakly simulating IQP with multiplicative error implies that the polynomial hierarchy collapses to the third level) to obtain a result about the approximation of the Ising model -namely that an FPRAS for the Ising model with parameter y = exp(iπ/8) would similarly entail collapse of the polynomial hierarchy. (As they mention, a similar result applies for other parameters that are universal for IQP.) This result is further discussed in Section 4.1. Other examples include (De las Cuevas et al. 2011, Result 2), (Iblisdir et al. 2014, Theorem 6.1), and (Matsuo et al. 2014, Theorems 2 and 3) which give BQPhardness of certain Ising model approximations, enabling the conclusion that certain efficient algorithms for approximating these partition function up to additive error are unlikely to exist. Iblisdir et al. (2014) point out that some instances that they prove hard do have multiplicative approximations, due to Jerrum & Sinclair (1993), emphasising the difference between additive and multiplicative approximation. (Matsuo et al. 2014, Theorem 4) also relate the simulation of IQP circuits to Ising model approximations with real parameters.
The motivation for our paper is to study more comprehen-sively the complexity of approximating the Ising partition function at complex parameters, and also to go the other way around, working from the combinatorial model to quantum computation. Partition functions are combinatorial in nature and classifying the difficulty of approximating these partition functions should not require a detailed understanding of quantum mechanics or quantum computation. Hence, we undertake a detailed classification of the complexity of the partition function problems, using combinatorial methods. We focus mainly on the Ising model since this model is particularly relevant in statistical physics (Section 3). This model is also connected to IQP (as explained in Section 4.1). We also consider the more general Tutte polynomial at any point (x, y) where x = −t and y = −t −1 for a root of unity t (this is connected to BQP, as will be explained in Section 5).
Our main result for the Ising model (Theorem 1.2) is a classification of the complexity of approximating the partition function with complex edge interactions. This result is illustrated in Section 1.3. As the figure shows, there are very few parameters (edge interactions) in the complex plane for which the approximation problem is tractable. For most edge interactions, it is extremely intractable (#P-hard to approximate the norm within any constant factor and to approximate the argument within ±π/3). Theorem 1.3 extends these results to a more relaxed setting in which approximation algorithms are unconstrained (allowed to output any rational number) if the correct output is zero. We emphasise that the goal of our work is to classify the difficulty of the problem for all fixed parameters in the complex plane. The proofs of our theorems are elementary and combinatorial. The main idea (see Lemma 3.6) is an extension of a bisection technique of Goldberg & Jerrum (2014) showing how to use an approximation for the norm of a function to get very close to a zero of the function. Our result for the Tutte polynomial (Theorem 1.7) is also proved using bisection. It shows that, for any relevant parameters, it is #P-hard to determine whether the sign of the polynomial is non-negative or non-positive (with an arbitrary answer being allowed when it is zero).
Using our classifications, we then revisit the connections to quantum computation, drawing conclusions that are a little different from (and incomparable to) the ones in the papers mentioned earlier, but along similar lines, as we now explain. Theorem 1.4 shows that strong simulation of IQP within any constant factor is #P-hard, even for the restricted class of circuits considered by Bremner et al. (2011). Our result is incomparable to their hardness result (Bremner et al. 2011, Corollary 3.3). Both results show hardness of multiplicative approximation. However, their result is for weak simulation (sampling from the output distribution of the circuit) whereas ours is for strong simulation (estimating the probability of a given output). In general, hardness results about weak simulation are more desirable, however multiplicative approximation is less appropriate for weak simulation, where total variation distance is more important. Also, our results (unlike those of Bremner et al. (2011)) are not sensitive to the behaviour of the algorithms when the correct value is zero. Moreover, our complexity assumption (that FP = #P) is implied by and therefore milder than theirs (that the polynomial hierarchy does not collapse to the third level). These results are discussed further in Section 4.1.
It seems that a result similar to our IQP result could also be obtained via Boson sampling Aaronson & Arkhipov (2013). In particular, (Aaronson & Arkhipov 2013, Theorem 4.3) have used a bisection technique similar to the one of Goldberg & Jerrum (2014) to show that approximating the square of the permanent of a realvalued input matrix within a constant factor is #P-hard. Any such input (Aaronson & Arkhipov 2013, Lemma 4.4) can be turned into a unitary matrix which can be viewed as a "Boson Sampling" input. The output of the Boson sampling problem is essentially the square of the permanent of the matrix (so is hard to approximate). Furthermore, the Boson sampling problem can be simulated by BQP circuits and adaptive IQP circuits (in the strong sense). Thus, while it is interesting to see that our Ising-model results have IQP applications, the important point concerning our result is the comprehensive classification of the Ising complexity, rather than the particular quantum applications.
As we explain in Section 1.5.2, classical simulation of the complexity class BQP is related to (but not directly a consequence of) determining the sign of the Tutte polynomial at a certain point (−t, −t −1 ). Theorem 1.7 shows that this problem is #P-hard (even when the algorithm is not required to handle the case in which the output is zero), answering a question raised by Bordewich et al. (2005). This is related to (but incomparable to) a result (Theorem 5.1) of Kuperberg (2015). These results are discussed further in Section 5.
Finally, we study Ising models with external fields. (De las Cuevas et al. 2011, Result 2) showed that with edge interaction i and external field e iπ/4 an additive approximation of the partition function is BQP-hard. Motivated by such connections, we focus on the problem of (multiplicatively) approximating the norm of the partition function when both the interaction parameter and the external field are roots of unity. We extend our hardness results to show that, for most such parameters, including the one studied by De las Cuevas et al., the approximation problem is #P-hard (for an exact statement, see Theorem 1.9). For the remaining parameters, the partition function can be evaluated exactly in polynomial time, and thus we get a complete dichotomy (Theorem 1.9). This extension relies on some lower bounds from transcendental number theory, which allow us to convert additive distances into multiplicative ones. The lower bound results are given in Section 6.1 and our hardness results are in Section 6.2.
As we have already mentioned, there are many papers encoding quantum simulations as Ising models, including especially the result of Fujii & Morimae (2013). We could use this encoding (along with our Theorem 1.3) to derive our quantum application (Theorem 1.4). In order to make the paper self-contained, and to make it accessible to readers from outside the area of quantum computation we instead give our own, more combinatorial, presentation of how to encode IQP circuits as Ising instances. This is given in Section 4.1.
1.1. The Ising model. The main partition function that we study is the partition function of the Ising model. Let y (called the edge interaction) and λ (called the external field ) be two parameters. The partition function is defined for a (multi)graph where m(σ) is the number of monochromatic edges under σ (that is, the number of edges (u, v) with σ(u) = σ(v)) and n 1 (σ) is the number of vertices v with σ(v) = 1. We write Z Ising (G; y) to denote Z Ising (G; y, 1).
We will consider complex parameters y and λ from the set Q of algebraic numbers. Thus, the real and imaginary parts of y and λ will be algebraic. We use arg(z) to denote the arg of a complex number z. For fixed y and λ, we study several computational problems. The first of them is approximating the norm of Z Ising (G; y, λ) within a factor K > 1.
Instance A (multi)graph G.
Output A rational number N such that N /K ≤ |Z Ising (G; y, λ)| ≤ K N .
We also consider the problem of approximating the argument of the partition function within an additive distance of ρ ∈ (0, 2π). Here we have to treat the zero case exceptionally since the argument is undefined.

Approximating Complex Numbers.
It makes sense that we approximate the norm of a complex number relatively, whereas we approximate the argument additively. This is natural because multiplying complex numbers multiplies norms and adds arguments, so it preserves the usual property that if you can approximate two numbers, you can approximate the product.
Other notions of approximation have been proposed. Most notably, Ziv (1982) has proposed that the distance between two complex numbers y and y should be measured as where d(0, 0) = 0. We also study the following approximation problem.
Instance A (multi)graph G and a positive integer R, in unary.
Output If |Z Ising (G; y, λ)| = 0 then the algorithm should output 0. Otherwise, it should output a complex number y such that d(y, Z Ising (G; y, λ)) ≤ 1 R .
As with the other problems, we use the notation Complex-Apx-Ising(y) for Complex-Apx-Ising(y, 1). We have specified the error R as an input of the problem, rather than as a parameter in order to emphasise the suitability of Complex-Apx-Ising(y, λ) as an appropriate notion of approximation for the Ising partition function when y is complex. The number R is expressed in unary so a polynomial time algorithm for Complex-Apx-Ising(y, λ) would give a so-called "fully polynomial time approximation scheme" for the norm of the partition function. For partition functions, it is well-known that approximating the norm within a factor that is an inverse polynomial in a unary input R is equivalent in difficultly to approximating the norm with any specific factor K > 1. We will return to this point later in Lemma 3.2.
1.3. Main results for the Ising model. The following theorem gives our main complexity results about the Ising model. The five white points correspond to the easy evaluations described in Item 1. The green line segment (1, ∞) corresponds to a region where approximation is in RP-See Item 2. The blue line segment (−∞, −1) corresponds to a region where approximation is equivalent to approximately counting perfect matchings. See Item 4. The red points on the axes (the imaginary axis and the segment (−1, 0)) and on the unit circle correspond to regions where approximation is #P-hard. See Items 5, 6, and 7. Elsewhere the points are coloured grey, and approximation is known to be NPhard (Items 3, 9 and 10) and sometimes to be #P-hard (Item 8, not pictured).
These results classify the problem of approximating the partition function over the entire complex plane. For every value of the parameter y, we either show that the problem is easy, in the sense that both the norm and the arg of the partition function can be well-approximated (and so can the conglomerate problem using the Ziv distance), or we show that approximating at least one these is hard (and so is the conglomerate problem using the Ziv distance).
The results for approximation of the norm are illustrated in Section 1.3.
(iv) If y < −1 is a real number then Factor-K-Norm-Ising(y) is equivalent in complexity to the problem of approximately counting perfect matchings in graphs and Complex-Apx-Ising(y) is as hard. However, Distance-(π/3)-Arg-Ising(y) is in FP.
1.4. Relaxed versions of the problems. A polynomial-time algorithm for any of the problems that we have defined is required to output 0 if it is given an input G such that Z Ising (G; y, λ) = 0. Theorem 1.2 gives hardness results for these problems. The hardness is not due to special difficulties which arise when the value of the partition function is zero. In order to demonstrate this point, (and in order to make certain reductions easier later on), we also consider the following, more relaxed versions of the problems, where the output is unconstrained if the value of the partition function is zero. As before, the parameter K is greater than 1 and the parameter ρ is in (0, 2π).
Instance A (multi)graph G.
Output If |Z Ising (G; y, λ)| = 0 then the algorithm may output any rational number. Otherwise, it must output a rational number N such that N /K ≤ |Z Ising (G; y, λ)| ≤ K N .
Instance A (multi)graph G.
Output If Z Ising (G; y, λ) = 0, then the algorithm may output any rational number. Otherwise, it must output a rational number A such that | A − arg(Z Ising (G; y, λ))| ≤ ρ.
Instance A (multi)graph G and a positive integer R, in unary.
Output If |Z Ising (G; y, λ)| = 0 then the algorithm may output any complex number. Otherwise, it must output a complex number z such that d(z, Z Ising (G; y, λ)) ≤ 1 R .
As in the un-relaxed versions of the problems, we drop the parameter "λ" from the problem name when it is 1. We give the following generalisation of Theorem 1.2. Theorem 1.3. All of the results in Theorem 1.2 extend to the relaxed case. That is, the results are still true with Factor-K-Norm-Ising(y), Distance-(π/3)-Arg-Ising(y) and Complex-Apx-Ising(y) replaced by Factor-K-Nonzero-Norm-Ising(y), Distance-(π/3)-Nonzero-Arg-Ising(y) and Complex-Apx-Nonzero-Ising(y), respectively.
1.5. Applications to quantum simulation.
1.5.1. IQP. IQP is characterised by a restricted class of quantum circuits Shepherd & Bremner (2009). We will give a formal definition in Section 4.1. There we will also discuss related work by Fujii & Morimae (2013), Bremner et al. (2011) and Jozsa & Van den Nest (2014). Here we give an informal description that enables us to state our theorem. Bremner et al. showed a hardness of a certain kind of "weak simulation" of a restricted class of circuits called IQP 1,2 (θ) circuits (see Definition 4.8). The qubits of the circuit travel along "lines" which go into (and out of) quantum gates. The output of such a circuit C is a random variable Y (over the qubits that get measured in the output). Given as input the string of all zero qubits and an output string y ∈ {0, 1} |I| on a set I -the set of qubits that are measured in the output, Pr C;I denotes the probability that Y = y. Strong simulation is the problem of (approximately) computing this probability. We consider the following problem where K > 1 is an error parameter.
Output A rational number p such that p/K ≤ Pr C;I (Y = y) ≤ Kp.
Our main result regarding this application is the following.
1.5.2. Connections between the Sign of the Tutte Polynomial and BQP. The partition function Z Ising (G; y, λ) is equivalent to a specialisation of the Tutte polynomial, which is a graph polynomial with two parameters, x and y, defined as follows, where n = |V (G)| and κ(A) is the number of connected components in the subgraph (V (G), A). If the quantity q = (x − 1)(y − 1) is a positive integer, then the Tutte polynomial with parameters x and y is closely related to the partition function of the Potts model, which includes the Ising model as the special case q = 2. In particular, when q = 2, Bordewich et al. (2005) raised the question "of determining whether the Tutte polynomial is greater than or equal to, or less than zero at a given point." As we will see, this question is relevant to the quantum complexity class BQP. We consider the following problems.
Name Sign-Real-Tutte(x, y) Instance A (multi)graph G.
Output Determine whether the sign of the real part of T (G; x, y) is positive, negative, or 0.
BQP is the class of decision problems solvable by a quantum computer in polynomial time with bounded error. The theorem (Bordewich et al. 2005, Theorem 6.1) shows that all of the problems in BQP can also be solved classically in polynomial time using an oracle that returns the sign of the real part of the Jones polynomial of a link, evaluated at the point t = exp(2πi/5). Thistlethwaite (1987) (see (Jaeger et al. 1990, (6.1))), showed that this problem is, in turn, related to the problem of evaluating the Tutte polynomial T (G; −t, −t −1 ), for a planar graph G. This inspired the question of Bordewich et al. about the complexity of determining the sign of the Tutte polynomial, particularly for the point (x, y) = (−t, −t −1 ). We show that problem is hard for values of t including the relevant value t = exp(2πi/5). Note that our result does not have direct implications for the simulation of BQP because we do not deal with planarity (though it does answer the question of Bordewich et al.). We give the details in Section 5, where we also discuss a related result of Kuperberg (2015). Our theorem is as follows.
1.6. Results about Ising models with fields. Our results in Section 1.3 are about the complexity of evaluating the Ising partition function in the absence of an external field (when λ = 1). This is appropriate for the application to IQP. Ising models with external fields are important for their own sake. Moreover, (De las Cuevas et al. 2011, Result 2) showed that with edge interaction i and external field e iπ/4 an additive approximation of the partition function is BQP-hard. Motivated by such quantum connections, we give the following extension.
Theorem 1.9. Let K > 1. Let y and z be two roots of unity. Then the following holds: (i) If y = ±i and z ∈ {1, −1, i, −i}, or y = ±1, then Z Ising (−; y, z) can be computed exactly in polynomial time.

Preliminaries
2.1. Facts about Approximating Complex Numbers. We will use the following technical lemma concerning Ziv's distance measure from Section 1.2.
so |z| is a suitable output to Factor-K-Norm-Ising(y, λ) or Factor-K-Nonzero-Norm-Ising(y, λ) with input G. On the other hand, if |Z Ising (G; y, λ)| = 0 then |z| is still suitable in both cases.
2.2. The multivariate Tutte polynomial. We will require the random cluster formulation of the multivariate Tutte polynomial. Given a (multi) graph G with edge weights γ : E(G) → Q and q ∈ Q, this is defined as where γ e is a shorthand for γ(e) for an edge e ∈ E(G).
To apply a technique from Goldberg & Jerrum (2014) we will require a multivariate version of the problem Factor-K-Nonzero-Norm-Ising(y, λ). We could do this for general q, but we will only use the following version, which is restricted to q = 2 and has two complex parameters, γ 1 and γ 2 .
Output If |Z Tutte (G; 2, γ)| = 0 then the algorithm may output any rational number. Otherwise, it should output a rational Suppose that s and t are two distinguished vertices of G. Let Z st (G; q, γ) be the contribution to Z Tutte (G; q, γ) from subgraphs A where s and t are in the same component of (V (G), A), that is, Similarly let Z s|t denote the contribution to Z Tutte (G; q, γ) from configurations A in which s and t are in different components.
2.3. Implementing new edge weights, series compositions, and parallel compositions . Our treatment of implementations, series compositions and parallel compositions is completely standard and is taken from (Goldberg & Jerrum 2012, Section 2.1). The reader who is familiar with this material can skip this section (which is included here for completeness).
Fix W ⊆ Q and q ∈ Q. Let w * ∈ Q be a weight (which may not be in W ) which we want to "implement". Suppose that there is a graph Υ, with distinguished vertices s and t and a weight function In this case, we say that Υ and γ implement w * (or even that W implements w * ). The purpose of "implementing" edge weights is this. Let G be a graph with weight function γ. Let f be some edge of G with weight γ f = w * . Suppose that W implements w * . Let Υ be a graph with distinguished vertices s and t with a weight function γ : E(Υ) → W satisfying (2.5). Construct the weighted graph G by replacing edge f with a copy of Υ (identify s with either endpoint of f (it doesn't matter which one) and identify t with the other endpoint of f and remove edge f ). Let the weight function γ of G inherit weights from γ and γ (so γ e =γ e if e ∈ E(Υ) and γ e = γ e otherwise). Then the definition of the multivariate Tutte polynomial gives So, as long as q = 0 and Z s|t (Υ; q, γ) is easy to evaluate, evaluating the multivariate Tutte polynomial of G with weight function γ is essentially the same as evaluating the multivariate Tutte polynomial of G with weight function γ.
Since the norm of the product of two complex numbers is the product of the norms, this reduces computing (or relatively approximating) the norm with weight function γ to the problem of computing (or relatively approximating) the norm with weight function γ . Also, since the argument of the product of two complex numbers is the sum of the arguments of the numbers, this reduces computing (or additively approximating) the argument with weight function γ to the problem of computing (or additively approximating) the argument with weight function γ .
Two especially useful implementations are series and parallel compositions. Parallel composition is the case in which Υ consists of two parallel edges e 1 and e 2 with endpoints s and t andγ e 1 = w 1 andγ e 2 = w 2 . It is easily checked from Equation (2.5) that w * = (1 + w 1 )(1 + w 2 ) − 1. Also, the extra factor in Equation (2.6) cancels, so in this case Z Tutte (G ; q, γ ) = Z Tutte (G; q, γ).
Series composition is the case in which Υ is a length-2 path from s to t consisting of edges e 1 and e 2 withγ e 1 = w 1 andγ e 2 = w 2 . It is easily checked from Equation (2.5) that w * = w 1 w 2 /(q + w 1 + w 2 ). Also, the extra factor in Equation (2.6) is q + w 1 + w 2 , so in this case Z Tutte (G ; q, γ ) = (q + w 1 + w 2 )Z Tutte (G; q, γ). It is helpful to note that w * satisfies We say that there is a "shift" from (q, α) to (q, α ) if there is an implementation of α consisting of some Υ and w : E(Υ) → W where W is the singleton set W = {α}. Taking y = α + 1 and y = α + 1 and defining x and x by q = (x − 1)(y − 1) = (x − 1)(y − 1) we equivalently refer to this as a shift from (x, y) to (x , y ). It is an easy, but important observation that shifts may be composed to obtain new shifts. So, if we have shifts from (x, y) to (x , y ) and from (x , y ) to (x , y ), then we also have a shift from (x, y) to (x , y ).
The k-thickening of Jaeger et al. (1990) is the parallel composition of k edges of weight α. It implements α = (1 + α) k − 1 and is a shift from (x, y) to (x , y ) where y = y k (and x is given by (x − 1)(y − 1) = q). Similarly, the k-stretch is the series composition of k edges of weight α. It implements an α satisfying It is a shift from (x, y) to (x , y ) where x = x k . (In the classical bivariate (x, y) parameterisation, there is effectively one edge weight, so the stretching or thickening is applied uniformly to every edge of the graph.) Thus, we have the following observation.
Observation 2.7. The k-thickening operation gives the following polynomial-time reductions.

Hardness results for the Ising model
In this section we prove Theorems Theorem 1.2 and Theorem 1.3.
3.1. Real weights. First we gather some known results regarding approximating the partition function Z Ising (G; y) of the Ising model when y is an algebraic real number. If y ∈ {−1, 0, 1}, then computing Z Ising (G; y) is trivial from the definition (1.1). A classical result by Jerrum and Sinclair Jerrum & Sinclair (1993) settles the complexity of approximating Z Ising (G; y) when y > 0. They show that there is a "fully polynomial randomised approximation scheme" (FPRAS) when y > 1 and that it is NP-hard to approximate the partition function when 0 < y < 1.
The negative case appears to be more complicated. Goldberg and Jerrum Goldberg & Jerrum (2008) showed that if −1 < y < 0, it is also NP-hard to approximate Z Ising (G; y), but if y < −1, the problem is equivalent to approximating the number of perfect matchings in a graph and it is not known whether there is an FPRAS. Technically, neither Jerrum and Sinclair nor Goldberg and Jerrum worked over the algebraic numbers. In order to avoid issues of real arithmetic, Jerrum and Sinclair used a computational model in which real arithmetic is performed with perfect accuracy, and Goldberg and Jerrum restricted attention to rationals. However, the operations in those papers are easily implemented over the algebraic real numbers. Using our notation, these results are summarised as follows. Suppose y ∈ Q and K > 1. Then Factor-K-Norm-Ising(y) • is equivalent in difficulty to approximately counting perfect matchings if y < −1.
Technically, the results in Goldberg & Jerrum (2008); Jerrum & Sinclair (1993) were not about the problem Factor-K-Norm-Ising(y) with fixed K. Instead, the accuracy parameter was viewed as part of the input as in the following problem.
Instance A (multi)graph G and a positive integer R, in unary.
Nevertheless, the hardness results in Lemma 3.1 follow easily from those papers using the following standard powering lemma.
Proof. The reduction from Factor-K-Norm-Ising(y, λ) to FPRAS-Norm-Ising(y, λ) is straightforward: Given an input G to Factor-K-Norm-Ising(y, λ), choose R so that K ≥ R/(R − 1) and run an algorithm for FPRAS-Norm-Ising(y, λ) with inputs G and R, returning the result. The other direction is almost as easy. Given an input (G, R) to FPRAS-Norm-Ising(y, λ), choose an integer k sufficiently large (which does not depend on the size of G) so that (1−1/R) k ≤ 1/K and (1+1/R) k ≥ K. Then form G k by taking k disjoint copies of G. Run an algorithm for Factor-K-Norm-Ising(y, λ) with input Note that the NP-hardness result for 0 < y < 1 in Lemma 3.1 is essentially best possible in the sense that the problem is not much harder than NP. As Goldberg & Jerrum (2008) observed, the problem can be solved in randomised polynomial time using an oracle for an NP predicate by applying the bisection technique of Valiant & Vazirani (1986). The situation is different for y < 0. (Goldberg & Jerrum 2014, Theorem 1, Region G) showed that it is #P-hard to determine the sign of Z Ising (G; y) if −1 < y < 0. Again, they stated their theorem for the case in which y is rational, but the proof applies equally well when y is an algebraic real number. In terms of our notation, they proved the following lemma.
If y is real then Z Ising (G; y) is real. Thus, either Z Ising (G; y) = 0, or arg(Z Ising (G; y)) ∈ {0, π}. Hence, approximating the argument within ±π/3 enables one to determine the sign of the real part. Using the connection (1.6) between the Tutte polynomial and the partition function of the Ising model and Lemma 2.2 we immediately obtain the following corollary.
In fact, we can extend Goldberg and Jerrum's #P-hardness interval-shrinking technique from Goldberg & Jerrum (2014) to also obtain #P-hardness for the relaxed version of the problems. We start with a general discussion of interval shrinking. Suppose that we have a linear function f (ε) = −εA + B for positive A and B and that we wish to find a valueε that is very close to the root ε * = B/A. Suppose that we also have an interval [ε , ε ] such that f (ε ) > 0 and f (ε ) < 0. Suppose that ε − ε = (so the interval has length ). Roughly, Goldberg and Jerrum had at hand an oracle for computing the sign of f (ε) (using an oracle for Sign-Real-Tutte(x, y)) and, using this, it is easy to bisect the interval, getting very close to ε * by binary search.
Using an oracle for the relaxed problem Sign-Real-Nonzero-Tutte(x, y) we can compute the sign whenever it is positive or negative, but we receive an unreliable answer for the sign of f (ε) if f (ε) = 0. Nevertheless, we observe that having a reliable answer in this case is not important for the progress of the binary search. If the binary search queries the value of f (ε) and f (ε) = 0 then the reply from the oracle is correct. Otherwise, the bisection technique described above recurses into a sub-interval that contains a zero of the function, as required. Thus, we have the following lemma. (We omit the formal proof since the lemma follows immediately from the observation that we have just made.) Lemma 3.5. For any algebraic real number y ∈ (−1, 0), Sign-Real-Nonzero-Tutte(x, y) is #P-hard, where x = 1 + 2/(y − 1). Also, the problems Distance-(π/3)-Nonzero-Arg-Ising(y) and Complex-Apx-Nonzero-Ising(y) are #P-hard.
We next show how to further extend the #P-hardness intervalshrinking technique to obtain #P-hardness for the problem Factor-K-Nonzero-Norm-Ising(y). This requires new ideas, so we will provide more details. Let us return to the discussion of interval shrinking. Let η = 1/21 (the exact value of η is not important, but we fix it for concreteness). Instead of having an oracle for the sign of f (ε) = −εA+B, we only will be able to assume that we have an oracle that, on input ε, returns a valuef (ε) satisfying except that again the valuef (ε) is completely unreliable if f (ε) = 0. Our strategy will be to divide the interval into 10 equal-length sub-intervals [ε i , ε i+1 ] for i ∈ {0, . . . , 9} with ε 0 = ε and ε 10 = ε . (The number 10 is not chosen to be optimal -however, it is easy to see that it suffices. Changing the number of sub-intervals would influence the choice of η above.) We then let s i be the sign (positive, negative, or zero) off (ε what the value of s i will be. However, this is true for at most two values of i. So either s 0 , s 1 , s 2 and s 3 are all positive (in which case ε 2 < ε * and we can recurse on the interval [ε 2 , ε 10 ]) or s 6 , s 7 , s 8 and s 9 are all negative (in which case ε 8 > ε * and we can recurse on the interval [ε 0 , ε 8 ]). Either way, the interval shrinks to 4/5 of its original length. Applying this idea in the proof of (Goldberg & Jerrum 2014, Lemma 1) yields the following.
Lemma 3.6. Suppose that γ 1 and γ 2 are algebraic reals with γ 1 ∈ (−2, −1) and γ 2 ∈ [−2, 0]. Then Factor-( 22 Proof. Apart from the interval shrinking idea discussed above, the proof is similar in structure to the proof of (Goldberg & Jerrum 2014, Lemma 1). We defer some calculations (which are unchanged) to Goldberg & Jerrum (2014) but we provide the rest of the proof to show how to get the stronger result. We use the fact that the following problem is #P-complete. This was shown by Provan & Ball (1983).
Output |{S ⊆ E : S is a minimum cardinality (s, t)-cut in G}|.
We will give a Turing reduction from #Minimum Cardinality (s, t)-Cut to the problem Factor-( 22 Let G, s, t be an instance of #Minimum Cardinality (s, t)-Cut. Assume without loss of generality that G has no edge from s to t. Let n = |V (G)| and m = |E(G)|. Assume without loss of generality that G is connected and that m ≥ n is sufficiently large. Let k be the size of a minimum cardinality (s, t)-cut in G and let C be the number of size-k (s, t)-cuts.
Let q = 2 and M * = 2 4m . Let h be the smallest integer such that (γ 2 + 1) h − 1 > M * and let M = (γ 2 + 1) h − 1. Note that we can implement M from γ 2 via an h-thickening, and h is at most a polynomial in m.
Let δ = 4 m /M . Let M be the constant weight function which gives every edge weight M . We will use the following facts: Fact (3.7) follows from the fact that each of the (at most 2 m ) terms in Z st (G; q, M ), other than the term with all edges in A, has size at most M m−1 q n and 2 m M m−1 q n ≤ δM m q. Fact (3.8) follows from the fact that all terms in Z s|t (G; q, M ) are complements of (s, t)-cuts. If more than k edges are cut then the term is at most M m−k−1 q n and For a parameter ε in the open interval (0, 1) which we will tune later, let γ = −1 − ε ∈ (−2, −1). We will discuss the implementation of γ later. Let G be the graph formed from G by adding an edge from s to t. Let γ be the edge-weight function for G that assigns weight M to every edge of G and assigns weight γ to the new edge. Using the definition of the (random cluster) Tutte polynomial, Goldberg and Jerrum noted that It is easily checked that Z Tutte (G ; 2, γ) is positive if ε is sufficiently small (ε = M −2m will do) and it is negative at ε = 1. Thus, viewing Z Tutte (G ; 2, γ) as a function of ε, we can perform interval shrinking (as discussed before the statement of the lemma) to find a value of ε for which Z Tutte (G ; 2, γ) is very close to 0. The interval shrinking uses an oracle for Factor-( 22 21 )-Nonzero-Norm-2Tutte(γ 1 , γ 2 ).
If we find an ε where Z Tutte (G ; q, γ) = 0, then for this value of ε, we have εZ st (G; q, M ) = Z s|t (G; q, M ) 1 − 1+ε 2 . Thus, using ε, we can calculate the fraction Z s|t (G; q, M )/Z st (G; q, M ). Plugging this (known) value into (3.7) and (3.8), we obtain Now, we don't know k, but C is an integer between 1 and 2 m , whereas M > 2 4m , so there is only one value of k that gives a solution C in the right range. Using the value of k, we can calculate C exactly. Technical issues arise both because we are somewhat constrained in what values ε we can implement and because we won't be able to discover the exact value of ε that we need (but we will be able to approximate it closely). These technical issues provide no more difficulty than they did in Goldberg & Jerrum (2014). Suppose first that we are able, for any given ε ∈ (M −2m , 1) to implement γ = −1−ε. Then our basic strategy is to do the interval shrinking, repeatedly sub-dividing the current interval Θ(log(M m 2 )) times, so eventually we'll get an interval of width at most M −m 2 that contains an ε where Z Tutte (G ; 2, γ) = 0. Goldberg & Jerrum (2014) have already shown that knowing such an interval enables the exact calculation of C (so having a small interval is OK -it is not necessary to know ε exactly).
The only issue, then, is implementing the weights γ = −1 − ε during the interval shrinking. As in Goldberg & Jerrum (2014) we cannot expect to implement any particular desired γ precisely. However, using stretching and thickening, we can implement a value that is within an additive error of M −m 2 /20 of any desired ε, and this suffices. The fact that we have algebraic, rather than rational, numbers is irrelevant since stretchings and thickenings can be computed on algebraic numbers.
Using stretching and thickening, we get the following corollary.
Using Lemma 2.2 and the trivial reduction from Factor-K-Nonzero-Norm-Ising(y) to Factor-K-Norm-Ising(y) and from Complex-Apx-Nonzero-Ising(y) to Complex-Apx-Ising(y) we get the following.
The following lemma enables us to determine the complexity of evaluating the Ising partition function when the complex edge interaction y ∈ Q is on the unit circle.
The hardness on the unit circle extends directly to the whole imaginary axis.
Lemma 3.15. Suppose y = ri and r = 0, ±1 where r is algebraic. There exists an algebraic real number y ∈ (−1, 0) that can be implemented by a sequence of stretchings and thickenings from y.
Otherwise suppose |y| > 1. We know that a k-stretch yields the weight z k = 1+2/(x k −1) where x = 1+2/(y −1) = (y +1)/(y −1). Re-arranging, we find that z k = (y+1) k +(y−1) k (y+1) k −(y−1) k . We will now argue that z k is purely imaginary. To see this, note that monomials in the numerator all have degrees of the same parity as k, whereas those in the denominator have degrees of the same parity as k − 1. Therefore, it must be the case that the numerator is real and the denominator is purely imaginary, or vice versa. In either case, z k is purely imaginary. Therefore, if we can find a positive integer k such that 0 < |z k | < 1 then we have reduced our problem to the previous case.
Combining Lemma 3.15 with Observation 2.7, Corollary 3.10, Lemma 3.5 and Corollary 3.11, we get the following corollary.
Finally, this hardness can be extended to some algebraic complex numbers off of the unit circle.
Lemma 3.17. Let y = re iθ be an algebraic complex number such that r > 0 and θ = aπ 2b , where a and b are two co-prime positive integers and a is odd. There exists an algebraic real number y ∈ (−1, 0) that can be implemented by a sequence of stretchings and thickenings from y.
Proof. If r = 1 then we are done by Lemma 3.13. Otherwise r = 1 and by a b-thickening it reduces to the case of Lemma 3.15.
Corollary 3.18. Let y = re iθ be an algebraic complex number such that r > 0 and θ = aπ 2b , where a and b are two co-prime positive integers and a is odd. Then for any K > 1, Factor-K-Nonzero-Norm-Ising(y), Distance-(π/3)-Nonzero-Arg-Ising(y) and Complex-Apx-Nonzero-Ising(y) are #P-hard. Hence, so are the un-relaxed versions of all three problems.
To obtain obtain NP-hardness results for other values of y, we start with the well-known NP-hard problem Max-Cut.
Name Max-Cut.
Instance A (multi)graph G and a positive integer b.
Output Is there a cut of size at least b.
Proof. We will reduce Max-Cut to Factor-K-Nonzero-Norm-Ising(y). Given a graph G and a constant b, we want to decide whether G has a cut of size at least b. We do a k-thickening on G, where k is the least positive integer such that 2 m |y| k < 1/4. Then the effective edge weight is y k = y k . Clearly |y k | = |y| k < 1.
Suppose the maximum cut of G has size c. Now rewrite (1.1) as where m is the number of edges in G and C i is the number of configurations under which there are exactly i bichromatic edges. Since the maximum cut of G has size c and G has m edges, m−c i=0 C i = 2 m . Also, since 2 m |y k | < 1, the i = c term dominates the sum, so Z Ising (G; y k ) is not equal to 0.
If c ≥ b, then our choice of k together with the triangle inequality implies that again by the triangle inequality and 2 m |y k | < 1/4. Therefore we could solve Max-Cut in polynomial time using an oracle for Factor-1.1-Nonzero-Norm-Ising(y k ). By Observation 2.7 it suffices to use an oracle for Factor-1.1-Nonzero-Norm-Ising(y). By Lemma 3.2, an oracle for Factor-K-Nonzero-Norm-Ising(y) will do. Finally, Lemma 2.2 gives the result for Complex-Apx-Nonzero-Ising(y).
The other case, when the norm of y is larger than 1, can be shown to be NP-hard by reduction from the previous case, unless the edge weight is real.
Lemma 3.20. Suppose K > 1. Let y be an algebraic complex number such that |y| > 1 and y ∈ R. Then Factor-K-Nonzero-Norm-Ising(y) is NP-hard and so is Complex-Apx-Nonzero-Ising(y).
Proof. We will prove that there exists a positive integer k such that the effective weight y k of a k-stretch satisfies |y k | < 1. Then we are done by Lemma 3.19.
3.3. Proof of Theorems Theorem 1.2 and Theorem 1.3. Theorems Theorem 1.2 and Theorem 1.3 follow from the following combined theorem. The hardness result in Item Theorem 1.2(iii) of Theorem 1.2 (and its counterpart in Theorem 1.3) follows from Item Theorem 3.21(ix) of the combined theorem.
(iv) If y < −1 is a real number then Factor-K-Nonzero-Norm-Ising(y) is equivalent in complexity to the problem of approximately counting perfect matchings in graphs and Complex-Apx-Nonzero-Ising(y) is as hard. However, Distance-(π/3)-Arg-Ising(y) is in FP.
Proof. Item (i) is from Jaeger et al. (1990). The randomised algorithm for Factor-K-Norm-Ising(y) referred to in Item (ii) is from Jerrum & Sinclair (1993). See also Lemma 3.1 and the surrounding text for a discussion of algebraic numbers and accuracy parameters. The same algorithm can be used for Complex-Apx-Ising(y) because Z Ising (G; y) is real and positive so an approximation N satisfing 1 − 1 R N ≤ Z Ising (G; y, λ) ≤ 1 + 1 R N also satisfies d( N , Z Ising (G; y, λ)) ≤ 1 R . The deterministic algorithm referred to in Items (ii) and (iii) is trivial because the argument of a positive real number is 0. The approximation equivalence in Item (iv) is from Goldberg & Jerrum (2008), since one can decide in polynomial time the existence of perfect matchings to lift the non-zero restriction. The hardness for Complex-Apx-Nonzero-Ising(y) follows from Lemma 2.2. The deterministic sign algorithm in Item (iv) is from Goldberg & Jerrum (2014). Item (v) is from Lemma 3.5 and Corollary 3.10 and Lemma 2.2. Item (vi) is from Corollary 3.14. Item (vii) is from Corollary 3.16. Item (viii) is from Corollary 3.18. Item (ix) is from Lemma 3.19. Finally, item (x) is from Lemma 3.20.

Quantum circuits and counting complexity
In this section we explain the connection between quantum computation and complex weighted Ising models. We begin with some basic notions about quantum circuits. We view qubits |0 and |1 as column vectors [ 1 0 ] and [ 0 1 ]. Similarly 0| and 1| are row vectors (1, 0) and (0, 1). For x ∈ {0, 1} n , let |x denote the tensor product ⊗ n j=1 |x j and x| is similar. Suppose C is a quantum circuit on n qubits and consists of m quantum gates U 1 , . . . , U m sequentially. A quantum gate is a function taking k input and k output variables and returning a value in C. Such a gate is called k-local and has a natural 2 k by 2 k square unitary matrix representation. In a circuit we also need to specify on which qubits the gate acts upon. To make the notation uniform we view unaffected qubits as simply copied and associate each quantum gate with the following 2 n by 2 n square unitary matrix. Let U be a quantum gate and x, y ∈ {0, 1} n two vectors specifying the input and output on all n qubits. Define the 2 n by 2 n matrix M U corresponding to gate U as M U ;x,y = U (x, y).
For example, let H be the Hadamard gate 1 √ 2 [ 1 1 1 −1 ] acting on the first qubit and suppose there are two qubits in total, illustrated as in Using this notation, given an input x ∈ {0, 1} n , the output of the quantum circuit C is a random variable Y subject to the distribution where y ∈ {0, 1} n . It is not necessary that we measure all qubits in the output. We may measure a subset I of all n qubits. Let y ∈ {0, 1} s where |I| = s. Then the output is a random variable Y subject to the distribution Pr C;I (Y = y ) = z∈{0,1} n such that z| I =y Pr C (Y = z). (4.2) Alternatively, we may treat such marginal probability in the counting perspective, as a partition function in the "sum of product" fashion. First let us consider composing two quantum gates, say U 1 and U 2 . Let the input variables of U 1 be x 1 , . . . , x n . Let z 1 , . . . , z n be the variables on the wires between U 1 and U 2 . Finally, let y 1 , . . . , y n be the outputs of U 2 . We use σ(x) to denote an assignment of values in {0, 1} to the variables x 1 , . . . , x n . We use σ(y) and σ(z) similarly. Then the composition U of U 1 followed by U 2 is given by U 1 (x, σ(z))U 2 (σ(z), y).
x 4 z 4 y 4 Figure 4.2: Two quantum gates U 1 and U 2 composed together. Figure 4.2 illustrates the composition of gate U 1 acting upon qubits 2, 3, 4 followed by U 2 acting upon 1, 2. In the matrix notation, it is easy to see that We now associate an intermediate variable z j,k to each edge on qubit k between gate U j and U j+1 for all 2 ≤ j ≤ m − 1 and 1 ≤ k ≤ n. Denote by z j the vector {z j,k | 1 ≤ k ≤ n} and z = ∪ m−1 j=2 z j . As the initial input and output of a quantum circuit are column vectors and row vectors respectively, they may be treated as function/gates with no output variables or no input variables. In particular, on the product input state |x input variables are set to {x k } where x ∈ {0, 1} n . Using (4.3) recursively we can rewrite (4.1) as follows: To simulate classically a quantum circuit, one can either (approximately) compute the probability Pr C (Y = y) -this is called "strong simulation" -or one can sample from a distribution that is sufficiently close to the one given by (4.1) or (4.4). This is called "weak simulation" 4.1. IQP and the Ising partition function. IQP, which stands for "instantaneous quantum polynomial time", is characterised by a restricted class of quantum circuits introduced by Shepherd & Bremner (2009). Bremner et al. (2011 showed that if IQP can be simulated classically in the sense of "weak simulation" with multiplicative error, then the polynomial hierarchy collapses to the third level. Fujii & Morimae (2013) showed that the marginal probabilities of possible outcomes of IQP circuits correspond to partition functions of Ising models with complex edge weights.
The key property of IQP is that all gates are diagonal in the |0 ±|1 basis. Therefore all gates are commutable. In other words, there is no temporal structure and hence it is called "instantaneous". Let H be the Hadamard gate 1 √ 2 [ 1 1 1 −1 ]. If a gate U is diagonal in the |0 ± |1 basis, there exists a diagonal matrix D such that M U = H ⊗n DH ⊗n . Moreover H is its own inverse; That is, HH = I 2 . Any two H's between each pair of gates cancel. This leads to an alternative view of IQP circuit in which each qubit line starts and ends with an H gate and all gates in between are diagonal.
Definition 4.5. An IQP circuit on n qubit lines is a quantum circuit with the following structure: each qubit line starts and ends with an H gate, and all other gates are diagonal.
We will focus particularly on 1, 2-local IQP, which means that every intermediate gate acts on 1 or 2 qubits. It was shown that a classical weak simulation of 1, 2-local IQP with multiplicative error implies the polynomial hierarchy collapse to the third level Bremner et al. other than H gates on two ends of each line. We will show that this class of IQP circuits corresponds to Ising models with complex edge interactions and that therefore the strong simulation of these circuits is #P-hard, even allowing an error of any factor K > 1.
To show the relationship between these circuits and Ising partition functions, it is convenient to use another set of gates. Let Hence we can replace every CZ gate on qubits j, k by 2 copies of R π/8 on j, k, 14 copies of P π/8 on qubit j, and 14 P π/8 on qubit k. It is easy to see that R π/8 can be replaced by CZ and P π/8 as well. We may therefore assume every gate is either P π/8 on 1 qubit or R π/8 on 2 qubits without changing the computational power of the circuit. In general we give the following definition.
Definition 4.8. An IQP 1,2 (θ) circuit on n qubit lines is a quantum circuit with the following structure: each qubit line starts and ends with an H gate, and every other gate is either P θ on 1 qubit or R θ on 2 qubits. We assume the input state is always |0 n .
An example IQP 1,2 (θ) circuit is given in Figure 4.3. The relationship between IQP 1,2 (θ) circuits and Ising models was first observed by Fujii & Morimae (2013). These connections will be shown next. For completeness we include our own proofs, which have a more combinatorial flavour than the original ones by Fujii & Morimae (2013). We introduce the following non-uniform Ising model which has been studied previously. See, for example Sokal (2005). Let G = (V, E) be a (multi)graph. The edge interaction is specified by a function ϕ : E → C and the external Figure 4.3: An IQP 1,2 (θ) circuit. We use two solid dots to denote R θ gate as it is diagonal and symmetric.
field is specified by a function τ : V → C. The partition function is defined as where δ(x, y) = 1 if x = y and δ(x, y) = 0 if x = y. We write Z Ising (G; y, τ ) when ϕ(e) = y is a constant function and similarly Z Ising (G; ϕ, λ) when τ (v) = λ. Notice that this notation is consistent with (1.1). We will show that the following problem is related to Factor-K-Strong-Sim-IQP 1,2 (θ) when e iθ is a root of unity.
Instance A (multi)graph G with an edge interaction function ϕ(−) taking value e iθ or e −iθ , and an external field function τ so that for each vertex v there are non-negative integers a v and b v so that τ (v) = (−1) av e iθ bv or τ (v) = (−1) av e −iθ bv .
We will first consider inputs to IQP 1,2 (θ) where I = [n] so all qubits are measured. Given an IQP 1,2 (θ) circuit C on n qubits and a string y ∈ {0, 1} n , we can construct a non-uniform Ising instance G C with edge interaction e i2θ and external field τ C;y such that The construction is as follows. The vertex set {v j } contains n vertices and each vertex corresponds to a qubit. For each gate R θ on two qubits j, k, add an edge (j, k) in G C . For qubit j, let p j be the number of gates P θ acting on qubit j in C. Let τ C;y (v j ) = e −i(2p j θ) (−1) y j . An example of the construction is given in Figure 4.4.
Lemma 4.11. Let C be an IQP 1,2 (θ) circuit on n qubits and y ∈ {0, 1} n be the output. Let G C and τ C;y be constructed as above. Then (4.10) holds.
Proof. Suppose C is composed sequentially by U 1 = H ⊗n , U 2 , . . . , U m−1 , U m = H ⊗n , where U j is either P θ on 1 qubit or R θ on 2 qubits for 2 ≤ j ≤ m − 1. Notice that U 1 (x, x ) = U m (x, x ) = 2 −n/2 n k=1 (−1) x k x k . As the input |x = |0 n , we can rewrite (4.4): Let Q denote the quantity inside the norm, that is, Since U j 's are diagonal for 2 ≤ j ≤ m − 1, any configuration σ with a non-zero contribution to Q must satisfy that for any k, σ(z 1,k ) = σ(z 2,k ) = · · · = σ(z m−1,k ). Therefore we may replace z j,k by a single variable v k for all 1 ≤ j ≤ m − 1 so that Moreover, if U j is the gate P θ on qubit k, then U j (σ(V ), σ(V )) = e iθ e −i2θ σ(v k ) . If U j is the gate R θ on qubits k 1 and k 2 , then U j (σ(V ), σ(V )) = e −iθ e i2θ δ(σ(v k 1 ),σ(v k 2 )) , where δ(x, y) = 1 if x = y and δ(x, y) = 0 if x = y. Recall that p k is the number of P θ gates on qubit k and τ C;y (v k ) = e −i(2p k θ) (−1) y k . Collecting all the contributions, we have where m j is the number of j qubit(s) gates for j ∈ {1, 2}, and, from (1.1), m(σ) is the number of monochromatic edges under σ. We get (4.10) by substituting (4.13) in (4.12).
Similar results hold when some qubits are not measured. To show it, we need the following fact. It can be viewed as an application of Parsevals's identity on the length-2 n vector {C z } indexed by z ∈ {0, 1} n over an orthonormal basis {e z } where basis element e z has value 2 − n 2 (−1) z·z in position z . We include a proof for completeness.
Claim 4.14. Let {C z } be 2 n complex numbers where z runs over {0, 1} n . Then we have Proof. Notice that for two complex numbers A and B, |A + B| 2 + |A − B| 2 = |A| 2 + |B| 2 − 2|A||B| cos θ where θ is the angle from A to B. Hence we have where in the last line we apply (4.15). The claim holds by induction.
We then have the following reduction.
Proof. If all qubits in the input to Factor-K-Strong-Sim-IQP 1,2 (θ) are measured, then the result follows from Lemma 4.11. Otherwise, without loss of generality we assume the first n − s qubits are measured. Let C, I = [n − s] and y ∈ {0, 1} n−s be the input to Factor-K-Strong-Sim-IQP 1,2 (θ). We use (4.2), (4.12), and the first line of (4.13): (4.17) where for z ∈ {0, 1} s , Q z is the contribution of assigning z l−n+s to v l without the possible −1 external field, that is, Apply Claim Claim 4.14 on (4.17): We construct the following instance of Factor-K 2 -Norm-IQP-Ising(2θ). We first construct G C = (V, E) with edge interaction e i2θ as before. The vertex set {v j } contains one vertex for each of the n qubits. For each gate R θ on two qubits j, k we add edge (j, k) with edge interaction e i2θ to G C . Now make a copy G C = (V , E ) such that the edge interaction is e i2θ = e −i2θ . Let ϕ C;I be this edge interaction function. Then we identify vertices v l with v l for all n − s + 1 ≤ l ≤ n. Let U be the set of these identified vertices and let V 1 = V − U and V 1 = V − U . The external field τ = τ C;I,y is defined as follows: for any v ∈ U , τ (v) = 1; for any v j ∈ V 1 , τ (v j ) = e −i(2p j θ) (−1) y j ; and for any v j ∈ V 1 , τ (v j ) = τ (v j ) = e i(2p j θ) (−1) y j . Informally, this instance was formed by putting G C and its complement together and identifying vertices that correspond to unmeasured qubits. Note that if two vertices in U are connected by an edge, then they are actually connected by two edges, and the product of the two edge interactions is 1. We therefore remove all edges with both endpoints in U . Call the resulting graph H C . One can verify that (H C , ϕ C;I , τ C;I,y ) is a valid instance of Factor-K 2 -Norm-IQP-Ising(2θ). An example of the construction is given in Figure 4.5.
Fix an assignment z ∈ {0, 1} s on U . The contribution Z z to Z Ising (H C ; ϕ C;I , τ C;I,y ) can be counted in two independent parts, V and V . Hence we have where given the configurations σ 1 (or σ 1 ), m * (σ 1 , z) (or m * (σ 1 , z)) is the number of monochromatic edges with at least one endpoint in V (or V ). Recall that τ (v j ) = e −i(2p j θ) (−1) y j . Comparing Z z to |Q z | 2 , the only difference is that in |Q z | 2 , e i2θ is raised to the number of monochromatic edges in the whole V instead of V 1 . However for any monochromatic edge in U , its contribution is independent from the configuration σ, and hence can be moved outside of the sum. All such terms are cancelled after taking the norm. This implies Z z = |Q z | 2 . Therefore (4.18) can be rewritten as Pr C;I (Y = y ) = 2 −2n+s z∈{0,1} s Z z = 2 −2n+s Z Ising (H C ; ϕ C;I , τ C;I,y ) = 2 −2n+s |Z Ising (H C ; ϕ C;I , τ C;I,y )|. The lemma follows from the above equation.
Remark 4.20. In fact, the construction of H C can be further simplified. If v ∈ V and v ∈ V connect to some u ∈ U , we can replace edges (u, v) and (u, v ) by a new edge (v, v ) with an Ising interaction 2 e i4θ +e −i4θ . (In case e i4θ + e −i4θ = 0 this interaction is equality and we identify v with v .) Therefore we can reduce an instance of Factor-K-Strong-Sim-IQP 1,2 (θ) to an Ising model of size linear in |I|, the number of measured qubits. If |I| = O(log n), then the reduced Ising instance is tractable and so is the simulation. This matches the strong simulation result by Shepherd (see (Bremner et al. 2011, Theorem 3.4) , the remark following that theorem and also Shepherd (2010).) The reduction also works in the other direction when e iθ is a root of unity.
Proof. Lemma 4.16 implies a reduction from the right hand side to the left hand side. In the rest of the proof we show the other direction. As e iθ is a root of unity, there exists a positive integer t such that e −i2θ = e i2tθ . Given an instance (G, ϕ, τ ) of Factor-K-Norm-IQP-Ising(2θ), we may replace each edge of interaction e −i2θ by t parallel edges of weight e i2θ . Moreover, we may assume the external field is of the form τ (v j ) = (−1) a j e −i2θ b j for the same reason. We construct an IQP 1,2 (θ) circuit C on n = |V | qubits. For each edge (v j , v k ) ∈ E, we add a quantum gate R θ on qubits j and k. For each 1 ≤ j ≤ n, we add b j many quantum gate P θ on qubits j and let the output y j = 1 on qubit j if a j is odd. By Lemma 4.11 we see that 2 2n Pr C (Y = y) = Z Ising (G; e i2θ , τ ) 2 .
Suppose the Ising instance in the proof of Theorem 4.21 has no external field and has a constant edge interaction e i2θ . Then it is not hard to see that the above construction does not rely on e iθ being a root of unity and works for general θ. Hence we have the following lemma.
We can now prove our main result about IQP.
Proof. This follows from Lemma 4.22 and Corollary 3.14.
In a related result, (Bremner et al. 2011, Corollary 3.3) showed that weakly simulating IQP with multiplicative error implies that the polynomial hierarchy collapses to the third level. More precisely, their result is the following. Suppose C is an IQP 1,2 (π/8) circuit on n qubits. If there exists a classical randomized polynomial time procedure to sample a binary string Z of length n, such that for every string y ∈ {0, 1} n and any constant 1 ≤ K < √ 2, then the polynomial hierarchy collapses to the third level. The usual measure for determining the quality of a sampling procedure is total variation distance. The notion of total variation distances is weaker than "multiplicative error" so the result in Bremner et al. (2011) does not rule out weak simulation with small variation distance. To see this, note that, if the multiplicative error is K, then obviously the total variation distance is at most K−1. On the other hand, consider two distributions supported by two n-bit Boolean strings. A sample from the first distribution is obtained uniformly choosing each of the n bits. A sample from the second distribution is obtained by uniformly choosing each of the first n − 1 bits. The last bit is 1 if all other bits are 0, and is chosen uniformly otherwise. The total variation distance is 2 −n , but the multiplicative error is infinity at the all 0 string. Note that the complexity implication "polynomial hierarchy collapses to the third level" is apparently weaker than the consequence of strong simulation from Theorem 1.4, which is FP = #P. Strong simulation is also studied with respect to other classes of quantum circuits, see for example Jozsa & Van den Nest (2014).
The allowable error is usually taken to be additive and exponentially small, instead of the constant factors that we have studied here. For example, Jozsa & Van den Nest (2014) requires that the output be computed with k bits of precision in an amount of time that is polynomial in both k and the size of the input. Additive error is quite different from multiplicative error. Also, the amount of accuracy is important. Lemma 3.2 shows that there is no difference between a constant factor and an FPRAS scenario, in which the error is allowed to be a factor of 1 ± 1/R for a unary input R. On the other hand, achieving a multiplicative error of 1±1/ exp(R) is an entirely different matter. Bordewich et al. (2005) raised the question "of determining whether the Tutte polynomial is greater than or equal to, or less than zero at a given point." Thus, they raised the question of determining the complexity of Sign-Real-Tutte(x, y). In fact, they were especially interested in the case x = −t, y = −t −1 where t = exp(2πi/5).

BQP and the Tutte polynomial
We next show that resolving this case is a simple corollary of our results. After that, we will discuss the motivation for considering this point (x, y) and its connection to the complexity class BQP. We will also briefly discuss a relevant general result of Kuperberg (2015), which resolves similar questions by using three results about quantum computation -the Solovay-Kitaev theorem, the FLW density theorem, and a result of Aaronson.
Motivated by connections to quantum computing, we consider the difficulty of the problem Sign-Real-Tutte(x, y) when xy = 1. In particular, we study the points (x, y) = (exp(−aπi/b), exp(aπi/b)), where a and b are positive integers. If a ∈ {0, b/2, b, 3b/2} then the problem is trivial since (x, y) is one of the so-called "special points" ((1, 1), (−1, −1), (−i, i) and (i, −i)) where evaluating the Tutte polynomial is in FP Jaeger et al. (1990). We can assume without loss of generality that a < 2b since adding 2π to the argument of a complex number doesn't change anything. We can now prove the main result of this section.
We implement (x , y ) using a b-thickening from (x, y). Then, since a is an odd positive integer, y = y b = exp(aπi) = −1.
So x = 1 + q/(y − 1) = 1 − q/2 = cos(aπ/b). Now since x < 11/27, (Goldberg & Jerrum 2014, Theorem 1, Region F) shows that computing the sign of Z Tutte (−; q, y − 1) is #P-hard. As we showed in the argument that established Lemma 3.5 (see the paragraph before the statement of the lemma), the same is true if the oracle returns any answer when the value is 0.
Since x and y are not 1, (2.4) shows that it is also hard to compute the sign of T (x , y ). The result now follows from Observation 2.7.
Theorem 1.7 is very close to a special case of the following result of Kuperberg. A link is a collection of smooth simple closed curves embedded in 3-dimensional space. V L (t) denotes the Jones polynomial of a link L evaluated at point t. We do not need the detailed definition of the Jones polynomial in order to state Kuperberg's theorem.
Theorem 5.1. (Kuperberg 2015, Theorem 1.2) Let V (L, t) be the Jones polynomial of a link L described by a link diagram, and let t be a principal root of unit other than exp(2πi/r) where r ∈ {1, 2, 3, 4, 5, 6}. Let 0 < A < B be two positive real numbers and assume as a promise that either |V (L, t)| < A or |V (L, t)| > B. Then it is #P-hard to decide which inequality holds. Moreover, it is still #P-hard when L is a knot.
The connection is as follows. There is a result of Thistlethwaite (1987) (see (Jaeger et al. 1990, (6.1))), showing that when L is an alternating link with associated planar graph G(L), then is an easily-computable factor which is plus or minus a half integer power of t. Thus, the evaluation of Jones polynomial of an alternating link is an easily-computable multiple of an evaluation of the Tutte polynomial along the hyperbola xy = 1 (where, for some value t, x = −t and y = −t −1 ), as in Theorem 1.7. The importance of these evaluations is established in (Bordewich et al. 2005, Theorem 6.1) which shows that all of the problems in the quantum complexity class BQP (consisting of those decisions problems that can be solved by a quantum computer in polynomial time) can also be solved classically in polynomial time using an oracle that returns the sign of the real part of the Jones polynomial of a link, evaluated at the point t = exp(2πi/5) (the point studied in Corollary 1.8).
Kuperberg's theorem (Theorem 5.1) is incomparable to Theorem 1.7. In some respects, Theorem 5.1 is more general -it does not have the restriction cos(aπ/b) < 11/27. Also, G(L) is always planar, which is essential for the connection to BQP, and it applies to a wide range of A and B. On the other hand, the most relevant case A = B = 0 (the one that relates to the BQP result of Bordewich et al. (2005)) is actually excluded from Theorem 5.1 since A and B must be different and positive. We are not sure whether Kuperberg's proof can be adapted to include this case, where the goal would be to determine whether |V (L, t)| ≥ 0 or |V (L, t)| ≤ 0. This is covered by Theorem 1.7.
In any case, it seems interesting to note that the proof of Theorem 1.7 is combinatorial (about Tutte polynomials only) whereas the proof of Theorem 5.1 is essentially about quantum computation. (Kuperberg describes it as "a mash-up of three standard theorems in quantum computation".) We refer the reader to Aharonov & Arad (2011) for more recent results giving BQP-hardness of multiplicative approximations of the Jones polynomial of the plat closure of a braid at roots of unity. Also, we note that other works such as Geraci & Lidar (2010) have suggested the idea of using tractable planar evaluations of these polynomials to give efficient classical simulations for special cases of quantum circuits.

Ising with a field
In Section 6.2, we will extend our Ising hardness results from Theorems Theorem 1.2 and Theorem 1.3 to the situation in which we have an external field λ = 1. To obtain our hardness results, we need a lower bound on the relevant partition functions.
6.1. Lower bounds on partition functions. Suppose we have two edge weights y and y that are close. It is easy to bound the distance between Z Ising (G; y) and Z Ising (G; y ) additively, but not multiplicatively. To convert an absolute error into a relative error, one needs some lower bound on the partition function. However, when the edge interaction y is negative or complex, it is possible that the partition function vanishes. Assuming that it doesn't vanish, we would like to know how close to zero could it get. When y is rational, an exponential lower bound is easy to obtain by a simple granularity argument, but the argument is more difficult when y is not rational. In this section we give an exponential lower bound which is valid when y is an algebraic number. The techniques that we use are standard in transcendental number theory, see e.g. Bugeaud (2004).
We begin with some basic definitions from Bugeaud (2004). For a polynomial with complex coefficients the (naive) height of P (x) is defined as H(P ) := max i {|a i |}. A more advanced tool, its Mahler measure, is defined as There is a standard inequality relating these two measures. It is proved for complex polynomials in (Bugeaud 2004, Lemma A.2). For completeness, we include the proof (following Bugeaud (2004)) for the case in which P (x) is a real polynomial, which is all that we require. a j a k (cos(j · 2πt) cos(k · 2πt) The claim holds as M(P ) ≤ ( n i=0 a 2 i ) 1/2 ≤ √ n + 1 H(P ).
Let y ∈ C be an algebraic number and its minimal polynomial over Z is P y (x). The degree of P y (x) is called the degree of y and H(P y ) is called the height of y, also denoted H(y).
We also need the following notion of resultants.
Definition 6.2. Let P (x) = a n n i=1 (x−α i ) and Q(x) = b m m i=1 (x − y i ) be two non-constant polynomials. The resultant of P (x) and Q(x) is defined as It is a standard result that Res(P, Q) is an integer polynomial in the coefficients of P (x) and Q(x). The resultant is also the determinant of the so-called Sylvester matrix. In particular, when P (x) and Q(x) are integer polynomials, Res(P, Q) is always an integer, as the Sylvester matrix is an integer matrix in this case. Moreover, we can rewrite the resultant as follows: Res(P, Q) = a m n 1≤i≤n Q(α i ) = (−1) mn b n m 1≤j≤m P (y j ). Now we are ready to give a lower bound for any integer polynomial evaluated at an algebraic number. It is a standard result in algebraic number theory. For completeness we provide a proof here and the treatment is from (Bugeaud 2004, Theorem A.1). Lemma 6.3. Let P (x) be an integer polynomial of degree n, and y ∈ C be an algebraic number of degree d. Then either P (y) = 0 or |P (y)| ≥ C −n y ((n + 1)H(P )) −d+1 .
where C y > 1 is an effectively computable constant that only depends on y.
Suppose there is an j = 1 such that P (y j ) = 0. As Q(x) is the minimal polynomial of y, none of y j could be a rational number. Hence there is an automorphism of the splitting field of Q(x) that maps y j to y. Applying this automorphism on both sides of P (y j ) = 0, we get P (y) = 0. Contradiction! Hence we have P (y i ) = 0 for all i and the resultant of P (x) and Q(x) is non-zero. Since Res(P, Q) is an integer, we have 1 ≤ |Res(P, Q)| = |b d | n 1≤i≤d |P (y i )|.
Lemma 6.4. Let G be a graph and y ∈ C a non-zero algebraic number of degree d. There exists a positive constant C > 1 depending only on y such that if Z Ising (G; y) = 0, then |Z Ising (G; y)| > C −m , where m is the number of edges in G.
Proof. Given a graph G, first suppose that G is not connected, G i 's are the components of G. Then Z Ising (G; y)= i Z Ising (G i ; y).
It is easy to see that if the claim holds for all components it hold for G as well. Therefore in the following we may assume G is connected. Then m ≥ n − 1 where n is the number of vertices. We can rewrite Z Ising (G; y) as a polynomial in y as follows, where C j is the number of configurations such that there are exactly j many monochromatic edges. Notice that m j=0 C j = 2 n , we have H(P ) ≤ 2 n . Assume P (y) = 0. Apply Lemma 6.3 and we obtain where C y > 1 is a constant depending only on y. As m ≥ n−1, the right hand side decays exponentially in m and the lemma follows.
Lemma 6.5. Let G be a graph and y, z ∈ C two roots of unity. Let n be the number of vertices in G and m the number of edges. There exists a positive constant C > 1 depending only on y and z such that if Z Ising (G; y, z) = 0, then |Z Ising (G; y, z)| > C −m .
Proof. As in the previous lemma we may assume G is connected and m ≥ n − 1. Suppose y is of order d 1 and z order d 2 . Let d be the least common multiple of d 1 and d 2 . Then there exists a root of unity w of order d such that y = w t 1 and z = w t 2 . Given a graph G, we can rewrite Z Ising (G; y, z) as a polynomial in y and z as follows, where C j,k is the number of configurations such that there are exactly j many monochromatic edges and k many 1 vertices. Let We have H(P ) ≤ 2 n . Assume P (w) = 0. Apply Lemma 6.3 and we obtain where C w > 1 is a constant depending only on w. As m ≥ n−1, the right hand side decays exponentially in m and the lemma follows.
6.2. Hardness results. In this section we will show hardness results when both the edge interaction and external field are roots of unity.
We first consider the external field −1. We describe the edge interaction by specifying an interaction matrix n 00 n 01 n 10 n 11 , where n ij is the weight when the two endpoints have spins i and j, respectively. In this notation, a binary equality is 1 0 0 1 , and an Ising interaction with weight y is y 1 1 y . Given a gadget with two distinguished vertices, we may view it as an edge and compute its effective interaction matrix M . Then we say the gadget implements M . Also, recall the definitions of k-stretch and k-thickening (Observation 2.7, for example).
Proof. We first argue that a binary equality can be implemented. Consider a 2-stretch with the edge interaction y and external field −1. It is easy to calculate that the (effective) interaction matrix is y 2 −1 0 0 1−y 2 . Then do a 2-thickening. The resulting matrix is (y 2 −1) 2 0 0 (1−y 2 ) 2 . Up to a constant of (y 2 − 1) 2 this is equality.
Suppose G = (V, E) is an input to Factor-K-Nonzero-Norm-Ising(y). We introduce a new vertex v for every vertex v ∈ V . Connect v and v via this equality gadget, that is, first a 2-stretch and then a 2-thickening. Hence the external field on v is cancelled with this construction. The reduction follows.
Next we consider the case when a real edge interaction can be implemented. If the norm of the interaction is less than 1, then we can cancel out the external field.
Lemma 6.7. Let K > 1 and K > 1. Let y and z be two roots of unity and z = ±1. Suppose some real number w ∈ (−1, 1) as an edge interaction is implementable for the Ising model with edge interaction y and external field z. Then we have Factor-K-Nonzero-Norm-Ising(y) ≤ T Factor-(KK )-Nonzero-Norm-Ising(y, z).
Suppose w = 0, which means we can implement inequality (see the remark above Lemma 6.6). For each vertex v i , we introduce a new vertex v i and connect v i and v i by the inequality. It is easy to verify that if v i is assigned 0, the weight from v i and v i together is z; when v i is assigned 1, the weight is also z. Hence the external field is effectively cancelled and the reduction follows.
Otherwise assume w = 0, that is w ∈ (−1, 0) ∪ (0, 1). For each vertex v i , we introduce a new vertex v i , and add 2t many new edges between v i and v i , where t is a positive integer which we will choose later. By assumption we can implement the edge interaction w and we put it on all new edges. Let V = {v i |1 ≤ i ≤ n} and we get a new graph G = (V ∪ V , E ).
For each vertex v i , the contribution of v i and v i (to the partition function) together is w 2t + z when v i is assigned 0 and z(1 + w 2t z) when v i is assigned 1. Let λ = z(1+w 2t z) w 2t +z . Notice that w 2t + z = 0 as |w| < 1 = |z|. We have Z Ising (G ; y, z) = (w 2t + z) n σ:V →{0,1} y m(σ) λ n 1 (σ) , where m(σ) is the number of monochromatic edges in E under σ and n 1 (σ) is the number of vertices in V that are assigned 1.
This finishes the proof.
A similar proof works when the implementable real field has a larger than 1 norm. Basically when this is the case we may power the external field z. If z is a root of unity then we could power it to 1. Lemma 6.12. Let K > 1 and K > 1. Let y and z be two roots of unity and z = ±1. Suppose some real number w ∈ (−∞, −1) ∪ (1, ∞) as an edge interaction is implementable for the Ising model with edge interaction y and external field z. Then we have Factor-K-Nonzero-Norm-Ising(y, z r ) ≤ T Factor-(KK )-Nonzero-Norm-Ising(y, z) for any positive integer r.
For each vertex v i , we introduce r − 1 many new vertices v i,j , and add 2t many new edges between v i and each v i,j , where j ∈ [r − 1] and t is a positive integer which we will choose later. By assumption we can implement the edge interaction w and we put it on all new edges. Let V = {v i,j |1 ≤ i ≤ n, 1 ≤ j ≤ r − 1} and we get a new graph G = (V ∪ V , E ).
For each vertex v i , the contribution of v i and all v i,j combined is (w 2t + z) r−1 when v i is assigned 0 and z (1 + w 2t z) r−1 when v i is assigned 1. Let λ = z(1+w 2t z) r−1 (w 2t +z) r−1 . Notice that w 2t + z = 0 as |w| > 1 = |z|. We have Z Ising (G ; y, z) = w 2t + z n(r−1) σ:V →{0,1} y m(σ) λ n 1 (σ) , where m(σ) is the number of monochromatic edges in E under σ and n 1 (σ) is the number of vertices in V that are assigned 1.
We will show how to implement a real edge interaction in the next lemma. Unless the norm of the new interaction is 1, the hardness holds due to the previous two lemmas. The failure cases are indeed polynomial-time computable.
Lemma 6.18. Let K > 1. Let y = ±i and z be a root of unity that is not one of {1, −1, i, −i}. Then Factor-K-Nonzero-Norm-Ising(y, z) is #P-hard.