I am implementing QAE (Quantum Amplitude Estimation), which is very similar to QPE (Quantum Phase Estimation) with a Grover Operator as the U matrix of QPE.
I want to check my results, in the outputs of the invQFT block, so I calculate the reduced density martix of the output of the circuit on the 3 iQFT qubits, and look on its diagonal to check the probabilities to each of the $2^3=8$ output combinations
On the other side, I compare it to the expected results by an analytical calculation of the expectation value of the 3 qubits to be in a state y (3 qubits state). So, as in QAE, the final state is : $$\frac{1}{\sqrt{2}}\frac{\mathrm{e}^{\mathrm{i}\theta}}{M}\sum_{y=0}^{M-1}\sum_{j=0}^{M-1}\mathrm{e}^{-\frac{2\pi\mathrm{i}j}{M}\left(y-\frac{\theta M}{\pi}\right)}|y\rangle\left|\Psi_+\right\rangle-\frac{1}{\sqrt{2}}\frac{\mathrm{e}^{-\mathrm{i}\theta}}{M}\sum_{y=0}^{M-1}\sum_{j=0}^{M-1}\mathrm{e}^{-\frac{2\pi\mathrm{i}j}{M}\left(y+\frac{\theta M}{\pi}\right)}|y\rangle\left|\Psi_-\right\rangle$$
The measurement operator is: $$M_y=\left|y_0\middle\rangle\middle\langle y_0\right|\otimes I$$
So the probability of the 3 qubits to be in y is calculated by: $$\mathrm{Pr}\left(y_0\right)=\langle\psi|M_y|\psi\rangle=\langle\psi|\left(\left|y_0\middle\rangle\middle\langle y_0\right|\otimes I\right)|\psi\rangle$$
And my final analytical result is: $$\mathrm{Prob}\left(y_0\right)=\frac{1}{2M^2}\left(\frac{1-\cos\left(\frac{2\pi M}{M}\left(y-\frac{\theta M}{\pi}\right)\right)}{1-\cos\left(\frac{2\pi}{M}\left(y-\frac{\theta M}{\pi}\right)\right)}+\frac{1-\cos\left(\frac{2\pi M}{M}\left(y+\frac{\theta M}{\pi}\right)\right)}{1-\cos\left(\frac{2\pi}{M}\left(y+\frac{\theta M}{\pi}\right)\right)}\right)$$
I am trying to figure out why my results are not exactly the same and have about a 30% difference in each of the 8 values of the probabilities.
Another metric I use is the fidelity like: $$F(\rho, \sigma)=\left[\mathrm{tr}\sqrt{\rho\sigma}\right]^2=\left(\sum_k\sqrt{p_kq_k}\right)^2=F(\mathbf{p},\mathbf{q})$$
And it is about $0.97$, which is still not close enough to 1
If anyone has an idea where are the differences coming from?
