1 Introduction↩︎

Quantum computing promises computational advantages in domains such as cryptography [1], finance [2], optimization [3], and the simulation of quantum physics [4] and chemistry [5]. For quantum computers to scale to practically useful tasks, various design challenges need to be solved, such as quantum circuit compilation, performance comparisons of quantum algorithms, and evaluating and mitigating the effects of noise. An important tool to tackle these challenges is the simulation of quantum circuits on classical computers [6]–[8].

Among the algorithms and data structures for quantum circuit simulation are decision diagrams (DD) [9], which have been proven to be successful in the simulation of quantum circuits and related tasks [10]–[19]. In quantum computing, states and transformations of states can be represented as pseudo-Boolean functions \(f: \{0, 1\}^n \rightarrow \mathbb{C}\), which are exponentially large in the number of qubits. Decision diagrams are a compact alternative that represents pseudo-Boolean functions as directed graphs. The nodes of these graphs represent sub-functions \(f_{\vec{a}} : \vec{x} \mapsto f(\vec{a}, \vec{x})\), for \(\vec{a} \in \{0, 1\}^m, m \leq n\). DDs achieve succinctness due to canonicity: by converting each node to a normal form, two nodes representing the same function are merged, avoiding storing multiple sub-functions if they are identical. DDs that offer stronger node merging also exist. For example, in quantum multiple-value decision diagrams (QMDDs [11], i.e., factored edge-valued DDs [20] with complex edge values), two nodes can be merged if they represent the same sub-function up to some scalar multiple.

Local-Invertible Map Decision Diagrams (LIMDDs) enable even stronger node merging: two nodes are merged if they represent the same sub-function up to some local invertible map (LIM), Kronecker products of invertible \(2\times 2\) matrices. Pauli-LIMDDs, where the LIMs consist of matrices from the Pauli group (a set of matrices well known in quantum physics), were proven to be exponentially faster for quantum-circuit simulation than QMDDs and Algebraic Decision Diagrams [21], [22], and can, for specific circuits, also outperform common techniques such as matrix product states and Clifford+T simulation [16], [23].

Practical implementations of LIMDDs have so far failed to fully realize their theoretical advantages, showcasing the difficulty of their (correct) implementation. The first cause is the complicated \(\mathcal{O}(n^3)\) run time merge rule algorithm, solving a group-theoretic task: in practice, a fast implementation is needed to outperform the \(\mathcal{O}(1)\) run time merging for QMDDs. In addition to this, canonicity, the crucial property for succinct and efficiently updatable decision diagrams, requires careful handling of many (edge) cases. To make matters worse, bugs in the canonicity implementation are hard to catch, as a non-canonical LIMDD can represent the correct quantum state (and thus pass all typical software tests) while still not being optimal in node number. Indeed, neither of the two existing proof-of-principle LIMDD implementations is fully canonical. First, MQT-LIMDD [24], built as an extension of the QMDD-simulator MQT-DDSIM [13], does not achieve minimal size in all cases. Second, LimTDD [18], which implements a generalization of Pauli-LIMDDs, deliberately does not include the canonicity property, voiding any guarantees against (exponential) blow-ups.

In this work, we present a novel algorithm for the Pauli-LIMDD merge rule, implemented as an open-source stand-alone library PauliLim [25]. The algorithm has less overhead and achieves a run time improvement from \(\mathcal{O}(n^3)\) to \(\mathcal{O}(n^2)\) for \(n\)-qubit nodes with a single child (one of the cases in which DDs achieve a succinctness advantage). We integrate this new algorithm into QolDDer, an open-source [26] from-scratch implementation of Pauli-LIMDDs. In an empirical evaluation, we compare QolDDer against the two previous LIMDD implementations [18], [24], as well as against the well-established QMDD implementation MQT-DDSIM [13]. These empirical results corroborate that LIMDDs are able to achieve exponential advantages over QMDDs and underline that LIMDDs require canonicity to realize these advantages. Our work will enable further application and development of LIMDD variants, not only for quantum design tasks but also for analysis of linear-algebra-based systems in general.

2 Preliminaries↩︎

2.1 Quantum computing↩︎

We briefly introduce quantum computing; see [27] for a full overview. The state of a single quantum bit (qubit) is a complex-valued 2-vector of unit norm: either \(\ket{0} = \begin{bmatrix} 1 & 0 \end{bmatrix}^{\intercal}\) or \(\ket{1} = \begin{bmatrix} 0 & 1 \end{bmatrix}^{\intercal}\), or a so-called superposition of both, \(\begin{bmatrix} \alpha_0 & \alpha_1 \end{bmatrix}^{\intercal} = \alpha_0\ket{0} + \alpha_1 \ket{1}\) where \(\alpha_0, \alpha_1\) are complex numbers satisfying \(|\alpha_0|^2 + |\alpha_1|^2 = 1\) and where \(|.|\) denotes the modulus of a complex number. The joint state of two qubits, \(\ket{\psi_A}\) and \(\ket{\psi_B}\), is obtained from the Kronecker product of their individual states; \(\begin{bmatrix} \alpha_0 & \alpha_1 \end{bmatrix}^\intercal\otimes \begin{bmatrix} \beta_0 & \beta_1 \end{bmatrix}^\intercal= \begin{bmatrix} \alpha_0 \beta_0 & ~\alpha_0 \beta_1 & ~\alpha_1 \beta_0 & ~\alpha_1 \beta_1 \end{bmatrix}^\intercal\). In general, the state of \(n\) qubits can be represented by a vector of unit norm, \(\ket{\psi} \in \mathbb{C}^{2^n}\), and can be written in the Shannon decomposition \(\ket{0} \otimes \ket{\psi_0} + \ket{1} \otimes \ket{\psi_1}\), where \(\ket{\psi_0}, \ket{\psi_1}\) are (unnormalized) \((n-1)\)-qubit states.

Quantum gates are (reversible) operations that map \(n\)-qubit states to \(n\)-qubit states, represented as \(2^n \times 2^n\) unitary matrices. Several common gates are listed below. Two sets of gates that are particularly relevant to this work are the Pauli gates \(\{I, X, Y, Z{=}P(\pi)\}\), and the Clifford gates, which are generated from \(\{H, P(\tfrac{\pi}{2}), \mathit{CX}\}\) under matrix multiplication and Kronecker products. \[\begin{align} {I}{=}\begin{bmatrix*}[r] 1 & 0\\ 0 & 1 \end{bmatrix*}~ X{=}\begin{bmatrix*}[r] 0 & 1\\ 1 & 0 \end{bmatrix*}~ Y{=}\begin{bmatrix*}[r] 0 & \text{-}i\\ i & 0 \end{bmatrix*}~ H{=}\frac{1}{\sqrt{2}} \begin{bmatrix*}[r] 1 & 1\\ 1 & \text{-}1 \end{bmatrix*}~ P(\theta){=}\begin{bmatrix*}[r] 1 & 0\\ 0 & e^{i\theta} \end{bmatrix*}~ \mathit{CX}{=}\begin{bmatrix} 1&0&0&0\\0&1&0&0\\0&0&0&1\\0&0&1&0\end{bmatrix} \end{align}\] The quantum state after the application of a gate is obtained through matrix-vector multiplication. When a gate is applied to a state with a greater number of qubits, the gate matrix is extended using Kronecker products with \(I\), so on a 3-qubit state \(\ket{\psi}\), applying \(X\) on the second qubit can be computed as \(I\otimes X \otimes I\ket{\psi}\).

