Expose seed in Estimator and StatevectorEstimator

[t-imamichi]

Summary

The Estimator implementations based on Statevector, i.e., Estimator and StatevectorEstimator may return expectation values in a stochastic way if input circuits include reset gates for a some subsystems.
This is because Statevector behaves in a stochastic way in such a situation and we need to set a random seed to make it reproducible. Here are the details.

qiskit/qiskit/quantum_info/states/statevector.py

Lines 611 to 616 in b75c2e0

	If all subsystems are reset this will return the ground state
	on all subsystems. If only a some subsystems are reset this
	function will perform a measurement on those subsystems and
	evolve the subsystems so that the collapsed post-measurement
	states are rotated to the 0-state. The RNG seed for this
	sampling can be set using the :meth:`seed` method.

This PR introduces a private utility to function to set a random number generator to make the output reprodicible.

Fixes #12519 (result with this PR #12519 (comment))

Details and comments

Based on #12862 (comment)

import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.primitives import Estimator, StatevectorEstimator
from qiskit.quantum_info import SparsePauliOp

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.reset(0)

op = SparsePauliOp("ZI")

n = 500

for seed in [123, None]:
    print("seed", seed)
    est1 = Estimator(options={"seed": seed})
    est2 = StatevectorEstimator(seed=seed)
    for _ in range(3):
        result = est1.run([qc for _ in range(n)], n * [op]).result()
        values = result.values
        print("Estimator\t\t", np.mean(values))

        result = est2.run([(qc, n * [op])]).result()
        values = result[0].data.evs
        print("StatevectorEstimator\t", np.mean(values))
    print()

main (the values are close to zero, but random)

seed 123
Estimator		 0.012
StatevectorEstimator	 0.016
Estimator		 0.008
StatevectorEstimator	 0.04
Estimator		 -0.008
StatevectorEstimator	 -0.008

seed None
Estimator		 0.012
StatevectorEstimator	 0.024
Estimator		 0.028
StatevectorEstimator	 -0.04
Estimator		 0.048
StatevectorEstimator	 -0.02

this PR (the values are close to zero and reproducible thanks to seed=123)

seed 123
Estimator		 0.04
StatevectorEstimator	 0.04
Estimator		 0.04
StatevectorEstimator	 0.04
Estimator		 0.04
StatevectorEstimator	 0.04

seed None
Estimator		 0.068
StatevectorEstimator	 0.012
Estimator		 -0.02
StatevectorEstimator	 0.06
Estimator		 0.04
StatevectorEstimator	 0.032

[qiskit-bot]

One or more of the following people are relevant to this code:

• @Qiskit/terra-core
• @ajavadia
• @levbishop
• @t-imamichi

[coveralls]

Pull Request Test Coverage Report for Build 10298480890

Details

• 8 of 8 (100.0%) changed or added relevant lines in 3 files are covered.
• 33 unchanged lines in 5 files lost coverage.
• Overall coverage decreased (-0.03%) to 89.72%

---

Files with Coverage Reduction	New Missed Lines	%
crates/qasm2/src/expr.rs	1	94.02%
qiskit/transpiler/passes/synthesis/unitary_synthesis.py	2	88.39%
crates/qasm2/src/lex.rs	5	92.23%
qiskit/synthesis/two_qubit/xx_decompose/decomposer.py	7	90.84%
crates/qasm2/src/parse.rs	18	96.69%

Totals
Change from base Build 10294759406:	-0.03%
Covered Lines:	67314
Relevant Lines:	75027

---

💛 - Coveralls

[Cryoris]

tl;dr: It seems good to expose the seed, but I think the user has the wrong expectation for what this class does.

---

While I would agree that we can expose the seed, I don't think it fixes the expectation the user describes in the issue. By setting the seed we're effectively fixing the state of the mixed state and averaging over multiple executions won't give the correct result. For example, creating a Bell state and then resetting qubit 1, we would expect qubit 2 to be in a mixture of $|0\rangle$ and $|1\rangle$, but if we set the seed we fix this outcome. E.g.

import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.primitives import Estimator
from qiskit.quantum_info import SparsePauliOp

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.reset(0)

op = SparsePauliOp("ZI")

n = 100
est = Estimator(options={"seed": 1})  
values = est.run([qc for _ in range(n)], n * [op]).result().values
print(np.mean(values))  # should be ~0, but with the seed it is either +1 or -1

If the user wants correct behavior they should probably average over multiple runs or use another primitive?

[SillieWous]

