Cannot replicate results in article on pricing financial derivatives on IBM Q

I am trying to implement a circuit for searching for the largest eigenvalue and respective eigenvector of an operator, i.e. phase estimation, introduced in article Towards Pricing Financial Derivatives with an IBM Quantum Computer, page 6, figure 3(a).

The authors of the article use the phase estimation for indentification of the most important component in principal component analysis. In particular, they look for the largest eigenvalue of this matrix

$$$$\rho= \begin{pmatrix} 0.6407 & 0.3288 \\ 0.3288 & 0.3593 \end{pmatrix}$$$$

The largest eigenvalue of the matrix is 0.8576 with respective eigenvector $$(0.8347; 0.5508)^T$$.

To search for that, the authors proposed following cicuit

Note that gates denoted by -1.57 and -0.79 are $$S^\dagger$$ and $$T^\dagger$$ respecitvelly. They act on "bottom" qubit and they are controlled by "top" qubit

$$U3$$ operators implement operator $$\mathrm{e}^{2\pi i \rho}$$ and its powers.

According to the article, results after measurement should be these:

Note that pink bars are results obtained on simulator, while blue ones on IBMQX2 processor.

Based on these results, the authors concluded that the largest eigenvalue is approximated by value $$0.111_2$$ (0.875 in decimal which is close to actual value 0.8576).

I tried to replicate their results with this circuit on IBM Q:

Note: please find a code in QASM at the end of the question to see parameters of gates

I received these results on IBM Q simulator.

According to my results, the largest eigenvalue should be zero which does not make sense.

So, my questions are these:

1. Where I did a mistake in circuit implementation? My circuit and the author's seems to be the same.
2. Why the qubit $$|q_3\rangle$$ is initialized by gates $$Ry$$ and $$Rz$$ with parameter $$\theta$$ equal to 1.00 and 0.33, respectively? I would expect only gate $$Ry$$ with $$\theta = 1.1665$$ as this gate produce the eigenvector $$(0.8347; 0.5508)^T$$. However, replacement of $$Ry(1.00)$$ and $$Rz(0.33)$$ with $$Ry(1.1665)$$ does not change resulting histogram significantly.

Here is a QASM code of my circuit

OPENQASM 2.0;
include "qelib1.inc";

qreg q[4];
creg c[4];

h q[0];
h q[1];
h q[2];
ry(1.00) q[3];
rz(0.33) q[3];
cu3(1.6,-1.12,2.03) q[2],q[3];
cu3(2.23,0.51,3.65) q[1],q[3];
cu3(0.8,-4.53,-1.39) q[0],q[3];
h q[0];
cu1(-pi/2) q[0],q[1];
cu1(-pi/4) q[0],q[2];
h q[1];
cu1(-pi/2) q[1],q[2];
h q[2];
measure q[3] -> c[3];
measure q[2] -> c[2];
measure q[1] -> c[1];
measure q[0] -> c[0];

• How does the outcome change as you vary the input state of $q_3$? You are supposed to be running it assuming you don't know what the eigenvector is, so you don't produce the eigenvector, but for the sake of debugging, you certainly can use the known eigenvector to see what answer it gives. Also, why do you measure $q_3$? You shouldn't need to. – DaftWullie Jan 10 at 16:40
• @DaftWullie: I know that I do not need to measure $q_3$ because the eigenvalue is a result of inverse Fourier transform on $q_0$ to $q_2$. Authors of the article proposed to put $\mathrm{H}$ on $q_3$ to measure it and iterate. For debuging purposes I tried to use $\mathrm{U3}$ producing exactly the eigenvector, however, the results are same. Firstly, I wanted to exactly replicate the circuit in the paper and see the results but they are different as you can see in figures in the question. I cannot find where is a difference between my and authors' circuit. – Martin Vesely Jan 10 at 17:15
• The QFT is wrong, you can check it. And I have a question, what are the matrix parameters in the controlled unitary in the figure 4(in the article) quantum circuit implementation for the 4×4 matrix? – user9774 Jan 13 at 7:01
• Thanks for reaction, what is wrong in particular with inverese QFT? Gates implementing QFT are simply upside down (as in the article), hence no swap gate before inverse QFT is needed. To be sure, I tried to implement the circuit with inverse QFT in convetional way but results are same. Regarding your question, you have to zoom it at 300 % to be readable. – Martin Vesely Jan 13 at 10:14

One thing that I noticed. If cu3 gate from $$q[2]$$ to $$q[0]$$ is some $$U$$, then the cu3 from $$q[2]$$ to $$q[0]$$ should be $$U^2$$ in the phase estimation algorithm, but the comparisons of operators with the help of numpy.array showed me that it's not true here. I tried to implement by replacing cu3 part of the QASM code with the following:

cu3(1.6, -1.12, 2.03) q[2] q[3]

cu3(1.6, -1.12, 2.03) q[1] q[3]
cu3(1.6, -1.12, 2.03) q[1] q[3]

cu3(1.6, -1.12, 2.03) q[0] q[3]
cu3(1.6, -1.12, 2.03) q[0] q[3]
cu3(1.6, -1.12, 2.03) q[0] q[3]
cu3(1.6, -1.12, 2.03) q[0] q[3]


And obtained a different result via IBM's 'qasm simulator':

{'0010': 39, '0101': 13, '1110': 16, '0110': 47, '1011': 4, '1010': 8, '1000': 92, '1101': 5, '1111': 143, '1001': 22, '0011': 10, '0001': 16, '1100': 3, '0100': 8, '0000': 235, '0111': 363}



And the sum of '1111' and '0111' outcomes is maximum (here 506 from 1024 measurements) as was obtained in the paper.

• Thanks, it works now. I though that $U3$ operators acting on $q_3$ controlled by $q_2$ and $q_1$ are $U^2$ and $U^4$, respectivelly. But it seems that their parameters are wrong, or at least too rounded. – Martin Vesely Jan 13 at 22:28
• Just note that I have same results as you but they are still a little bit different from those in the paper, especially frequency of state $|0000\rangle$. But it again seems that this is a matter of rounding gate parameters. – Martin Vesely Jan 13 at 22:30
• Yes, I also think that it can be because of rounding. Also, maybe they are ignoring somewhere a global phase. One can ignore a global phase, but one shouldn't ignore the controlled global phase. – Davit Khachatryan Jan 13 at 22:33