Quantum Phase Estimation Tutorial¶

Quantum phase estimation is one of the most important subroutines in quantum computing. In this document, we'll go through its implementation in PsiQDK Algorithms and how the code mirrors the design that was outlined in the Quantum Phase Estimation Domain Model ⧉. This document has a counterpart in the how-to guide on using the framework.

The domain model can be seen here (NOTE (SEAN): This model needs updating to be a bit more in line with the implementation's language):

We're going to focus on the quantum subroutines that define the model, how they are defined in Workbench algorthms and how they fit together.

QPEBuilder and presets¶

Most applications of QPE fit into a small number of common patterns. To facilitate the easy implementation of these patterns, the QPE framework in PsiQDK Algorithms exposes a qpe.QPEBuilder class that generates QPE programs by hooking up different components under the hood. In addition, a number of common presets are made available in qpe.presets that make it even easier to instantiate QPE programs. For a thorough overview of the available presets, see the how-to guide on using the QPE framework.

As an illustration of how the presets work, here's a quick setup of a QPE simulation whereby the phases are directly applied onto the time-frequency register. The result is post-processed using classical "cheat" methods to extract the probabilities of obtaining different results (this is useful for performing state vector simulations of QPE programs, though we won't be able to use it on a real QPU).

from psiqdk.algorithms.experimental import quantum_phase_estimation as qpe

qpu = QPU(num_qubits=4)
program = (
    qpe.presets.phase_debug(angle=180, bits_of_precision=4, window_function=qpe.RectWindowV2())
    .with_debug_resolver()
    .build()
)
result = program(qpu)
print(result.get_most_probable_phase())

QPE Program¶

The top level entry points for quantum phase estimation are QPEPrograms, classes that can be called to apply the full QPE subroutine, from Qubits allocation, to sampling the output of the time-frequency register and parsing the results.

For QFT-based QPE, the quantum program we'll use is the following:

class QFTBasedQPE:
    """Quantum program for running QFT based QPE."""
    def __call__(self, qpu: QPU):
        """Run the program."""
        for _ in self.frequency_resolver.iterate():
            qpu.reset(qpu.num_qubits)
            qpu.release_all_rotation_catalyst_qubits()
            signal_source_reg = self.signal_source_state_factory.get_signal_source_state(qpu)
            time_frequency_reg = Qubits(self.config.bits_of_precision, "time_frequency_reg", qpu)
            self.qpe_subroutine.compute(time_frequency_reg, signal_source_reg, self.config)
            self.frequency_resolver.resolve(time_frequency_reg)
        return self.frequency_resolver

The window, signal-generation, and inverse time–frequency transformation are delegated to QPESubroutine.compute (shown here as a single call for clarity).

This loop shows the steps in QPE:

• For each shot in our algorithm (controlled by the frequency resolver):
  • We reset the QPU and drop any rotation catalyst state — each shot starts from a clean slate
  • We use a SignalSourceStateFactory to allocate and prepare the signal source state, and we allocate the time–frequency register
  • We call QPESubroutine.compute, which applies the window, generates the signal onto the time–frequency register (via the time–signal sampling loop), and applies the inverse time–frequency transformation (e.g. inverse QFT)
  • We resolve the time–frequency register with the FrequencyResolver
  • (Register lifetimes follow your Workbench conventions; ancillary scratch used inside Qubricks is released by those routines)
• We then return the output of the shots

We'll go through each step in this algorithm in detail in this notebook.

Signal source state preparation factory¶

The first step in the "standard" QPE algorithm is the preparation of the state we're going to use to generate our signal – we refer to this as the "signal source state" (as opposed to "eigenstate", since it does not necessarily need to be an eigenstate of our signal unitary).

We could pass a Qubits object here, but there are a couple of issues with this:

1. This register needs to be re-initialized on every shot of the QPE algorithm.
2. We don't necessarily want to have to pass such a register – e.g. if we're probing the properties of window functions, we can do this by directly applying phases onto the time-frequency register, saving qubits that would otherwise hold the signal source state (and, importantly, halving the memory requirements for state vector simulation).

