How to get exact measurement probabilities when having intermediate measurements with Qiskit?

Question

Suppose we have a circuit with two qubits, A and B. Both are initialized to $|0\rangle$. Over qubit A we apply a single rotation gate (e.g. $R_y$) with an angle given by $x_0$, and then we entangle the two qubits by an arbitrary unitary operator $U$. Then we measure qubit A while B is never measured. Qubit A is reset to $|0\rangle$.We repeat iteratively this process, only changing the initial rotation angle over qubit A, $x_1,...x_n$. $U$ is always the same. The full circuit is represented in the figure.

Qiskit, by default, performs the simulation by making the intermediate measurements and returns the distribution of all the possible bit strings given by the classical register in which we store all the intermediate measurements. After that, we can recover the probability of getting 0 or 1 in each individual measurement. But there seems not to exist a method for computing the probability in each individual measurement in exact form (like e.g. computing the modulus of the statevector amplitudes when the measurement is only at the end).

Is there any form of getting the exact individual probabilities $p_0,p_1,...p_n$ with Qiskit? It is also not possible when using another method on AerSimulator, like density_matrix.

Does there exist any other open source kit different from qiskit that solves this?

Answer 1

To do what you want in Qiskit, you could replace each measurement operation in your circuit by a simple call to the Statevector.probabilities(qargs) method, where the list qargs contains the indices of the subsystems you need to measure (in this case qargs=[0] for the first qubit).

Here is a complete working example (I generate random angles $x_i \in [0, 2\pi[$ and a random $4 \times 4$ unitary $U$):

import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import random_unitary, Statevector

reps = 3
angles = np.random.rand(reps) * (2*np.pi)
u = random_unitary(4)

probs = []
qc = QuantumCircuit(2)
for x in angles:
    qc.reset(0)
    qc.ry(x, 0)
    qc.unitary(u, [0, 1])
    p0_p1 = Statevector(qc).probabilities(qargs=[0]) 
    probs.append(p0_p1)

print(probs)

EDIT: If you run the code above multiple times, you will get different probabilities outcomes even though you set a fixed array angles and a fixed unitary u. This is because, in Qiskit, the reset instruction is implemented as a measurement instruction in the computational basis (making the state collapse to either $| 0 \rangle$ or $| 1 \rangle$) followed by an $X$-gate conditioned on the measure outcome (flipping the qubit back to the $| 0 \rangle$ if $1$ was measured):

Since the measurement is implemented sampling random numbers (Qiskit calls the numpy.random.default_rng internally), and the given quantum state is a two-qubits entangled state, the statevector of the full system after the measurement will depend on the actual measure outcome and state collapse, yielding different state amplitudes (and so different probabilities).

So, to run the same "deterministic" statevector evolution multiples times and get exactly the same probabilities, you should set the Qiskit random seed by the method Statevector.seed, as shown in the following code:

import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import random_unitary, Statevector

reps = 3
angles = np.random.rand(reps) * (2*np.pi)
u = random_unitary(4)

qc = QuantumCircuit(2)
statevec = Statevector(qc)
statevec.seed(value=42)

probs = []
for x in angles:
    qc.reset(0)
    qc.ry(x, 0)
    qc.unitary(u, [0, 1])
    p0_p1 = statevec.evolve(qc).probabilities(qargs=[0]) 
    probs.append(p0_p1)

print(probs)

Answer 2

Alternatively, you can use Quirk for small circuits. Especially it's chance and amplitude displays are useful for visualizing measurement probabilities.

Answer 3

You can use save_probabilities() function to save the measurement outcome probabilities anywhere in your quantum circuit when using AerSimulator.

First, let's create the circuit:

from qiskit import QuantumCircuit, transpile
from qiskit.quantum_info import random_unitary
from qiskit_aer import AerSimulator
import numpy as np

ROUNDS = 5
np.random.seed(0)
x = np.random.uniform(low=0, high=np.pi, size=(ROUNDS,))

circ = QuantumCircuit(2, ROUNDS)
for m in range(ROUNDS):
    circ.reset(0)
    circ.ry(x[m], 0)
    circ.unitary(random_unitary(4, seed = m), [0, 1], label = '  $U_' + str(m) + '$  ')
    circ.save_probabilities([0], label='prob-' + str(m)) # <== here
    circ.measure(0, m)

circ.draw('mpl', fold = -1)

Now, we simulate it and read the measurement outcome probabilities:

simulator = AerSimulator()
tr_circ = transpile(circ)
job = simulator.run(tr_circ)
sim_data = job.result().data()
for m in range(ROUNDS):
    print('Probabilities {}: {}'.format(m, sim_data['prob-' + str(m)]))