Primitive estimator with shots number fails with complex valued averages

[AlexandreFoley]

Environment

• Qiskit version: 1.0.1
• Python version: 3.12.2
• Operating system: Mac OS 14.3.1

What is happening?

The qiskit.primitives.Estimator fails to compute the average when supplied with an operator with non-real average and a shots number, with the following error:

Cannot cast scalar from dtype('complex128') to dtype('float64') according to the rule 'safe' ComplexWarning: Casting complex values to real discards the imaginary part expectation_value_with_error = rng.normal(expectation_value, standard_error)

The problem is due to the way the noise is added onto the exact expectation, it cannot deal with complex values. The random number generator is unsuitable for that case, and the actual variance is not computed. In fact, the current implementation computes the pseudovariance in the complex case.

How can we reproduce the issue?

from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import Estimator
from qiskit import QuantumCircuit

A = SparsePauliOp('XX')+1j*SparsePauliOp('ZZ')
circuit = QuantumCircuit(2)
circuit.h(0)
circuit.cx(0,1)
estimator = Estimator(options={'shots':4000}) 
aver = estimator.run(circuit,A).result().values[0]

What should happen?

A noisy complex valued average should have been computed.

Any suggestions?

Treating the complex and imaginary part as the average of two correlated random variables is likely the most rigorous thing to do. This would require computing the covariance matrix, and generating the real and imaginary part from a bivariate normal distribution.

But treating them as independent random variables might be good enough. In that case, one would compute the variance of the real and imaginary part, and sample a normal distribution with the relevant variance for each parts.

[ikkoham]

Estimator does not support non Hermitian operator. It would be nice to add more validation.

[AlexandreFoley]

The source code for primitives.Estimator suggest otherwise: it attempts to treat complex valued averages, just in a way that is defective.

If truly estimator must not compute averages of non-hermitian operators, then this must be enforced as a pre-condition on the arguments, and documented.

There's little reasons not to accept non-hermitian operators in estimator because whether Hermitian or not, the operators has to be decomposed on a basis of observables that the backend can measure: The supplied operator is never measured directly and therefore does not need to be an observable itself. And simulators are not limited by such physical constraints.