Implementing Quantum Phase Estimation (QPE) for Eigenvalue Evaluation in Multi-Qubit Systems
Quantum Phase Estimation is a circuit construction used to estimate the phase associated with an eigenvector of a unitary operator. It is a core primitive in quantum algorithms for eigenvalue estimation, Hamiltonian simulation analysis, and spectral decomposition. This document defines a formal implementation pattern for constructing QPE circuits in a Python-based Qiskit-compatible environment.
Let \(U\) be a unitary operator and let \(|\psi\rangle\) be an eigenvector of \(U\). The eigenvalue relation is:
\[ U|\psi\rangle = e^{2\pi i \phi}|\psi\rangle \]
The purpose of QPE is to estimate \(\phi \in [0,1)\). If the evaluation register contains \(m\) qubits, the estimated phase is represented on a grid with spacing:
\[ \Delta \phi = 2^{-m} \]
The algorithm uses two registers. The first register stores the binary phase estimate. The second register stores the eigenstate of the target unitary. The evaluation register is placed into uniform superposition, controlled powers of \(U\) are applied, and an inverse Quantum Fourier Transform maps relative phase information into computational-basis amplitudes.
Immediately before the inverse Fourier transform, the ideal state is:
\[ \frac{1}{2^{m/2}} \sum_{k=0}^{2^m-1} e^{2\pi i k\phi}|k\rangle|\psi\rangle \]
When \(\phi\) has an exact \(m\)-bit binary representation, the measurement result is deterministic. If the observed integer is \(y\), the decoded phase is:
\[ \hat{\phi}=\frac{y}{2^m} \]
The corresponding unitary eigenvalue is reconstructed as:
\[ \hat{\lambda}=e^{2\pi i\hat{\phi}} \]
Reference Implementation
from __future__ import annotations
import math
from typing import Dict
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit.circuit.library import QFT
from qiskit_aer import AerSimulator
def controlled_phase_power(theta: float, power: int) -> QuantumCircuit:
"""
Return a controlled-U^power circuit for:
U = diag(1, exp(2*pi*i*theta))
"""
circuit = QuantumCircuit(2, name=f"cU^{power}")
circuit.cp(2.0 * math.pi * theta * power, 0, 1)
return circuit
def build_qpe_circuit(
phase: float,
evaluation_qubits: int
) -> QuantumCircuit:
evaluation = QuantumRegister(evaluation_qubits, "eval")
system = QuantumRegister(1, "sys")
classical = ClassicalRegister(evaluation_qubits, "c")
circuit = QuantumCircuit(evaluation, system, classical)
# Prepare eigenstate |1> of the phase unitary.
circuit.x(system[0])
# Initialize the evaluation register.
for qubit in evaluation:
circuit.h(qubit)
# Apply controlled powers U^(2^k).
for k in range(evaluation_qubits):
power = 2 ** k
gate = controlled_phase_power(phase, power).to_gate()
circuit.append(gate, [evaluation[k], system[0]])
inverse_qft = QFT(
num_qubits=evaluation_qubits,
inverse=True,
do_swaps=True
).to_gate(label="IQFT")
circuit.append(inverse_qft, list(evaluation))
circuit.measure(evaluation, classical)
return circuit
def decode_phase(counts: Dict[str, int], evaluation_qubits: int) -> float:
dominant_bitstring = max(counts, key=counts.get)
integer_value = int(dominant_bitstring, 2)
return integer_value / (2 ** evaluation_qubits)
def run_reference_qpe() -> None:
exact_phase = 0.3125
evaluation_qubits = 5
shots = 4096
circuit = build_qpe_circuit(
phase=exact_phase,
evaluation_qubits=evaluation_qubits
)
simulator = AerSimulator()
compiled = transpile(circuit, simulator, optimization_level=2)
result = simulator.run(compiled, shots=shots).result()
counts = result.get_counts()
estimated_phase = decode_phase(counts, evaluation_qubits)
estimated_eigenvalue = complex(
math.cos(2.0 * math.pi * estimated_phase),
math.sin(2.0 * math.pi * estimated_phase)
)
print("counts:", counts)
print("exact phase:", exact_phase)
print("estimated phase:", estimated_phase)
print("estimated eigenvalue:", estimated_eigenvalue)
print(compiled.draw(output="text"))
if __name__ == "__main__":
run_reference_qpe()The controlled powers are the essential mechanism of the algorithm. The \(k\)-th evaluation qubit controls \(U^{2^k}\). This exponential schedule encodes the binary expansion of the phase. If \(U\) is defined through a Hamiltonian simulation operator \(U=e^{-iHt}\), then the phase can be mapped to an energy eigenvalue. If \(H|\psi\rangle=E|\psi\rangle\), then:
\[ e^{-iEt}=e^{2\pi i\phi} \]
Therefore:
\[ E=-\frac{2\pi\phi}{t} \pmod{\frac{2\pi}{t}} \]
This expression shows that QPE requires careful spectral interpretation. The phase is periodic, so the selected simulation time \(t\) must be compatible with the expected energy range. If the range is not constrained, two different energies may produce indistinguishable phases.
If the input state is not an exact eigenstate but a linear combination \( |\eta\rangle=\sum_j \alpha_j|\psi_j\rangle \), QPE samples eigenphases according to the squared overlap:
\[ \Pr(\phi_j)\approx|\alpha_j|^2 \]
For this reason, a complete QPE implementation must record the state-preparation routine, controlled-unitary construction, evaluation register size, shot count, transpilation settings, and bit-ordering convention. Increasing \(m\) improves phase precision but increases circuit depth because higher controlled powers require more applications of the unitary. In noisy execution environments, the practical precision may be limited by decoherence and two-qubit gate errors rather than by register size alone.