(Commenting as the user who filed the issue)
@Cryoris From Estimator docs

shots (None or int) – The number of shots. If None, it calculates the exact expectation values

Since default is None, and I do not set it, I do not see how 'exact' can result in stochastic behaviour.

[t-imamichi]

@Cryoris I see. Using seed to make Statevector could be confusing for your example. It might be nice to provide a separate option to make Statevector deterministic.

Updated: the next commit seems to work fine.

[t-imamichi]

I updated this PR to make the result of this case #12862 (comment) stochastic but reproducible by fixing a random seed.

output of this PR

0.06

It also works #12519 case too if we set seed.
I noticed that the case #12519 requires reproducible results.
If results are reproducible, the stochastic results due to reset works fine for the case.
Here is a plot #12519 (comment)

I updated the document of Estimator. As #12862 (comment) mentioned, "exact expectation values" is not an accurate description with resets of subsystems. I copied a phrase from StatevectorEstimator.

If None, it calculates the exact expectation values

[t-imamichi]

updated description #12862 (comment)

[Cryoris]

Since default is None, and I do not set it, I do not see how 'exact' can result in stochastic behaviour.

tl;dr: The statevector class cannot represent classical mixtures (states that require a density matrix representation) and randomly picks one of the valid statevectors.

---

The statevector simulation is only that -- simulation of statevectors by means of unitary evolution 🙂
Adding non-unitary operations such as resets results in mixed states, which cannot generally be represented by a statevector but would need a density matrix representation. In your case, for example, we can write the reset operation $\mathcal{R}$ using Kraus operators as

$$\mathcal{R}: \rho \mapsto K_1\rho K_1^\dagger + K_2\rho K_2^\dagger$$

with

$$K_1 = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix},~ K_2 = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}.$$

If we start from a Bell state $|\psi\rangle = (|00\rangle + |11\rangle)/\sqrt{2}$ and apply a reset on the second qubit, we perform the operation

$$|\psi\rangle\langle\psi| \mapsto \rho = (\mathbb 1 \otimes \mathcal{R})(|\psi\rangle\langle\psi|) = \frac{1}{2} \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

which is mixed ($Tr(\rho^2) = 1/2$) and can be understood as a equiprobable classical mixture of the two pure states $|00\rangle$ and $|01\rangle$. Since the Statevector class cannot represent that, it picks either of the two possible states, which makes the "exact" statevector calculation stochastic.

So to get an exact solution to your circuits that include resets, you might want to use a density matrix simulator (Aer should have one). Alternatively you can run a large number of statevector simulations, and average over the results, which equals a manual average over the classical mixture.

[SillieWous]

Thank you for taking the time to explain it. I will be reading up on the Kraus operators.

Let me start by saying that I am and was not questioning your expertise. Quite the opposite, you definitely understand the code and physics much better than me. But due to that expertise your expectations are based on the implementation, whereas user expectations come from the documentation.

Therefore when the documentation states 'exact' it is not strange that the user thinks it is deterministic. Apparently that is not possible due to the way it is implemented (which probably make a lot of sense), but in my opinion it is not fair to assume the user knows implementation details.

As someone with more of a background in EE and software, I see the reset gate as just a (non-unitary) matrix multiplication. My knowledge reaches as far as that non-unitary operations are forbidden in a real circuit, but running locally is not a real circuit. So when I read 'exact', that makes me think no such limitations apply (running locally).

PS
I'm sorry if I come across as rude or something, that is not my intention. I just get triggered slightly at things that come even remotely close to 'blaming the user'.

[t-imamichi]

Thank you for the detailed explanation, @Cryoris. Yes, since the statevector class cannot represent the classical mixture, I suggest rewording "exact" in the documentation in this PR.

As for the issue #12519, the script expects the same statevector for all observables and this PR helps do it. Of course, they still need to conduct the experiments multiple times with different seeds to deal with the stochastic statevector behavior due to reset for subsystems.

[Cryoris]

If this runs very fast, could we increase n (maybe 1000)? If I did the math right this has still a ~1% chance of failing, which could be quite a lot in the long run. If it is not very fast I think we can also skip this test, as the behavior of reset should be tested in the Statevector class.

[t-imamichi]

n = 1000 takes around 0.3 sec on mac M1 (while n = 250 takes around 0.09 sec).
I will increase the number to n = 1000.

[Cryoris]

We discussed with the team and the failure likelihood (now ~0.1%) is still too high as it would likely fail the CI once a day. Could you move that test to test/randomized? Then CI will not be blocked on the failures but we can still get a notice 🙂