A quantum circuit is a sequence of gates (and measurements) performed on a register of qubits. A straightforward way to simulate quantum circuits is to compute the final state vector by performing matrix-vector multiplication for all gates in the circuit, starting from the initial state. Desired information, such as measurement outcomes, can then be extracted from this final state. The limiting factor of this approach is the exponential (\(2^n\)) size of \(n\)-qubit state vectors.

We finish with notions from Pauli stabilizer algebra [27], [28]. By an \(n\)-qubit Pauli word \(P\in \mathrm{\small Pauli}_n\), we mean a length-\(n\) Kronecker product of the Pauli matrices, e.g. \(X\otimes Y\otimes Z\) or \(Z\otimes I\), abbreviated by omitting ‘\(\otimes\)’ as \(XYZ\) and \(ZI\). We will denote the \(2^n \times 2^n\) identity matrix \(I_2^{\otimes n}\) as \(I_{n}\) and write \(\mathbb{I}\) when \(n\) is clear from the context. For this work, a useful representation of a Pauli gate is a binary tuple \((x, z) \in \{0, 1\}^2\), defined as \(I\simeq (0, 0), Z\simeq (0, 1), X \simeq (1, 0), Y \simeq (1, 1)\). We extend this to Pauli words using a reverse order than is common [28], defining the binary vector of \(P\) as \(\textsf{revbin}\mathopen{}\mathclose{\left(P_1 \otimes P_2 \otimes \dots P_n}\right) = (x_n, x_{n-1}, \dots, x_1, z_{n}, z_{n-1}, \dots, z_1)\). For example, \(\textsf{revbin}\mathopen{}\mathclose{\left(YZ}\right) = (0, 1, 1, 1)\).

The \(n\)-qubit Pauli group is the closure of \(n\)-qubit Pauli words under matrix multiplication. A stabilizer group of a state \(\ket{\phi}\) is the commutative subgroup \(\text{Stab}(\ket{\phi}) = \{\lambda P \mid \lambda P \ket{\phi} = \ket{\phi}, P\in\mathrm{\small Pauli}_n, \lambda\in\{\pm 1\}\}\) of the Pauli group. A generating set for a stabilizer group is a set \(S\) of signed Pauli words whose closure under multiplication, \(\langle S\rangle\), is that stabilizer group. An example is \(\{I, -Z\} = \langle -Z\rangle\), with generating set \(\{-Z\}\) as \((-Z)^2 = I\). For a signed Pauli word \((-1)^s \cdot P\) with \(s\in \{0, 1\}\), we define \(\textsf{revbin}\mathopen{}\mathclose{\left(sP}\right)\) as the concatenation of \(\textsf{revbin}\mathopen{}\mathclose{\left(P}\right)\) and the bit \(s\), e.g., \(\textsf{revbin}\mathopen{}\mathclose{\left(-YZ}\right) = (0, 1, 1, 1, 1)\). Stabilizer generating sets \(\{g_1, g_2, \dots, g_k\}\) can be written as binary matrix where the \(r\)-th row is the binary vector \(\textsf{revbin}\mathopen{}\mathclose{\left(g_r}\right)\). Such matrices can be brought in row-echelon form (REF) using Gaussian elimination, an elementary linear algebra algorithm, by a slight adjustment of vector addition (see the ‘rowsum’ algorithm in [28]). We will say that a stabilizer generating set is in REF if its binary representation is.

2.2 Decision diagrams for quantum computing↩︎

Decision diagrams (DDs) [9] are a data structure that heuristically compresses the representation of (pseudo-)Boolean functions, i.e., functions of the form \(\{0, 1\}^n \to \mathbb{D}\) for some domain \(\mathbb{D}\). Quantum states \(\ket{\psi}\) and gates \(U\) on \(n\) qubits can be viewed as pseudo-Boolean functions \(\psi : \{0,1\}^n \to \mathbb{C}\) and \(U : \{0,1\}^{2n} \to \mathbb{C}\), where \(\psi(\vec{x})\) returns the element of the vector \(\ket{\psi}\) at index \(\vec{x}\), which can be encoded in decision diagrams. A crucial aspect of DDs is that the functions they represent can be manipulated without uncompressing the data structure [9]. It is also well known how to do matrix-vector multiplication (i.e., quantum circuit simulation in our context) in the DD representation [29].

In the context of quantum computing, various types of decision diagrams have been used, including binary decision diagrams (BDDs) [14], [30], multi-terminal BDDs (MTBDDs) [10], [12], quantum multiple-valued DDs (QMDDs) [11], [13], [19], (weighted) context-free-language ordered BDDs (CFLOBDDs) [17], [31], and local invertible map DDs (LIMDDs) [16], [18], [24]. Because we focus on the effect of canonicity of LIMDDs (2.3), we limit our scope to (canonical and non-canonical) LIMDDs as well as QMDDs, of which LIMDDs are a generalization.

Structurally, a DD that encodes a vector \(\ket{\psi}\) is a directed acyclic graph whose internal nodes have two outgoing edges. An example of a QMDD or LIMDD node is shown in 1. Such a node represents the Shannon decomposition, \(\ket{\psi} = \ket{0} \otimes A\cdot \ket{\psi_{0}} + \ket{1} \otimes B\cdot \ket{\psi_{1}}\), where the edge labels \(A, B\) are complex numbers for QMDDs, or local invertible maps (LIMs), which are Kronecker products of \(2\times 2\) invertible matrices, for LIMDDs. For a node \(v\), we use \(\textsf{low}(v)\) and \(\textsf{high}(v)\) to refer to nodes corresponding to subvectors \(\ket{\psi_0}\) and \(\ket{\psi_1}\) respectively. A DD achieves succinctness by recognizing identical or equivalent subvectors and storing only one node for each unique subvector. This recognition of equivalent subvectors is done efficiently by storing DD nodes in a canonical form (details in 2.3). While in general every node has one or two child nodes, there is a single unique terminal node that represents the scalar 1.

