07_Variational_Circuits.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:28:53.754490Z",
     "start_time": "2019-01-30T17:28:52.741902Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/gianlucafugante/anaconda3/lib/python3.8/site-packages/qiskit/aqua/__init__.py:86: DeprecationWarning: The package qiskit.aqua is deprecated. It was moved/refactored to qiskit-terra For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide>\n",
      "  warn_package('aqua', 'qiskit-terra')\n"
     ]
    },
    {
     "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": {},
   "source": [
    "# Quantum approximate optimization algorithm\n",
    "\n",
    "QAOA is a shallow-circuit variational algorithm that is easy to understand if you already grasped quantum annealing. It is, in fact, just a particular type of a larger family of algorithms called variational quantum eigensolvers. Approximating the adiabatic pathway as QAOA does is just one option of how to find the eigenvalues of a system.\n",
    "\n",
    "Even then, QAOA has many moving elements. Let us import some handy packages and define some functions that we are going to use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:28:54.708430Z",
     "start_time": "2019-01-30T17:28:54.700061Z"
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "from functools import partial, reduce\n",
    "from scipy.optimize import minimize\n",
    "np.set_printoptions(precision=3, suppress=True)\n",
    "\n",
    "# Functions useful if you're using Qiskit\n",
    "def pauli_x(qubit, coeff):\n",
    "    eye = np.eye((n_qubits))\n",
    "    return PauliOp(Pauli((np.zeros(n_qubits), eye[qubit])), coeff)\n",
    "\n",
    "def pauli_z(qubit, coeff):\n",
    "    eye = np.eye((n_qubits))\n",
    "    return PauliOp(Pauli((eye[qubit], np.zeros(n_qubits))), coeff)\n",
    "\n",
    "def product_pauli_z(q1, q2, coeff):\n",
    "    eye = np.eye((n_qubits))\n",
    "    return PauliOp(Pauli((eye[q1], np.zeros(n_qubits))) & Pauli((eye[q2], np.zeros(n_qubits))), coeff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1** (2 points). Define a mixing Hamiltonian on two qubits. Store it in an object called `Hm`. If you're doing the assignement in PyQuil, remember that the Hamiltonian should be a list of PauliTerms (and not a PauliSum) in order to be exponentiated easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:28:56.353161Z",
     "start_time": "2019-01-30T17:28:56.349370Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/gianlucafugante/anaconda3/lib/python3.8/site-packages/qiskit/aqua/operators/operator_base.py:46: DeprecationWarning: The package qiskit.aqua.operators is deprecated. It was moved/refactored to qiskit.opflow (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide>\n",
      "  warn_package('aqua.operators', 'qiskit.opflow', 'qiskit-terra')\n"
     ]
    }
   ],
   "source": [
    "n_qubits = 2\n",
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "\n",
    "\n",
    "Hm = pauli_x(0, 1.0) + pauli_x(1, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:28:57.312028Z",
     "start_time": "2019-01-30T17:28:57.305208Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise1",
     "locked": true,
     "points": "2",
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "if isinstance(Hm, SummedOp):\n",
    "    assert np.alltrue(Hm.to_matrix() == np.array([[0., 1., 1., 0.],\n",
    "                                                  [1., 0., 0., 1.],\n",
    "                                                  [1., 0., 0., 1.],\n",
    "                                                  [0., 1., 1., 0.]]))\n",
    "elif isinstance(Hm, list):\n",
    "    assert len(Hm) == n_qubits\n",
    "    assert all([isinstance(Hm[i], PauliTerm) for i in range(n_qubits)])\n",
    "    assert all([Hm[i].compact_str() == '(1+0j)*X{}'.format(i) for i in range(n_qubits)])\n",
    "else:\n",
    "    raise ValueError(\"Unknown type for Hamiltonian!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2** (2 points). Define the cost Hamiltonian: $$H_c = -\\sigma^Z_1\\sigma^Z_2-0.5\\sigma^Z_1$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:07.809767Z",
     "start_time": "2019-01-30T17:29:07.805940Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.5+0.j, -0. +0.j, -0. +0.j, -0. +0.j],\n",
       "       [-0. +0.j,  1.5+0.j, -0. +0.j, -0. +0.j],\n",
       "       [-0. +0.j, -0. +0.j,  0.5+0.j, -0. +0.j],\n",
       "       [-0. +0.j, -0. +0.j, -0. +0.j, -0.5+0.j]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "\n",
    "Hc = product_pauli_z(0, 1, -1.0) + pauli_z(0, -0.5)\n",
    "Hc.to_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:08.286947Z",
     "start_time": "2019-01-30T17:29:08.280462Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise2",
     "locked": true,
     "points": "2",
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "if isinstance(Hc, SummedOp):\n",
    "    assert np.alltrue(Hc.to_matrix() == np.diag([-1.5, 1.5, 0.5, -0.5]))\n",
    "    \n",
    "elif isinstance(Hc, list):\n",
    "    assert len(Hc) == 2\n",
    "    assert all([isinstance(Hc[i], PauliTerm) for i in range(n_qubits)])\n",
    "    assert Hc[0].compact_str() == '(-1+0j)*Z0Z1'\n",
    "    assert Hc[1].compact_str() == '(-0.5+0j)*Z0'\n",
    "else:\n",
    "    raise ValueError(\"Unknown type for Hamiltonian!\")    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3** (2 points). We know that the ground state of the mixer Hamiltonian is the uniform superposition. Create a circuit `circuit_init` that will contain this initial state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:20.701673Z",
     "start_time": "2019-01-30T17:29:20.698795Z"
    }
   },
   "outputs": [],
   "source": [
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "\n",
    "init_state = StateFn([0.5, 0.5, 0.5, 0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:22.069567Z",
     "start_time": "2019-01-30T17:29:21.927383Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise3",
     "locked": true,
     "points": "2",
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "amplitudes = init_state.to_matrix()\n",
    "assert np.allclose(amplitudes, np.array([0.5, 0.5, 0.5, 0.5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set p=2 and initialize the beta_i and gamma_i parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:33.696855Z",
     "start_time": "2019-01-30T17:29:33.693611Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4.82  4.952] [4.388 1.354]\n"
     ]
    }
   ],
   "source": [
    "p = 2\n",
    "beta = np.random.uniform(0, np.pi*2, p)\n",
    "gamma = np.random.uniform(0, np.pi*2, p)\n",
    "\n",
    "print(beta, gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to create the complete variational circuit, made of $$e^{-\\beta H} \\space and \\space e^{-\\gamma H}$$ We will use a function `create_circuit` that takes `gamma` and `beta` as argument, and the state preparation circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:34.599407Z",
     "start_time": "2019-01-30T17:29:34.587974Z"
    }
   },
   "outputs": [],
   "source": [
    "def evolve(hamiltonian, angle):\n",
    "    return (angle * hamiltonian).exp_i()\n",
    "\n",
    "def create_circuit(init_state, beta, gamma):\n",
    "    if isinstance(init_state, qiskit.aqua.operators.StateFn):\n",
    "        circuit = ComposedOp([evolve(Hm, beta[i]) @ evolve(Hc, gamma[i]) for i in range(p)] + [init_state])\n",
    "    elif isinstance(circuit_init, pyquil.quil.Program):\n",
    "        exp_Hm = []\n",
    "        exp_Hc = []\n",
    "        for term in Hm:\n",
    "            exp_Hm.append(exponential_map(term))\n",
    "        for term in Hc:\n",
    "            exp_Hc.append(exponential_map(term))\n",
    "        circuit = Program()\n",
    "        circuit += circuit_init\n",
    "        for i in range(p):\n",
    "            for term_exp_Hm in exp_Hm:\n",
    "                circuit += term_exp_Hm(-beta[i])\n",
    "            for term_exp_Hc in exp_Hc:\n",
    "                circuit += term_exp_Hc(-gamma[i])\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we need a function `evaluate_circuit` to compute the average energy of the circuit, i.e. compute: $$\\langle\\psi(\\beta, \\gamma)|H_c|\\psi(\\beta, \\gamma)\\rangle$$ where: $$|\\psi(\\beta, \\gamma)\\rangle$$ is the circuit built above. This function should take a unique argument `beta_gamma` (concatenation of the lists `beta` and `gamma`) in order to be used directly by optimizers, and return a real value corresponding to the expectancy of *Hc*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:35.802611Z",
     "start_time": "2019-01-30T17:29:35.797170Z"
    }
   },
   "outputs": [],
   "source": [
    "def evaluate_circuit(beta_gamma):\n",
    "    n = len(beta_gamma)//2\n",
    "    circuit = create_circuit(init_state, beta_gamma[:n], beta_gamma[n:])\n",
    "    if isinstance(circuit, qiskit.aqua.operators.OperatorBase):\n",
    "        meas = OperatorStateFn(Hc, is_measurement=True)\n",
    "        circuit_meas = PauliTrotterEvolution(trotter_mode=Suzuki(order=3, reps=1)).convert(meas @ circuit)\n",
    "        circuit_expectation = PauliExpectation().convert(circuit_meas)\n",
    "        \n",
    "        sampler = CircuitSampler(backend=Aer.get_backend('statevector_simulator'))\n",
    "        sampled_trotter_exp_op = sampler.convert(circuit_expectation)\n",
    "        return sampled_trotter_exp_op.eval().real\n",
    "    elif isinstance(circuit, pyquil.quil.Program):\n",
    "        qvm = pyquil.api.QVMConnection(endpoint=fc.sync_endpoint, compiler_endpoint=fc.compiler_endpoint)\n",
    "        return np.real(qvm.pauli_expectation(circuit, sum(Hc)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4** (2 points). The *p* parameter defines the number of steps in the Trotterization. The real question from here is how we optimize the beta_i and gamma_i parameters. If we can find a method that makes fewer evaluations to arrive at the same result, that is a win, since we have to execute fewer loops on the quantum computer. Try various methods for minimizing the evaluate function. We used L-BFGS-B before. Try another one and write the outcome in an object called `result`. You will see that the number of function evaluation (`nfev`) differs and so does the function value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:29:39.475428Z",
     "start_time": "2019-01-30T17:29:38.504544Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/gianlucafugante/anaconda3/lib/python3.8/site-packages/qiskit/aqua/quantum_instance.py:135: DeprecationWarning: The class qiskit.aqua.QuantumInstance is deprecated. It was moved/refactored to qiskit.utils.QuantumInstance (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide>\n",
      "  warn_class('aqua.QuantumInstance',\n"
     ]
    }
   ],
   "source": [
    "###\n",
    "### YOUR CODE HERE\n",
    "###\n",
    "result = minimize(evaluate_circuit, np.concatenate([beta, gamma]), method='Nelder-Mead')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-30T17:30:02.493851Z",
     "start_time": "2019-01-30T17:30:02.199109Z"
    },
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "grade": true,
     "grade_id": "exercise4",
     "locked": true,
     "points": "2",
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Function evaluations: 100 versus 207\n",
      "Function values: -1.191382 versus -1.050576\n"
     ]
    }
   ],
   "source": [
    "result2 = minimize(evaluate_circuit, np.concatenate([beta, gamma]), method='L-BFGS-B')\n",
    "import scipy\n",
    "assert isinstance(result, scipy.optimize.optimize.OptimizeResult)\n",
    "assert result2.nfev != result.nfev\n",
    "print(\"Function evaluations: %d versus %d\" % (result2.nfev, result.nfev))\n",
    "print(\"Function values: %f versus %f\" % (result2.fun, result.fun))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the circuit had an infinite capacity, you should obtain a minimum of -1.5 (minimum value of the Hamiltonian we defined above). However, with *p=2*, you might have a bigger value. Try increasing *p* to see the effect on the minimum reached by the circuit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variational circuits are actually very similar to neural networks: we do some kind of gradient descent over a parameter space. Automatic differentiation of neural networks has been a major step in scaling them up and we can expect that similar techniques in variational quantum circuit can be useful. This is exactly what projects like [QuantumFlow](https://github.com/rigetti/quantumflow) and [PennyLane](https://pennylane.ai/) try."
   ]
  }
 ],
 "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
}