LIMDD: A Decision Diagram for Simulation of
Quantum Computing Including Stabilizer States
Lieuwe Vinkhuijzen∗,1 , Tim Coopmans∗,1,2 , David Elkouss2,3 , Vedran Dunjko1 , andAlfons Laarman1
1 Leiden University, The Netherlands
2 Delft University of Technology, The Netherlands

arXiv:2108.00931v5 [quant-ph] 6 Sep 2023

3 Networked Quantum Devices Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan

Efficient methods for the representation and simulation of quantum states and quantum operations are crucial for the optimization of quantum circuits. Decision diagrams(DDs), a well-studied data structure originally used to represent Boolean functions,have proven capable of capturing relevant aspects of quantum systems, but their limits are not well understood. In this work, we investigate and bridge the gap between
existing DD-based structures and the stabilizer formalism, an important tool for simulating quantum circuits in the tractable regime. We first show that although DDs
were suggested to succinctly represent important quantum states, they actually require exponential space for certain stabilizer states. To remedy this, we introduce amore powerful decision diagram variant, called Local Invertible Map-DD (LIMDD). We
prove that the set of quantum states represented by poly-sized LIMDDs strictly contains the union of stabilizer states and other decision diagram variants. Finally, there
exist circuits which LIMDDs can efficiently simulate, while their output states cannot
be succinctly represented by two state-of-the-art simulation paradigms: the stabilizer
decomposition techniques for Clifford + T circuits and Matrix-Product States. By
uniting two successful approaches, LIMDDs thus pave the way for fundamentally more
powerful solutions for simulation and analysis of quantum computing.

∗ These authors contributed equally.
Lieuwe Vinkhuijzen∗, : l.t.vinkhuijzen@liacs.leidenuniv.nl

Tim Coopmans∗, : t.j.coopmans@liacs.leidenuniv.nl, https://timcp.github.io/
David Elkouss: d.elkousscoronas@tudelft.nl, https://www.davidelkouss.com
Vedran Dunjko: v.dunjko@liacs.leidenuniv.nl, https://liacs.leidenuniv.nl/∼dunjkov
Alfons Laarman: a.w.laarman@liacs.leidenuniv.nl, https://alfons.laarman.com

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

1

Contents
1

Introduction

3

2Preliminaries

4

2.1

Decision diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Pauli operators and stabilizer states . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.3

Matrix product states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

3

4

Local Invertible Map Decision Diagrams

7

3.1

The LIMDD data structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

3.2

Succinctness of LIMDDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

3.3

Pauli-LIMDD manipulation algorithms for simulation of quantum computing . . . .

14

3.4

Comparing LIMDD-based simulation with other methods . . . . . . . . . . . . . . .

21

Canonicity: Reduced LIMDDs with efficient MakeEdge algorithm

23

4.1

LIMDD canonical form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4.2

The MakeEdge subroutine: Maintaining canonicity during simulation . . . . . . .

28

5Related work

29

6

Discussion

30

7

Acknowledgements

31

A Linear-algebra algorithms for Pauli operators

36

B Proof that cluster states and coset states need exponentially-large QMDDs

37

C How to write graph states, coset states and stabilizer states as Tower-LIMDDs

40

D Efficient algorithms for choosing a canonical high label

44

D.1 Constructing the stabilizer subgroup of a LIMDD node . . . . . . . . . . . . . . . .

47

D.2 Efficiently finding a minimal LIM by multiplying with stabilizers . . . . . . . . . .

50

E Measuring an arbitrary qubit

52

F LIMDDs prepare the W state efficiently

54

G Numerical search for the stabilizer rank of Dicke states

59

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

2

1 Introduction
Classical simulation of quantum computing is useful for circuit design [1, 2], verification [3, 4] andstudying noise resilience in the era of Noisy Intermediate-Scale Quantum (NISQ) computers [5].
Moreover, identifying classes of quantum circuits that are classically simulatable, helps in excluding
regions where a quantum computational advantage cannot be obtained. For example, circuits
containing only Clifford gates (a non-universal quantum gate set), using an all-zero initial state,
only compute the so-called ‘stabilizer states’ and can be simulated in polynomial time [6–10].
Stabilizer states, and associated formalisms for expressing them, are fundamental to many quantum
error correcting codes [8] and play a role in measurement-based quantum computation [11]. In fact,
simulation of universal quantum circuits is fixed-parameter tractable in the number of non-Cliffordgates [12], which is why many modern simulators are based on stabilizer decomposition [12–17].
Another method for simulating universal quantum computation
is based on decision diagrams (DDs) [18–21], including Algebraic
DDs [22–25], Affine Algebraic DDs [26], Quantum Multi-valued
DDs [27, 28], and Tensor DDs [29]. A DD is a directed acyclic
graph (DAG) in which each path represents a quantum amplitude,
enabling the succinct (and exact) representation of many quantum states through the combinatorial nature of this structure. A
DD can also be thought of as a homomorphic (lossless) compression scheme, since various manipulation operations for DDs exist
which implement any quantum gate operation, including measurement (without requiring decompression). Strong simulationis therefore easily implemented using a DD data structure [27–29].
Indeed, DD-based simulation was empirically shown to be competitive with state-of-the-art simulators [21, 28, 30] and is used inseveral simulator and circuit verification implementations [31, 32].
DDs and the stabilizer formalism are introduced in Sec. 2.

State space
poly-size LIMDD
poly-size
QMDD

Stabilizer
statespoly-size
MPS(pseudo) cluster states

Figure 1: The set of stabilizer
states and states representable
as poly-sized: (Pauli-)LIMDDs(this work), QMDDs and MPS.

QMDDs are currently the most succinct DD supporting quantum simulation, but in this paper
we show that they require exponential size to represent a type of stabilizer state called a cluster
state [33]. In order to unite the strengths of DDs and the stabilizer formalism and inspired by
SLOCC (Stochastic Local Operations and Classical Communication) equivalence of quantum states[34, 35], in Sec. 3, we propose LIMDD: a new DD for quantum computing simulation using local
invertible maps (LIMs). Specifically, LIMDDs eliminate the need to store multiple states which
are equivalent up to LIMs, allowing more succinct DD representations. For the local operations in
the LIMs, we choose Pauli operations, creating a Pauli-LIMDD, which we will simply refer to as
LIMDD. We prove that there is a family of quantum states —called pseudo cluster states— thatcan be represented by poly-sized (Pauli-)LIMDDs but that require exponentially-sized QMDDs andcannot be expressed in the stabilizer formalism. We also show the same separation for matrix
product states (MPS) [36, 38? ]. Fig. 1 visualizes the resulting separations.
Further, we give algorithms for simulating quantum circuits using Pauli-LIMDDs. We continue
by investigating the power of these algorithms compared to state-of-the-art simulation algorithms
based on QMDD, MPS and stabilizer decomposition. We find circuit families which Pauli-LIMDD
can efficiently simulate, which stands in stark contrast to the exponential space needed by QMDDbased, MPS-based and a stabilizer-decomposition-based simulator (the latter result is conditionalon the exponential time hypothesis). This is the first analytical comparison between decision
diagrams and matrix product states.
Efficient decision diagram operations for both classical [39] and quantum [2] applications cruciallyrely on dynamic programming (storing the result of each intermediate computation) and canonicity(each quantum state has a unique, smallest representative as a LIMDD) [40–42]. We provide algorithms for both in Sec. 4. Indeed, the main technical contribution of this paper is the formulationof a canonical form for Pauli-LIMDDs together with an algorithm which brings a Pauli-LIMDD
into this canonical form. By interleaving this algorithm with the circuit simulation algorithms, weensure that the algorithms act on LIMDDs that are canonical and as small as possible.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

3

The canonicity algorithm effectively determines whether two n-qubit quantum states |φ⟩ , |ψ⟩, each
represented by a LIMDD node φ, ψ, are equivalent up to a Pauli operator P , i.e, |φ⟩ = P |ψ⟩, which
we call an isomorphism between |φ⟩ and |ψ⟩. Here P = Pn ⊗ · · · ⊗ P1 consists of single qubit Pauli
operators Pi (ignoring scalars for now). In general, there are multiple choices for P , so the goal is
to make a deterministic selection among them, to ensure canonicity of the resulting LIMDD. To doso, we first take out one qubit and write the states as, e.g., |φ⟩ = c0 |0⟩ |φ0 ⟩ + c1 |1⟩ |φ1 ⟩ for complex
coefficients c0 , c1 . We then realize that Prest = Pn−1 · · · ⊗ P1 must map the pair (|φ0 ⟩ , |φ1 ⟩) toeither (|ψ0 ⟩ , |ψ1 ⟩) or (|ψ1 ⟩ , |ψ0 ⟩) (in case Pn is a diagonal or antidiagonal, respectively). Hence
Prest is a member of the intersection of the two sets of isomorphisms. Next, we realize that the setof all isomorphisms, e.g. mapping |φ0 ⟩ to |ψ0 ⟩, is a coset π · G of the stabilizer group G of |φ0 ⟩(i.e. the set of isomorphisms mapping |φ0 ⟩ to itself) where π is a single isomorphism |φ0 ⟩ → |ψ0 ⟩.
Thus, to find a (canonical) isomorphism between n-qubit states |φ⟩ → |ψ⟩ (or determine no suchisomorphism exists), we need three algorithms: to find (a) an isomorphism between (n − 1)-qubit
states, (b) the stabilizer group of an (n − 1)-qubit state (in fact: the group generators, which form
an efficient description), (c) the intersection of two cosets in the Pauli group (solving the Paulicoset intersection problem). Task (a) and (b) are naturally formulated as recursive algorithms
on the number of qubits, which invoke each other in the recursion step. For (c) we provide aseparate algorithm which first rotates the two cosets such that one is a member of the Pauli
Z group, hence isomorphic to a binary vector space, followed by using existing algorithms forbinary coset (hyperplane) intersection. Having characterized all isomorphisms |φ⟩ → |ψ⟩, we selectthe lexicographical minimum to ensure canonicity. We emphasize that the algorithm works forarbitrary quantum states, not only stabilizer states.

2 Preliminaries
Here, we briefly introduce two methods to manipulate and succinctly represent quantum states:
decision diagrams, which support universal quantum computing, and the stabilizer formalism, inwhich a subset of all quantum computations is supported which can however be efficiently classically simulated. Both support strong simulation, i.e. the probability distribution of measurement
outcomes can be computed (through weak simulation one only samples measurement outcomes).

2.1 Decision diagrams
An n-qubit quantum state |φ⟩ can be represented as a 2n -dimensional vector of complex numbers(modeling amplitudes) and can thus be described by a pseudo-Boolean function f : {0, 1}n → C
where
X
|φ⟩ =
f (xn , . . . , x1 ) |xn ⟩ ⊗ · · · ⊗ |x1 ⟩ .
(1)x1 ,...,xn ∈{0,1}

The Quantum Multi-valued Decision Diagram (QMDD) [27] is a data structure which can succinctly
represent functions of the form f : {0, 1}n → C, and thus can represent any quantum state per Eq. 1.
A QMDD is a rooted DAG with a unique leaf node 1 , representing the value 1. Fig. 2 (d) shows
an example (and its construction from a binary tree). Each node has two outgoing edges, called its
low edge (dashed line) and its high edge (solid line). The diagram has levels as each node is labeledwith (the index of) a variable; the root has index n, its children n − 1, etc, until the leaf withindex 0 (the set of nodes with index k form level k). Hence each path from root to leaf visits nodes
representing the variables x3 , x2 , x1 in order. The value f (xn , . . . , x1 ) = ⟨xn . . . x1 |φ⟩ is computedby traversing the diagram, starting at the root edge and then for each node at level i following thelow edge (dashed line) when xi = 0, and the high edge (solid line) when xi√= 1, while multiplying
the edge weights (shown in boxes) along the path, e.g., f (1, 1, 0) = 12 · 1 · − 2 · 1 = − √12 in Fig. 2.
A path from the root to a node v with index k (on level k) thus corresponds to a partial assignment(xn = an , . . . , xk−1 = ak−1 ), which induces subfunction f⃗a (xk , . . . , x1 ) ≜ f (an , . . . , ak−1 , xk , . . . , x1 ).
The node v represents this subfunction up to a complex factor γ, which is stored on the edge incoming to v along that path. This allows any two nodes which represent functions equal up
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

4

a)

b)

c)

1/2

d)

Level:
3

2

1
2

1
2

0
v

1

w
0

0
0

0

1
2

0

1
2

1

0 − √2 0

1
2

1
2

0

s

t

0 -

√1
2

- √1

1

√
− 2

0

2

u

u
0

0

0
1

1

1

Figure 2: Different decision diagrams representing the 3-qubit state [0, 0, 12 , 0, 21 , 0, − √12 , 0]⊤ , evolving into a
QMDD (right). Left, (a) shows the exponential binary tree, where a node on level i represents xi (see Eq. 1) and
its outgoing arrows xi = 0 (dashed) and xi = 1 (solid). The leaf contains the complex amplitude f (x1 , x2 , x3 )(see Eq. 1) for the assignment of (x1 , x2 , x3 ) corresponding to the path from the root, e.g. f (1, 1, 0) = − √12 .
Next (b), the leaves are merged by dividing out common factors, putting these as weights (shown in boxes) on
the edges going out of level-1 nodes (note in particular that we can suppress a separate 0 leaf, as 0 = 0·1). Thenthe same trick is applied to level-1 nodes in (c). In this example, all level-1 nodes v, w, s, t become isomorphic
and can be merged into a new node u, representing the vector |u⟩ = [1, 0]⊤ . This can be done because thelevel-1 nodes v, w, s, t respectively represent the vectors [0, 0]⊤ , [ 21 , 0]⊤ , [ 12 , 0]⊤ , [ √12 , 0]⊤ , which can be written

as c · |u⟩ = c · [1, 0]⊤ for respective weights c = 0, 12 , 12 , √12 . Finally, (d) shows the resulting QMDD, applying
the same tactic to nodes on levels 2 and 3. Note that a QMDD requires a root edge. Merging (isomorphic)
nodes makes QMDDs succinct. By convention, unlabelled edges have label 1.

to a complex factor to be merged. For instance, the node u on level 1 in Fig. 2 represents−1
f01 = f10 = √f = 0 · f00 . When all eligible nodes have been merged, the QMDD is reduced. A
2 11
reduced QMDD is a canonical representation: a given function has a unique reduced QMDD.
Canonicity ensures that the QMDD is always as small as possible as redundant nodes are merged.
But more importantly, canonicity allows for quick equality checks: two diagrams represent thesame state if and only if their root edges are the same (i.e., have the same label and point tothe same root node). This allows for efficient QMDD manipulation algorithms (i.e. updating the
QMDD upon performing a gate or measurement) through dynamic programming, which avoids
traversing all paths (exponentially many in the size of the diagram in the worst case). For all
quantum gates, there are algorithms to update the QMDD accordingly and measurement is also
supported (even efficiently). Therefore, QMDDs can simulate any quantum circuit, although they
may become exponentially large (in the number of qubits) already after applying part of the gatesfrom the circuit. The resulting simulator is strong, as the complete final state is computed as
QMDD (and computing measurement outcome probabilities on QMDD is tractable).
Finally, we can also define the semantics of a node v recursively, overloading Dirac notation: |v⟩.
For convenience, we denote an edge to node v labeled with ℓ pictographically asαℓv . Now a nodeβv0 and high edge
v1 , represents the state: |v⟩ ≜ α |0⟩ ⊗ |v ⟩ + β |1⟩ ⊗ |v ⟩,v with low edge
01where in the base case | 1 ⟩ ≜ 1 as defined above. We later define LIMDD semantics similarly.

2.2 Pauli operators and stabilizer states
In contrast to decision diagrams, the stabilizer formalism [6] forms a classically simulatable subset
of quantum computation. Instead of explicitly representing the (exponential) amplitude vector, thestabilizer formalism describes the symmetries a quantum state using so-called stabilizers. A unitary
operator U stabilizes a state |φ⟩ if |φ⟩ is a +1 eigenvector of U , i.e., U |φ⟩ = |φ⟩. The formalism
considers stabilizers U made up of the single-qubit Pauli operators Pauli ≜ {I, X, Y, Z} as defined
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

5below. In fact, a stabilizer is taken from the n-qubit Pauli group, defined as Paulin ≜ Pauli⊗n ,
i.e. it is the group generated by all n-qubit Pauli strings Pn ⊗ · · · ⊗ P1 with Pi ∈ Pauli. Here
we used the notation ⟨G⟩ = H to denote that G ⊆ H is a generator set for a group H. One
can check that Paulin = {ic Pn ⊗ · · · ⊗ P1 | P1 , . . . , Pn ∈ Pauli, c ∈ {0, 1, 2, 3}}, so in particular
we have Pauli1 = {±P, ±iP | P ∈ Pauli} (the Pauli set with a factor ±1 or ±i).

I ≜

1
0



0
0
,X ≜
1
1






1
0 −i
1 0
,Y ≜
,Z ≜
0
i 0
0 −1

The set of Pauli stabilizers Stab(|φ⟩) ⊂ Paulin of an n-qubit quantum state |φ⟩ necessarily forms
a subgroup of Paulin , since the identity operator I⊗n is a stabilizer of any n-qubit state andmoreover if U and V stabilize |φ⟩, then so do U V, V U and U −1 . Furthermore, any Pauli stabilizer
group G is abelian, i.e. A, B ∈ G implies AB = BA. The reason for this is that elements ofPaulin either commute (AB = BA) or anticommute (AB = −BA) under multiplication andanticommuting elements can never be stabilizers of the same state |φ⟩, because if A, B ∈ Stab(|φ⟩)
and AB = −BA then |φ⟩ = AB |φ⟩ = −(BA) |φ⟩ = − |φ⟩, a contradiction. Finally, note that−I⊗n can never be a stabilizer. In fact, these conditions are necessary and and sufficient: the classof abelian subgroups G of Paulin , not containing −I⊗n , are precisely all n-qubit stabilizer groups.
For clarity, we adopt the convention that we denote Pauli strings without phase using the symbols
P, Q, R, . . . and we use the symbols A, B, C, . . . for Pauli operators including phase; e.g., we may
write A = λP . The phase λ of any stabilizer λP ∈ Pauli can only be λ = ±1, derived as∀λP ∈ Stab(|φ⟩) :|φ⟩ = (λP ) |φ⟩ = (λP )2 |φ⟩ = λ2 I |φ⟩ = λ2 |φ⟩=⇒λ = ±1.

(2)

The number of generators k for a n-qubit stabilizer group S can range from 1 to n, and S has 2k
elements. If k = n, then there is only a single quantum state |φ⟩ (a single vector up to complex⊗nscalar) which is stabilized by S; such a state is called a stabilizer state. Equivalently, |φ⟩ = C |0⟩
where C is a circuit composed of only Clifford unitaries, a group generated by the Clifford gates:
1 0 0 0
1 1 1
1 0, (phase gate) S ≜, and CN OT ≜ 00 10 00 01 .
(Hadamard gate) H ≜ √0 −i
2 1 −10010
In the stabilizer formalism, an n-qubit stabilizer state is succinctly represented through n independent generators of its stabilizer group, each of which is represented by O(n) bits to encode
the Pauli string (plus factor), yielding O(n2 ) bits in total. Examples of (generators of) stabilizer
groups are ⟨Z⟩ for |0⟩ and ⟨X ⊗ X, Z ⊗ Z⟩ for √12 (|00⟩ + |11⟩). Updating a stabilizer state’s generators after application of a Clifford gate or a single-qubit computational-basis measurement can bedone in polynomial time in n [6, 7]. Various efficient algorithms exist for manipulating stabilizer(sub)groups S, including testing membership (is A ∈ Paulin a member of S?) and finding agenerating set of the intersection of two stabilizer (sub)groups. These algorithms predominantly
use standard linear algebra, e.g., Gauss-Jordan elimination, as described in App. A in detail.
In this work, we also consider states which are not stabilizer states and which therefore have anonmaximal stabilizer group (i.e. < n generators). To emphasize that a stabilizer group neednot be maximal, i.e. it is a subgroup of maximal stabilizer groups, we will use the term stabilizer
subgroup. Such objects are also studied in the context of simulating mixed states [43] and quantum
error correction [8]. Examples of stabilizer subgroups are {I} for √12 (|0⟩ + eiπ/4 |1⟩), ⟨−Z⟩ for |1⟩
and ⟨X ⊗ X⟩ for √16 (|00⟩ + |11⟩) + √13 (|01⟩ + |10⟩). In contrast to stabilizer states, in general astate is not uniquely identified by its stabilizer subgroup.
1
Graph states on n qubits are the output states of circuits with input state 2n/2(|0⟩+|1⟩)⊗n followed
by only CZ ≜ |00⟩⟨00| + |01⟩⟨01| + |10⟩⟨10| − |11⟩⟨11| gates, and form a strict subset of all stabilizer
states that is also important in error correction and measurement-based quantum computing [44].
By the (two-dimensional) cluster state on n2 qubits, we mean the graph state whose graph is then × n grid.

GivenPa vector space V ⊆ {0, 1}n and a length-n bitstring s, the corresponding coset state is√1x∈V |x + s⟩ where ‘+’ denotes bitwise xor-ing [45]. Each coset state is a stabilizer state.
|V |Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

6

Stabilizer decomposition-based methods [12–17] extend the stabilizer formalism to families of Clifford circuits with arbitrary input states |φn ⟩, enabling the simulation of universal quantum computation [46]. By decomposing the n-qubit state |φn ⟩ as linear combination of χ stabilizer states
followed by simulating the circuit on each of the χ stabilizer states, the measurement outcomescan be computed in time O(χ2 · poly(n)), where the least χ is referred to as the stabilizer rankof |φn ⟩. Therefore, stabilizer-rank based methods are efficient for any family of input states |φn ⟩whose stabilizer rank grows polynomially in n.
A specific method for obtaining a stabilizer decomposition of the output state of an n-qubit circuit
is by rewriting the circuit into Clifford gates and T = |0⟩⟨0|+eiπ/4 |1⟩⟨1| gates (a universal gate set).
Next, each of the T gates can be converted into an ancilla qubit initialized to the state T |+⟩ where|+⟩ = √12 (|0⟩ + |1⟩); thus, an n-qubit circuit containing t T gates will be converted into an n + t⊗n⊗t

qubit Clifford circuit with input state |φ⟩ = |0⟩ ⊗ (T |+⟩) [12]. We will refer to the resulting
specific stabilizer-rank based simulation method as the ‘Clifford + T simulator,’ whose simulation⊗truntime scales with χt = χ((T |+⟩) ), the number of stabilizer states in the decomposition of |φ⟩.
t
Trivially, we have χt ≤ 2 , and although recent work [12, 13] has found decompositions smaller
than 2t terms based on weak simulation methods, the scaling of χt remains exponential in t. Weemphasize that the Clifford + T decomposition is not necessarily optimal, in the sense that the⊗tintermediate states of the circuit might have lower stabilizer rank than |T ⟩ does. Consequently,
if a given circuit contains t = Ω(n) T -gates, then the Clifford + T simulator requires exponential
time (in n) for simulating this n-qubit circuit, even if there exist polynomially-large stabilizer
decompositions of each of the circuit’s intermediate and output states (i.e., in principle, there
might exist another stabilizer rank-based simulator that can simulate this circuit efficiently).

2.3 Matrix product states
Representing quantum states as matrix product states (MPS) has proven successful for solving alarge range of many-body physics problems [36, 37]. For qubits, an n-qubit MPS M is formallydefined as a series of 2n matrices Axk ∈ CDk ×Dk+1 where k ∈ [n], x ∈ {0, 1}, Dk ∈ N≥1 andD1 = Dn+1 = 1. Here, Dk+1 is the matrix dimension over the k-th bond. The interpretation |M ⟩is determined as ⟨x1 x2 . . . xn |M ⟩ = Ax1 1 Ax2 2 · · · Axnn for x1 , . . . , xn ∈ {0, 1}. If the bond dimension
may scale exponentially in the number of qubits, any family of quantum states can be represented
exactly by an MPS.
The Schmidt rank of a state |φ⟩ on n qubits, relative to a bipartition of the qubits into two sets
A and B, is the smallest integer m ≥ 1 such that |φ⟩ can be expressed as the superposition |φ⟩ =Pm
j=1 cj |aj ⟩A |bj ⟩B for complex coefficients cj , where the states |aj ⟩A (|bj ⟩B ) form an orthonormal
basis for the Hilbert space of the A register (B register). The relation with MPS is that the
maximum Schmidt rank with respect to any bipartition A = {x1 , . . . , xk }, B = {xk + 1, . . . , xn } isprecisely the smallest possible bond dimension Dk+1 required to exactly express a state in MPS.
Vidal [38] showed that MPS-based circuit simulation is possible in time O(n · poly(χ)) per elementary operation, where n is number of qubits and χ the maximum Schmidt rank for all intermediate
states computed.

3 Local Invertible Map Decision Diagrams
Sec. 3.1 introduces a LIMDD definition parameterized with different local operations. We mainly
consider the Pauli-LIMDD and refer to it simply as LIMDD. We show how LIMDDs generalize
QMDDs and can represent arbitrary quantum states, normalized or not. We then use this definitionin Sec. 3.2 to show how LIMDDs succinctly —i.e., in polynomial space— represent graph states (inparticular cluster states), coset states and, more generally, stabilizer states. On the other hand,
QMDDs and MPS require exponential size to represent two-dimensional cluster states.
We translate this exponential advantage in quantum state representation to (universal and strong)
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

7quantum circuit simulation in Sec. 3.3 by giving algorithms to update and query the LIMDD data
structure. These LIMDD manipulation algorithms take a LIMDD φ, representing some state |φ⟩,and return another LIMDD ψ that represents the state |ψ⟩ = U |φ⟩ for standard gates U and also
for arbitrary unitaries U (by preparing U in LIMDD form first; we show how). The measurement
algorithm we give returns the outcome in linear time in size of the LIMDD representation of thequantum state.
For many quantum operations, we show that our manipulation algorithms are efficient on all
quantum states, i.e., take polynomial time in the size of the LIMDD representation of the state.
Algorithms for certain other operations are efficient for certain classes of states, e.g., all Cliffordgates can be applied in polynomial time to a LIMDD representing a stabilizer state. We show that
LIMDDs can be exponentially faster than QMDDs, while they are never slower by more than amultiplicative factor O(n3 ). These algorithms use a canonical form of LIMDDs, such that for eachstate there is a unique LIMDD. We defer this subject to Sec. 4, which introduces reduced LIMDDs
and efficient algorithms to compute them.
With these algorithms, a quantum circuit simulator can be engineered by applying the circuit’s
gates one by one on the representation of the state as LIMDD. Prop. 1 provides the bottom line ofthis section by comparing simulator runtimes. In Sec. 3.4, we prove Prop. 1.
+T
Proposition 1. Let QSimClifford
denote the runtime of the Clifford + T simulator on circuit CC(allowing for weak simulation as in [13]). Let QSimD
C denote the runtime of strong simulation of
circuit C using method D = (Pauli−)LIMDD, QMDD, QMDD ∪ Stab, MPS, QMDD ∪ Stab.∗Here, the latter is an (imaginary) ideal combination of QMDD (not tractable for all Clifford circuits)
and the stabilizer formalism (tractable for Clifford circuits), i.e., one that always inherits the best
worst-case runtime from either method.
The following holds, where Ω∗ discards polynomial factors, i.e., Ω∗ (f (n)) ≜ Ω(nO(1) f (n)).

There is a family of circuits C such that:1.

Clifford + T
LIMDD is exponentially faster than Clifford + T : QSimC= Ω∗ (2n · QSimLIMDD
),†
C

2.

LIMDD is exponentially faster than MPS: QSimMPS= Ω∗ (2n · QSimLIMDD), andCC

3.

LIMDD is exponentially faster than QMDD: QSimQMDD= Ω∗ (2n · QSimLIMDD).
CC

4. For all C, LIMDD is at worst cubically slower than QMDD: QSimLIMDD= O(n3 · QSimQMDD).
CC
5. Item 3 and 4 hold when replacing QMDD with QMDD ∪ Stab.

3.1 The LIMDD data structure
Where QMDDs only merge nodes representing the same complex vector up to a constant factor,
the LIMDD data structure goes further by also merging nodes that are equivalent up to local
operations, called Local Invertible Maps (LIMs) (see Def. 1). As a result, LIMDDs can be exponentially more succinct than QMDDs, for example in the case of stabilizer states (see Sec. 3.2). We
will call nodes which are equivalent under LIMs, (LIM-) isomorphic. This definition generalizes
SLOCC equivalence (Stochastic Local Operations and Classical Communication); if we choose theparameter G to be the linear group, then the two notions coincide (see [34, App. A] and [35, 47]).
Definition 1 (G-LIM, G-Isomorphism). An n-qubit G-Local Invertible Map (LIM) is an operator O of the form O = λOn ⊗ · · · ⊗ O1 , where G is a group of invertible 2 × 2 matrices, Oi ∈ G
and λ ∈ C \ {0}. A G-isomorphism between two n-qubit quantum states |φ⟩ , |ψ⟩ is a LIM O suchthat O |φ⟩ = |ψ⟩, denoted |φ⟩ ≃G |ψ⟩. Note that G-isomorphism is an equivalence relation.
∗ We are not aware of any (potentially better) weak D-based simulation approaches and do not consider them.
† Assuming the exponential time hypothesis (ETH). See Sec. 3.4.3 for details.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

