{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Quantum phase estimation\n", "\n", "In this tutorial, we will look at various methods for estimating eigenenergies $E_j$ of a Hamiltonian $H$, or equivalently (for sufficiently small $t$), eigenphases $\\phi_j=-E_jt$ of a unitary operator $U=e^{-iHt}$. \n", "\n", "1. Multi-ancilla QPE\n", " 1. Controlled unitary circuits\n", " 1. Unitaries on multiple ancillas\n", " 1. The quantum Fourier transform\n", "1. Single-ancilla QPE\n", " 1. Robust phase estimation\n", " 1. Prony estimation\n", " 1. Scaling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Multi-ancilla QPE\n", "Let us take a very simple Hamiltonian to begin:\n", "$$ H = H_0 ZI + H_1 IZ + H_2 ZZ $$\n", "\n", "Though not particularly interesting, this Hamiltonian has the advantage of being easy to solve and characterize\n", "\n", "**Problem 1:** Write down the eigenvalues and corresponding eigenvectors of $H$ in terms of the three parameters $H_0$, $H_1$, and $H_2$. What are the equivalent eigenvalues and eigenvectors of $U$?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def get_eigenvalues(H0, H1, H2):\n", " return [H0+H1+H2, H0-H1-H2, -H0+H1-H2, -H0-H1+H2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Controlled unitary circuits\n", "Quantum phase estimation requires that we perform $U=e^{-iHt}$ conditional on an ancilla qubit. That is, that we apply the unitary\n", "$$ I\\oplus U = \\left(\\begin{array}{cc}1 & 0 \\\\ 0 & U\\end{array}\\right) $$\n", "Let's first try implementing $U$ by itself. Recall from previous lectures that we can apply $U$ via Trotterization --- separately applying $e^{-iH_0ZIt}$, $e^{-iH_1IZt}$, and $e^{-iH_2ZZt}$.\n", "\n", "**Problem 2:** Write circuits to implement the above three unitary operations using rz rotations and CNOT gates" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import cirq\n", "\n", "qubits = [cirq.GridQubit(0,0), cirq.GridQubit(0,1)]\n", "\n", "def time_evolution_H0(H0, t):\n", " # Note the factor of -2, as cirq.rz(theta) implements exp(-i*Z*theta/2).\n", " yield cirq.rz(2*H0*t).on(qubits[0])\n", " \n", "def time_evolution_H1(H1, t):\n", " ### Insert your code here!\n", " \n", "def time_evolution_H2(H2, t):\n", " ### Insert your code here!\n", "\n", "# Plug in some values for our parameters \n", "H0 = 0.7\n", "H1 = 1.2\n", "H2 = 0.8\n", "t = 0.5\n", "circuit = cirq.Circuit([time_evolution_H0(H0, t),\n", " time_evolution_H1(H1, t),\n", " time_evolution_H2(H2, t)])\n", "print(circuit)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Problem 3:** Use cirq.unitary to check the values of the generated unitary against your answers to problem 1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "### Insert your code here!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The simplest method to implement a circuit controlled by an ancilla qubit is simply to control each and every one of the gates. However, this is relatively high-cost, as we must then decompose three-qubit gates in terms of two-qubit gates.\n", "\n", "A more economical method is to realize that, in the above circuit, chains of CNOT gates will cancel out if the Z-rotation gate between them is not implemented. This implies that we need only control the Z gates on the above qubits for QPE (a much more practical implementation).\n", "\n", "One technical point of the above - the native gate in superconducting qubit hardware is the controlled-phase gate cirq.CZ**\\phi (with a diagonal $(1,1,1,e^{i\\phi*\\pi})$), and not the controlled-rz gate (with a diagonal $(1,1,e^{-i\\phi},e^{i\\phi})$). We need to make sure we translate between these so as to not pick up spurious phases.\n", "\n", "**Problem 4:** Write a function that implements controlled-rz in terms of the cirq.CZ gate" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def crz(control, target, phase):\n", " ### Insert your code here!\n", " \n", "test_unitary = cirq.unitary(cirq.Circuit(crz(qubits[0],qubits[1],np.pi/2)))\n", "# What do you expect to see here? Does your code match your expectations?\n", "print((test_unitary / test_unitary[0,0]).round(2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's use this to make our circuit controlled\n", "\n", "**Problem 5:** Implement a controlled version of the above circuits." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ancilla = cirq.GridQubit(1,0)\n", "\n", "def controlled_time_evolution_H0(ancilla, H0, t):\n", " # Check factors of 2 here again!\n", " yield crz(ancilla, qubits[0], H0*t)\n", " \n", "def controlled_time_evolution_H1(ancilla, H1, t):\n", " ### Insert your code here!\n", " \n", "def controlled_time_evolution_H2(ancilla, H2, t):\n", " ### Insert your code here!\n", " \n", "circuit = cirq.Circuit([controlled_time_evolution_H0(ancilla, H0, t),\n", " controlled_time_evolution_H1(ancilla, H1, t),\n", " controlled_time_evolution_H2(ancilla, H2, t)])\n", "print(circuit)\n", "\n", "# Check angles again\n", "unitary = cirq.unitary(circuit)\n", "print()\n", "print(np.round([np.angle(unitary[j,j]/unitary[0,0]) / -t for j in range(8)], 2))\n", "print(np.round(get_eigenvalues(H0, H1, H2), 2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multiple-ancilla evolution\n", "\n", "Traditionally QPE is achieved by using a register of $L\\sim\\log(1/\\epsilon)$ ancilla qubits labeled by an index $l=0,\\ldots,L-1$, and performing controlled-$e^{i2^lHt}$ with each ancilla.\n", "\n", "Luckily for us, the above code is written so as to make the extension quite easy.\n", "\n", "**Problem 6:** write a function using the above code to implement controlled $e^{i2^lHt}$ for a register of $l=0,\\ldots,L-1$ ancilla qubits." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "L=2\n", "ancillas = [cirq.GridQubit(1,l) for l in range(L)]\n", "\n", "def controlled_time_evolution_multi_ancilla(ancillas, H0, H1, H2, t):\n", " for l, ancilla in enumerate(ancillas):\n", " ### insert your code here!\n", " \n", "circuit = cirq.Circuit(controlled_time_evolution_multi_ancilla(ancillas, H0, H1, H2, t))\n", "print(circuit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note here that for ease of checking the code / reading the result you might have and implemented $e^{i2^lHt}$ in a single shot. This is known as fast-forwarding, and for a general Hamiltonian is known to not be possible --- to implement $e^{i2^lHt}$ with constant error in both $l$ and $t$ will typically take time $2^lt$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The quantum Fourier transform\n", "\n", "Traditionally to process QPE data, we would perform the inverse quantum Fourier transform on the ancilla qubits, and then read them out in the computational basis. To do this, we need to implement the quantum Fourier transform.\n", "\n", "**Problem 7:** Write a function to implement the quantum Fourier transform on a register of qubits" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def iqft(register):\n", " ### Insert your code here\n", " \n", "L=5\n", "ancillas = [cirq.GridQubit(1,l) for l in range(L)]\n", "circuit = cirq.Circuit(iqft(ancillas))\n", "print(circuit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we put everything together, we need to prepare the ancilla and system states.\n", "\n", "The easiest ancilla state preparation is to prepare $\\frac{1}{\\sqrt{2}^L}\\sum_l|l\\rangle$.\n", "\n", "For the system preparation, we have the freedom to prepare whatever we want - let's prepare a superposition of the $|01\\rangle$ and $|10\\rangle$ states.\n", "\n", "**Problem 7:** Write a circuit to prepare the above states" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def state_preparation_ancilla(register):\n", " ### Insert your code here\n", " \n", "def state_preparation_system():\n", " ### Insert your code here\n", "\n", "circuit = cirq.Circuit(state_preparation_system())\n", "print(circuit)\n", "print()\n", " \n", "circuit = cirq.Circuit(state_preparation_ancilla(ancillas))\n", "print(circuit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we go any further, let's check that our QFT returns the right result on the prepared state (when no phase is accumulated). If the following does not return only 0, we have made a mistake!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "circuit = cirq.Circuit(state_preparation_ancilla(ancillas))\n", "circuit.append(iqft(ancillas))\n", "circuit.append(cirq.measure(*ancillas, key='phase'))\n", "simulator = cirq.Simulator()\n", "results = simulator.run(circuit, repetitions = 1000)\n", "print(results.histogram(key='phase'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's put it all together!\n", "\n", "**Problem 8:** combine the above code and simulate quantum phase estimation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "L=5\n", "H0 = 0.125 * 2 * np.pi\n", "H1 = 0.4 * 2 * np.pi\n", "H2 = 0.3 * 2 * np.pi\n", "ancillas = [cirq.GridQubit(1,l) for l in range(L)] \n", "\n", "### Insert your code here\n", "\n", "simulator = cirq.Simulator()\n", "result = simulator.run(circuit, repetitions=10000)\n", "hist = result.histogram(key='phase')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's plot the result!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "plt.figure(figsize=(11,7))\n", "histx = [x / 2**L * 2 * np.pi for x in hist]\n", "histy = [hist[x] for x in hist]\n", "plt.plot(histx, histy, 'k+', markersize=15,markeredgewidth=3)\n", "plt.vlines([-x % (2 * np.pi) for x in get_eigenvalues(H0,H1,H2)],0, 10000, linestyle='dashed')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Single-ancilla QPE\n", "\n", "Though the above performs quite well, the use of multiple ancilla qubits adds extra overhead we can do away with. In effect, we have no need to perform the quantum Fourier transform - if we can extract the phase data accumulated on a single ancilla qubit, we can convert this into a signal classically.\n", "\n", "To see this, let's go back to the single-ancilla phase estimation circuit, and investigate the phase accumulated between the two states of the ancilla when our system is in an eigenstate - this can be found as the expectation value of $X-iY$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "H0 = 1\n", "H1 = 0.5\n", "H2 = 0.8\n", "\n", "M = np.kron(np.array([[0, 0], [2, 0]]),np.identity(4))\n", "\n", "\n", "phase_values = []\n", "time_values = np.linspace(0, 20, 101)\n", "for t in time_values:\n", " circuit = cirq.Circuit([cirq.H(ancilla),\n", " controlled_time_evolution_H0(ancilla, H0, t),\n", " controlled_time_evolution_H1(ancilla, H1, t),\n", " controlled_time_evolution_H2(ancilla, H2, t)])\n", " \n", " #simulator.simulate extracts the wavefunction at the end of the circuit\n", " result = simulator.simulate(circuit, qubit_order = [ancilla, *qubits])\n", " phase_values.append(result.final_state.conj().T @ M @ result.final_state)\n", "\n", "plt.figure(figsize=(11,7))\n", "plt.plot(time_values, np.real(phase_values), 'k+')\n", "plt.plot(time_values, np.imag(phase_values), 'r+')\n", "plt.plot(time_values, np.cos(time_values*(H0+H1+H2)), 'k--')\n", "plt.plot(time_values, np.sin(time_values*(H0+H1+H2)), 'r--')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we need a piece of code that estimates $\\langle X-iY\\rangle$. However, this is not a Hermitian operator, so we need to write $\\langle X-iY\\rangle = \\langle X\\rangle - i\\langle Y\\rangle$, and estimate these quantities separately via partial state tomography.\n", "\n", "**Problem 9:** Write code to estimate $\\langle X\\rangle$ and $\\langle Y\\rangle$ using two separate circuits." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# We will replace this function soon\n", "def state_preparation():\n", " return []\n", "\n", "def get_expectation_X(state_preparation, H0, H1, H2, t, num_repetitions):\n", " ### Insert your code here\n", " \n", "def get_expectation_Y(state_preparation, H0, H1, H2, t, num_repetitions):\n", " ### Insert your code here\n", "\n", "num_repetitions = 100\n", "new_phase_values = [get_expectation_X(state_preparation, H0, H1, H2, t, num_repetitions) -\n", " 1j * get_expectation_Y(state_preparation, H0, H1, H2, t, num_repetitions)\n", " for t in time_values]\n", "\n", "plt.figure(figsize=(11,7))\n", "plt.plot(time_values, np.real(new_phase_values), 'k+')\n", "plt.plot(time_values, np.imag(new_phase_values), 'r+')\n", "plt.plot(time_values, np.cos(time_values*(H0+H1+H2)), 'k--')\n", "plt.plot(time_values, np.sin(time_values*(H0+H1+H2)), 'r--')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Robust phase estimation\n", "\n", "If you look closely, you can see some deviations between the expected signal and the obtained points. This is sampling noise --- it's proportional to the square root of the number of repetitions made. This noise will propagate through to our estimation of the phase.\n", "\n", "As we started in an eigenstate, our signal is just a single wave - this means we can estimate the phase from a single point!\n", "\n", "**Problem 10:** Write a function to estimate the phase from a single point of the signal." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def phase_estimation_single_point(time, signal):\n", " ### Insert your code here\n", "\n", "print(\"Estimating phase from point at time {}\".format(time_values[1]))\n", "print(\"Estimated phase = {}\".format(phase_estimation_single_point(time_values[1], new_phase_values[1])))\n", "print(\"True phase = {}\".format(H0+H1+H2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Probably you're seeing quite a large error here. This is because with only 100 repetitions, our estimate of the signal has an error on the order of $\\sqrt{1/100}=0.1$, and this error is divided by the value of the time that we estimated at --- so the standard error of the estimate is on the order of $0.5$.\n", "\n", "Let's investigate what happens when we tune the knobs of the number of measurements and the time-step that we estimate at:\n", "\n", "**Problem 11:** Repeat the above procedure for different values of the number of repetitions and of the time. Plot the error in the estimation for both on a log-log scale." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "num_repetitions_range = [100, 200, 500, 1000, 2000, 5000, 10000, 100000]\n", "error_increasing_repetitions = []\n", "true_phase = H0+H1+H2\n", "t = 0.2\n", "for num_repetitions in num_repetitions_range:\n", " ### Insert your code here\n", "\n", "num_repetitions = 10000\n", "time_range = [0.1, 0.2, 0.3, 0.5, 1, 1.5, 2, 3, 5, 6, 10]\n", "error_increasing_time = []\n", "for t in time_range:\n", " ### Insert your code here\n", " \n", "plt.figure(figsize=(11,7))\n", "plt.plot(num_repetitions_range, error_increasing_repetitions, 'k+', markersize=15, markeredgewidth=3, label='data')\n", "plt.plot(num_repetitions_range, np.array(num_repetitions_range)**(-1/2)*2, 'r--', linewidth=2, label='approximate fit')\n", "plt.xlabel('number of repetitions', fontsize=18)\n", "plt.ylabel('error', fontsize=18)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.legend(fontsize=18)\n", "\n", "plt.figure(figsize=(11,7))\n", "plt.plot(time_range, error_increasing_time, 'k+', markersize=15, markeredgewidth=3, label='data')\n", "plt.plot(time_range, np.array(time_range)**(-1)*0.01, 'r--', linewidth=2, label='approximate fit')\n", "plt.xlabel('time', fontsize=18)\n", "plt.ylabel('error', fontsize=18)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.legend(fontsize=18)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I've thrown in some guide lines for the eye (to confirm this you should repeat the calculation of each point many times and average the error). Hopefully you see that the scaling with the number of repetitions works as expected... But there is a horrible jump in the scaling with the time!\n", "\n", "What is going on here?\n", "\n", "The answer is that our estimate of $E_jt$ is only accurate modulo $2\\pi$ - whenever $E_jt>2\\pi$, we don't see that $e^{iE_jt}$ is wrapped around the complex unit circle. This places a limit on the value of $t$ we can use for this estimation. This is a bit of a shame, as our scaling with $t$ is twice as good as our scaling with the number of repetitions.\n", "\n", "There are many methods to circumvent this issue. The basic idea is to use a small value of $t$ to form a rough estimate of $E_jt$, then use this to keep track of the modulo $2\\pi$ as $t$ grows larger and larger. In this tutorial, we will use the method known as 'robust phase estimation'.\n", "\n", "Robust phase estimation requires that we perform separate estimates $\\tilde{F_j}^{(l)}$ of $E_jt\\mod 2\\pi$ for $t=2^lt_0$ for $l=0,\\ldots,L$. Then, at each order $j>0$, we estimate $\\tilde{E_j}^{(l)}$ as the number closest to $\\tilde{E_j}^{(l-1)}$ such that \n", "\n", "$$\\tilde{E_j}^{(l)}\\mod (2\\pi / 2^lt_0) = \\tilde{F_j}^{(l)}.$$\n", "Let's try this!\n", "\n", "**Problem 12:** Write a function to generate a set of estimators $\\tilde{F_j}^{(l)}$ from a set of input experiments.\n", "\n", "**Problem 13:** Write an estimator for $\\tilde{E_j}^{(l)}$ using the above set of estimators $\\tilde{F_j}^{(l)}$. Test the scaling you obtain!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def phase_estimator_modulo_2pi(exponentially_sparse_signal, t0):\n", " ### Insert your code here\n", " \n", "def robust_phase_estimation(exponentially_sparse_signal, t0):\n", " ### Insert your code here\n", "\n", "num_repetitions = 10000\n", "t0 = 0.5\n", "J = 10\n", "time_range = [t0 * 2**j for j in range(0, J+1)]\n", "exponentially_sparse_signal = [get_expectation_X(state_preparation, H0, H1, H2, t, num_repetitions) -\n", " 1j * get_expectation_Y(state_preparation, H0, H1, H2, t, num_repetitions)\n", " for t in time_range]\n", "\n", "error_increasing_time = robust_phase_estimation(exponentially_sparse_signal, t0)\n", "\n", "plt.figure(figsize=(11,7))\n", "plt.plot(time_range, error_increasing_time, 'k+', markersize=15, markeredgewidth=3, label='data')\n", "plt.plot(time_range, np.array(time_range)**(-1)*0.01, 'r--', linewidth=2, label='approximate fit')\n", "plt.xlabel('time', fontsize=18)\n", "plt.ylabel('error', fontsize=18)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.legend(fontsize=18)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It turns out that, as stated, the above method is not completely optimal in time complexity. We need to be careful to choose the number of repetitions at later values of $j$ to be low enough that we don't completely swamp the total number of $e^{iHt}$. But this is probably beyond what we have time for in this tutorial." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prony phase estimation\n", "\n", "QPE signal processing gets more complicated when we use an initial state that is not an eigenstate of our system:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "num_repetitions = 1000\n", "time_values = np.linspace(0, 100, 101)\n", "new_phase_values = [get_expectation_X(state_preparation_system, H0, H1, H2, t, num_repetitions) -\n", " 1j * get_expectation_Y(state_preparation_system, H0, H1, H2, t, num_repetitions)\n", " for t in time_values]\n", "\n", "plt.figure(figsize=(11,7))\n", "plt.plot(time_values, np.real(new_phase_values), 'k+-')\n", "plt.plot(time_values, np.imag(new_phase_values), 'r+-')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is going on here? We have two frequencies superimposed, like two notes of a chord. Robust phase estimation will fail here, as the estimators $\\tilde{F_j}^{(l)}$ explicitly assume that the input signal is a single frequency.\n", "\n", "One way to extract multiple frequencies from this wave is via Prony's method, which can be found in the openfermion package." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from openfermion import prony\n", "amplitudes, phases = prony(new_phase_values)\n", "plt.figure(figsize=(11,7))\n", "\n", "plt.plot(np.angle(phases) % (2*np.pi), np.abs(amplitudes), 'k+', markersize=15, markeredgewidth=2)\n", "plt.vlines([x % (2 * np.pi) for x in get_eigenvalues(H0,H1,H2)],0, 1, linestyle='dashed')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is not the only way of estimating single frequencies from a more complicated signal; applying these methods to quantum phase estimation is very much a work in progress.\n", "\n", "But, to finish this tutorial see if you can do better than yours truly at using Prony's method to estimate eigenvalues of a Hamiltonian.\n", "\n", "**Problem 14:** Investigate the scaling of Prony's method with the number of repetitions and maximum value of $t$. How well can you do?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.7.5rc1" } }, "nbformat": 4, "nbformat_minor": 2 }