MQT Core DD¶
MQT Core represents and manipulates quantum states and operations with decision
diagrams (DDs). The C++ library and mqt.core.dd Python module support
simulation, synthesis, and verification. Start with the quickstart for Python
usage or the introduction below for the data structure and its limits.
Quickstart¶
The MQT Compiler Collection uses the DD package to simulate
QCOProgram objects. The simulator supports
mid-circuit measurements, resets, and classically controlled operations. This
example compiles and samples a Bell-state program:
1from mqt.core.mlir import sample
2
3bell_qasm = """OPENQASM 3.0;
4include "stdgates.inc";
5qubit[2] q;
6bit[2] result;
7h q[0];
8cx q[0], q[1];
9result = measure q;
10"""
11
12counts = sample(bell_qasm, shots=1024, seed=1)
13print(counts)
{'00': 506, '11': 518}
The sample(), simulate(), and
build_functionality() functions accept source text,
paths, Qiskit circuits, and typed compiler programs. They lower each input
directly to QCO. The corresponding QCOProgram methods
provide the DD-native interface for reusable compiled programs, custom initial
states, and dynamic simulation. The top-level simulate function starts a
closed program in the all-zero state. It and build_functionality manage the DD
package internally and materialize their results directly into NumPy arrays.
The QCOProgram methods avoid constructing exponentially large dense arrays
unless the result is explicitly converted with
get_vector() or
get_matrix().
1import numpy as np
2from mqt.core.dd import DDPackage
3from mqt.core.mlir import QCOProgram, build_functionality, simulate
4
5unitary_program = QCOProgram.from_mlir_str("""
6module {
7 func.func @main() attributes {mqt.entry_point} {
8 %q0 = qco.static 0 : !qco.qubit
9 %q1 = qco.static 1 : !qco.qubit
10 %q0_h = qco.h %q0 : !qco.qubit -> !qco.qubit
11 %q0_out, %q1_out = qco.ctrl(%q0_h) targets(%target = %q1) {
12 %target_out = qco.x %target : !qco.qubit -> !qco.qubit
13 qco.yield %target_out : !qco.qubit
14 } : ({!qco.qubit}, {!qco.qubit}) -> ({!qco.qubit}, {!qco.qubit})
15 qco.sink %q0_out : !qco.qubit
16 qco.sink %q1_out : !qco.qubit
17 return
18 }
19}
20""")
21
22vec = simulate(unitary_program)
23unitary = build_functionality(unitary_program)
24
25dd = DDPackage(2)
26zero_state_dd = dd.zero_state(2)
27out_state_dd = unitary_program.simulate(zero_state_dd, dd)
28vec = np.array(out_state_dd.get_vector(), copy=False)
29with np.printoptions(precision=3, suppress=True):
30 print(vec)
31
32functionality_dd = unitary_program.build_functionality(dd)
33unitary = np.array(functionality_dd.get_matrix(2), copy=False)
34with np.printoptions(precision=3, suppress=True):
35 print(unitary)
[0.707+0.j 0. +0.j 0. +0.j 0.707+0.j]
[[ 0.707+0.j 0.707+0.j 0. +0.j 0. +0.j]
[ 0. +0.j 0. +0.j 0.707+0.j -0.707+0.j]
[ 0. +0.j 0. +0.j 0.707+0.j 0.707+0.j]
[ 0.707+0.j -0.707+0.j 0. +0.j 0. +0.j]]
Use to_svg() to render a decision diagram as SVG.
It uses
PyGraphviz 2
or later when installed, or the dot command otherwise. IPython can display the
resulting file in a notebook. DOT exports use unique node IDs assigned in
traversal order, so ordinary exports do not depend on memory addresses. The
memory=True option includes addresses as debugging information.
1from IPython.display import SVG
2
3out_state_dd.to_svg("bell_state.svg")
4SVG(filename="bell_state.svg")
See DDPackage for the full API.
Memory and Cache Growth¶
Node unique tables start with 64 buckets per qubit level and grow independently as levels fill. This avoids reserving large tables for sparsely populated levels. The matrix-vector cache starts with 16,384 entries. After garbage collection, it can grow to accommodate surviving vector nodes when past cache hits indicate reuse. Automatic growth stops at 1,048,576 entries. Other compute caches keep their initial capacities.
This policy trades memory for less recomputation; it can slow workloads whose
useful cache entries already fit. The cache ceiling does not bound total package
memory. C++ callers can still set initial capacities with dd::DDPackageConfig;
an initial matrix-vector capacity at or above the growth ceiling stays fixed.
Collection and reset retain grown table capacities. Node pools retain their
slabs until reset or destruction and zero fresh entries on acquisition, so
unused reserved slots need not occupy resident memory.
How do Quantum Decision Diagrams Work?¶
Decision diagrams were introduced in the 1980s as a data structure for the efficient representation and manipulation of Boolean functions [3]. This led to the emergence of a wide variety of decision diagrams, including BDDs, FBDDs, KFDDs, MTBDDs, and ZDDs (see, for example, [4, 5, 6, 7, 8, 9]), which made them a crucial tool in the development of modern circuits and systems. Because of their previous success, decision diagrams have been proposed for application in the realm of quantum computing [10, 11, 12, 13, 14, 15, 16]. Particularly for design tasks like simulation [16, 17, 18, 19, 20, 21, 22, 23, 24, 25], synthesis [26, 27, 28, 29, 30, 31], and verification [32, 33, 34, 35, 36, 37] of quantum circuits, they recently attracted great attention.
The following sections explain how decision diagrams represent quantum states and operations, and how computations act on those representations.
Representation of Quantum States¶
First, we review how quantum states are represented using decision diagrams. To this end, we consider the simple case of a single-qubit system. The state \(\ket{\Psi}\) of such a system is described by two complex-valued, normalized amplitudes \(\alpha_0\) and \(\alpha_1\), that is,
which is commonly represented as a statevector
The vector in (1) splits into the contributions of the \(\ket{0}\) state (\(\alpha_0\)) and the \(\ket{1}\) state (\(\alpha_1\)):
This decomposition is the core of the decision-diagram formalism. The decision diagram representing \(\ket{\Psi}\) has the structure
It consists of a single node with one incoming edge that represents the entry point in the decision diagram, as well as two successors that represent the split shown in (2) and end in a terminal node (the black box). The state’s amplitudes are annotated at the respective edges. Edges without annotations correspond to an edge weight of 1.
Example (Single-Qubit States)
Consider the computational basis states \(\ket{0}\) and \(\ket{1}\). Then, the corresponding decision diagrams have the structures
and
In each of the cases, one of the successors ends in the terminal node, while the other ends in a zero stub (indicated by a black dot)—uncannily resembling the corresponding vector descriptions.
Building off the intuition of a single-qubit state, we can move to larger systems.
Example (Multi-Qubit States)
Consider the following statevector of a three-qubit system:
Then, \(\ket{\Psi}\) can be recursively split into equally-sized parts similar to (2), i.e.,
where \(q_2, q_1, q_0 \in \{0, 1\}\). This directly translates to the decision-diagram formalism:
Each level of the decision diagram consists of decision nodes with corresponding left and right successor edges. These successors represent the path that leads to an amplitude where the local quantum system (corresponding to the level of the node, annotated here with the labels) is in the \(\ket{0}\) (left successor) or the \(\ket{1}\) state (right successor).
The diagrams above represent each part of the statevector separately. Merging redundant subgraphs makes the representation compact.
Example (Redundancy in Decision Diagrams)
Observe how, as in the previous example, the left and right successors of the top-level node (labeled \(q_2\)) lead to exactly the same structure (highlighted by dashed rectangles in the unreduced diagram). As a result, the whole sub-diagram does not need to be represented twice, i.e.,
From a memory perspective, this reduction alone has compressed the overall memory required to represent the state by 50%.
Identifying redundancies in these kinds of representations heavily depends on the use of what is referred to as a normalization scheme for the decision diagram nodes [13]. Such a normalization scheme makes sure two decision diagram nodes that represent the same functionality do indeed have the same numerical structure. In computer science, this property is called canonicity.
The most widely used and practically relevant normalization scheme is to normalize the outgoing edges of a node by dividing both weights by the norm of the vector containing both edge weights and extracting a common phase into the incoming edge [19]. This normalizes the sum of the squared magnitudes of the outgoing edge weights to \(1\) and is consistent with quantum semantics, where basis states \(\ket{0}\) and \(\ket{1}\) are observed after measurement with probabilities that are squared magnitudes of the respective weights. MQT Core selects a maximum-magnitude edge (preferring the left edge when squared magnitudes agree within relative tolerance) and makes its normalized weight real and nonnegative. The incoming edge retains its complex phase. Normalization proceeds bottom-up; complex-number comparisons use the package tolerance. Recursive addition divides both operands by their largest component before visiting child nodes and restores that scale on return. This keeps small basis-state amplitudes from being discarded before the normalized parent is reconstructed, while making the dominant component exactly one for subgraph sharing. The same rule applies to magnitude addition.
Cached vector normalization projects nearly equal or opposite coefficients onto the corresponding balanced pair before dividing by their norm. For unit-norm child states, this changes the local vector by at most the absolute tolerance in Euclidean norm, apart from roundoff. The incoming weight also compensates for reuse of a stored dominant weight. These steps limit amplification of small coefficient differences; they do not bound accumulated circuit error.
Floating-point arithmetic makes this canonicity approximate. Ordinary real
components reuse the nearest stored value within the absolute tolerance,
preferring the smaller magnitude on a tie; zero, one, and \(1/\sqrt{2}\) have
priority. The default tolerance is \(2^{-42}\) (1024 times double-precision
machine epsilon). The numeric index hashes binary intervals without rounding
stored values. C++ callers can set a positive, normal global tolerance with
dd::ComplexNumbers::setTolerance. A smaller tolerance can reduce
error amplification in small subproblems, but may also prevent sharing of nearly
equal subgraphs. Neither tolerance choice guarantees polynomial DD size for a
circuit. Set the tolerance before creating packages: changing it does not
recanonicalize existing decision diagrams.
Example (Normalization of Decision Diagrams)
Considering the decision diagram from the previous example, this results in the following normalized and reduced decision diagram:
The first two levels (\(q_2\) and \(q_1\)) of the above diagram naturally encode that the respective qubits have a \(50/50\) chance to be in \(\ket{0}\) and \(\ket{1}\) (since \(\vert1/\sqrt{2}\vert^2 = 0.5\)). Meanwhile, the bottom level (\(q_0\)) encodes that the probability of \(q_0\) depends on the state of \(q_1\). If \(q_1\) is in the \(\ket{0}\) state (following the left successor), then \(q_0\) has probability \(0.5\) for both \(\ket{0}\) and \(\ket{1}\). If \(q_1\) is in the \(\ket{1}\) state (following the right successor), it is guaranteed that the remaining qubit is in the \(\ket{0}\) state.
A statevector DD recursively halves the vector and shares redundant subgraphs. This representation has the following properties:
Decision diagrams can be initialized in their compact form (as, for example, shown in the last example above). There is no need to create the maximally large decision diagram (as shown, for example, in the unreduced diagram) at any point in a calculation.
Determining a particular amplitude of the represented state corresponds to multiplying the edge weights along a single-path traversal from the top edge of the decision diagram (called its root) to a terminal node.
The efficiency of decision diagrams is commonly measured by their size, that is, the number of nodes in the decision diagram—the smaller the number of nodes, the higher the compaction achieved by the data structure. Note that the terminal (node) is typically not counted towards the size of a decision diagram.
Any product state naturally has a decision diagram consisting of a single node per site. However, a compact DD does not correlate with the state being trivial. Even entangled states such as the GHZ state or the W state have decision diagrams whose size (that is, the number of nodes) is linear in the number of qubits.
The worst-case size, for states without redundancy, is exponential in the number of qubits. More specifically, a maximally large decision diagram has \(1+2^1+2^2+\dots+2^{n-1} = 2^n-1\) nodes.
To reduce visual clutter in illustrations of decision diagrams, edge weights are commonly not explicitly annotated, but their magnitude and phase are reflected in the thickness and the color of the respective edge. In addition, to make the correspondence of the individual levels in a decision diagram to a system’s qubits more explicit, the nodes are frequently annotated with the qubit’s index as an identifier. See [38] for further details on common techniques for visualization of decision diagrams.
Representation of Quantum Operations¶
Quantum operations are fundamentally described by complex-valued matrices. Matrix decision diagrams are a natural extension to vector decision diagrams by an additional dimension. To this end, consider the base case of a \(2\times 2\) matrix \(U\), that is,
Then, the decision diagram representing this matrix has the structure
which again resembles the general structure of the matrix. Note that \(U_{ij}\) can be interpreted as the transformation of \(\ket{j}\) to \(\ket{i}\).
Example (Single-Qubit Operations)
The following shows decision diagram representations for selected single-qubit operations:
The last equivalence demonstrates how a common factor between the edge weights can be pulled out and attached to the incoming (root) edge.
The generalization to larger matrices works analogously to the vector case. To construct the decision diagram representing a matrix, the matrix is recursively divided into quarters, and the four elements correspond to the four successors of the node to represent that split. As for vector decision diagrams, a normalization scheme makes the representation canonical so equivalent subgraphs can be shared. Each node’s outgoing edge weights are divided by the weight with the highest magnitude, selecting the leftmost one in a tie. The normalized outgoing weights have magnitude at most \(1\); the extracted factor moves to the incoming edge.
The matrix root carries the global scale. Its real components use a separate exact index, retaining tolerance-based priority for nonzero special constants. This preserves roots below the ordinary absolute tolerance, such as the \(2^{-64}\) root of \(H^{\otimes 128}\). Internal normalized coefficients still use the ordinary tolerance. Matrix normalization removes a common power-of-two scale before squared magnitudes and division, then retains the original root weight. Small local matrix entries remain subject to zero tolerance, and intermediate and final values must still fit the floating-point representation.
Example (Matrix Decision Diagrams)
Consider the maximally-entangling two-qubit \(R_{xx}\) rotation represented by the matrix
This matrix is equivalent to blocks of \(2 \times 2\) matrices corresponding to the identity \(I\) and the Pauli-\(X\) matrix, i.e.,
The corresponding (already reduced) decision diagram has the following structure:
Notice how the decision diagram naturally resembles the structure of the matrix. The nodes at the bottom represent the identity and the \(X\) matrix while the node at the top encodes the redundancy of the upper left quadrant and the bottom right quadrant, as well as the upper right and lower left quadrant in (3). Similarly to the vector example above, exploiting redundancy has halved the overall memory requirement.
Again, some interesting properties to point out:
Just as in the vector case, it is always possible to work with the reduced form of matrix decision diagrams right away, that is, without ever constructing the exponentially-sized, maximally-large diagram.
A maximally-large matrix decision diagram for \(n\) qubits has \(\sum_{i=1}^n 4^{i-1} = \frac{(4^n -1)}{3}\) nodes.
Decision diagrams are not limited to local interactions. Even long-range interactions between arbitrary qubits typically produce compact representations as decision diagrams. For example, any two-qubit interaction between arbitrary qubits can be represented as a decision diagram with at most \(1+4(n-1)\) nodes—an exponential reduction.
Decision diagrams are not limited to two-qubit interactions either. For example, controlled quantum gates with arbitrarily many controls (such as the multi-controlled Toffoli gate) give rise to decision diagrams with a linear number of nodes.
Fundamental Operations on Decision Diagrams¶
DD operations recursively split computations along the graph structure and cache shared subproblems. Their cost depends on the distinct subproblems visited and the size of the result, as described below. The examples use vectors; the same recursive approach extends to matrices.
Kronecker Product¶
The Kronecker product is necessary to create product states and to chain together local operations. For vectors, it can be expressed as
The DD Kronecker product replaces the nonzero terminal edges of the first diagram with the root edge of the second, multiplying their weights. For the example above:
As such, its complexity is linear in the number of nodes of the first decision diagram.
Addition¶
Standard vector addition can be recursively broken down according to
where \(w\) and \(w'\) are common factors of the terms in \(\ket{\Psi}\) and \(\ket{\Phi}\), respectively.
In the decision-diagram formalism, this corresponds to a simultaneous traversal of both decision diagrams from their roots to the terminal (multiplying edge weights along the way until the individual amplitudes are reached) and back again (accumulating the results of the recursive computations). More precisely,
where the dashed nodes represent the respective successor decision diagrams. The cost depends on the distinct weighted subproblems and the resulting DD. Even two compact inputs can produce an exponentially large sum; input node counts alone do not give a linear time bound.
Matrix-Vector Multiplication¶
Matrix-vector multiplication can be handled in a very similar fashion as addition. Standard matrix-vector multiplication can be expressed as
This implies that a multiplication boils down to four smaller multiplications and two additions. In the decision-diagram formalism, this has the form
where the dashed nodes again represent the respective successor decision diagrams. Runtime depends on the distinct weighted subproblems, intermediate additions, and output size. Cache reuse can reduce repeated work; compact input DDs alone do not guarantee a compact result.
Inner Product¶
Computing the inner product of two vectors can be recursively broken down according to
This implies that the inner product boils down to two smaller inner product computations and adding the results. As with the matrix-vector multiplication, this is done recursively for each level of the decision diagram. In the decision-diagram formalism, this has the following form
The recursion visits pairs of subdiagrams and reuses cached results. Its cost depends on the pairs visited and cache reuse, rather than only the size of the larger input.
Check the algebra with complex amplitudes¶
The same operations can be compared directly with NumPy. A nonsymmetric matrix makes row/column mistakes visible, while complex amplitudes exercise conjugation and phase handling.
1matrix = np.array([[0.6, -0.8], [0.8, 0.6]], dtype=complex)
2state = np.array([1, 1j], dtype=complex) / np.sqrt(2)
3other = np.array([0, 1], dtype=complex)
4package = DDPackage(1)
5state_dd = package.from_vector(state)
6other_dd = package.from_vector(other)
7matrix_dd = package.from_matrix(matrix)
8product = package.matrix_vector_multiply(matrix_dd, state_dd)
9np.testing.assert_allclose(product.get_vector(), matrix @ state)
10np.testing.assert_allclose(package.inner_product(state_dd, other_dd), np.vdot(state, other))
11np.testing.assert_allclose(package.vector_add(state_dd, other_dd).get_vector(), state + other)
12with np.printoptions(precision=3, suppress=True):
13 print("Matrix-vector product:", np.asarray(product.get_vector()))
14 print("Inner product:", package.inner_product(state_dd, other_dd))
Matrix-vector product: [0.424-0.566j 0.566+0.424j]
Inner product: -0.7071067811865476j
Compact inputs can have a large sum¶
Both inputs below are product states with one nonterminal node per qubit. Their
sum needs many more nodes. The public size() includes the terminal; subtract
one when reporting nonterminal nodes.
1print("Qubits | Left nodes | Right nodes | Sum nodes")
2for n in (4, 6, 8):
3 package = DDPackage(n)
4 left = np.ones(1, dtype=complex)
5 right = left.copy()
6 for j in range(n):
7 left = np.kron(left, [1, 1]) / np.sqrt(2)
8 angle = 0.2 + 0.031 * j
9 right = np.kron(right, [np.cos(angle), np.sin(angle)])
10 left_dd = package.from_vector(left)
11 right_dd = package.from_vector(right)
12 result_dd = package.vector_add(left_dd, right_dd)
13 np.testing.assert_allclose(result_dd.get_vector(), left + right, atol=1e-10)
14 print(n, left_dd.size() - 1, right_dd.size() - 1, result_dd.size() - 1, sep=" | ")
Qubits | Left nodes | Right nodes | Sum nodes
4 | 4 | 4 | 15
6 | 6 | 6 | 63
8 | 8 | 8 | 255