QMDDs merge nodes which represent subvectors \(\ket{\phi}, \ket{\psi}\) when these satisfy \(\ket{\phi} = M \cdot \ket{\psi}\), where \(M\) is a complex number. For LIMDDs, this works likewise, but \(M\) is a LIM. This broader merging, based on Kronecker products of local (\(2 \times 2\)) matrices, allows for even stronger compression of the representation. The original work on LIMDDs [16] mainly focuses on Pauli-LIMDDs, where LIMs are restricted to Pauli-LIMs: a complex number times a Pauli word.

Examples of both data types are found in 2. Here, to simplify the notation, edges with value 1 or \(\mathbb{I}\) are left unlabeled. Note that in QMDDs, state vector values can be retrieved by following dashed/solid edges according to the vector index and multiplying all values on the path. Similarly, retrieving values from a LIMDD is done by multiplying the subvectors with the LIMs on the encountered edges. Observe that in this example, the LIMDD can merge the two \(x_2\) nodes because \(\begin{bmatrix} 1 & 2 & 3 & 0 \end{bmatrix}^\intercal= (Z \otimes I)\begin{bmatrix} 1 & 2 & \text{-}3 & 0 \end{bmatrix}^\intercal\).

We note that QMDDs can be seen as LIMDDs with LIMs restricted to \(\lambda I_{2^n}\), with \(\lambda \in \mathbb{C}\). Hence, LIMDDs are strictly more succinct than QMDDs. A more succinct diagram directly benefits quantum-circuit simulation run time as well, as the run time of many DD operations scales with the number of nodes.

Pauli-LIMDDs, similar to the stabilizer formalism [28] and unlike QMDDs, can efficiently simulate Clifford circuits (i.e., circuits composed of Clifford gates). Moreover, it has been shown that there exist classes of circuits where LIMDDs are exponentially faster than both QMDDs and the extended stabilizer formalism for the universal Clifford+T gate set, such as W-state preparation [16]. LIMDDs have also been shown to be exponentially faster on certain classes of circuits than other simulation methods, such as matrix product states (MPS)  [16], [23].

2.3 Canonicity for Pauli-LIMDDs↩︎

While simulating quantum circuits, the DD algorithms create new diagram nodes. As mentioned in 2.2, DDs obtain their succinctness and resulting run time improvements for simulation by recognizing whether an equivalent node is already present in the diagram, so that the new node need not be added. Although one could implement this by checking the equivalence of the newly created node with each existing node in the diagram, the common (and faster) approach is to bring each node into canonical form. This reduces the equivalence check with existing nodes in the diagram to an identity check, which is typically implemented using a hash table for \(\mathcal{O}(1)\) lookup in practice [9].

To define a canonical form for LIMDDs, the authors of [16] realized that two LIMDDs in general could represent Pauli-equivalent states, have all identical nodes merged, and still be different, due to the degrees of freedom listed below and visualized in 3.

1. Because \(0\cdot \vec{\psi} = 0 \cdot \vec{\varphi}\), edges with a 0 label can in principle point to any node. This degree of freedom is made canonical by requiring \(\textsf{low}(v) = \textsf{high}(v)\) whenever one or both has edge label 0.

2. Because \(\ket{0}\otimes\ket{\psi_0} + \ket{1}\otimes\ket{\psi_1} = (X \otimes I^{k}) \left(\ket{1}\otimes\ket{\psi_0} + \ket{0}\otimes\ket{\psi_1}\right)\), there is a degree of freedom regarding the order of \(\textsf{low}(v)\) and \(\textsf{high}(v)\) whenever \(X\) is in the LIM group (which it is in Pauli-LIMs). This order is made canonical by requiring \(\textsf{low}(v) \preccurlyeq\textsf{high}(v)\) according to some total order over all nodes.

3. Because \((I\otimes L) \left(\ket{0}\otimes\ket{\psi_0} + \ket{1}\otimes\ket{\psi_1}\right) = \ket{0}\otimes L\ket{\psi_0} + \ket{1}\otimes L\ket{\psi_1}\), any LIM \(L\) can be factored out of the children’s edge labels. This is brought into a canonical form by requiring that the label of a low edge pointing to a node representing an \(n\)-qubit state is \(I_n\) (whenever it is not zero).

4. (The high determinism rule) Two nodes can have different high labels and still be Pauli-equivalent (i.e., equivalent up to a Pauli word). For example, given a node \(v\) that represents \(\ket{0}\),\(~\raisebox{-.mm}{ \pgftext{\settowidth{\global\pgfcalcparm}{\scriptsize \,\,\,I_2\,\,\,}} \pgftext{\settowidth{\global\pgfcalcparmm}{\scriptsize \,\,\,L\,\,\,}} \tikz{ \node[state,inner sep=0pt,minimum size=10pt] (v){\scriptsize }; \node[state,inner sep=0pt,minimum size=10pt,left=\pgfcalcparm of v](v0){\scriptsize v}; \draw[dotted, ->] (v) to node[above,pos=.45]{\scriptsize I_2} (v0); \node[state,inner sep=0pt,minimum size=10pt,right=\pgfcalcparmm of v](v1){\scriptsize v}; \draw[->] (v) to node[above,pos=.45]{\scriptsize L} (v1) } }\) represents the (unnormalized) state \(\ket{0}\otimes \ket{0} + \ket{1} \otimes \ket{0}\) for both \(L=I_2\) and \(L=Z\). The high label is brought into a canonical form by choosing the minimum high label according to a defined total order on Pauli-LIMs.

A reduced LIMDD is one that does not represent the vector \((0, 0, \dots, 0)^{\intercal}\) and adheres to the above rules as well as the merge rule: it has no duplicate nodes. It is semi-reduced if it satisfies all except high determinism (d). It was proven in [16] that the above rules eliminate all degrees of freedom: two reduced LIMDD have the same root node if and only if they represent equivalent states.

These rules thus establish a canonical form for nodes: given a newly created node \(v\) and edge with label \(L\) pointing to it, we want to find a reduced node \(v'\) and edge label \(L'\) such that \(v'\) satisfies the rules (a)-(d) above, and \(L\ket{v} = L'\ket{v'}\), i.e., they represent the same quantum state. (If we were performing quantum circuit simulation, we would next check if \(v'\) is already present in the LIMDD.)

Given any node \((v,L)\), finding a semi-reduced \((v', L')\), i.e., that satisfies the first three rules (a)-(c), is straightforward and can be done in time \(\mathcal{O}(n)\); see 3 for how to update \((v, L)\). For Pauli-LIMDDs, for the final rule (d), high determinism, an \(\mathcal{O}(n^3)\) algorithm exists [16], which we improve upon in 3.

= [state, ellipse, minimum size=4mm, inner sep=0.5mm, node distance=1cm]

