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

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?

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)

• Thank you for you proposal! However, try to make a fixed array of angles and fix the seed of random_unitary, and then run the circuit many times. The values of the statevector change, so we have no exact probabilities. Jan 30 at 9:32
• @dviqu You are right.. And this is weird! Setting a fixed array of angles and a fixed unitary, the statevector evolution should be deterministic and lead exactly to the same probabilities if you repeat the full experiment multiple times. The problem looks to be hidden in how the reset operation works in Qiskit: it seems that at each qc.reset(0) call the statevector evolution splits in 2 possible ways, so that after $N$ repetitions of your block $U R_y(x_i)$ the statevector has $2^N$ possible configurations... Jan 30 at 11:12
• I figured out what the problem with the reset operation is and edited my original answer with all the details Jan 30 at 13:53
• You are right! The reset is implemented as a measurement plus a c_if, and this leads to two possible statevectors after each reset. Jan 31 at 8:45

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

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)]))

• Thank you! Your proposal is useful for my case because it separates the results from each measurement. However, try to run your code more than once and compare the outputs. If you look at them carefully, apart from the first, they are slightly different. As we have intermediate measurements, the system collapses to one state or another depending on whether the measurement outcome was 0 or 1. I think that Qiskit actually does a lot of shots and then computes the probabilities based on these results, but there are still some statistical fluctuations and it is not computing the exact probability. Jan 30 at 9:23
• PD. To run the code several times, you have to fix the x array and the seed for unitary to have reproducibility. Jan 30 at 9:34
• The code is just to demonstrate the solution. I used randomness in two places: (1) generating rotation angles, (2) generating 2-qubit unitaries. If you replace these two things with yours, you should get same probabilities after each run. Also, I modified my answer to pass a seed to random functions to make sure the same values will be returned each time. Jan 30 at 12:13