8

1
4

a)

1
4

b)

q3′
ω

-ω
q2

q2′

p2
i

p′2
i

-ω

-ω
q1′

q1

1

−1

q2

q2′

p2
i

q3
ω

-ω

Z -ω Zℓ1

1iZ

Level: 4

−1

q3′

q3

1
·I⊗I⊗I⊗X
4

d)

−1

-1
q3

1
4

c)

p′2

-ωT ZX ⊗ I

I⊗Z
i·I

-ωZ

q3′

ℓ′2

ℓ2

3

ℓ3

ωT ⊗ X

I⊗I
ℓ2

Z⊗I⊗Z

ω·T ⊗Z
i·I

Z

ℓ′2
Z

ℓ1

ℓ1

1

1

2

1

Figure 3: A QMDD (a) representing the state 14 [1, 1, i, i, −ω, ω, i, i, −1, 1, −i, i, −ω, −ω, i, −i]⊤ with ω = eiπ/4 ,evolving into a LIMDD (d). As in Fig. 2, diagram nodes are horizontally ordered in ‘levels’ with qubit indices4, 3, 2, 1. Low edges are dashed, high edges solid. See the text for an explanation.
By convention, unlabelled edges have label 1 (for QMDD) or I⊗k (for LIMDD nodes at level k).

We define PauliLIMn ≜ ⟨Pauli⟩-LIM, i.e., the group of Pauli operators P ∈ Paulin with
arbitrary complex factor λ ∈ C \ {0} (λ can absorb the factor γ = ±1, ±i in P = γPn ⊗ · · · ⊗ P1 .
Note λ = ±1 still for PauliLIMn operators which are stabilizers, by eq. (2)).
Before we give the formal definition of LIMDDs in Def. 2, we give a motivating example in Fig. 3,
which uses ⟨X, Y, Z, T ⟩-LIMs to demonstrate how the use of isomorphisms can yield small diagrams
for a four-qubit state. This figure shows how to merge nodes in four steps, shown in subfigures(a)-(d), starting with a large QMDD (a) and ending with a small LIMDD (d). In the QMDD (a),⊤⊤the nodes labeled q1 and q1′ represent the single-qubit states |q1 ⟩ = [1, 1] and |q1′ ⟩ = [1, −1] ,′respectively. By noticing that these two vectors are related via |q1 ⟩ = Z |q1 ⟩, we merge nodes q1 , q1′
into node ℓ1 in (b), storing the isomorphism Z on all incoming edges that previously pointed to
q1′ . From step (b) to (c), we first merge q2 , q2′ into ℓ2 , observing that |q2′ ⟩ = I ⊗ Z |q2 ⟩. Second,
we create a node ℓ′2 such that |p2 ⟩ = T ZX ⊗ I |ℓ′2 ⟩ and |p′2 ⟩ = T ⊗ X |ℓ′2 ⟩. So we can merge nodes
p2 , p′2 into ℓ′2 , placing these isomorphisms on the respective edges. To go from (c) to (d), we merge
nodes q3 , q3′ into node ℓ3 by noticing that |q3′ ⟩ = (Z ⊗ I ⊗ Z) |q3 ⟩. This isomorphism Z ⊗ I ⊗ Z isstored on the high edge out of the root node. We have |q3 ⟩ = I ⊗ I ⊗ X |ℓ3 ⟩, so we propagate theisomorphism I ⊗ I ⊗ X upward, and store it on the root edge. Therefore, the final LIMDD has the
LIM 41 I ⊗ I ⊗ I ⊗ X on its root edge.
The resulting data structure in Fig. 3 is a LIMDD of only six nodes instead of ten, but requires
additional storage for the LIMs. Sec. 3.2 shows that merging isomorphic nodes sometimes leads toexponentially smaller diagrams, while the additional cost of storing the isomorphisms results only
costs a linear factor of space (linear in the number of qubits).
The transformation presented above (for Fig. 3) only considers particular choices for LIMs. For
instance, it would be equally valid to select LIM I ⊗ Z instead of −I ⊗ XZ for mapping q2′ onto
q2 . In fact, efficient algorithms to select LIMs in such a way that a canonical LIMDD is obtained
are a cornerstone for the LIMDD manipulation algorithms presented in Sec. 3.3. Sec. 4 provides asolution for ⟨Pauli⟩-LIMs (the basis for all results presented in the current article), which is basedon using the stabilizers of each node, e.g., the group generated by {I ⊗ X, Y ⊗ I} for q2 .
Definition 2. An n-G-LIMDD is a rooted, directed acyclic graph (DAG) representing an n-qubit
quantum state. Formally, it is a 6-tuple (Node ∪ {Leaf}, idx, low, high, label, eroot ), where:• Leaf (a sink) is a unique leaf node with qubit index idx(Leaf) = 0;
• Node is a set of nodes with qubit indices idx(v) ∈ {1, . . . , n} for v ∈ Node;Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

9• eroot is a root edge without source pointing to the root node r ∈ Node with idx(r) = n;
• low, high : Node → Node ∪ {Leaf} indicate the low and high edge functions, respectively.
We write lowv (or highv ) to obtain the edge (v, w) with w = low(v) (or w = high(v)). For allv ∈ Node it holds that idx(low(v)) = idx(high(v)) = idx(v) − 1 (no qubits are skipped‡ );• label : low∪high∪{eroot } → k−G-LIM∪{0} is a function labeling edges ( . , w) with k-G-LIMs
or 0, where k = idx(w)
A

B

u
w
We will find it convenient to write v
for a node u with low and high edges to nodesAv for a (root) edge to v
v and w labeled with A and B, respectively. We will also denotev, the LIM should be interpreted as I⊗idx(v) .
labeled with A. When omitting A or B, e.g.,

We define the semantics of a leaf, node v and an edge e to node v by overloading the Dirac notation:|Leaf⟩ ≜ 1|e⟩ ≜ label(e) · |v⟩|v⟩ ≜ |0⟩ ⊗ |lowv ⟩ + |1⟩ ⊗ |highv ⟩It follows from this definition that a node v with idx(v) = k represents a quantum state on k qubits.
This state is however not necessarily normalized: For instance, a normalized state α |0⟩ + β |1⟩,
α

β

β/α

α

v
v
1 or a LIMDD
1
1 with root edgev .
can be represented as a LIMDD 1
So the node v represents a state up to global scalar. But, in general, any scalar can be applied to
the root edge, or any other edge for that matter. So LIMDDs can represent any complex vector.
A
Bv and e1w can be
The tensor product |e0 ⟩ ⊗ |e1 ⟩ of the G-LIMDDs with root edges e0α1computed just like for QMDDs [27]: Take all edges
pointing to the leaf in the LIMDD e0
α·B

w pointing to the e1 root node w. The result is an n + m
and replace them with edges
level LIMDD if e0 has n levels and e1 has m. In addition, the LIMs C on the other edges in the
LIMDD e0 should be extended to C ⊗ I⊗m .

We can now consider various instantiations of the above general LIMDD definition for different
LIM groups G. A G-LIMDD with G = {I} yields precisely all QMDDs by definition, i.e., alledges labels effectively only contain scalars. As all groups G contain the identity operator I, theuniversality of G-LIMDDs (i.e., all quantum states can be represented) follows from the universality
of QMDDs. It also follows that any state that can efficiently be represented by QMDD, can be
efficiently represented by a G-LIMDD for any G. Similarly, we can consider ⟨Z⟩ and ⟨X⟩, which aresubgroups of the Pauli group, and define a ⟨Z⟩-LIMDD and a ⟨X⟩-LIMDD; instances that we will
study for their relation to graph states and coset states in Sec. 3.2. Finally, and most importantly,⟨Pauli⟩-LIMDDs can represent all stabilizer states in polynomial space, which is a feature thatneither QMDDs nor matrix product states (MPS) posses, as shown in Sec. 3.2.
In what follows, we only consider ⟨X⟩-, ⟨Z⟩-, and ⟨Pauli⟩-LIMDDs, or Pauli-LIMDD for short.
For Pauli-LIMDDs, we now illustrate how to find the amplitude of a computational basis state⟨x|ψ⟩ for a bitstring x ∈ {0, 1}n by traversing the LIMDD of the state |ψ⟩ from root to leaf, as
follows. Suppose that this diagram’s root edge eroot points to node r and is labeled with the LIM
label(eroot ) = A = λPn ⊗ · · · ⊗ P1 ∈ PauliLIMn . First, we substitute |r⟩ = |0⟩ |lowr ⟩ + |1⟩ |highr ⟩,where lowr , highr are the low and high edges going out of r, thus obtaining ⟨x|ψ⟩ = ⟨x|eroot ⟩ =⟨x| A (|0⟩ |lowr ⟩ + |1⟩ |highr ⟩). Next, we notice that ⟨x| A = λ(⟨xn | Pn ) ⊗ · · · ⊗ (⟨x1 | P1 ) = γ ⟨y| forsome γ ∈ C and a computational basis state ⟨y|. Therefore, letting y ′ = yn−1 . . . y1 , it suffices tocompute ⟨yn | ⟨y ′ | (|0⟩ |lowr ⟩ + |1⟩ |highr ⟩), which reduces to computing either ⟨y ′ |lowr ⟩ if yn = 0, or⟨y ′ |highr ⟩ if yn = 1. By applying this simple rule repeatedly, one walks from the root to the leaf,
encountering one node on each level. The amplitude ⟨x|ψ⟩ is then found by multiplying together
the scalars γ found along this path. Alg. 1 formalizes this. Its runtime is O(n2 ).
‡ Decision diagram definitions [18, 19, 48] often allow to skip (qubit) variables, interpreting them as ‘don’t cares.’
We disallow this here, since it complicates definitions and proofs, while at best it yields linear size reductions [41].

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

10

Algorithm 1 Read the amplitude for basis state |xn . . . x1 ⟩ from n-qubit state |e⟩ = A · |v⟩ withA = λPn ⊗ · · · ⊗ P1 ∈ PauliLIMn .
1: procedure ReadAmplitude(Edge eAv , xn , . . . , x1 ∈ {0, 1} with n = idx(v))

2:

if n = 0 then return λ3:γ ⟨yn . . . y1 | := ⟨xn . . . x1 | λPn ⊗ · · · ⊗ P1

4:

if yn = 0 then

▷ yn = 0

5:6:

return γ · ReadAmplitude(lowv , yn−1 , . . . , y1 )else

▷ yn = 1

7:

return γ · ReadAmplitude(highv , yn−1 , . . . , y1 )▷ O(n)-computable LIM operation

3.2 Succinctness of LIMDDs
Succinctness is crucial for efficient simulation, as we show later. In this section, we show exponential advantages for representing states with LIMDDs over two other state-of-the-art data structures:
QMDDs and Matrix Product States (MPS) [36, 37]. Specifically, QMDDs and MPS require exponential space in the number of qubits to represent specific stabilizer states called (two-dimensional)
cluster states. We also show that an ad-hoc combination of QMDD with the stabilizer formalism
still requires exponential space for ‘pseudo-cluster states.’ These results are visualized in Fig. 1.
3.2.1

LIMDDs are exponentially more succinct than QMDDs (union stabilizer states)

Fig. 4 visualizes succinctness relations between
different quantum state representations, as
proved in Prop. 2. In particular, G-LIMDDswith G = ⟨Pauli⟩ can be exponentially more
succinct than QMDDs, and retain this exponential advantage even with G = ⟨Z⟩ , ⟨X⟩.
In Corollary 3, we show the strongest result,namely that LIMDDs are also more succinct thanthe union of QMDDs and stabilizer states, written QMDD ∪ Stab, which can be thought ofa structure that switches between QMDD andthe stabilizer formalism depending on its content (stabilizer or non-stabilizer state). This
demonstrates that ad-hoc combinations of existing formalisms do not make LIMDDs obsolete.
Proposition 2. The inclusions and separationsin Fig. 4 hold.

Pauli-LIMDD

QMDD ∪ Stab
⟨X⟩-LIMDD

stabilizer states⟨Z⟩-LIMDD
graph states

coset states
QMDD

2D cluster states

Figure 4: Relations between non-universal classes ofquantum states (gray) and decision diagrams, where
we consider a diagram as the set of states that it canrepresent in polynomial size. Solid arrows denote set
inclusion. Dashed arrows D1 99K D2 signify an exponential separation between two classes, i.e., somequantum states have polynomial-size representationin D1 , but only exponential-size in D2 .
By transitivity, QMDD is exponentially separated
from all representations (not drawn for clarity).

Proof. The inclusions between the sets of states
shown in gray are well known [44, 45]. The inclusions between decision diagrams hold because,e.g., a QMDD is a G-LIMDD with G = {I}, i.e., each label is of the form λIn with λ ∈ C, as
discussed in Sec. 3.1. The relations between coset, graph, stabilizer states and G-LIMDD withG = ⟨X⟩ , ⟨Z⟩ , ⟨Pauli⟩ are proven in Th. 1 and App. C (which also shows that poly-sized LIMDD
includes QMDD ∪ Stab). Corollary 3 shows that there is family of a non-stabilizer states (withsmall LIMDD) for which QMDD is exponential, hence the separation between QMDD ∪ Stab. Th. 2shows the separation√ with QMDDs by demonstrating that the so-called (two-dimensional) cluster
state, requires 2Ω( n) nodes as QMDD. Finally, App. C proves the same for coset states.
Th. 1 shows that any stabilizer state can be represented as a ⟨Pauli⟩-Tower LIMDD (Def. 3).

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

11

Definition 3. A n-qubit G-Tower-LIMDD, is a G-LIMDD with exactly one node on each level.
Edges to nodes on level k are labeled as follows: low edges are labeled with I⊗k , high edges withP ∈ G⊗k ∪ {0} and the root edge is labeled with λ · P with P ∈ G⊗k and λ ∈ C \ {0} (i.e.,
in contrast to high edges, the root edge can have an arbitrary scalar). Fig. 6 depicts a n-qubit
G-Tower LIMDD.
Theorem 1. Let n > 0. Each n-qubit stabilizer state is represented up to normalization by a⟨Pauli⟩-Tower LIMDDs of Def. 3, e.g., where the scalars λ of the PauliLIMs λP on high edges
are restricted as λ ∈ {0, ±1, ±i}. Conversely, every such LIMDD represents a stabilizer state.
Proof sketch of Th. 1. (Full proof in App. C) The n = 1 case: the six single-qubit states |0⟩ , |1⟩ , |0⟩±|1⟩ and |0⟩ ± i |1⟩ are all represented by a ⟨Pauli⟩-Tower LIMDD with a single node on top of
the leaf. The induction step: Let |ψ⟩ be an n-qubit stabilizer state. First, consider the casethat |ψ⟩ = |a⟩ |ψ ′ ⟩ where |a⟩ = α |0⟩ + β |1⟩ (with α, β ∈ {0, ±1, ±i}) and |ψ ′ ⟩ are stabilizer
states on respectively 1 and n − 1 qubits. Then |ψ⟩ is represented by the ⟨Pauli⟩-Tower-LIMDDβI

αI. In the remaining case, |ψ⟩ = √12 (|0⟩ |ψ0 ⟩ + |1⟩ |ψ1 ⟩), where both |ψ0 ⟩ and |ψ1 ⟩are stabilizer states. Moreover, since |ψ⟩ is a stabilizer state, there is always a set of single-qubit
Pauli gates P1 , . . . , Pn and a λ ∈ {±1, ±i} such that |ψ1 ⟩ = λPn ⊗ · · · ⊗ P1 |ψ0 ⟩. That is, in our
terminology, the states |ψ0 ⟩ and |ψ1 ⟩ are isomorphic. Hence |ψ⟩ can be written as
ψ′

ψ′

1
|ψ⟩ = √ [|0⟩ |ψ0 ⟩ + λ |1⟩ ⊗ (Pn ⊗ · · · ⊗ P1 |ψ0 ⟩)]
2
I

(3)

λPn ⊗ · · · ⊗ P1

ψ0Hence |ψ⟩ is represented by the Tower Pauli-LIMDD ψ0. In both cases, |ψ ′ ⟩is represented by a Tower Pauli-LIMDDs (up to normalization) by the induction hypothesis.
√
1/ 2 · I⊗3√
1/ 8 · I⊗3

eGHZ

e+++

√
1/ 8 · I⊗3eφ

λ · LnLevel n
Ln−1

X ⊗X

X ⊗X

Level n − 1

v
0

Y

0

−1

..
.
Level 11

11Figure 5:
Example ⟨Pauli⟩-Tower LIMDDs for three stabilizer states:
the GHZ state |eGHZ ⟩ = √12 (|000⟩ + |111⟩), for|e+++ ⟩ = |+ + +⟩ where |+⟩ = √12 (|0⟩ + |1⟩), and the state |eφ ⟩ =√1 (|000⟩ − |001⟩ + i |010⟩ + i |011⟩ + i |100⟩ + i |101⟩ − |110⟩ + |111⟩)8with stabilizer group generators {X ⊗ X ⊗ X, −Z ⊗ Z ⊗ X, I ⊗ Y ⊗ Z}.

L0
1

Figure 6: An n-qubit GTower LIMDD. We let Li ∈G⊗i ∪ {0} and λ ∈ C \ {0}(only root edges have an arbitrary scalar).

We stress that obtaining the LIMs for the Pauli Tower-LIMDD of a stabilizer state is not immediatefrom the stabilizer generators; specifically, the edge labels in the Pauli-LIMDD are not directly thestabilizers of the state. For example, the GHZ state √12 (|000⟩ + |111⟩) is represented by |eGHZ ⟩ =
√1
2

I

X⊗X

with |v⟩ = |00⟩ in Fig. 5, but X ⊗ X is not a stabilizer of |00⟩. Nonetheless,Th. 1 implicitly contains an algorithm that constructs a ⟨Pauli⟩-Tower LIMDD stabilizer state.
Sec. 4 also provides the inverse construction, which we use to make LIMDDs (representing any
quantum state) canonical in time O(mn3 ) (using Alg. 3).
vv

We also note that Th. 1 demonstrates that for any n-qubit stabilizer state |φ⟩, the (n − 1)-qubit
states (⟨0| ⊗ I2n−1 ) |φ⟩ and (⟨1| ⊗ I2n−1 ) |φ⟩ are not only stabilizer states, but also PauliLIMAccepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

12

isomorphic. While we believe this fact is known in the community,§ we have not found thisstatement written down explicitly in the literature. More importantly for this work, to the bestof our knowledge, the resulting recursive structure (which DDs are) has not yet been exploited in
the context of classical simulation.
Next, Th. 2 shows the separation
with QMDDs by demonstrating that the so-called (two-dimensional)√
cluster state, requires 2Ω( n) nodes as QMDD. Corollary 3 shows that a trivial combination withstabilizer formalism does not solve this issue.
Theorem 2. Denote by |Gn ⟩ the two-dimensional cluster state, defined as a graph state on then × n lattice. Each QMDD representing |Gn ⟩ has at least 2⌊n/12⌋ nodes.
Proof sketch. Consider a partition of the vertices of the n × n lattice into two sets S and T of
size 12 n2 , corresponding to the first 12 n2 qubits under some variable order. Then there are at least⌊n/3⌋ vertices in S that are adjacent to a vertex in T [49, Th. 11]. Because the degree of thevertices is small, many vertices on this boundary are not connected and therefore influence theamplitude function independently of one another. From this independence, it follows that, for
1 2nany variable order, the partial assignments ⃗a ∈ {0, 1} 2 induce 2⌊n/12⌋ different subfunctions f⃗a ,2nwhere f : {0, 1} → C is the amplitude function of |Gn ⟩. The lemma follows by noting that a
QMDD has a single node per unique subfunction modulo phase. For details see App. B.
Corollary 3 (Exponential separation between Pauli-LIMDD versus QMDD union stabilizer states).
There is a family of non-stabilizer states, which we call pseudo cluster states, that have polynomialsize Pauli-LIMDD but exponential-size QMDDs representation.
Proof. Consider the pseudo cluster state |φ⟩ = √12 (|0⟩ + eiπ/4 |1⟩) ⊗ |Gn ⟩ where |Gn ⟩ is the graphstate on the n × n grid. This is not a stabilizer state, because each computational-basis coefficientof a stabilizer state is of the form z · √1 k for z ∈ {±1, ±i} and some integer k ≥ 1 [9], while2 n2 +1⊗n2
⟨1| ⊗ ⟨0||φ⟩ = eiπ/4 · √12
is not of this form. Its canonical QMDD and Pauli-LIMDD
1

eiπ/4

I

eiπ/4 I

Gn
Gn
have root nodes Gn
and Gn, where the respective diagram for Gn isexponentially large (Th. 2) and polynomially small (Th. 1).

3.2.2

LIMDDs are exponentially more succinct than matrix product states

Th. 4 states that matrix product states (MPS) require large bond dimension for representing thetwo-dimensional cluster states, which follows directly from the well-known results that these states
have large Schmidt rank.
Theorem 4. To represent the graph state on the n × n grid (the two-dimensional cluster state onn2 qubits), an MPS requires bond dimension 2Ω(n) .
Proof. Van den Nest et al. [50] consider spanning trees over the complete graph where each nodecorresponds to a qubit and define the Schmidt-rank width: the largest encountered base-2 logarithmof the Schmidt rank between the two connected components resulting from removing an edge fromthe spanning tree, minimized over all possible spanning trees. It then follows from the relation
between bond dimension and Schmidt rank (see Sec. 2) that any quantum state with Schmidt-rank
width w requires bond dimension 2w for representation by an MPS. Van den Nest et al. also showedthat for graph states, the Schmidt-rank width equals the so-called rank width of the graph, which
for n × n grid graphs was shown to equal n − 1 by Jelinek [51]. This proves the theorem.
§ For instance, this fact can be observed (excluding global scalars) by executing the original algorithm forsimulating single-qubit computational-basis measurement on the first qubit, as observed in [6]. Similarly, thecharacterization in Prop. 2 of ⟨Z⟩-Tower-LIMDDs as representing precisely the graph states, is immediate by defining
graph states recursively (see App. C). The fact that ⟨X⟩-Tower LIMDDs represent coset states is less evident andrequires a separate proof, which we also give in App. C.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

13

In contrast, the Pauli-LIMDD efficiently represents cluster states, and more generally all stabilizer
states (Th. 1).

3.3 Pauli-LIMDD manipulation algorithms for simulation of quantum computing
In this section, we give all algorithms that are necessary to simulate a quantum circuit with PauliLIMDDs (referred to simply as LIMDD from now on). We provide algorithms which update the
LIMDD after an arbitrary gate and after a single-qubit measurement in the computational basis.
In addition, we give efficient specialized algorithms for applying a Clifford gate to a stabilizer state(represented by a ⟨Pauli⟩-Tower LIMDD) and computing a measurement outcome. We also showthat many (Clifford) gates can in fact be applied to an arbitrary state in polynomial time. Table 1
provides an overview of the LIMDD algorithms and their complexities compared to QMDDs.
Central to the speed of many DD algorithms is keeping the diagram canonical throughout thecomputation. Recall from Sec. 3.1, that a G-LIMDD can merge isomorphic nodes v ≃G w, i.e.,
if there exists a G-LIM C such that |w⟩ = C |v⟩. To achieve this, we require a ‘MakeEdge’
A

B

C

v1
w
v with C |v⟩ = |w⟩, where v is the, returns
subroutine which, given the node v0
unique, canonical node in the diagram that is G-isomorphic to node w. Sec. 4.2 provides a O(n3 )MakeEdge algorithm for ⟨Pauli⟩-LIMDDs satisfying this specification. For now, the reader may
assume the provisional implementation in Alg. 2, which does not yet merge LIM-isomorphic nodes
and hence does not yield canonical diagrams.

In line with other existing efficient decision-diagram algorithms, we use dynamic programming in
our algorithms to avoid traversing all paths (possibly exponentially many) in the LIMDD. In thisapproach, the decision diagram is manipulated and queried using recursive algorithms, which store
intermediate results for each recursive call to avoid unnecessary recomputations. For instance,
Alg. 3 makes any LIMDD canonical using dynamic programming and the (real) O(n3 ) MakeEdge
algorithm from Sec. 4.2. It recursively traverses child nodes at Line 3, reconstructing the diagram
bottom up in the backtrack at Line 4. By virtue of dynamic programming it visits each node only
once: The table CanonicalCache : Node → Edge stores for each node its canonical counterpart
as soon as it is computed at Line 4. The algorithm therefore runs in time O(n3 m) where m is thenumber of nodes in the original diagram.
This recursive algorithmic structure that uses dynamic programming and reconstructs the diagram in the backtrack, is typical for all DD manipulation algorithms. Note that constant-timecache lookups (using a hash table) therefore require the canonical nodes produced by MakeEdge.
LIMDDs additionally require the addition of LIMs to the caches; Sec. 3.3.3 shows how we do this.
Table 1: Worst-case complexity of currently best-known algorithms for applying specific operations, in terms of
the size of the input diagram size m (i.e., the number of nodes in the DD) and the number of qubits n. Althoughaddition (Add) of quantum states is not, strictly speaking, a quantum operation, we include it because it is asubroutine of gate application. Note that several of the LIMDD algorithms invoke MakeEdge and therefore
inherit its cubic complexity (as a factor).

Operation \ input:
Single |0⟩ / |1⟩-basis measurement
Single Pauli gate
Single Hadamard gate / Add()
Clifford gate on stabilizer state
Multi-qubit gate
MakeEdge
Checking state equality

QMDD
O(m)
O(m)
O(2n ) ¶O(2n )O(4n )
O(1)
O(1)

LIMDD
O(m)
O(1)
O(n3 2n ) ¶O(n4 )O(n3 4n )O(n3 )O(n3 )

SectionSec. 3.3.1
Sec. 3.3.2
Sec. 3.3.2
Sec. 3.3.4
Sec. 3.3.3
Sec. 4.2
Sec. 4.2.2¶ The worst-case of QMDDs and LIMDDs is caused by the vector addition introduced by the Hadamard gate [52,Table 2, +BC, +SLDD]. See Fig. 9 for an example.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

14

Algorithm 2 Provisionary algorithm MakeEdge for creating a new node/edge. Given two edges
representing states A |v⟩ , B |w⟩, it returns an edge representing the state |0⟩ A |v⟩ + |1⟩ B |w⟩. The
real MakeEdge algorithm (Sec. 4.2) returns a canonical node, assuming v, w are already canonical.
1: procedure MakeEdge(Edge
A

2:

u :=

3:

return Edge

v

u

B

A

v , Edge

B

w)

w
I⊗k

▷ Where k = idx(u)uAlgorithm 3 Make any LIMDD canonical using MakeEdge.
1: procedure MakeCanonical(Edge

A

v )

2:
3:
4:

if v ∈/ CanonicalCache then
▷ Compute result once for v and store in cache:e0 , e1 := MakeCanonical(low(v)), MakeCanonical(high(v))
CanonicalCache[v] := MakeEdge(e0 , e1 )5:

return A · CanonicalCache[v]▷ Retrieve result from cache

Finally, in this section, we often decompose LIMS using A = λPn ⊗ P ′ . Here λ ∈ C is a nonzero scalar, P ′ a Pauli string and Pn ∈ {I, X, Z, Y } = Pauli. Our algorithms will use the Follow
procedure from Alg. 4 to easily navigate diagrams according to edge semantics. Provided with a bit
string xn . . . x1 , the procedure is the same as ReadAmplitude. If however fewer bits are supplied,
it returns a LIMDD root edge representing a subvector. For instance, the subvector for |10⟩ of the1
4 I ⊗ XZ

ℓ2 ⟩ =
LIMDD root edge er in Fig. 3 (d) is computed by taking |follow10 (er )⟩ = |
1
n−ℓ
·
[−1,
1,
−i,
i].
So,
we
can
specifyitas|follow
(e)⟩
=
(⟨b|
⊗
I
)
|e⟩,
i.e.,
select
the
bth
block
b
4
of size 2n−ℓ from the vector |e⟩ (or rather, return a LIMDD edge representing that block).

Algorithm 4 Follow: a generalization of ReadAmplitude, returning edges.
1: procedure Follow(Edge e

λPn ⊗ · · · ⊗ P1

λPn ⊗ · · · ⊗ P1

v , xn , . . . , xk ∈ {0, 1} with n = idx(v) and k ≥ 1)

2:

if k > n then returnv▷ End of bit string3:

γ ⟨yn . . . yk | := ⟨xn . . . xk | λPn ⊗ · · · ⊗ Pk

▷ O(n)-computable LIM operation

4:

if yn = 0 then

▷ yn = 0

5:6:

return γ · Follow(lowv , yn−1 , . . . , yk )else

▷ yn = 1

7:

return γ · Follow(highv , yn−1 , . . . , yk )3.3.1

Performing a measurement in the computational basis

