2.6. Simulators#

geqo has several simulators for simulating the execution of gates and quantum circuits.

  • ensembleSimulatorSymPy

  • mixedStateSimulatorSymPy

  • simulatorUnitarySymPy

  • ensembleSimulatorCuPy

  • mixedStateSimulatorCuPy

  • unitarySimulatorCuPy

  • statevectorSimulatorCuPy

  • simulatorStatevectorNumpy

Simulators based on SymPy can process symbolic and numeric values for angles and matrix entries. These simulators do not have a high performance and they are only suitable for rather small quantum systems.

Simulators based on NumPy can execute numeric simulations and have a higher performance. They are suitable for simulating systems with 10 to 20 qubits.

Simulators based on CuPy offer either GPU-accelerated simulations or vectorized numpy-based simulations in the absence of a GPU. These simulators have the highest performance and can implement larger circuits.

2.6.1. ensembleSimulatorSymPy#

This simulator keeps track of quantum states, measurement results and their corresponding probabilities. It supports unitary operations and other operations. For instance, it supports the non-unitary operations DropQubits and SetDensityMatrix as well as post-measurement operations and multiple rounds of measurement. Each round of measurement creates several branches of possible measurement outcomes, all of which are collected in an ensemble.

The output of the simulator is a dictionary. The keys of the dictionary are all values of the classical bits of the circuit. Different values for the classical bits can be obtained with Measure or SetBits in a quantum circuit. If there are no classical bits, then the list of bit values is empty. The value for each key is a pair. The first component of the pair is the probability for obtaining the key and the second component is the corresponding density matrix of the quantum part of the system.

The following example shows the result of the measurement of an EPR pair followed by a Hadamard operation and a second measurement.

from geqo.simulators.sympy import ensembleSimulatorSymPy
from geqo.gates.fundamental_gates import Hadamard
from geqo.gates.multi_qubit_gates import CNOT
from geqo.operations.measurement import Measure

sim = ensembleSimulatorSymPy(
    2, 2
)  # We consider a system with 2 classical bits and 2 qubits.
sim.apply(Hadamard(), [0])  # We apply the Hadamard transform to the first qubit.
sim.apply(CNOT(), [0, 1])  # We entangle both qubits.
sim.apply(
    Measure(2), [0, 1], [0, 1]
)  # We measure the 2 qubits and we store the results in the two bits.
print("first measurement", sim.ensemble)

sim.apply(Hadamard(), [0])
sim.apply(Measure(2), [0, 1], [0, 1])
print("second measurement", sim.ensemble)
first measurement {(0, 0): (1/2, Matrix([
[1, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]])), (1, 1): (1/2, Matrix([
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 1]]))}
second measurement {(0, 0): (1/4, Matrix([
[1, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]])), (1, 0): (1/4, Matrix([
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 0]])), (0, 1): (1/4, Matrix([
[0, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]])), (1, 1): (1/4, Matrix([
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 1]]))}

The result is a dictionary with the keys 00, 01, 10 and 11, which are the possible measurement outcomes of this quantum circuit. The probabilities are 1/4 for each case and the corresponding density matrices are also in the dictionary.

The operation SetBits can be used to reset the classical bits, leading to mixing both cases together.

from geqo.initialization.state import SetBits

sim.setValue("a", [0, 1])  # Set the name a to a list of bit values.
sim.apply(SetBits("a", 2), [], [0, 1])  # Set the classical bits in the simulator.
sim.ensemble
{(0, 1): (1,
  Matrix([
  [1/4,   0,   0,   0],
  [  0, 1/4,   0,   0],
  [  0,   0, 1/4,   0],
  [  0,   0,   0, 1/4]]))}

The operation SetBits sets all bits to a fixed value. As a consequence, the measurement results are lost and the states corresponding to the measurement results turn into a mixed state.

2.6.2. ensembleSimulatorCuPy#

This is the CuPy version of the ensembleSimulatorSymPy. It has almost identical functionalities except that it’s much more efficient for numerical simulations.

The following example shows the result of the measurement of a GHZ state followed by a Hadamard operation and a second measurement.

