spectral_measurement.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Algorithm for Spectral Measurement with Lower Gate Count\n",
    "\n",
    "This tutorial shows how to implement the algorithm introduced in the following paper:\n",
    "\n",
    "**Quantum Algorithm for Spectral Measurement with Lower Gate Count**\n",
    "by David Poulin, Alexei Kitaev, Damian S. Steiger, Matthew B. Hastings, Matthias Troyer\n",
    "[Phys. Rev. Lett. 121, 010501 (2018)](https://doi.org/10.1103/PhysRevLett.121.010501)\n",
    "([arXiv:1711.11025](https://arxiv.org/abs/1711.11025))\n",
    "\n",
    "For details please see the above paper. The implementation in ProjectQ is discussed in the PhD thesis of Damian S. Steiger (soon available online). A more detailed discussion will be uploaded soon.\n",
    "Here we only show a small part of the paper, namely the implementation of W and how it can be used with iterative phase estimation to obtain eigenvalues and eigenstates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "import math\n",
    "\n",
    "import scipy.sparse.linalg as spsl\n",
    "\n",
    "import projectq\n",
    "from projectq.backends import Simulator\n",
    "from projectq.meta import Compute, Control, Dagger, Uncompute\n",
    "from projectq.ops import All, H, Measure, Ph, QubitOperator, R, StatePreparation, X, Z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Let's use a simple Hamiltonian acting on 3 qubits for which we want to know the eigenvalues:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_qubits = 3\n",
    "\n",
    "hamiltonian = QubitOperator()\n",
    "hamiltonian += QubitOperator(\"X0\", -1/12.)\n",
    "hamiltonian += QubitOperator(\"X1\", -1/12.)\n",
    "hamiltonian += QubitOperator(\"X2\", -1/12.)\n",
    "hamiltonian += QubitOperator(\"Z0 Z1\", -1/12.)\n",
    "hamiltonian += QubitOperator(\"Z0 Z2\", -1/12.)\n",
    "hamiltonian += QubitOperator(\"\", 7/12.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this quantum algorithm, we need to normalize the hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "hamiltonian_norm = 0.\n",
    "for term in hamiltonian.terms:\n",
    "    hamiltonian_norm += abs(hamiltonian.terms[term])\n",
    "normalized_hamiltonian = deepcopy(hamiltonian)\n",
    "normalized_hamiltonian /= hamiltonian_norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1.** We implement a short helper function which uses the ProjectQ simulator to numerically calculate some eigenvalues and eigenvectors of Hamiltonians stored in ProjectQ's `QubitOperator` in order to check our implemenation of the quantum algorithm. This function is particularly fast because it doesn't need to build the matrix of the hamiltonian but instead uses implicit matrix vector multiplication by using our simulator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_eigenvalue_and_eigenvector(n_sites, hamiltonian, k, which='SA'):\n",
    "    \"\"\"\n",
    "    Returns k eigenvalues and eigenvectors of the hamiltonian.\n",
    "    \n",
    "    Args:\n",
    "        n_sites(int): Number of qubits/sites in the hamiltonian\n",
    "        hamiltonian(QubitOperator): QubitOperator representating the Hamiltonian\n",
    "        k: num of eigenvalue and eigenvector pairs (see spsl.eigsh k)\n",
    "        which: see spsl.eigsh which\n",
    "    \n",
    "    \"\"\"\n",
    "    def mv(v):\n",
    "        eng = projectq.MainEngine(backend=Simulator(), engine_list=[])\n",
    "        qureg = eng.allocate_qureg(n_sites)\n",
    "        eng.flush()\n",
    "        eng.backend.set_wavefunction(v, qureg)\n",
    "        eng.backend.apply_qubit_operator(hamiltonian, qureg)\n",
    "        order, output = deepcopy(eng.backend.cheat())\n",
    "        for i in order:\n",
    "            assert i == order[i]\n",
    "        eng.backend.set_wavefunction([1]+[0]*(2**n_sites-1), qureg)\n",
    "        return output\n",
    "\n",
    "    A = spsl.LinearOperator((2**n_sites,2**n_sites), matvec=mv)\n",
    "\n",
    "    eigenvalues, eigenvectormatrix = spsl.eigsh(A, k=k, which=which)\n",
    "    eigenvectors = []\n",
    "    for i in range(k):\n",
    "        eigenvectors.append(list(eigenvectormatrix[:, i]))\n",
    "    return eigenvalues, eigenvectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use this function to find the 4 lowest eigenstates of the normalized hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.29217007 0.36634371 0.5        0.57417364]\n"
     ]
    }
   ],
   "source": [
    "eigenvalues, eigenvectors = get_eigenvalue_and_eigenvector(\n",
    "    n_sites=num_qubits,\n",
    "    hamiltonian=normalized_hamiltonian,\n",
    "    k=4)\n",
    "print(eigenvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the eigenvalues are all positive as required (otherwise increase identity term in hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2.** Let's define the W operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def W(eng, individual_terms, initial_wavefunction, ancilla_qubits, system_qubits):\n",
    "    \"\"\"\n",
    "    Applies the W operator as defined in arXiv:1711.11025.\n",
    "    \n",
    "    Args:\n",
    "        eng(MainEngine): compiler engine\n",
    "        individual_terms(list<QubitOperator>): list of individual unitary\n",
    "                                               QubitOperators. It applies\n",
    "                                               individual_terms[0] if ancilla\n",
    "                                               qubits are in state |0> where\n",
    "                                               ancilla_qubits[0] is the least\n",
    "                                               significant bit.\n",
    "        initial_wavefunction: Initial wavefunction of the ancilla qubits\n",
    "        ancilla_qubits(Qureg): ancilla quantum register in state |0>\n",
    "        system_qubits(Qureg): system quantum register\n",
    "    \"\"\"\n",
    "    # Apply V:\n",
    "    for ancilla_state in range(len(individual_terms)):\n",
    "        with Compute(eng):\n",
    "            for bit_pos in range(len(ancilla_qubits)):\n",
    "                if not (ancilla_state >> bit_pos) & 1:\n",
    "                    X | ancilla_qubits[bit_pos]\n",
    "        with Control(eng, ancilla_qubits):\n",
    "            individual_terms[ancilla_state] | system_qubits\n",
    "        Uncompute(eng)\n",
    "    # Apply S: 1) Apply B^dagger\n",
    "    with Compute(eng):\n",
    "        with Dagger(eng):\n",
    "            StatePreparation(initial_wavefunction) | ancilla_qubits\n",
    "    # Apply S: 2) Apply I-2|0><0|\n",
    "    with Compute(eng):\n",
    "        All(X) | ancilla_qubits\n",
    "    with Control(eng, ancilla_qubits[:-1]):\n",
    "        Z | ancilla_qubits[-1]\n",
    "    Uncompute(eng)\n",
    "    # Apply S: 3) Apply B\n",
    "    Uncompute(eng)\n",
    "    # Could also be omitted and added when calculating the eigenvalues:\n",
    "    Ph(math.pi) | system_qubits[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.** For testing this algorithm, let's initialize the qubits in a superposition state of the lowest and second lowest eigenstate of the hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "eng = projectq.MainEngine()\n",
    "system_qubits = eng.allocate_qureg(num_qubits)\n",
    "\n",
    "# Create a normalized equal superposition of the two eigenstates for numerical testing:\n",
    "initial_state_norm =0.\n",
    "initial_state = [i+j for i,j in zip(eigenvectors[0], eigenvectors[1])]\n",
    "for amplitude in initial_state:\n",
    "    initial_state_norm += abs(amplitude)**2\n",
    "normalized_initial_state = [amp / math.sqrt(initial_state_norm) for amp in initial_state]\n",
    "\n",
    "#initialize system qubits in this state:\n",
    "StatePreparation(normalized_initial_state) | system_qubits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4.** Split the normalized_hamiltonian into individual terms and build the wavefunction for the ancilla qubits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "individual_terms = []\n",
    "initial_ancilla_wavefunction = []\n",
    "for term in normalized_hamiltonian.terms:\n",
    "    coefficient = normalized_hamiltonian.terms[term]\n",
    "    initial_ancilla_wavefunction.append(math.sqrt(abs(coefficient)))\n",
    "    if coefficient < 0:\n",
    "        individual_terms.append(QubitOperator(term, -1))\n",
    "    else:\n",
    "        individual_terms.append(QubitOperator(term))\n",
    "\n",
    "# Calculate the number of ancilla qubits required and pad\n",
    "# the ancilla wavefunction with zeros:\n",
    "num_ancilla_qubits = int(math.ceil(math.log(len(individual_terms), 2)))\n",
    "required_padding = 2**num_ancilla_qubits - len(initial_ancilla_wavefunction)\n",
    "initial_ancilla_wavefunction.extend([0]*required_padding)\n",
    "\n",
    "# Initialize ancillas by applying B\n",
    "ancillas = eng.allocate_qureg(num_ancilla_qubits)\n",
    "StatePreparation(initial_ancilla_wavefunction) | ancillas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.** Perform an iterative phase estimation of the unitary W to collapse to one of the eigenvalues of the `normalized_hamiltonian`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Semiclassical iterative phase estimation\n",
    "bits_of_precision = 8\n",
    "pe_ancilla = eng.allocate_qubit()\n",
    "\n",
    "measurements = [0] * bits_of_precision\n",
    "\n",
    "for k in range(bits_of_precision):\n",
    "    H | pe_ancilla\n",
    "    with Control(eng, pe_ancilla):\n",
    "        for i in range(2**(bits_of_precision-k-1)):\n",
    "            W(eng=eng,\n",
    "              individual_terms=individual_terms,\n",
    "              initial_wavefunction=initial_ancilla_wavefunction,\n",
    "              ancilla_qubits=ancillas,\n",
    "              system_qubits=system_qubits)\n",
    "\n",
    "    #inverse QFT using one qubit\n",
    "    for i in range(k):\n",
    "        if measurements[i]:\n",
    "            R(-math.pi/(1 << (k - i))) | pe_ancilla\n",
    "\n",
    "    H | pe_ancilla\n",
    "    Measure | pe_ancilla\n",
    "    eng.flush()\n",
    "    measurements[k] = int(pe_ancilla)\n",
    "    # put the ancilla in state |0> again\n",
    "    if measurements[k]:\n",
    "        X | pe_ancilla\n",
    "\n",
    "est_phase = sum(\n",
    "    [(measurements[bits_of_precision - 1 - i]*1. / (1 << (i + 1)))\n",
    "     for i in range(bits_of_precision)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We measured 0.203125 corresponding to energy 0.290284677254\n"
     ]
    }
   ],
   "source": [
    "print(\"We measured {} corresponding to energy {}\".format(est_phase, math.cos(2*math.pi*est_phase)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6.** We measured the lowest eigenstate. You can verify that this happens with 50% probability as we chose our initial state to have 50% overlap with the ground state. As the paper notes, the `system_qubits` are not in an eigenstate and one can easily test that using our simulator to get the energy of the current state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.33236578253447085"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eng.backend.get_expectation_value(normalized_hamiltonian, system_qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**7.** As explained in the paper, one can change this state into an eigenstate by undoing the `StatePreparation` of the ancillas and then by measuring if the ancilla qubits in are state 0. The paper says that this should be the case with 50% probability. So let's check this (we require an ancilla to measure this):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5004522593645913"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with Dagger(eng):\n",
    "    StatePreparation(initial_ancilla_wavefunction) | ancillas\n",
    "measure_qb = eng.allocate_qubit()\n",
    "with Compute(eng):\n",
    "    All(X) | ancillas\n",
    "with Control(eng, ancillas):\n",
    "    X | measure_qb\n",
    "Uncompute(eng)\n",
    "eng.flush()\n",
    "eng.backend.get_probability('1', measure_qb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, we would measure 1 (corresponding to the ancilla qubits in state 0) with probability 50% as explained in the paper. Let's assume we measure 1, then we can easily check that we are in an eigenstate of the `normalized_hamiltonian` by numerically calculating its energy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.29263140625433537"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eng.backend.collapse_wavefunction(measure_qb, [1])\n",
    "eng.backend.get_expectation_value(normalized_hamiltonian, system_qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed we are in the ground state of the `normalized_hamiltonian`. Have a look at the paper on how to recover, when the ancilla qubits are not in state 0."
   ]
  },
  {
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}