We discuss algorithms for measuring, sampling and updating after measurement of the top qubit.
App. E gives general algorithms with the same worst-case runtimes.
The procedure MeasurementProbability in Alg. 5 computes the probability p of observing theoutcome |0⟩ for state |e⟩. If the quantum state can be written as |e⟩ = |0⟩ |e0 ⟩ + |1⟩ |e1 ⟩, then theprobability is p = ⟨e0 |e0 ⟩ / ⟨e|e⟩, where we have ⟨e|e⟩ = ⟨e0 |e0 ⟩ + ⟨e1 |e1 ⟩. Hence we compute thesquared norms of ex = followx (e) using the SquaredNorm subroutine. The total runtime isdominated by the subroutine SquaredNorm, which computes the quantity ⟨e|e⟩ given a LIMDDv by traversing the entire LIMDD. We have ⟨e|e⟩ = |λ|2 ⟨v| P † P |v⟩ = |λ|2 ⟨v|v⟩,
edge e =†because P P = I for Pauli matrices. Therefore, to this end, it computes the squared norm of |v⟩.
Since ⟨v|v⟩ = ⟨lowv |lowv ⟩+⟨highv |highv ⟩, this is accomplished by recursively computing the squared
norm of the node’s low and high edges. This subroutine visits each node at most once by virtue ofdynamic programming, which stores intermediate results in a cache SNormCache : Node → R
for all recursive calls (Line 7, 8). Therefore, it runs in time O(m) for a diagram with m nodes.
λP

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

15

