{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Quantum Phase Estimation\n","\n","**Download this notebook - {nb-download}`phase_estimation.ipynb`**\n","\n","When constructing circuits for quantum algorithms it is useful to think of higher level operations than just individual quantum gates. In `pytket` we can construct circuits using box structures which abstract away the complexity of the underlying circuit. This notebook is intended to complement the [boxes section](https://tket.quantinuum.com/user-guide/manual/manual_circuit.html#boxes) of the user manual which introduces the different box types.\n","\n","To demonstrate boxes in `pytket` we will consider the Quantum Phase Estimation algorithm (QPE). This is an important subroutine in several quantum algorithms including Shor's algorithm and fault-tolerant approaches to quantum chemistry.\n","\n","## Overview of Phase Estimation\n","\n","The Quantum Phase Estimation algorithm can be used to estimate the eigenvalues of some unitary operator $U$ to some desired precision.\n","\n","The eigenvalues of $U$ lie on the unit circle, giving us the following eigenvalue equation\n","\n","$$\n","\\begin{equation}\n","U |\\psi \\rangle = e^{2 \\pi i \\theta} |\\psi\\rangle\\,, \\quad 0 \\leq \\theta \\leq 1\n","\\end{equation}\n","$$\n","\n","Here $|\\psi \\rangle$ is an eigenstate of the operator $U$. In phase estimation we estimate the eigenvalue $e^{2 \\pi i \\theta}$ by approximating $\\theta$.\n","\n","\n","The circuit for Quantum phase estimation is itself composed of several subroutines which we can realise as boxes.\n","\n","![](images/phase_est.png \"Quantum Phase Estimation Circuit\")"]},{"cell_type":"markdown","metadata":{},"source":["QPE is generally split up into three stages\n","\n","1. Firstly we prepare an initial state in one register. In parallel we prepare a uniform superposition state using Hadamard gates on some ancilla (measurement) qubits. The number of ancilla qubits determines how precisely we can estimate the phase $\\theta$.\n","\n","2. Secondly we apply successive controlled $U$ gates. This has the effect of \"kicking back\" phases onto the ancilla qubits according to the eigenvalue equation above.\n","\n","3. Finally we apply the inverse Quantum Fourier Transform (QFT). This essentially plays the role of destructive interference, suppressing amplitudes from \"undesirable states\" and hopefully allowing us to measure a single outcome (or a small number of outcomes) with high probability.\n","\n","\n","There is some subtlety around the first point. The initial state used can be an exact eigenstate of $U$ however this may be difficult to prepare if we don't know the eigenvalues of  $U$ in advance. Alternatively we could use an initial state that is a linear combination of eigenstates, as the phase estimation will project into the eigenspace of $U$."]},{"cell_type":"markdown","metadata":{},"source":["We also assume that we can implement $U$ with a quantum circuit. In chemistry applications $U$ could be of the form $U=e^{-iHt}$ where $H$ is the Hamiltonian of some system of interest. In the textbook algorithm, the number of controlled unitaries we apply scales exponentially with the number of measurement qubits. This allows more precision at the expense of a larger quantum circuit."]},{"cell_type":"markdown","metadata":{},"source":["## The Quantum Fourier Transform"]},{"cell_type":"markdown","metadata":{},"source":["Before considering the other parts of the QPE algorithm, lets focus on the Quantum Fourier Transform (QFT) subroutine.\n","\n","Mathematically, the QFT has the following action.\n","\n","$$\n","\\begin{equation}\n","QFT : |j\\rangle\\ \\longmapsto \\frac{1}{\\sqrt{N}} \\sum_{k=0}^{N - 1} e^{2 \\pi ijk/N}|k\\rangle, \\quad N= 2^n\n","\\end{equation}\n","$$\n","\n","This is essentially the Discrete Fourier transform except the input is a quantum state $|j\\rangle$.\n","\n","We can build the circuit for the $n$ qubit QFT using $n$ Hadamard gates $\\lfloor{\\frac{n}{2}}\\rfloor$ swap gates and $\\frac{n(n-1)}{2}$ controlled unitary rotations $\\text{CU1}$.\n","\n","$$\n"," \\begin{equation}\n","U1(\\phi) =\n"," \\begin{pmatrix}\n"," 1 & 0 \\\\\n"," 0 & e^{i \\pi \\phi}\n"," \\end{pmatrix}\\, , \\quad\n"," CU1(\\phi) =\n"," \\begin{pmatrix}\n"," 1 & 0 & 0 & 0 \\\\\n"," 0 & 1 & 0 & 0 \\\\\n"," 0 & 0 & 1 & 0 \\\\\n"," 0 & 0 & 0 & e^{i \\pi \\phi}\n"," \\end{pmatrix}\n"," \\end{equation}\n","$$\n","\n","The circuit for the Quantum Fourier transform on three qubits is the following\n","\n","![](images/qft.png \"QFT Circuit\")\n","\n","We can build this circuit in `pytket` by adding gate operations manually:"]},{"cell_type":"markdown","metadata":{},"source":["lets build the QFT for three qubits"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.circuit import Circuit\n","from pytket.circuit.display import render_circuit_jupyter"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qft3_circ = Circuit(3)\n","qft3_circ.H(0)\n","qft3_circ.CU1(0.5, 1, 0)\n","qft3_circ.CU1(0.25, 2, 0)\n","qft3_circ.H(1)\n","qft3_circ.CU1(0.5, 2, 1)\n","qft3_circ.H(2)\n","qft3_circ.SWAP(0, 2)\n","render_circuit_jupyter(qft3_circ)"]},{"cell_type":"markdown","metadata":{},"source":["We can generalise the quantum Fourier transform to $n$ qubits by iterating over the qubits as follows"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def build_qft_circuit(n_qubits: int) -> Circuit:\n","    circ = Circuit(n_qubits, name=\"QFT\")\n","    for i in range(n_qubits):\n","        circ.H(i)\n","        for j in range(i + 1, n_qubits):\n","            circ.CU1(1 / 2 ** (j - i), j, i)\n","    for k in range(0, n_qubits // 2):\n","        circ.SWAP(k, n_qubits - k - 1)\n","    return circ"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qft4_circ: Circuit = build_qft_circuit(4)\n","render_circuit_jupyter(qft4_circ)"]},{"cell_type":"markdown","metadata":{},"source":["Now that we have the generalised circuit we can wrap it up in a `CircBox` which can then be added to another circuit as a subroutine."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.circuit import CircBox"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qft4_box: CircBox = CircBox(qft4_circ)\n","qft_circ = Circuit(4).add_gate(qft4_box, [0, 1, 2, 3])\n","render_circuit_jupyter(qft_circ)"]},{"cell_type":"markdown","metadata":{},"source":["Note how the `CircBox` inherits the name `QFT` from the underlying circuit."]},{"cell_type":"markdown","metadata":{},"source":["Recall that in our phase estimation algorithm we need to use the inverse QFT.\n","\n","$$\n","\\begin{equation}\n","\\text{QFT}^† : \\frac{1}{\\sqrt{N}} \\sum_{k=0}^{N - 1} e^{2 \\pi ijk/N}|k\\rangle \\longmapsto |j\\rangle\\,, \\quad N= 2^n\n","\\end{equation}\n","$$\n","\n","\n","Now that we have the QFT circuit we can obtain the inverse by using `CircBox.dagger`. We can also verify that this is correct by inspecting the circuit inside with `CircBox.get_circuit()`."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["inv_qft4_box = qft4_box.dagger\n","# Explicitly set the name of the `CircBox` to \"QFT†\"\n","inv_qft4_box.circuit_name = \"QFT†\"\n","qft_inv_circ = Circuit(4)\n","qft_inv_circ.add_gate(inv_qft4_box, [0, 1, 2, 3])\n","render_circuit_jupyter(qft_inv_circ)"]},{"cell_type":"markdown","metadata":{},"source":["## Building the Phase Estimation Circuit"]},{"cell_type":"markdown","metadata":{},"source":["We can now define a function to build our entire QPE circuit. We can make this function take a state preparation circuit and a unitary circuit as input as well. The function also has the number of measurement qubits as input which will determine the precision of our phase estimate."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.circuit import QControlBox"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def build_phase_estimation_circuit(\n","    n_measurement_qubits: int, state_prep_circuit: Circuit, unitary_circuit: Circuit\n",") -> Circuit:\n","    # Define a Circuit with a measurement and prep register\n","    qpe_circ: Circuit = Circuit()\n","    n_state_prep_qubits = state_prep_circuit.n_qubits\n","    measurement_register = qpe_circ.add_q_register(\"m\", n_measurement_qubits)\n","    state_prep_register = qpe_circ.add_q_register(\"p\", n_state_prep_qubits)\n","    qpe_circ.add_circuit(state_prep_circuit, list(state_prep_register))\n","\n","    # Create a controlled unitary with a single control qubit\n","    unitary_circuit.name = \"U\"\n","    controlled_u_gate = QControlBox(CircBox(unitary_circuit), 1)\n","\n","    # Add Hadamard gates to every qubit in the measurement register\n","    for m_qubit in measurement_register:\n","        qpe_circ.H(m_qubit)\n","\n","    # Add all (2**n_measurement_qubits - 1) of the controlled unitaries sequentially\n","    for m_qubit in range(n_measurement_qubits):\n","        control_index = n_measurement_qubits - m_qubit - 1\n","        control_qubit = [measurement_register[control_index]]\n","        for _ in range(2**m_qubit):\n","            qpe_circ.add_gate(\n","                controlled_u_gate, control_qubit + list(state_prep_register)\n","            )\n","\n","    # Finally, append the inverse qft and measure the qubits\n","    qft_box = CircBox(build_qft_circuit(n_measurement_qubits))\n","    inverse_qft_box = qft_box.dagger\n","    inverse_qft_box.circuit_name = \"QFT†\"\n","    qpe_circ.add_gate(inverse_qft_box, list(measurement_register))\n","    qpe_circ.measure_register(measurement_register, \"c\")\n","    return qpe_circ"]},{"cell_type":"markdown","metadata":{},"source":["## Phase Estimation with a Trivial Eigenstate\n","\n","Lets test our circuit construction by preparing a trivial $|1\\rangle$ eigenstate of the $\\text{U1}$ gate. We can then see if our phase estimation circuit returns the expected eigenvalue."]},{"cell_type":"markdown","metadata":{},"source":["$$\n","\\begin{equation}\n","U1(\\phi)|1\\rangle = e^{i \\pi \\phi}|1\\rangle = e^{2 \\pi i \\theta} |1\\rangle \\implies \\theta = \\frac{\\phi}{2}\n","\\end{equation}\n","$$\n","\n","So we expect that our ideal phase $\\theta$ will be half the input angle $\\phi$ to our $U1$ gate."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["prep_circuit = Circuit(1).X(0)  # prepare the |1> eigenstate of U1"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["input_angle = 0.73  # angle as number of half turns"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["unitary_circuit = Circuit(1).U1(input_angle, 0)  # Base unitary for controlled U ops"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qpe_circ_trivial = build_phase_estimation_circuit(\n","    4, state_prep_circuit=prep_circuit, unitary_circuit=unitary_circuit\n",")"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["render_circuit_jupyter(qpe_circ_trivial)"]},{"cell_type":"markdown","metadata":{},"source":["Lets use the noiseless `AerBackend` simulator to run our phase estimation circuit."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.extensions.qiskit import AerBackend"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["backend = AerBackend()"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["compiled_circ = backend.get_compiled_circuit(qpe_circ_trivial)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["n_shots = 1000\n","result = backend.run_circuit(compiled_circ, n_shots)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["print(result.get_counts())"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.backends.backendresult import BackendResult\n","import matplotlib.pyplot as plt"]},{"cell_type":"markdown","metadata":{},"source":["plotting function for QPE Notebook"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def plot_qpe_results(\n","    sim_result: BackendResult,\n","    n_strings: int = 4,\n","    dark_mode: bool = False,\n","    y_limit: int = 1000,\n",") -> None:\n","    \"\"\"\n","    Plots results in a barchart given a BackendResult. the number of stings displayed\n","    can be specified with the n_strings argument.\n","    \"\"\"\n","    counts_dict = sim_result.get_counts()\n","    sorted_shots = counts_dict.most_common()\n","    n_most_common_strings = sorted_shots[:n_strings]\n","    x_axis_values = [str(entry[0]) for entry in n_most_common_strings]  # basis states\n","    y_axis_values = [entry[1] for entry in n_most_common_strings]  # counts\n","    if dark_mode:\n","        plt.style.use(\"dark_background\")\n","    fig = plt.figure()\n","    ax = fig.add_axes((0, 0, 0.75, 0.5))\n","    color_list = [\"orange\"] * (len(x_axis_values))\n","    ax.bar(\n","        x=x_axis_values,\n","        height=y_axis_values,\n","        color=color_list,\n","    )\n","    ax.set_title(label=\"Results\")\n","    plt.ylim([0, y_limit])\n","    plt.xlabel(\"Basis State\")\n","    plt.ylabel(\"Number of Shots\")\n","    plt.show()"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["plot_qpe_results(result, y_limit=int(1.2 * n_shots))"]},{"cell_type":"markdown","metadata":{},"source":["As expected we see one outcome with high probability. Lets now extract our approximation of $\\theta$ from our output bitstrings.\n","\n","suppose the $j$ is an integer representation of our most commonly measured bitstring."]},{"cell_type":"markdown","metadata":{},"source":["$$\n","\\begin{equation}\n","\\theta_{estimate} = \\frac{j}{N}\n","\\end{equation}\n","$$"]},{"cell_type":"markdown","metadata":{},"source":["Here $N = 2 ^m$ where $m$ is the number of measurement qubits."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.backends.backendresult import BackendResult"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def single_phase_from_backendresult(result: BackendResult) -> float:\n","    # Extract most common measurement outcome\n","    basis_state = result.get_counts().most_common()[0][0]\n","    bitstring = \"\".join([str(bit) for bit in basis_state])\n","    integer_j = int(bitstring, 2)\n","\n","    # Calculate theta estimate\n","    return integer_j / (2 ** len(bitstring))"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["theta = single_phase_from_backendresult(result)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["print(theta)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["print(input_angle / 2)"]},{"cell_type":"markdown","metadata":{},"source":["Our output is close to half our input angle $\\phi$ as expected. Lets calculate our error $E$ to three decimal places."]},{"cell_type":"markdown","metadata":{},"source":["$$\n","\\begin{equation}\n","E = |\\phi - 2 \\, \\theta_{estimate}|\n","\\end{equation}\n","$$"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["error = round(abs(input_angle - (2 * theta)), 3)\n","print(error)"]},{"cell_type":"markdown","metadata":{},"source":["## Phase Estimation with Time Evolution"]},{"cell_type":"markdown","metadata":{},"source":["In the phase estimation algorithm we repeatedly perform controlled unitary operations. In the textbook variant of QPE presented here, the number of controlled unitaries will be $2^m - 1$ where $m$ is the number of measurement qubits."]},{"cell_type":"markdown","metadata":{},"source":["In the example above we've shown a trivial instance of QPE where we know the exact phase in advance. For more realistic applications of QPE we will have some non-trivial state preparation required.\n","\n","For chemistry or condensed matter physics $U$ typically be the time evolution operator $U(t) = e^{- i H t}$ where $H$ is the problem Hamiltonian.\n","Suppose that we had the following decomposition for $H$ in terms of Pauli strings $P_j$ and complex coefficients $\\alpha_j$.\n","\n","$$\n","\\begin{equation}\n","H = \\sum_j \\alpha_j P_j\\,, \\quad \\, P_j \\in \\{I, \\,X, \\,Y, \\,Z\\}^{\\otimes n}\n","\\end{equation}\n","$$\n","\n","Here the term Pauli strings refers to tensor products of Pauli operators. These strings form an orthonormal basis for $2^n \\times 2^n$ matrices."]},{"cell_type":"markdown","metadata":{},"source":["If we have a Hamiltonian in the form above, we can then implement $U(t)$ as a sequence of Pauli gadget circuits. We can do this with the [PauliExpBox](https://tket.quantinuum.com/api-docs/circuit.html#pytket.circuit.PauliExpBox) construct in pytket. For more on `PauliExpBox` see the [user manual](https://tket.quantinuum.com/user-manual/manual_circuit.html#pauli-exponential-boxes)."]},{"cell_type":"markdown","metadata":{},"source":["Once we have a circuit to implement our time evolution operator $U(t)$, we can construct the controlled $U(t)$ operations using [QControlBox](https://tket.quantinuum.com/api-docs/circuit.html#pytket.circuit.QControlBox). If our base unitary is a sequence of `PauliExpBox`(es) then there is some structure we can exploit to simplify our circuit. See this [blog post](https://tket.quantinuum.com/blog/posts/controlled_gates/) on [ConjugationBox](https://tket.quantinuum.com/api-docs/circuit.html#pytket.circuit.ConjugationBox) for more."]},{"cell_type":"markdown","metadata":{},"source":["As an exercise, try to use phase estimation to calculate the ground state of diatomic hydrogen $H_2$."]},{"cell_type":"markdown","metadata":{},"source":["## Suggestions for further reading\n","\n","* Quantinuum paper on Bayesian phase estimation -> https://arxiv.org/pdf/2306.16608.pdf\n","* Blog post on `ConjugationBox` (efficient circuits for controlled gates) -> https://tket.quantinuum.com/blog/posts/controlled_gates/"]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.4"}},"nbformat":4,"nbformat_minor":2}