=[draw, rectangle,minimum size=4.mm, inner sep=3pt] =[circle,draw=black!70,solid,thick,minimum size=6mm] =[regular polygon, regular polygon sides=3, draw=black!70,solid,thick,inner sep=0.5mm] =[->,loosely dashed,thick] =[->,solid,thick] =[->,densely dashdotted,ultra thick] =[draw,fill=white,inner sep=2pt, minimum size=0cm,line width=.5pt]

2.4 Related work: implementations↩︎

One implementation of these Pauli-LIMDDs is MQT-LIMDD  [24], which is an extension of the QMDD implementation MQT-DDSIM [13]. However, MQT-LIMDD implements canonicity partially, and does not achive minimal size LIMDDs in all cases. Aside from this implementation, LimTDDs [18] have extended LIMDDs to a more expressive group of LIMs, specifically all matrices generated under Kronecker and matrix products from \(\{X, P(\tfrac{2\pi}{N})\}\) (and a complex scalar), where \(N \in \mathbb{N}\) can be chosen freely. This XP group is a generalization of the one generated by the Pauli matrices, as setting \(N=2\) gives \(P(\pi) = Z\), \(P(\pi)X = iY\), and \(XX=I\) which yields Pauli-LIMs. However, the implementation accompanying [18] does, by design, not guarantee canonicity and thus possibly yields LIMDDs that are not minimal in size.

3 Faster algorithm for High Determinism↩︎

Here, we derive a new algorithm for transforming a semi-reduced Pauli-LIMDD node into one that satisfies the high determinism rule (2.3) with the following advantages over the original algorithm of [16]: (a) For the case that the node has a single child (\(n\) qubits, \(m\) stabilizer generators), a run time of \(\mathcal{O}(nm)\) instead of \(\mathcal{O}(nm^2)\); (b) For the case of two distinct children, the same worst-case \(\mathcal{O}(n^3)\) run time but with less overhead. Our advantage is gained by adjusting a well-known algorithm for computing intersections and unions of vector spaces.

3.1 High determinism algorithm: problem statement↩︎

Given \(n\)-qubit semi-reduced node \(u := ~\raisebox{-.mm}{ \pgftext{\settowidth{\global\pgfcalcparm}{\scriptsize \,\,\,I\,\,\,}} \pgftext{\settowidth{\global\pgfcalcparmm}{\scriptsize \,\,\,\beta B\,\,\,}} \tikz{ \node[state,inner sep=0pt,minimum size=10pt] (v){\scriptsize }; \node[state,inner sep=0pt,minimum size=10pt,left=\pgfcalcparm of v](v0){\scriptsize v_0}; \draw[dotted, ->] (v) to node[above,pos=.45]{\scriptsize I} (v0); \node[state,inner sep=0pt,minimum size=10pt,right=\pgfcalcparmm of v](v1){\scriptsize v_1}; \draw[->] (v) to node[above,pos=.45]{\scriptsize \beta B} (v1) } }\) with \(\beta \in \mathbb{C}, B \in \mathrm{\small Pauli}_n\), it was shown in [16] that the node \(u_\text{can} := ~\raisebox{-.mm}{ \pgftext{\settowidth{\global\pgfcalcparm}{\scriptsize \,\,\,I\,\,\,}} \pgftext{\settowidth{\global\pgfcalcparmm}{\scriptsize \,\,\,\alpha A\,\,\,}} \tikz{ \node[state,inner sep=0pt,minimum size=10pt] (v){\scriptsize }; \node[state,inner sep=0pt,minimum size=10pt,left=\pgfcalcparm of v](v0){\scriptsize v_0}; \draw[dotted, ->] (v) to node[above,pos=.45]{\scriptsize I} (v0); \node[state,inner sep=0pt,minimum size=10pt,right=\pgfcalcparmm of v](v1){\scriptsize v_1}; \draw[->] (v) to node[above,pos=.45]{\scriptsize \alpha A} (v1) } }\) (\(\alpha \in \mathbb{C}, A \in \mathrm{\small Pauli}_n\)) is a valid choice for the reduced version of \(u\), where \[\label{eq:ucan} \alpha A = \min_{\substack{s, x \in \{0, 1\},\\ g\in \text{Stab}(v_0), h \in \text{Stab}(v_1)}} \{(-1)^s \cdot \beta^{(-1)^x} g B h \mid x \neq 1 \textrm{ if } v_0 \neq v_1\} .\tag{1}\] Indeed, \(u\) and \(u_\text{can}\) are Pauli-equivalent as [16] \[\label{eq:xzgmin} \ket{u_\text{can}} = (X \otimes \beta B)^x \cdot (Z^s \otimes g_\text{min}) \ket{u}\tag{2}\] where \(g_\text{min} \in \text{Stab}(v_0)\) achieves the minimum in 1 . Canonicity is enforced by the minimization, defined through the following total order [16]. In the definition below, for a complex number \(r\cdot e^{i\theta}\) in polar form, i.e. \(r \in \mathbb{R}_{\geq 0}\) and \(\theta \in [0, 2\pi)\) and \(P\) a Pauli word, we denote the vector \(\vec{re^{i\theta} P}\) as the concatenation of \(\textsf{revbin}\mathopen{}\mathclose{\left(P}\right)\) and the tuple \((r, \theta)\), for example \(\vec{3\cdot e^{\pi/2} \cdot YZ} = (0, 1, 1, 1, 3, \pi/2)\).

For example, \(\vec{-\tfrac{1}{2}X} = \vec{\tfrac{1}{2} e^{i\pi} X} = (1, 0, \tfrac{1}{2}, \pi) > (0, 1, 2, 0) = \vec{2 e^{0} Z} = \vec{2 Z}\).

The goal of the high determinism algorithm is to evaluate eq. 1 , i.e. to find \(\alpha A\) given \(\beta B\) and \(\text{Stab}(v_0), \text{Stab}(v_1)\). In general, there are multiple tuples \((g, h, s, x)\) that achieve the minimum of 1 . However we show in 6 that any two such tuples \((g, h, s, x), (g', h', s', x')\) satisfy \(g\cdot B\cdot h = \pm g'\cdot B \cdot h'\). Because the Pauli-LIM ordering is lexicographic, i.e., the scalars only matter if the Pauli words are equal, if we found \(g,h\) that is off by a factor of \(-1\) of the true minimum, we correct for this by choosing \(s\) accordingly to get the right factor \((-1)^s\). This is why we may greedily perform the minimization in eq. 1 : We first find the Pauli word \(A\) as \(\min_{g \in \text{Stab}(v_0), h \in \text{Stab}(v_0)} \textsf{word}\mathopen{}\mathclose{\left(gBh}\right)\), where \(\textsf{word}\mathopen{}\mathclose{\left(L}\right)\) is the Pauli word of Pauli LIM \(L\), e.g. \(\textsf{word}\mathopen{}\mathclose{\left(+XY}\right) = XY\) and \(\textsf{word}\mathopen{}\mathclose{\left(-iZX}\right) = ZX\). We give our algorithm for doing so in 3.2. Next, we find the minimal scalar \(\alpha\) by trying the four combinations of \(x, s \in \{0, 1\}\).