Algorithm 5 Algorithms MeasurementProbability and UpdatePostMeas for respectively
computing the probability of observing outcome |0⟩ when measuring the top qubit of a Pauli
LIMDD in the computational basis and converting the LIMDD to the post-measurement state after
outcome m ∈ {0, 1}. The subroutine SquaredNorm takes as input a Pauli LIMDD edge e, andreturns ⟨e|e⟩. It uses a cache to store the value s of a node v.
1: procedure MeasurementProbability(Edge e)2:
s0 := SquaredNorm(follow0 (e))
s1 := SquaredNorm(follow1 (e))
3:4:return s0 /(s0 + s1 )5: procedure SquaredNorm(Edge
6:
7:
8:

λP

v

with λ ∈ C, P ∈ Pauliidx(v) )if idx(v) = 0 then return |λ|2
if v ∈/ SNormCache thenI
SNormCache[v] := SquaredNorm(follow0 (▷ Compute result once for v and store in cache:Iv ))v )) + SquaredNorm(follow1 (return |λ|2 · SNormCache[v]▷ Retrieve result for v from cache and multiply with |λ|2λP
v , measurement outcome m ∈ {0, 1})
10: procedure UpdatePostMeas(Edge e
11:
if m = 0 then12:
er := MakeEdge(follow0 (e), 0 · follow0 (e))
else
13:
14:
er := MakeEdge(0 · follow0 (e), follow1 (e))
p15:
return 1/ SquaredNorm(er ) · er
9:The outcome m ∈ {0, 1} can then be chosen by flipping a p-biased coin. The corresponding
state update is implemented by the procedure UpdatePostMeas. In order to update the state|e⟩ = |0⟩ |e0 ⟩+|1⟩ |e1 ⟩ after the top qubit is measured to be m, we simply construct an edge |m⟩ |em ⟩using the MakeEdge subroutine. This state is finally normalized by multiplying (the scalar on)
the resulting root edge with a normalization constant computed using squared norm.
To sample from a quantum state in the computational basis, simply repeat the measurement prov with k = idx(v), throw a p-biased coin to determine xk , use followx (v )cedure for edgekto go to level k − 1 and repeat the process.
3.3.2

Gates with simple LIMDD algorithms

As a warm up, before we give the algorithm for arbitrary gates and Clifford gates, we first givealgorithms for several gates that have a relatively simple and efficient LIMDD manipulation operation. In the case of a controlled gate, we distinguish two cases, depending whether the control or
the target qubit comes first; we call these a downward and an upward controlled gate, respectively.
Here, we let Lk denote the unitary applying local gate L on qubit k, i.e., Lk ≜ I⊗n−k ⊗ L ⊗ I⊗k−1 .
Applying a single-qubit Pauli gate Q to qubit k of a LIMDD, by updating the diagram’s rootedge from A to Qk A, i.e., change A = λPn ⊗ · · · ⊗ P1 to λPn ⊗ · · · ⊗ Pk+1 ⊗ QPk ⊗ Pk−1 ⊗ · · · ⊗ P1 .
Since only nodes —and not root edges— need be canonical, this can be done in constant time,
provided that the LIMDD is stored in the natural way (uncompressed with objects and pointers).
Applying any diagonal or antidiagonal single-qubit gate to the top qubit can be doneBrootv , we can construct
efficiently, e.g., applying the T -gate to the top qubit. For root edge e =ex = followx(e), which propagates the root edge’s LIM to the root’s two children. Then, for a
diagonal node α0 β0 , we construct a new root node MakeEdge(α ·e0 , β ·e1 ). For the anti-diagonal






gate α0 β0 , it is sufficient to note that α0 β0 = X · α0 β0 ; thus, we can first apply a diagonal gate,
and then an X gate, as described above.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

16

BrootBroot Sk Broot S †v
Applying a phase gate (S = [ 10 0i ]) to qubit with index k on
kis also efficient. Alg. 6 gives a recursive procedure. If k < n = idx(v)′⇝ vv(top qubit), then note Sk Broot |v⟩ = (Sk Broot Sk† )Sk |v⟩ where Sk Broot Sk† isAB
A iB
the new (O(n)-computable) root ⟨Pauli⟩-LIM because Sk is a Clifford gate.
Hence, we can ‘push’ Sk through the LIMs down the recursion, rebuilding
v0 v1
v0 v1
the LIMDD in the backtrack with MakeEdge on Line 6 and 7. To apply
Sk to v when k = n = idx(v), we finally multiply the high edge label with i on Line 4. Dynamic
programming, using table SGateCache, ensures a linear amount of recursive calls in the numberof nodes m. The total runtime is therefore O(mn3 ), as MakeEdge’s is cubic (see Sec. 4).

Algorithm 6 Apply gate S to qubit k for Pauli-LIMDD
1: procedure SGate(Edge
2:
3:
4:
5:
6:
7:

A

v

if v ∈/ SGateCache thenAv . We let n = idx(v).

with A ∈ Pauli-LIM, k ∈ {1, . . . , idx(v)})▷ Compute result once for v and store in cache:

if idx(v) = k thenSGateCache[v] := MakeEdge(lowv , i · highv )else
SGateCache[v] := MakeEdge(SGate(lowv , k), SGate(highv , k))return Sk ASk† · SGateCache[v]▷ Retrieve result from cache

Applying a Downward Controlled-Pauli gate CQct , where Q is a
idx(v) = c :single-qubit Pauli gate, c the control qubit and t the target qubit with t < c,⇝ v′
v
to a node v can also be done recursively. If idx(v) > c, then since CQct is
A B
A Qt Ba Clifford gate, we may push it through the node’s root label, and apply
it to the children low(v) and high(v), similar to the S gate. Otherwise, if
v0 v1
v0 v1
idx(v) = c, then update v’s high edge label as B 7→ Qt B, and do not recurse.
Alg. 7 shows the recursive procedure, which is similar to Alg. 6 and also has O(mn3 ) runtime.
Av .
Algorithm 7 Apply gate CX with control qubit c and target qubit t for Pauli-LIMDD
We let n = idx(v). We can replace CX, with CY, CZ. modifying Line 4 accordingly (i.e. to Yt , Zt ).

1: procedure CPauliGate(Edge
2:
3:
4:
5:
6:
7:

if v ∈/ CPauliCache thenAvwith A ∈ Pauli-LIM, c, t with 1 ≤ c < t ≤ n)
▷ Compute result once for v and store in cache:

if idx(v) = k thenCPauliCache[v] := MakeEdge(lowv , Xt · highv )else
CPauliCache[v] := MakeEdge(CPauliGate(lowv , c, t), CPauliGate(highv , c, t))return CXtc · A · CXtc † · CPauliCache[v]▷ Retrieve result from cacheSec. 3.3.4 shows that all Clifford gates (including Hadamard and upward CNOT) have runtimeO(n4 ) when applied to a stabilizer state represented as a LIMDD. We first show how to apply
general gates, in Sec. 3.3.3, as this yields some machinery required for Hadamards (specifically, apointwise addition operation).
3.3.3

Applying a generic multi-qubit gate to a state

We use a standard approach [24] to represent quantum gates (2n ×2n unitary matrices) as LIMDDs.
Here a matrix U is interpreted as a function u(r1 , c1 , . . . , rn , cn ) ≜ ⟨r| U |c⟩ on 2n variables, which
returns the entry of U on row r and column c. The function u is then represented using a LIMDD
of 2n levels. The bits of r and c are interleaved to facilitate recursive descent on the structure. In
particular, for x, y ∈ {0, 1}, the subfunction uxy represents a quadrant of the matrix, namely the

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

17

submatrix uxy (r2 , c2 , . . . , rn , cn ) ≜ u(x, y, r2 , c2 , . . . , rn , cn ), as follows:
u0∗

u=

z

u00
u10

}| {
u01
u∗1
u11

(4)

Def. 4 formalizes this idea. Fig. 7 shows a few examples of gates represented as LIMDDs.
Definition 4 (LIMDDs for gates). A LIMDD edge e =Au

can represent a (unitary) 2n × 2n maAu )trix U iff idx(u) = 2n. The value of the matrix cell Ur,c is defined as followr1 c1 r2 c2 ...rn cn (where r, c are the row and column index, respectively, with binary representation r1 , . . . , rn andc1 , . . . , cn . The semantics of a LIMDD edge e as a matrix is denoted [e] ≜ U (as opposed to its
semantics |e⟩ as a vector).

The procedure ApplyGate (Alg. 8) applies a gate U to a state |φ⟩, represented by LIMDDs eU
and eφ . It outputs a LIMDD edge representing U |φ⟩. It works similar to well-known matrix-vector
product algorithms for decision diagrams [24, 27], except that we also handle edge weights with
LIMs (see Fig. 8 for an illustration). Using the followx (e) procedure, we write |φ⟩ and U as|φ⟩ = |0⟩ |φ0 ⟩ + |1⟩ |φ1 ⟩(5)

U = |0⟩ ⟨0| ⊗ U00 + |0⟩ ⟨1| ⊗ U01 + |1⟩ ⟨0| ⊗ U10 + |1⟩ ⟨1| ⊗ U11(6)

Then, on Line 6, we compute each of the fourterms Urc |φc ⟩ for row/column bits r, c ∈ {0, 1}.
We do this by constructing four LIMDDs fr,crepresenting the states |fr,c ⟩ = Ur,c |φc ⟩, using four recursive calls to the ApplyGate algorithm. Next, on Line 7 and 8, the appropriate states are added, using Add (Alg. 9), producing LIMDDs e0 and e1 for the states |e0 ⟩ =U00 |φ0 ⟩ + U01 |φ1 ⟩ and for |e1 ⟩ = U10 |φ0 ⟩ +U11 |φ1 ⟩. The base case of ApplyGate is thecase where n = 0, which means U and |v⟩ aresimply scalars, in which case both eU and eφ areedges that point to the leaf.

−X

X

X ⊗I⊗X

Z

0
0

0
X

1

1

1

I gateZ gateH gate0

CNOT gate:

1

Figure 7: LIMDDs representing various gates.

Algorithm 8 Applies the gate [eU ] to the state |eφ ⟩. Here eU and eφ are LIMDD edges. The
output is a LIMDD edge ψ satisfying |ψ⟩ = [eU ] |eφ ⟩.
1: procedure ApplyGate(Edge eU =
λ·γ

λP

2:

if idx(v) = 0 then return3:

P , Q := RootLabel(

4:

if (P , u, Q , v) ∈/ Apply-cache then

5:
6:
7:
8:
9:
10:
11:

′

′

′

P

u , Edge eφ =γQvwith idx(u) = 2 · idx(v))

1u ), RootLabel(▷ P =Q=1
Qv )▷ Get canonical root labels′▷ Compute result for the first time:

for r, c ∈ {0, 1} do
Edge fr,c := ApplyGate(followrc (P′u ), followc (

Q′

v ))

Edge e0 := Add(f0,0 , f0,1 )Edge e1 := Add(f1,0 , f1,1 )ApplyCache[(P ′ , u, Q′ , v)] := MakeEdge(e0 , e1 )e′ψ := Apply-cache[(P ′ , u, Q′ , v)]
return λγ · e′ψ▷ Store in cache▷ Retrieve from cache

Caching in ApplyGate. A straightforward way to implement dynamic programming would be
to simply store all results of ApplyGate in the cache, i.e., when ApplyGate(
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

λP

u ,

γQ

v ) is

18

eU

eφ

u

λA

u00 u01 u10 u11

v0
×

U

e′ψ

v

γB

e0

e1

f0,0 + f0,1

f1,0 + f1,1

v1
|φ⟩

=

U |φ⟩ = λγ · |e′ψ ⟩ = |eψ ⟩Figure 8: An illustration of ApplyGate (Alg. 8), where matrix U is applied to state B |v⟩, both represented
as Pauli-LIMDDs. The edges f0,0 , f0,1 , etc. are the edges made on Line 6. The dotted box indicates thatthese states are added, using Add, producing edges e0 , e1 , which are then passed to MakeEdge, producing
the result edge. For readability, not all edge labels are shown.

called, store an entry with key (P, u, Q, v) in the cache. This would allow us to retrieve the result
the next time ApplyGate is called with the same parameters. However, we can do much better,in such a way that we can retrieve the result from the cache also when the procedure is called with
B

A

λP

γQ

A

B

y ) satisfying [y ⟩.
u ] = [x ] and |v ⟩ = |x ,parameters ApplyGate(
This can happen even when λP ̸= A or γQ ̸= B; therefore this may prevent many recursive calls.

To this end, we store not just an edge-edge tuple from the procedure’s parameters, but a canonical edge-edge tuple. To obtain canonical edge labels, our algorithms use the function RootLabelABv ) = RootLabel(v )which returns a canonically chosen LIM, i.e., it holds that RootLabel(
whenever A |v⟩ = B |v⟩. A specific choice for RootLabel is the lexicographic minimum of all possible root labels. In Alg. 17, we give an O(n3 )-time algorithm for computing the lexicographicallyminimal root label, following the same strategy as the MakeEdge procedure in Sec. 4.2. As alast optimization, we opt to not store the scalars λ, γ in the cache (they are “factored out”), sothat we can retrieve this result also when ApplyGate is called with inputs that are equal up to
a complex phase. These scalars are then factored back in on Line 11 and 9.
The subroutine Add (Alg. 9) adds two quantum states, i.e., given two LIMDDs representing |e⟩and |f ⟩, it returns a LIMDD representing |e⟩ + |f ⟩. It proceeds by simple recursive descent on thechildren of e and f . The base case is when both edges point to the diagram’s leaf. In this case,these edges are labeled with scalars A, B ∈ C, so we return the edge

A+B

1 .

Algorithm 9 Given two n-LIMDD edges e, f , constructs a new LIMDD edge a with |a⟩ = |e⟩+|f ⟩.
1: procedure Add(Edge e =
2:
3:

A

if idx(v) = 0 then returnif v ̸≼ w then return Add(
A−1 B

v , Edge f =

A+B
B

4:

C := RootLabel(5:if (v, C, w) ∈/ Add-Cache then

6:
7:
8:
9:

B

w

with idx(v) = idx(w))

1w,▷ A, B ∈ C
Av )▷ Normalize for cache lookupw)▷ Compute result for the first time:C

Edge a0 := Add(follow0 ( v ), follow0 (C
Edge a1 := Add(follow1 ( v ), follow1 (Add-Cache[(v, C, w)] := MakeEdge(a0 , a1 )

w ))
w ))

return A · Add-Cache[(v, C, w)]▷ Store in cache▷ Retrieve from cache

Caching in Add. A straightforward way to implement the cache would be to store a tuple with key
BAv ,w ). However, we can do much better; namely, we remark(A, v, B, w) in the call Add(
that we are looking to construct the state A |v⟩+B |w⟩, and that this is equal to A·(|v⟩+A−1 B |w⟩).
This gives us the opportunity to “factor out” the LIM A, and only store the tuple (v, A−1 B, w).
A−1 Bw ) (on Line 4)
We can do even better by finding a canonically chosen LIM C = RootLabel(and storing (v, C, w) (on line Line 8). This way, we get a cache hit at Line 5 upon the call

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

19

D

E

v ,
w ) whenever A−1 B |w⟩ = D −1 E |w⟩. This happens of course in particular when
Add(
(A, v, B, w) = (D, v, E, w), but can happen in exponentially more cases; therefore, this technique
works at least as well as the “straightforward” way outlined above. Finally, on Line 3, we takeadvantage of the fact that addition is commutative; therefore it allows us to pick a preferred orderin which we store the nodes, thus improving possible cache hits by a factor two. We also use C in
the recursive call at Line 6 and 7.

The worst-case runtime of Add is O(n3 2n ) (exponential as expected), where n is the number of qubits. Thiscan happen when the resulting LIMDD is exponential
in the input sizes (bounded by 2n ), as identified for
QMDDs in [52, Table 2]. The reason for this is that addition may remove any common factors, as illustrated
in Fig. 9. However, the Add algorithm is polynomialtime when v = w and v is a stabilizer state, which is
sufficient to show that the Hadamard gate can be efficiently applied to stabilizers represented as LIMDD, as
we demonstrate next in Sec. 3.3.4.
3.3.4

1

4

1

0

1
1

=
1

2

1

2

+
2

3 1
1

1

4
1

Figure 9: Adding two states (0, 1, 0, 4) and(1, 2, 2, 4) as QMDDs can cause an exponentially larger result QMDD (1, 3, 2, 8) due tothe loss of common factors.

LIMDD operations for Clifford gates are polynomial time on stabilizer states

We give an algorithm for the Hadamard gate and then show that it can be applied to a stabilizer
state in polynomial time. Together with the results of Sec. 3.3.2, this shows that all Clifford gatescan be applied to stabilizer states in polynomial time. The key ingredient is Th. 7, which describes
situations in which the Add procedure looks up the same tuples in the cache in both its recursive
calls (modulo ±1). Th. 5 gives the final result.
Theorem 5. Any Clifford gate (H, S, CNOT) can be applied in O(n4 ) time to any (combinationof) qubits to a LIMDD representing a stabilizer state.
Proof. Let |ψ⟩ be an n qubit stabilizer state, represented by a LIMDD with root edgeTh. 1, this LIMDD is a ⟨Pauli⟩-Tower-LIMDD with m = n nodes apart from the leaf.

Av . BySec. 3.3.2 shows that any S-gate can be applied in time O(n3 m), so we get O(n4 ).
Th. 6 shows that any Hadamard gate can be applied on any qubit in time O(n4 ).
Sec. 3.3.2 shows that any downward CNOT-gate can be applied in time O(n3 m), so in this caseO(n4 ). By applying Hadamard to the target and control qubits, before and after the downward
CNOT, we obtain an upward CNOT, i.e., CXct = (H ⊗ H)CXtc (H ⊗ H), still in time O(n4 ).
 1 
√
1/ 2
To apply a Hadamard gate (H = √12 11 −1) to the first qubit,⇝vv′
we first construct edges representing the states |a0 ⟩ = |e0 ⟩ + |e1 ⟩ and
a0
a1
e1 e0
|a1 ⟩ = |e0 ⟩ − |e1 ⟩, using the Add procedure (Alg. 9 and multiplying
e0 + e1e0 − e1
the root edge with −1). Then we construct an edge representing thestate |0⟩ |a0 ⟩ + |1⟩ |a1 ⟩ using MakeEdge. Lastly, the complex factor
on the new edge’s root label is multiplied by √12 . Since the Hadamard
is also a Clifford gate, we can apply this operation to any qubit in the LIMDD by pushing it throughthe LIMs, as we saw in Sec. 3.3.2. Alg. 10 shows the complete algorithm.
Theorem 6. Let e be an n-qubit ⟨Pauli⟩-Tower-LIMDD. HGate(e, k) of Alg. 10 takes O(n4 )time.
Proof. By virtue of the cache, HGate is called at most once per node. Since the LIMDD is a
Tower, there are only n nodes; so HGate is called at most n times. For the node at level k,
HGate makes two calls to Add on Line 4. By applying induction over the qubits 1, . . . , k, using
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

20

Algorithm 10 Apply gate H to qubit k for Pauli-LIMDD
1: procedure HGate(Edge2:v

if v ∈/ HGateCache thenv . We let n = idx(v).

with A ∈ Pauli-LIM, k ∈ {1, . . . , idx(v)})▷ Compute result once for v and store in cache:

if idx(v) = k then√
HGateCache[v] := 1/ 2 · MakeEdge(Add(low(v), high(v)), Add(low(v), −high(v)))
else
HGateCache[v] := MakeEdge(HGate(low(v), k), HGate(high(v), k))

3:
4:
5:
6:
7:

A

A

return Hk AHk† · HGateCache[v]▷ Retrieve result from cacheTh. 7, it is easy to see that at each level, the cache in Alg. 9 is consulted at Line 5 with a tuple(vk , P (k) , vk ) or (vk , −P (k) , vk ). This tells us that Add performs at most 2k recursive calls. Each
recursive call to Add may invoke the MakeEdge procedure, which runs in time O(n3 ), yielding
a total worst-case running time of O(n4 ), when k = Ω(n).
Theorem 7. If Alg. 9 is called on two edges pointing to the same ⟨Pauli⟩-Tower-LIMDD nodev with low(v) = high(v), then the recursive Add calls at Line 6, 7 both lookup the same LIM incache up to a factor ±1.
Qw
be a node on which the algorithm was
Proof. Assume the algorithm is at Line 6. Let w v
called. Let C = Pn ⊗ P be the n qubit Pauli-LIM computed at Line 4 with Pn ∈ Pauli and P an
n−1 qubit Pauli-LIM. At Line 6 and 7, Add makes two recursive calls computing ax for x ∈ {0, 1},
as listed in the header of the below table. The followx (e) semantics yield four cases cases forthe parameters in a recursive Add calls, depending on x and Pn . The following table shows thetuples computed for cache normalization at Line 4 in the recursive call, ignoring the RootLabel()function for now. E.g., if ⇝ denotes cache normalization, thensince QR = ±RQ for P, Q ∈ PauliLIM:
a0 := Add(follow0 (Pn = I:Pn = X:Pn = Y :Pn = Z:

w,
w,
w,
w,

P

w

PQ

w

−iP Q
P

w

w

v ), follow0 (

C

v ))

γR

v ,

γQ−1 R

a1 := Add(follow1 (⇝ (w, P, w)

Q⇝ (w, P Q, w)

Q⇝ (w, −iP Q, w)

Q⇝ (w, P, w)

Q

w,

PQ

w,

P

w,

iP

w,

w

⇝ (v, ±Q, w)v ), follow1 (C⇝ (w, ±P, w)ww⇝ (w, ±P Q, w)w⇝ (w, ±iP Q, w)

−P Q

v ))

w

⇝ (w, ±P , w)

In all four cases, the cache-normalized LIMs computed in both recursive calls are equivalent up
to a factor ±1. Finally, our RootLabel() function from Sec. 4.2.1, which selects the lexicographic
smallest label, satisfies RootLabel(completes the proof.

Qw ) = −RootLabel(−Qw ) for any PauliLIM Q.

This

Since stabilizer states are closed under Clifford gates, ⟨Pauli⟩-Tower-LIMDDs should also be closed
under the respective LIMDD manipulation operations. We show this in Lemma 4 (App. C).

3.4 Comparing LIMDD-based simulation with other methods
Prop. 1 shows exponential advantages of (Pauli-)LIMDDs over three state-of-the-art classical quantum circuit simulators: those based on QMDDs and MPS [36, 37], and the Clifford + T simulator.
In this section we prove the proposition, mainly using results from the current section: To show
the separation between simulation with LIMDDs and Clifford + T , we present Th. 8.
Our proofs often rely on the fact that LIMDDs are exponentially more succinct representations of acertain class of quantum states S that are generated by circuits with a certain (non-universal) gate
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

21

set G. For instance, the stabilizer states that are generated by the Clifford gate set. LIMDD-based
simulation —similar to MPS [38] and QMDD-based [28] simulation— proceeds by representing astate |φt ⟩ at time step t as a LIMDD φt . It then applies the gate Ut ∈ G in the circuit corresponding
to this time step to obtain a LIMDD φt+1 with |φt+1 ⟩ = Ut |φt ⟩, thus yielding strong simulation at
the final time step as reading amplitudes from the final LIMDD is easy (see Sec. 3.1).
It follows that LIMDD-based simulation is efficient provided that it can execute all gates Ut inpolynomial time (in the size of the LIMDD representation), at least for the states in S. Note inparticular that since the execution stays in S, i.e., |φt ⟩ ∈ S =⇒ |φt+1 ⟩ ∈ S, the representation size
can not grow to exponential size in multiple steps (S can be considered an inductive invariant in
the style of Floyd [53] and de Bakker & Meertens [54]). On the other hand, since MPS and QMDD
are exponentially sized for cluster states, they necessarily require exponential time on circuits
computing this family of states.
3.4.1

LIMDD is exponentially faster than QMDD-based simulation

As state set S, we select the stabilizer states and for G the Clifford gates. Th. 1 shows that LIMDDs
for stabilizers are always quadratic in size in the number of qubits n, as the diagram contains nnodes and n + 1 LIMs, each of size at most n (see Def. 3). Sec. 3.3.2 shows that LIMDD can execute
all Clifford gates on stabilizer states in time O(n4 ).
On the other hand, Th. 2 shows that QMDDs for cluster states are exponentially sized. It follows
that in simulation also, there is an exponential separation between QMDD and LIMDD, proving) (Prop. 1 Item 3).
that QSimQMDD= Ω∗ (2n · QSimLIMDDCCFor the other direction, we now show that LIMDDs are at most a factor O(n3 ) slower than QMDDs
on any given circuit. First, a LIMDD never contains more nodes than a QMDD representing thesame state (because QMDD is by definition a specialization of LIMDD, see Sec. 3.1). The LIMDD
additionally uses O(n) memory per node to store two Pauli LIMs; thus, the total memory usage is
at most a factor O(n) worse than QMDDs for any given state. The ApplyGate and Add algorithms
introduced in Sec. 3.3.3 are very similar to the ones used for QMDDs in [1, 24]. In particular, ourApplyGate and Add algorithms never make more recursive calls than those for QMDDs. However,
one difference is that our MakeEdge algorithm runs in time O(n3 ) instead of O(1). Therefore,
in the worst case these LIMDD algorithms make the same number of recursive calls to ApplyGate
and Add, in which case they are slower by a factor O(n3 ).
Finally, Corollary 3 shows that the pseudo-cluster state |φ⟩ has a polynomial representation in
LIMDD. By definition of the pseudo-cluster state, post-selecting (constraining) the top qubit to
0 (or 1) yields the cluster state |Gn ⟩. Therefore, QMDD for the pseudo-cluster state must haveexponential size, as constraining can never increase the size of DD [55, Th 2.4.1]. Together withthe universal simulation discussed above, this proves that the above also holds for for a simulator
based on the combination QMDD ∪ Stab (Prop. 1 Item 5).
3.4.2

LIMDD is exponentially faster than MPS

In Sec. 3.4.1, we saw that LIMDD can simulate the cluster state in polynomial time. On the other
hand, Th. 4 shows that MPS for cluster states are exponentially sized. It follows that in simulationalso, there is an exponential separation between MPS and LIMDD, proving Prop. 1 Item 2.
3.4.3

LIMDD is exponentially faster than Clifford + T

In this section, we consider a circuit family that LIMDDs can efficiently simulate, but which isdifficult for the Clifford+T simulator because the circuit contains many T gates, assuming theExponential Time Hypothesis (ETH, a standard complexity-theoretic assumption which is widely

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

22

believed to be true). This method decomposesa given quantum circuit into a circuit consisting0
only of Clifford gates and the T = 10 eiπ/4gate, as explained in Sec. 2.
⊗nThe circuit family, given my McClung [56], maps the input state |0⟩to the n-qubit W state|Wn ⟩, which is the equal superposition over computational-basis states with Hamming weight 1,
1|Wn ⟩ = √ (|100 . . . 00⟩ + |010 . . . 00⟩ + · · · + |000 . . . 01⟩)n
Arunachalam et al. showed that, assuming ETH, any circuit which deterministically produces the|Wn ⟩ state in this way requires Ω(n) T gates [57]. Consequently, the Clifford + T simulator cannot
efficiently simulate the circuit family, even when one allows for preprocessing with a compilation
algorithm aiming to reduce the T -count of the circuit (such as the ones developed in [58, 59]).
Th. 8 now shows that the exponential separation between simulation with LIMDD and Clifford + T ,+Ti.e., that QSimClifford= Ω(2n · QSimLIMDD) (Prop. 1 Item 1). App. F gives its proof.
CC⊗nTheorem 8. There exists a circuit family Cn such that Cn |0⟩ = |Wn ⟩, that Pauli-LIMDDs canefficiently simulate. Here simulation means that it constructs representations of all intermediate
states, in a way which allows one to, e.g., efficiently simulate any single-qubit computational-basis
measurement or compute any computational basis amplitude on any intermediate state and theoutput state.
We note that we could have obtained a similar result using the simpler scenario where one applies a
T gate to each qubit of the (|0⟩+|1⟩)⊗n input state. However, our goal is to show that LIMDDs cannatively simulate scenarios which are relevant to quantum applications, such as the stabilizer statesfrom the previous section. The W state is a relevant example, as several quantum communication
protocols use the W state [60–62]. In contrast, the circuit with only T gates yields a product state,hence it is not relevant unless we consider it as part of a larger circuit which includes multi-qubit
operations.
Lastly, it would be interesting to analytically compare LIMDD with general stabilizer rank based
simulation (without assuming ETH). However, this would require finding a family of states with
provably superpolynomial stabilizer rank, which is a major open problem. Instead, we implemented
a heuristic algorithm by Bravyi et al. [14] to empirically find upper bounds on the stabilizer rank
and applied it to a superset of the W states, so-called Dicke states, which can be represented aspolynomial-size LIMDD. The O(n2 )-size LIMDD can be obtained via a construction by Bryant [19],since the amplitude function of a Dicke state is a symmetric function. The results hint at a possible
separation but are inconclusive due to the small number of qubits which the algorithm can feasibly
investigate in practice. See App. G for details.

4 Canonicity: Reduced LIMDDs with efficient MakeEdge algorithm
Unique representation, or canonicity, is a crucial property for the efficiency and effectiveness ofdecision diagrams. In the first place, it allows for circuit analysis and simplification [20, 27], byfacilitating efficient manipulation operations through dynamic programming efficiently, as discussed
in Sec. 3.3. In the second place, a reduced diagram is smaller than an unreduced diagram because
it merges nodes with the same semantics. For instance, Pauli-LIMDDs allow all states in the same
≃Pauli equivalence class to be merged. Here, we define a reduced Pauli-LIMDD, which is canonical.
In general, many different LIMDDs can represent a given quantum state, as illustrated in Fig. 10.
However, by imposing a small number of constraints on the diagram, listed in Def. 5 and visualized
in Fig. 11, we ensure that every quantum state is represented by a unique ‘reduced ’ Pauli-LIMDD.
We present a MakeEdge algorithm (Alg. 11 in Sec. 4.2) that computes a canonical node assuming
its children are already canonical. The algorithms for quantum circuit simulation in Sec. 3.3 ensurethat all intermediate LIMDDs are reduced by creating nodes exclusively through this subroutine.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

23

1
√
I⊗I
2

1
√
X⊗I
2

1
√
I⊗I
2

1
√
Z⊗I
2

−X

X
v′

v
0

0
1

v′

v
0

v

0
1

v
0

0

1

1

Figure 10: Four different Pauli-LIMDDs representing the Bell state √12 (|00⟩ + |11⟩). From left to right: as
I-LIMDD, swapping high and low nodes v, v ′ by placing an X on the root LIM, merging v ′ into v by observing
that |v⟩ = X |v ′ ⟩ and selecting a different high LIM −X together with changing the root LIM. This section
shows that selecting a unique high LIM is the most challenging, as in general many LIMs can be chosen.

4.1 LIMDD canonical form
The main insight used to obtain canonical decision diagrams is that a canonical form can becomputed locally for each node, assuming its children are already canonical. In other words,
if the diagram is constructed bottom up, starting from the leaf, it can immediately be madecanonical. (This is why decision diagram manipulation algorithms always construct the diagramin the backtrack of the recursion using a typical ‘MakeNode’ procedure for constructing canonical
α

β

1

β/α

w
with α, β ∈ C \ {0} cannodes [24], like in Sec. 3.3.) For instance, a QMDD node v
be reduced into a canonical node by dividing out a common factor α and placing it on the rootw
edge. Assuming that v, w are canonical, the resulting node v
can be stored as a tuple
β(1, v, /α, w) in a hash table. Moreover, any other node that is equal to this node up to a scalar is
reduced to the same tuple with this strategy [27] and thus merged in the hash table.

For LIMDD, we use a similar approach of dividing out ‘common LIM factors.’ However, we need to
do additional work to obtain a unique high edge label (β/α in the example above), as the PauliLIM
group is more complicated than the group of complex numbers (scalars).
Def. 5 gives reduction rules for LIMDDs and Fig. 11 illustrates them. The merge (1) and low
factoring (4) rules fulfill the same purpose as in the QMDD case discussed above. In a PauliLIMDD, we may always swap high and low edges of a node v by multiplying the root edge LIM
with X ⊗ I, as illustrated in Fig. 10. The low precedence rule (3) makes this choice deterministic,
but only in case low(v) ̸= high(v). Next, the zero edges (2) rule handles the case when α ór β
are zero in the above, as in principle a edge e with label 0 could point to any node on the nextlevel k, as this always yields a 0 vector of length 2k (see semantics below Def. 2). The rule forces
low(v) = high(v) in case either edge has a zero label. We explain the interaction among the zeroedges (2), low precedence (3) and low factoring (4) rules below. Finally, the high determinism rule(5) defines a deterministic function to choose LIMs on high edges, solving the most challenging
problem of uniquely selecting a LIM on the high edge. We give an O(n3 ) algorithm for this function
in Sec. 4.2.
Definition 5 (Reduced LIMDD). A Pauli-LIMDD is reduced when it satisfies the following constraints. It is semi-reduced if it satisfies all constraints except possibly high determinism.
1. Merge: No two nodes are identical: We say two nodes v, w are identical if low(v) = low(w),high(v) = high(w), label(low(v)) = label(low(w)), label(high(v)) = label(high(w)).
2. (Zero) edge: For any edge (v, w) ∈ high ∪ low, if label(v, w) = 0, then both edges outgoingfrom v point to the same node, i.e., high(v) = low(v) = w.
3. Low precedence: Each node v has low(v) ≼ high(v), where ≼ is a total order on nodes.
4. Low factoring: The label on every low edge to a node v is the identity I⊗idx(v) .
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

24

C

D

u

v

C
⇝

A

B
A
w

C

w

(1) Merge u, v into v(I ⊗ A) · C

C

0

B

A

A

B

A

B

u

w

w

u

u

w

u

w

v
B

(X ⊗ Ik ) · C

C

⇝ v′
B 0 B

v

A
B

C

D

⇝ v′

v

w

(2) Zero edges(3) Low precedencewith u ≼ w

v′

⇝

v

w

(4) Low factoring

v′

⇝

v
A

A−1 B
u

C · (AH −1 )

C

w

u

H
w

u

(5) High determinism
H = HighLabel(v)Figure 11: Illustrations of the reduction rules from Def. 5 applied at level k + 1 (i.e., k + 1 = idx(v) = idx(v ′ )).
Note that, in general, the top edges are not necessarily root edges, but could be high and low edges for nodes
on level k + 2. So, in general, there can be multiple such incoming edges (dashed and solid).

5. High determinism: The label on the high edge of any node v is Bhigh = HighLabel(v),where HighLabel is a function that takes as input a semi-reduced n-Pauli-LIMDD node v, andoutputs an (n − 1)-Pauli-LIM Bhigh satisfying |v⟩ ≃Pauli |0⟩ |low(v)⟩ + |1⟩ ⊗ Bhigh |high(v)⟩.
Moreover, for any other semi-reduced node w with |v⟩ ≃Pauli |w⟩, it satisfies HighLabel(w) =Bhigh . In other words, the function HighLabel is constant within an isomorphism class.
We make several observations about reduced LIMDDs. First, let us apply this definition to a state|0⟩ ⊗ A |φ⟩ + |1⟩ ⊗ B |ψ⟩ with |φ⟩ ̸≃Pauli |ψ⟩, where A, B ∈ PauliLIM. Assume we already havecanonical LIMDDs for φ and ψ (note that necessarily φ ̸= ψ). We will transform this node sothat it satisfies all the reduction rules above. There is a choice between representing this state
A

B

B

A

φ
ψ
as either φor ψ, as these are related by the isomorphism X ⊗ I. The
low precedence rule resolves this choice here. Assuming φ ≺ ψ, low factoring can now be realized

by dividing out the LIM A, yielding a node

I

φ
I

A−1 B
B

−1

ψ

(with root edge I ⊗ A as in Fig. 11Aφ(4)). Otherwise, if ψ ≺ φ, we obtain node ψwith incoming edge X ⊗ B. Finally,−1since there might be other LIMs Bhigh not equal to B A that yield the same state, the high
I

Bhigh

φ
as shown in Fig. 12.
determinism rule is finally needed to obtain a canonical node ψThis last step turns a semi-reduced node into a (fully) reduced node. Sec. 4.2 discusses it in detail.

Now, let us apply the definition to a state |1⟩ ⊗ A |φ⟩. First, notice that the zero edges rule forces
low(v) = high(v) = φ in this case. There is a choice between representing this state as either
φ

A

0

φ

0

A

φ
or φ
, which denote the states |0⟩ ⊗ A |φ⟩ and |1⟩ ⊗ A |φ⟩, as these are
related by the isomorphism X ⊗ I. The low factoring rule requires that the low edge label isA0φI, yielding a node of the form φwith root label X ⊗ A: In other words, this ruleenforces swapping high and low edges, placing a X on the root label, and dividing out the LIM
A. Consequently, the high edge must be labeled with 0, and therefore, semi-reduction, in thiscase, coincides with (full) reduction (no high determinism is required). Notice also that there is
no reduced LIMDD for the 0-vector, because low factoring requires low edges with label I. This is
not a problem, since the 0-vector is not a quantum state.

The rules in Def. 5 are defined only for Pauli-LIMDDs, to which our results pertain (except forthe brief mention of ⟨X⟩ and ⟨Z⟩-LIMDDs in Sec. 3.2). We briefly discuss alternative groups here.
If G is a group without the element X ̸∈ G, the reduced G-LIMDD based on the same rules is notuniversal (does not represent all quantum states), because the low precedence rule cannot always
be satisfied, since it requires that v0 ≼ v1 for every node. Hence, in this case, reduced G-LIMDD
cannot represent a state |0⟩ |v0 ⟩ + |1⟩ |v1 ⟩ when v1 ≺ v0 . However, it is not difficult to formulate
rules to support these groups G; for instance, when G = {I}, we recover the QMDD and may use
its reduction rules [28].
Nodes and edges in a reduced LIMDD need not represent normalized quantum states, just like in(unreduced) LIMDDs as explained in Sec. 3.1. Consider, e.g., node ℓ2 in Fig. 3, which represents
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

25

⇝

v′

0

⇝

v ′′

v

A−1 B
vL

w

⇝

v′

B
A

0

A

(I ⊗ A)Broot

I⊗A

I⊗A

≼

v
Bhigh

vR

Figure 12: Reduced node construction in case |φ1 ⟩ = 0 (left), and |φ0 ⟩ , |φ1 ⟩ ̸= 0 and vL ≼ vR (right). Not
shown: for cases |φ0 ⟩ = 0 and vR ≼ vL , we take instead root edge X ⊗ A and swap low/high edges.

state [1, 1, i, i]⊤ . Because the normalization constant was divided out (see factor 1/4 on the rootedge), this state is not normalized. In fact, the root node does not need to be normalized, as evenreduced LIMDDs can represent any vector (except for the zero vector).
Lastly, the literature on other decision diagrams [18, 19, 48] often considers a “redundant test”
or “deletion” rule to remove nodes with the same high and low child. This would introduce theskipping of qubit levels, which our syntactic definition disallows, as already discussed in Footnote ‡.
However, if needed Def. 2 could be adapted and a deletion rule could be added to Def. 5.
We now give a proof of Th. 9, which states that reduced LIMDDs are canonical.
Theorem 9 (Node canonicity). For each n-qubit quantum state |φ⟩, there exists a unique reduced
Pauli-LIMDD L with root node vL such that |vL ⟩ ≃ |φ⟩.
Proof. We use induction on the number of qubits n to show universality (the existence of an
isomorphic LIMDD node) and uniqueness (canonicity).
Base case. If n = 0, then |φ⟩ is a complex number λ. A reduced Pauli-LIMDD for this state is
the leaf node representing the scalar 1. To show it is unique, consider that nodes v other than theleaf have an idx(v) > 0, by the edges rule, and hence represent multi-qubit states. Since the leaf
node itself is defined to be unique, the merge rule is not needed and canonicity follows.
Finally, |φ⟩ is represented by root edge

λ1 .

Inductive case. Suppose n > 0. We first show existence, and then show uniqueness.
Part 1: existence. We use the unique expansion of |φ⟩ as |φ⟩ = |0⟩ ⊗ |φ0 ⟩ + |1⟩ ⊗ |φ1 ⟩ where |φ0 ⟩and |φ1 ⟩ are either (n − 1)-qubit state vectors, or the all-zero vector. We distinguish three cases
based on whether |φ0 ⟩ , |φ1 ⟩ = 0.
Case |φ0 ⟩ , |φ1 ⟩ = 0: This case is ruled out because |φ⟩ ≠ 0.
Case |φ0 ⟩ = 0 or |φ1 ⟩ = 0: In case |φ0 ⟩ ̸= 0, by the induction hypothesis, there exists a PauliLIMDD with root node w satisfying |w⟩ ≃ |φ0 ⟩. By definition of ≃, there exists an n-qubit Pauli
isomorphism A such that |φ0 ⟩ = A |w⟩. We construct the following reduced Pauli-LIMDD for |φ⟩:
I

0

I⊗A

v as illustrated in Fig. 12 (left). In case |φ1 ⟩ =, adding a root edge er =
̸ 0,
we do the same for root node In case |φ1 ⟩ ≠ 0, we do the same for root |w⟩ ≃ |φ1 ⟩ = A |w⟩, butX⊗Av (similar to Fig. 11 (3)). In
switch the high and the low edge by instead a root edge er =both cases, it is easy to check that the root node v is reduced as it can be represented by a tuple(I, w, 0, w), where w is canonical because of the induction hypothesis. Also in both cases, we also
have |φ⟩ = |er ⟩ because either |φ⟩ = I ⊗ A |v⟩ or |φ⟩ = X ⊗ A |v⟩.
w

v

w

Case |φ0 ⟩ , |φ1 ⟩ ̸= 0: By applying the induction hypothesis twice, there exist Pauli-LIMDDs L
and R with root nodes |vL ⟩ ≃ |φ0 ⟩ and |vR ⟩ ≃ |φ1 ⟩. The induction hypothesis implies only a ‘local’
reduction of LIMDDs L and R, but not automatically a reduction of their union. For instance, Lmight contain a node v and R a node w such that v ≃ w. While the other reduction rules ensurethat v and w will be structurally the same, the induction hypothesis only applies the merge rule L
and M in isolation, leaving two copies of identical nodes v, w. We can solve this by applying merge

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

26

on the union of nodes in L and M , to merge any equivalent nodes, as they are already structurally
equivalent by the induction hypothesis. This guarantees that (also) vL , vR are identical nodes.
By definition of ≃, there exist n-qubit Pauli isomorphisms A and B such that |φ0 ⟩ = A |vL ⟩and |φ1 ⟩ = B |vR ⟩. In case vL ≼ vR , we construct the following reduced Pauli-LIMDD for |φ⟩:
I

E

I

A−1 B

vR
vR
the root node is vL v
, where E is the LIM computed by HighLabel( vL).
Otherwise, if vR ≼ vL , then we construct the following reduced Pauli-LIMDD for |φ⟩: the root nodeI

FI

B −1 A

vL
vR
v
is vR
, where F = HighLabel( vL). It is straightforward to check that, inboth cases, this Pauli-LIMDD is reduced. Moreover, |v⟩ isomorphic to |φ⟩ as illustrated in Fig. 12(right).

Part 2: uniqueness. To show uniqueness, let L and M be reduced LIMDDs with root nodesvL , vM such that |vL ⟩ ≃ |φ⟩ ≃ |vM ⟩, as follows,
0
vL

AL

vL

BL

0
vM

1
vL

AM

vM

BM

1
vM

(7)

The fact that these nodes are isomorphic means that there is a Pauli isomorphism P such that
P |vL ⟩ = |vM ⟩. We write P = λPtop ⊗ Prest ̸= 0 where Ptop is a single-qubit Pauli matrix and Prestan (n − 1)-qubit Pauli LIM. Expanding the semantics of vL and vM , we obtain,01
01
λPtop ⊗ Prest (|0⟩ ⊗ AL |vL⟩ + |1⟩ ⊗ BL |vL⟩) = |0⟩ ⊗ AM |vM⟩ + |1⟩ ⊗ BM |vM⟩.

(8)

We distinguish two cases from here on: where Ptop ∈ {I, Z} or Ptop ∈ {X, Y }.
Case Ptop = I, Z. If Ptop = [ 10 z0 ] for z ∈ {1, −1}, then Eq. 8 gives:
00λPrest AL |vL⟩ = AM |vM
⟩

and

1
1
zλPrest BL |vL⟩ = BM |vM
⟩

(9)

0
0
0
⟩ is isomorphic⟩. Hence |vL⟩ = |vM
By low factoring, we have AL = AM = I, so we obtain λPrest |vL000with |vM ⟩, so by the induction hypothesis, we have vL = vM . We now show that also vL = vM byconsidering two cases.
11
11
represent isomorand vM⟩, so the nodes vL⟩ = BM |vM
BL ̸= 0 and BM ̸= 0: then zλPrest BL |vL
11
phic states, so by the induction hypothesis we have vL = vM . We already noticed by the low
factoring rule that vL and vM have I as low edge label. By the high edge rule, their high edge
labels are HighLabel(vL ) and HighLabel(vM ), and since the nodes vL and vM are semi-reduced
and |vL ⟩ ≃ |vM ⟩, we have HighLabel(vM ) = HighLabel(vL ) by definition of HighLabel.
1⟩. Since the state vector
BL = 0 or BM = 0: In case BL = 0, we see from Eq. 9 that 0 = BM |vM1
̸ 0 by the observation that a reduced node does not represent the zero vector, it follows|vM ⟩ =
1
⟩ = 0. We have zλ ̸= 0,
that BM = 0. Otherwise, if BM = 0, then Eq. 9 yields zλPrest BL |vL
1
Prest ̸= 0 by definition, and we observed |vL ⟩ ≠ 0 above. Therefore BL = 0. In both cases,
BL = BM .

We conclude that in both cases vL and vM have the same children and the same edge labels, sothey are identical by the merge rule.
 ∗Case Ptop = X, Y . If Ptop = z0 z0 for z ∈ {1, i}, then Eq. 8 gives:
01
λzPrest AL |vL⟩ = BM |vM
⟩

and

1
0
λz ∗ Prest BL |vL⟩ = AM |vM
⟩.

0
1
1
0
By low factoring, AL = AM = I, so we obtain zλPrest |vL⟩ = BM |vM⟩ and λz ∗ Prest BL |vL⟩ = |vM⟩.
To show that vL = vM , we consider two cases.
0110BL ̸= 0 and BM ̸= 0: we find |vL⟩ ≃ |vM⟩ and |vL
⟩ ≃ |vM⟩, so by the induction hypothesis,
0
1
1
0
1
1
0
0
vL = vM and vL = vM . By low precedence, it must be that vL
= vM
= vL
= vM
. Now usehigh determinism to infer that BL = BM as in the Ptop = I, Z case.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

27

BL = 0 or BM = 0: This case leads to a contradiction and thus cannot occur. BL cannot be zero,0
because then |vM⟩ is the all-zero vector, which we excluded. The other case: if BM = 0, then00
it must be that λzPrest AL |vL⟩ is zero. Since λzPrest ̸= 0 and AL = I, it follows that |vL⟩ isthe all-zero vector, which is again excluded.
We conclude that vL and vM have the same children and the same edge labels for all choices ofPtop , so they are identical by the merge rule.

4.2 The MakeEdge subroutine: Maintaining canonicity during simulation
To construct new nodes and edges, our algorithms use the MakeEdge subroutine (Alg. 11),
as discussed in Sec. 4.1. MakeEdge produces a reduced parent node (with root edge) given two
reduced children, so that the LIMDD representation becomes canonical. Here we give the algorithm
for MakeEdge and show that it runs in time O(n3 ) (assuming the input nodes are reduced).
The MakeEdge subroutine distinguishes two cases, depending on whether both children are nonzero vectors, which both largely follow the discussion below Def. 5. It works as follows:• First it ensures low precedence, switching e0 and e1 if necessary at Line 3. This is also done
if e0 ’s label A is 0 to allow for low factoring (avoiding divide by zero).
• Low factoring, i.e., dividing out the LIM A, placing it on the root node, is visualized inFig. 12 for the cases e1 = 0/e1 ̸= 0, and done in the algorithm at Line 6,7 / 9,11.
• The zero edges rule is enforced in the B = 0 branch by taking v1 := v0 .
• The canonical high label Bhigh is computed by GetLabels, discussed below, for the semiIÂ

v1
with v0 ̸= v1 . With the resulting high label, it now satisfies the
reduced node v0 w
high determinism rule of Def. 5 with HighLabel(w) = Bhigh .

• Finally, we merge nodes by creating an entry (v0 , Bhigh , v1 ) in a table called the uniquetable [40] at Line 13.
All steps except for GetLabels have complexity O(1) or O(n) (for checking low precedence, we
use the nodes’ order in the unique table). The algorithm GetLabels, which we sketch below inSec. 4.2.1 and fully detail in App. D, has runtime O(n3 ) if both input nodes are reduced, yielding
an overall complexity O(n3 ).
4.2.1

Choosing a canonical high-edge label

In order to choose the canonical high edge label of node v, the MakeEdge algorithm calls GetLabels (Line 10 of Alg. 11). The function GetLabels returns a uniquely chosen LIM Bhigh among
all possible high-edge labels which yield LIMDDs representing states that are Pauli-isomorphic
to |v⟩. We sketch the algorithm for GetLabels here and provide the algorithm in full detailin App. D (Alg. 12). First, we characterize the eligible high-edge labels. That is, given a semireduced node
I

v0

I

v

Â v
1

, we characterize all C such that the node

v0

I

C v
1

is isomorphic

Â v
v
1

to v0
. Our characterization shows that, modulo some complex factor, the eligible labels
C are of the formC ∝ g0 · Â · g1 , for g0 ∈ Stab(|v0 ⟩), g1 ∈ Stab(|v1 ⟩)(10)
where Stab(|v0 ⟩) and Stab(|v1 ⟩) are the stabilizer subgroups of |v0 ⟩ and |v1 ⟩, i.e., the already
reduced children of our input node v. Note that the set of eligible high-edge labels might be
exponentially large in the number of qubits. Fortunately, eq. (10) shows that this set has apolynomial-size description by storing only the generators of the stabilizer subgroups.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

28

Algorithm 11 Algorithm MakeEdge takes two root edges to (already reduced) nodes v0 , v1 , thechildren of a new node, and returns a reduced node with root edge. It assumes that idx(v0 ) =idx(v1 ) = n. We indicate which lines of code are responsible for which reduction rule in Def. 5.
1: procedure MakeEdge(Edge e0
2:
3:
4:
5:

if v0 ̸≼ v1 or A = 0 then

9:
10:
11:
12:
13:
14:

B v
1 , with v , v
0 1 reduced, A ̸= 0 or B ̸= 0)▷ Enforce low precedence and enable factoring

return (X ⊗ I⊗n ) · MakeEdge(e1 , e0 )if B = 0 thenv1 := v0

▷ Enforce zero edges

I⊗n

6:
7:
8:

A v
0 , e10 v
0v := v0
Broot := I ⊗ A
else
Â := A−1 B
Bhigh , Broot := GetLabels(Â, v0 , v1 )I⊗n▷ Enforce low factoring▷ Enforce high determinism

Bhigh

v1
v :=
Broot := (I ⊗ A)Brootv0▷ Enforce low factoring▷ Broot |v⟩ = |0⟩ ⊗ A |v0 ⟩ + |1⟩ ⊗ B |v1 ⟩▷ Broot |v⟩ = |0⟩ ⊗ |v0 ⟩ + |1⟩ ⊗ A−1 B |v1 ⟩▷ (I ⊗ A)Broot |v⟩ = |0⟩ ⊗ A |v0 ⟩ + |1⟩ ⊗ B |v1 ⟩vr := Find or create unique table entry Unique[v] = (v0 , Bhigh , v1 )returnBroot▷ Enforce mergevr

Our algorithm chooses the lexicographically smallest eligible label, i.e., the smallest C of the formC ∝ g0 Âg1 (the definition of ‘lexicographically smallest’ is given in App. A). To this end, we use
two subroutines: (1) an algorithm which finds (a generating set of) the stabilizer group Stab(|v⟩)
of a LIMDD node v; and (2) an algorithm that uses these stabilizer subgroups of the children nodes
to choose a unique representative of the eligible-high-label set from eq. (10).
For (1), we use an algorithm which recurses on the children nodes. First, we note that, if the Pauli
LIM A stabilizes both children, then I ⊗ A stabilizes the parent node. Therefore, we compute (agenerating set for) the intersection of the children’s stabilizer groups. Second, our method finds
out whether the parent node has stabilizers of the form Pn ⊗ A for Pn ∈ {X, Y, Z}. This requires
us to decide whether certain cosets of the children’s stabilizer groups are empty. These groups are
relatively simple, since, modulo phase, they are isomorphic to a binary vector space, and cosets
are hyperplanes. We can therefore rely in large part on existing algorithms for linear algebra invector spaces. The difficult part lies in dealing with the non-abelian aspects of the Pauli group.
We provide the full algorithm, which is efficient, also in App. D.
Our algorithm for (2) applies a variant of Gauss-Jordan elimination to the generating sets ofStab(|v0 ⟩) and Stab(|v1 ⟩) to choose g0 and g1 in eq. (10) which, when multiplied with Â as ineq. (10), yield the smallest possible high label C. (We recall that Gauss-Jordan elimination,a standard linear-algebra technique, is applicable here because the stabilizer groups are groupisomorphic to binary vector spaces, see also App. A). We explain the full algorithm in App. D.
4.2.2

Checking whether two LIMDDs are Pauli-equivalent

To check whether two states represented as LIMDDs are Pauli-equivalent, it suffices to check
whether they have the same root node. Namely, due to canonicity, and in particular the Merge
rule (in Def. 5), there is a unique LIMDD representing a quantum state up to phase and local Pauli
operators.

5 Related work
We mention related work on classical simulation formalisms and decision diagrams other than QMDD.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

29

The Affine Algebraic Decision Diagram, introduced by Tafertshofer and Pedam [63], and by Sannerand McAllister [26], is akin to a QMDD except that its edges are labeled with a pair of real numbers(a, b)⊗nv represents the state vector a |v⟩ + b |+⟩(a, b), so that an edge
(i.e., here b is added toeach element of the vector a |v⟩). To the best of our knowledge, this diagram has not been applied
to quantum computing.

Context-Free-Language Ordered Binary Decision Diagrams (CFLOBDDs) [64, 65] extend BDDs
with insights from visibly pushdown automata [66]. An extension of CFLOBDD to the complex
domain [67] shows good performance for various simulation of quantum computing benchmarks.
Sentential Decision Diagrams [68] generalize BDDs by replacing their total variable order with avariable tree (vtree). Although Kisa et al. [69] introduced an SDD which represents probability
distributions, SDDs have not yet been used to simulate quantum computing, to the best of ourknowledge. The Variable-Shift SDD (VS-SDD) [70] improves on the SDD by merging isomorphic
vtree nodes. We remark that CFLOBDDs are similar to VS-SDD with a balanced vtree.
Günther and Drechsler introduced a BDD variant [71] which, in LIMDD terminology, has a labelon the root node only. To be precise, this diagram’s root edge is labeled with an invertible matrix
A ∈ Fn×n. If the root node represents the function r : Bn → B, then the diagram represents
2
the function f (⃗x) = r(A · ⃗x). (This concepts extends trivially to the domain of pseudo-Booleanfunctions, by replacing the BDD with an ADD.) In contrast, LIMDDs allow a label on every edgein the diagram, not only the root edge. We show that this is essential to capture stabilizer states.
A multilinear arithmetic formula is a formula over +, × which computes a polynomial in which
no variable appears raised to a higher power than 1. Aaronson showed that some stabilizer statesrequire superpolynomial-size multilinear arithmetic formulas [33, 45].

6 Discussion
We have introduced LIMDD, a novel decision diagram-based method to simulate quantum circuits,
which enables polynomial-size representation of a strict superset of stabilizer states and the states
represented by polynomially large QMDDs. To prove this strict inclusion, we have shown thefirst lower bounds on the size of QMDDs: they require exponential size for certain families ofstabilizer states. Our results show that these states are thus hard for QMDDs. We also give thefirst analytical comparison between simulation based on decision diagrams, and matrix productstates, and the Clifford + T simulator.
LIMDDs achieve a more succinct representation than QMDDs by representing states up to local
invertible maps which uses single-qubit (i.e., local) operations from a group G. We have investigatedthe choices G = Pauli, G = ⟨Z⟩ and G = ⟨X⟩, and found that any choice suffices for an exponential
advantage over QMDDs; notably, the choice G = Pauli allows us to succinctly represent any
stabilizer state. Furthermore, we showed how to simulate arbitrary quantum circuits, encoded asPauli-LIMDDs. The resulting algorithms for simulating quantum circuits are exponentially fasterthan for QMDDs in the best case, and never more than a polynomial factor slower. In the case of
Clifford circuits, the simulation by LIMDDs is in polynomial time (in contrast to QMDDs).
We have shown that Pauli-LIMDDs can efficiently simulate a circuit family outputting the W states,in contrast to the Clifford + T simulator which requires exponential time to do so (assuming thewidely believed ETH), even when allowing for preprocessing of the circuit with a T -count optimizer.
Since we know from experience that implementing a decision diagram framework is a major endeavor, we leave an implementation of the Pauli-LIMDD, in order to observe its runtimes in practiceon relevant quantum circuits, to future work. We emphasize that from the perspective of algorithm
design, we have laid all the groundwork for such an implementation, including the key ingredient for
the efficiency of many operations for existing decision diagrams: the existence of a unique canonicalrepresentative of the represented function, combined with a tractable MakeEdge algorithm to findit.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

30

Regarding extensions of the LIMDD data structure, an obvious next step is to investigate otherchoices of G. Of interest are both the representational capabilities of such diagrams (do they
represent interesting states?), and the algorithmic capabilities (can we still find efficient algorithms
which make use of these diagrams?). In this vein, an important question is what the relationship is
between G-LIMDDs (for various choices of G) and existing formalisms for the classical simulationof quantum circuits, such as those based on match gates [72–74] and tensor networks [29, 75]. It
would also be interesting to compare LIMDDs to graphical calculi such as the ZX calculus [76],following similar work for QMDDs [77].
Lastly, we note that the current definition of LIMDD imposes a strict total order over the qubits
along every path from root to leaf. It is known that the chosen order can greatly influence the sizeof the DD [55, 78], making it interesting to investigate variants of LIMDDs with a flexible ordering,for example taking inspiration from the Sentential Decision Diagram [68? ].

7 Acknowledgements
We thank the anonymous reviewers at Quantum for their valuable suggestions, which helped
to greatly improve the presentation of the manuscript. We thank Dan Browne for help withestablishing stabilizer ranks. We thank Marie Anastacio, Jonas Helsen, Yash Patel and Matthijs
Rijlaarsdam for their feedback on early versions of the manuscript. We thank Patrick Emonts andAdriàn Pérez-Salinas for discussions on MPS, and Kenneth Goodenough for useful discussions ingeneral. The second author acknowledges the QIA project (funded by European Union’s Horizon
2020, Grant Agreement No. 820445). The third and fourth author are funded by the Netherlands
Organization for Scientific Research (NWO/OCW), as part of the Quantum Software Consortium
program (Project No. 024.003.037/3368). The last author is funded by the research program VENI
with project number 639.021.649 of the Netherlands Organization for Scientific Research (NWO).

References[1] Alwin Zulehner and Robert Wille. “One-pass design of reversible circuits: Combining embedding and synthesis for reversible logic”. IEEE Transactions on Computer-Aided Design ofIntegrated Circuits and Systems 37, 996–1008 (2017).
[2] Lukas Burgholzer and Robert Wille. “Improved DD-based equivalence checking of quantum
circuits”. In 2020 25th Asia and South Pacific Design Automation Conference (ASP-DAC).
Pages 127–132. IEEE (2020).
[3] Lukas Burgholzer, Richard Kueng, and Robert Wille. “Random stimuli generation for theverification of quantum circuits”. In Proceedings of the 26th Asia and South Pacific Design
Automation Conference. Pages 767–772. (2021).
[4] Lukas Burgholzer and Robert Wille. “Advanced equivalence checking for quantum circuits”.
IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 40, 1810–
1824 (2020).
[5] John Preskill. “Quantum computing in the NISQ era and beyond”. Quantum 2, 79 (2018).
[6] Daniel Gottesman.
“The Heisenberg representation of quantum computers” (1998).
url: arxiv.org/abs/quant-ph/9807006.
[7] Scott Aaronson and Daniel Gottesman. “Improved simulation of stabilizer circuits”. PhysicalReview A70 (2004).
[8] Daniel Gottesman. “Stabilizer codes and quantum error correction”. PhD thesis. California
Institute of Technology. (1997).
[9] Maarten Van den Nest, Jeroen Dehaene, and Bart De Moor. “Local unitary versus local
Clifford equivalence of stabilizer states”. Phys. Rev. A 71, 062323 (2005).
[10] Matthias Englbrecht and Barbara Kraus. “Symmetries and entanglement of stabilizer states”.
Phys. Rev. A 101, 062302 (2020).
[11] Robert Raussendorf and Hans J. Briegel. “A one-way quantum computer”. Phys. Rev. Lett.
86, 5188–5191 (2001).
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

31

[12] Sergey Bravyi, Graeme Smith, and John A. Smolin. “Trading classical and quantum computational resources”. Phys. Rev. X 6, 021043 (2016).
[13] Sergey Bravyi and David Gosset. “Improved classical simulation of quantum circuits dominated by Clifford gates”. Phys. Rev. Lett. 116, 250501 (2016).
[14] Sergey Bravyi, Dan Browne, Padraic Calpin, Earl Campbell, David Gosset, and MarkHoward. “Simulation of quantum circuits by low-rank stabilizer decompositions”. Quantum
3, 181 (2019).
[15] Yifei Huang and Peter Love. “Approximate stabilizer rank and improved weak simulation of
Clifford-dominated circuits for qudits”. Phys. Rev. A 99, 052307 (2019).
[16] Lucas Kocia and Peter Love. “Stationary phase method in discrete Wigner functions andclassical simulation of quantum circuits”. Quantum 5, 494 (2021).
[17] Lucas Kocia and Mohan Sarovar. “Classical simulation of quantum circuits using fewer gaussian eliminations”. Physical Review A 103, 022603 (2021).
[18] Sheldon B. Akers. “Binary decision diagrams”. IEEE Computer Architecture Letters 27,
509–516 (1978).
[19] Randal E. Bryant. “Graph-based algorithms for Boolean function manipulation”. IEEE Trans.
Computers 35, 677–691 (1986).
[20] Randal E Bryant and Yirng-An Chen. “Verification of arithmetic circuits with binary moment
diagrams”. In 32nd Design Automation Conference. Pages 535–541. IEEE (1995).
[21] G.F. Viamontes, I.L. Markov, and J.P. Hayes. “High-performance QuIDD-based simulationof quantum circuits”. In Proceedings Design, Automation and Test in Europe Conference andExhibition. Volume 2, pages 1354–1355 Vol.2. (2004).
[22] R. I. Bahar, E. A. Frohm, C. M. Gaona, G. D. Hachtel, E. Macii, A. Pardo, and F. Somenzi.
“Algebraic decision diagrams and their applications”. In Proceedings of 1993 International
Conference on Computer Aided Design (ICCAD). Pages 188–191. (1993).
[23] George F Viamontes, Igor L Markov, and John P Hayes. “Improving gate-level simulation ofquantum circuits”. Quantum Information Processing 2, 347–380 (2003).
[24] Masahiro Fujita, Patrick C. McGeer, and JC-Y Yang. “Multi-terminal binary decision diagrams: An efficient data structure for matrix representation”. Formal methods in systemdesign 10, 149–169 (1997).
[25] E. M. Clarke, K. L. McMillan, X Zhao, M. Fujita, and J. Yang. “Spectral transforms forlarge boolean functions with applications to technology mapping”. In Proceedings of the
30th International Design Automation Conference. Pages 54–60. DAC ’93New York, NY,
USA (1993). Association for Computing Machinery.
[26] Scott Sanner and David McAllester.
“Affine algebraic decision diagrams (AADDs)
and their application to structured probabilistic inference”.
In Proceedings ofthe 19th International Joint Conference on Artificial Intelligence.
Pages 1384–
1390.
IJCAI’05San Francisco, CA, USA (2005). Morgan Kaufmann Publishers Inc.
url: www.ijcai.org/Proceedings/05/Papers/1439.pdf.
[27] D Michael Miller and Mitchell A Thornton. “QMDD: A decision diagram structure for reversible and quantum circuits”. In 36th International Symposium on Multiple-Valued Logic(ISMVL’06). Pages 30–30. IEEE (2006).
[28] Alwin Zulehner and Robert Wille. “Advanced simulation of quantum computations”.
IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 38, 848–
859 (2018).
[29] Xin Hong, Xiangzhen Zhou, Sanjiang Li, Yuan Feng, and Mingsheng Ying. “A tensor network
based decision diagram for representation of quantum circuits”. ACM Trans. Des. Autom.
Electron. Syst.27 (2022).
[30] Stefan Hillmich, Richard Kueng, Igor L. Markov, and Robert Wille. “As accurate as needed,
as efficient as possible: Approximations in DD-based quantum circuit simulation”. In Design, Automation & Test in Europe Conference & Exhibition, DATE 2021, Grenoble, France,
February 1-5, 2021. Pages 188–193. IEEE (2021).
[31] George F Viamontes, Igor L Markov, and John P Hayes. “Quantum circuit simulation”.
Springer Science & Business Media. (2009).

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

32

[32] Xin Hong, Mingsheng Ying, Yuan Feng, Xiangzhen Zhou, and Sanjiang Li. “Approximate
equivalence checking of noisy quantum circuits”. In 2021 58th ACM/IEEE Design Automation
Conference (DAC). Pages 637–642. (2021).
[33] Hans J. Briegel and Robert Raussendorf. “Persistent entanglement in arrays of interacting
particles”. Phys. Rev. Lett. 86, 910–913 (2001).
[34] Wolfgang Dür, Guifre Vidal, and J Ignacio Cirac. “Three qubits can be entangled in two
inequivalent ways”. Physical Review A 62, 062314 (2000).
[35] Eric Chitambar, Debbie Leung, Laura Mančinska, Maris Ozols, and Andreas Winter. “Everything you always wanted to know about LOCC (but were afraid to ask)”. Communicationsin Mathematical Physics 328, 303–326 (2014).
[36] Steven R White. “Density matrix formulation for quantum renormalization groups”. Physicalreview letters 69, 2863 (1992).
[37] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac. “Matrix product state representations”. Quantum Information & Computation 7, 401–430 (2007).
[38] Guifré Vidal. “Efficient classical simulation of slightly entangled quantum computations”.
Physical review letters 91, 147902 (2003).
[39] Adnan Darwiche and Pierre Marquis. “A knowledge compilation map”. Journal of Artificial
Intelligence Research 17, 229–264 (2002).
[40] Karl S Brace, Richard L Rudell, and Randal E Bryant. “Efficient implementation of a BDD
package”. In Proceedings of the 27th ACM/IEEE design automation conference. Pages 40–45.
(1991).
[41] Donald Ervin Knuth. “The art of computer programming. volume 4, fascicle 1”. AddisonWesley. (2005).
[42] Fabio Somenzi. “Efficient manipulation of decision diagrams”. International Journal on Software Tools for Technology Transfer 3, 171–181 (2001).
[43] Koenraad M R Audenaert and Martin B Plenio. “Entanglement on mixed stabilizer
states: normal forms and reduction procedures”. New Journal of Physics 7, 170 (2005).
url: http://stacks.iop.org/1367-2630/7/i=1/a=170.
[44] Marc Hein, Wolfgang Dür, Jens Eisert, Robert Raussendorf, M Nest, and H. J. Briegel.
“Entanglement in graph states and its applications”. In Proceedings of the International
School of Physics ”Enrico Fermi”. Volume Volume 162: Quantum Computers, Algorithms
and Chaos. IOS Press (2006).
[45] Scott Aaronson. “Multilinear formulas and skepticism of quantum computing”. In Proceedingsof the Thirty-Sixth Annual ACM Symposium on Theory of Computing. Page 118–127. STOC
’04New York, NY, USA (2004). Association for Computing Machinery.
[46] Sergey Bravyi and Alexei Kitaev. “Universal quantum computation with ideal Clifford gatesand noisy ancillas”. Phys. Rev. A 71, 022316 (2005).
[47] Charles H Bennett, Herbert J Bernstein, Sandu Popescu, and Benjamin Schumacher. “Concentrating partial entanglement by local operations”. Physical Review A 53, 2046 (1996).
[48] David Y Feinstein and Mitchell A Thornton. “On the skipped variables of quantum multiplevalued decision diagrams”. In 2011 41st IEEE International Symposium on Multiple-Valued
Logic. Pages 164–169. IEEE (2011).
[49] Richard J Lipton, Donald J Rose, and Robert Endre Tarjan. “Generalized nested dissection”.
SIAM journal on numerical analysis 16, 346–358 (1979).
[50] M. Van den Nest, W. Dür, G. Vidal, and H. J. Briegel. “Classical simulation versus universality
in measurement-based quantum computation”. Phys. Rev. A 75, 012337 (2007).
[51] Vı́t Jelı́nek. “The rank-width of the square grid”. Discrete Applied Mathematics 158, 841–
850 (2010).
[52] Hélene Fargier, Pierre Marquis, Alexandre Niveau, and Nicolas Schmidt. “A knowledge compilation map for ordered real-valued decision diagrams”. In Proceedings of the AAAI Conferenceon Artificial Intelligence. Volume 28. (2014).
[53] Robert W Floyd. “Assigning meanings to programs”. In Program Verification. Pages 65–81.
Springer (1993).
[54] JW De Bakker and Lambert G. L. T. Meertens. “On the completeness of the inductive
assertion method”. Journal of Computer and System Sciences 11, 323–357 (1975).

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

33

[55] Ingo Wegener. “Branching programs and binary decision diagrams: theory and applications”.
SIAM. (2000).
[56] James McClung. “Constructions and applications of W-states”. PhD thesis. Worcester
Polytechnic Institute. (2020).
[57] Srinivasan Arunachalam, Sergey Bravyi, Chinmay Nirkhe, and Bryan O’Gorman. “The parameterized complexity of quantum verification” (2022).
[58] Aleks Kissinger and John van de Wetering. “Reducing T-count with the ZX-calculus” (2019).
arXiv:1903.10477.
[59] Himanshu Thapliyal, Edgard Munoz-Coreas, TSS Varun, and Travis S Humble. “Quantum
circuit designs of integer division optimizing T-count and T-depth”. IEEE Transactions on
Emerging Topics in Computing 9, 1045–1056 (2019).
[60] Wang Jian, Zhang Quan, and Tang Chao-Jing. “Quantum secure communication scheme with
W state”. Communications in Theoretical Physics 48, 637 (2007).
[61] Wen Liu, Yong-Bin Wang, and Zheng-Tao Jiang. “An efficient protocol for the quantum
private comparison of equality with W state”. Optics Communications 284, 3160–3163 (2011).
[62] Victoria Lipinska, Gláucia Murta, and Stephanie Wehner. “Anonymous transmission in anoisy quantum network using the W state”. Phys. Rev. A 98, 052320 (2018).
[63] Paul Tafertshofer and Massoud Pedram. “Factored edge-valued binary decision diagrams”.
Formal Methods in System Design 10, 243–270 (1997).
[64] Meghana Sistla, Swarat Chaudhuri, and Thomas Reps. “CFLOBDDs: Context-free-language
ordered binary decision diagrams” (2023). arXiv:2211.06818.
[65] Meghana Sistla, Swarat Chaudhuri, and Thomas Reps. “Symbolic quantum simulation
with quasimodo”. In Constantin Enea and Akash Lal, editors, Computer Aided Verification.
Pages 213–225. Cham (2023). Springer Nature Switzerland.
[66] Rajeev Alur and P. Madhusudan. “Visibly pushdown languages”. In Proceedings of theThirty-Sixth Annual ACM Symposium on Theory of Computing. Pages 202–211. STOC
’04New York, NY, USA (2004). Association for Computing Machinery.
[67] Meghana Sistla, Swarat Chaudhuri, and Thomas Reps. “Weighted context-free-language ordered binary decision diagrams” (2023). arXiv:2305.13610.
[68] Adnan Darwiche. “SDD: a new canonical representation of propositional knowledge bases”.
In Proceedings of the Twenty-Second international joint conference on Artificial IntelligenceVolume Volume Two. Pages 819–826. AAAI Press (2011).
[69] Doga Kisa, Guy Van den Broeck, Arthur Choi, and Adnan Darwiche. “Probabilistic sentential
decision diagrams”. In Proceedings of the Fourteenth International Conference on Principlesof Knowledge Representation and Reasoning. Pages 558–567. KR’14. AAAI Press (2014).
url: cdn.aaai.org/ocs/8005/8005-36908-1-PB.pdf.
[70] Kengo Nakamura, Shuhei Denzumi, and Masaaki Nishino. “Variable shift SDD: A more
succinct sentential decision diagram”. In Simone Faro and Domenico Cantone, editors, 18th
International Symposium on Experimental Algorithms (SEA 2020). Volume 160 of Leibniz International Proceedings in Informatics (LIPIcs), pages 22:1–22:13. Dagstuhl, Germany (2020).
Schloss Dagstuhl–Leibniz-Zentrum für Informatik.
[71] Wolfgang Gunther and Rolf Drechsler. “Minimization of bdds using linear transformations
based on evolutionary techniques”. In 1999 IEEE International Symposium on Circuits andSystems (ISCAS). Volume 1, pages 387–390. IEEE (1999).
[72] Barbara M. Terhal and David P. DiVincenzo. “Classical simulation of noninteracting-fermion
quantum circuits”. Phys. Rev. A 65, 032325 (2002).
[73] Richard Jozsa and Akimasa Miyake. “Matchgates and classical simulation of quantum
circuits”. Proceedings: Mathematical, Physical and Engineering SciencesPages 3089–
3106 (2008).
[74] Martin Hebenstreit, Richard Jozsa, Barbara Kraus, and Sergii Strelchuk. “Computational
power of matchgates with supplementary resources”. Physical Review A 102, 052604 (2020).
[75] Román Orús. “A practical introduction to tensor networks: Matrix product states and projected entangled pair states”. Annals of Physics 349, 117–158 (2014).
[76] Bob Coecke and Ross Duncan. “Interacting quantum observables: categorical algebra anddiagrammatics”. New Journal of Physics 13, 043016 (2011).

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

34

[77] Renaud Vilmart. “Quantum multiple-valued decision diagrams in graphical calculi” (2021).
arXiv:2107.01186.
[78] Richard Rudell. “Dynamic variable ordering for ordered binary decision diagrams”. In Proceedings of 1993 International Conference on Computer Aided Design (ICCAD). Pages 42–47.
IEEE (1993).
[79] Ewout van den Berg and Kristan Temme. “Circuit optimization of Hamiltonian simulation
by simultaneous diagonalization of Pauli clusters”. Quantum 4, 322 (2020).
[80] Eugene M Luks, Ferenc Rákóczi, and Charles RB Wright. “Some algorithms for nilpotent
permutation groups”. Journal of Symbolic Computation 23, 335–354 (1997).
[81] Pavol Ďuriš, Juraj Hromkovič, Stasys Jukna, Martin Sauerhoff, and Georg Schnitger. “On
multi-partition communication complexity”. Information and computation 194, 49–75 (2004).
[82] Hector J. Garcia, Igor L. Markov, and Andrew W. Cross. “Efficient inner-product algorithm
for stabilizer states” (2012). arXiv:1210.6646.
[83] “Stabranksearcher: code for finding (upper bounds to) the stabilizer rank of a quantum state”.
https://github.com/timcp/StabRankSearcher (2021).
[84] Padraic Calpin. “Exploring quantum computation through the lens of classical simulation”.
PhD thesis. UCL (University College London). (2020).

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

35

A Linear-algebra algorithms for Pauli operators
In Sec. 2, we defined the stabilizer group for an n-qubit state |φ⟩ as the group of Pauli operators
A ∈ Paulin which stabilize |φ⟩, i.e. A |φ⟩ = |φ⟩. Here, we explain existing efficient algorithms
for solving various tasks regarding stabilizer groups (whose elements commute with each other).
We also outline how the algorithms can be extended and altered to work for general PauliLIMs,
which do not necessarily commute. For sake of clarity, in the explanation below we first ignore
the scalar λ of a PauliLIM or Pauli element λP . At the end, we explain how the scalars can betaken into account when we use these algorithms as subroutine in LIMDD operations.
Any n-qubit Pauli string can (modulo factor ∈ {±1, ±i}) be written as (X xn Z zn ) ⊗ · · · ⊗ (X x1 Z z1 )for bits xj , zj , 1 ≤ j ≤ n. We can therefore write an n-qubit Pauli string P as a length-2n binary
vector as follows [7],(xn , xn−1 , . . . x1 | zn , zn−1 , . . . , z1 ),
{z
} |
{z
}
|
X block

Z block

where we added the horizontal bar (|) only to guide the eye. We will refer to such vectors as check
vectors. For example, we have X ∼ (1, 0) and Z ⊗ Y ∼ (0, 1|1, 1) . This equivalence induces anordering on Pauli strings following the lexicographic ordering on bit strings. For example, X < Ybecause (1|0) < (1|1) and Z ⊗ I < Z ⊗ X because (00|10) < (01|10).
A set of k Pauli strings thus can be written as 2n × k binary matrix, often called check matrix, as
the following example shows.
 X ⊗ X ⊗ X
1 1 1 | 0 0 0∼.
I ⊗ Z ⊗ Y
0 0 1 | 0 1 1
Furthermore, if P, Q are Pauli strings corresponding to binary vectors (⃗xP , ⃗zP ) and (⃗xQ , ⃗zQ ), then
P ·Q∝

n 
O
j=1

P

P

X xj Z zjn   Q Q O

Q
Q
P
P
X xj Z zj =X xj ⊕xj Z zj ⊕zj
j=1

and therefore the group of n-qubit Pauli strings with multiplication (disregarding factors) is groupisomorphic to the vector space {0, 1}2n (i.e., F2n2 ) with bitwise addition ⊕ (i.e., exclusive or; ‘xor’).
Consequently, many efficient algorithms for linear-algebra problems carry over to sets of Paulistrings. In particular, if G = {g1 , . . . , gk } are length−2n binary vectors (/ n-qubit Pauli strings)
with k ≤ n, then we can efficiently perform the following operations.
RREF: bring G into a reduced-row echelon form (RREF) using Gauss-Jordan elimination (bothare standard linear algebra notions) where each row in the check matrix has strictly moreleading zeroes than the row above. The RREF is achievable by O(k 2 ) row additions (/ multiplications modulo factor) and thus O(k 2 · n) time (see [79] for a similar algorithm). In the
RREF, the first 1 after the leading zeroes in a row is called a ‘pivot’.
Construct Minimal-size Generator Set convert G to a (potentially smaller) set G′ by performing the RREF procedure and discarding resulting all-zero rows. It holds that ⟨G⟩ = ⟨G′ ⟩,i.e., these sets generate the same group modulo phase.
Membership: determining whether a given a vector (/ Pauli string) h has a decomposition inelements of G. This can be done by obtaining minimal-size generating sets H1 , H2 for thesets G and G ∪ {h}, respectively. Then the generating sets have the same number of elements(i.e., rows) if and only if h ∈ ⟨G⟩; otherwise, if h ̸∈ ⟨G⟩, it holds that |H2 | = |H1 | + 1.
Intersection: determine all Pauli strings which, modulo a factor, are contained in both GA and
GB , where GA , GB are generator sets for n-qubit stabilizer subgroups. More specifically, we
obtain the generator set of this group, i.e., we obtain a set GC such that ⟨GC ⟩ = ⟨GA ⟩∩⟨GB ⟩.
This can be achieved using the Zassenhaus algorithm [80] for computing the intersection of
two subspaces of a vector space, in time O(n3 ).
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

36

Division remainder: given a vector h (/ Pauli string h), determine hrem := ming∈⟨G⟩ {g ⊕ h}(minimum in the lexicographic ordering). We do so in the check matrix picture by bringing
G into RREF, and then making the check vector of h contain as many zeroes as possible byadding rows from G:
1: for column index j = 1 to 2n do
2:if hj = 1 and G has a row gi with its pivot at position j then h := h ⊕ gi

The resulting h is hrem . This algorithm’s runtime is dominated by the RREF step; O(n3 ).
To include the scalar into the representation, we remark that Pauli LIMs that appear as labels on
diagrams may have λ ∈ C, i.e., any complex number is allowed. Therefore, to store LIMs, we use
a minor extension to the check vector form introduced above, in order to also include the phase.
Specifically, the phase is stored using two real numbers, by writing λ = r · eiθ with r ∈ R>0 andθ ∈ [0, 2π). Consequently, the check vector has 2n + 2 entries, where the last entries store r and θ,
e.g.:
 


1 1 1 | 0 0 0 | 3 0
3X ⊗ X ⊗ X∼
0 0 1 | 0 1 1 | 12 3π− 21 iI ⊗ Z ⊗ Y
2where we used 3 = 3 · ei·0 and − 12 i = 12 · e3πi/2 . This extended check vector also easily allows
a total ordering, namely, we simply use the ordering on real numbers for r and θ. For example,(1, 1, |0, 0|2, 12 ) < (1, 1|1, 0|3, 0). Let us stress that the factor encoding (r, θ) is less significant than
the Pauli string encoding (xn , . . . , x1 |zn , . . . , z1 ). As a consequence, we can greedily determine theminimum of two Pauli operators, by reading their check vectors from left to right.
Finally, we emphasize that the algorithms above rely on bitwise xor-ing, which is a commutative operation. Since conventional (i.e., factor-respecting) multiplication of Pauli operators is notcommutative, the algorithms above are not straightforwardly applicable to arbitrary PauliLIMn
input. (When the input consist of pairwise commuting Pauli operators, such as stabilizer subgroups [7], the algorithms can be made to work by adjusting row addition to keep track of thescalar.) Fortunately, since Pauli strings either commute or anti-commute, row addition may only
yield factors up to the ± sign, not the resulting Pauli strings. This feature, combined with thestipulated order assigning least significance to the factor, enables us to invoke the algorithms above
as subroutine. We do so in Sec. 4.2.1 and Sec. D.1.

B Proof that cluster states and coset states need exponentially-large
QMDDs
In this appendix, we show that QMDDs which represent both clusters states, and coset states,
are exponentially large in the worst case (respectively, Th. 2 and Corollary 11). On the other
hand, in App. C, we will show that these states can be represented using only O(n) nodes by ⟨X⟩LIMDDs, showing that they are exponentially more succinct than QMDDs. We first fix notationand definitions, after which we prove the theorem using two lemmas.
Let G be an undirected graph with vertices VG = {v1 , . . . , vn } and edge set EG ⊆ VG × VG . For asubset of vertices S ⊆ VG , the S-induced subgraph of G has vertices S and edge set (S × S) ∩ E.
Given G, its graph state |G⟩ is expressed asX|G⟩ =(−1)fG (⃗x) |⃗x⟩(11)⃗
x∈{0,1}nwhere fG (⃗x) is the number of edges in the S-induced subgraph of G.
For a function f : {0, 1}n → C and bit string ⃗a = a1 · · · ak ∈ {0, 1}k , we denote by f⃗a thesubfunction of f restricted to ⃗a:f⃗a (xk+1 , . . . , xn ) := f (a1 , . . . , ak , xk+1 , . . . , xn )Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

(12)
37

We also say that f⃗a is a subfunction of f of order |⃗a| = k.
We will also need the notions of boundary and strong matching.
Definition 6 (Boundary). For a set S ⊆ VG of vertices in G, the boundary of S is the set ofvertices in S adjacent to a vertex outside of S.
Definition 7 (Strong Matching). Let G = (V, E) be an undirected graph. A strong matching is
a subset of edges M ⊆ E that do not share any vertices (i.e., it is a matching) and no two edgesof M are incident to the same edge of G, i.e., an edge in E \ M . Alternatively, a strong matchingis a matching M s.t. G[V (M )] = M . We say that M is an (S, T )-strong matching for two sets ofvertices S, T ⊂ V if M ⊆ S × T . For a strong matching M and a vertex v ∈ V (M ), we let M (v)denote the unique vertex to which v is matched by M .
Using these definitions and notation, we prove Th. 2.
Theorem 2. Denote by |Gn ⟩ the two-dimensional cluster state, defined as a graph state on then × n lattice. Each QMDD representing |Gn ⟩ has at least 2⌊n/12⌋ nodes.
Proof. Let G = lattice(n, n) be the undirected graph of the n × n lattice, with vertex set V ={v1 , . . . , vn2 }. Let σ = v1 v2 · · · vn2 be a variable order, and let S = {v1 , v2 , . . . , v 12 n2 } ⊂ V be thefirst 12 n2 vertices in this order.
The proof proceeds broadly as follows. First, in Lemma 1, we show that any (S, S)-strong matching
M effects 2|M | different subfunctions of fG . Second, Lemma 2 shows that the lattice contains alarge (S, S)-strong matching for any choice of S. Put together, this will prove the lower bound on
the number of QMDD nodes as in Th. 2 by the fact that a QMDD for the cluster state G has anode per unique subfunction of the function fG . Fig. 13 illustrates this setup for the 5 × 5 lattice.
Lemma 1. Let M be a non-empty (S, S)-strong matching for the vertex set S chosen above. If
σ = v1 v2 · · · vn2 is a variable order where all vertices in S appear before all vertices in S, thenfG (x1 , . . . , xn2 ) has 2|M | different subfunctions of order |S|.
Proof. Let SM := S ∩ V (M ) and S M := S ∩ M be the sets of vertices that are involved in thestrong matching. Write χ(x1 , . . . , xn ) for the indicator function for vertices: χ(x1 , . . . , xn ) :={vi | xi = 1, i ∈ [n]}. Choose two different subsets A, B ⊆ SM and let ⃗a = χ−1 (A) and ⃗b = χ−1 (B)be the corresponding length-|S| bit strings. These two strings induce the two subfunctions fG,⃗aand fG,⃗b . We will show that these subfunctions differ in at least one point.
First, if fG,⃗a (0, . . . , 0) ̸= fG,⃗b (0, . . . , 0), then we are done. Otherwise, take a vertex s ∈ A ⊕ B andsay w.l.o.g. that s ∈ A \ B. Let t = M (s) be its partner in the strong matching. Then we have,|E[A ∪ {t}]| = |E[A]| + 1 but |E[B ∪ {t}]| = |E[B]|. Therefore we havefG,⃗a (0, . . . , 0, xt = 0, 0, . . . , 0) ̸= fG,⃗a (0, . . . , 0, xt = 1, 0, . . . , 0)(13)fG,⃗b (0, . . . , 0, xt = 0, 0, . . . , 0) = fG,⃗b (0, . . . , 0, xt = 1, 0, . . . , 0)

(14)

We see that each subset of SM corresponds to a different subfunction of fG . Since there are 2|M |subsets of M , fG has at least that many subfunctions.
We now show that the n × n lattice contains a large enough strong matching.
Lemma 2. Let S = {v1 , . . . , v 12 n2 } be a set of 12 n2 vertices of the n × n lattice, as above. Then1 
the graph contains a (S, S)-strong matching of size at least 12n .
Proof. Consider the boundary BS of S. This set contains at least n/3 vertices, by Theorem 11in 1 [49]. Each vertex of the boundary of S has degree at most 4. It follows
 1 that there 1 is a set of|B|vertices
which
share
no
neighbors.
In
particular,
there
is
a
set
of
|B
|
≥
S
S
4
4
12 n vertices
in BS which share no neighbors in S.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

38

Put together, every choice of half the vertices in the latticea set with a boundary of at least 1 yields
n/3 nodes, which yields a strong matching of at least 12
n edges, which shows that fG has at
1
least 2⌊ 12 n⌋ subfunctions of order 1 n2 .
2Figure 13: The 5 × 5 lattice, partitioned in a vertex set S and its complement S. A strong matching between
S and S is indicated by black edges.

Proof that coset states need exponentially large QMDDs. We now show that QMDDs
which represent coset states are exponentially large in the worst case. We will use the following
result by Ďuriš et al. on binary decision diagrams (BDDs), which are QMDDs with codomain {0, 1}.
This result concerns vector spaces, but of course, every vector space of {0, 1}n is, in particular, a
coset.
Theorem 10 (Ďuriš et al.[81]). The characteristic function fV : {0, 1}n → {0, 1} of a randomly
chosen vector space V in {0, 1}n , defined as fV (x) = 1 if x ∈ V and 0 otherwise, needs a BDD ofsize 2Ω(n) /(2n) with high probability.
Our result
P follows by noting that if f has codomain {0, 1} as above, then the QMDD of the state|f ⟩ = x f (x) |x⟩ has the same structure as the BDD of f . Consequently, in particular the BDD
and QMDD have the same number of nodes.
Corollary 11. For a random vector space V ⊆ {0, 1}n , the coset state |V ⟩ requires QMDDs of
size 2Ω(n) /(2n) with high probability.
Proof. We will show that the QMDD has the same number of nodes as a BDD. A BDD encodes afunction f : {0, 1}n → {0, 1}. In this case, the BDD encodes fV , the characteristic function of V .
A BDD is a graph which contains one node for each subfunction of f . (In the literature, such a
BDD is sometimes called a Full BDD, so that the term BDD is reserved for a variant where thenodes are in one-to-one correspondence with the subfunctions f which satisfy f0 ̸= f1 ).
P
Similarly, a QMDD representing a state |φ⟩ = x f (x) |x⟩ can be said to represent the function
f : {0, 1}n → C, and contains one node for each subfunction of f modulo scalars. We will show that,two distinct subfunctions of fV are never equal up to a scalar. To this end, let fV,a , fV,b be distinct
subfunctions of fV induced by partial assignments a, b ∈ {0, 1}k . We will show that there is no
λ ∈ C∗ such that fV,a = λfV,b . Since the two subfunctions are not pointwise equal, say that the two
subfunctions differ in the point x ∈ {0, 1}n−k , i.e., fV,a (x) ̸= fV,b (x). Say without loss of generalitythat fV,a (x) = 0 and fV,b (x) = 1. Then, since λ ̸= 0, we have λ = λfS,b (x) ̸= fV,a (x) = 0, so
fV,a ̸= λfB,b .

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

39

Because distinct subfunctions of fV are not equal up to a scalar, the QMDD of |V ⟩ contains a nodefor every unique subfunction of fV . We conclude that, since by Th. 10 with high probability the
BDD representing fV has exponentially many nodes, so does the QMDD representing |V ⟩.

C How to write graph states, coset states and stabilizer states as TowerLIMDDs
In this appendix, we prove that the families of ⟨Z⟩-, ⟨X⟩-, and ⟨Pauli⟩-Tower-LIMDDs correspondto graph states, coset states, and stabilizer states, respectively, in Th. 12, Th. 13 and Th. 1 below.
Def. 5 for reduced Pauli-LIMDDs requires modification for G = ⟨Z⟩-LIMDDs because of the absenceof X as discussed below the definition. Note that the proofs do not rely on the specialized definitionof reduced LIMDDs, but only on Def. 2 which allows parameterization of the LIM G. They only
rely on the Tower LIMDD in Def. 3.
Before we give the proof, we remark that graph states present an interesting special case because
the LIMDD’s edge labels contain meaningful information. Namely, the labels on the high edges
of a graph state’s LIMDD are precisely the edges in the original graph. Specifically, suppose agraph G gives rise to a graph state |φG ⟩ represented by a LIMDD. Let P = Pk−1 ⊗ · · · ⊗ P1 be thelabel on the high edge out of the LIMDD node at level k. Then G contains an edge (vk , vj ) if andonly if Pj = Z (with the roles of k and j reversed if k < j). These edge labels come about in astraightforward manner during the construction of the graph state. Namely, the graph state |φG ⟩⊗nis produced by starting from the state |+⟩ , and applying controlled-Z gates to qubit pairs (u, v)
for every edge (u, v) in the graph. Applying such a controlled-Z gate to qubit pair (u, v) has theeffect of setting Pv to Z in the high edge outgoing from the vertex at level u. In general, however,the labels on the high edges cannot be easily inferred from the stabilizer state.
A G-Tower-LIMDD representing an n-qubit state is a LIMDD which has n nodes, not counting theleaf. It has G-LIMs on its high edges. Def. 3 gives an exact definition.
Theorem 12 (Graph states are ⟨Z⟩-Tower-LIMDDs). Let n ≥ 1. Denote by Gn the set of n-qubit
graph states and write Zn for the set of n-qubit quantum states which are represented byN⟨Z⟩Tower-LIMDDs a defined in Def. 3, i.e, a tower with low-edge-labels I and high-edge labels λ j Pj
with Pj ∈ {I, Z} and λ = 1, except for the root edge where λ ∈ C \ {0}. Then Gn = Zn .
Proof. We establish Gn ⊆ Zn by providing a procedure to convert any graph state in Gn to a
⟨Z⟩-Tower-LIMDD in Zn . See Fig. 14 for an example of a 4-qubit graph state. We describe theprocedure by induction on the number n of qubits in the graph state.
Base case: n = 1. We note that √there is only one single-qubit graph state by definition (seeEq. 11), which is |+⟩ := (|0⟩ + |1⟩)/ 2 and can be represented as LIMDD by a single node (inaddition to the leaf node): see Fig. 14(a).
Induction case. We consider an (n + 1)-qubit graph state |G⟩ corresponding to the graph G. We
isolate the (n + 1)-th qubit by decomposing the full state definition from Eq. 11:








O
1 

|G⟩ = √ |0⟩ ⊗ |G1..n ⟩ + |1⟩ ⊗ 
Zj  |G1..n ⟩2
(n+1,j)∈E


|
{z
}

(15)

Isomorphism B

where E is the edge set of G and G1..n is the induced subgraph of G on vertices 1 to n. Thus,|G1..n ⟩ is an n-qubit graph state on qubits 1 to n. Since |G1..n ⟩ is a graph state on n qubits, by
the induction hypothesis, we have a procedure to convert it to a ⟨Z⟩-Tower-LIMDD ∈ Zn . Now
we construct a ⟨Z⟩-Tower-LIMDD for |G⟩ as follows. The root node has two outgoing edges, both

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

40

a)

A

b)

