Solving SAT by VQE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational quantum eigensolver\n",
    "\n",
    "Any qubit Hamiltonian can be decomposed into a sum of Pauli operators:\n",
    "$$\n",
    "H = \\sum_{\\boldsymbol{\\alpha} \\in \\{0 ... 3\\}^n} J_{\\boldsymbol{\\alpha}} \\   \n",
    "\\sigma_{\\alpha_1} \\otimes \\dots \\otimes \\sigma_{\\alpha_n}, \n",
    "\\quad J_{\\boldsymbol{\\alpha}} \\in \\mathbb{R}\n",
    "$$\n",
    "\n",
    "If the number of terms scales polynomially with the number of qubits, then \n",
    "$\\langle \\psi(\\boldsymbol{\\theta}) \\ H \\ \\psi(\\boldsymbol{\\theta}) \\rangle$ can be efficiently measured on a quantum computer: prepare the state $\\psi(\\boldsymbol{\\theta})$ and measure each term separately. Then minimize the total energy over $\\boldsymbol{\\theta}$ [Peruzzo, A. et al.  Nat. Comm. 2014, 5 (1)]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alexey/quantumenv/lib/python3.6/site-packages/marshmallow/schema.py:364: ChangedInMarshmallow3Warning: strict=False is not recommended. In marshmallow 3.0, schemas will always be strict. See https://marshmallow.readthedocs.io/en/latest/upgrading.html#schemas-are-always-strict\n",
      "  ChangedInMarshmallow3Warning\n"
     ]
    }
   ],
   "source": [
    "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister\n",
    "from qiskit import Aer, execute\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from skopt import gp_minimize\n",
    "from scipy.optimize import minimize\n",
    "from itertools import cycle\n",
    "from functools import reduce\n",
    "\n",
    "import warnings\n",
    "import sys\n",
    "\n",
    "if not sys.warnoptions:\n",
    "    warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tensor network states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TensorNetwork:\n",
    "    \"\"\"\n",
    "    Tensor Network states prepared on qubits\n",
    "    \"\"\"\n",
    "    backend = \"qasm_simulator\"\n",
    "    \n",
    "    def __init__(self, q, c):\n",
    "        self.backend = \"qasm_simulator\"\n",
    "        self.q = q\n",
    "        self.c = c\n",
    "        self.n_qubits = len(q)\n",
    "    \n",
    "    def set_params(self, params):\n",
    "        if self.n_params == len(params):\n",
    "            self.params = params\n",
    "        else:\n",
    "            raise ValueError('Incorrect number of parameters!')\n",
    "            \n",
    "class TNWithEntangler(TensorNetwork):\n",
    "    '''Tensor network that requires an entangler'''\n",
    "    def __init__(self, q, c, entangler):\n",
    "        super().__init__(q, c)\n",
    "        if not isinstance(entangler, Entangler):\n",
    "            raise TypeError('Should supply an entangler')\n",
    "        self.entangler = entangler            \n",
    "            \n",
    "class RankOne(TensorNetwork):\n",
    "    '''Tensor network for building rank-1 (unentangled) states'''\n",
    "    def __init__(self, q, c):\n",
    "        super().__init__(q, c)\n",
    "        self.n_params = len(q) * 2\n",
    "        self.params = [0] * self.n_params\n",
    "        self.n_tensors = self.n_qubits\n",
    "\n",
    "    def construct_circuit(self, params=None):\n",
    "        ps = params if params is not None else self.params\n",
    "        assert len(ps) == self.n_params, 'Incorrect number of parameters!'\n",
    "        pc = cycle(ps)\n",
    "        circ = QuantumCircuit(self.q, self.c)\n",
    "        for i in range(self.n_qubits):\n",
    "            circ.ry(next(pc), self.q[i])\n",
    "            circ.rz(next(pc), self.q[i])\n",
    "        return circ\n",
    "    \n",
    "class Checkerboard(TNWithEntangler):\n",
    "    '''\n",
    "    Tensor network where gates are applied in a checkerboard pattern\n",
    "    '''\n",
    "\n",
    "    def __init__(self, q, c, entangler, params=None, depth=4):\n",
    "        super().__init__(q, c, entangler)\n",
    "        self.depth = depth\n",
    "        \n",
    "        if self.n_qubits % 2:\n",
    "            raise NotImplementedError\n",
    "        if self.n_qubits % 2 == 0:\n",
    "            self.n_tensors = (self.n_qubits // 2) * depth\n",
    "        else:\n",
    "            self.n_tensors = depth * (self.n_qubits // 2)\n",
    "        self.n_params = self.n_tensors * entangler.n_params\n",
    "        \n",
    "        if params is None:\n",
    "            self.params = [0] * self.n_params\n",
    "        else:\n",
    "            self.set_params(params)\n",
    "\n",
    "    def construct_circuit(self, params=None):\n",
    "        ps = params if params is not None else self.params\n",
    "        assert len(ps) == self.n_params, 'Incorrect number of parameters!'\n",
    "        pc = cycle(ps)\n",
    "        circ = QuantumCircuit(self.q, self.c)\n",
    "        def bulknext(pc):\n",
    "            return ([next(pc) for j in range(self.entangler.n_params)])\n",
    "        for layer in range(self.depth):\n",
    "            remainder = layer % 2\n",
    "            n_gates = (self.n_qubits) // 2\n",
    "            for gate in range(n_gates):\n",
    "                first_qubit = gate * 2 + remainder\n",
    "                second_qubit = (first_qubit + 1) % self.n_qubits\n",
    "                self.entangler.apply(self.q, circ, first_qubit, second_qubit, bulknext(pc))\n",
    "                circ.barrier()\n",
    "            circ.barrier()\n",
    "        return circ    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Entangler:\n",
    "    '''\n",
    "    Two-qubit gate used as an element of a TN state \n",
    "    '''\n",
    "    def __init__(self, params=None):\n",
    "        self.n_params = 2 \n",
    "        if params is None:\n",
    "            self.params = [0] * self.n_params\n",
    "        else:\n",
    "            self.set_params(params)\n",
    "\n",
    "    def set_params(self, params):\n",
    "        '''Set the parameters of the gate'''\n",
    "        if len(params) == self.n_params:\n",
    "            self.params = params\n",
    "        else:\n",
    "            raise ValueError('Incorrect amount of parameters!')\n",
    "        \n",
    "    def apply(self, q, circ, params=None):\n",
    "        '''\n",
    "        Apply with internal parameters if None\n",
    "        otherwise use those supplied'''\n",
    "        print('Applying parent method')\n",
    "        \n",
    "class IsingEntangler(Entangler):\n",
    "    '''\n",
    "    Entangler that uses the operators found in \n",
    "    the Ising model\n",
    "    <| (exp(X) \\otimes exp(X)) exp(ZZ) (exp(Z) \\otimes exp(Z)) \n",
    "    '''\n",
    "    def __init__(self, params=None):\n",
    "        self.n_params = 5\n",
    "        if params is None:\n",
    "            self.params = [0] * self.n_params\n",
    "        else:\n",
    "            self.set_params(params)\n",
    "    \n",
    "    def apply(self, q, circ, i, j, params=None):\n",
    "        '''\n",
    "        q - QuantumRegister\n",
    "        circ - QuantumCircuit\n",
    "        i, j - qubits to act on\n",
    "        '''\n",
    "        if params is None:\n",
    "            pc = cycle(self.params)\n",
    "        else:\n",
    "            if len(params) == self.n_params:\n",
    "                pc = cycle(params)\n",
    "            else:\n",
    "                raise ValueError('Incorrect number of parameters!')\n",
    "        circ.rx(next(pc), q[i])\n",
    "        circ.rx(next(pc), q[j])\n",
    "        circ.cx(q[i], q[j])\n",
    "        circ.u1(next(pc), q[j])\n",
    "        circ.cx(q[i], q[j])\n",
    "        circ.rz(next(pc), q[i])\n",
    "        circ.rz(next(pc), q[j])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hamiltonians\n",
    "\n",
    "In this tutorial, we try VQE on two examples: the Ising model with transverse field and Hamiltonians embedding SAT problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array([[0, 1], [1, 0]], dtype=np.complex64)\n",
    "Y = np.array([[0, -1j], [1j, 0]])\n",
    "Z = np.diag([1, -1])\n",
    "I = np.eye(2)\n",
    "\n",
    "def ising_model(n_spins, J, hx):\n",
    "    '''Ising model with transverse field'''\n",
    "    ham = {}\n",
    "    line = 'Z' + 'Z' + 'I' * (n_spins - 2)\n",
    "    for i in range(n_spins):\n",
    "        term = line[-i:] + line[:-i]\n",
    "        ham[term] = J\n",
    "    line = 'X' + 'I' * (n_spins - 1)\n",
    "    if hx != 0:\n",
    "        for i in range(n_spins):\n",
    "            term = line[-i:] + line[:-i]\n",
    "            ham[term] = hx\n",
    "    return ham\n",
    "\n",
    "def explicit_hamiltonian(ham_dict):\n",
    "    n_qubits = len(list(ham_dict.keys())[0])\n",
    "    #I = np.eye(2)\n",
    "    #X = np.array([[0, 1], [1, 0]])\n",
    "    #Y = np.array([[0, -1j], [1j, 0]])\n",
    "    #Z = np.diag([1, -1])\n",
    "    pauli={}\n",
    "    pauli['I'] = I\n",
    "    pauli['X'] = X\n",
    "    pauli['Y'] = Y\n",
    "    pauli['Z'] = Z\n",
    "    H = np.zeros((2**n_qubits, 2**n_qubits), dtype='complex128')\n",
    "    for term, energy in ham_dict.items():\n",
    "        matrices=[]\n",
    "        for sym in term:\n",
    "            matrices.append(pauli[sym])\n",
    "        total_mat = energy * reduce(np.kron, matrices)\n",
    "        H +=total_mat\n",
    "    return H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_clause(num_variables, k=3):\n",
    "    '''Generate a random clause for k-SAT. Returns a list of three\n",
    "    integers\n",
    "    '''\n",
    "    clause = np.random.choice(num_variables, k, replace=False) + 1\n",
    "    clause = [c * ((-1)**np.random.randint(2)) for c in clause]\n",
    "    return clause\n",
    "\n",
    "def random_SAT_instance(num_variables, num_clauses, k=3):\n",
    "    '''\n",
    "    Generate a random instance of K-SAT\n",
    "    '''\n",
    "    return [random_clause(num_variables, k) for i in range(num_clauses)]\n",
    "\n",
    "def make_SAT_Hamiltonian(sat_instance, num_bits, pauli=Z):\n",
    "    '''Creates a Hamiltonian from a SAT instance. pauli specifies which\n",
    "    basis to use\n",
    "    '''    \n",
    "    H = np.zeros((2**num_bits, 2**num_bits))\n",
    "    for clause in sat_instance:\n",
    "        term_list = [I] * num_bits\n",
    "        for i in clause:\n",
    "            term_list[abs(i) - 1] = 0.5 * (I + pauli * np.sign(i))\n",
    "        H += reduce(np.kron, term_list)\n",
    "    return H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variational quantum eigensolver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv_backend = Aer.get_backend(\"statevector_simulator\")\n",
    "\n",
    "def measure_ham(circ, ham_dict=None, explicit_H=None, shots=1000):\n",
    "    '''Measures the expected value of a Hamiltonian passed either as\n",
    "    ham_dict (uses QASM backend) or as explH (uses statevector backend)\n",
    "    '''\n",
    "    if ham_dict is not None:\n",
    "        raise NotImplementedError(\"Use statevector backend for now\")\n",
    "    elif explicit_H is not None:\n",
    "        job = execute(circ, sv_backend)\n",
    "        result = job.result()\n",
    "        state = result.get_statevector(circ)\n",
    "        state = np.array(state).reshape((len(state), 1))\n",
    "        E = (state.T.conj() @ explicit_H @ state)[0,0].real\n",
    "        return E\n",
    "    else:\n",
    "        raise TypeError('pass at least something!')\n",
    "        \n",
    "def any_order_VQE(TN: TensorNetwork, params_order: list,\n",
    "                    init_vals=None,\n",
    "                    ham_dict=None, explicit_H=None,\n",
    "                    initial_circuit=None,\n",
    "                    n_calls=100, verbose=True):\n",
    "    '''\n",
    "    Performs VQE by optimizing the tensor network in the order\n",
    "    supplied in params_order.\n",
    "    '''\n",
    "    if init_vals:\n",
    "        vals = [u for u in init_vals]\n",
    "    else:\n",
    "        vals = [0] * TN.n_params\n",
    "    \n",
    "    for free_parameters in params_order:\n",
    "        print('Optimizing parameters ', free_parameters, end=' ... ')\n",
    "        f = restrained_objective(TN, free_parameters, vals,\n",
    "                                   explicit_H=explicit_H, ham_dict=ham_dict,\n",
    "                                   initial_circuit=initial_circuit)\n",
    "        suggested_point = [vals[i] for i in free_parameters]\n",
    "\n",
    "        ## Supposedly the unitary exp(i F \\theta) yields a pi-periodic\n",
    "        ## cost function if F is a Pauli operator or any such that F**2 = 1\n",
    "        ## Which is sadly not always the case\n",
    "        res = gp_minimize(f, [(0,  2 * np.pi)] * len(free_parameters), n_calls=n_calls,\n",
    "                          x0=suggested_point, verbose=verbose)\n",
    "\n",
    "        \n",
    "        #print(res.x)\n",
    "        for i, n in enumerate(free_parameters):\n",
    "            vals[n] = res.x[i]\n",
    "        print('E = {0:0.4f}'.format(res.fun))\n",
    "        #print(['{0:0.6f}'.format(v) for v in vals])\n",
    "        #print(measure_ham_2(TN.construct_circuit(vals), explicit_H=explicit_H))\n",
    "\n",
    "    return res.fun, vals\n",
    "        \n",
    "def restrained_objective(TN, free_parameters, default_vals, ham_dict=None, explicit_H=None, initial_circuit=None):\n",
    "    '''Makes an objective function good for minimization. Locks most\n",
    "    parameters, while leaving those listed in free_parameters as free\n",
    "    '''\n",
    "\n",
    "    def f(x):\n",
    "        assert(len(x) == len(free_parameters)), 'Free parameters qty mismatch!'\n",
    "        params = [u for u in default_vals]  ## this is probably bad\n",
    "        for i, n in enumerate(free_parameters):\n",
    "            params[n] = x[i]\n",
    "        #print(params)\n",
    "        circ = TN.construct_circuit(params)\n",
    "        if initial_circuit:\n",
    "            circ = initial_circuit + circ\n",
    "        return measure_ham(circ, ham_dict=ham_dict, explicit_H=explicit_H)\n",
    "    return f        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running VQE on the examples\n",
    "\n",
    "The variational quantum eigensolver prepares a quantum state that depends on a number of parameters, samples energy of that state and uses the result as a value of the cost function. A classical optimizer then finds the next point to evaluate. In this example, the energy of the state is measured by explicitly using the Hamiltonian matrix, so it ignores the effects of finite measurement. The minimizer we use implements Bayesian optimization using Gaussian Processes.\n",
    "\n",
    "Here we use two ansatz states: unentangled state (<tt>RankOne</tt>) and checkerboard tensor network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_qubits = 6\n",
    "\n",
    "\n",
    "## SAT Hamiltonian\n",
    "sat_instance = random_SAT_instance(n_qubits, 23)\n",
    "H = make_SAT_Hamiltonian(sat_instance, n_qubits) \n",
    "\n",
    "## Ising Hamiltonian\n",
    "H_dict = ising_model(n_qubits, 1, -1)\n",
    "H_2 = explicit_hamiltonian(H_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimizing parameters  [0, 1] ... E = 1.7653\n",
      "Optimizing parameters  [2, 3] ... E = 0.6719\n",
      "Optimizing parameters  [4, 5] ... E = -1.6164\n",
      "Optimizing parameters  [6, 7] ... E = -2.9814\n",
      "Optimizing parameters  [8, 9] ... E = -4.9873\n",
      "Optimizing parameters  [10, 11] ... E = -5.2911\n",
      "Optimizing parameters  [0, 1] ... E = -5.4414\n",
      "Optimizing parameters  [2, 3] ... E = -6.1748\n",
      "Optimizing parameters  [4, 5] ... E = -6.2143\n",
      "Optimizing parameters  [6, 7] ... E = -6.9470\n",
      "Optimizing parameters  [8, 9] ... E = -7.0682\n",
      "Optimizing parameters  [10, 11] ... E = -7.0682\n"
     ]
    }
   ],
   "source": [
    "q, c = QuantumRegister(n_qubits), ClassicalRegister(n_qubits)\n",
    "circ = QuantumCircuit(q, c)\n",
    "\n",
    "ent = IsingEntangler()\n",
    "\n",
    "TN = RankOne(q, c)\n",
    "n_block_pars = 2\n",
    "\n",
    "#TN = Checkerboard(q, c, ent, depth=4)\n",
    "#n_block_pars = ent.n_params\n",
    "\n",
    "params_order = [list(range(k * n_block_pars, (k + 1) * n_block_pars))\n",
    "                        for k in range(TN.n_tensors)] * 2\n",
    "\n",
    "E, vals = any_order_VQE(TN, params_order, \n",
    "                          explicit_H=H_2, n_calls=20, \n",
    "                          verbose=False, init_vals=None)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3 (quantumenv)",
   "language": "python",
   "name": "quantumenv"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}