graphix.sim.statevec module¶

Numba JIT-compiled MBQC statevector backend simulator based on Ref. [1].

[1] McGuffin, M. J., Robert J-M., and Ikeda K. “How to Write a Simulator for Quantum Circuits from Scratch: A Tutorial.”, 2025 (arXiv:2506.08142).

NUM_QUBIT_PARALLEL = 15¶: This compilation constant determines the number of qubits above which matrix operations are multi-threaded. For lower counts, the advantage from parallelization is lost due to the overhead of initalizing multiple threads. This number was determined empirically and may be platform dependent.

class Statevector(data: Data = graphix.states.PlanarState(Plane.XY, 0), nqubit: int | None = None, max_qubits: int | None = None)[source]¶

Bases: DenseState

Statevector object.

_psi¶

Complex-valued 1-dimensional array representing the quantum statevector. Only the first 2**nqubit complex values have meaning.

Type:: npt.NDArray[np.complex128]

_nqubit¶

Number of active qubits at any given time.

Type:: int

_max_qubits¶

Maximum Hilbert space size allowed for internal computations. It determines the size of self._psi. For circuit simulations, it corresponds to the number of qubits, while for pattern simulations it corresponds to the pattern’s maximum space. The method Statevector.ensure_capacity() allows to increase this number.

Type:: int

Notes

The internal representation of the quantum state is guaranteed to be normalized after initialization, and it is assumed to remain normalized thereafter.

Using evolve_single(), expectation_single(), evolve(), or expectation_value() with non-unitary operators does not preserve the norm of the statevector and may lead to unexpected behavior.

In pattern simulation, node measurements call evolve_single() with a projector (a non-unitary operator). However, the measured qubit is immediately removed via remove_qubit(), which restores the unit norm of the internal quantum state. See graphix.sim.base_backend.DenseStateBackend.measure() for additional details.

__init__(data: Data = graphix.states.PlanarState(Plane.XY, 0), nqubit: int | None = None, max_qubits: int | None = None) → None[source]¶

Initialize a statevector object.

data can be: - a single graphix.states.State (classical description of a quantum state) - an iterable of graphix.states.State objects - an iterable of scalars (a \(2^n\) numerical statevector) - a single graphix.sim.statevec.Statevector

If nqubit is not provided, it is inferred from data. If max_qubits is not provided, it is set to match the inferred max_qubits. If only one graphix.states.State is provided and nqubit is a valid integer, the statevector is initialized in the tensor product state. If a class:graphix.sim.statevec.Statevector is provided, a copy is returned. Consistency between provided nqubit, max_qubits and data is checked.

Parameters:

data (Data, optional) – Input data to prepare the state. Can be a classical description or a numerical input, defaults to graphix.states.BasicStates.PLUS.
nqubit (int | None, optional) – Number of qubits to prepare. If None (default), it’s inferred from data.
max_qubits (int | None, optional.) – Maximum Hilbert space size for array preallocation. If None (default), it’s set equal to nqubit.

Raises:

ValueError – If nqubit, max_qubits or data are not consistent with each other.

property psi: NDArray[complex128]¶

View of the meaningful elements in self._psi.

These are the first 2**self.nqubit elements.

property nqubit: int¶: Number of qubits.

property max_qubits: int¶: Preallocated number of qubits.

ensure_capacity(required_qubits: int) → None[source]¶

Extend the state vector if the required qubit capacity exceeds the current one.

It copies the full state vector if required_qubits > self.max_qubits, otherwise does nothing.

Parameters:: required_qubits (int) – Minimum number of qubits the state vector must support. If expansion is needed, self._psi is extended to size 2**required_qubits.

flatten() → NDArray[object_ | complex128][source]¶

Return flattened state vector.

A view of only the first 2**self.nqubit elements of self._psi is returned.

add_nodes(nqubit: int, data: Data) → None[source]¶

Add nodes (qubits) to the state vector and initialize them in a specified state.

Previously existing nodes remain unchanged.

Parameters:

nqubit (int) – The number of qubits to add to the state vector.
data (Data) –
The state in which to initialize the newly added nodes.
- If a single basic state is provided, all new nodes are initialized in that state.
- If a list of basic states is provided, it must match the length of nodes, and each node is initialized with its corresponding state.
- A single-qubit state vector will be broadcast to all nodes.
- A multi-qubit state vector of dimension \(2^n\), where \(n = \mathrm{len}(nodes)\), initializes the new nodes jointly.

Notes

This method can extend the size of self._psi for convenience, but this requires allocating a full new array.

entangle(qubits: tuple[int, int]) → None[source]¶