3.2 Finding the minimal Pauli word \(A\)↩︎

We here assume that we are given \(n\)-qubit stabilizer generating sets \(S_0, S_1\) and an \(n\)-qubit Pauli word \(B\). We will describe an algorithm to compute \[\label{eq:minimal-objective} (g_\text{min}, h_\text{min}) = \operatornamewithlimits{argmin}_{g\in \langle S_0\rangle, h \in \langle S_1 \rangle} \textsf{word}\mathopen{}\mathclose{\left(gBh}\right) .\tag{3}\] In the context of the high determinism rule, \(\langle S_0\rangle\) (\(\langle S_1 \rangle\)) is the stabilizer group of child node \(v_0\) (\(v_1\)), but our algorithm is agnostic to the interpretation of \(S_0\) and \(S_1\), except for whether \(v_0\) and \(v_1\) are identical reduced nodes.

The algorithm from [16] to solve eq. 3 consists of many individual steps, which include the calculation of intersections of stabilizer groups and cosets of stabilizer groups. We here sidestep having to calculate these by adjusting one part of the algorithm, the Zassenhaus procedure, to directly solve eq. 3 . Below, in 3.3, we will also re-use the adjusted algorithm’s output to derive a faster and simpler Stabilizer Finding algorithm.

Given two sets of basis vectors \(V, W \subset \mathbb{R}^m\) for some positive integer \(m\), the Zassenhaus algorithm [32] finds bases for the vector spaces \(V \cap W\) and \(\langle V \cup W \rangle\) (where \(\langle.\rangle\) denotes closure under vector addition). To do so, it defines a matrix of \(2m\) columns, which we call the Zassenhaus matrix, with rows \(\{v \circ v \mid v \in V\} \cup \{w\circ 000\dots0 \mid w \in W\}\), where \(\circ\) denotes vector concatenation and ‘\(000\dots0\)’ is the zero vector. The algorithm then brings the left register of the matrix (the leftmost \(m\) columns) into row-echelon form. To read off a basis for the union \(\langle V \cup W \rangle\) after this, one takes all vectors in the left register that are not equal to the zero vector, while a basis for the intersection \(V \cap W\) is the set of all (nonzero) vectors in the right register where the left vectors are zero.