B

A

c)

d)

B

A

C

A

B

C

D

Z⊗I⊗Z
−1

Z

Z⊗I

Z⊗I

Z

Z

−1
−1

0
1

1

1

1

1

−1

Figure 14: Construction of the ⟨Z⟩-Tower LIMDD for the 4-qubit cluster state, by iterating over the vertices
in the graph, as described in the proof of Th. 12. (a) First, we consider the single-qubit graph state, whichcorresponds to a the subgraph containing only vertex A. (b) Then, we add vertex B, which is connected to
A by an edge. The resulting LIMDD is constructed from the LIMDD from (a) by adding a new root node. In
the figure, the isomorphism is ZB ⊗ IA , since vertex C is connected to vertex B (yielding the Z operator) but
not to A (yielding the identity operator I). (c) This process is repeated for a third vertex C until we reach the
LIMDD of the full 4-qubit cluster state (d). For comparison, (d) also depicts a regular QMDD for the samegraph state, which has width 4 instead of 1 for the LIMDD.

going to the node representing |G1..n ⟩. The node’s low edge has label I, and the node’s high edge
has label B, as follows,
O
B=
Zj
(16)
(n+1,j)∈E

Thus the root node represents the state |0⟩ |G1..n ⟩ + |1⟩ B |G1..n ⟩, satisfying Eq. 15.
To prove Zn ⊆ Gn , we show how to construct the graph corresponding to a given ⟨Z⟩-Tower
LIMDD. Briefly, we simply run the algorithm outlined above in reverse, constructing the graphone node at a time. Here we assume without loss of generality that the low edge of every node islabeled I.
Base case. The LIMDD node above the Leaf node, representing the state |+⟩, always represents
the singleton graph, containing one node.
Induction case. Suppose that the LIMDD node k + 1 levels above the Leaf has a low edge labeledI, and a high edge labeled Pk ⊗ · · · ⊗ P1 , with Pj = Z aj for j = 1 . . . k. Here by Z aj we mean
Z 0 = I and Z 1 = Z. Then we add a node labeled k + 1 to the graph, and connect it to those nodes
j with aj = 1, for j = 1 . . . k. The state represented by this node is of the form given in Eq. 15, so
it represents a graph state.
A simple counting argument based on the above construction shows that |Zn | = |Gn | = 2( 2 ) , sonthe conversion is indeed a bijection. Namely, there are 2( 2 ) graphs, since there are n2 edges tonchoose, and there are 2( 2 ) ⟨Z⟩-Tower-LIMDDs,
because the total number of single-qubit operators
of the LIMs on the high edges is n2 , each of which can be independently chosen to be either I orZ.
n

