How to estimate noisy expectation value in limit of infinite shots?

I suspect this is an obvious question, but it has proven surprisingly elusive!

Suppose I have a quantum circuit, as well as a pre-defined initial state. After performing this circuit, I evaluate an observable which is bounded from below, e.g. at 0 -- in particular, I expect this output state to be vanishingly close to 0.

However, I want to model this circuit in the presence of hardware noise. Of course I can follow the usual procedure of sampling from a distribution of circuits which models this noise, and evaluate some large number of shots.

This brings me to my question: is there some analytical or computational way I can predict the average deviation of my expectation value , in the limit of infinite shots? Naively I might have expected this deviation to be 0, but of course it cannot be since the observable is near its lower bound -- some proportion of the shots will deliver erroneous values of , and those values will nearly always be too high.

• Well, it will surely depend on the noise model you're considering! If I understand correctly, your question boils down to: if the observable $O$ is positive semi-definite and $|\psi\rangle$ is an eigenvector of $0$ associated to the value $0$, how can we estimate the mean error $|\langle\psi'|O|\psi'\rangle|^2$, with $|\psi'\rangle$ being the state after some potential error happened, right? Commented Jul 5 at 17:12
• @TristanNemoz, that's a great summary, thank you -- and yes, it will certainly depend on the model. I suppose what I'm hoping for it a method (or outline of one) for how to form an effective infinite-shot limit of a noise model. Commented Jul 5 at 22:21

Depending on your noise model, you should be able to compute the final density matrix that you obtain after your computations, let us denote it $$\rho$$. Let us assume that the state you intended to prepare was $$|\psi\rangle$$. The quantity you're looking for is then: $$\mathrm{Tr}[\rho O]-\langle\psi|O|\psi\rangle$$ Note that you could also put an absolute value or a squared one, but fundamentally, this will depend on this quantity.

Let us take a small example to understand where this is coming from. Let us assume that you apply two $$X$$ gates to $$|0\rangle$$ and then measure in the computational basis. That is, we have $$|\psi\rangle=|0\rangle$$ and $$\langle\psi|O|\psi\rangle=1$$.

The noise model we'll consider is that after each gate that we apply, a bit flip happens with probability $$p_1$$. Furthermore, when we measure, a bit-flip happens with probability $$p_2$$. You can convince your self that the final density matrix is: $$\rho=\left[2p_1\left(1-p_1\right)p_2+\left(1-p_1\right)^2\left(1-p_2\right)+p_1^2\left(1-p_2\right)\right]|0\rangle\left[2p_1\left(1-p_1\right)\left(1-p_2\right)+\left(1-p_1\right)^2p_2+p_1^2p_2\right]|1\rangle$$ You can see this as a classical incertitude over the state you're getting. With probability $$2p_1\left(1-p_1\right)p_2+\left(1-p_1\right)^2\left(1-p_2\right)+p_1^2\left(1-p_2\right)$$ you get the state you expect, and with probability $$2p_1\left(1-p_1\right)\left(1-p_2\right)+\left(1-p_1\right)^2p_2+p_1^2p_2$$ you don't. You then have to compute the expectation value over all possibilities, weighted by their respective probability. In our case, you'll find that you'll measure $$1$$ with probability $$2p_1\left(1-p_1\right)p_2+\left(1-p_1\right)^2\left(1-p_2\right)+p_1^2\left(1-p_2\right)$$ and $$-1$$ with probability $$2p_1\left(1-p_1\right)\left(1-p_2\right)+\left(1-p_1\right)^2p_2+p_1^2p_2$$, which gives you an expected value of $$2p_1\left(1-p_1\right)\left(2p_2-1\right)+\left(1-p_1\right)^2\left(1-2p_2\right)+p_1^2\left(1-2p_2\right)$$.

More generally, you can express the final state as: $$\rho=\sum_ip_i\left|\psi_i\middle\rangle\!\middle\langle\psi_i\right|$$ Which means that you will get $$\left|\psi_i\right\rangle$$ with probability $$p_i$$. Thus, the expectation value that you get is: $$\sum_ip_i\left\langle\psi_i\middle|O\middle|\psi_i\right\rangle=\sum_ip_i\mathrm{Tr}\left[O\left|\psi_i\middle\rangle\!\middle\langle\psi_i\right|\right]=\mathrm{Tr}\left[O\sum_ip_i\left|\psi_i\middle\rangle\!\middle\langle\psi_i\right|\right]=\mathrm{Tr}[O\rho].$$

You then only have to compare that with your desired expectation value.

Note that this method may not be the most convenient though, as it requires you to compute $$\rho$$, which may be tedious depending on your noise model.

I'm not sure whether that answers your question, if it doesn't tell me and I'll happily delete it!