We can solve these two issues by using a factory pattern – rather than passing the object itself to the QPE program, we pass a class that can be used to construct the signal source state dynamically. Note that we could pass a Callable here, which would be simpler, but we use a class in anticipation of the fact that for flagship applications, the state preparation is likely to be quite complicated involving (at the very least!) Qubricks and temporary ancilla allocations/deallocations.

To see how this works, here's an example for a factory that allocates a Qubits register and pushes the signal source state into it, just to illustrate that we can do things that aren't just Qubrick calls (although we can do those too!):

class PushStateVectorFactory:
    def __init__(self, state: Iterable[float]):
        # process the state such that we have a power of two number of amplitudes
        state = np.pad(state, [0, (1 << (len(state) - 1).bit_length()) - len(state)])
        self.state = state

    def get_signal_source_state(self, qpu):
        num_qubits = len(self.state).bit_length() - 1
        signal_source_reg = Qubits(num_qubits, "signal_source", qpu)
        signal_source_reg.push_state(self.state)
        return signal_source_reg

We can also do more exotic things like

class FourierBasisStateFactory:
    def __init__(self, n:int, basis: int):
        self.n = n 
        self.basis = basis

    def get_signal_source_state(self, qpu: QPU):
        signal_source_reg = Qubits(self.n, "signal_source", qpu)
        signal_source_reg.write(self.basis)
        signal_source_reg.QFT()
        return signal_source_reg

We won't say anything more about the state factories since the problem of initial state preparation is so problem specific (and well known!). The main takeaway is that the pattern is designed to be flexible enough that any Workbench program can be included in the get_signal_source_state call and the QPE will use it as intended.

Quantum window function¶

The window functions are straightforward – they're just Qubricks and get computed like any other Qubrick. The only noteworthy thing to mention is that they now take the number of elements in the window to facilitate non-power-of-two QPE.

Signal generator¶

The signal generator is where the bulk of the complexity of QPE lies. The signal generator itself is straightforward, but its constituents are where the magic lies. Conceptually, all we're doing is looping over the times we're sampling (controlled by the time sampler) and sampling the signal at that time (using the signal sampler). The implementation reflects this concept -- both of these samplers are implemented as iterators and we simply loop over them:

class SignalGenerator(Qubrick, Generic[T]):
    """Qubrick for generating the signal for QPE."""

    def __init__(self, sampling_strategy: SamplingStrategy, **kwargs):
        self.sampling_strategy = sampling_strategy
        super().__init__(**kwargs)

    def _compute(
        self,
        time_frequency_register: Qubits,
        signal_source_register: T,
        num_time_steps: int | None = None,
        ctrl: Qubits | None = None,
    ):
        for time_sample_register, signal_func in self.sampling_strategy.sample_from(
            time_frequency_register, num_time_steps
        ):
            signal_func(signal_source_register, time_sample_register)

Note that unlike most components we're talking about here, users are typically expected to use this class as-is, rather than implementing a new version themselves.

Eagle-eyed readers might notice that the SignalGenerator actually accepts a sampling strategy rather than time and signal samplers. This is an implementation detail that just makes the signal generator class a bit nicer. The sampling strategy class looks like this:

@dataclass(frozen=True)
class SamplingStrategy(Generic[T]):
    """Utility class for combining time and signal sampling iterators.

    Users should be able to take this class and use it without modification with any compatible time_sampler and
    signal_sampler.
    """

    time_sampler: TimeSampler
    signal_sampler: SignalSampler

    def sample_from(
        self,
        time_frequency_register: Qubits,
        num_time_steps: int | None = None,
    ) -> Generator[tuple[Qubits, Callable[[T, Qubits], None]], None, None]:
        yield from zip(
            self.time_sampler.sample_from(time_frequency_register, num_time_steps),
            self.signal_sampler.sample_signal(num_time_steps),
        )

Again, this should work with any compatible time and signal sampler, users shouldn't need to write their own.

