Suppose I have a random unitary matrix, known eigenvectors and eigenvalues. I know that exact eigenvalue for the given matrix is $$0.5491617699847768+0.835716070437315j$$. From here, if I'm not mistaken the phase should be equal to $$5.039153255477287$$ - something between $$5$$ and $$6$$. I apply the QPE algorithm for $$5$$ qubit scheme for eigenvalue estimation. The distribuiton of the estimated phases is in the image below.

Question: why the distribution is not concentrated around 5 and 6 only? How to interpret this distibution? If I take the math expectation (equals to $$13.380000000000003$$), devide it on $$2^n$$, where $$n$$ in the number of qubits for eigenvalue etimation, add $$2\pi i$$ and exponentiate, again I get the value (equals to $$-0.87057047+0.49204375j$$) not close to the real eigenvalue

UPD

• added matrix and the eigenvector;
• compressed scheme to 3 qubits for eigenvalue estimation;
• added script for the circuit execution

• This does not look like a correct distribution for phase estimation algorithm. The phase for your value is actually around $0.9894$, or $0.15747\cdot 2\pi$. Are you sure your eigenvector is correct? Is there noise in your computation? Commented Feb 16 at 20:47
• In the question I used my own implementation of QPE. Now I checked using qiskit and got the same distribution as a result. By the way, how did you get $0.9894$? If the eigenvalue is $eig = e^{2\cdot \pi \cdot i \cdot \frac{p}{2^n}}$, where $p$ is a phase and $n$ - number of qubits for the eigenvalue estimation, then the phase equals to $p = \frac{\ln{eig}\cdot 2^n}{2\cdot \pi i} = 5.039153255477287$ (updated value in the question) Commented Feb 16 at 23:37
• 5.039... is correct, I was referring to the phase in radians. I'll try to check the distribution. Commented Feb 16 at 23:42

Doublecheck your inverse QFT circuit. I get a distribution similar to yours if I put QFT instead of its inverse in the second part of the phase estimation algorithm. If the correct inverse QFT is used, the distribution is very concentrated around 5.

Here is Qiskit code:

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from numpy import pi

qreg_phi = QuantumRegister(1, 'phi')
qreg_q = QuantumRegister(5, 'q')
creg_c = ClassicalRegister(5, 'c')
circuit = QuantumCircuit(qreg_phi, qreg_q, creg_c)

theta = 5.039153255477287 * pi / 16

# Prepare |1> in qreg_phi
# |1> is the eigenvector of p(theta)
circuit.x(qreg_phi[0])

# Start phase estimation for U = p(theta)

circuit.h(qreg_q[0])
circuit.h(qreg_q[1])
circuit.h(qreg_q[2])
circuit.h(qreg_q[3])
circuit.h(qreg_q[4])

# Apply controlled U's
circuit.cp(theta, qreg_q[0], qreg_phi[0])
circuit.cp(2 * theta, qreg_q[1], qreg_phi[0])
circuit.cp(4 * theta, qreg_q[2], qreg_phi[0])
circuit.cp(8 * theta, qreg_q[3], qreg_phi[0])
circuit.cp(16 * theta, qreg_q[4], qreg_phi[0])

# Uncompute QFT
circuit.h(qreg_q[4])

circuit.cp(-pi / 2, qreg_q[4], qreg_q[3])
circuit.h(qreg_q[3])

circuit.cp(-pi / 4, qreg_q[4], qreg_q[2])
circuit.cp(-pi / 2, qreg_q[3], qreg_q[2])
circuit.h(qreg_q[2])

circuit.cp(-pi / 8, qreg_q[4], qreg_q[1])
circuit.cp(-pi / 4, qreg_q[3], qreg_q[1])
circuit.cp(-pi / 2, qreg_q[2], qreg_q[1])
circuit.h(qreg_q[1])

circuit.cp(-pi / 16, qreg_q[4], qreg_q[0])
circuit.cp(-pi / 8, qreg_q[3], qreg_q[0])
circuit.cp(-pi / 4, qreg_q[2], qreg_q[0])
circuit.cp(-pi / 2, qreg_q[1], qreg_q[0])
circuit.h(qreg_q[0])

# Swap qubits
circuit.swap(qreg_q[0], qreg_q[4])
circuit.swap(qreg_q[1], qreg_q[3])

circuit.reset(qreg_phi)
circuit.measure(qreg_q, creg_c)


UPD. Here is an example with arbitrary number of qubits

def QPE(U, phi, n):
p = U.num_qubits
qreg_phi = QuantumRegister(p, 'phi')
qreg_q = QuantumRegister(n, 'q')

circuit = QuantumCircuit(qreg_phi, qreg_q)

circuit.initialize(phi, qreg_phi)

circuit.h(qreg_q)

for i in range(n):
circuit.append(U.power(2**i).control(1), [qreg_q[i], qreg_phi])

for i in range(n):
for j in range(i):
circuit.cp(-pi / 2**(i - j), qreg_q[n - 1 - j], qreg_q[n - 1 - i])
circuit.h(qreg_q[n - i - 1])

circuit.swap(qreg_q[:n//2], qreg_q[-1:-(n//2)-1:-1])

return circuit

theta = 5.039153255477287 * pi / 16
U = PhaseGate(theta)

n = 5

circuit = QPE(U, [0.0, 1.0], n)
qreg_q = circuit.qregs[1]
creg_x = ClassicalRegister(n, 'x')
circuit.measure(qreg_q, creg_x)

• Yep, answer from this script gives more intuitive results. I added script that I used in my approach (from yours one), egenvector and entry matrix. The script still does not provide me answer close to what I see using your approach. Commented Feb 18 at 23:43
• As far as I see, your vector is not an eigenvector. Commented Feb 19 at 7:33
• np.linalg.eig gives the eigenvectors as columns of a matrix, so you need vecs_u[:,0] instead of vecs_u[0]. Commented Feb 19 at 8:38
• I used phase rotation with theta*16, which is the same as 16-fold repetition of phase rotation with theta. So in my code it also goes 1,2,4,8,16. I added an example of a general construction. Commented Feb 20 at 13:07
• Note that in your current code you are basically estimating phase of $U^4$ and find that it's around $5/4 \pi$. If you run it for $U$, you'll get the distriibution around $1$, which means the phase is around $1/4 \pi$. If you increase to 5 bits, you get a better estimate of $5/16 \pi$. Commented Feb 20 at 13:11