We now prove that coset states are represented by ⟨X⟩-Tower-LIMDDs.
Theorem 13 (coset states are ⟨X⟩-Tower-LIMDDs). Let n ≥ 1. Denote by Vn the set of n-qubit
coset states and write Xn for the set of n-qubit quantum states which are represented
N by ⟨X⟩Tower-LIMDDs as per Def. 3, i.e., a tower with low edge labels I and high edge labels λ j Pj with
Pj ∈ {I, X} and λ ∈ {0, 1}, except for the root edge where λ ∈ C \ {0}. Then Vn = Xn .

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

41

Proof. We first prove Vn ⊆ Xn by providing a procedure for constructing a Tower-LIMDD for a
coset state. We prove the statement for the case when C is a group rather than a coset; the result
will then follow by noting that, by placing the label X an ⊗ · · · ⊗ X a1 on the root edge, we obtainthe coset state |C + a⟩. The procedure is recursive on the number of qubits.
√
Base case: n = 1. In this case, there are two coset states: |0⟩ and (|0⟩ + |1⟩)/ 2, which are
represented by a single node which has a low and high edge pointing to the leaf node with low/highedge labels 1/0 and 1/1, respectively.
Induction case. Now consider an (n + 1)-qubit coset state |S⟩ for a group S ⊆ {0, 1}n+1 forsome n ≥ 1 and assume we have a procedure to convert any n-qubit coset state into a TowerLIMDD in Xn . We consider two cases, depending on whether the first bit of each element of S iszero:(a) The first bit of each element of S is 0. Thus, we can write S = {0x | x ∈ S0 } for some set
S0 ⊆ {0, 1}n . Then 0a, 0b ∈ S =⇒ 0a ⊕ 0b ∈ S implies a, b ∈ S0 =⇒ a ⊕ b ∈ S0 and thus S0
is an length-n bit string vector space. Thus by assumption, we have a procedure to convert itto a Tower-LIMDD in Xn . Convert it into a Tower-LIMDD in Xn+1 for |S⟩ by adding a fresh
node on top with low edge label I⊗n and high edge label 0, both pointing to the the root S.
(b) There is some length-n bit string u such that 1u ∈ S. Write S as the union of the sets
{0x | x ∈ S0 } and {1x | x ∈ S1 } for sets S0 , S1 ⊆ {0, 1}n . Since S is closed under element-wise
XOR, we have 1u ⊕ 1x = 0(u ⊕ x) ∈ S for each x ∈ S1 and therefore u ⊕ x ∈ S0 for eachx ∈ S1 . This implies that S1 = {u ⊕ x | x ∈ S0 } and thus S is the union of {0x | x ∈ S0 } and{1u ⊕ 0x | x ∈ S0 }. By similar reasoning as in case (a), we can show that S0 is a vector spaceon length-n bit strings.
We build a Tower-LIMDD for |S⟩ as follows. By the induction hypothesis, there is a TowerLIMDD with root node v which represents |v⟩ = |S0 ⟩. We construct a new node whose twooutgoing edges both go to this node v. Its low edge has label I⊗n and its high edge has label
P = Pn ⊗ · · · ⊗ P1 where Pj = X if uj = 1 and Pj = I if uj = 0.
We now show Vn ⊆ Xn , also by induction.
Base case: n = 1. There are only two Tower-LIMDDs on 1 qubit satisfying the description above,
namely(1) A node whose two edges point to the leaf. Its low edge has label 1, and its high edge has label 0.
This node represents the coset state |0⟩, corresponding to the vector space V = {0} ⊆ {0, 1}1 .
(2) A node whose two edges point to the leaf. Its low edge has label 1 and its high edge also
has label 1. This node represents the coset state |0⟩ + |1⟩, corresponding to the vector space
V = {0, 1}.
Induction case. Let v be the root node of an n + 1-qubit Tower ⟨X⟩-LIMDD as described above.
We distinguish two cases, depending on whether v’s high edge has label 0 or not.
(a) The high edge has label 0. Then |v⟩ = |0⟩ |v0 ⟩ for a node v0 , which represents a coset state|v0 ⟩ corresponding to a coset V0 ⊆ {0, 1}n , by the induction hypothesis. Then v corresponds
to the coset {0x | x ∈ V0 }.
(b) the high edge has label P = Pn ⊗ · · · ⊗ P1 with Pj ∈ {I, X}. Then |v⟩ = |0⟩ |v0 ⟩ + |1⟩ ⊗ P |v0 ⟩.
By the observations above, this is a coset state, corresponding to the vector space V = {0x|x ∈V0 } ∪ {1(ux)|x ∈ V0 } where u ∈ {0, 1}n is a string whose bits are uj = 1 if Pj = X and uj = 0
if Pj = I, and V0 is the vector space corresponding to the coset state |v0 ⟩.
Lastly, we prove the stabilizer-state case, showing that they are exactly equivalent to the ⟨Pauli⟩Tower-LIMDD, as defined in Def. 3. For this, we first need Lemma 3 and Lemma 4, which state that,
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

42

if one applies a Clifford gate to a ⟨Pauli⟩-Tower-LIMDD, the resulting state is another ⟨Pauli⟩Tower-LIMDD. First, Lemma 3 treats the special case of applying a gate to the top qubit; thenLemma 4 treats the general case of applying a gate to an arbitrary qubit.
Lemma 3. Let |φ⟩ be an n-qubit stabilizer state which is represented by a ⟨Pauli⟩-Tower-LIMDD
as defined in Def. 3. Let U be either a Hadamard gate or S gate on the top qubit (n-th qubit),or a downward CNOT with the top qubit as control. Then U |φ⟩ is still represented by a ⟨Pauli⟩Tower-LIMDD.
Proof. The proof is on the number n of qubits.

√
Base case: n = 1. For n = 1, there are six single-qubit stabilizer states |0⟩ , |1⟩ and (|0⟩+α |1⟩)/ 2for α ∈ {±1, ±i}. There are precisely represented by Pauli-Tower-LIMDDs with high edge label
factor ∈ {0, ±1, ±i} as follows:• for |0⟩:

1

1

1

• for |1⟩: A ·

0
1

1
0

√
• for (|0⟩ + α |1⟩)/ 2):

1
1

where A ∝ X or A ∝ Y1α1

Since the H and S gate permute these six stabilizer states, U |φ⟩ is represented by a ⟨Pauli⟩Tower-LIMDD if |φ⟩ is.
Induction case. For n > 1, we first consider U = S and U = CNOT. Let R be the labelof the root edge. If U = S, then the high edge of the top node is multiplied with i, while adownward CNOT (target qubit with index k) updates the high edge label A 7→ Xk A. Next, theroot edge label is updated to U RU † , which is still a Pauli string, since U is a Clifford gate. Since
the high labels of the top qubit in the resulting diagram is still a Pauli string, and the high edge’s
weights are still ∈ {0, ±1, ±i}, we conclude that both these gates yield a ⟨Pauli⟩-Tower-LIMDD.
Finally, for the Hadamard, we decompose |φ⟩ = |0⟩ ⊗ |ψ⟩ + α |1⟩ ⊗ P |ψ⟩ for some (n − 1)-qubit
stabilizer state |ψ⟩, α ∈ {0, ±1, ±i} and P is an (n − 1)-qubit Pauli string. Now we note thatH |φ⟩ ∝ |0⟩ ⊗ |ψ0 ⟩ + |1⟩ ⊗ |ψ1 ⟩ where |ψx ⟩ := (I + (−1)x αP ) |ψ⟩ with x ∈ {0, 1}. Now we consider
two cases, depending on whether P commutes with all stabilizers of |ψ⟩:(a) There exist a stabilizer g of |ψ⟩ which anticommutes with P . We note two things. First,⟨ψ|P |ψ⟩ = ⟨ψ|P g|ψ⟩ = ⟨ψ|g · (−P )|ψ⟩ = −⟨ψ|P |ψ⟩, hence ⟨ψ|P |ψ⟩ = 0. It follows from
Lemma 15 of [82] that |ψx ⟩ is a stabilizer state, so by the induction hypothesis it can be
written as a ⟨Pauli⟩-Tower-LIMDD. Let v be the root node of this LIMDD. Next, we note thatg |ψ0 ⟩ = g(I + αP ) |ψ⟩ = (I − αP )g |ψ⟩ = |ψ1 ⟩. Hence,⟨Pauli⟩-Tower-LIMDD for H |φ⟩.

v

g

I

v

is the root node of a(b) All stabilizers of |ψ⟩ commute with P . Then (−1)y P is a stabilizer of |ψ⟩ for either y = 0 or
y = 1. Hence, |ψx ⟩ = (I +(−1)x αP ) |ψ⟩ = (1+(−1)x+y α) |ψ⟩. Therefore, |φ⟩ = |a⟩⊗|ψ⟩ where|a⟩ := (1 + (−1)y α) |0⟩ + (1 + (−1)y+1 α |1⟩). It is not hard to see that |a⟩ is a stabilizer state forall choices of α ∈ {0, ±1, ±i}. By the induction hypothesis, both |a⟩ and |ψ⟩ can be represented
as ⟨Pauli⟩-Tower-LIMDDs. We construct a ⟨Pauli⟩-Tower-LIMDD for H |φ⟩ by replacing theleaf of the LIMDD of |a⟩ by the root node of the LIMDD of |ψ⟩, and propagating the root edge
label of |ψ⟩ upwards. Specifically, if the root edge of |a⟩ isand if the root edge of |ψ⟩ is
I

w

βI

w

B

A

v

with v =

1

1

β

1 ,

w , then a ⟨Pauli⟩-Tower-LIMDD for H |φ⟩ has root nodeand has root edge label A ⊗ B.

Lemma 4. Let |φ⟩ be an n-qubit state state represented by a ⟨Pauli⟩-Tower-LIMDD, as definedin Def. 3. Let U be either a Hadamard gate, an S gate or a CNOT gate. Then U |φ⟩ is a state
which is also represented by a ⟨Pauli⟩-Tower-LIMDD.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

43

Proof. The proof is by induction on n. The case n = 1 is covered by Lemma 3. Suppose that theinduction hypothesis holds, and let |φ⟩ be an n + 1-qubit state represented by a ⟨Pauli⟩-TowerLIMDD. First, we note that a CNOT gate CXct can be written as CXct = (H ⊗ H)CXtc (H ⊗ H),
so without loss of generality we may assume that c > t. We treat two cases, depending on whether
U affects the top qubit or not.
(a) U affects the top qubit. Then U |φ⟩ is represented by a ⟨Pauli⟩-Tower-LIMDD, according toLemma 3.
(b) U does not affect the top qubit. Suppose |φ⟩ = |0⟩ ⊗ |φ0 ⟩ + |1⟩ ⊗ αP |φ0 ⟩ (with P a Pauli stringand α ∈ {0, ±1, ±i}). Then U |φ⟩ = |0⟩ ⊗ U |φ0 ⟩ + |1⟩ ⊗ (αU P U † )U |φ0 ⟩. Since U is either
a Hadamard, S gate or CNOT, and |φ0 ⟩ is an n-qubit state, the induction hypothesis statesAv be the root edgethat the state U |φ0 ⟩ is represented by a ⟨Pauli⟩-Tower-LIMDD. Letof this ⟨Pauli⟩-Tower-LIMDD, representing U |φ0 ⟩. Then U |φ⟩ is represented by the root edgeI⊗Aw , where w is the nodeI

αA−1 U P U † A. The label αA−1 U P U † A is a Pauli
LIM, and may therefore be used as the label on the high edge of w.
v

v

Finally, we show that stabilizer states are precisely the ⟨Pauli⟩-Tower-LIMDDs.
Theorem 1. Let n > 0. Each n-qubit stabilizer state is represented up to normalization by a⟨Pauli⟩-Tower LIMDDs of Def. 3, e.g., where the scalars λ of the PauliLIMs λP on high edges
are restricted as λ ∈ {0, ±1, ±i}. Conversely, every such LIMDD represents a stabilizer state.
Proof. We first prove that each stabilizer state is represented by a ⟨Pauli⟩-Tower-LIMDD. We
recall that each stabilizer state can be obtained as the output state of a Clifford circuit on input⊗nstate |0⟩ . Each Clifford circuit can be decomposed into solely the gates H, S and CNOT. The⊗nstate |0⟩is represented by a ⟨Pauli⟩-Tower-LIMDD. According to Lemma 4, applying an H, S
or CNOT gate to a ⟨Pauli⟩-Tower-LIMDD results a state represented by another ⟨Pauli⟩-TowerLIMDD. One can therefore apply the gates of a Clifford circuit to the initial state |0⟩, and obtain a⟨Pauli⟩-Tower-LIMDD for every intermediate state, including the output state. Therefore, every
stabilizer state is represented by a ⟨Pauli⟩-Tower-LIMDD.
For the converse direction, the proof is by induction on n. We only need to note that a state
represented by a ⟨Pauli⟩-Tower-LIMDD can be written as |φ⟩ = |0⟩ ⊗ |φ0 ⟩ + |1⟩ ⊗ αP |φ0 ⟩ =C(P )(|0⟩ + α |1⟩) ⊗ |φ0 ⟩ where C(P ) := |0⟩⟨0| ⊗ I + |1⟩⟨1| ⊗ P is the controlled-(P ) gate. Using therelations Z = HXH, Y = SXS † and S = Z 2 , we can decompose C(P ) as CNOT, H and S, henceC(P ) is a Clifford gate. Since both |0⟩ + α |1⟩ and |φ0 ⟩ can be written as ⟨Pauli⟩-Tower-LIMDDs,
they are stabilizer states by the induction hypothesis. Therefore, the state |ψ⟩ = (|0⟩ + α |1⟩) ⊗ |φ0 ⟩is also a stabilizer state. Thus, the state |φ⟩ = C(P ) |ψ⟩ is obtained by applying the Clifford gate
C(P ) to the stabilizer state |φ⟩. Therefore, |φ⟩ is a stabilizer state.

D Efficient algorithms for choosing a canonical high label
I

λP

v1
Here, we present an efficient algorithm which, on input Pauli-LIMDD node v0 w
, returns
a canonical choice for the high label Bhigh (algorithm GetLabels, in Alg. 12). By canonical, wemean that it returns the same high label for any two nodes in the same isomorphism equivalenceclass, i.e., for any two nodes v, w for which |v⟩ ≃Pauli |w⟩.

We first characterize all eligible labels Bhigh in terms of the stabilizer subgroups of the children
nodes v0 , v1 , denoted as Stab(v0 ) and Stab(v1 ) (see Sec. 2 for the definition of stabilizer subgroup).
Then, we provide the algorithm GetLabels which correctly finds the lexicographically minimal
eligible label (and corresponding root label), and runs in time O(n3 ) where n is the number of
qubits.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

44

Fig. 15 illustrates this process. In the figure, the left node w summarizes the status of the MakeEdge algorithm on Line 10, when this algorithm has enough information to construct the
I⊗n

λP

v1
w
semi-reduced node v0, shown on the left. The node v r , on the right, is the canonicalnode, and is obtained by replacing w’s high edge’s label by the canonical label Bhigh . This labelxis chosen by minimizing the expression Bhigh = (−1)s λ(−1) g0 P g1 , where the minimization is overs, x ∈ {0, 1}, g0 ∈ Stab(|v0 ⟩), g1 ∈ Stab(|v1 ⟩), subject to the constraint that x = 0 if v0 ̸= v1 . We
have |w⟩ ≃Pauli |v r ⟩ by construction as intended, namely, they are related via |w⟩ = Broot |v r ⟩.
Th. 14 shows that this way to choose the high label indeed captures all eligible high labels, i.e., a

node

I

v0

vr

Bhigh

v1

is isomorphic to |w⟩ if and only if Bhigh is of this form.
I⊗n

λP

v1
w
Theorem 14 (Eligible high-edge labels). Let v0
be a semi-reduced n-qubit node
in a Pauli-LIMDD, where v0 , v1 are reduced, P is a Pauli string and λ ̸= 0. For all nodes v =
v0

I⊗n

v

Bhigh

v1

, it holds that |w⟩ ≃ |v⟩ if and only if
Bhigh = (−1)s · λ(−1) g0 P g1x

(17)

for some g0 ∈ Stab(v0 ), g1 ∈ Stab(v1 ), s, x ∈ {0, 1} and x = 0 if v0 ̸= v1 . An isomorphism mapping|w⟩ to |v⟩ is
Broot = (X ⊗ λP )x · (Z s ⊗ (g0 )−1 ).
(18)
Proof. It is straightforward to verify that the isomorphism Broot in eq. (18) indeed maps |w⟩ to|v⟩ (as x = 1 implies v0 = v1 ), which shows that |w⟩ ≃ |v⟩. For the converse direction, suppose
there exists an n-qubit Pauli LIM C such that C |w⟩ = |v⟩, i.e.,
C (|0⟩ ⊗ |v0 ⟩ + λ |1⟩ ⊗ P |v1 ⟩) = |0⟩ ⊗ |v0 ⟩ + |1⟩ ⊗ Bhigh |v1 ⟩ .

(19)

We show that if Bhigh satisfies eq. (19), then it has a decomposition as in eq. (17). We write
C = Ctop ⊗ Crest where Ctop is a single-qubit Pauli operator and Crest is an (n − 1)-qubit Pauli
LIM (or a complex number ̸= 0 if n = 1). We treat the two cases Ctop ∈ {I, Z} and Ctop ∈ {X, Y }separately:0y(a) Case Ctop ∈ {I, Z}. Then Ctop = 10 (−1)
for y ∈ {0, 1}. In this case, Eq. 19 implies
Ctop |0⟩ Crest |v0 ⟩ = |0⟩ |v0 ⟩, so Crest |v0 ⟩ = |v0 ⟩, in other words Crest ∈ Stab(|v0 ⟩). Moreover,−1
Eq. 19 implies (−1)y λCrest P |v1 ⟩ = Bhigh |v1 ⟩, or, equivalently, (−1)−y λ−1 P −1 Crest
Bhigh ∈
Stab(v1 ). Hence, by choosing s = y and x = 0, we computey (−1)0(−1) λCrest P (−1)| {z }|−y −1∈Stab(v0 )