Time sampler¶

The signal generation step in QPE consists of generating samples of the signal over times $t$ and storing them in the phases of the time-frequency register at the corresponding basis state $|t\rangle$, i.e.

$$ |t\rangle \mapsto e^{2 \pi i \phi t}|t\rangle $$

There are many strategies for doing this, but they all involve sampling some set of time basis states and then generating the signal for those times.

The most common approach (the one used in textbook QPE) is to perform a binary expansion of the signal by sampling times with a set bit in each position of their binary representation – in more familiar terms, to control over each qubit of the time-frequency register individually!

class BitwiseTimeSampler(TimeSampler):
    """Yields samples of the time-frequency register where each bit in the register is set.

    That is, on the nth iteration, yields the nth bit such that all time steps t for which
    `t & (1 << n) == 1` are sampled in superposition.

    Note:
        This is the "standard" sampling protocol that you see in, e.g. Nielsen and Chuang.
    """

    def _sample_from(self, time_frequency_register, num_time_steps: int | None = None) -> Generator[Qubits, Any, Any]:
        # TODO: in a final implementation, we would also handle the adders for
        # non-power-of-two QPE here, but for now we just support power-of-two.
        for qubit in time_frequency_register:
            yield qubit

We can also sample time in a unary manner, controlling on each time individually:

class UnaryTimeSampler(TimeSampler):
    def _sample_from(self, time_frequency_register, num_time_steps: int | None = None) -> Generator[Qubits, Any, Any]:
        for t in range(num_time_steps):
            yield time_frequency_register == t

This version is simple, but we can also implement the unary iteration circuit from Encoding Electronic Spectra in Quantum Circuits with Linear T Complexity ⧉ for more efficient compilation of this circuit.

Signal sampler¶

This will be the main entry point for integrating new circuits into the QPE framework. Signal samplers inherit from the SignalSampler base class:

class SignalSampler(ABC, Generic[SigReg]):
    """Base class for signal samplers in QPE."""

    MAX_SAMPLE_SIZE = 10_000_000

    @abstractmethod
    def _sample_signal(
        self, signal_source_register: SigReg, time_sample_register: Qubits, sample_step: int, num_time_steps: int | None = None
    ):
        raise NotImplementedError

    def sample_signal(self, num_time_steps: int | None) -> Generator[Callable[[SigReg, Qubits], None], None, None]:
        sample_step = 0
        # generator bound will typically be set by the TimeSampler
        # this check is just to avoid infinite loops
        while sample_step < self.MAX_SAMPLE_SIZE:
            yield partial(self._sample_signal, sample_step=sample_step, num_time_steps=num_time_steps)
            sample_step += 1

Just like Qubricks, there is a method to call when using the class, sample_signal and a private method to use when writing the class, _sample_signal. The _sample_signal method is what users should overwrite, and can contain any python code (just like _compute in Qubricks).

This acts as an adaptor for QPE, taking the signal unitary (which in most cases comprises some Qubrick or composition of multiple Qubricks) and wrapping it such that it is compatible with the rest of the QPE machinery.

In some cases, this may involve a needless wrapping of a Qubrick in an additional class that needs to be added, which would be annoying to use. For flagship routines like Trotterization and Qubitization, this is less of an issue, since we can define top-level signal samplers for those and just fill in the details for specific applications as you would with a Qubrick, but for generic unitaries this becomes a problem.

To resolve this, we also added a generic signal sampler that can be used to incorporate any Qubrick, QPU op or function into QPE without needing to define an additional class. This is primarily for prototyping and debugging rather than flagship application development – we can't take advantage of bidirectional phase kickback using this method, for example.

The sampler accepts an op of exactly one shape: something that takes the signal source register and the time sample register and nothing else. Anything with a wider signature is adapted first by bind_generic_unitary_arguments, which freezes the remaining arguments and works out which parameter each of the two registers belongs in. The signal source register goes to the first parameter that was not frozen; the time sample register goes to ctrl, cond or condition_mask if the op has one of those, and to the next free parameter otherwise. A Qubrick has its parameter names read from _compute (or from _estimate, if that is all it implements), but is still invoked through compute, so temporary register allocation and control reduction happen as usual.