Apply a CZ gate on two qubits.

Parameters:: qubits (tuple[int, int]) – (control, target) qubit indices.

evolve_single(op: NDArray[object_ | complex128], qubit: int) → None[source]¶

Apply a single-qubit operator.

Parameters:

op (npt.NDArray[np.complex128]) – Complex-valued matrix of shape \((2, 2)\) representing the operator to apply.
qubit (int) – Target qubit index.

expectation_single(op: NDArray[object_ | complex128], qubit: int) → complex[source]¶

Return the expectation value of a single-qubit operator.

Parameters:

op (npt.NDArray[np.complex128]) – Complex-valued matrix of shape \((2, 2)\) representing the operator to measure.
qubit (int) – Target qubit index.

Returns:

Expectation value.

Return type:

complex

Notes

This method assumes that quantum state represented by self.psi is normalized. See the class docstring for details.

evolve(op: Matrix, qubits: Sequence[int]) → None[source]¶

Apply a multi-qubit operator.

Parameters:

op (npt.NDArray[np.complex128]) – Complex-valued matrix of shape \((2^n, 2^n)\) representing the operator to apply.
qubits (Sequence[int]) – Target qubit indices.

Notes

This method is a fallback for circuit simulation and it’s not required for pattern simulation. It does not have an efficient JIT-compiled implementation.

expectation_value(op: Matrix, qubits: Sequence[int]) → complex[source]¶

Return the expectation value of a multi-qubit operator.

Parameters:

op (npt.NDArray[np.complex128]) – Complex-valued matrix of shape \((2^n, 2^n)\) representing the operator to measure.
qubits (Sequence[int]) – Target qubit indices.

Notes

This method assumes that quantum state represented by self.psi is normalized. See the class docstring for details.

remove_qubit(qubit: int) → None[source]¶

Remove a separable qubit from the system and assemble the statevector of the remaining qubits.

This method is deprecated. See self.project_qubit().

Parameters:: qubit (int) – Target qubit index.

project_qubit(op: NDArray[object_ | complex128], qubit: int) → None[source]¶

Project out a qubit from the system and assemble the statevector of the remaining qubits.

This method evolves the statevector with op and removes qubit. It assumes that after the application of op, qubit is a separable qubit.

Let \(P \ket{\psi} = \sum c_i \ket{i}\) be the statevector after the application of the projector \(P\) (op), with sum taken over \(i \in [ 0 \dots 00,\ 0\dots 01,\ \dots,\ 1 \dots 11 ]\), this method returns

