Keywords

1 Introduction

Physicists and chemists regularly deal with ‘quantum \(\textsf{NP}\)’-hard problems, for example when finding the ground state (energy) of a physical system [30] or assessing the consistency of local density matrices (the quantum analog of deciding the consistency of marginal probability distributions) [32]. Quantum computing not only holds the potential to provide a matching computational resource for tackling these challenges but also serves as a bridge to incorporate classical reasoning techniques for tackling nature’s hardest problems. Quantum circuits, in particular, offer a precise view into these problems, because the quantum circuit equivalence checking problem is also ‘quantum \(\textsf{NP}\)’-hard.

Circuit equivalence [2, 4, 8, 22, 23, 52, 55, 56, 61] also has many important applications. Since quantum computers are highly affected by noise, it is necessary to optimize the circuits to maximize the performance when running them on a real device. Furthermore, many devices can only handle shallow-depth circuits and are subject to various constraints such as connectivity, topology, and native gate sets. An essential aspect of designing and optimizing quantum circuits isverifying whether two quantum circuits implement the same quantum operation.

Equivalence checking for so-called Clifford circuits is tractable [52], which is surprising considering their wide applicability, e.g. in quantum error correction [9, 48, 49]. Extending the Clifford gate set with any non-Clifford gate, however, e.g. with a T or Toffoli gate, makes the problem immediately ‘quantum \(\textsf{NP}\)’-hard, that is: \(\textsf{NQP}\)-hard to compute exactly [51] and \(\textsf{QMA}\)-hard to approximate [24], even for constant-depth circuits [25].Footnote 1 The exact formulation of equivalence checking allows its discretization [29], exposing the underlying combinatorial problem that classical reasoning methods excel in. Indeed, exact reasoning methods based on decision diagrams are even used to compute the approximate version of the problem (see e.g. [23, 57]).

Our aim is to use reasoning tools based on satisfiability (SAT) for exact equivalence checking of universal quantum circuits. Like SAT solvers [7, 16], model counters, or #SAT solvers, can handle complex constraints from industrial-scale applications [40, 47], despite the #P-completeness of the underlying problem.

We propose a new equivalence-checking algorithm based on weighted model counting (WMC). To do so, we generalize the WMC encoding of quantum circuit simulation from [34], showing that it essentially only relies on expressing quantum information in the so-called Pauli basis [18], thus obviating the need for the arguably more complex stabilizer theory [20, 63]. In addition, we extend the encoding with support for the (non-Clifford) Toffoli gate, allowing more efficient encodings for many circuits. We then prove that a folklore theorem on quantum circuit equivalence checking [52] enables the reduction of the problem to a sequence of weighted Boolean formulas that can be solved using existing weighted model counters (provided they support negative weights [34]).

We show how the WMC encoding satisfies the conditions of the theorem from [52] and implement the proposed equivalence checking algorithm in the open-source tool ECMC, which uses the weighted model-counting tool GPMC [50].Footnote 2 To assess the scalability and practicality of ECMC, we conduct experimental evaluations using random Clifford+T circuits which closely resemble quantum chemistry applications [59] and various quantum algorithms from the MQT benchmark [43], which includes important quantum algorithms such as QAOA, W-state, and VQE among others. We compare the results of our method against that of the state-of-the-art circuit equivalence checker QCEC  [8], showing that in several cases the WMC approach used by our ECMC tool is competitive.

In summary, this paper provides a many-to-many reduction of (universal) quantum circuit equivalence to weighted model counting (WMC). As a consequence, we contribute additional new benchmarks for the WMC competition: basically, each pair of universal quantum circuits can be reduced to a sequence of weighted CNF encodings that need to be solved to (dis)prove equivalence. This opens up numerous possibilities and challenges to better adapt model counters for this new application area in quantum computing.

2 General Background

We only provide the necessary background. For a more complete description see the full version of this paper [35].

Quantum Computing. We fix n as the number of qubits in the circuit(s) under consideration and write [m] for the set \(\left\{ 1, \dots , m \right\} \). Qubits are numbered as [n]. We represent an n-qubit quantum state \(|{\varphi }\rangle \in \mathbb {C} ^{2^n}\) as its density matrix \(\vert \varphi \rangle \!\langle {\varphi }\vert \in \mathbb {C} ^{2^n} \times \mathbb {C} ^{2^n}\), where \(\langle \varphi \vert \) represents the conjugate transpose \(|{\varphi }\rangle ^\dagger \) of \(|{\varphi }\rangle \) [38].

A quantum gate G on n qubits can be expressed by a \(2^n\times 2^n\) complex matrix \(U_G\) which is unitary, i.e. \(U_G\) is invertible and satisfies \(U_G^\dag = U_G^{-1}\). If a quantum state is represented by a density matrix \(\vert \varphi \rangle \!\langle {\varphi }\vert \), then the density matrix after applying G is given by conjugation of \(\vert \varphi \rangle \!\langle {\varphi }\vert \), i.e. \(U_G\vert \varphi \rangle \!\langle {\varphi }\vert U_G^{\dagger }\). For an n-qubit quantum system, applying a single-qubit gate U on the j-th qubit is represented by

$$\begin{aligned} U_j = I^{\otimes j -1}\otimes U \otimes I^{\otimes n- j}, \end{aligned}$$
(1)

where I is the single-qubit identity matrix and \(\otimes \) denotes the Kronecker product. A circuit in our text is simply a list of n-qubit unitaries, i.e., \(C = (G^0, \dots , G^{m-1})\) where C can in turn be understood as unitary itself \(U_C = U_{G^{m-1}} \cdot U_{G^{m-2}} \cdots U_{G^0}\). We will sometimes refer to a gate or circuit as its unitary, and vice versa, because it is clear from context which is meant.

The gates

$$ H = \frac{1}{\sqrt{2}} \begin{bmatrix} 1&{}1\\ 1&{}-1\end{bmatrix}, \quad S = \begin{bmatrix} 1&{} 0\\ 0&{}i\end{bmatrix}, \quad CZ = \left[ {\begin{matrix} 1&{} 0&{}0 &{}0\\ 0&{} 1&{} 0 &{} 0\\ 0 &{} 0 &{} 1 &{} 0\\ 0 &{} 0 &{} 0 &{} -1\end{matrix}}\right] $$

form the so-called Clifford (generating) set.

Though non-universal and classically simulatable [1], Clifford circuits, i.e., circuits composed of Clifford gates only, are expressive enough to describe entanglement, teleportation and superdense coding, and are used in quantum error-correcting codes [9, 48, 49] and in measurement-based quantum computation [46]. Nonetheless, even equivalence checking of Clifford circuits is in \(\textsf{P}\) [52]. By extending the Clifford gate set with any non-Clifford gate, such as the \(T = \sqrt{S}\), Toffoli or arbitrary rotation gates \(R_X,R_Y,R_Y\), we immediately obtain a universal gate set, in the sense that arbitrary unitaries can be approximated [14, 30, 31].

