Phase estimation into an ancilla register

What it fills

• Estimate the eigenphase of a unitary

  Given a circuit whose controlled powers you can apply, and a routine preparing a state with non-negligible overlap on one of its eigenvectors, return that eigenvector's phase as a number with an error bar. The phase is read out of an ancilla, never out of the system register — the system is only ever the thing the controlled powers act on.

When it applies

Kitaev's paper is where measuring an eigenvalue of a unitary is introduced — "Our method is based on a procedure for measuring an eigenvalue of a unitary operator" — and it is introduced in the service of something else, the Abelian Stabilizer Problem, which is the shape ADR-0026 permits. The construction is the ancilla-and-controlled-U one: with $S$ the Hadamard, "the operator $\Xi(U) = S \Lambda(U) S$ is a measurement operator for the observable $\varphi$", and precision comes from having the controlled powers available — "the situation is different if we have in our disposal the operators $\Lambda(U^k)$ for all $k$".

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

given  U a unitary on a subspace N of a register A, eigenvalues of the form
       lambda(phi) = exp(2*pi*i*phi), phi a real number mod 1 ("the
       observable phi"); E(U,phi) the corresponding eigenspaces
                                                   (Sec. "Important example")
       a state |xi> in N -- need NOT be an exact eigenvector: Xi(U) below is
           a measurement operator in the paper's sense, and for a general
           |xi> the composite probability formula weights each eigenspace by
           P(|xi>, V), so the readout acts correctly on a superposition
                                                       (Def. 5; Lemma 8, Eq. 18)
       l, the number of precision bits wanted;  r = 2^l - 1
       eps > 0, the total error probability that may be tolerated
       the basis Q = U(B1) union {tau, controlled-tau}          (Sec. 2.4)
# NB on case: Lambda = the controlled-U operator (the paper's capital Lambda);
#   lambda = the eigenvalue exp(2*pi*i*phi) (the paper's lower-case lambda).

requires  ONE extra gate beyond Q -- the register gate
    U^[0,r] : |a,xi> -> |a> (x) U^a|xi>,  a in {0,...,r} = B^l          (Eq. 20)
    Each controlled power Lambda(U^(2^j)), j = 0,...,l-1, is NOT assumed
    separately: "By Lemma 7, the operators Lambda(U^(2^j)) ... can be
    represented in terms of U^[0,r]." So U^[0,r] is the one extra primitive
    this construction assumes                             (Lemma 7; remark after Eq. 20)

# --- the single-level readout operator -------------------------------------
S = (1/sqrt(2)) * [[1, 1], [1, -1]]                (Sec. "Important example")
Xi(U)[A,1] = S[1] . Lambda(U)[1,A] . S[1]                              (Eq. 16)
    # a measurement operator for the observable phi. On |xi,0> with |xi> in
    #   E(U,phi) it leaves the extra bit in
    #   |eta> = ((1+lambda(phi))/2)|0> + ((1-lambda(phi))/2)|1>, so
    #   P(phi,0) = (1+cos(2*pi*phi))/2,  P(phi,1) = (1-cos(2*pi*phi))/2
    #                                                        (Eq. 16, Eq. 17)

# --- Lemma 9: turn one such bit into a fixed-precision estimate -------------
apply Xi(U) s times to the same A with fresh extra bits 1,...,s; measure them
    and let y = number of 1s. Then for constant delta > 0,
    Prob[ |y/s - P(phi,1)| > delta ] <= 2*exp(-c(delta)*s), c(delta) > 0
    # substituting iU for U turns cos(2*pi*phi) into -sin(2*pi*phi), so both
    #   cos and sin are measurable, which determines phi
    # => precision delta at error probability <= eps in O(log(1/eps))
    #   operations over Q union {Lambda(U)}. NOTE delta is a CONSTANT here:
    #   Lemma 9 buys cheap confidence, not arbitrary precision      (Lemma 9)

# --- per-level localization: pin 2^j*phi one register level at a time -------
for j = 0 to l-1:
    using Lambda(U^(2^j)) obtained from U^[0,r] by Lemma 7, localize
        2^j*phi (mod 1) into one of the 8 overlapping width-1/4 intervals
        [(s-1)/8, (s+1)/8],  s = 0,...,7,  spaced 1/8 apart
    # cost O(log(l/eps)) operations at error probability <= eps/l on this
    #   level. This is the paper's own stated figure. It is stated flatly:
    #   the paper does NOT spell out the mechanism, and the natural reading
    #   -- that it is Lemma 9's repeat-and-count trick run per level with
    #   eps/l in place of eps -- is an inference, not a cited step
    #                                                (remark after Eq. 20)

# --- combine the l level-estimates -----------------------------------------
phi~ = COMBINE(interval estimates for j = 0,...,l-1)
    # "by a polynomial algorithm" -- classical, and left unnamed and
    #   unnumbered in the paper; COMBINE is this listing's placeholder, not
    #   a name Kitaev gives                              (remark after Eq. 20)

return  phi~, with |phi~ - phi| <= (1/8)*2^(-(l-1)) = 2^(-l-2), total error
    probability <= eps, in O(l*log(l/eps)) + poly(l) operations over
    Q union {U^[0,r]}; U^[0,r] itself is used at most O(l*log(l/eps))
    times                                                          (Lemma 10)
    # Lemma 10 states eps as the achieved TOTAL directly; the paper does not
    #   write out a union bound over the l per-level eps/l failures

# l is fixed once, by the precision wanted (2^-l-2): the register does not
#   grow as eps shrinks. Only the operation count grows with 1/eps, and only
#   logarithmically -- that is what makes "having Lambda(U^k) at our
#   disposal for all k" strictly better than iterating Lemma 9 alone, whose
#   cost blows up as the precision is tightened.
#   Caution on the source's wording: the paper writes that the cost "grows
#   polynomially in delta", which as written is backwards -- the cost grows
#   as delta SHRINKS. The usable form is: polynomial in 1/delta at fixed
#   error probability                                    (remark before Eq. 20)

# NOT part of this construction: an inverse QFT on the register followed by a
#   computational-basis measurement. That is the standard textbook framing of
#   phase estimation; it does not appear in Kitaev 1995, which reads the
#   phase out by the interval localization above.
# Also NOT part of it: Coppersmith's approximate Fourier transform AFFT_m.
#   It prices a DIFFERENT transform -- the ordinary FFT inside Shor's
#   factoring circuit, decomposed into L one-qubit gates P_J and L(L-1)/2
#   two-qubit gates Q_JK on an L-qubit register (Coppersmith's L, not this
#   record's l). AFFT_m deletes every Q_JK with K >= J+m, leaving about L*m
#   two-qubit gates, which multiplies each output amplitude by a phase
#   exp(i*eps') with |eps'| <= 2*pi*L*2^(-m) -- about 3/1000 in his own
#   worked case L=500, m=20, an instance, not a general guarantee
#                            (Coppersmith, "Approximate Fourier Transform";
#                             "Calculating the AFFT")

Cost, as the source states it

Coppersmith prices only the Fourier read-out of Shor's factoring construction — not a generic controlled-$U$ ladder, and not Kitaev's interval-narrowing scheme priced above. On a register of $L$ qubits his exact transform uses $L$ single-qubit gates $P_J$ and $L(L-1)/2$ two-qubit gates $Q_{JK}$ acting on qubits $J$ and $K$. Deleting every $Q_{JK}$ with $K \ge J+m$, for an integer truncation parameter $m$, leaves about $Lm$ two-qubit operations. That truncation multiplies each output amplitude by a phase $\exp(i\epsilon)$ with $|\epsilon| \le 2\pi L 2^{-m}$; his own worked case $L=500$, $m=20$ evaluates to about $3/1000$, an instance rather than a guarantee. In his factoring variables ($n$ the modulus bit-length, $q = 2^L \approx 5n^2$) he estimates $(\log n)(\log\log n + \log(1/\epsilon))$ operations for a required final precision $\epsilon$, against $(\log n)^2$ for the exact transform and $(\log n)^3$ for Shor's mixed-radix one. All three are the paper's own informal estimates, not proved bounds.

Implementations

• Qiskit's register QPE: circuit-library primitive plus algorithm runner

  Qiskit implements the coherent-register-plus-inverse-QFT circuit as a first-class primitive in its circuit library, and the separate qiskit-algorithms package runs that circuit end to end against a sampler primitive and decodes the measured bitstring into a phase. Both name the same construction: Hadamards on an evaluation register, controlled powers of the unitary into the target register, then an inverse quantum Fourier transform on the register.

  `qiskit.circuit.library.phase_estimation(num_evaluation_qubits, unitary, name="QPE")` (qiskit/circuit/library/phase_estimation.py) builds an evaluation register of `num_evaluation_qubits`, puts it in Hadamard superposition, composes `unitary.power(2**j).control(annotated=False)` into the target register for each evaluation qubit j, then appends the inverse of `QFTGate(num_evaluation_qubits)` followed by a bit-reversal `PermutationGate`. `unitary` may be a `QuantumCircuit` or a `Gate`. The older class-based constructor, `qiskit.circuit.library.PhaseEstimation.__init__`, built the same circuit against `QFT(num_evaluation_qubits, inverse=True, do_swaps=False).reverse_bits()`; it is now decorated `@deprecate_func(since="2.1", additional_msg='Use qiskit.circuit.library.phase_estimation instead.', removal_timeline="in Qiskit 3.0")` on Qiskit's main branch. The class-based `qiskit_algorithms.PhaseEstimation` (qiskit_algorithms/phase_estimators/phase_estimation.py in the separate qiskit-algorithms package) is described in its own docstring as running "QPE with a multi-qubit register for reading the phases"; its `construct_circuit` composes a caller-supplied state-preparation circuit onto `circuit.library.PhaseEstimation(num_evaluation_qubits, unitary)`, and the whole thing runs on a `BaseSamplerV2` sampler supplied at construction (the constructor raises if none is given). The returned `PhaseEstimationResult.phase` property returns "the most likely phase": the bit string with the highest measured weight, decoded to a float in [0,1). `PhaseEstimation` sits in the same `qiskit_algorithms.phase_estimators` subpackage as `IterativePhaseEstimation`, the single-ancilla route, and both are re-exported at the package's top level.

  No dataset — the unitary and its eigenstate are supplied by the caller as circuits or gates. The circuit-library function's own docstring example runs it on a `SparsePauliOp`-derived `PauliEvolutionGate` Hamiltonian evolution; the qiskit-algorithms test suite exercises single-qubit Pauli eigenproblems (Z, X, and an RZ(pi/2) rotation).

  https://github.com/Qiskit/qiskit, qiskit/circuit/library/phase_estimation.py — function `phase_estimation` (current) and deprecated class `PhaseEstimation`. https://github.com/qiskit-community/qiskit-algorithms, qiskit_algorithms/phase_estimators/phase_estimation.py — class `PhaseEstimation`, exported as `qiskit_algorithms.PhaseEstimation`.

  The qiskit-algorithms test suite (test/test_phase_estimator.py, class TestPhaseEstimation) runs `PhaseEstimation(num_evaluation_qubits=6, sampler=StatevectorSampler(seed=42))` — Qiskit's statevector-simulator sampler primitive — against a Z-gate unitary prepared in the |1> eigenstate by an X gate, and asserts the returned `.phase` equals 0.5 exactly (test_qpe_Z_sampler); the same six-qubit register recovers phase 0.125 for an RZ(pi/2) unitary on the X-prepared state and 0.875 on the untouched |0> state (test_qpe_RZ_sampler), all on the same noiseless statevector sampler.

• Qualtran's TextbookQPE bloq

  Qualtran, Google Quantum AI's Python library for building and resource-counting quantum algorithms as composable bloqs, records the same Hadamard-register/controlled-power/inverse-QFT circuit as a single bloq, attributed in its own docstring to "Chapter 5.2 of Neilson & Chuang" [sic in the source] rather than to Kitaev.

  `TextbookQPE` (qualtran/bloqs/phase_estimation/text_book_qpe.py), an `attrs.frozen` `GateWithRegisters`, takes a `unitary: Bloq`, a configurable `ctrl_state_prep: QPEWindowStateBase` for the register's initial state (uniform Hadamard superposition by default, via `RectangularWindowState`), and a configurable `qft_inv`, defaulting to `QFTTextBook(self.m_bits, with_reverse=True).adjoint()`. Its `decompose_from_registers` applies `ctrl_state_prep` to the phase register `qpe_reg`, then for each phase qubit applies `cirq.pow(unitary_op.controlled_by(qbit), 2**i)`, then applies `qft_inv` to the register — the same construction as Qiskit's, expressed as bloqs rather than gates. The class documents one closed-form choice of register size: for `n` bits of precision at failure probability $\delta$, its Eq. 5.35 of Nielsen & Chuang gives $m = n + \lceil \log_2(2 + 2/\delta) \rceil$ evaluation qubits, reachable by constructing `TextbookQPE(unitary, RectangularWindowState.from_precision_and_delta(precision, delta))`.

  No dataset. The module's own tests run it on a single-qubit `ZPowGate` with a literal rotation, on a qubitization walk operator built from small (2-4 term) synthetic Pauli Hamiltonians, and separately on a walk operator for the Hubbard model wrapped by `HamiltonianSimulationByGQSP`.

  https://github.com/quantumlib/Qualtran, Python, Apache License 2.0. qualtran/bloqs/phase_estimation/text_book_qpe.py, class `TextbookQPE`.

  The bloq's own test file (qualtran/bloqs/phase_estimation/text_book_qpe_test.py) decomposes `TextbookQPE(ZPowGate(exponent=2*theta), RectangularWindowState(3))` to a Cirq circuit and simulates it on `cirq.Simulator()`, a state-vector simulator, recovering theta in {0.234, 0.78, 0.54} to within an asserted error bound of 0.1 at 3 bits of register precision (test_textbook_phase_estimation_zpow_theta). A separate test, at 5 bits of precision on qubitization walk operators built from small (2-4 term) synthetic Pauli Hamiltonians, simulates the same way and checks the recovered eigenvalue against each Hamiltonian's own `np.linalg.eigh` decomposition to within `atol=0.05`, on circuits kept under 23 qubits total (test_textbook_phase_estimation_qubitized_walk). A third, separate test runs `TextbookQPE` at 3 bits of precision against a walk operator for the Hubbard model wrapped by `HamiltonianSimulationByGQSP`, but checks only that the bloq decomposes validly (`qlt_testing.assert_valid_bloq_decomposition`) — not a numeric eigenvalue recovery (test_qpe_of_gqsp).

• Cirq's phase-estimator worked example

  Cirq, Google Quantum AI's circuit-construction and simulation library, does not ship a reusable phase-estimation class the way Qiskit does; the register construction instead lives as a runnable worked example in the repository's own examples directory, built directly from Cirq's own `cirq.qft` gate.

  examples/phase_estimator.py's `run_estimate(unknown_gate, qnum, repetitions)` allocates `qnum` `LineQubit`s plus one ancilla `LineQubit`, applies `cirq.H.on_each(*qubits)`, then for each of the `qnum` register qubits applies `unknown_gate.on(ancilla).controlled_by(qubits[i]) ** (2**i)`, then applies `cirq.qft(*qubits, without_reverse=True) ** -1` as the inverse transform, then measures all `qnum` register qubits into a classical key — the same Hadamard/controlled-power/inverse-QFT structure as the other two artefacts, run through `cirq.sample`.

  No dataset. The example's own `example_gate(phi)` is a synthetic single-qubit `MatrixGate` with eigenvector |0> and eigenvalue $\exp(2\pi i \phi)$, swept over $\phi \in \{0.0, 0.1, \ldots, 0.9\}$.

  https://github.com/quantumlib/Cirq, examples/phase_estimator.py, function `run_estimate` (circuit construction) and `experiment`/`main` (the parameter sweep and reporting).

  Run via `cirq.sample`, Cirq's own circuit sampler, at register sizes `qnum` in (2, 4, 8) with 100 repetitions per phase value. The file's own recorded example output for `qnum=8` reports an RMS error of 0.0011 across the ten swept phases, e.g. target=0.5000 recovered as estimate=0.5000=128/256 and target=0.7000 recovered as estimate=0.6992=179/256.

What it needs

Nobody has taken this apart yet. That is a gap in this graph, not a claim that the method has no parts.

Other ways to fill the same slot

Different approaches

• Iterative phase estimation on one ancilla

  Use one ancilla and measure it, over and over, least significant bit first. Each measured bit is fed back classically as a rotation angle on the next round, so the register the other route holds in superposition is replaced by a classical string that grows one bit at a time.

In the Atlas

• Quantum phase estimation

  The eigenphase-estimation primitive behind chemistry, simulation, and period-finding workflows.

• Quantum Fourier Transform

  The unitary that maps computational basis states to Fourier-basis phase patterns, and the phase-estimation primitive behind Shor's and HHL's speedups.

Sources

• Quantum measurements and the Abelian Stabilizer Problem↗A. Yu. Kitaev · 1995theory yes · simulation ? · hardware no
• An approximate Fourier transform useful in quantum factoring↗D. Coppersmith · 2002theory yes · simulation ? · hardware no