λ

−1

P
{z

0
(−1)y λ(−1)−1Crest Bhigh =Bhigh = Bhighy}∈Stab(v1 )(−1) λx−1
Broot = (λX ⊗ P ) · Z s ⊗ (g0 )
⇝

w

I⊗n

λP
I⊗n
g0 ∈ Stab(v0 )vrBhigh = (−1)s λ(−1) g0 P g1
x

v0

≼

v1

Stab(v1 ) ∋ g1

Choose s, x ∈ {0, 1}, g0 ∈ Stab(v0 ), g1 ∈ Stab(v1 ) s.t. Bhigh is minimal and x = 0 if v0 ̸= v1 .
Figure 15: Illustration of finding a canonical high label for a semi-reduced node w, yielding a reduced node v r .
The chosen high label is the minimal element from the set of eligible high labels based on stabilizers g0 , g1 ofxv0 , v1 (drawn as self loops). The minimal element holds a factor λ(−1) for some x ∈ {0, 1}. There are twocases: if v0 ̸= v1 or x = 0, then the factor is λ and the root edge should be adjusted with an I or Z on theroot qubit. The other case, x = 1, leads to an additional multiplication with an X on the root qubit.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

45

(b) Case Ctop ∈ {X, Y }. Write Ctop =

0 z −1z 0

where z ∈ {1, i}. Now, eq. (19) implies

zCrest |v0 ⟩ = Bhigh |v1 ⟩andz −1 λCrest P |v1 ⟩ = |v0 ⟩ .

(20)

From Eq. 20, we first note that |v0 ⟩ and |v1 ⟩ are isomorphic, so by Corollary 9, and because
the diagram has merged these two nodes, we have v0 = v1 . Consequently, we find from Eq. 20−1that z −1 Crest
Bhigh ∈ Stab(v0 ) and z −1 λCrest P ∈ Stab(v1 ). Now choose x = 1 and choose−1s such that (−1)s · z −2 Crest
Bhigh Crest = Bhigh (recall that Pauli LIMs either commute or
anticommute, so Bhigh Crest = ±Crest Bhigh ). This yields:−1−1
Bhigh Crest = Bhigh(−1)s λ−1 · z −1 Crest
Bhigh ·P · z −1 λP Crest = λ−1 · λ · (−1)s z −2 · Crest
|
{z
}
{z
}
|
∈Stab(v1 )∈Stab(v0 )where we used the fact that P 2 = I⊗(n−1) because P is a Pauli string.

Corollary 15. As a corollary of Th. 14, we find that taking, as in Fig. 15,
HighLabel( v0

I

v

λP

v1

)=

o
n
x
( (−1)s · λ(−1) · g0 · P · g1 x ̸= 1 if v0 ̸= v1 )mini,s,x∈{0,1},gi ∈Stab(vi )yields a proper implementation of HighLabel as required by Def. 5, because it considers all possible
Bhigh such that |v⟩ ≃Pauli |0⟩ |v0 ⟩ + |1⟩ ⊗ Bhigh |v1 ⟩.
A naive implementation for GetLabels would follow the possible decompositions of eligible LIMs
(see Eq. 17) and attempt to make this LIM smaller by greedy multiplication, first with stabilizers
of g0 ∈ Stab(v0 ), and then with stabilizers g1 ∈ Stab(v1 ). To see why this does not work, consider
the following example: the high edge label is Z and the stabilizer subgroups are Stab(v0 ) = ⟨X⟩
and Stab(v1 ) = ⟨Y ⟩. Then the naive algorithm would terminate and return Z because X, Y > Z,
which is incorrect since the high-edge label X · Z · Y = −iI is smaller than Z.
Algorithm 12 Algorithm for finding LIMs Bhigh and Broot required by MakeEdge. Its parameters represent a semi-reduced node

v0

I

v

λP

v1

and it returns LIMs Bhigh , Broot such that

Bhigh

I

v1
|v⟩ = Broot |w⟩ with v0 w
. The LIM Bhigh is chosen canonically as the lexicographically smallest from the set characterized in Th. 14. It runs in O(n3 )-time (with n the number of
qubits), provided GetStabilizerGenSet has been computed for children v0 , v1 (an amortizedcost).

1: procedure GetLabels(PauliLim λP , Node v0 , v1 with λ ̸= 0 and v0 , v1 reduced)

Output: canonical high label Bhigh and root label Broot
G0 , G1 := GetStabilizerGenSet(v0 ), GetStabilizerGenSet(v1 )3:(g0 , g1 ) := ArgLexMin(G0 , G1 , λP )if v0 = v1 then
4:
n
o
x
5:
(x, s) := arg min (−1)s λ(−1) g0 P g12:(x,s)∈{0,1}2

6:
7:
8:

else
x := 0s := arg min {(−1)s λg0 P g1 }s∈{0,1}9:10:Bhigh := (−1)s · λ(−1) · g0 · P · g1
Broot := (X ⊗ λP )x · (Z s ⊗ (g0 )−1 )11:

return (Bhigh , Broot )x

To overcome this, we consider the group closure of both Stab(v0 ) and Stab(v1 ). See Alg. 12 for
the O(n3 )-algorithm for GetLabels, which proceeds in two steps. In the first step (Line 3), we
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

46

use the subroutine ArgLexMin for finding the minimal Pauli LIM A such that A = λP · g0 · g1
for g0 ∈ Stab(v0 ), g1 ∈ Stab(v1 ). We will explain and prove correctness of this subroutine below
in Sec. D.2. In the second step (Line 4-8), we follow Corollary 15 by also minimizing over x
and s. Finally, the algorithm returns Bhigh , the minimum of all eligible edge labels according toCorollary 15, together with a root edge label Broot which ensures the represented quantum stateremains the same.
Below, we will explain O(n3 )-time algorithms for finding generating sets for the stabilizer subgroupof a reduced node and for ArgLexMin. Since all other lines in Alg. 12 can be performed in linear
time, its overall runtime is O(n3 ).

D.1 Constructing the stabilizer subgroup of a LIMDD nodeIn this section, we give a recursive subroutine GetStabilizerGenSet to construct the stabilizer
subgroup Stab(|v⟩) = {A ∈ PauliLIMn | A |v⟩ = |v⟩} of an n-qubit LIMDD node v (see Sec. 2).
This subroutine is used by the algorithm GetLabels to select a canonical label for the high edge
and root edge. If the stabilizer subgroup of v’s children have been computed already, GetStabilizerGenSet’s runtime is O(n3 ). GetStabilizerGenSet returns a generating set for thegroup Stab(|v⟩). Since these stabilizer subgroups are generally exponentially large in the numberof qubits n, but they have at most n generators, storing only the generators instead of all elements
may save an exponential amount of space. Because any generator set G of size |G| > n can bebrought back to at most n generators in time O(|G| · n2 ) (see App. A), we will in the derivation
below show how to obtain generator sets of size linear in n and leave the size reduction implicit.
We will also use the notation A · G and G · A to denote the sets {A · g | g ∈ G} and {g · A | g ∈ G},
respectively, when A is a Pauli LIM.
We now sketch the derivation of the algorithm. The base case of the algorithm is the Leaf node of
the LIMDD, representing the number 1, which has stabilizer group {1}. For the recursive case, we
I

Bhigh

v1
v
. If Bhigh = 0,
wish to compute the stabilizer group of a reduced n-qubit node v = v0′
then it is straightforward to see that λPn ⊗ P |v⟩ = |v⟩ implies Pn ∈ {I, Z}, and further that
Stab(|v⟩) = ⟨{Pn ⊗ g | g ∈ G0 , Pn ∈ {I, Z}}⟩, where G0 is a stabilizer generator set for v0 .

Otherwise, if Bhigh ̸= 0, then we expand the stabilizer equation λP |v⟩ = |v⟩:λPn ⊗ P ′ (|0⟩ ⊗ |v0 ⟩ + |1⟩ ⊗ Bhigh |v1 ⟩) = |0⟩ ⊗ |v0 ⟩ + |1⟩ ⊗ Bhigh |v1 ⟩ , which implies:
λP ′ |v0 ⟩ = |v0 ⟩ and zλP ′ Bhigh |v1 ⟩ = Bhigh |v1 ⟩y ∗ λP ′ Bhigh |v1 ⟩ = |v0 ⟩and λP ′ |v0 ⟩ = y ∗ Bhigh |v1 ⟩if Pn = [ 10 z0 ] with z ∈ {1, −1}h ∗i
if Pn = y0 y0 , with y ∈ {1, i}

(21)
(22)

The stabilizers can therefore be computed according to Eq. 21 and 22 as follows.
[
Stab(|v⟩) =[ 10 z0 ] ⊗ (Stab(|v0 ⟩) ∩ z · Stab(Bhigh |v1 ⟩))z=∈{1,−1},y∈{1,i}

∪

h

0 y
y∗ 0

i


⊗ Iso(y ∗ Bhigh |v1 ⟩ , |v0 ⟩) ∩ Iso(|v0 ⟩ , y ∗ Bhigh |v1 ⟩)(23)

where Iso(v, w) denotes the set of Pauli isomorphisms A which map |v⟩ to |w⟩ and we have denoted
π·G := {π · g | g ∈ G} for a set G and a single operator π. Lemma 5 shows that such an isomorphism
set can be expressed in terms of the stabilizer group of |v⟩.
Lemma 5. Let |φ⟩ and |ψ⟩ be quantum states on the same number of qubits. Let π be a Pauli
isomorphism mapping |φ⟩ to |ψ⟩. Then the set of Pauli isomorphisms mapping |φ⟩ to |ψ⟩ is
Iso(|v⟩ , |w⟩) = π · Stab(|φ⟩). That is, the set of isomorphisms |φ⟩ → |ψ⟩ is a coset of the stabilizer
subgroup of |φ⟩.
Proof. If P ∈ Stab(|φ⟩), then π · P is an isomorphism since π · P |φ⟩ = π |φ⟩ = |ψ⟩. Conversely,if σ is a Pauli isomorphism which maps |φ⟩ to |ψ⟩, then π −1 σ ∈ Stab(|φ⟩) because π −1 σ |φ⟩ =π −1 |ψ⟩ = |φ⟩. Therefore σ = π(π −1 σ) ∈ π · Stab(|φ⟩).
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

47

With Lemma 5 we can rewrite eq. (23) as
Stab(|v⟩) =I ⊗ (Stab(|v0 ⟩) ∩ Stab(Bhigh |v1 ⟩))|{z}
stabilizer subgroup∪ Z ⊗ (I · Stab(|v0 ⟩) ∩ −I · Stab(Bhigh |v1 ⟩)){z}|
isomorphism set∪[
y∈{1,i}

h

0 y
y∗ 0

i

⊗ π · Stab(y ∗ Bhigh · |v1 ⟩) ∩ π −1 · Stab(|v0 ⟩)
|
{z
}



(24)

isomorphism set

where π denotes a single isomorphism y ∗ Bhigh |v1 ⟩ → |v0 ⟩.
Given generating sets for Stab(v0 ) and Stab(v1 ), evaluating eq. (24) requires us to:• Compute Stab(A |w⟩) from Stab(w)(as generating sets) for Pauli LIM A and node w. It isstraightforward to check that AgA† g ∈ G , with ⟨G⟩ = Stab(w), is a generating set for
Stab(A |w⟩).
• Find a single isomorphism between two edges, pointing to reduced nodes. In a reduced
LIMDD, edges represent isomorphic states if and only if they point to the same nodes. This
results in a straightforward algorithm, see Alg. 16.
• Find the intersection of two stabilizer subgroups, represented as generating sets G0 and G1(IntersectStabilizerGroups, Alg. 15). First, it is straightforward to show that the intersectionof two stabilizer subgroups is again a stabilizer subgroup (namely, it is abelian and does notcontain −I. It is never empty since I is a stabilizer of all states).‖ The algorithm proceeds in
two steps: first, we compute the intersection of G0 and G1 modulo phase; second, we “correctfor” the fact that the phases play a role.
Very broadly speaking, we use the following algebraic properties of Pauli groups. First, when
considering a Pauli string λP modulo phase, it is convenient to think of it as simply the Paulistring P with phase equal to +1. This allows us to take an element P ∈ G0 from a Paulistabilizer group G0 modulo phase, and, by abuse of language, multiply it by a phase λ ∈ C toobtain λ · P ∈ G0 . Second, for any Pauli string P ∈ G0 ∩ G1 , i.e., in the group modulo phase,
there exists a unique λ such that λ·P ∈ ⟨G0 ⟩; and a unique ω such that ω·P ∈ ⟨G1 ⟩. Moreover,
we have ω = ±λ in this case due to anti-commutativity. Consequently, if S = {P1 , . . . , Pℓ }is a generating set for the group ⟨S⟩ = G0 ∩ G1 modulo phase, then we can divide thesegenerators Pj into two sets, S = S0 ∪ S1 , where each P ∈ S0 satisfies λP ∈ G0 ∩ G1 for someλ ∈ C and each P ∈ S1 satisfies λP ∈ G0 and −λP ∈ G1 . The algorithm finds these sets
S0 and S1 , including phase, in the loop in Line 5-11. Given such sets S0 , S1 , any element in
G0 ∩ G1 (i.e., the set we are interested in) can be written as a product of elements of S0 andan even number of elements from S1 . The set {e1 · e2 , . . . , e1 · em }, found by the algorithm
in Line 14-16, generates precisely this set of elements generated from an even number of
elements of S1 .
All of the above steps can be performed in O(n3 ) time, where n is the number of qubits. In
particular, a generating set for the intersection of G0 and G1 modulo phase is simply theintersection of two vector spaces over F2 , which is constructed in time O(n3 ) on Line 3 using
the Zassenhaus algorithm. On line Line 7, checking whether λP ∈ G0 for given λ, P can bedone in O(n2 ) time; this happens at most O(|S|) = O(n) times, so in total this operation
takes up O(n3 ) time. Lastly, the loop in Line 14-16 runs in O(n2 ) time, as there are at most|S| − 1 = O(n) multiplications of Pauli strings, each of which takes O(n) time. We remarkthat the Zassenhaus algorithm cannot be straightforwardly applied to find the intersection of
the groups G0 and G1 directly, since the elements of G0 may not commute with those of G1 .
‖ To be clear, here we consider the stabilizers including their phase, i.e., we are not considering the groups modulo
phase. Indeed, computing the intersection of two groups modulo phase is relatively easy, as shown in App. A.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

48• IntersectIsomorphismSets: Find the intersection of two isomorphism sets, represented
as single isomorphism (π0 , π1 ) with a generator set of a stabilizer subgroup (G0 , G1 ), see
Lemma 5. This is the coset intersection problem for the PauliLIMn group. Isomorphism
sets are coset of stabilizer groups (see Lemma 5) and it is not hard to see that that theintersection of two cosets, given as isomorphisms π0/1 and generator sets G0/1 , is either
empty, or a coset of ⟨G0 ⟩ ∩ ⟨G1 ⟩ (this intersection is computed using Alg. 15). Therefore,
we only need to determine an isomorphism π ∈ π0 ⟨G0 ⟩ ∩ π1 ⟨G1 ⟩, or infer that no suchisomorphism exists.
We solve this problem in O(n3 ) time in two steps (see Alg. 14 for the full algorithm). First,
we note that that π0 ⟨G0 ⟩ ∩ π1 ⟨G1 ⟩ = π0 [⟨G0 ⟩ ∩ (π0−1 π1 )⟨G1 ⟩], so we only need to find anelement of the coset S := ⟨G0 ⟩ ∩ (π0−1 π1 )⟨G1 ⟩. Now note that S is nonempty if and only ifthere exists g0 ∈ ⟨G0 ⟩, g1 ∈ ⟨G1 ⟩ such that g0 = π0−1 π1 g1 , or, equivalently, π0−1 π1 ·g1 ·g0−1 = I.
We show in Lemma 6 that such g0 , g1 exist if and only if I is the smallest element in theset Sπ0−1 π1 ⟨G1 ⟩ · ⟨G0 ⟩. Hence, for finding out if S is empty we may invoke the LexMin
algorithm we have already used before in GetLabels and we will explain below in Sec. D.2.
If it is not empty, then we obtain g0 , g1 as above using ArgLexMin, and output π0 · g0 as
an element in the intersection. Since Lexmin and ArgLexMin take O(n3 ) time, so does
Alg. 14.
Lemma 6. The coset S := ⟨G0 ⟩ ∩ π1−1 π0 · ⟨G1 ⟩ is nonempty
if and only if the lexicographicallysmallest element of the set S = π0−1 π1 ⟨G1 ⟩ · ⟨G0 ⟩ = π0−1 π1 g1 g0 g0 ∈ G0 , g1 ∈ G1 is 1 · I.
Proof. (Direction ⇒) Suppose that the set ⟨G0 ⟩ ∩ π0−1 π1 ⟨G1 ⟩ has an element a. Then a = g0 =π0−1 π1 g1 for some g0 ∈ ⟨G0 ⟩ , g1 ∈ ⟨G1 ⟩. We see that I = π0−1 π1 g1 g0−1 ∈ π0−1 π1 ⟨G1 ⟩ · ⟨G0 ⟩, i.e.,
I ∈ S. Note that I is, in particular, the lexicographically smallest element, since its check vector isthe all-zero vector (⃗0|⃗0|00).
(Direction ⇐) Suppose that I ∈ π0−1 π1 ⟨G1 ⟩ · ⟨G0 ⟩. Then I = π0−1 π1 g1 g0 , for some g0 ∈ ⟨G0 ⟩ , g1 ∈⟨G1 ⟩, so we get g0−1 = π0−1 π1 g1 ∈ ⟨G0 ⟩ ∩ π0−1 π1 ⟨G1 ⟩, as promised.
The four algorithms above allow us to evaluate each of the four individual terms in eq. (24). To
finish the evaluation of eq. (24), one would expect that it is also necessary that we find the unionof isomorphism sets. However, we note that if πG is an isomorphism set, with π an isomorphism
and G an stabilizer subgroup, then Pn ⊗ (πg) = (Pn ⊗ π)(I ⊗ g) for all g ∈ G. Therefore, wewill evaluate eq. (24), i.e. find (a generating set) for all stabilizers of node v in two steps. First,
we construct the generating set for the first term, i.e. I ⊗ (Stab(|v0 ⟩) ∩ Stab(Bhigh |v1 ⟩)), using
the algorithms above. Next, for each of the other three terms Pn ⊗ (πG), we add only a single
stabilizer of the form Pn ⊗ π for each Pn ∈ {X, Y, Z}. We give the full algorithm in Alg. 13 andprove its efficiency below.
Lemma 7 (Efficiency of function GetStabilizerGenSet). Let v be an n-qubit node. Assume
that generator sets for the stabilizer subgroups of the children v0 , v1 are known, e.g., by an earlier
call to GetStabilizerGenSet, followed by caching the result (see Line 28 in Alg. 13). Then
Alg. 13 (function GetStabilizerGenSet), applied to v, runs in time O(n3 ).
Proof. If n = 1 then Alg. 13 only evaluates Line 2–4, which run in constant time. For n > 1,
the algorithm performs a constant number of calls to GetIsomorphism (which only multiplies
two Pauli LIMs and therefore runs in time O(n)) and four calls to IntersectIsomorphismSets.
Note that the function IntersectIsomorphismSets from Alg. 14 invokes O(n3 )-runtime external
algorithms:• the Zassenhaus algorithm [80] to calculate a basis for the intersection of two subspaces of avector space,• the RREF algorithms mentioned in App. A, andAccepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

49

• Algorithm 2 from [82] to synthesize a circuit that transforms any stabilizer state to a basis
state. Specifically, this algorithm receives as input a stabilizer subgroup G and outputs a
Clifford circuit U such that U GU † = {Z1 , . . . , Z|G| }. We remark that Garcı́a et al. assume
in their work that G is the stabilizer group of a stabilizer state, i.e., |G| = n, but in fact thealgorithm works also without that assumption, i.e., in the more general case when G is any
abelian group of Pauli operators not containing −I. Our algorithms use this more generaluse case.
Therefore, GetStabilizerGenSet has runtime is O(n3 ).
Algorithm 13 Algorithm for constructing the Pauli stabilizer subgroup of a Pauli-LIMDD node.
The algorithm always returns a set in reduced row echelon form (see App. A), which is accomplishedin line 27. In particular, the set always returns at most n elements for n-qubit states.
1: procedure GetStabilizerGenSet(Edge e0
2:
3:
4:
5:

I⊗n

v 0 , e1

Bhigh

v1

with v0 , v1 reduced)

if n=1 thenif there exists P ∈ ±1 · {X, Y, Z} such that P |v⟩ = |v⟩ then return P
else return None
else

6:

if v ∈ StabCache[v] then return StabCache[v]7:G0 := GetStabilizerGenSet(v0 )8:

if Bhigh = 0 then

9:
10:
11:

return {I2 ⊗ g | g ∈ G0 } ∪ {Z ⊗ I⊗n−1 }else
G := ∅no†G1 := Bhigh · g · Bhighg ∈ GetStabilizerGenSet(v1 )

12:
13:
14:
15:
16:

(π, B) := IntersectIsomorphismSets((I⊗n−1 , G0 ), (I⊗n−1 , G1 ))G := G ∪ {I2 ⊗ g | g ∈ B}
▷ Add all stabilizers of the form I ⊗ . . .

17:

π0 , π1 := I⊗n−1 , GetIsomorphism(e1 , −1 · e1 )(π, B) := IntersectIsomorphismSets((π0 , G0 ), (π1 , G1 ))18:

if π ̸= None then G := G ∪ {Z ⊗ π}

▷ Add stabilizer of form Z ⊗ . . .

19:
20:
21:

π0 , π1 := GetIsomorphism(e0 , e1 ), GetIsomorphism(e1 , e0 ))(π, B) := IntersectIsomorphismSets((π0 , G0 ), (π1 , G1 ))22:

if π ̸= None then G := G ∪ {X ⊗ π}

23:
24:
25:

▷ Add stabilizer of form X ⊗ . . .

π0 , π1 := GetIsomorphism(e0 , −i · e1 ), GetIsomorphism(−i · e1 , e0 ))(π, B) := IntersectIsomorphismSets((π0 , G0 ), (π1 , G1 ))26:

if π ̸= None then G := G ∪ {Y ⊗ π}27:

G := RREF (G)28:

StabCache[v] := G29:

return G▷ Add stabilizer of form Y ⊗ . . .

▷ Bring G to reduced row echelon form, potentially pruning some rows

D.2 Efficiently finding a minimal LIM by multiplying with stabilizers
Here, we give O(n3 ) subroutines solving the following problem: given generators sets G0 , G1 ofstabilizer subgroups on n qubits, and an n-qubit Pauli LIM A, determine min(g0 ,g1 )∈⟨G0 ,G1 ⟩ A·g0 ·g1 ,and also find the g0 , g1 which minimize the expression. We give an algorithm for finding both the
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

50

Algorithm 14 An O(n3 ) algorithm for computing the intersection of two sets of isomorphisms,
each given as single isomorphism with a stabilizer subgroup (see Lemma 5).
output: Pauli LIM π, stabilizer subgroup generating set G s.t. π⟨G⟩ = π0 ⟨G0 ⟩ ∩ π1 ⟨G1 ⟩1: procedure IntersectIsomorphismSets(stabilizer subgroup generating sets G0 , G1 ,Pauli-LIMs π0 , π1 )2:
π := LexM in(G0 , G1 , π1−1 π0 )
3:
4:
5:
6:

if π = I then(g0 , g1 ) = ArgLexM in(G0 , G1 , π1−1 π0 )
π := π0 · g0
G := IntersectStabilizerGroups(G0 , G1 )7:
8:

return (π, G)
else

9:

return None

Algorithm 15 Algorithm for finding the intersection ⟨G0 ⟩ ∩ ⟨G1 ⟩ of the groups generated by two
stabilizer subgroup generating sets G0 and G1 .
output: a generating set for ⟨G0 ⟩ ∩ ⟨G1 ⟩1: procedure IntersectStabilizerGroups(stabilizer subgroup generating sets G0 , G1 )2:
Use the Zassenhaus algorithm to compute the intersection modulo phase3:
S := IntersectGroupsModuloPhase(G0 , G1 )4:
J, S0 , S1 := ∅
5:
6:
7:
8:
9:
10:
11:

for P ∈ S do
By abuse of language, we treat P as a Pauli string with phase +1
Find λ ∈ {±1, ±i} such that λ · P ∈ ⟨G0 ⟩if λ · P ∈ ⟨G1 ⟩ thenS0 := S0 ∪ {λ · P }else if −λ · P ∈ ⟨G1 ⟩ thenS1 := S1 ∪ {λ · P }12:J := S013:

if ∃e ∈ S1 then

14:
15:
16:
17:

for e′ ∈ S1 \ {e} do
q := e′ · e
J := J ∪ {q}return J

Algorithm 16 Algorithm for constructing a single isomorphism between the quantum states
represented by two Pauli-LIMDD edges, each pointing to a canonical node.
1: procedure GetIsomorphism(Edge
2:
3:
4:

A

v ,

B

w

with v, w reduced, A ̸= 0 ∨ B ̸= 0)

if v = w and A, B ̸= 0 thenreturn B · A−1
return None

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

51

minimum (LexMin) and the arguments of the minimum (ArgLexMin) in Alg. 17. The intuition
behind the algorithms are the following two steps: first, the lexicographically minimum Pauli LIM
modulo scalar can easily be determined using the scalar-ignoring DivisionRemainder algorithm
from App. A. Since in the lexicographic ordering, the scalar is least significant (App. A), theresulting Pauli LIM has the same Pauli string as the the minimal Pauli LIM including scalar. We
show below in Lemma 8 that if the scalar-ignoring minimization results in a Pauli LIM λP , thenthe only other eligible LIM, if it exists, is −λP . Hence, in the next step, we only need to determine
whether such LIM −λP exists and whether −λ < λ; if so, then −λP is the real minimal Pauli LIM∈ ⟨G0 ∪ G1 ⟩.
Lemma 8. Let v0 and v1 be LIMDD nodes, R a Pauli string and ν, ν ′ ∈ C. Define G = Stab(v0 ) ∪Stab(v1 ). If νR, ν ′ R ∈ ⟨G⟩, then ν = ±ν ′ .
Proof. We prove g ∈ ⟨G⟩ and λg ∈ ⟨G⟩ implies λ = ±1. Since Pauli LIMs commute or anticommute, we can decompose both g and λg as g = (−1)x g0 g1 and λg = (−1)y h0 h1 for some x, y ∈ {0, 1}
and g0 , h0 ∈ Stab(v0 ) and g1 , h1 ∈ Stab(v1 ). Combining these yields λ(−1)x g0 g1 = (−1)y h0 h1 .
We recall that, if g ∈ Stab(|φ⟩) is a stabilizer of any state, then g 2 = I. Therefore, squaring bothsides of the equation, we get λ2 (g0 g1 )2 = (h0 h1 )2 , so λ2 I = I, so λ = ±1.
The central procedure in Alg. 17 is ArgLexMin, which, given a LIM A and sets G0 , G1 which
generate stabilizer groups, finds g0 ∈ ⟨G0 ⟩ , g1 ∈ ⟨G1 ⟩ such that A · g0 · g1 reaches its lexicographic
minimum over all choices of g0 , g1 . It first performs the scalar-ignoring minimization (Line 5) to find
g0 , g1 modulo scalar. The algorithm LexMin simply invokes ArgLexMin to get the arguments
g0 , g1 which yield the minimum and uses these to compute the actual minimum.
The subroutine FindOpposite finds an element g ∈ G0 such that −g ∈ G1 , or infers that nosuch g exists. It does so in a similar fashion as IntersectStabilizerGroups from Sec. D.1:
by conjugation with a suitably chosen unitary U , it maps G1 to {Z1 , Z2 , . . . , Z|G1 | }. Analogously
to our explanation of IntersectStabilizerGroups, the group generated by U G1 U † contains
precisely all Pauli LIMs which satisfy the following three properties: (i) the scalar is 1; (ii) its
Pauli string has an I or Z at positions 1, 2, . . . , |G1 |; (iii) its Pauli string has an I at positions|G1 | + 1, . . . , n. Therefore, the target g only exists if there is a LIM in ⟨U G0 U † ⟩ which (i’) has
scalar −1 and satisfies properties (ii) and (iii). To find such a g, we put U G0 U † in RREF formand check all resulting generators for properties (i’), (ii) and (iii). (By definition of RREF, it
suffices to check only the generators for this property) If a generator h satisfies these properties,
we return U † hU and None otherwise. The algorithm requires O(n3 ) time to find U , the conversion
G 7→ U GU † can be done in time O(n3 ), and O(n) time is required for checking each of the O(n2 )generators. Hence the runtime of the overall algorithm is O(n3 ).

