Canonical Quantum Phase estimation¶
Download this notebook - canonical-qpe.ipynb
In this example notebook we will implement quantum phase estimation (\(QPE\)) in Guppy. We will define Guppy functions to implement controlled unitary and quantum Fourier transform (\(QFT\)) parts of the algorithm. We will also use this as a demonstration of Guppy modifiers which allow us to easily control and dagger quantum operations.
Based on this stack exchange post. See also the old pytket example notebook on QPE.
from guppylang import guppy
from guppylang.std.angles import pi
from guppylang.std.quantum import qubit, h, crz, x, measure_array, discard_array, cx, rz, collect_measurements
from guppylang.std.builtins import array, output, control, dagger, nat
Background¶
Phase estimation is an important quantum algorithm for estimating the eigenvalues of a unitary operator \(U\) to some precision. Quantum phase estimation appears as an important subroutine in Shor’s algorithm and various fault tolerant approaches to quantum chemistry. In this notebook we will consider the “canonical” phase estimation variant which is implemented by a pure unitary circuit.
If \(U\) is a unitary matrix, its eigenvalues must lie on the unit circle.
Here \(|\psi\rangle\) is an eigenstate of \(U\).
We estimate the eigenvalues by approximating \(\theta\) in the equation above.
We will consider a very simplified version of phase estimation wherein \(U\) is a diagonal matrix. This means the true eigenvalues can be read off the diagonal. This will allow us to clearly see that our implementation is correct.
Trivial State Preparation¶
Clearly the matrix has the eigenvalue \(e^{i \frac{\pi}{8}}\) corresponding to the eigenstate \(|11\rangle = (0, 0, 0, 1)^T\). We can prepare this trivial eigenstate with two Pauli \(X\) gates.
@guppy
def prepare_trivial_eigenstate() -> array[qubit, 2]:
q0, q1 = qubit(), qubit()
x(q0)
x(q1)
return array(q0, q1)
Controlled-\(U\) Operations¶
Next we need to create a subroutine applying a Controlled-\(U\) operation. This will be repeatedly applied and will kick back phase factors of \(e^{i \frac{\pi}{8}}\) onto the ancilla qubits. We will implement this diagonal \(U\) operator with two \(CX\) gates and an \(Rz\) rotation.
@guppy(controllable=True)
def u(q0: qubit, q1: qubit) -> None:
cx(q0, q1)
rz(q1, -pi / 4)
cx(q0, q1)
Note that as we will be controlling the \(U\) operation later, we mark the u function as controllable.
Quantum Fourier Transform¶
The final subroutine we need is the inverse quantum Fourier transform (\(QFT^\dagger\)). This has the effect of inducing destructive interference at the end of our circuit. This means that we are more likely to measure a single basis state (or a small set of basis states).
We can define a generalised Guppy function over \(n\) qubits. We do this by making our Guppy program for the \(QFT\) polymorphic over \(n\). We can implement the \(QFT\) with Hadamard and controlled-\(Rz\) gates followed by a layer of qubit swaps.
For the purposes of this example, we can implement a swap with three \(CX\) gates. Using this swap function as part of \(QFT^\dagger\), we mark the function as daggerable.
@guppy(daggerable=True)
def swap(q0: qubit, q1: qubit) -> None:
cx(q0, q1)
cx(q1, q0)
cx(q0, q1)
We can define a generalised Guppy function over \(n\) qubits. We do this by making our Guppy program for the QFT polymorphic over \(n\).
@guppy.comptime(daggerable=True)
def qft[n: nat](qs: array[qubit, n]) -> None:
for i in range(n):
h(qs[i])
for j in range(i + 1, n):
crz(qs[j], qs[i], pi / 2.0 ** (j - i))
# Reverse qubit order with swaps
for k in range(n // 2):
swap(qs[k], qs[n - k - 1])
Note that here we have defined a function implementing the QFT which we will then invert with the dagger modifier to implement the inverse QFT. The fact that qft is a comptime function allows us to apply the inverse even though the function contains loops. We also mark the qft function as daggerable as we did with swap.
For more on how comptime functions work, consult the comptime section of the Guppy language guide.
Note that we can have a QFT subroutine for any number of qubits that we like by adjusting the size of the input array.
The \(QPE\) Program¶
Now that we have defined subroutines for the state preparation, controlled unitaries and IQFT steps, we can combine these into a single function to perform quantum phase estimation.
First we define an array of measurement qubits of size \(m\). The more measurement qubits we have the more precise our estimate of the phase \(\theta\) will be.
Now we can define a function to implement phase estimation using our qft and u functions. Here we will have \(n\) measurement qubits. We fix the size of the initial state to have only two qubits.
The \(QPE\) construction can be generalised in a similar manner to the inverse \(QFT\) function. A larger value of \(n\) will mean that we can estimate the eigenphase of \(U\) to greater precision.
@guppy
def phase_estimation[n: nat](measured: array[qubit, n], state: array[qubit, 2]) -> None:
for i in range(n):
h(measured[i])
# Add 2^n - 1 controlled unitaries sequentially
for n_index in range(n):
control_index: int = n - n_index - 1
for _ in range(2**n_index):
with control(measured[control_index]):
u(state[0], state[1])
with dagger:
qft(measured)
Execution on Selene¶
Let’s execute this \(QPE\) program on the Selene emulator for 500 shots.
We can define a main function which includes our six qubit phase estimation subroutine and measurements. This can then be compiled for execution on the Selene simulator.
@guppy
def main() -> None:
state = prepare_trivial_eigenstate()
qubits_to_measure = array(qubit() for _ in range(4))
# Apply phase estimation subroutine.\n
phase_estimation(qubits_to_measure, state)
# Measure the qubits encoding the phase.
measurements = measure_array(qubits_to_measure)
# State prep qubits are not measured so have to be explicitly discarded.
discard_array(state)
# Create an output from the measured array.
output("c", collect_measurements(measurements))
n_shots = 500
sim_result = main.emulator(n_qubits=6).with_seed(5).with_shots(n_shots).run()
Now that we have executed our phase estimation instance on the Selene emulator we can analyse our results.
Let’s look at our measurement outcomes. In this highly idealised \(QPE\) instance we expect all of our measurement outcomes to be \(|0001\rangle\).
result_counter = sim_result.register_counts()["c"]
assert result_counter["0001"] == n_shots
print(result_counter)
Counter({'0001': 500})
We see that all of our measurements yield the basis state \(|0001\rangle\) which encodes the integer \(j=1\) in four bits.
Here \(n\) is the number of evaluation qubits (4 in our case). The value of \(j\) is given by the decimal value of the most frequent measurement outcome.