\[\begin{split}\ket{\psi}' =& c_{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k}}0_{\mathrm{k+1}} \dots 00} \ket{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k+1}} \dots 00} \\ & + c_{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k}}0_{\mathrm{k+1}} \dots 01} \ket{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k+1}} \dots 01} \\ & + c_{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k}}0_{\mathrm{k+1}} \dots 10} \ket{0 \dots 0_{\mathrm{k-1}}0_{\mathrm{k+1}} \dots 10} \\ & + \dots \\ & + c_{1 \dots 1_{\mathrm{k-1}}0_{\mathrm{k}}1_{\mathrm{k+1}} \dots 11} \ket{1 \dots 1_{\mathrm{k-1}}1_{\mathrm{k+1}} \dots 11},\end{split}\]

(after normalization) for \(k =\) qubit. If the \(k\) th qubit is in the \(\ket{1}\) state, all the amplitudes above will be zero. In that case the returned state will be the one above with \(0_{\mathrm{k}}\) replaced with \(1_{\mathrm{k}}\).

Warning

This method assumes that op is a projector, and therefore that qubit is a separable qubit after applying op. It is implemented as a significantly faster alternative for partial trace to be used after single-qubit measurements. Separability is not checked.

Parameters:

op (npt.NDArray[np.complex128]) – Complex-valued matrix of shape \((2, 2)\) representing the projector to apply.
qubit (int) – Target qubit index.

Notes

This method combines evolve_single() and remove_qubit().

First, the relevant entries of self.psi are partitioned into two disjoint subvectors. The partition depends on the value of q. For a three-qubit statevector [a0, a1, a2, a3, a4, a5, a6, a7]:

q = 0: sv0 = [a0, a1, a2, a3], sv1 = [a4, a5, a6, a7]
q = 1: sv0 = [a0, a1, a4, a5], sv1 = [a2, a3, a6, a7]
q = 2: sv0 = [a0, a2, a4, a6], sv1 = [a1, a3, a5, a7]

The statevector is then iterated over as in evolve_single(), updating self.psi in place while computing the norm of each partition separately.

If both resulting subvectors have nonzero norm, they are guaranteed to be identical. The subvector sv0 is normalized and written to the first 2**(nqubit - 1) entries of self.psi.

If exactly one subvector has nonzero norm, that subvector is selected. If both have zero norm, a RuntimeError is raised.

swap(qubits: tuple[int, int]) → None[source]¶

Apply SWAP gate between two qubits.

Parameters:: qubits (tuple[int, int]) – (control, target) qubit indices.

tensor(other: Statevector) → None[source]¶

Tensor product state with other qubits.

Results in self \(\otimes\) other.

Parameters:: other (graphix.sim.statevec.Statevector) – Statevector to be tensored with self.

Notes

This method is used internally by add_nodes().

permute(permutation: Sequence[int]) → None[source]¶

Reorder the qubits.

Parameters:: permutation (Sequence[int]) – The permutation to apply. For each position in the resulting order, the value gives the index of the qubit in the original ordering.

fidelity(other: Statevector | NDArray[object_ | complex128]) → float[source]¶

Calculate the fidelity against another statevector.

The fidelity is defined as \(|\langle\psi_1|\psi_2\rangle|^2\).

Parameters:: other (graphix.sim.statevec.Statevector) – Statevector to compare with.
Returns:: Fidelity between the two statevectors.
Return type:: float

isclose(other: Statevector | NDArray[object_ | complex128], *, rtol: float = 1e-09, atol: float = 0.0) → bool[source]¶

Check if two quantum states are equal up to global phase.

Two states are considered close if their fidelity is close to 1.

Parameters:

other (graphix.sim.statevec.Statevector) – Statevector to compare with.
rtol (float) – Relative tolerance for math.isclose().
atol (float) – Absolute tolerance for math.isclose().

Returns:

True if the states are equal up to global phase.

Return type:

bool

to_dict(encoding: _ENCODING = 'MSB', *, rtol: float = 0.0, atol: float = 1e-08) → dict[str, np.object_ | np.complex128][source]¶

Convert the statevector to dictionary form.

This dictionary representation uses a ket-like notation where the dictionary keys are qubit strings for the basis vectors and values are the corresponding complex amplitudes. Amplitudes below a certain threshold are filtered out.

Parameters:

encoding (Literal["LSB", "MSB"], default="MSB") – Encoding for the basis kets. See notes for additional information.
rtol (float, default=0.0) – Relative tolerance used when deciding whether a coefficient should be treated as zero. Values whose magnitude is within this relative tolerance of zero are omitted from the resulting dictionary.
atol (float, default=1e-8) – Absolute tolerance used when deciding whether a coefficient should be treated as zero. Values whose magnitude is within this relative tolerance of zero are omitted from the resulting dictionary.

Returns:

The statevector in dictionary form.

Return type:

dict[str, complex]

Notes

The encoding determines the bit ordering convention used when mapping basis states to dictionary keys. Consider a tensor product of three qubits:

\[\]

lvertpsirangle = q_0 otimes q_1 otimes q_2.

If encoding == "MSB" the first qubit is represented in the Most Significant Bit -> q0q1q2. This is the default representation in Graphix. If encoding == "LSB" the first qubit is represented in the Least Significant Bit -> q2q1q0. This is the default representation in other software packages such as Qiskit.

Example

>>> from graphix.states import BasicStates
>>> from graphix.sim.statevec import Statevector
>>> sv = Statevector(data=[BasicStates.ZERO, BasicStates.ONE])
>>> sv.to_dict()
{'01': np.complex128(1+0j)}
>>> sv.to_dict(encoding="LSB")
{'10': np.complex128(1+0j)}

to_prob_dict(encoding: _ENCODING = 'MSB', *, rtol: float = 0.0, atol: float = 1e-08) → dict[str, np.object_ | np.float64][source]¶

Convert the statevector to a probability distirbution in a dictionary form.

This dictionary representation uses a ket-like notation where the dictionary keys are qubit strings for the basis vectors and values are the corresponding probabilities. Basis vector whose amplitude is below a certain threshold are filtered out.

Parameters:

encoding (Literal["LSB", "MSB"], default="MSB") – Encoding for the basis kets. See to_dict() for additional information.
rtol (float, default=0.0) – Relative tolerance used when deciding whether a coefficient should be treated as zero. Values whose magnitude is within this relative tolerance of zero are omitted from the resulting dictionary.
atol (float, default=1e-8) – Absolute tolerance used when deciding whether a coefficient should be treated as zero. Values whose magnitude is within this relative tolerance of zero are omitted from the resulting dictionary.

Returns:

The probability distribution associated to the statevector in dictionary form.

Return type:

dict[str, float]