E Measuring an arbitrary qubit
Alg. 18 allows one to measure a given qubit. Specifically, given a quantum state |e⟩ represented by a
LIMDD edge e, a qubit index k and an outcome b ∈ {0, 1}, it computes the probability of observing|b⟩ when measuring the k-th significant qubit of |e⟩. The algorithm proceeds by traversing the
LIMDD with root edge e at Line 7. Like Alg. 5, which measured the top qubit, this algorithm finds
the probability of a given outcome by computing the squared norm of the state when the k-th qubit
is projected onto |0⟩, or |1⟩. The case that is added, relative to Alg. 5, is the case when n > k,in which case it calls the procedure SquaredNormProjected. On input e, y, k, the procedure
SquaredNormProjected outputs the squared norm of Πyk |e⟩, where Πyk = In−k ⊗ |y⟩ ⟨y| ⊗ Ik−1
is the projector which projects the k-th qubit onto |y⟩.
After measurement of a qubit k, a quantum state is typically projected to |0⟩ or |1⟩ (b = 0 or
b = 1) on that qubit, depending on the outcome. Alg. 19 realizes this. It does so by traversing the
LIMDD until a node v with idx(v) = k is reached. It then returns an edge to a new node by calling
MakeEdge(follow0 (e), 0) to project onto |0⟩ or MakeEdge(0, follow1 (e)) to project onto |1⟩,
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

52

Algorithm 17 Algorithms LexMin and ArgLexMin for computing the minimal element fromthe set A · ⟨G0 ⟩ · ⟨G1 ⟩ = {Ag0 g1 | g0 ∈ G0 , g1 ∈ G1 }, where A is a Pauli LIM and G0 , G1 are
generating sets for stabilizer subgroups. The algorithms make use of a subroutine FindOpposite
for finding an element g ∈ ⟨G0 ⟩ such that −g ∈ ⟨G1 ⟩. A canonical choice for the Rootlabel (seeSec. 3.3.3) of an edge e pointing to a node v is LexMin(G, {I}, label(e)) where G is a stabilizer
generator group of Stab(v).
output: min(g0 ,g1 )∈⟨G0 ∪G1 ⟩ A · g0 · g1
1: procedure LexMin(stabilizer subgroup generating sets G0 , G1 and Pauli LIM A)
2:(g0 , g1 ) := ArgLexMin(G0 , G1 , A)
return A · g0 · g1
3:
output: arg ming0 ∈G0 ,g1 ∈G1 A · g0 · g1

4: procedure ArgLexMin(stabilizer subgroup generating sets G0 , G1 and Pauli LIM A)
5:(g0 , g1 ) :=arg min{h | h ∝ A · g0 · g1 }(g0 ,g1 )∈⟨G0 ∪G1 ⟩6:▷ Using the scalar-ignoring DivisionRemainder algorithm from App. A,
g ′ := FindOpposite(G0 , G1 , g0 , g1 )7:

if g ′ is None then

8:
9:
10:

return (g0 , g1 )else
h0 , h1 := g0 · g ′ , (−g ′ ) · g1

▷ g0 g1 = −h0 h1

11:

if A · h0 · h1 <lex A · g0 · g1 then return (h0 , h1 )12:else return (g0 , g1 )output: g ∈ G0 such that −g ∈ G1 , or None if no such g exists

13: procedure FindOpposite(stabilizer subgroup generating sets G0 , G1 )14:
Compute U s.t. U G1 U † = {Z1 , Z2 , . . . , Z|G1 | }, using Algorithm 2 from [82]▷ Zj is the Z gate applied to qubit with index j
15:
16:

H0 := U G0 U †H0RREF := H0 in RREF form

17:
18:

for h ∈ H0RREF do
if h satisfies all three of the following: (i) h has scalar −1; the Pauli string of h (ii)contains only I or Z at positions 1, 2, . . . , |G1 |, and (iii) only I at positions |G1 | + 1, . . . , n then

19:
20:

return U † hU
return None

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

53

Algorithm 18 Compute the probability of observing |y⟩ when measuring the k-th qubit of the state|e⟩. Here e is given as LIMDD on n qubits, y is given as a bit, and k is an integer index. For example,
to measure the top-most qubit, one calls Measure(e, 0, n). The procedure SquaredNorm(e, y, k)
computes the scalar ⟨e| (I ⊗ |y⟩ ⟨y| ⊗ I) |e⟩, i.e., computes the squared norm of the state |e⟩ after
the k-th qubit is projected to |y⟩. For readability, we omit calls to the cache, which implement
dynamic programming.
1: procedure MeasurementProbability(Edge e

λPn ⊗ P ′v , y ∈ {0, 1}, k ∈ [1 . . . idx(v)])

3:
4:

if n = k thenp0 := SquaredNorm(follow0 (e))
p1 := SquaredNorm(follow1 (e))

5:
6:
7:
8:

return pj /(p0 + p1 ) where j = 0 if Pn ∈ {I, Z} and j = 1 if Pn ∈ {X, Y }else
p0 := SquaredNormProjected(follow0 (e), y, k)
p1 := SquaredNormProjected(follow1 (e), y, k)2:

9:return (p0 + p1 )/SquaredNorm(e)10: procedure SquaredNorm(Edge

λP

v )

2

11:

if n = 0 then return |λ|12:

s := Add(SquaredNorm(follow0 (13:Iv )), SquaredNorm(follow1 (

16:
17:
18:
19:
20:
21:
22:
23:

v )))

2

return |λ| s

14: procedure SquaredNormProjected(Edge e
15:

I

λPn ⊗ P ′b := (Pn ∈ {X, Y })if n = 0 thenv , y ∈ {0, 1}, k ∈ [1 . . . idx(v)])▷ i.e., b = 1 iff Pn is Anti-diagonal

return |λ|2else if n = k thenreturn SquaredNorm(followb⊕y (e))
elseIα0 := SquaredNormProjected(follow0 (I
α1 := SquaredNormProjected(follow1 (v ), b ⊕ y, k)v ), b ⊕ y, k)return |λ|2 · (α0 + α1 )on Line 6, recreating a node on level k in the backtrack on Line 8. The projection operator Πbk
commutes with any LIM P when Pk is a diagonal operator (i.e., Pk ∈ {I2 , Z}). Otherwise, if Pk(1−b)
is an antidiagonal operator (i.e, Pk ∈ {X, Y }), have Πbk · P = P Πk
. The algorithm applies thiscorrection on Line 2. The resulting state should still be normalized as shown in Sec. 3.3.1.

F LIMDDs prepare the W state efficiently
In this section, we prove Th. 8. To this end, we show that LIMDDs can efficiently simulate acircuit family given by McClung [56], which prepares the |W ⟩ state when initialized to the |0⟩state. We thereby show a separation between LIMDD and the Clifford+T simulator, as explainedin Sec. 3.4.3. Figure Fig. 16 shows the circuit for the case of 8 qubits.
⊗n

Theorem 8. There exists a circuit family Cn such that Cn |0⟩ = |Wn ⟩, that Pauli-LIMDDs canefficiently simulate. Here simulation means that it constructs representations of all intermediate
states, in a way which allows one to, e.g., efficiently simulate any single-qubit computational-basis
measurement or compute any computational basis amplitude on any intermediate state and theoutput state.
Proof. The proof outline is as follows. First, we establish that the LIMDD of each intermediate state
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

54

Av to state |b⟩ for qubit k, i.e., produce a
Algorithm 19 Project the state given by LIMDD
LIMDD representing the state (In−k ⊗ |b⟩ ⟨b| ⊗ Ik−1 ) · A |v⟩, with A = λPn ⊗ · · · ⊗ P1 .
λPn ⊗ · · · ⊗ P1

v , k ∈ [1 . . . idx(v)], b ∈ {0, 1})
1: procedure UpdatePostMeas(Edgeb′ := x ⊕ b where x = 0 if Pk ∈ {I, Z} and x = 1 if Pk ∈ {X, Y } ▷ flip b if Pk is anti-diagonal
2:

3:
4:

if (v, k, b′ ) ∈ Cache then return Cache[v, k, b′ ]n := idx(v)

5:
6:
7:
8:

if n = k thene := MakeEdge((1 − b′ ) · lowv , b′ · highv )▷ Project |v⟩ to |b′ ⟩ ⟨b′ | ⊗ I⊗n−12else▷ n ̸= k:e := MakeEdge(UpdatePostMeas(lowv , k, b′ ), UpdatePostMeas(highv , k, b′ ))9:Cache[v, k, b′ ] := e10:

return e(Lemma 11), as well as of each gate (Lemma 12), has polynomial size. Second, we establish that thealgorithms presented in Sec. 3.3 can apply each gate to the intermediate state in polynomial time
(Lemma 16). To this end, we observe that the circuit only produces relatively simple
Pn intermediate
states. Specifically, each intermediate and output state is of the form |ψt ⟩ = √1n k=1 |xk ⟩ where
the xk ∈ {0, 1}n are computational basis vectors (Lemma 10). For example, the output state hask−1⊗n−k|xk ⟩ = |0⟩|1⟩ |0⟩. The main technical tool we will use to reason about the size of the
LIMDDs of these intermediate states, are the subfunction rank and computational basis rank ofa state. Both these measures are upper bounds of the size of a LIMDD (in Lemma 9), and also
allow us to upper bound the time taken by the ApplyGate and Add algorithms (in Lemma 13
for ApplyGate and Lemma 14 Add).
The theorem follows from Lemma 16 and Corollary 16.
Fig. 16 shows the circuit for the case of n = 8 qubits. For convenience and without loss of generality,
we only treat the case when the number of qubits is a power of 2, since the circuit is simplest inthat case. In general, the circuit works as follows. The qubits are divided into two registers;
register A, with log n qubits, and register B, with the remaining n − log n qubits. First, thecircuit applies a Hadamard gate to each qubit in register A, to bring the state to the superposition⊗ log nn−log n|+⟩|0⟩
. Then it applies n − log n Controlled-X gates, where, in each gate, each qubit

Figure 16: Reproduced from McClung [56]. A circuit on eight qubits (n = 8) which takes as input the |0⟩⊗8
state and outputs the |W8 ⟩ state. In the general case, it contains log n Hadamard gates, and its Controlled-X
gates act on one target qubit and at most log n control qubits.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

55

of register A acts as the control qubits and one qubit in register B is the target qubit. Lastly, itapplies n − log n Controlled-X gates, where, in each gate, one qubit in register B is the control
qubit and one or more qubits in register A are the target qubits. Each of the three groups of gates⊗nis highlighted in a dashed rectangle in Fig. 16. On input |0⟩ , the circuit’s final state is |Wn ⟩. We
emphasize that the Controlled-X gates are permutation gates (i.e., their matrices are permutation
matrices). Therefore, these gates do not influence the number of non-zero computational basis
state amplitudes of the intermediate states. We refer to the t-th gate of this circuit as Ut , and thet-th intermediate state as |ψt ⟩, so that |ψt+1 ⟩ = Ut |ψt ⟩ and |ψ0 ⟩ = |0⟩ is the initial state.
We refer to the number of computational basis states with nonzero amplitude as a state’s computational basis rank, denoted χcomp (|ψ⟩).
P
Definition 8. (Computational basis rank) Let |ψ⟩ =x∈{0,1}n α(x) |x⟩ be a quantum state
defined by the amplitude function α : {0, 1}n → C. Then the computational basis rank of |ψ⟩ is
χcomp (|ψ⟩) = | {x | α(x) ̸= 0} |, the number of nonzero computational basis amplitudes.
Recall that, for a given function α : {0, 1}n → C, a string a ∈ {0, 1}ℓ induces a subfunction
αy : {0, 1}n−ℓ → C, defined as αy (x) = α(y, x). We refer to the number of subfunctions of a state’s
amplitude function as its subfunction rank. The following definition makes this more precise.
P
Definition 9. (Subfunction rank) Let |ψ⟩ = x∈{0,1}n αψ (x) |x⟩ be a quantum state defined by
the amplitude function αψ : {0, 1}n → C, as above. Let χsub (|ψ⟩ , ℓ) be the number of unique
non-zero subfunctions induced by strings of length ℓ, as follows,χsub (|ψ⟩ , ℓ) = | αyψ : {0, 1}n−ℓ → C αy ̸= 0, y ∈ {0, 1}ℓ |(25)
We define the subfunction rank of |ψ⟩ as χsub (|ψ⟩) = max
P ℓ=0,...n χsub (|ψ⟩ , ℓ). We extend thesedefinitions in the natural way for an n-qubit matrix U = r,c∈{0,1}n αU (r, c) |r⟩ ⟨c| defined by thefunction αU : {0, 1}2n → C.
It is easy to check that χsub (|ψ⟩) ≤ χcomp (|ψ⟩) holds for any state.
For the next lemma, we use the notion of a prefix of a LIMDD node. This lemma will serve as atool which allows us to show that a LIMDD is small when its computational basis rank is low. We
apply this tool to the intermediate states of the circuit in Lemma 11.
Definition 10 (Prefix of a LIMDD node). For a given string x ∈ {0, 1}ℓ , consider the path traversedRr ) subroutine, which starts at the diagram’s root edge and ends at a nodeby the followx (v on level ℓ. We will say that x is a prefix of the node v. We let Labels(x) be the product of the
LIMs on the edges of this path (i.e., including the root edge). The set of prefixes of a node v isdenoted pre(v).

Lemma 9. If a LIMDD represents the state |φ⟩, then its width at any given level (i.e., the numberof nodes at that level) is at most χcomp (|φ⟩).
Proof. For notational convenience, let us number the levels so that the root node is on level 0, its
children are on level 1, and so on, with the Leaf on level n (contrary to Fig. 3). Let r be the root
node of the LIMDD, and R the root edge’s label. By construction of a LIMDD, the state represented
by the LIMDD can be expressed as follows, for any level ℓ ≥ 0,
R |r⟩ =X|x⟩ ⊗ followx (

R

r )

(26)

x∈{0,1}ℓ

Since

R

r

is the root of our diagram, if x is a prefix of v, thenfollowx (Rr ) = Labels(x) · |v⟩

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

(27)

56

A string x ∈ {0, 1}ℓ can be a prefix of only one node; consequently, the prefix sets of two nodes on
the same level are disjoint, i.e., pre(vp )∩pre(vq ) = ∅ for p ̸= q. Moreover, each string x is a prefix ofRr ) subroutine arrives).
some node on level ℓ (namely, simply the node at which the followx (Say that the ℓ-th level contains m nodes, v1 , . . . , vm . Therefore, the sets pre(v1 ), . . . , pre(vm )partition the set {0, 1}ℓ . Therefore, by putting Eq. 27 and Eq. 26 together, we can express theroot node’s state in terms of the nodes v1 , . . . , vm on level ℓ:

R |r⟩ =m
X

X

|x⟩ ⊗ followx (

R

r )

(28)

k=1 x∈pre(vk )

=

m
X

X

|x⟩ ⊗ Labels(x) · |vk ⟩(29)

k=1 x∈pre(vk )P
We now show that each term x∈pre(vk ) |x⟩ ⊗ Labels(x) · |vk ⟩ contributes a non-zero vector. It
then follows that the state has computational basis rank at least m, since these terms are vectors
with pairwise disjoint support, since the sets pre(vk ) are pairwise disjoint. Specifically, we showthat each node has at least one prefix x such that Labels(x) · |v⟩ is not the all-zero vector. In
principle, this can fail in one of three ways: either v has no prefixes, or all prefixes x ∈ pre(vk )have Labels(x) = 0 because the path contains an edge labeled with the 0 LIM, or the node v
represents the all-zero vector (i.e., |v⟩ = ⃗0). First, we note that each node has at least one prefix,
since each node is reachable from the root, as a LIMDD is a connected graph. Second, due to thezero edges rule (see Def. 5), for any node, at least one of its prefixes has only non-zero LIMs on theedges. Namely, each node v has at least one incoming edge labeled with a non-zero LIM, since, if
it has an incoming edge from node w labeled with 0, then this must be the high edge of w and by
the zero edges rule the low edge of w must also point to v and moreover must be labeled with I
by the low factoring rule. Together, via a simple inductive argument, there must be at least one
non-zero path from v to the root. Lastly, no node represents the all-zero vector, due to the low
factoring rule (in Def. 5). Namely, if v is a node, then by the low factoring rule, the low edge hasA

v1 , then the node v
label I. Therefore, if this edge points to node v0 , and the high edge isrepresents |v⟩ = |0⟩ |v0 ⟩ + |1⟩ A |v1 ⟩ with possibly A = 0, so, if |v0 ⟩ ≠ ⃗0, then |v⟩ ≠ ⃗0. An argument
by induction now shows that no node in the reduced LIMDD represents the all-zero vector.
Rr ) ̸= ⃗Therefore, each node has at least one prefix x such that followx (0. We conclude thatthe equation above contains at least m non-zero contributions. Hence m ≤ χcomp (R |r⟩), at any
level 0 ≤ ℓ ≤ n.

Lemma 10. Each intermediate state in the circuit in Fig. 16 (with n = 2c ) has χcomp (|ψ⟩) ≤ n.
⊗n

Proof. The initial state is |ψ0 ⟩ = |0⟩ , which is a computational basis state, so χcomp (ψ0 ) = 1.
The first log n gates are Hadamard gates, which produce the state⊗ log n|ψlog n ⟩ = H ⊗ log n ⊗ In−log n |0⟩ = |+⟩

n−1

⊗ |0⟩

⊗n−log n

1 X⊗n−log n=√|x⟩ |0⟩
n x=0

(30)

This is a superposition of n computational basis states, so we have χcomp (|ψlog n ⟩) = n. All
subsequent gates are controlled-X gates; these gates permute the computational basis states, butthey do not increase their number.
Lemma 11. The reduced LIMDD of each intermediate state in the circuit in Fig. 16 has polynomial
size.
Proof. By Lemma 9, the width of a LIMDD representing |φ⟩ is at most χcomp (|φ⟩) at any level.
Since there are n levels, the total size is at most nχcomp (|φ⟩). By Lemma 10, the intermediate
states in question have polynomial χcomp , so the result follows.
Lemma 12. The LIMDD of each gate in the circuit in Fig. 16 (with n = 2c ) has polynomial size.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

57

Proof. Each gate acts on at most k = log n + 1 qubits. Therefore, the width of any level of the
LIMDD is at most 4k = 4n2 . The height of the LIMDD is n by definition, so the LIMDD has at
most 4n3 nodes.
The ApplyGate procedure handles the Hadamard gates efficiently, since they apply a single-qubit
gate to a product state. The difficult part is to show that the same holds for the controlled-X
gates. To this end, we show a general result for the speed of LIMDD operations (Lemma 13).
Although this worst-case upper bound is tight, it is exponentially far removed from the best case,e.g., in the case of Clifford circuits, in which case the intermediate states can have exponentialχsub , yet the LIMDD simulation is polynomial-time, as shown in Sec. 3.3.4.
Lemma 13. The number of recursive calls made by subroutine ApplyGate(U, |ψ⟩), is at mostnχsub (U )χsub (|ψ⟩), for any gate U and any state |ψ⟩.
Proof. Inspecting Alg. 8, we see that every call to ApplyGate(U, |ψ⟩) produces four new recursive calls, namely ApplyGate(followrc (U ), followc (|ψ⟩)) for r, c ∈ {0, 1}. Therefore,
the set of parameters in all recursive calls of ApplyGate(U, |ψ⟩) is precisely the set of tuples(followrc (U ), followc (|ψ⟩)), with r, c ∈ {0, 1}ℓ with ℓ = 0 . . . n. The terms followrc (U ) andfollowc (|ψ⟩) are precisely the subfunctions of U and |ψ⟩, and since there are at most χsub (U )and χsub (|ψ⟩) of these, the total number of distinct parameters passed to ApplyGate in recursive
calls at level ℓ, is at most χsub (U, ℓ) · χsub (|ψ⟩ , ℓ) ≤ χsub (U ) · χsub (|ψ⟩). Summing over the nlevels of the diagram, we see that there are at most nχsub (U )χsub (|ψ⟩) distinct recursive calls in
total. As detailed in Sec. 3.3.3, the ApplyGate algorithm caches its inputs in such a way that
it will achieve a cache hit on a call ApplyGate(U ′ , |ψ ′ ⟩) when it has previously been called withparameters U, |ψ⟩ such that U = U ′ and |ψ⟩ = |ψ ′ ⟩. Therefore, the total number of recursive callsthat is made, is equal to the number of distinct calls, and the result follows.
In our case, both χsub (U ) and χsub (|ψ⟩) are polynomial, so a polynomial number of recursive calls
to ApplyGate is made. We now show that also the Add subroutine makes only a small numberof recursive calls every time it is called from ApplyGate. First, Lemma 14 shows expresses aworst-case upper bound on the number of recursive calls to Add in terms of χsub . Then Lemma 15
uses this result to show that, in our circuit, the number of recursive calls is polynomial in n.
Lemma 14. The number of recursive calls made by the subroutine Add(|α⟩ , |β⟩) is at mostnχsub (|α⟩) · χsub (|β⟩), if |α⟩ , |β⟩ are n-qubit states.
Proof. Inspecting Alg. 9, every call to Add(|α⟩ , |β⟩) produces two new recursive calls, namely
Add(follow0 (|α⟩), follow0 (|β⟩)) and Add(follow1 (|α⟩), follow1 (|β⟩)). It follows that theset of parameters on n−ℓ qubits with which Add is called is the set of tuples (followx (|α⟩), followx (|β⟩)),for x ∈ {0, 1}ℓ . This corresponds precisely to the set of subfunctions of α and β induced by length-ℓstrings, of which there are χsub (|α⟩ , ℓ) and χsub (|β⟩ , ℓ), respectively. Because the results of previous computations are cached, as explained in Sec. 3.3.3, the total number of recursive calls is
the number of distinct recursive calls. Therefore, we get the upper bound of χsub (|α⟩) · χsub (|β⟩)
for each level of the LIMDD. Since the LIMDD has n levels, the upper bound nχsub (|α⟩) · χsub (|β⟩)
follows.
Lemma 15. The calls to Add(|α⟩ , |β⟩) that are made by the recursive calls to ApplyGate(Ut , |ψt ⟩),satisfy χsub (|α⟩), χsub (|β⟩) = poly(n).
Proof. We have established that the recursive calls to ApplyGate are all called with parameters of
the form ApplyGate(followr,c (Ut ), followc (|ψt ⟩)) for some r, c ∈ {0, 1}ℓ . Inspecting Alg. 8, we
see that, within such a call, each call to Add(|α⟩ , |β⟩) has parameters which are both of the form|α⟩ , |β⟩ = ApplyGate(followrx,cy (Ut ), followcy (|ψt ⟩)) for some x, y ∈ {0, 1}; therefore, theparameters |α⟩ , |β⟩ are of the form |α⟩ , |β⟩ = followr,c (Ut )·followr (|ψt ⟩). Here followcy (|ψt ⟩)is a quantum state on n − (ℓ + 1) qubits.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

58

The computational basis rank of a state is clearly non-increasing under taking subfunctions; thatis, for any string x, it holds that, χcomp (followx (|ψ⟩)) ≤ χcomp (|ψ⟩). In particular, we haveχcomp (followcy (|ψt ⟩)) ≤ χcomp (|ψt ⟩) = O(n). The matrix followrx,cy (Ut ) is a subfunction of a
permutation gate, and applying such a matrix to a vector cannot increase its computational basis
rank, so we haveχsub (|α⟩) =χsub (followrx,cy (Ut ) · followcy (|ψt ⟩))(31)≤χcomp (followrx,cy (Ut ) · followcy (|ψt ⟩)) ≤ χcomp (followcy (|ψt ⟩))(32)≤χcomp (|ψt ⟩) = O(n)(33)

This proves the lemma.
Lemma 16. Each call to ApplyGate(Ut , |ψt ⟩) runs in polynomial time, for any gate Ut in thecircuit in Fig. 16 (with n = 2c ).
Proof. If Ut is a Hadamard gate, then LIMDDs can apply this in polynomial time by Th. 6, since|ψt ⟩ is a stabilizer state. Otherwise, Ut is one of the controlled-X gates. In this case there area polynomial number of recursive calls to ApplyGate, by Lemma 13. Each recursive call to
ApplyGate makes two calls to Add(|α⟩ , |β⟩), where both α and β are states with polynomial
subfunction rank, by Lemma 15. By Lemma 14, these calls to Add all complete in time polynomialin the subfunction rank of its arguments.
Corollary 16. The circuit in Fig. 16 (with n = 2c ) can be simulated by LIMDDs in polynomial
time.

G Numerical search for the stabilizer rank of Dicke states
Given the separation between the Clifford + T simulator —a specific stabilizer-rank based simulator—and Pauli-LIMDDs, it would be highly interesting to theoretically compare Pauli-LIMDDs and general stabilizer-rank simulation. However, proving an exponential separation would require us tofind a family of states for which we can prove its stabilizer rank scales exponentially, which is amajor open problem. Instead, we here take the first steps towards a numerical comparison bychoosing a family of circuits which Pauli-LIMDDs can efficiently simulate and using Bravyi et al.’sheuristic algorithm for searching the stabilizer rank of the circuits’ output states [12]. If the stabilizer rank is very high (specifically, if it grows superpolynomially in the number of qubits), thenwe have achieved the goal of showing a separation. We cannot use W states for showing thisseparation because the n-qubit W state |Wn ⟩ has linear stabilizer rank, since it is a superpositionnof only n computational basis states. Instead we focus on their generalization, Dicke states |Dw⟩,
which are equal superpositions of all n-qubit computational-basis status with Hamming weight w(note |Wn ⟩ = |D1n ⟩),
n
|Dw
⟩= q

1

X

n

w

|x⟩

(34)

x:|x|=w

We implemented the algorithm by Bravyi et al.: see [83] for our open-source implementation.
Unfortunately, the algorithm’s runtime grows significantly in practice, which we believe is due tothe fact that it acts on sets of quantum state vectors, which are exponentially large in the number of
qubits. Our implementation allowed us to go to at most 9 qubits using the SURF supercomputing
cluster. We believe this is a limitation of the algorithm and not of our implementation, since
Bravyi et al. do not report beyond 6 qubits while Calpin uses the same algorithm and reaches atmost 10 qubits [84]. Table 2 shows the heuristically found stabilizer ranks of Dicke states with our
implementation. Although we observe the maximum found rank over w to grow quickly in n, thefeasible regime (i.e. up to 9 qubits) is too small to draw a firm conclusion on the stabilizer ranks’scaling.
Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

59

nTable 2: Heuristically-found upper bounds on the stabilizer rank χ of Dicke states |Dw⟩ (eq. (34)) using theheuristic algorithm from Bravyi et al. [12], see text in App. G for details. We investigated up to 9 qubits using
the SURF supercomputing cluster (approximately the same as the number of qubits reached in the literature as
described in the text). Empty cells indicate non-existing or not-investigated states. In particular, we have notnnnn
investigated w > ⌊ n2 ⌋ since |Dw⟩ and |Dn−w⟩ have identical stabilizer rank because X ⊗n |Dw⟩ = |Dn−w⟩. For|D38 ⟩ and |D49 ⟩, we have run the heuristic algorithm to find sets of stabilizers up to size 11 (theoretical upperbound) and 10, respectively, but the algorithm has not found sets in which these two Dicke states could be
decomposed. We emphasize that the algorithm is heuristic, so even if there exists a stabilizer decomposition ofa given rank, the algorithm might not find it.

#qubits n
1
2
3
4
5
6
7
8
9

0
1
1
1
1
1
1
1
1

Hamming weight w
1
2
3
1
2
2
3
3
4
4

2
2
4
7
8

2
4
≤ 11

4

5
> 10?

Since our heuristic algorithm finds only an upper bound on the stabilizer rank, and not a lower
bound, by construction we cannot guarantee any statement on the scaling of the rank itself. However, our approach could have found only stabilizer decompositions of very low rank, thereby
providing evidence that Dicke states have very slowly growing rank, meaning that stabilizer-rank
methods can efficiently simulate circuits which output Dicke states. This is not what we observe;at the very least we can say that, if Dicke states have low stabilizer rank, then the current stateof-the-art method by Bravyi et al. does not succeed in finding the corresponding decomposition.
Further research is needed for a conclusive answer.
We now explain the heuristic algorithm by Bravyi et al. [12], which has been explained in moredetail in [84]. The algorithm follows a simulated annealing approach: on input n, w and χ, it
performs a random walk through sets of χ stabilizer states. It starts with a random set V of χ
stabilizer states on n qubits. In a single ‘step’, the algorithm picks one of these states |ψ⟩ ∈ V
at random, together with a random n-qubit Pauli operator P , and replaces the state |ψ⟩ with|ψ ′ ⟩ := c(I + P ) |ψ⟩ with c a normalization constant (or repeats if |ψ ′ ⟩ = 0), yielding a new setnn⟩ | with ΠV| ΠV |Dw
V ′ . The step is accepted with certainty if FV < FV ′ , where FV := | ⟨Dw
the projector on the subspace of the n-qubit Hilbert space spanned by the stabilizer states in V .
Otherwise, it is accepted with probability exp(−β(FV ′ − FV )), where β should be interpreted asnthe inverse temperature. The algorithm terminates if it finds FV = 1, implying that |Dw⟩ can bewritten as linear combination of V , outputting the number χ as (an upper bound on) the stabilizer
rank of |ψ⟩. For a fixed χ, we use identical values to Bravyi et al. [12] and vary β from 1 to 4000
in 100 steps, performing 1000 steps at each value of β.

Accepted in Quantum 2023-07-17, click title to verify. Published under CC-BY 4.0.

60