Solving SAT by QAOA.ipynb

{
 "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",
    "pip3 install projectq\n",
    "```\n",
    "\n",
    "or, if it fails,\n",
    "\n",
    "```\n",
    "pip3 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, 3, 4], [1, 2, 3], [-1, 2, 4], [-2, -3, 4], [-1, -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": 28,
   "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": 16,
   "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.125 Z0 +\n",
      "0.375 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": 17,
   "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": 18,
   "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": 19,
   "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": 20,
   "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": 21,
   "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": 22,
   "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": 23,
   "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": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables_number = 4\n",
    "clauses_number = 8\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": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time: 163.4611372947693\n",
      "Approximate eigenvalue: 0.008393808893097142\n",
      "Exact eigenvalue: 0.0\n",
      "Overlap between approxiate and exact eigenvector: 0.9920214114585775\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": "myenv",
   "language": "python",
   "name": "myenv"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}