In this work, we express matrices not in the standard basis but in the Pauli basis. We define the \(2 \times 2\) Pauli matrices X, Y, Z, together with identity, as:

$$ \begin{aligned} \sigma [00] \equiv \,I\, \equiv \left[ {\begin{matrix} 1 &{} 0 \\ 0 &{} 1 \end{matrix}}\right] ,~ \sigma [01] \equiv Z \equiv \left[ {\begin{matrix} 1 &{} 0 \\ 0 &{} -1 \end{matrix}}\right] ,~ \sigma [10] \equiv X \equiv \left[ {\begin{matrix} 0 &{} 1 \\ 1 &{} 0 \end{matrix}}\right] ,~ \sigma [11] \equiv Y \equiv \left[ {\begin{matrix} 0 &{} -\dot{\imath } \\ \dot{\imath } &{} 0 \end{matrix}}\right] \end{aligned} $$

For n qubits, we define the set of “Pauli strings” \(\mathcal {\hat{P}}_n \,\triangleq \,\{P_1\otimes P_2\otimes \ldots \otimes P_n \mid P_j \in \{ I, X, Y, Z \} \}\). Inheriting the properties of Pauli matrices, Pauli strings are unitary, involutory and Hermitian. It is well-known that the scaled Pauli strings \(\{ \frac{1}{\sqrt{2}^n} \cdot P \mid P \in \mathcal {\hat{P}}_n \}\) form an orthonormal basis for \(2^n \times 2^n\) complex matrices [27]. Hence, we can decompose any \(2^n \times 2^n\) complex matrix M as \(M = \sum _{P \in \mathcal {\hat{P}}_n} \gamma _P \cdot P\) where the Pauli coefficient \(\gamma _P = \frac{1}{2^n} \textrm{Tr}(P^{\dagger } \cdot M)\).

In general, the coefficients \(\gamma _P\) are complex numbers, but for Hermitian matrices, they are real [18, 35].

Example 1

The matrix \(M = \left[ {\begin{matrix} 1&{}4+i\\ 4-i&{}-5\end{matrix}}\right] \) is Hermitian. We calculate the coefficients:

$$\begin{aligned} \frac{1}{2^1} \textrm{Tr}( I^{\dagger } M ) = -2 , {} & {} \frac{1}{2^1} \textrm{Tr}( Z^{\dagger } M ) = 3, {} & {} \frac{1}{2^1} \textrm{Tr}( X^{\dagger } M ) =4, {} & {} \frac{1}{2^1} \textrm{Tr}( Y^{\dagger } M ) = -1 \end{aligned}$$

It is straightforward to verify that these are M’s Pauli real coefficients:

$$ -2 I + 4 X -1 Y + 3 Z = -2 \cdot \left[ {\begin{matrix}1&{}0\\ 0&{}1\end{matrix}}\right] + 4 \cdot \left[ {\begin{matrix}0&{}1\\ 1&{}0\end{matrix}}\right] -1 \cdot \left[ {\begin{matrix}0&{}-i\\ i&{}0\end{matrix}}\right] + 3 \cdot \left[ {\begin{matrix}1&{}0\\ 0&{}-1\end{matrix}}\right] = \left[ {\begin{matrix} 1&{}4+i\\ 4-i&{}-5\end{matrix}}\right] . $$

   \(\square \)

Weighted Model Counting (WMC). In this work, we will encode the Pauli coefficients of specific matrices as weighted model counting: a sum of weights over all satisfying assignments of a boolean formula. We here formally describe WMC.

For boolean variables \(x,y \in \mathbb {B} = \left\{ 0,1 \right\} \), we define a literal as e.g. x and \(\overline{x}\) and write conjunctions of literals (cubes) as products, e.g., \(x\overline{y} = x \wedge \overline{y}\). A clause is a disjunction of literals, e.g., \(\overline{x} \vee y\). A formula in conjunctive normal form (CNF) is a conjunction of clauses.

Let \(F :\mathbb {B} ^{\vec x} \rightarrow \mathbb {B} \) be a propositional formula over boolean variables \(\vec x \in \mathbb {B} ^n\). We assign weights to literals using a weight function \(W:\left\{ \overline{x}, x \mid x\in \vec x \right\} \rightarrow \mathbb {R}\). Given an assignment \(\alpha \in \mathbb {B} ^{\vec x}\), let \(W(\alpha ) = \prod _{x\in \vec x} W(x = \alpha (x))\). We define weighted model counting [7, 10, 19, 21] as follows.

$$ MC_W(F) \,\triangleq \,\sum _{\alpha \in \mathbb {B} ^{\vec {x}}} F(\alpha )\cdot W(\alpha ) $$

Example 2

An example, consider a formula \( F = b \wedge c \) over \( \vec x = \left( a,b,c \right) \). There exist two satisfying assignments: \( \alpha _1 = abc \) and \( \alpha _2 = \overline{a} bc \). Suppose a weight function \( W \) is defined as follows: \( W(a) = -2 \), \( W(\overline{a}) = 3 \), \( W(b) = 1/2 \), \( W(\overline{b}) = 2 \), while \( c \) remains unbiased, i.e., \(W(c) = W(\overline{c}) = 1\). The weighted model counting for \( F \) with respect to \( W \) is computed as follows. \( MC_W(F) = F(abc) \cdot W( a bc) + F(\overline{a} bc) \cdot W(\overline{a} bc) = (-2 \cdot \frac{1}{2} \cdot 1) + (3 \cdot \frac{1}{2} \cdot 1) = \frac{1}{2}. \)    \(\square \)

3 Equivalence Checking Circuits in the Pauli Basis

In this section, we introduce (exact) equivalence checking [2, 4, 8, 22, 52, 55, 56, 61] in Definition 1, the task we set out to solve. In this work, we will only consider circuits which consist of gates, and do not contain measurements (this is without loss of generality since measurements be deferred to the end of the circuit [38]).

Definition 1

Given two n-qubit circuits U and V where \(n\in \mathbb {N}^+\), U is equivalent to V, written \(U\equiv V\), if there exists a complex number c (the global phase [38]) such that for all input states \(|{\psi }\rangle \), we have \(U|{\psi }\rangle = cV|{\psi }\rangle \).

At first sight, one might think that Definition 1 requires iterating over all quantum states. However, although the n-qubit quantum state space is continuous, it is a complex vector space of dimension \(2^n\), so it suffices to only consider \(2^n\) basis vectors for proving U and V equivalent. In fact, the novice approach to equivalence checking is to decompose U and V in the standard basis; that is, to find U and V each by writing each of their individual gates in the standard basis and determining the full unitaries U and V by matrix multiplication, and finally checking whether the matrix entries of U equal those of V, modulo a uniform constant c. One could also perform such an approach when the individual gates in U and V are specified in a different basis, such as the Pauli basis (see Sect. 2), but this would have no a priori advantage over the use of the standard basis. Instead, we will use the following folklore result (for proof see e.g. [52]).

