# Quantum Phase Estimation answers distribution

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? 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) Feb 16 at 23:37
• 5.039... is correct, I was referring to the phase in radians. I'll try to check the distribution. 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')

• np.linalg.eig gives the eigenvectors as columns of a matrix, so you need vecs_u[:,0] instead of vecs_u[0]. Feb 19 at 8:38
• 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$. Feb 20 at 13:11