We here expand the algorithm to signed Pauli words (recall that stabilizers can only have scalar \(\pm 1\)) to create the Zassenhaus matrix \(M\) consisting of three registers: \(M = \{ \textsf{word}\mathopen{}\mathclose{\left(g}\right) \circ g \circ (+I_{2^n}) \mid g \in S_0\} \cup \{\textsf{word}\mathopen{}\mathclose{\left(h}\right)\circ (+I_{2^n}) \circ h \mid h \in S_1\}\). Thus, the left register contains length-\(2n\) binary vectors, the middle only elements of \(\langle S_0\rangle\), and the right register only elements of \(\langle S_1\rangle\). We next bring the left register in REF by performing Gaussian elimination on \(M\), where addition (\(+\)) of rows is defined using Pauli-LIM multiplication (\(\cdot\)) as \((\textsf{word}\mathopen{}\mathclose{\left(k}\right), g, h) + (\textsf{word}\mathopen{}\mathclose{\left(k'}\right), g', h') = (\textsf{word}\mathopen{}\mathclose{\left(k \cdot k'}\right), g' \cdot g, h \cdot h')\). Because \(\langle S_0\rangle\) and \(\langle S_1 \rangle\) are stabilizer groups, whose elements commute, this addition is commutative, as Gaussian elimination requires. As a consequence, the resulting matrix \(M'\) has the following properties:

1. A basis of minimal size for the binary vector space \(\langle \textsf{word}\mathopen{}\mathclose{\left(g \cdot h}\right) \mid g \in \langle S_0 \rangle, h \in \langle S_1 \rangle \rangle\) is found as the set of vectors in the left register that are nonzero. This follows directly from the fact that the left register contains (the Pauli words of) all elements of both \(S_0\) and \(S_1\)

2. Each left vector is the product of the middle and right Pauli stabilizers: each row \((k, g, h)\) satisfies \(k = \textsf{word}\mathopen{}\mathclose{\left(g \cdot h}\right)\). This is an invariant of the Gaussian elimination algorithm and is hence straightforwardly proven by induction on the algorithm’s steps.

To find \(g_\text{min}\) and \(h_\text{min}\) from 2 , we now perform the standard algorithm for finding the minimal vector between a vector space \(V\) and a hyperplane \(\{v_0 + v \mid v \in V\}\) for a given vector \(v_0\) [16]. First, append a row \(\textsf{word}\mathopen{}\mathclose{\left(B}\right) \circ (+I_{2^n}) \circ (+I_{2^n})\) at the bottom of \(M\). Next, add rows downward to minimize \(\textsf{word}\mathopen{}\mathclose{\left(B}\right)\) in the left register: iterating a column index \(c\) from \(1\) to \(2n\), if the value of the bottom row \(b\) at \(c\) contains the value \(1\) and there is a row \(r\) which contains a leading entry \(1\) at column \(c\), then replace \(b\) by \(b+r\). Finally, the bottom row’s left register will contain \(\textsf{word}\mathopen{}\mathclose{\left(g_{\textrm{min}} \cdot B \cdot h_{\textrm{min}}}\right)\) by property [item:first-property] above, while by property [item:second-property], \(g_\text{min}\) and \(h_\text{min}\) are found in the appended row’s middle and right registers, respectively.

The run time of the algorithm above is dominated by the Zassenhaus procedure, taking time \(\mathcal{O}(n \cdot (|S_0| + |S_1|)^2)\), i.e., at most \(\mathcal{O}(n^3)\), identically to the original algorithm in [16]. However, in the case that \(u\) only has a single child node, i.e. \(v_0 = v_1\) and hence \(\langle S_0 \rangle = \langle S_1\rangle\), the algorithm for finding \(A\) is simpler: note that \(\min_{g, h \in \langle S_0\rangle} \textsf{word}\mathopen{}\mathclose{\left(gBh}\right) = \min_{g, h \in \langle S_0\rangle} \textsf{word}\mathopen{}\mathclose{\left(\pm B gh}\right) =\min_{h \in \langle S_0\rangle} \textsf{word}\mathopen{}\mathclose{\left(Bh}\right)\) hence we only need to minimize \(B\) by multiplying with elements of \(\langle S_0\rangle\), i.e. performing the minimal-vector algorithm only, an \(\mathcal{O}(n\cdot |S_0|)\) procedure.

As \(|S_0| \leq n\), this improves the original algorithm’s \(\mathcal{O}(n^3)\) run time to \(\mathcal{O}(n^2)\).

3.3 Finding the stabilizer group of a reduced node↩︎

So far, we have shown how to reduce a node by minimizing 1 , provided that we have generator sets for the stabilizer groups of the node’s children, in row-echelon form (REF). Here, we give a recursive algorithm to find the latter: find an REF generating set for the stabilizer group of a reduced node. Compared to the original algorithm [16], our algorithm has less overhead as it avoids computing group intersections multiple times by re-using the Zassenhaus result from the high determinism algorithm in 3.2, and has a speedup from \(\mathcal{O}(n^3)\) to \(\mathcal{O}(n^2)\) for reduced-nodes with a single child.

The base case for our recursive algorithm is the leaf node, which has a stabilizer group consisting only of the number \(1\). For the inductive case, suppose that we have an \(n\)-qubit reduced node \(~\raisebox{-.mm}{ \pgftext{\settowidth{\global\pgfcalcparm}{\scriptsize \,\,\,I\,\,\,}} \pgftext{\settowidth{\global\pgfcalcparmm}{\scriptsize \,\,\,\alpha A\,\,\,}} \tikz{ \node[state,inner sep=0pt,minimum size=10pt] (v){\scriptsize }; \node[state,inner sep=0pt,minimum size=10pt,left=\pgfcalcparm of v](v0){\scriptsize \phi_0}; \draw[dotted, ->] (v) to node[above,pos=.45]{\scriptsize I} (v0); \node[state,inner sep=0pt,minimum size=10pt,right=\pgfcalcparmm of v](v1){\scriptsize \phi_1}; \draw[->] (v) to node[above,pos=.45]{\scriptsize \alpha A} (v1) } }\) where \(\alpha \in \mathbb{C}, A \in \mathrm{\small Pauli}_{n-1}\) and \(\phi_0, \phi_1\) are reduced nodes. Denote \(S_0\) (\(S_1\)) for an REF (row-echelon form) generating set of node \(\phi_0\) (\(\phi_1\)). Our goal is to find all \(\lambda P\) such that \(\lambda P \ket{\phi} = \ket{\phi}\), where \(\lambda \in \{\pm 1\}\) and \(P \in \mathrm{\small Pauli}_n\).

We first expand \(\lambda P = P_1 \otimes P'\) where \(P_1 \in \{I, X, Y, Z\}\) (absorbing \(\lambda\) into the Pauli-LIM \(P'\)) and distinguish cases \(P_1 \in \{I, Z\}\) and \(P_1\in \{X, Y\}\). For the case \(P_1 \in \{I, Z\}\) we write \(P_1 = { \left(\begin{smallmatrix}1&0\\0&x\end{smallmatrix}\right)}\) for \(x \in \{1, -1\}\). Because \(\ket{\phi} = \ket{0} \otimes \ket{\phi_0} + \ket{1}\otimes \alpha A \ket{\phi_1}\), in this case \(P_1 \otimes P' \ket{\phi} = \ket{\phi}\) is equivalent to \[\label{eq:iz-derivation} P' \ket{\phi_0} = \ket{\phi_0} \quad \text{and} \quad P' \cdot \alpha A \ket{\phi_1} = x \alpha A \ket{\phi_1}\tag{4}\] If \(\alpha = 0\), then the solutions \(P'\) are precisely all stabilisers of \(\ket{\phi_0}\). That is, the stabilizers of \(\phi\) are \(\{P_1 \otimes s \mid P_1 \in \{I, Z\}, s \in S_0\}\); we return \(\{Z \otimes I_{n-1}\} \cup \{I \otimes s \mid s \in S_0\}\) as generating set. If \(\alpha \neq 0\), we rewrite the latter equation of 4 as \(x \sigma(A, P') P' \ket{\phi_1} = \ket{\phi_1}\) using the fact that \(x = 1/x\) and where \(\sigma(\cdot, \cdot) = 1\) if the arguments commute and \(-1\) otherwise. So, a generating set for \(\phi\)’s stabilizers is \[\label{eq:zcase} \left\{\begin{pmatrix}1&0\\0&\sigma(g, A) \cdot \gamma_g\end{pmatrix} \otimes g \mid g \in S'\right\}\tag{5}\] where \(S'\) is an REF generating set for the set of stabilisers \(g\) of \(\ket{\phi_0}\) such that \(g \ket{\phi_1} = \gamma_g \cdot \ket{\phi_1}\) for some \(\gamma_g \in \{\pm 1\}\).

If \(v_0 = v_1\), then \(S' = S_0\), otherwise \(S\) is found as \(S' = \{g \mid (000\dots0, g, \gamma \cdot g) \in M'\}\), where \(M'\) is the Zassenhaus matrix of \(S_0\) and \(S_1\) (see 3.2).

Second, for the case \(P_1 \in \{X, Y\}\) we write \(P_1 = { \left(\begin{smallmatrix}0&y^*\\y&0\end{smallmatrix}\right)}\) with \(y \in \{1, i\}\). Then \(P_1 \otimes P' \ket{\phi} = \ket{\phi}\) is expanded as \[\label{eq:y-derivation} y^*\alpha P' A \ket{\phi_1} = \ket{\phi_0} \quad \text{and} \quad y P'\ket{\phi_0} = \alpha A \ket{\phi_1}\tag{6}\] If \(\alpha = 0\), then there are no solutions because that would imply that \(\ket{\phi_0}=\ket{\phi_1}=0\), which is not a reduced node by definition [16].

If \(\alpha \neq 0\), then \(\ket{\phi_0}\) and \(\ket{\phi_1}\) are Pauli-isomorphic. But by assumption, \(\phi_0\) and \(\phi_1\) are reduced, hence \(\phi_0 = \phi_1\).

Using the fact and multiplying both sides of the latter equation with \(y^* P'\), 6 yields the same equation twice \[\label{eq:y-derivation-two} y^* \alpha P' A \ket{\phi_0} = \ket{\phi_0}\tag{7}\] where we used that \(y^* = 1/y\) and that \(P'\) and \(A\) are their own inverse. Since all Pauli-LIMs which are stabilizers have scalar either \(1\) or \(-1\) [27], while \(y\in \{1, i\}\) and the multiplication of \(P'\) and \(A\) yields a scalar \(\in \{\pm 1\}\), it follows that a solution in this case is only possible if \(\alpha = \pm y\), which in particular implies that \(\alpha \in \{\pm 1, \pm i\}\). Provided that that is the case, a solution to 7 is \(P' = y/\alpha A\). That is, the second case yields \[\label{eq:xcase} \begin{pmatrix}0&\alpha^* \\ \alpha &0\end{pmatrix} \otimes A\tag{8}\] if \(\alpha \in \{\pm 1, \pm i\}\), and no stabilizer generators otherwise. Since any generating set of a stabilizer group can be brought back to a form where only a single generator has \(P_1 \in \{X, Y\}\) as the first Pauli entry [33], it suffices to return this one only.

The above algorithm returns a generating set for a node’s stabilizer group. To return an REF generating set, we call Gaussian elimination on the output in case \(\alpha \neq 0\). In the case of distinct children with \(\alpha \neq 0\), the algorithm’s run time is dominated by the Zassenhaus procedure, taking \(\mathcal{O}(n \cdot (|S_0| + |S_1|)^2)\) time. In the case of a single child with \(\alpha \neq 0\), the procedure takes \(\mathcal{O}(n\cdot |S_0|)\) time to calculate \(\sigma(s, A)\) for all \(s \in S_0\). As \(S_0\) is already in REF by assumption and our REF definition has the most significant qubit at the right (see \(\textsf{revbin}\mathopen{}\mathclose{\left(}\right)\) definition in 2), the set in 5 is in REF already. Hence, Gaussian elimination only updates the single row appended in 8 , taking \(\mathcal{O}(n \cdot |S_0|)\) time. If \(\alpha = 0\), then the output is already in REF, as the last row is \(\textsf{revbin}\mathopen{}\mathclose{\left(Z\otimes I_{n-1}}\right) = (0000\dots001)\), i.e., it is majorized by \(\textsf{revbin}\mathopen{}\mathclose{\left(I\otimes s}\right)\) for any stabilizer \(s \neq I_{2^n}\). As \(|S_0|, |S_1| \leq n\), the run times are at most \(\mathcal{O}(n^3), O(n^2)\) and \(\mathcal{O}(n^2)\).

4 Empirical evaluation↩︎

We implemented our high determinism algorithm (3) as an open-source C library called PauliLim [25], and we have created a from-scratch LIMDD implementation in C/C++ named QolDDer [26]. Together, these allow us to investigate the effects of canonicity on the simulation of quantum circuits. The implementations are separated to improve modularity and extensibility for future work. Below, we highlight two important design decisions of QolDDer.

As is typical for a DD implementation, to efficiently recognize equivalent nodes, QolDDer’s nodes are stored in a hash table after bringing them into a canonical form [9]. We resolve hash collisions through chaining. Each node is assigned a unique ID based on its index in the table. These IDs are then used as canonical representations of the children when hashing a parent node.

While most of the information that uniquely defines a node (qubit number, children nodes, Pauli words of edge labels) is discrete and can be hashed directly; scalars need to be treated separately. Floating-point computations are subject to small rounding errors that can prevent nodes from merging (e.g., 1/sqrt(2) yields 0.707...5, while cos(pi/4) yields 0.707...6). Most DD implementations that include floating points handle this by setting a threshold \(\delta\) and consider two scalars \(\alpha\) and \(\alpha'\) equivalent when \(|\alpha - \alpha'| \leq \delta\) [34]–[36]. In our case, with \(\alpha,\alpha' \in \mathbb{C}\), we consider them equivalent when \(|\alpha_{\text{real}} - \alpha'_{\text{real}}| \leq \delta\) and \(|\alpha_{\text{imag}} - \alpha'_{\text{imag}}| \leq \delta\). We leave \(\delta\) configurable, and set the default \(\delta = 10^{-13}\); for 64-bit floats, \(\delta\) should generally be greater than \(10^{-15}\) [37]. To find \(\delta\)-close scalars, we store them in a search tree. Additionally, in our implementation, we make a minor change to the LIMDD normalization rules that has been shown to mitigate numerical errors in QMDDs [19] while preserving LIMDD canonicity. Details are in App. 7.

The correctness of a LIMDD implementation has two components: representing the correct state vector at all times and representing it using reduced nodes only. We rigorously tested QolDDer on the first aspect by verifying that QolDDer’s output, converted to a state vector, is identical to the output of MQT-DDSIM and LimTDD on the entire MQT-Bench benchmark set [38] (see 4.1) up to 10 qubits, and for the first 10 vector entries up to 20 qubits. We tested the second aspect by verifying PauliLim’s correctness against a simple brute-force implementation that solves 1 on \(10^5\) randomly-generated 4 to 8-qubit instances of high edge labels and stabilizer groups.

4.1 Empirical results↩︎

We investigated the effect of canonicity on simulation run time using two benchmark sets. First, we generate a set of random Clifford gates, RandCliff, as Pauli-LIMDDs are known to outperform QMDDs on such circuits [16], [23]. It consists of randomly generated \(n\)-qubit Clifford circuits with \(100n\) gates, ranging from \(n=10\) to \(30\) qubits, each with \(10\) instances. Gates are drawn uniformly from the Clifford generating set \(\{H,P(\tfrac{\pi}{2}),\mathit{CX}\}\). Next, MQT-Bench [38] is an existing benchmark set comprised of a broad range of circuit types.

All benchmarks were run on an AMD Ryzen 7 9800X3D CPU, with 64 GB of available memory. Each run uses a 10-minute timeout and a 16 GB memory limit. Reproducible benchmarks are available online [39].

None

By turning the use of high determinism on or off, we investigate the effect of simulation on RandCliff with reduced versus semi-reduced Pauli-LIMDDs (see 2.3). 4 shows a run time comparison. We observe that full canonicity results in significantly faster run times. On this log-log scale, the data points follow a trend that appears linear, or even sub-linear, indicating a polynomial run time improvement or better. To verify that the effect of canonical simulation is due to a reduction in node count, we show the number of nodes in the DD of the output state of the random Clifford circuits in 4b. Here we see that, indeed, the node count is significantly larger for semi-reduced Pauli-LIMDDs.

None

Next, in 5, we compare QolDDer’s run time against the QMDD simulator MQT-DDSIM and two LIMDD simulators MQT-LIMDD and LimTDD (see also 2.4). Based on the observed sublinear scaling in the log-log plot of 5, we see that QolDDer is exponentially faster than QMDDs on RandCliff. This is as expected, as LIMDDs can represent the output of Clifford circuits in polynomial space, but QMDDs cannot succinctly represent highly entangled states [23], which random Clifford circuits produce [40]. We also observe exponential speedup compared to MQT-LIMDD. We expect that this is due to the fact that their canonicity is only partially implemented; in theory, our algorithm should only show polynomial improvement (see 3). In 5c, we observe that QolDDer is a constant 25 times faster than LimTDD on some instances, while on larger instances LimTDD times out while QolDDer does not. This difference is likely due to QolDDer’s canonicity, as we observed (not depicted) that LimTDD has much larger node numbers on this benchmark set.

In 5d-f, we observe that the advantages of canonicity do not hold on all circuits in MQT-Bench. In 5d and f, we grouped the benchmarks into sets A = {Quantum Phase Estimation (QPE) exact, QPE inexact}, B = {Quantum Neural Networks (QNN), Quantum random walks, Grover}, and C = {all other} to highlight the observed trends. In 5e, only the Greenberger–Horne–Zeilinger state preparation (GHZ) and Deutsch-Josza (DJ) circuits have been highlighted.

In 5d, we observe that QolDDer significantly outperforms MQT-DDSIM on QPE circuits. However, we also observe that the latter is faster on QNN, Quantum random walks, and Grover algorithms (“Set B”). We speculate this is because the MQT-Bench benchmark set contains many circuits for which LIMs do not achieve more merging than QMDDs. In that case, the LIMs in the DD only add overhead. Nonetheless, we observe roughly a linear scaling on Sets B and C, indicating that QolDDer is only slower by a constant factor of \(20\).

In 5e, we see that QolDDer’s performance is similar to MQT-LIMDD on MQT-Bench, with two types of outliers. First, on GHZ and DJ circuits, QolDDer shows improved scaling over MQT-LIMDD, although due to their sub-1-second run times, we draw no conclusions from these. Second, on various circuits, MQT-LIMDD times out or terminates with an error, indicating possible issues in the implementation. In 5f, we observe a constant-factor speedup of QolDDer over LimTDD on QPE (Set A; \(30\times\) faster). For others, this factor decreases as circuit size increases. We speculate this is due to the larger LIM group that LimTDD has (XP, a supergroup of the Pauli group), enabling more merging, despite not being canonical.

None

To investigate the hypothesis that LIMDD canonicity does not speed up run times for MQT-Bench, we compared semi-reduced and reduced Pauli-LIMDDs in 6. We see little performance difference (6a), which is explained by a limited gain in succinctness (6b).

5 Discussion↩︎

We presented a novel, more efficient algorithm to bring LIMDD nodes into a canonical form, and implemented this algorithm as part of our LIMDD implementation QolDDer (both open source [25], [26]). In our empirical evaluation, we have compared both run times and node counts for quantum circuit simulation with canonical and non-canonical LIMDDs. Our results corroborate that LIMDDs can achieve exponential advantages over QMDDs. Also, they show that, without canonicity, LIMDDs cannot always realize these advantages. On some benchmarks, however, ensuring canonicity is unnecessary overhead.

In the future, we will extend our canonicity algorithm beyond the Pauli group. Specifically, we expect that our algorithm is extendable to the XP group, which LimTDD uses [18]. We will also apply LIMDDs to a wider range of circuits, such as the universal Clifford+\(T\) circuits, for which LIMDDs have scaling guarantees [41]. Moreover, we will look into applications and tasks beyond quantum circuit simulation, such as quantum circuit equivalence checking, a task on which QMDDs have shown success  [19], [42]. In the classical domain, we will investigate the use of LIMDDs for linear-algebra tasks such as stochastic reasoning and Markov-chain analysis, investigating whether relevant tasks possess structure that the strong merge rule of LIMDDs could explicitly make use of.

5.0.1 ↩︎

SB is supported by the NWO Quantum Technology program (project number NGF.1582.22.035). AQ acknowledges the support received though the NWO NWA‑ORC program, project Divide & Quantum (NWA.1389.20.241). TC is supported by an NWO Veni grant (VI.Veni.232.381).

6 Proof of High Determinism minimizer form↩︎

Fix the number of qubits \(n \geq 0\). Let \(G_0, G_1\) be \(n\)-qubit stabilizer groups and let \(B\) be an \(n\)-qubit Pauli word. Suppose that both \((g, h, s, x), (g', h', s', x') \in G_0 \times G_1 \times \{0, 1\} \times \{0, 1\}\) achieve the minimum 1 . We here show that \(g\cdot B \cdot h = \pm g' \cdot B \cdot h'\).

First, we note that by definition of the total order on PauliLIMs in 1, \(\textsf{word}\mathopen{}\mathclose{\left(gBh}\right) = \textsf{word}\mathopen{}\mathclose{\left(g'Bh'}\right)\). Hence \(gBh = \alpha \cdot g'Bh'\) for some \(\alpha \in \mathbb{C}\). Now we use the fact that stabilizers are their own inverse (because Pauli words square to identity, and the scalars of stabilizers are \(\pm 1\)), we reorder factors to find that \[(g' g) B (h h') = \alpha B\] hence \[(g' g) B (h h') B = \alpha I_{2^n} .\] As Pauli words either commute or anticommute, we find \[\pm (g' g) B B (h h') = \alpha I_{2^n}\] so \[\label{eq:pmgprimeg} \pm (g'g) (h h') = \alpha I_{2^n} \text{.}\tag{9}\] As \(g'\cdot g\) and \(h\cdot h'\) are stabilizers, their scalars are \(\pm 1\). Moreover, following 9 we see that \(g'g\) and \(hh'\) multiply to a multiple of \(I_{2^n}\), which is only possible for PauliLIMs if \(\textsf{word}\mathopen{}\mathclose{\left(g' \cdot g}\right) = \textsf{word}\mathopen{}\mathclose{\left(h \cdot h'}\right)\). Hence, we may write \(g'\cdot g = \pm W\) and \(h \cdot h' = \pm W\) for some Pauli word \(W\). Plugging back into 9 , we find \[\pm (\pm W) \cdot (\pm W) = \alpha I_{2^n} .\] Hence \(\alpha = \pm 1\), which finishes our proof.

7 Scalar normalization↩︎

Given a node with low edge label \(A\) and a high edge label \(B\), the low-factoring rule (see 3) normalizes the edge labels by multiplying the root edge label with \(A\), setting the high edge label to \(A^{-1} B\), and leaving the low edge an identity label. This means that, for \(\alpha\) and \(\beta\) being the scalar components of \(A\) and \(B\) respectively, the normalized scalars will be \(\alpha' = \alpha/\alpha = 1\) and \(\beta' = \beta/\alpha\). However, for QMDDs, it has been shown that the issues that arise from floating-point imprecision can be mitigated by using a different normalization strategy [19]. Specifically, instead of dividing \(\alpha\) and \(\beta\) always by \(\alpha\), dividing by \(\max(\alpha, \beta)\) yields fewer numerical errors.

In LIMDDs, this can be implemented by either modifying the low factoring rule (3c) or by modifying the low precedence rule (3b). We opt for the latter as it is simpler in implementation. The modified low precedence rule, which we shall call the max precedence rule, is given in 2. In 3c, the low factoring rule has “\(u \preccurlyeq w\)” in the caption. This condition would be replaced by “max precedence holds”.

We now give an alternative definition of a LIMDD to be reduced, where the Low precedence rule (Fig. 3b) is replaced by Max precedence.

What remains is to show that max-reduced LIMDD is canonical, which is the equivalent of the proof that reduced LIMDD is canonical in [16].

References↩︎