from geqo.simulators.cupy import ensembleSimulatorCuPy
from geqo.gates.fundamental_gates import Hadamard
from geqo.gates.multi_qubit_gates import CNOT
from geqo.operations.measurement import Measure

n = 10  # 10-qubit GHZ state
sim = ensembleSimulatorCuPy(n, n)
sim.apply(Hadamard(), [0])  # We apply the Hadamard transform to the first qubit.
for i in range(n - 1):
    sim.apply(CNOT(), [i, i + 1])  # We entangle all the qubits.
sim.apply(Measure(n), [*range(n)], [*range(n)])
print("first measurement", sim.ensemble)

sim.apply(Hadamard(), [0])
sim.apply(Measure(n), [*range(n)], [*range(n)])
print("second measurement", sim.ensemble)
first measurement {(0, 0, 0, 0, 0, 0, 0, 0, 0, 0): (np.float64(0.4999999999999999), array([[1.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       ...,
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]],
      shape=(1024, 1024))), (1, 1, 1, 1, 1, 1, 1, 1, 1, 1): (np.float64(0.4999999999999999), array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       ...,
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 1.+0.j]],
      shape=(1024, 1024)))}
second measurement {(0, 0, 0, 0, 0, 0, 0, 0, 0, 0): (np.float64(0.2499999999999999), array([[1.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       ...,
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]],
      shape=(1024, 1024))), (1, 0, 0, 0, 0, 0, 0, 0, 0, 0): (np.float64(0.2499999999999999), array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       ...,
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]],
      shape=(1024, 1024))), (0, 1, 1, 1, 1, 1, 1, 1, 1, 1): (np.float64(0.2499999999999999), array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       ...,
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]],
      shape=(1024, 1024))), (1, 1, 1, 1, 1, 1, 1, 1, 1, 1): (np.float64(0.2499999999999999), array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       ...,
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 1.+0.j]],
      shape=(1024, 1024)))}

2.6.3. mixedStateSimulatorSymPy#

This simulator is similar to ensembleSimulatorSymPy, but the output of measurement results and states is more compact. Instead of tracking every possible measuremnt outcome, all branches are compactified into a single mixed-state density matrix after the measurement. Hence, it is also more efficient.

from geqo.simulators.sympy import mixedStateSimulatorSymPy
from geqo.gates.fundamental_gates import Hadamard
from geqo.gates.multi_qubit_gates import CNOT
from geqo.operations.measurement import Measure

sim = mixedStateSimulatorSymPy(
    2, 2, return_density=True
)  # We consider a system with 2 classical bits and 2 qubits.
sim.apply(Hadamard(), [0])  # We apply the Hadamard transform to the first qubit.
sim.apply(CNOT(), [0, 1])  # We entangle both qubits.
sim.apply(
    Measure(2), [0, 1], [0, 1]
)  # We measure the 2 qubits and we store the results in the two bits.

sim.apply(Hadamard(), [0])
sim.apply(Measure(2), [0, 1], [0, 1])

print("first measurement", sim.measureHistory[0])
print("second measurement", sim.measureHistory[1])
first measurement {(0, 0): 1/2, (1, 1): 1/2, 'mixed_state': Matrix([
[1/2, 0, 0,   0],
[  0, 0, 0,   0],
[  0, 0, 0,   0],
[  0, 0, 0, 1/2]])}
second measurement {(0, 0): 1/4, (0, 1): 1/4, (1, 0): 1/4, (1, 1): 1/4, 'mixed_state': Matrix([
[1/4,   0,   0,   0],
[  0, 1/4,   0,   0],
[  0,   0, 1/4,   0],
[  0,   0,   0, 1/4]])}

2.6.4. mixedStateSimulatorCuPy#

This is the CuPy version of the mixedStateSimulatorSymPy. It has almost identical functionalities except that it’s much more efficient for numerical simulations.

from geqo.simulators.cupy import mixedStateSimulatorCuPy
from geqo.gates.fundamental_gates import Hadamard
from geqo.gates.multi_qubit_gates import CNOT
from geqo.operations.measurement import Measure

n = 10  # 10-qubit GHZ state
sim = mixedStateSimulatorCuPy(n, n, return_density=True)
sim.apply(Hadamard(), [0])  # We apply the Hadamard transform to the first qubit.
for i in range(n - 1):
    sim.apply(CNOT(), [i, i + 1])  # We entangle all the qubits.
