prepare_gaussian_state.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Hubbard Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.utils import HubbardSquareLattice \n",
    "\n",
    "x_n = 2 \n",
    "y_n = 2 \n",
    "n_dofs = 1\n",
    "periodic = 0\n",
    "spinless = 0 \n",
    "\n",
    "lattice = HubbardSquareLattice(x_n, y_n, n_dofs=n_dofs, periodic=periodic, spinless=spinless)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.hamiltonians import FermiHubbardModel\n",
    "from openfermion.utils import SpinPairs\n",
    "tunneling = [('neighbor', (0, 0), 1.)] \n",
    "interaction = [('onsite', (0, 0), 2., SpinPairs.DIFF)] \n",
    "potential = None\n",
    "mag_field = 0. \n",
    "particle_hole_sym = False \n",
    "\n",
    "hubbard = FermiHubbardModel(lattice , tunneling_parameters=tunneling, interaction_parameters=interaction, \n",
    "                            potential_parameters=potential, magnetic_field=mag_field, \n",
    "                            particle_hole_symmetry=particle_hole_sym)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find tunneling eigenvectors with highest fidelity with Hubbard ground state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvector v_tun[:, 11] has ovelap 0.12622062025255504.\n",
      "Eigenvector v_tun[:, 19] has ovelap 0.1985620731386243.\n",
      "Eigenvector v_tun[:, 34] has ovelap 0.10964243588147657.\n"
     ]
    }
   ],
   "source": [
    "from openfermion.utils import inner_product \n",
    "from openfermion.transforms import get_sparse_operator\n",
    "\n",
    "def overlap(a, b):\n",
    "    \"\"\"Calculates the overlap between vectors a and b. This metric is also known as fidelity. \"\"\"\n",
    "    inner = inner_product(a, b)\n",
    "    return (np.conjugate(inner) * inner).real\n",
    "\n",
    "# Lowest eigenvalue/vector of Hubbard Hamiltonian\n",
    "# The column v_hub[:, i] is the eigenvector corresponding to the eigenvalue w_hub[i]\n",
    "hub_sparse = get_sparse_operator(hubbard.hamiltonian())\n",
    "w_hub, v_hub = scipy.sparse.linalg.eigsh(hub_sparse, k=1, which='SA') \n",
    "\n",
    "# From looking at the spectrum,  I noticed the ground state was degenerate. \n",
    "# There were 16 eigenvectors with the lowest eigenvalue. \n",
    "# Actually, the initial state might not even be a ground state. Let's try \n",
    "# as many of the eigenvectors as we can. \n",
    "tun_sparse = get_sparse_operator(hubbard.tunneling_terms())\n",
    "# k can be at most n-2 where n is the dimension of the matrix\n",
    "w_tun, v_tun = scipy.sparse.linalg.eigsh(tun_sparse, k=np.shape(tun_sparse)[0]-2, which='SA')\n",
    "\n",
    "# Let's check overlap of each of the ground states. Hopefully, there's an obvious winner. \n",
    "for i in range(len(w_tun)):\n",
    "    fid = overlap(v_hub[:, 0], v_tun[:, i]) \n",
    "    if fid > 0.10:\n",
    "        print('Eigenvector v_tun[:, {}] has ovelap {}.'.format(i, fid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "None of these stand out as the obvious best choice, so we'll perturb tunneling term with the interaction term and hope we have better luck. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find perturbed eigenvectors with highest fidelity with Hubbard ground state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvector v_per[:, 0] has overlap 0.9770315592604026.\n",
      "Eigenvector v_per[:, 1] has overlap 3.658337908110008e-20.\n",
      "Eigenvector v_per[:, 2] has overlap 5.060794393010509e-25.\n",
      "Eigenvector v_per[:, 3] has overlap 8.521614327774115e-26.\n",
      "Eigenvector v_per[:, 4] has overlap 2.4103963324261132e-25.\n"
     ]
    }
   ],
   "source": [
    "s = 1e-4\n",
    "int_sparse = get_sparse_operator(hubbard.interaction_terms())\n",
    "perturbed_sparse = tun_sparse + s * int_sparse\n",
    "\n",
    "w_per,v_per = scipy.sparse.linalg.eigsh(perturbed_sparse, k=5, which='SA') \n",
    "\n",
    "for i in range(len(w_per)):\n",
    "    fid = overlap(v_hub[:, 0], v_per[:, i])\n",
    "    print('Eigenvector v_per[:, {}] has overlap {}.'.format(i, fid))\n",
    "    \n",
    "# WOW, there's a clear winner here v_per[:, 0]\n",
    "per_state_most_overlap = v_per[:, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find tunneling eigenvector with highest fidelity with best perturbed eigenvector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvector v_tun[:, 19] had the maximum overlap of 0.2032311164794801 with the best perturbed state.        It has eigenvalue -4.000000000000009\n"
     ]
    }
   ],
   "source": [
    "max_overlap = 0\n",
    "index_max_overlap = 0\n",
    "\n",
    "for i in range(len(w_tun)):\n",
    "    fid = overlap(per_state_most_overlap, v_tun[:, i])\n",
    "    if fid > max_overlap: \n",
    "        max_overlap = fid \n",
    "        index_max_overlap = i \n",
    "print(\"Eigenvector v_tun[:, {}] had the maximum overlap of {} with the best perturbed state. \\\n",
    "       It has eigenvalue {}\".format(index_max_overlap, max_overlap, w_tun[index_max_overlap]))\n",
    "initial_state = v_tun[:, index_max_overlap]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use `prepare_gaussian_state` and `final_wavefunction` to find circuit that produces state with highest fidelity with above tunneling eigenvector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2.00000000e+00 -2.00000000e+00 -1.91765100e-16  1.64346022e-32\n",
      "  8.93943564e-17  5.82889146e-16  2.00000000e+00  2.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "from openfermion.transforms import get_quadratic_hamiltonian \n",
    "\n",
    "# Convert to QuadraticHamiltonian instance \n",
    "tun_quad = get_quadratic_hamiltonian(hubbard.tunneling_terms())\n",
    "\n",
    "orbital_energies, constant = tun_quad.orbital_energies()\n",
    "print(orbital_energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Orbital energies [0, 1] resulted in a state with 0.2032310999224089 overlap with our desired state.\n"
     ]
    }
   ],
   "source": [
    "# Get all possible combinations of orbital energies which have a sum of -4 \n",
    "import itertools \n",
    "zero_indices = [2, 3, 4, 5] # Indices of 0's in orbital_energies\n",
    "orbital_energies_combs = [] \n",
    "for r in range(len(zero_indices) + 1):\n",
    "    for subset in itertools.combinations(zero_indices, r):\n",
    "        # Include indices for both -2's and then indices of some 0's\n",
    "        orbital_energies_combs.append([0, 1] + list(subset))\n",
    "        \n",
    "# Now, generate each eigenvector from a combination of orbital_energies \n",
    "# and compare it to our desired tunneling eigenvector. \n",
    "# This way we know which orbital_energies combination yields our \n",
    "# desired eigenvector. \n",
    "\n",
    "# from openfermion.utils import gaussian_state_preparation_circuit\n",
    "from openfermioncirq import prepare_gaussian_state\n",
    "from cirq import Circuit, final_wavefunction, LineQubit\n",
    "\n",
    "overlaps = []\n",
    "for comb in orbital_energies_combs:\n",
    "    state = final_wavefunction(Circuit(prepare_gaussian_state(\n",
    "            LineQubit.range(8), # We have an 8 qubit Hamiltonian \n",
    "            tun_quad, \n",
    "            occupied_orbitals=comb)))\n",
    "    overlaps.append(overlap(state, initial_state))\n",
    "\n",
    "# Find state with top overlap and save the Circuit object that created it \n",
    "best_state_index = np.argmax(overlaps)\n",
    "\n",
    "print(\"Orbital energies {} resulted in a state with {} overlap with our desired state.\".format(\n",
    "    orbital_energies_combs[best_state_index], overlaps[best_state_index]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Shouldn't we have a much higher overlap here? We use `prepare_gaussian_state` to generate ground states of the tunneling term and calculate fidelity with another ground state of the tunneling term?** We should hit an eigenvector with fidelity close to 1, right?\n",
    "\n",
    "If I instead compute fidelity with the best *perturbed* state, I get an eigenvector with fidelity close to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Orbital energies [0, 1] resulted in a state with 0.999999917691482 overlap with our desired state.\n"
     ]
    }
   ],
   "source": [
    "overlaps = []\n",
    "for comb in orbital_energies_combs:\n",
    "    state = final_wavefunction(Circuit(prepare_gaussian_state(\n",
    "            LineQubit.range(8), # We have an 8 qubit Hamiltonian \n",
    "            tun_quad, \n",
    "            occupied_orbitals=comb)))\n",
    "    overlaps.append(overlap(state, per_state_most_overlap))\n",
    "\n",
    "best_state_index = np.argmax(overlaps)\n",
    "print(\"Orbital energies {} resulted in a state with {} overlap with our desired state.\".format(\n",
    "    orbital_energies_combs[best_state_index], overlaps[best_state_index]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is good, but it seems we just got lucky with the basis that `occupied_orbitals` chose. A degenerate ground state means any vector in the span of the ground states is a valid ground state. Is there any way we can specify the basis that `occupied_orbitals` chooses? Or at least specify linear combinations of the basis vectors? "
   ]
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}