08_Sampling_a_Thermal_State.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T22:01:50.475321Z",
     "start_time": "2019-01-29T22:01:49.695728Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Available frameworks:\n",
      "Qiskit\n",
      "D-Wave Ocean\n"
     ]
    }
   ],
   "source": [
    "%run -i \"assignment_helper.py\"\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Quantum annealing\n",
    "\n",
    "**Exercise 1** (2 points). The thermal state: $$\\rho=\\frac{1}{Z} \\sum_n e^{-E_n/T} |n\\rangle \\langle n|$$ is a Boltzmann distribution over the energy eigenvalues *En* of some Hamiltonian *H*. An open quantum system equilibrates with the environment to this state. Equilibration is hard to simulate classically and therefore this process could be exploited for calculations. Using dimod, create a random Ising model over 5 spins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T21:52:31.512988Z",
     "start_time": "2019-01-29T21:52:31.507436Z"
    }
   },
   "outputs": [],
   "source": [
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "\n",
    "n_spins = 5\n",
    "\n",
    "h = {i: 1. for i in range(n_spins)}\n",
    "J = {(i, i+1): 0. for i in range(n_spins-1)}\n",
    "\n",
    "model = dimod.BinaryQuadraticModel(h, J, 0., dimod.SPIN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T21:52:32.433092Z",
     "start_time": "2019-01-29T21:52:32.429882Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise1",
     "locked": true,
     "points": "2",
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "assert isinstance(model, dimod.binary_quadratic_model.BinaryQuadraticModel), \"Wrong model type\"\n",
    "assert model.vartype == dimod.SPIN, \"Wrong variables: binary model instead of spin system\"\n",
    "assert len(model.variables) == 5, \"Wrong number of spins\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2** (2 points). Sample the model at near-zero temperature (0.01) and at a high temperature (100). Write the list of energies in two arrays, `energies_low` and `energies_high`. Note that beta denotes the inverse temperature *1/T*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T21:52:46.422181Z",
     "start_time": "2019-01-29T21:52:45.192862Z"
    }
   },
   "outputs": [],
   "source": [
    "n_samples = 100\n",
    "sampler = dimod.SimulatedAnnealingSampler()\n",
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "T_low  = 0.01\n",
    "T_high = 100.\n",
    "\n",
    "sampleset_low  = sampler.sample(model, beta_range=[1./T_low, 1./T_low], num_reads=n_samples)\n",
    "sampleset_high = sampler.sample(model, beta_range=[1./T_high, 1./T_high], num_reads=n_samples)\n",
    "\n",
    "energies_low  = [sampleset_low.record[i][1]  for i in range(len(sampleset_low.record))]\n",
    "energies_high = [sampleset_high.record[i][1] for i in range(len(sampleset_high.record))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-5.0]\n",
      "[-5.0, -3.0, -1.0, 1.0, 3.0, 5.0]\n"
     ]
    }
   ],
   "source": [
    "print(sorted(list(dict.fromkeys(energies_low))))\n",
    "print(sorted(list(dict.fromkeys(energies_high))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T21:52:46.455503Z",
     "start_time": "2019-01-29T21:52:46.452117Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise2",
     "locked": true,
     "points": "2",
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "assert np.isclose(min(energies_low), max(energies_low)), \"At near-zero temperature, we should only see one energy value, the ground state energy\"\n",
    "assert np.isclose(min(energies_low), min(energies_high)), \"The ground state should appear at both temperatures\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as we saw in the introduction, the distribution flattens out at a high temperature (*T3*). On the other hand, the energy is peaked for a low temperature, and we do not even have samples for high-energy states."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Approximate Thermalization\n",
    "\n",
    "**Exercise 3** (3 points). In QAOA, we approximated the ground state of a target Hamiltonian, starting from the ground state of a mixer Hamiltonian. We can actually perform the exact same optimization for approximating the thermal state of a target system, starting from the thermal state of some system.\n",
    "\n",
    "We exploit that if we trace out a subsystem of an entangled system, we end up with a mixed state.\n",
    "\n",
    "It can be shown that: $$| \\psi \\rangle =1/\\sqrt{2 \\cosh \\frac{1}{2T}} \\sum_{z \\in {-1,1}} e^{- z/T} |z \\rangle_{\\mathcal{H_1}} \\otimes | z \\rangle_{\\mathcal{H_2}}$$ purifies: $$\\rho=\\frac{1}{Z}e^{- H_m/T}$$ where *Hm* is the mixing Hamiltonian in QAOA. This state can be built with a circuit composed uniquely of RX gates and CNOT gates. Since the overall structure of the algorithm resembles QAOA closely, it is less important to go through the details. Most of it is more about forcing the quantum computing packages to do what we want than about the problem itself. The biggest difference to QAOA is the preparation of the initial state. Prepare this state at inverse temperature *beta=0* in your preferred framework. Note that *beta=1/T*, so this corresponds to the infinite temperature. The state preparation should be in an object called `circuit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T22:01:59.135749Z",
     "start_time": "2019-01-29T22:01:59.125142Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">      ┌───┐┌───┐┌───┐     \n",
       "q0_0: ┤ H ├┤ Z ├┤ S ├──■──\n",
       "      └───┘└───┘└───┘┌─┴─┐\n",
       "q0_1: ───────────────┤ X ├\n",
       "                     └───┘</pre>"
      ],
      "text/plain": [
       "      ┌───┐┌───┐┌───┐     \n",
       "q0_0: ┤ H ├┤ Z ├┤ S ├──■──\n",
       "      └───┘└───┘└───┘┌─┴─┐\n",
       "q0_1: ───────────────┤ X ├\n",
       "                     └───┘"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_qubits = 1\n",
    "β = 0\n",
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "\n",
    "circuit = QuantumCircuit()\n",
    "qr = QuantumRegister(2)\n",
    "circuit.add_register(qr)\n",
    "\n",
    "\n",
    "circuit.h(qr[0])\n",
    "circuit.z(qr[0])\n",
    "\n",
    "circuit.s(qr[0])\n",
    "circuit.cx(qr[0],qr[1])\n",
    "\n",
    "circuit.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T22:02:01.252611Z",
     "start_time": "2019-01-29T22:02:00.773007Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise3",
     "locked": true,
     "points": "3",
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "amplitudes = get_amplitudes(circuit)\n",
    "assert np.allclose(amplitudes, np.array([1/np.sqrt(2), 0, 0, -1j/np.sqrt(2)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4** (1 points).  Create the density matrix for this quantum state in a variable called `dm`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T22:02:12.122651Z",
     "start_time": "2019-01-29T22:02:12.104843Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.000000e+00+0.j , 3.061617e-17-0.5j],\n",
       "       [0.000000e+00+0.j , 0.000000e+00+0.j ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "\n",
    "amplitudes = amplitudes.reshape(2,2)\n",
    "\n",
    "dm = np.outer(amplitudes[0], amplitudes[1])\n",
    "dm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T22:02:12.656345Z",
     "start_time": "2019-01-29T22:02:12.650757Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise4",
     "locked": true,
     "points": "1",
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "###\n",
    "### AUTOGRADER TEST - DO NOT REMOVE\n",
    "###\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you look at the structure of this density matrix and trace out the ancilla system, you get: $$\\begin{bmatrix}0.5&0\\\\0&0.5\\end{bmatrix}$$ That is the maximally mixed state, which corresponds to the thermal state at infinite temperature, which is exactly what we wanted."
   ]
  },
  {
   "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.8.8"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}