def bind_generic_unitary_arguments(
    op: Qubrick | Callable,
    *,
    skip_params: Iterable[str] = (),
    **kwargs,
) -> Callable[[BaseQubits, BaseQubits], None]:
    if isinstance(op, Qubrick):
        signature, target = inspect.signature(_qubrick_compute_method(op)), op.compute
    else:
        signature, target = inspect.signature(op), op

    signal_source_param, time_sample_param = _infer_and_validate_param_names(signature, set(skip_params), kwargs)

    def _partial_func(signal_source_register: BaseQubits, time_sample_register: BaseQubits, /):
        target(
            **kwargs,
            **{signal_source_param: signal_source_register, time_sample_param: time_sample_register},
        )

    return _partial_func

Which leaves the sampler itself with very little to do:

class GenericUnitarySignalSampler(SignalSampler[Qubits]):
    """This is an adaptor for handling generic, arbitrary unitaries for QPE.

    The op is either a callable taking the signal source register and the time sample register and nothing else, or a
    Qubrick that already has that shape: its `_compute` needs a `ctrl` parameter for the time sample register, and one
    other parameter for the signal source register.

    Any other unitary has to have its remaining arguments frozen first, with
    `bind_generic_unitary_arguments`. The `generic_unitary` preset does this for you.

    While this is pretty generic and flexible, there are unitaries that can't
    easily be fit into this structure. The suggested API is to write a custom
    adaptor for whatever specific unitary you have in mind.

    Note:
        - For high performance applications, this is probably not the way to go.
        - We are not doing bidirectional phase kickback, as there is no generic
          way to do so.
        - This is only intended to support textbook (e.g. Nielsen and Chuang) QPE,
          with the exponentiation achieved by applying the unitary exponentially
          many times.
    """

    def __init__(self, op: Qubrick | Callable[[BaseQubits, BaseQubits], None]) -> None:
        if isinstance(op, Qubrick):
            _validate_qubrick_has_ctrl_param(op)
            self.fn = bind_generic_unitary_arguments(op)
        elif callable(op):
            self.fn = op
        else:
            raise TypeError("GenericUnitarySignalSampler op must be a Qubrick or callable.")
        super().__init__()

    def _sample_signal(
        self,
        signal_source_register: Qubits,
        time_sample_register: Qubits,
        sample_step: int,
        qpe_config: QPEConfig,
    ):
        for _ in range(qpe_config.get_num_iterations_per_time_step(sample_step)):
            self.fn(signal_source_register, time_sample_register)

Users going through the generic_unitary preset do not have to call the binder themselves. The preset takes the op directly and takes its other arguments as a dictionary, e.g. bound_arguments={"theta": 45.0}.

Time-frequency transformation¶

Much like the window function, the time-frequency transformation is simply a Qubrick, no additional functionality is added.

Frequency resolver¶

When we talk about QPE, we typically leave the conversation at the previous step – we apply our QPE circuit, and then perhaps we have some statement about needing to repeat it a certain number of times to actually get our answer with a satisfactory degree of certainty. But this is a key aspect of running the algorithm! Skipping over this detail is problematic, as it means that part of our algorithm will be untested.

In the new framework, this iterative nature is built into the implementation. But what if we're just benchmarking, say, a window function and we don't want to have to iterate over a bunch of trials because we want to inspect the state vector to ensure we have produced the correct probability distribution? We have a trick up our sleeves here – we can use the same API to do both of these things! The way we do this is with a FrequencyResolver class.

What is a FrequencyResolver?¶

A FrequencyResolver is a class that implements two methods:

1. iterate, which returns an Iterable
2. resolve, which takes in the time_frequency_register and returns nothing

These two methods give us a dynamic way to control how the iteration in QPE proceeds. At a high level, what we're going to do is:

1. Iterate over the elements yielded by FrequencyResolver.iterate and for each of them
2. Run the QPE subroutine, the result of which is the frequencies in the signal loaded into the time_frequency register
3. Which can then be interpreted by the FrequencyResolver.resolve method.

Why not just a simple loop?¶

The above might seem a little over-engineered – why not just pass a number of iterations to the QPE program and loop over that, storing up measurement results? There are a few reasons:

1. For full applications, it's likely the case that the number of iterations we expect is an upper bound on what we actually get in practice. We want some functionality to analyze the statistical properties of the outcomes and break early if we are satisfied that we've got the correct result.
2. We may want to do some more complicated post-processing of the outcomes, for example throwing away measurements from some of the qubits in order to boost the success probability of getting the desired outcome.
3. When debugging or unit testing some parts of the subroutine (e.g. the window functions) we may not want to measure the time-frequency register at all, but rather peek the read probabilities of the different amplitudes within it.

By using a class to handle this, we abstract away these details such that the QPE program itself doesn't need to know about any of them.

Parsing frequencies¶

One other role that we can delegate to a dedicated class is the interpretation of measurement outcomes in QPE. This is something that is notoriously error-prone (Do I need to include a factor of $2\pi$? Are my energies normalized? Did I remember the negative sign in the conversion from qubitization?) and which needs to be handled every time we use QPE. Since we have a class dedicated to handling frequencies, we can do this parsing there, allowing us to (for example) print out the results as "raw" fixed point numbers, as binary strings and as post-processed values all using the same class.

An example of how this can be extremely useful can be seen in the following class that allows us to extract out the most probable result from the QPE to get a sense of the behavior in the limit of taking an infinite number of shots. This is extremely useful for debugging/testing QPE and subroutines that constitute it.

class DebugQFTQPEFrequencyResolver:
    """Class to handle QPE data accumulation.

    This class specifically handles the cases where we use
    classical "cheat" methods to pull out information about
    probabilities etc rather than using measurements as you'd
    have to do for a real QPU.
    """

    frequency_interpreter: FrequencyInterpreter = default_frequency_interpreter

    def iterate(self):
        """Iteration for building up QFT QPE statistics.

        Since we're debugging, this is effectively just a dummy
        method that just yields None once. The statistics we're
        going to get are going to come from read probabilities
        instead of repeated measurements.
        """
        yield

    @property
    def measurement_outcomes(self):
        return {
            self.frequency_interpreter(phase): probability
            for phase, probability in zip(self.phases, self.probabilities)
        }

    @property
    def fixed_point_measurement_outcomes(self):
        """Get the measurement outcomes interpreted as raw fixed point numbers."""
        return {phase: probability for phase, probability in zip(self.phases, self.probabilities)}

    def resolve(self, time_frequency_register: Qubits):
        num_frequencies = 1 << time_frequency_register.num_qubits
        phases = np.arange(num_frequencies) / num_frequencies
        probabilities = [time_frequency_register.peek_read_probability(i) for i in range(1 << time_frequency_register.num_qubits)]
        self.phases = phases
        self.probabilities = probabilities

    def calculate_holevo_variance(self):
        # the calculation of the Holevo variance assumes we're just interpreting the values as phases - we thus skip
        # the frequency_interpreter call and do the multiplication here
        return calculate_holevo_variance(2 * np.pi * self.phases, self.probabilities)

    def get_most_probable_phase(self):
        return list(self.measurement_outcomes.keys())[np.argmax(list(self.measurement_outcomes.values()))]

Different frequency resolvers can be accessed via the QPEBuilder as methods (e.g. QPEBuilder.with_debug_frequency_resolver()).

Note that the task of interpreting the phase is delegated to a callable – for example, here's a factory to generate the interpreter callable for qubitization:

def get_qubitization_frequency_interpreter(norm: float):
        def qubitization_frequency_resolver(phase: float):
            return -1 * np.cos(2 * np.pi * phase) * norm
        return qubitization_frequency_resolver