# MQT Core DD

MQT Core represents and manipulates quantum states and operations with decision
diagrams (DDs). The C++ library and [`mqt.core.dd`](api/mqt/core/dd/index.html.md#module-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`](api/mqt/core/mlir/index.html.md#mqt.core.mlir.QCOProgram) objects. The simulator supports
mid-circuit measurements, resets, and classically controlled operations. This
example compiles and samples a Bell-state program:

```ipython3
from mqt.core.mlir import sample

bell_qasm = """OPENQASM 3.0;
include "stdgates.inc";
qubit[2] q;
bit[2] result;
h q[0];
cx q[0], q[1];
result = measure q;
"""

counts = sample(bell_qasm, shots=1024, seed=1)
print(counts)
```

```myst-ansi
{'00': 506, '11': 518}
```

The [`sample()`](api/mqt/core/mlir/index.html.md#mqt.core.mlir.sample), [`simulate()`](api/mqt/core/mlir/index.html.md#mqt.core.mlir.simulate), and
[`build_functionality()`](api/mqt/core/mlir/index.html.md#mqt.core.mlir.build_functionality) functions accept source text,
paths, Qiskit circuits, and typed compiler programs. They lower each input
directly to QCO. The corresponding [`QCOProgram`](api/mqt/core/mlir/index.html.md#mqt.core.mlir.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()`](api/mqt/core/dd/index.html.md#mqt.core.dd.VectorDD.get_vector) or
[`get_matrix()`](api/mqt/core/dd/index.html.md#mqt.core.dd.MatrixDD.get_matrix).

```ipython3
import numpy as np
from mqt.core.dd import DDPackage
from mqt.core.mlir import QCOProgram, build_functionality, simulate

unitary_program = QCOProgram.from_mlir_str("""
module {
  func.func @main() attributes {mqt.entry_point} {
    %q0 = qco.static 0 : !qco.qubit
    %q1 = qco.static 1 : !qco.qubit
    %q0_h = qco.h %q0 : !qco.qubit -> !qco.qubit
    %q0_out, %q1_out = qco.ctrl(%q0_h) targets(%target = %q1) {
      %target_out = qco.x %target : !qco.qubit -> !qco.qubit
      qco.yield %target_out : !qco.qubit
    } : ({!qco.qubit}, {!qco.qubit}) -> ({!qco.qubit}, {!qco.qubit})
    qco.sink %q0_out : !qco.qubit
    qco.sink %q1_out : !qco.qubit
    return
  }
}
""")

vec = simulate(unitary_program)
unitary = build_functionality(unitary_program)

dd = DDPackage(2)
zero_state_dd = dd.zero_state(2)
out_state_dd = unitary_program.simulate(zero_state_dd, dd)
vec = np.array(out_state_dd.get_vector(), copy=False)
with np.printoptions(precision=3, suppress=True):
  print(vec)

functionality_dd = unitary_program.build_functionality(dd)
unitary = np.array(functionality_dd.get_matrix(2), copy=False)
with np.printoptions(precision=3, suppress=True):
  print(unitary)
```

```myst-ansi
[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()`](api/mqt/core/dd/index.html.md#mqt.core.dd.VectorDD.to_svg) to render a decision diagram as SVG.
It uses
[PyGraphviz](https://pygraphviz.github.io/documentation/stable/install.html) 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.

```ipython3
from IPython.display import SVG

out_state_dd.to_svg("bell_state.svg")
SVG(filename="bell_state.svg")
```

```myst-ansi
Fontconfig warning: "/usr/share/fontconfig/conf.avail/05-reset-dirs-sample.conf", line 6: unknown element "reset-dirs"
```

See [`DDPackage`](api/mqt/core/dd/index.html.md#mqt.core.dd.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](references.html.md#id8)]. This led to the emergence of a
wide variety of decision diagrams, including BDDs, FBDDs, KFDDs, MTBDDs, and
ZDDs (see, for example,
[[4](references.html.md#id9), [5](references.html.md#id39), [6](references.html.md#id17), [7](references.html.md#id16), [8](references.html.md#id7), [9](references.html.md#id26)]),
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](references.html.md#id41), [11](references.html.md#id42), [12](references.html.md#id25), [13](references.html.md#id30), [14](references.html.md#id45), [15](references.html.md#id22), [16](references.html.md#id37)].
Particularly for design tasks like *simulation*
[[16](references.html.md#id37), [17](references.html.md#id36), [18](references.html.md#id44), [19](references.html.md#id20), [20](references.html.md#id12), [21](references.html.md#id19), [22](references.html.md#id14), [23](references.html.md#id18), [24](references.html.md#id23), [25](references.html.md#id31)],
*synthesis*
[[26](references.html.md#id29), [27](references.html.md#id5), [28](references.html.md#id34), [29](references.html.md#id46), [30](references.html.md#id6), [31](references.html.md#id24)],
and *verification*
[[32](references.html.md#id11), [33](references.html.md#id13), [34](references.html.md#id15), [35](references.html.md#id38), [36](references.html.md#id33), [37](references.html.md#id21)]
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,

<a id="equation-ssstate"></a>
$$
\ket{\Psi} = \alpha_0 \ket{0} + \alpha_1 \ket{1},
$$

which is commonly represented as a statevector

$$
\ket{\Psi}\equiv \begin{bmatrix} \alpha_0 & \alpha_1	\end{bmatrix}^\top.
$$

The vector in [(1)](#equation-ssstate) splits into the contributions of the $\ket{0}$ state
($\alpha_0$) and the $\ket{1}$ state ($\alpha_1$):

<a id="equation-splitting"></a>
$$
\bigl(
\overbrace{
\overset{\ket{0}}{\begin{bmatrix} \alpha_0
\end{bmatrix}}
\ \ \
\overset{\ket{1}}{\begin{bmatrix} \alpha_1
\end{bmatrix}}
}^{\ket{\Psi}}
\bigr)^\top.
$$

This decomposition is the core of the decision-diagram formalism. The decision
diagram representing $\ket{\Psi}$ has the structure

![One qubit with outgoing edges weighted by its zero and one amplitudes.](_static/dd-figure-01.svg)

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)](#equation-splitting) 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.

Building off the intuition of a single-qubit state, we can move to larger
systems.

The diagrams above represent each part of the statevector separately. Merging
redundant subgraphs makes the representation compact.

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](references.html.md#id30)]. 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](references.html.md#id20)]. 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`](cpp/classdd_1_1ComplexNumbers.html#ac59b134735021aebca3264ba52fe15f0). 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.

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](#dd-three-qubits)) 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](references.html.md#id43)] 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,

$$
U &= \begin{bmatrix}
U_{00} & U_{01} \\ U_{10} & U_{11}
\end{bmatrix} = U_{00} \ket{0}\!\bra{0} + U_{01} \ket{0}\!\bra{1} + U_{10} \ket{1}\!\bra{0} + U_{11} \ket{1}\!\bra{1} .
$$

Then, the decision diagram representing this matrix has the structure

![Matrix DD with four successors in row-major order.](_static/dd-figure-07.svg)

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}$.

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.

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

<!-- prettier-ignore -->

<a id="equation-kronecker"></a>
$$
\ket{\Psi} \otimes \ket{\Phi} = \begin{bmatrix}
\Psi_{0} \ket{\Phi} \\
\Psi_{1} \ket{\Phi}
\end{bmatrix}
= \begin{bmatrix}
\Psi_{0} \begin{bmatrix} \Phi_0 \\ \Phi_1 \end{bmatrix} \\
\Psi_{1} \begin{bmatrix} \Phi_0 \\ \Phi_1 \end{bmatrix}
\end{bmatrix}.
$$

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:

![Kronecker product replacing terminal edges with the second DD.](_static/dd-figure-10.svg)

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

<!-- prettier-ignore -->

<a id="equation-addition"></a>
$$
\ket{\Psi} + \ket{\Phi} = \begin{bmatrix} \Psi_0 \\ \Psi_1 \end{bmatrix} + \begin{bmatrix} \Phi_0 \\ \Phi_1 \end{bmatrix} = w \begin{bmatrix} \alpha_0  \\ \alpha_1 \end{bmatrix} + w' \begin{bmatrix} \alpha'_0 \\ \alpha'_1 \end{bmatrix} = \begin{bmatrix} w \alpha_0 + w' \alpha'_0 \\ w \alpha_1 + w' \alpha'_1 \end{bmatrix},
$$

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,

![Addition recursively combining corresponding weighted successors.](_static/dd-figure-11.svg)

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

<a id="equation-multiplication"></a>
$$
U\ket{\Psi} = \begin{bmatrix} U_{00} & U_{01} \\
U_{10} & U_{11} \end{bmatrix} \begin{bmatrix} \Psi_0 \\ \Psi_1 \end{bmatrix}
 = w \begin{bmatrix} u_{00} & u_{01} \\
u_{10} & u_{11} \end{bmatrix} w' \begin{bmatrix} \alpha_0 \\ \alpha_1 \end{bmatrix} = ww' \begin{bmatrix} u_{00} \cdot \alpha_0 + u_{01} \cdot \alpha_1 \\
u_{10} \cdot \alpha_0 + u_{11} \cdot \alpha_1 \end{bmatrix}.
$$

This implies that a multiplication boils down to four smaller multiplications
and two additions. In the decision-diagram formalism, this has the form

![Matrix-vector product combining each matrix row with the vector.](_static/dd-figure-12.svg)

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

<!-- prettier-ignore -->

<a id="equation-innerproduct"></a>
$$
\langle\Psi \vert \Phi\rangle = \begin{bmatrix} \Psi^*_0 & \Psi^*_1 \end{bmatrix} \begin{bmatrix} \Phi_0 \\ \Phi_1 \end{bmatrix}
= w^* \begin{bmatrix} \alpha^*_0 & \alpha^*_1 \end{bmatrix} w' \begin{bmatrix} \alpha'_0 \\ \alpha'_1 \end{bmatrix} = w^*w' (\alpha^*_0 \alpha'_0 + \alpha^*_1 \alpha'_1)
$$

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

![Inner product conjugating the first vector and summing paired branches.](_static/dd-figure-13.svg)

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.

```ipython3
matrix = np.array([[0.6, -0.8], [0.8, 0.6]], dtype=complex)
state = np.array([1, 1j], dtype=complex) / np.sqrt(2)
other = np.array([0, 1], dtype=complex)
package = DDPackage(1)
state_dd = package.from_vector(state)
other_dd = package.from_vector(other)
matrix_dd = package.from_matrix(matrix)
product = package.matrix_vector_multiply(matrix_dd, state_dd)
np.testing.assert_allclose(product.get_vector(), matrix @ state)
np.testing.assert_allclose(package.inner_product(state_dd, other_dd), np.vdot(state, other))
np.testing.assert_allclose(package.vector_add(state_dd, other_dd).get_vector(), state + other)
with np.printoptions(precision=3, suppress=True):
    print("Matrix-vector product:", np.asarray(product.get_vector()))
    print("Inner product:", package.inner_product(state_dd, other_dd))
```

```myst-ansi
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.

```ipython3
print("Qubits | Left nodes | Right nodes | Sum nodes")
for n in (4, 6, 8):
    package = DDPackage(n)
    left = np.ones(1, dtype=complex)
    right = left.copy()
    for j in range(n):
        left = np.kron(left, [1, 1]) / np.sqrt(2)
        angle = 0.2 + 0.031 * j
        right = np.kron(right, [np.cos(angle), np.sin(angle)])
    left_dd = package.from_vector(left)
    right_dd = package.from_vector(right)
    result_dd = package.vector_add(left_dd, right_dd)
    np.testing.assert_allclose(result_dd.get_vector(), left + right, atol=1e-10)
    print(n, left_dd.size() - 1, right_dd.size() - 1, result_dd.size() - 1, sep=" | ")
```

```myst-ansi
Qubits | Left nodes | Right nodes | Sum nodes
4 | 4 | 4 | 15
6 | 6 | 6 | 63
8 | 8 | 8 | 255
```