sim.apply(Measure(n), [*range(n)], [*range(n)])

sim.apply(Hadamard(), [0])
sim.apply(Measure(n), [*range(n)], [*range(n)])

print("first measurement", sim.measureHistory[0])
print("second measurement", sim.measureHistory[1])
first measurement {(0, 0, 0, 0, 0, 0, 0, 0, 0, 0): np.float64(0.4999999999999999), (1, 1, 1, 1, 1, 1, 1, 1, 1, 1): np.float64(0.4999999999999999), 'mixed_state': array([[0.5+0.j, 0. +0.j, 0. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
       [0. +0.j, 0. +0.j, 0. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
       [0. +0.j, 0. +0.j, 0. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
       ...,
       [0. +0.j, 0. +0.j, 0. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
       [0. +0.j, 0. +0.j, 0. +0.j, ..., 0. +0.j, 0. +0.j, 0. +0.j],
       [0. +0.j, 0. +0.j, 0. +0.j, ..., 0. +0.j, 0. +0.j, 0.5+0.j]],
      shape=(1024, 1024))}
second measurement {(0, 0, 0, 0, 0, 0, 0, 0, 0, 0): np.float64(0.24999999999999992), (0, 1, 1, 1, 1, 1, 1, 1, 1, 1): np.float64(0.24999999999999992), (1, 0, 0, 0, 0, 0, 0, 0, 0, 0): np.float64(0.24999999999999992), (1, 1, 1, 1, 1, 1, 1, 1, 1, 1): np.float64(0.24999999999999992), 'mixed_state': array([[0.25+0.j, 0.  +0.j, 0.  +0.j, ..., 0.  +0.j, 0.  +0.j, 0.  +0.j],
       [0.  +0.j, 0.  +0.j, 0.  +0.j, ..., 0.  +0.j, 0.  +0.j, 0.  +0.j],
       [0.  +0.j, 0.  +0.j, 0.  +0.j, ..., 0.  +0.j, 0.  +0.j, 0.  +0.j],
       ...,
       [0.  +0.j, 0.  +0.j, 0.  +0.j, ..., 0.  +0.j, 0.  +0.j, 0.  +0.j],
       [0.  +0.j, 0.  +0.j, 0.  +0.j, ..., 0.  +0.j, 0.  +0.j, 0.  +0.j],
       [0.  +0.j, 0.  +0.j, 0.  +0.j, ..., 0.  +0.j, 0.  +0.j, 0.25+0.j]],
      shape=(1024, 1024))}

2.6.5. simulatorStatevectorNumpy#

This simulator is based on NumPy and can handle numeric values only. It simulates the state vector evolution of the system. The two parameters are the number of classical bits and of qubits.

The state vector of the system can be accessed by the member variable state. In the following example, the state vector of the simulator is shown after the application of a single Hadamard gate.

from geqo.simulators.numpy import simulatorStatevectorNumpy
from geqo.gates.fundamental_gates import Hadamard

gate = Hadamard()

sim = simulatorStatevectorNumpy(
    1, 0
)  # Create a state vector simulator with 0 classical bits and 1 qubit.

sim.apply(gate, [0])  # Apply the Hadamard gate on qubit 0.
sim.state  # get the state vector
array([[0.70710678+0.j],
       [0.70710678+0.j]])

The simulator provides access to all possible measurement results with their respective probabilities. These values can be accessed by the member variable measurementResult. The result is a dictionary with the possible results as keys and the probabilities as values.

This simulator does currently not support the simulation of individual measurements including the subsequent wave function reduction. Therefore, it is not possible to apply further gates after a measurement using this simulator.

The following example shows the measurement results and their probabilities after applying a Hadamard gate and measuring in the standard basis.

from geqo.simulators.numpy import simulatorStatevectorNumpy
from geqo.gates.fundamental_gates import Hadamard
from geqo.operations.measurement import Measure

gate = Hadamard()
meas = Measure(1)

sim = simulatorStatevectorNumpy(1, 1)

sim.apply(gate, [0])
sim.apply(meas, [0], [0])

res = sim.measurementResult
for r in res:
    print("result", r, "measured with probability", res[r])
result (0,) measured with probability 0.4999999999999999
result (1,) measured with probability 0.4999999999999999

2.6.6. statevectorSimulatorCuPy#

This is the CuPy version of the simulatorStatevectorNumpy. It has almost identical functionalities except that it’s much more efficient.

from geqo.simulators.cupy import statevectorSimulatorCuPy
from geqo.gates.fundamental_gates import Hadamard

gate = Hadamard()

sim = statevectorSimulatorCuPy(
    10, 0
)  # Create a state vector simulator with 0 classical bits and 10 qubit.

for i in range(10):
    sim.apply(gate, [i])  # Apply the Hadamard gate on all the qubits.
sim.state  # get the state vector
array([[0.03125+0.j],
       [0.03125+0.j],
       [0.03125+0.j],
       ...,
       [0.03125+0.j],
       [0.03125+0.j],
       [0.03125+0.j]], shape=(1024, 1))
from geqo.simulators.cupy import statevectorSimulatorCuPy
from geqo.gates.fundamental_gates import Hadamard
from geqo.operations.measurement import Measure

gate = Hadamard()
meas = Measure(10)

sim = statevectorSimulatorCuPy(10, 10)

sim.apply(gate, [0])
sim.apply(meas, [*range(10)], [*range(10)])

res = sim.measurementResult
for r in res:
    if res[r] > 1e-3:
        print("result", r, "measured with probability", res[r])
result (0, 0, 0, 0, 0, 0, 0, 0, 0, 0) measured with probability 0.4999999999999999
result (1, 0, 0, 0, 0, 0, 0, 0, 0, 0) measured with probability 0.4999999999999999

2.6.7. simulatorUnitarySymPy#

This simulator is specifically designed for obtaining the unitary matrix corresponding to a gate or a sequence of gates. It only supports unitary gates. After applying all gates, the corresponding unitary can be accessed via the u member variable.

from geqo.simulators.sympy import simulatorUnitarySymPy
from geqo.gates.fundamental_gates import PauliX

gate = PauliX()
sim = simulatorUnitarySymPy(1)
sim.apply(gate, [0])

sim.u
\[\begin{split}\displaystyle \left[\begin{matrix}0 & 1\\1 & 0\end{matrix}\right]\end{split}\]

2.6.8. unitarySimulatorCuPy#

This is the CuPy version of the simulatorUnitarySymPy. It has almost identical functionalities except that it’s much more efficient for numerical simulations.

from geqo.simulators.cupy import unitarySimulatorCuPy
from geqo.gates.fundamental_gates import PauliX
from geqo.algorithms import QFT

gate = QFT(10)
sim = unitarySimulatorCuPy(10)
sim.apply(gate, [*range(10)])

sim.u
array([[0.03125   +0.j        , 0.03125   +0.j        ,
        0.03125   +0.j        , ..., 0.03125   +0.j        ,
        0.03125   +0.j        , 0.03125   +0.j        ],
       [0.03125   +0.j        , 0.03124941+0.00019175j,
        0.03124765+0.00038349j, ..., 0.03124471-0.00057521j,
        0.03124765-0.00038349j, 0.03124941-0.00019175j],
       [0.03125   +0.j        , 0.03124765+0.00038349j,
        0.03124059+0.00076691j, ..., 0.03122882-0.00115023j,
        0.03124059-0.00076691j, 0.03124765-0.00038349j],
       ...,
       [0.03125   +0.j        , 0.03124471-0.00057521j,
        0.03122882-0.00115023j, ..., 0.03120236+0.00172485j,
        0.03122882+0.00115023j, 0.03124471+0.00057521j],
       [0.03125   +0.j        , 0.03124765-0.00038349j,
        0.03124059-0.00076691j, ..., 0.03122882+0.00115023j,
        0.03124059+0.00076691j, 0.03124765+0.00038349j],
       [0.03125   +0.j        , 0.03124941-0.00019175j,
        0.03124765-0.00038349j, ..., 0.03124471+0.00057521j,
        0.03124765+0.00038349j, 0.03124941+0.00019175j]],
      shape=(1024, 1024))