[t-imamichi]

Because I fixed seed, it does not fail as far as there is no breaking change in np.random or Statevector.reset.

I printed result.values.mean() several times and the values are consistent.

(dev-qiskit) ima@mbp3 ~/t/4/q/terra (primitives/estimator-reset)> python -m unittest test/python/primitives/test_estimator.py test/python/primitives/test_statevector_estimator.py
.....0.006
0.032
........................0.006
0.03223976655569642
......
----------------------------------------------------------------------
Ran 35 tests in 1.942s

OK
(dev-qiskit) ima@mbp3 ~/t/4/q/terra (primitives/estimator-reset)> python -m unittest test/python/primitives/test_estimator.py test/python/primitives/test_statevector_estimator.py
.....0.006
0.032
........................0.006
0.03223976655569642
......
----------------------------------------------------------------------
Ran 35 tests in 1.838s

OK
(dev-qiskit) ima@mbp3 ~/t/4/q/terra (primitives/estimator-reset)> python -m unittest test/python/primitives/test_estimator.py test/python/primitives/test_statevector_estimator.py
.....0.006
0.032
........................0.006
0.03223976655569642
......
----------------------------------------------------------------------
Ran 35 tests in 1.834s

OK

diff to display mean values.

diff --git a/test/python/primitives/test_estimator.py b/test/python/primitives/test_estimator.py
index 2c251af65..f17b0a7d7 100644
--- a/test/python/primitives/test_estimator.py
+++ b/test/python/primitives/test_estimator.py
@@ -370,7 +370,8 @@ class TestEstimator(QiskitTestCase):
                 estimator = Estimator(options={"seed": seed})
                 result = estimator.run([qc for _ in range(n)], [op] * n).result()
             # expectation values should be stochastic due to reset for subsystems
-            np.testing.assert_allclose(result.values.mean(), 0, atol=1e-1)
+            np.testing.assert_allclose(z:=result.values.mean(), 0, atol=1e-1)
+            print(z)

             with self.assertWarns(DeprecationWarning):
                 result2 = estimator.run([qc for _ in range(n)], [op] * n).result()
@@ -383,7 +384,8 @@ class TestEstimator(QiskitTestCase):
                 estimator = Estimator(options={"seed": seed})
                 result = estimator.run([qc for _ in range(n)], [op] * n, shots=shots).result()
             # expectation values should be stochastic due to reset for subsystems
-            np.testing.assert_allclose(result.values.mean(), 0, atol=1e-1)
+            np.testing.assert_allclose(z:=result.values.mean(), 0, atol=1e-1)
+            print(z)

             with self.assertWarns(DeprecationWarning):
                 result2 = estimator.run([qc for _ in range(n)], [op] * n, shots=shots).result()
diff --git a/test/python/primitives/test_statevector_estimator.py b/test/python/primitives/test_statevector_estimator.py
index 0cd76550f..2913e8c53 100644
--- a/test/python/primitives/test_statevector_estimator.py
+++ b/test/python/primitives/test_statevector_estimator.py
@@ -321,7 +321,8 @@ class TestStatevectorEstimator(QiskitTestCase):
         with self.subTest("precision=0"):
             result = estimator.run([(qc, [op] * n)]).result()
             # expectation values should be stochastic due to reset for subsystems
-            np.testing.assert_allclose(result[0].data.evs.mean(), 0, atol=1e-1)
+            np.testing.assert_allclose(z:=result[0].data.evs.mean(), 0, atol=1e-1)
+            print(z.item())

             result2 = estimator.run([(qc, [op] * n)]).result()
             # expectation values should be reproducible due to seed
@@ -331,7 +332,8 @@ class TestStatevectorEstimator(QiskitTestCase):
             precision = 0.01
             result = estimator.run([(qc, [op] * n)], precision=precision).result()
             # expectation values should be stochastic due to reset for subsystems
-            np.testing.assert_allclose(result[0].data.evs.mean(), 0, atol=1e-1)
+            np.testing.assert_allclose(z:=result[0].data.evs.mean(), 0, atol=1e-1)
+            print(z.item())

             result2 = estimator.run([(qc, [op] * n)], precision=precision).result()
             # expectation values should be reproducible due to seed

[Cryoris]

@SillieWous no worries, happy to discuss 🙂

@t-imamichi maybe we could add a note to the StatevectorEstimator docstring that explaining that it is exact if the circuit contains only unitary operations?