Theorem 1

Let U, V be two circuits on \(n \in \mathbb {N}^{+}\) qubits. Then U is equivalent to V if and only if the following condition holds (for notation \(P_j\) see Eq. 1):

For all \(j\in [n]\) and \(P\in \left\{ X, Z \right\} \), we have \(U P_j U^{\dagger } = V P_j V^{\dagger }\).

The main advantage of using Theorem 1 instead of directly computing the (matrix entries of the) unitaries U and V is that for Clifford gates \(G, G'\) it is computationally easy to update the Pauli coefficients of \(GP_j G^{\dagger }\) to those of \((GG') P_j (GG')^{\dagger } = G \left( G' P_j G'^{\dagger }\right) G^{\dagger }\). This feature forms the basis for efficient simulation of Clifford circuits and has lead to efficient Clifford circuit equivalence checking [52]. Here, we will include T gates, Toffoli, and Pauli rotation gates, enabling equivalence checking of universal quantum computing (lifting the hardness of equivalence checking to quantum analogs of \(\textsf{NP}\), see Sect. 1). Another advantage of Theorem 1 is that, since U is a unitary, \(U P_j U^{\dagger }\) is Hermitian, so that its Pauli coefficients are real numbers as noted in Sect. 2, relieving us from the need to use complex numbers.

figure a

Example 3

Choose \(V = S_1\) and \(U = T_1 T_1\). In order to determine whether \(U\equiv V\), we compute the Pauli coefficients of \(UXU^{\dagger }, UZU^{\dagger }, VXV^{\dagger }\) and \(VZV^{\dagger }\) as follows using Table 1. By Theorem 1, this implies that U and V are equivalent, which we verify by writing their unitaries in the standard basis as follows.

$$ U = S = \left[ {\begin{matrix} 1&{}0\\ 0 &{} i\end{matrix}}\right] , \quad V = T\cdot T = \left[ {\begin{matrix} 1&{}0\\ 0&{}\sqrt{i}\end{matrix}}\right] \cdot \left[ {\begin{matrix} 1&{}0\\ 0&{}\sqrt{i}\end{matrix}}\right] = \left[ {\begin{matrix} 1&{}0\\ 0&{}i\end{matrix}}\right] $$

Finally, we remark that represents both constructive (Y terms add up) as well as destructive interference (X terms cancel).

We will finish this section by explaining the intuition behind Theorem 1, by rephrasing its proof from [52]. The first step in the proof is to realize that Definition 1 is equivalent to the following in density matrix representation.

Lemma 1

Given two n-qubit circuits U and V where \(n\in \mathbb {N}^+\), U is equivalent to V iff for all n-qubit quantum states \(|{\varphi }\rangle \), we have \(U\vert \varphi \rangle \!\langle {\varphi }\vert U^{\dagger } = V \vert \varphi \rangle \!\langle {\varphi }\vert V^{\dagger }\).

Recall that for any unitary U, with \(|{\psi }\rangle = U|{\varphi }\rangle \), the corresponding operation on the density matrix \(\vert {\varphi }\rangle \!\langle {\varphi }\vert \) is conjugation, i.e., \(\vert {\psi }\rangle \!\langle {\psi }\vert = U \vert {\varphi }\rangle \!\langle {\varphi }\vert U^\dagger \). Density matrices are \(2^n \times 2^n\) Hermitian matrices and can thus be expressed as a (real-weighted) linear combination of Pauli strings. For this reason, we observe that if \(UPU^{\dagger } = VPV^{\dagger }\) for each Pauli string P, i.e. U and V coincide on all Pauli strings by conjugation, then U and V must also coincide on all density matrices by conjugation, and thus they are equivalent by Lemma 1.

The final step in proving Theorem 1 is to realize that for a unitary matrix, the conjugation action is completely determined by fixing its conjugation action on only all \(X_j\) and \(Z_j\) for \(j\in [n]\). This insight relies on two parts: First, each Pauli string can be written as the product of \(X_j\) and \(Z_j\) modulo a factor \(\in \{\pm 1, \pm i\}\). Second, for a unitary M, we have \(M^{\dagger } M = I\), which implies that instead of first multiplying \(X_j\)s and \(Z_j\)s to construct a Pauli string, followed by conjugation, one can first conjugate and subsequently multiply to arrive at the same result.Footnote 3 For example, \(M X_j M^{\dagger } \cdot M Z_j M^{\dagger } = M X_j IZ_j M^{\dagger } = M X_j Z_j M^{\dagger }\).

We observe that in Table 1, the last two non-Clifford gates yield a linear combination of Pauli strings [34] for each Pauli string (matrix). This potentially causes an explosion of the number of Pauli strings when conjugating multiple non-Clifford gates. To handle this, we will exploit the strength of model counters in Sect. 4 by representing Pauli strings \({\hat{P}} \) as satisfying assignments which are weighted by the coefficient \(\gamma _{\hat{P}} \), as explained next in Sect. 4.

4 Encoding Quantum Circuit Equivalence in SAT

The previous section Sect. 3, centered around Theorem 1, explained that equivalence checking can be done by conjugating Pauli strings with unitaries, and that the required calculations for this approach are the same as in simulation of quantum circuits using a density matrix representation of the quantum state. In this section, we show how we reduce equivalence checking of universal quantum circuits to weighted model counting, which is formalized in Corollary 1 below. Our approach is based on the \(\mathcal O(n +m)\)-length encoding for quantum circuit simulation provided in [34]. Finally, our encoding in this work extends [34] with Toffoli gates. For the rest of the paper, we use P for an unweighted Pauli string and we use \(\textbf{P}\) for a summation of weighted Pauli strings, e.g. \(\frac{1}{\sqrt{2}}X + \frac{1}{\sqrt{2}}Y\).

Table 1. Lookup table for conjugating Pauli gates by Clifford+T+\(R_X\) gates. The subscripts “c” and “t” stand for “control” and “target”. Adapted from [34].
Full size table

To simplify notation, we will solve a rephrased version of the equivalence checking problem from Definition 1 in Sect. 3: to check whether a unitary A is equivalent to the identity unitary \(I\), which leaves every input unchanged. By choosing \(A \,\triangleq \,V^{\dagger }U\), we see that \(U \equiv V\) precisely if \(A \equiv I\). If U and V consist of gates \(U = (U_0, U_1, \dots , U_{m-1})\) and \(V = (V_0, V_1, \dots , V_{\ell -1})\) for \(m, \ell \in \mathbb {N}^+\), then a circuit for A is given as the \(m + \ell \) gates \(A = (U_0, U_1, \dots , U_{m-1}, V_{\ell -1}^{\dagger }, V_{\ell - 2}^{\dagger }, \dots , V_{0}^{\dagger })\). Following Theorem 1, our task will be as follows: Given a circuit \(A =\left( G^0, \dots , G^{m-1} \right) \in \left\{ H_j,S_j, CZ_{jk}, T_j, \text {Toffoli}_{jkl}, R_X(\theta )_j, \dots \mid j,k,l \in [n] \right\} ^m\), we need to obtain \(\textbf{P}^m = A \textbf{P}^0 A^{\dagger }\) from an initial \(\textbf{P}^0\in \{+X_i, +Z_i \mid i\in [n]\}\), showing that \(\textbf{P}^m = \textbf{P}^0\). Since \({\textbf {P}}^0\) is a Pauli string and thus Hermitian, so is \({\textbf {P}}^m\). Our approach is to construct a boolean formula whose weighted model counts represent the terms in the Pauli decomposition of \({\textbf {P}}^m\).

4.1 Encoding Pauli Coefficients as Weighted Model Counts

We first explain the encoding for circuit simulation from [34], where we encode the real-weighted sum of Pauli operators \(\textbf{P}\) and the update rules of the circuit A as weighted boolean formulas. We start with the simplest case—a Pauli string, then consider how to encode a single summand, i.e., a single weighted Pauli operator, and in the end extend this to a weighted sum of Pauli operators.

Given a Pauli string \(P=\bigotimes _{i\in [n]}\sigma [a_i, b_i]\) with \(a_i, b_i \in \{0,1\}\), the corresponding encoding is denoted as \(F_P\), which is the boolean formula which only has \(\{x_1\leftarrow a_1,\cdots , x_n\leftarrow a_n, z_1\leftarrow b_1,\cdots , z_n\leftarrow b_n\}\) as satisfying assignment, for example \(F_{Z \otimes X} = F_{\sigma [01]\otimes \sigma [10]} = \overline{x}_1 x_2 z_1 \overline{z}_2\). When it comes to weighted Pauli string, although the weights are never imaginary in case of a Hermitian matrix, they can still have a ± sign. A weighted Pauli operator can be therefore encoded by \(2n+1\) boolean variables: two bits \(x_i,z_i\) for each of the n Pauli matrices and one sign bit r, such that \(\textbf{P}=(-1)^r\sigma [x_{1},z_{1}] \otimes \ldots \otimes \sigma [x_{n},z_{n}]\). For example, consider boolean formula \(F_{\textbf{P}} = r\overline{x}_1 z_1x_2z_2\) where \(\textbf{P} = -Z\otimes Y\). Its one satisfying assignment is \(\{r\leftarrow 1, x_1 \leftarrow 0, z_1\leftarrow 1, x_2 \leftarrow 1, z_2\leftarrow 1\} \equiv -Z\otimes Y\). We later introduce weights \(W(r) = -1\) and \(W(\overline{r}) = 1\) to interpret the sign. So for a formula \(F(x_1,z_1, \dots ,x_n,z_n,r)\), we let the satisfying assignment represent a set (linear combination) of Pauli strings. The base case is the formula \(F_{\textbf{P}^0} = F_{P}\) for a Pauli string \(P \in \{X_j, Z_j \mid j \in [n]\}\). Next, we need to encode how sums of Pauli operators evolve when conjugating with the gates of the circuit, one by one. For this, our encoding duplicates the variables for all m gates (each time step) as follows (which is similar to encodings for bounded model checking [6]).

$$\begin{aligned} \vec w^t = \{x_j^t, z_j^t, r^t \mid j\in [n]\}\text { for } t\in \left\{ 0,1, \dots , m \right\} \text { and } \vec v^t = \bigcup _{i\in [t] \cup \{0\}} \vec w^i. \end{aligned}$$
(2)

For example, \(\textbf{P}^0 = X_1\) is encoded as \(\overline{r}^0 x_1^0 \overline{z}^0_1 x_2^0 \overline{z}^0_2 \dots x_n^0 \overline{z}^0_n \). Also, the satisfying assignments of a boolean formula \(F_A(\vec v^m)\) projected to variables \(\vec w^t\) represent the sum of Pauli operators after conjugating the initial t gates \(G_0, G_1, \dots , G_{t-1}\) of the circuit A, written:

$$ F_A(\vec v^m)[\vec w^t] = \sum _{\alpha \in \left\{ 0,1 \right\} ^{\vec v^m}} F_A(\alpha ) \cdot (-1)^{\alpha (r^t)}\,\cdot \bigotimes _{j\in [n]}\sigma [\alpha (x^t_{j}),\alpha (z^t_{j})] $$

The next question is how to encode gate semantics, i.e., define a constraint to get \(\textbf{P}^1\) by conjugating gate \(G^0\) to \(\textbf{P}^0\), etc. Note that since \(\textbf{P}^0 \in \left\{ X_j, Z_j \mid j \in [n] \right\} \) consists of a sum of only one Pauli operator. For Clifford circuits C, there will only be a single satisfying assignment \(\alpha \) for all time steps \(t\in [m]\), since e.g. \(H X H^\dag = Z\) (and not e.g. \(Z + Y\)). Non-Clifford gates, like T or Toffoli, will add satisfying assignments representing summands with different weights (e.g. sums of accumulated weights of \(\nicefrac 1{\sqrt{2}}\) for the T gate as discussed above). To encode these weights, we introduce new variables \(u^t\), but only at time steps t with a T gate (i.e., \(G^t= T\)).

When a gate \(T_j\) is performed and there is a satisfying assignment with \(x^t_{j}=1\), it means that we are conjugating a T gate on the j-th qubit set to \(\pm X\) or \(\pm Y\) and the result should be either \(TXT^\dag = \frac{1}{\sqrt{2}}(X+Y)\) or \(TYT^\dag =\frac{1}{\sqrt{2}}(Y-X)\) (modulo sign). To achieve this the encoding should let \(z^{t+1}\) unconstrained and set \(u^t\Leftrightarrow x^t_j\). Accordingly, we set the weights \(W(u^t) = \frac{1}{\sqrt{2}}\) and \(W(\overline{u}^t) = 1\). Table 2 illustrates how the boolean variables \(\vec w^t\) and \(\vec w^{t+1}\) relate for a T gate (derived by computing \(T P T^\dagger \) for Pauli gate P).

The encoding of gate semantics can be derived similarly. For example the boolean constraint for \(H^t_j\) follows from Table 1 and is given by

$$ F_{H^t_j}(\vec w^{t}, \vec w^{t+1}) \,\triangleq \,r^{t+1} \Leftrightarrow r^{t} \oplus x^t_{j}z^t_{j} \wedge ~ z^{t+1}_{j} \Leftrightarrow x^t_{j} \wedge ~ x^{t+1}_{j} \Leftrightarrow z^t_{j} $$

Here we omit additional constraints \(a^{t+1} \Leftrightarrow a^t\) for all unconstrained time-step-\(t+1\) variables a, i.e., for \(a = x^{t+1}_l,z_l^{t+1}\) with \(l \ne j\). Similarly, by abbreviating \(F_{G^t}(\vec w^t, \vec w^{t+1})\) as \(G^t\), the encoding for other Clifford+T gates are as follows:

$$\begin{aligned}\ S^t_j \,\triangleq \,~& r^{t+1} \Leftrightarrow r^{t} \oplus x^t_{j}z^t_{j} ~\wedge ~ z^{t+1}_{j} \Leftrightarrow x^t_{j} \oplus z^t_{j}, \\ \ T^t_j \,\triangleq \,~& r^{t+1} \Leftrightarrow r^t \oplus x^t_{j} z^t_{j} \lnot z^{t+1}_{j} ~\wedge ~ x^{t+1}_{j} \Longleftrightarrow x^t_{j} ~\wedge ~ x^{t}_{j} \vee (z^{t+1}_{j} \Leftrightarrow z^t_{j}) ~\wedge ~ u^t \Leftrightarrow x^t_{j},\\ CZ^t_{jk} \,\triangleq \,~&r^{t+1} \Leftrightarrow r^{t} \oplus x^t_{j}x^t_{k} (z^t_{k} \oplus z^t_{j}) ~\wedge ~ z^{t+1}_{k} \Leftrightarrow z^t_{j} \oplus x^t_{k} ~\wedge ~ z^{t+1}_{j} \Leftrightarrow z^t_{k} \oplus x^t_{j}. \end{aligned}$$
Table 2. Boolean variables under the action of conjugating one T gate. Here we omit the sign \((-1)^{r^{t}}\) for all P and sign \((-1)^{r^{t+1}}\) for all \({TPT^\dag }\).
Full size table

To this end, we can inductively define boolean constraints for each time step as \(F_{\textbf{P}^{t}}(\vec v^{t}) = F_{\textbf{P}^0}(\vec w^0)\wedge \bigwedge _{i\in [t-1] \cup \{0\}} G^i (\vec w^{i}, \vec w^{i+1})\) for \(t\ge 1\), where \(G_i\) denotes the gate at time step i and \(F_{\textbf{P}^0}(\vec v^0)\) encodes \(\textbf{P}^0\).

Example 4

Reconsider the circuit \(U = T \cdot T\) from Example 3. Starting with \(\textbf{P}^0 = X\), the formulas are \(F_{\textbf{P}^0} = x^0_1 \overline{z^0_1} \overline{r^0}\), \(F_{\textbf{P}^1} = F_{\textbf{P}^0} \wedge F_{T^0_1}\), i.e.

$$\begin{aligned} F_{\textbf{P}^0} \wedge F_{T^0_1} = &~ x^0_1 \overline{z^0_1} \overline{r^0} ~\wedge ~ x^{1}_{1}\Leftrightarrow x^0_{1} ~\wedge ~ x^{0}_{1} \vee (z^{1}_{1} \Leftrightarrow z^0_{1}) ~\wedge ~ \\ &~r^{1} \Leftrightarrow r^0 \oplus x^0_{1} z^0_{1} \lnot z^{1}_{1} ~\wedge ~ u^0 \Leftrightarrow x^0_{1}, \end{aligned}$$

and similarly \(F_{\textbf{P}^2} = F_{\textbf{P}^1} \wedge F_{T^1_1}\).    \(\square \)

Formalizing the explanation above as induction over the gates proves Proposition 1, relating weighted model counting the Pauli coefficients (see Sect. 2).

Proposition 1

(WMC computes the Pauli coefficients). Let \(C_A =(G^0, \) \( \dots , G^{m-1})\) be an n-qubit circuit, \(A = G^0 \cdots G^{m-1}\) the corresponding unitary and \(\textbf{P}^0\) a Pauli string, so that the encoding of \(\textbf{P}^m \,\triangleq \,A \textbf{P}^0 A^\dag \) is given by \(F_{\textbf{P}^m} \,\triangleq \,F_{\textbf{P}^0} \wedge \bigwedge _{i=0}^{m-1} F_{G^i}\) with according weight function W. For any \(\textbf{P}^0 \in \{+X_j, +Z_j \mid j \in [n]\}\) and \(P \in \mathcal {\hat{P}}_n \), the weighted model count of \(F_{\textbf{P}^m}\wedge F_{P}\) equals the Pauli coefficient \(\gamma _{P}\) of \(\textbf{P}^m\). That is, \(MC_W(F_{\textbf{P}^m}\wedge F_{P}) = \frac{1}{2^n} \cdot \textrm{Tr}(P^{\dagger } \cdot A\textbf{P}^0 A^{\dagger })\) for all \(P \in \mathcal {\hat{P}}_n \).

We emphasize the necessity for using negative weights. For example, in Example 3, we have for \(\textbf{P}^0 = X\), where the terms X and \(-X\) cancel each other out, while the Y terms add up. This is why weighted model counting with negative weights is required; to reason about such constructive and destructive interference, ubiquitous to quantum computing.

Example 5

Following Example 4, we have the satisfying assignments for \(F_{\textbf{P}^0}\), \(F_{\textbf{P}^1}\) and \(F_{\textbf{P}^2}\) as:

$$\begin{aligned} SAT(F_{\textbf{P}^0}) &~= \{x^0_1\overline{z^0_1}\overline{r^0_1}\}, \\ SAT(F_{\textbf{P}^1}) &~= \{x^0_1\overline{z^0_1}\overline{r^0_1} \ x^1_1 \overline{z^1_1} \overline{r^1_1} \ u^0, \ x^0_1\overline{z^0_1}\overline{r^0_1} \ x^1_1z^1_1\overline{r^1_1} \ u^0\}, \\ SAT(F_{\textbf{P}^2}) &~= \{ x^0_1\overline{z^0_1}\overline{r^0_1} \ x^1_1 \overline{z^1_1} \overline{r^1_1} \ x^2_1 \overline{z^2_1} \overline{r^2_1} \ u^0 u^1, \ x^0_1\overline{z^0_1}\overline{r^0_1} \ x^1_1 \overline{z^1_1} \overline{r^1_1} \ x^2_1 {z^2_1} \overline{r^2_1} \ u^0 u^1, \\ & \qquad x^0_1\overline{z^0_1}\overline{r^0_1} \ x^1_1 {z^1_1} \overline{r^1_1} \ x^2_1 {z^2_1} \overline{r^2_1} \ u^0 u^1, \ x^0_1\overline{z^0_1}\overline{r^0_1} \ x^1_1 {z^1_1} \overline{r^1_1} \ x^2_1 \overline{z^2_1} {r^2_1} \ u^0 u^1\}, \end{aligned}$$

with the weight function \(W(r^2_1) = -1\), \(W(\overline{r^2_1}) = 1\), \(W(u^0) = W(u^1) = \frac{1}{\sqrt{2}}\) and \(W(\overline{u^0}) = W(\overline{u^1}) = 1\). Each of the satisfying assignments corresponds to a term in the Pauli decomposition of \(\textbf{P}^2\), which we recall from Example 3 to be

$$\begin{aligned} \textbf{P}^2 = \frac{1}{2}X+\frac{1}{2}Y+\frac{1}{2}Y-\frac{1}{2}X = (\frac{1}{2} - \frac{1}{2}) X + (\frac{1}{2} + \frac{1}{2}) Y = Y . \end{aligned}$$
(3)

For example, the term \(-\frac{1}{2} X\) is encoded by \(x^0_1\overline{z^0_1}\overline{r^0_1} \ x^1_1 {z^1_1} \overline{r^1_1} \ x^2_1 \overline{z^2_1} {r^2_1} \ u^0 u^1\) because it contains \(x^2_1 \overline{z^2_1}\) (corresponding to X) and its weight is \(W(r^2_1) \cdot W(u^0) \cdot W(u^1) = (-1) \cdot \frac{1}{\sqrt{2}} \cdot \frac{1}{\sqrt{2}} = -\frac{1}{2}\). We verify that the constructive interference of the Y terms in (3) (i.e. they add up) results in an aggregate Pauli coefficient \(\gamma _Y\) of \(\textbf{P}^2\) of 1:

$$ MC_W(F_{\textbf{P}^2} \wedge F_Y) = \frac{1}{\sqrt{2}}\cdot \frac{1}{\sqrt{2}} + \frac{1}{\sqrt{2}}\cdot \frac{1}{\sqrt{2}} = 1 = \frac{1}{2}\textrm{Tr}(Y\cdot \textbf{P}^2). $$

Similarly, we verify that destructive interference of the X terms in (3) (i.e. they cancel) results in the coefficient \(\gamma _X\) being 0:

$$ MC_W(F_{\textbf{P}^2} \wedge F_X) = \frac{1}{\sqrt{2}}\cdot \frac{1}{\sqrt{2}} - \frac{1}{\sqrt{2}}\cdot \frac{1}{\sqrt{2}} = 0 = \frac{1}{2}\textrm{Tr}(X\cdot \textbf{P}^2). $$

   \(\square \)

Toffoli Gate. Similar to the way gate encodings of other non-Clifford gates were derived, we can encode the Toffoli gate. To this end, we brute forced the Toffoli gate behavior in the Pauli domain. To keep things readable, we will only present a lookup table in the Pauli basis in Table 3, like Table 1. The corresponding boolean constraint can easily be derived. To subsequently obtain a minimal (weighted) CNF formula, we applied the Quine-McCluskey algorithm [33, 44].

Table 3. An partial lookup table for the Toffoli gate for in/output Pauli operators P and Q. The extended version of this paper [35] includes the full table.
Full size table

4.2 WMC-Based Algorithm for Equivalence Checking

The previous subsection explains how to encode the Pauli coefficients of \(A P A^{\dagger }\), where A is a unitary and P a Pauli string, in a boolean formula together with a weight function. We here connect this encoding to Theorem 1, which expresses that determining whether a unitary A is equivalent to the identity circuit can be done by checking if \(APA^{\dagger } {\mathop {=}\limits ^{?}} P\) for Pauli strings \(P \in \{X_j, Z_j \mid j \in [n]\}\). We use the following lemma, which expresses that for any unitary A and Pauli string P, the P-Pauli coefficient of \(AP A^{\dagger }\) can only become 1 if \(A P A^{\dagger }\) equals P.

Lemma 2

Let A be a unitary and \(P \in \mathcal {\hat{P}}_n \) be a Pauli string. Then \(AP A^{\dagger } = P\) if and only if \(\frac{1}{2^n} \textrm{Tr}(A P A^{\dagger } \cdot P) = 1\).

Proof

If \(AP_j A^{\dagger } = P_j\), then \(\textrm{Tr}(A P_j A^{\dagger } \cdot P_j) = \textrm{Tr}(P_j \cdot P_j) = \textrm{Tr}(I^{\otimes n}) = 2^n\). For the converse direction, we observe that \(\textrm{Tr}(A P_j A^{\dagger } \cdot P_j)\) is the Frobenius inner product \(\langle U, V \rangle \,\triangleq \,\textrm{Tr}(U^{\dagger } V)\) for \(U \,\triangleq \,A P_j A^{\dagger }\) and \(V \,\triangleq \,P_j\). It now follows from the Cauchy-Schwarz inequality \(|\langle U, V\rangle |^2 \le \langle U,U \rangle \cdot \langle V ,V \rangle \) that

$$\begin{aligned} |\textrm{Tr}(A P_j A^{\dagger } \cdot P_j)|^2 &\le \textrm{Tr}((A P_j A^{\dagger })^{\dagger } \cdot A P_j A^{\dagger }) \cdot \textrm{Tr}(P_j^{\dagger } \cdot P_j) \\ &= \textrm{Tr}(A P_j A^{\dagger } \cdot A P_j A^{\dagger }) \cdot \textrm{Tr}(P_j^{\dagger } \cdot P_j) & \\ &= \textrm{Tr}(A P_j \cdot P_j A^{\dagger }) \cdot \textrm{Tr}(I^{\otimes n}) \qquad \quad \text {( }A \text { and } P_j \text { are unitary)} \\ &= \textrm{Tr}(A A^{\dagger }) \cdot \textrm{Tr}(I^{\otimes n}) \qquad \qquad \qquad \text {( } P^2 = I \text { for all } P \in \mathcal {\hat{P}}_n \text { )} \\ &= \textrm{Tr}(I^{\otimes n}) \cdot \textrm{Tr}(I^{\otimes n})= 2^n \cdot 2^n = 4^n \qquad \qquad \text {( }A \text { is unitary)} \end{aligned}$$

and therefore \(|\textrm{Tr}(A P_j A^{\dagger } \cdot P_j)| \le 2^n\). Since \(\textrm{Tr}(A P_j A^{\dagger } \cdot P_j) = 2^n\) by assumption, the Cauchy-Schwarz inequality is tight, which only happens if \(U=AP_jA^{\dagger }\) and \(V=P_j\) are linearly dependent. Thus, there exists a complex number \(\lambda \) such that \(AP_jA^{\dagger } = \lambda P_j\). Substituting this expression in \(\textrm{Tr}(AP_jA^{\dagger } \cdot P_j)\) yields \(\textrm{Tr}(\lambda P_j \cdot P_j) = \lambda \cdot \textrm{Tr}(I^{\otimes n}) = \lambda 2^n\), hence \(\lambda = 1\) and \(A P_j A^{\dagger } = P_j\).   \(\square \)

Combining Lemma 2 and Proposition 1 with Theorem 1 yields Corollary 1 below, which in turn implies correctness of Algorithm 1 which reduces equivalence checking to WMC.

Corollary 1

Let A be an n-qubit circuit with m gates and \(P \in \left\{ X_j, Z_j \mid j\in [n] \right\} \), which are encoded by \(F_{A}\) and \(F_{P}\) respectively, with according weight function W. We have \(A \equiv I\) if and only if \(MC_W(F_{P}(\vec {w}^0) \wedge F_A\left( \vec {v}^t \right) \wedge F_{P}(\vec {w}^m)) = 1\) for all \(P \in \left\{ X_j, Z_j \mid j\in [n] \right\} \), where \(\vec {w}^{t+1}\) are boolean variables encoding the quantum state in circuit A after the t-th gate of A (\(0\le t \le m -1\)) and \(\vec {v}^t = \bigcup _{t \in [m] \cup \{0\}} \vec {w}^t\) as defined in Eq. (2).

figure d

Example 6

Consider \(A = V^{\dagger } U\) where =\(U = (T, T)\) and \(V = (S)\) as in Example 3. We show how to reduce the equivalence check \(A{\mathop {\equiv }\limits ^{?}} I\) to weighted model counting. First, we encode the check \(A X A^{\dagger } {\mathop {=}\limits ^{?}} X\) using \(F_1 \,\triangleq \,F_{AXA^\dagger }\wedge F_{X}\):

figure e

The satisfying assignments of \(F_1\) are

$$\begin{aligned} SAT{(F_1)} = \{ & x^1_1 \overline{z^1_1} \overline{r^1_1} \ x^2_1 {z^2_1} \overline{r^2_1} \ x^3_1 \overline{z^3_1} \overline{r^3_1} \ u^0 u^1, \ x^1_1 \overline{z^1_1} \overline{r^1_1} \ x^2_1 {z^2_1} \overline{r^2_1} \ x^3_1 \overline{z^3_1} \overline{r^3_1} \ u^0 u^1\}. \end{aligned}$$

so \(MC_W(F_1) = \sum _{\sigma \in SAT(F_1)}W(\sigma (r^3_1))W(\sigma (u^0))W(\sigma (u^1)) = \frac{1}{\sqrt{2}}\cdot \frac{1}{\sqrt{2}} + \frac{1}{\sqrt{2}}\cdot \frac{1}{\sqrt{2}} = 1\).

Now we turn to the check \(AZA^{\dagger } {\mathop {=}\limits ^{?}} Z\), obtaining the formula \(F_2 \,\triangleq \,\ F_{AZA^\dagger } \wedge F_Z\), where \(F_{AZA^\dagger }\) is the same formula from \(F_1\) and \(F_Z = \overline{x^3_1} z^3_1\). The satisfying assignments of \(F_2\) are \(SAT(F_2) = \{\overline{x^0_1} {z^0_1} \overline{r^0_1}\overline{x^1_1} {z^1_1} \overline{r^1_1} \overline{x^2_1} {z^2_1} \overline{r^2_1}\overline{x^3_1} {z^3_1} \overline{r^3_1} \overline{u^0} \overline{u^1}\}\), and \(MC_W(F_2)=W(\overline{r^3_1})W(\overline{u^0})W(\overline{u^1}) = 1\). Since both weighted model counts evaluate to 1, we conclude that \(A\equiv I\).    \(\square \)

5 Implementation: The ECMC Tool

We implemented our method in an open-source tool called ECMC, available at https://github.com/System-Verification-Lab/Quokka-Sharp. ECMC takes two quantum circuits in QASM format [13] as input. It encodes these circuits to a sequence of 2n weighted conjunctive normal form (CNF) formulas as explained in Sect. 4, and then uses the weighted model counter GPMC [50] to solve these constraints in parallel, terminating as soon as one returns a negative result. Here we set the number of parallel cores to be 16 as it is shown to be the optimal number of cores for our task.

We choose GPMC as it supports the negative weights in our encoding and performs the best among solvers with that capability in the model counting competition 2023 [21]. To demonstrate the effectiveness of our method, we conducted a set of broad experiments as discussed in the following.

We performed equivalence checking of quantum circuits comparing our method against the state-of-the-art tool QCEC  [8], which runs different algorithms and heuristics based on ZX calculus and decision diagrams (shorted as DD) in portfolio with 16 parallel threads [60]. Similar to ECMC, QCEC also terminates earlier when one thread returns “non-equivalent”. Since the ZX-calculus based method is still incomplete for universal quantum circuits, in the sense that it is only capable of proving equivalence, we use this tool under two settings: one is the default setting which uses DD and ZX calculus in portfolio; the other is to exclusively enable DD [8]. We use two families of circuits: (i) random Clifford+T circuits, which mimic hard problems arising in quantum chemistry [59] and quantum many-body physics [17]; (ii) all benchmarks from the public benchmark suite MQT Bench [43], which includes many important quantum algorithms like QAOA, VQE, QNN, Grover, etc. All experiments have been conducted on a 3.5 GHz M2 Machine with MacOS 13 and 16 GB RAM. We set the time limit to be 5 min (300 s) and include the time to read a QASM file, construct the weighted CNF and perform the model counting in all reported runtimes.

Fig. 1.
figure 1

Equivalence check of typical random Clifford+T circuits against their optimized circuits (equivalent cases, Fig 1 & Fig 2) and optimized circuits with one random gate missing (non-equivalent cases, Fig 3 & Fig 4). (Both vertical axes are on a logarithmic scale.)

Full size image

Results. First, to show the scalability of both methods on checking equivalence, we consider random circuits that resemble typical oracle implementations—random quantum circuits with varying qubits and depths, which comprise the CX, H, S, and T gates with appearing ratio 10%, 35%, 35%, 20% [41]. We use a ZX-calculus tool PyZX [28] to generate optimized circuits, to construct equivalent, yet very different, counterparts. To construct non-equivalent instances, we inject an error by removing one random gate from the corresponding optimized circuits. So by construction, we know the correct answer for all equivalence checking instances in advance. The resulting runtimes can be seen in Fig. 1.

Table 4. Results of verifying equivalence of circuits from MQT bench against optimized circuits. For cases within time limit, we give runtime (sec), while \(>300\) represents a timeout (5 min) and ✕ means that the result was ‘unknown’.
Full size table

In addition to random circuits, to test structural quantum circuits, we empirically evaluated our method on the MQTBench benchmark set [43]. We also generate the optimized circuits of the circuits from MQT-bench using PyZX [28]. To generate non-equivalent instances, three kinds of errors are injected into the optimized circuits: one with a random gate removed, one where a random CNOT gate is flipped, switching control and target qubits, and one where the phase of the angle of a random rotation gate is shifted. For the last error, since many optimizations on rotation gates involve phase shifts in the rotation angles, we consider two sizes of phase shift: one with the angle of a random rotation gate added by \(10^{-4}\), one with the angle added by \(10^{-7}\). We note that this experimental setup is stronger than the one used in [41], where only two errors are considered: bit flip and phase shift without giving the shifting scale. We present a representative subset of equivalence checking results in Table 4. The complete results can be found in the extended version of this paper [35]. The first three columns list the number of qubits n and gates |G| in original circuits, and the number of gates \(|G'|\) in optimized circuits. Then we give the runtime of the weighted model counting tool ECMC, the decision diagram-based QCEC (DD) and the default setting of QCEC respectively. For the non-equivalent cases, we show the flipped-CNOT and one-gate-missing error in Table 5. The first three columns are the same as Table 4 and then the performance of all three tools on CNOT flipped error and one-gate-missing error respectively. Finally, Table 6 shows the performance of phase shift errors, where Shift-\(10^{-4}\) (resp. Shift-\(10^{-4}\)) denotes adding \(10^{-4}\) (resp. \(10^{-7}\)) to the phase of a random rotation gate.

Table 5. Results of verifying non-equivalence of circuits from MQT bench against optimized circuits with flipped CNOT gate (Flipped) and one missing gate (1 Gate Missing). For cases within time limit, we give runtime (sec), while \(>300\) represents a timeout (5 min).
Full size table

Discussion. For random circuits, Fig. 1 shows that the runtime of ECMC exhibits a clear correlation with the size of the circuits. While QCEC and QCEC (DD) are very fast for small size circuits, for non-equivalent cases, both of them are less scalable and reach time limit much earlier than ECMC. For the equivalent cases, QCEC benefits from ZX calculus and outperforms the other two methods. We suspect that QCEC (DD) shows poor performance when solving random circuits because these circuits don’t contain the structure found in quantum algorithms, which decision diagrams can typically exploit.

Table 6. Results of verifying non-equivalence of circuits from MQT bench against optimized circuits with \(10^{-4}\) size and \(10^{-7}\) size phase shift in one random rotation gate. For cases within time limit, we give runtime (sec), while \(>300\) represents a timeout (5 min), “wrong” a wrong result and ✕ that the results was ‘unknown’.
Full size table

When considering structural quantum circuits, the results vary between equivalent and non-equivalent instances. For equivalent instances, QCEC (DD) significantly surpasses ECMC on Grover, QFT and QNN, primarily due to the decision diagram-based method’s proficiency in handling circuits featuring repeated structures and oracles. While for those circuits featuring a large number of rotation gates with various rotation angles, like graphstate and wstate, ECMC demonstrates clear advantages. Moreover, the default QCEC is much faster than QCEC (DD) on all cases while it reports “no information” for many cases as ZX calculus method and decision diagram method give different answers.

For non-equivalent instances, since ECMC can terminate when a single out of 2n WMC calls returns a negative result, it shows better performance than checking equivalence. For example, in the case of QPE, where both tools face time constraints when checking equivalent instances, ECMC can efficiently demonstrate non-equivalence and resolve the majority of cases within the time limit, while both QCEC and QCEC (DD) still get timeout in most instances.

In all instances, ECMC outperforms both QCEC and QCEC (DD) on graph state and wstate, each featuring many rotation gates. When dealing with rotation gates, decision diagrams might suffer from numerical instability [39, 41], as can be clearly observed in Table 6 for the instances with errors in the phase shift, where both QCEC and QCEC (DD) get wrong results for many benchmarks. In contrast, the WMC approach—also numerical in nature—iteratively computes a sum of products, which we think avoids numerical instability. Table 6 also demonstrates this point as ECMC yields the correct answer for most benchmarks with \(10^{-4}\) and \(10^{-7}\)-size error. In contrast, the default QCEC gives no answer for a large amount of cases.

6 Related Work

Bauer et al. [3] verify quantum programs by encoding the verification problem in SMT, using an undecidable theory of nonlinear real arithmetic with trigonometric expressions. An SMT theory for quantum computing was proposed in [11]. Berent et al. [4] realize a Clifford circuit simulator and equivalence checker based on a SAT encoding. The equivalence checker was superseded by the deterministic polynomial-time algorithm proposed and implemented in [52]. Using weighted model counting, universal quantum circuit simulation is realized in [34], which we extend by providing encodings for the CZ and Toffoli gates and which we apply to circuit equivalence checking according to the approach of [52]. Amy [2] uses path integrals to check equivalence of circuits, which is complete for Clifford circuits and can prove equivalence of Clifford+T and Clifford+R circuits.

Yu and Palsberg [62] use an abstract interpretation to simulate quantum circuits. Abstraqt [5] improves upon this by using the stabilizer basis. SAT solvers have proven successful in quantum compilation [53], e.g., for reversible simulation of circuits [58] and optimizing space requirements of quantum circuits [36, 45].

The ZX calculus [12] offers a diagrammatic approach to manipulate and analyze quantum circuits. A circuit is almost trivially expressible as a diagram, but the diagram language is more powerful and circuit extraction is consequently #\(\textsf{P}\)-complete [15]. It has proven enormously successful in applications from equivalence checking [41, 42], to circuit optimization [28] and simulation [29].

Decision diagrams [37] have been used for simulating quantum circuits, checking their equivalence [8] and synthesis [64]. Jimenez et al. use bisimulation for circuit reduction, reducing simulation time compared to DDs in some cases [26].

7 Conclusions

We have shown circuit equivalence checking reduces to weighted model counting by considering quantum states in the Pauli basis, which allows for an efficient reduction of the equivalence checking problem to weighted model counting. We extended a linear-length encoding with the three-qubit Toffoli gate, so that most common non-Clifford gates are supported (previously the T, phase shift and rotation gates were already supported).

Given two n-qubit quantum circuits, their equivalence (up to global phase) can be decided by 2n calls to a weighted model counter, each with an encoding that is linear in the circuit size. Our open source implementation demonstrates that this technique is competitive to state-of-the-art methods based on a combination of decision diagrams and ZX calculus. This result demonstrates the strength of classical reasoning tools can transfer to the realm of quantum computing, despite the general ‘quantum-hardness’ of these problems. In future work, we plan to extract diagnostics for non-equivalent circuits from the satisfying assignments of the model counter.