Canonical Quantum Phase estimation¶
Download this notebook - canonical-qpe.ipynb
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
from guppylang.std.builtins import result, array, mem_swap
from pytket.circuit import DiagonalBox, QControlBox
from pytket.passes import AutoRebase
from pytket import OpType
import numpy as np
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” QPE 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 \(U\) 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 can generate this using a pytket DiagonalBox where we apply an additional control.
d_box = DiagonalBox(np.array([1, 1, np.exp(1j * np.pi / 4), np.exp(1j * np.pi / 8)]))
controlled_u_op = QControlBox(d_box, 1)
Circuits created using pytket can be loaded into Guppy as functions using the load_pytket function (or alternatively by specifying a function stub with the same signature as the circuit and annotating it with @pytket(circ)). However, the circuit being loaded can only contain gates that exist in the Guppy quantum library, so we need to rebase the generated circuit to remove a TK1 gate.
rebase = AutoRebase({OpType.CX, OpType.Rz, OpType.H, OpType.CCX})
circ = controlled_u_op.get_circuit()
rebase.apply(circ)
controlled_u = guppy.load_pytket("controlled_u_circuit", circ, use_arrays=False)
Inverse Quantum Fourier Transform¶
The final subroutine we need is the inverse quantum fourier transform (IQFT). 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 defining a natural number variable \(n\). Our Guppy program for the IQFT is then polymorphic over \(n\).
n = guppy.nat_var("n")
@guppy
def inverse_qft(qs: array[qubit, n]) -> None:
# Reverse qubit order with swaps
for k in range(n // 2):
mem_swap(qs[k], qs[n - k - 1])
for i in range(n):
h(qs[n - i - 1])
for j in range(n - i - 1):
crz(qs[n - i - 1], qs[n - i - j - 2], -pi / 2 ** (j + 1))
Note that we can have a IQFT 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 inverse_qft and controlled_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(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):
controlled_u(measured[control_index], state[0], state[1])
inverse_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()
measured = array(qubit() for _ in range(4))
phase_estimation(measured, state)
# state qubits are not measured so have to be explicitly discarded
discard_array(state)
# Create a result from the measured array
result("c", measure_array(measured))
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.