From 3944d31c05293de42491be2cb22a17a478da41a2 Mon Sep 17 00:00:00 2001 From: kardashin Date: Thu, 21 Feb 2019 15:45:53 +0300 Subject: [PATCH] Add files via upload --- Solving SAT by QAOA.ipynb | 516 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 516 insertions(+) create mode 100644 Solving SAT by QAOA.ipynb diff --git a/Solving SAT by QAOA.ipynb b/Solving SAT by QAOA.ipynb new file mode 100644 index 0000000..88bfba1 --- /dev/null +++ b/Solving SAT by QAOA.ipynb @@ -0,0 +1,516 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Solving the SAT problem using QAOA\n", + "\n", + "In order to embed an instance of the SAT problem into a Hamiltonian $H_c$, we use the following mapping:\n", + "\n", + "$$\n", + " x_j \\rightarrow P^j_0 \\\\ \n", + " \\overline{x}_j \\rightarrow P^j_1 \\\\ \n", + " \\lor \\rightarrow \\otimes \\\\\n", + " \\land \\rightarrow +\n", + "$$\n", + "\n", + "In QAOA, we prepare the state\n", + "$$\n", + " |\\psi(\\boldsymbol{\\alpha}, \\boldsymbol{\\beta})\\rangle = e^{i \\beta_k H_c} e^{i \\alpha_k H_m} \\ldots e^{i \\beta_1 H_c} e^{i \\alpha_1 H_m} |\\boldsymbol{0}\\rangle,\n", + "$$\n", + "\n", + "where *the mixing operator* is chosen to be $H_m = -\\sum_j{\\sigma_x^j},\\, \\boldsymbol{\\alpha} = \\{\\alpha_i\\},\\, \\boldsymbol{\\beta} = \\{\\beta_i\\}$ and $\\alpha_i, \\beta_i \\in [0, 2\\pi)$.\n", + "\n", + "We minimize the expectation value of $H_c$ over the sets of parameters $\\boldsymbol{\\alpha}, \\boldsymbol{\\beta}$, e.g. we calculate:\n", + "\n", + "$$\n", + " E_{min} = \\min_{\\boldsymbol{\\alpha}, \\boldsymbol{\\beta}}{\\langle \\psi(\\boldsymbol{\\alpha}, \\boldsymbol{\\beta})| H_c |\\psi(\\boldsymbol{\\alpha}, \\boldsymbol{\\beta})\\rangle}\n", + "$$\n", + "\n", + "We also search for the sets of parameters\n", + "$$\n", + " \\widetilde{\\boldsymbol{\\alpha}}, \\widetilde{\\boldsymbol{\\beta}} = \\arg\\min_{\\boldsymbol{\\alpha}, \\boldsymbol{\\beta}}{\\langle \\psi(\\boldsymbol{\\alpha}, \\boldsymbol{\\beta})| H_c |\\psi(\\boldsymbol{\\alpha}, \\boldsymbol{\\beta})\\rangle}\n", + "$$\n", + "such that \n", + "$$\n", + " \\langle \\psi(\\widetilde{\\boldsymbol{\\alpha}}, \\widetilde{\\boldsymbol{\\beta}})| H_c |\\psi(\\widetilde{\\boldsymbol{\\alpha}}, \\widetilde{\\boldsymbol{\\beta}})\\rangle = E_{min}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this tutorial, we use the *ProjectQ* package for QAOA calculations.\n", + "\n", + "After we found the smallest eigenvalue, we compare the corresponding eigenvector with the exact solution for a SAT instance." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To install ProjectQ, run\n", + "\n", + "```\n", + "pip install projectq\n", + "```\n", + "\n", + "or, if it fails,\n", + "\n", + "```\n", + "pip install projectq --global-option=--without-cppsimulator\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import projectq\n", + "from projectq.ops import All, Measure, QubitOperator, TimeEvolution\n", + "from projectq.ops import X, Y, Z, H, Rx, Ry, Rz, CNOT, Swap" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from numpy import kron\n", + "import scipy\n", + "from scipy import linalg as la\n", + "from scipy.optimize import minimize\n", + "import random\n", + "from functools import reduce\n", + "import time" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "P_0 = np.array([[1,0],[0,0]]) # |0><0|\n", + "P_1 = np.array([[0,0],[0,1]]) # |1><1|\n", + "X_np = np.array([[0,1],[1,0]]) # X Pauli matrix\n", + "Z_np = np.array([[1,0],[0,-1]]) # Z Pauli matrix\n", + "I_np = np.array([[1,0],[0,1]]) # 2x2 identity matrix " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def cnf_unique_generator(sat, variables_number, clauses_number):\n", + " '''\n", + " Generates a list of unique CNFs\n", + " ''' \n", + " cnf_list = []\n", + " \n", + " while len(cnf_list) <= clauses_number-1:\n", + " a = np.array(sorted(random.sample(range(1, variables_number+1), sat)))\n", + " b = np.array([(-1)**random.randint(0, 1) for i in range(sat)])\n", + " clause = list(a*b)\n", + " \n", + " if clause not in cnf_list:\n", + " cnf_list.append(clause)\n", + " \n", + " return cnf_list" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "An example of a 3SAT instance: \n", + "\n", + " [[1, 2, -4], [-2, -3, 4], [2, -3, -4], [-1, 3, 4], [2, -3, 4]]\n" + ] + } + ], + "source": [ + "print(\"An example of a 3SAT instance: \\n\\n\", cnf_unique_generator(sat=3, variables_number=4, clauses_number=5))" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "def cnf_to_hamiltonian(cnf, sat):\n", + " '''\n", + " Transforms a CNF to a Hamiltonian for ProjectQ\n", + " ''' \n", + " def projector(qubit, index):\n", + " if index == 1:\n", + " return 0.5 * (QubitOperator(\"\") + QubitOperator(\"Z\"+str(qubit)))\n", + " if index == -1: \n", + " return 0.5 * (QubitOperator(\"\") - QubitOperator(\"Z\"+str(qubit)))\n", + " \n", + " sat = len(cnf[0])\n", + "\n", + " hamiltonian = projector(abs(cnf[0][0])-1, np.sign(cnf[0][0]))\n", + " for i in range(1, sat):\n", + " hamiltonian *= projector(abs(cnf[0][i])-1, np.sign(cnf[0][i]))\n", + "\n", + " for clause in cnf[1:]:\n", + " \n", + " ham = projector(abs(clause[0])-1, np.sign(clause[0]))\n", + " for i in range(1, sat):\n", + " ham *= projector(abs(clause[i])-1, np.sign(clause[i]))\n", + "\n", + " hamiltonian += ham\n", + "\n", + " return hamiltonian " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "An example of a 3SAT Hamiltonian:\n", + "\n", + " 0.375 I +\n", + "0.125 Z1 Z2 +\n", + "-0.375 Z0 +\n", + "-0.125 Z0 Z1 Z2 +\n", + "-0.125 Z2 +\n", + "0.125 Z1 +\n", + "0.125 Z0 Z2 +\n", + "-0.125 Z0 Z1\n" + ] + } + ], + "source": [ + "cnf = cnf_unique_generator(sat=3, variables_number=3, clauses_number=3)\n", + "print(\"An example of a 3SAT Hamiltonian:\\n\\n\", cnf_to_hamiltonian(cnf=cnf, sat=3))" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "def cnf_to_matrix(sat_instance, num_bits, pauli=Z_np):\n", + " '''\n", + " Creates a SAT-Hamiltonian as a NumPy array. pauli specifies which basis to use\n", + " ''' \n", + " H = np.zeros((2**num_bits, 2**num_bits))\n", + " for clause in sat_instance:\n", + " term_list = [I_np] * num_bits\n", + " for i in clause:\n", + " term_list[abs(i) - 1] = 0.5 * (I_np + pauli * np.sign(i))\n", + " H += reduce(kron, term_list)\n", + " return H" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def mixing_operator(variables_number):\n", + " '''\n", + " Creates the mixing operator H_m for ProjectQ\n", + " ''' \n", + " operator = QubitOperator('X0')\n", + " for k in range(1, variables_number):\n", + " operator += QubitOperator('X' + str(k))\n", + " \n", + " return -operator" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "An example of a mixing operator: \n", + "\n", + " -1.0 X0 +\n", + "-1.0 X1 +\n", + "-1.0 X2\n" + ] + } + ], + "source": [ + "print(\"An example of a mixing operator: \\n\\n\", mixing_operator(variables_number=3))" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "def mixing_operator_matrix(variables_number):\n", + " '''\n", + " Creates the mixing operator H_m as a NumPy array \n", + " ''' \n", + " sum_operator = np.zeros((2**variables_number,2**variables_number))\n", + " \n", + " for i in range(variables_number):\n", + " \n", + " operator_ini = -X_np\n", + " \n", + " for j in range(i):\n", + " operator_ini = np.kron(I_np, operator_ini)\n", + " \n", + " for k in range(i, variables_number-1):\n", + " operator_ini = np.kron(operator_ini, I_np)\n", + " \n", + " sum_operator += operator_ini\n", + " \n", + " return sum_operator" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "def optimize_sat_instance(sat, hamiltonian, mixer, variables_number, clauses_number, sequence_length, method, engine):\n", + " '''\n", + " Minimizes the expectation value of a given Hamiltonian\n", + " ''' \n", + " def function(x):\n", + " \n", + " # allocate a quantum register in the |00...0> state\n", + " quantum_state = engine.allocate_qureg(variables_number)\n", + " \n", + " # transform to the |++...+> state\n", + " for i in range(variables_number):\n", + " H | quantum_state[i]\n", + " \n", + " # parameters \\alpha and \\beta for QAOA\n", + " parameters = x\n", + " \n", + " # apply the sequence of e^{-i*\\beta*H_m} e^{-i*\\alpha*H_c} operators\n", + " for i in range(sequence_length): \n", + " # apply unitary evolution e^{-i*parameter*operator}\n", + " TimeEvolution(parameters[2*i], hamiltonian) | quantum_state\n", + " TimeEvolution(parameters[2*i+1], mixer) | quantum_state\n", + " \n", + " # send all the gates to the backend\n", + " engine.flush()\n", + " \n", + " # calculate the expectation value\n", + " expectated_value = engine.backend.get_expectation_value(hamiltonian, quantum_state)\n", + " \n", + " # deallocate the quantum register\n", + " All(Measure) | quantum_state\n", + "\n", + " return expectated_value\n", + " \n", + " initial_parameters = [0.25*np.pi for i in range(2*sequence_length)]\n", + " \n", + " bounds = [(0, 2*np.pi) for i in range(2*sequence_length)]\n", + " \n", + " result = minimize(function, initial_parameters, method=method, bounds=bounds)\n", + " \n", + " return result" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This function takes the sets $\\boldsymbol{\\alpha}, \\boldsymbol{\\beta}$ and constructs $|\\psi(\\boldsymbol{\\alpha}, \\boldsymbol{\\beta})\\rangle = e^{i \\beta_k H_c} e^{i \\alpha_k H_m} \\ldots e^{i \\beta_1 H_c} e^{i \\alpha_1 H_m} |\\boldsymbol{0}\\rangle$ as a NumPy array:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "def construct_vector(variables_number, mixing_operator_matrix, hamiltonian_matrix, sequence_length, parameters):\n", + " '''\n", + " Creates the approximate ground eigensate as a NumPy vector\n", + " ''' \n", + " # |+> state\n", + " vector = np.array([1]*(2**variables_number)) / np.sqrt(2**variables_number)\n", + " \n", + " # apply the sequence of operators\n", + " for i in range(sequence_length):\n", + " vector = la.expm(-1j*parameters[2*i+1]*mixing_operator_matrix).dot(la.expm(-1j*parameters[2*i]*hamiltonian_matrix)).dot(vector)\n", + " \n", + " return vector" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Suppose that the SAT instance is satisfiable, then the solution vector lies in the span of all ground eigenvectors.\n", + "\n", + "The *check_overlap(...)* function finds the set of ground eigenvectors of a given Hamiltonian, $\\{|\\varphi_1\\rangle, \\ldots, |\\varphi_l\\rangle \\}$, and then calculates the overlap between the approximate eigenvector $|\\psi(\\widetilde{\\boldsymbol{\\alpha}}, \\boldsymbol{\\widetilde{\\beta}})\\rangle$ and this set as\n", + "\n", + "$$\n", + " overlap = \\sum_j{|\\langle\\varphi_j|\\psi(\\widetilde{\\boldsymbol{\\alpha}}, \\boldsymbol{\\widetilde{\\beta}})\\rangle|^2}\n", + "$$\n", + "\n", + "This quantity shows how far our solution is from the exact one." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "def check_overlap(target_vector, hamiltonian_matrix, variables_number):\n", + " '''\n", + " Calculates the overlap between the approximate eigenstate and the span of exact eigenstates\n", + " ''' \n", + " eigenvalues, eigenvectors = la.eig(hamiltonian_matrix)\n", + " ground_eigenvectors_indices = np.where(eigenvalues == min(eigenvalues))\n", + " \n", + " overlap = 0\n", + " for l in ground_eigenvectors_indices[0]:\n", + " overlap += abs(target_vector.conjugate().transpose().dot(eigenvectors[l]/np.linalg.norm(eigenvectors[l])))**2 \n", + " \n", + " return overlap" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Running an Example " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If the approximate eigenvalue is smaller than $1$, then we consider that the SAT instance is satisfiable." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(Note: This is the (slow) Python simulator.)\n" + ] + } + ], + "source": [ + "variables_number = 3\n", + "clauses_number = 7\n", + "sequence_length = 5\n", + "sat = 3\n", + "\n", + "# optimization method\n", + "method = 'L-BFGS-B'\n", + "\n", + "# set the default backend for ProjectQ\n", + "engine = projectq.MainEngine()\n", + "\n", + "# generate a SAT instance\n", + "cnf = cnf_unique_generator(sat, variables_number, clauses_number)\n", + "\n", + "# generate a Hamiltonian and a mixing operator for ProjectQ\n", + "hamiltonian = cnf_to_hamiltonian(cnf, sat)\n", + "mixer = mixing_operator(variables_number)\n", + "\n", + "# generate a Hamiltonian and a mixing operator as the NumPy arrays\n", + "hamiltonian_matrix = cnf_to_matrix(cnf, variables_number)\n", + "mixer_matrix = mixing_operator_matrix(variables_number)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time: 29.50817584991455\n", + "Approximate eigenvalue: 1.3029854972756993e-12\n", + "Exact eigenvalue: 0.0\n", + "Overlap between approxiate and exact eigenvector: 0.999999999998697\n" + ] + } + ], + "source": [ + "# optimize the Hamiltonian, obtain the ground eigenvalue and the sets of parameters \\alpha and \\beta\n", + "start_time = time.time()\n", + "result = optimize_sat_instance(sat, hamiltonian, mixer, variables_number, clauses_number, sequence_length, method, engine)\n", + "finish_time = time.time() - start_time\n", + "\n", + "# construct the ground eigenvector using the \\alpha and \\beta parameters\n", + "parameters = result.x\n", + "vec_opt = construct_vector(variables_number, mixer_matrix, hamiltonian_matrix, sequence_length, parameters)\n", + "\n", + "print(\"Time:\", finish_time)\n", + "print(\"Approximate eigenvalue:\", result.fun)\n", + "print(\"Exact eigenvalue:\", la.eigh(hamiltonian_matrix)[0][0])\n", + "print(\"Overlap between approxiate and exact eigenvector:\", check_overlap(vec_opt, hamiltonian_matrix, variables_number))\n", + "#print(\"Parameters:\", parameters)" + ] + }, + { + "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.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}