Knapsack.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The 0-1 Knapsack Problem solved with a VQE\n",
    "\n",
    "[The Knapsack Problem](https://en.wikipedia.org/wiki/Knapsack_problem) is an optimzation problem that goes like this: You are given a list of items that each has a weight and a value. You also have a knapsack that can hold a maximum weigh. Determine which items to take in the knapsack so as to maximize the total value taken without going over the maiximum weight.\n",
    "\n",
    "The most efficient approach would be a greedy approach, but that is not guaranteed to give the best result. For example let's take the these items\n",
    " - a tv of weight 12 and value 75\n",
    " - a watch of weight 3 and value 50\n",
    " - a vase of weight 10 and value 40\n",
    " - a ball of weight 5 and value 10\n",
    " \n",
    "Taking the most valuable product (the tv and the watch) would give us a total value of 125. Getting the lightest objects (the watch, the ball and the vase) gives us a value 100. Getting the items with the best value to weight ratio (the watch and the tv) gives us 125. The best choice is to take the tv, the watch and the ball with a total value of 135.\n",
    "\n",
    "## The VQE approach\n",
    "\n",
    "The VQE algorithm can find the ground state of an Ising Hamiltonian. In order to solve an optimization problem using this algorithm we need to \n",
    " 1. first express the solution to our problem as a binary string. \n",
    " 2. Then we need to find a cost function that acts on binary strings and that reachs it minimal value at the point wehere the binary string is the solution to the optimization problem. I.e. If we minimize the cost function and get the input to the minimum value that input is our solution. \n",
    " 3. We then need to map our function to an Ising Hamiltonian\n",
    " 4. Use the VQE Algorithm to find the ground state of that Hamiltonian\n",
    "\n",
    "We are done - the ground state that we found is the solution to our problem.\n",
    "\n",
    "Fortunately for step 3 there is a cookie-cutter technique to mapping a cost function to an Ising Hamiltonian and step 4 is built into Qiskit. So we really only need to worry about step 1 and 2.\n",
    "\n",
    "### The cost function\n",
    "\n",
    "For the 0-1 Knapsack Problem the solution can be expressed as an array $X =[x_{0}..x_{n-1}]$ of 0s and 1s in the following way: Each item is represented at an index of the array. If the value is 1 then we take the item in the knapsac. If it's 0 we don't take it.\n",
    "\n",
    "We can represent our weights as an array $W =[w_{0}..w_{n-1}]$ and the values as an array $V =[v_{0}..v_{n-1}]$ and name the maximum weight $ W_{max} $. \n",
    "\n",
    "To build the cost function we must keep in mind that we want to maximize the sum $ \\sum_{i} x_{i}v_{i} $ while comlying with the constraint that $ W_{max} >= \\sum_{i} x_{i}w_{i} $. This cost function is difficult to express so we use a mathematical trick. We add another term  $ S $ that will vary with our solution. In order to make it change wit the solution we enhance $ X $ with a number of aditional bits $ X' = [x_0,..x_{n-1},y_{n}..y_{n+m-1}]  $. The bytes $ y_i $ will be the binary representation of $ S $. In this way the optimizer will be able to optimize $ S $ as well as $ X $\n",
    "\n",
    "We can now write out cost function\n",
    "\n",
    "$$ C(X') = M(W_{max} - \\sum_{i=0}^{n-1} x_{i}w_{i} - S)^2 - \\sum_{i}^{n-1} x_{i}v_{i} $$ where $ S = \\sum_{j=n}^{n+m} 2^jy_{j} $ and $M$ is a number large enought to dominate the sum of values. \n",
    "\n",
    "We see that, because $S$ can only be positive, when $W_{max} >= \\sum_{i} x_{i}w_{i}$ the optimizer can find an $S$ (or better the $y_j$ that compose $S$) so that it will take the first term to 0. This way the function is dominated by the sum of values. If $W_{max} < \\sum_{i} x_{i}w_{i}$ then the first term can never be 0 and, multiplied by a large $M$, will always dominate the function.\n",
    "The minimum value of the function will be that where the constraint is respected and the sum of values is maximized.\n",
    "\n",
    "Note that because $[y_n..y_{n_m-1}]$ is the binary representation of $S$ and $S$ is at most $W_{max}$ them $m = \\log(W_{max})$\n",
    "\n",
    "\n",
    "### Mapping to an Ising Hamiltonian\n",
    "\n",
    "To map this function to a Hamiltonian we use a known technique. We first expand the first term. We then ignore the constant terms since they cannot be minimized. We then take each term of the remaining sum and apply the mapping $x_i\\rightarrow (I-Z_i)/2$, where $I$ is the identity matrix of order $n+m$ and $Z_i$ is the matrix that results from tensoring $n+m$ matrixes of which at every positin other than #i# the factors of the tensor are identity matrixes ans at position $i$ is the Pauly $Z$ matrix.\n",
    "\n",
    "This is the cookie-cutter technique I mentioned earlier. And luckily we don't even need to do the tensor products ourselves. Qiskit has an object that takes the inde of the Pauly $Z$ in the tensor and does the rest.\n",
    "\n",
    "Afer doing this replacement we expand again and remove the constant terms. The rest we need to return as a list of terms to the VQE algorithm.\n",
    "\n",
    "### Getting the solution\n",
    "\n",
    "The VQE minimizes the function and gives us $X'$. All we need to do is thake the first $n$ terms and we have our solution\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# useful additional packages \n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.axes as axes\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "\n",
    "from qiskit import BasicAer\n",
    "from qiskit.tools.visualization import plot_histogram\n",
    "from qiskit.aqua import run_algorithm\n",
    "from qiskit.aqua.input import EnergyInput\n",
    "from qiskit.aqua.translators.ising import max_cut, tsp\n",
    "from qiskit.aqua.algorithms import VQE, ExactEigensolver\n",
    "from qiskit.aqua.components.optimizers import SPSA\n",
    "from qiskit.aqua.components.variational_forms import RY\n",
    "from qiskit.aqua import QuantumInstance\n",
    "\n",
    "from qiskit.quantum_info import Pauli\n",
    "from qiskit.aqua.operators import WeightedPauliOperator\n",
    "from collections import OrderedDict\n",
    "import math\n",
    "\n",
    "# setup aqua logging\n",
    "import logging\n",
    "from qiskit.aqua import set_qiskit_aqua_logging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining the terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_knapsack_operator(values, weights, max_weight):\n",
    "    if len(values) != len(weights):\n",
    "        raise ValueError(\"The values and weights must have the same length\")\n",
    "\n",
    "    if any(v < 0 for v in values) or any(w < 0 for w in weights):\n",
    "        raise ValueError(\"The values and weights cannot be negative\")\n",
    "\n",
    "    if all(v == 0 for v in values):\n",
    "        raise ValueError(\"The values cannot all be 0\")\n",
    "\n",
    "    if max_weight < 0:\n",
    "        raise ValueError(\"max_weight cannot be negative\")\n",
    "\n",
    "    y_size = int(math.log(max_weight, 2)) + 1 if max_weight > 0 else 1\n",
    "    print(y_size)\n",
    "    n = len(values)\n",
    "    num_values = n + y_size\n",
    "    pauli_list = []\n",
    "    shift = 0\n",
    "\n",
    "    M = 2000000  #10 * np.sum(values)\n",
    "\n",
    "    # term for sum(x_i*w_i)**2\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            coefficient = -1 * 0.25 * weights[i] * weights[j] * M\n",
    "\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[j] = not zp[j]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[i] = not zp[i]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "            coefficient = -1 * coefficient\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[i] = not zp[i]\n",
    "            zp[j] = not zp[j]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "    # term for sum(2**j*y_j)**2\n",
    "    for i in range(y_size):\n",
    "        for j in range(y_size):\n",
    "            coefficient = -1 * 0.25 * (2 ** i) * (2 ** j) * M\n",
    "\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[n + j] = not zp[n + j]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[n + i] = not zp[n + i]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "            coefficient = -1 * coefficient\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[n + i] = not zp[n + i]\n",
    "            zp[n + j] = not zp[n + j]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "    # term for -2*W_max*sum(x_i*w_i)\n",
    "    for i in range(n):\n",
    "        coefficient = max_weight * weights[i] * M\n",
    "\n",
    "        xp = np.zeros(num_values, dtype=np.bool)\n",
    "        zp = np.zeros(num_values, dtype=np.bool)\n",
    "        zp[i] = not zp[i]\n",
    "        pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "        shift -= coefficient\n",
    "\n",
    "    # term for -2*W_max*sum(2**j*y_j)\n",
    "    for j in range(y_size):\n",
    "        coefficient = max_weight * (2 ** j) * M\n",
    "\n",
    "        xp = np.zeros(num_values, dtype=np.bool)\n",
    "        zp = np.zeros(num_values, dtype=np.bool)\n",
    "        zp[n + j] = not zp[n + j]\n",
    "        pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "        shift -= coefficient\n",
    "\n",
    "    for i in range(n):\n",
    "        for j in range(y_size):\n",
    "            coefficient = -1 * 0.5 * weights[i] * (2 ** j) * M\n",
    "\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[n + j] = not zp[n + j]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[i] = not zp[i]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "            coefficient = -1 * coefficient\n",
    "            xp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp = np.zeros(num_values, dtype=np.bool)\n",
    "            zp[i] = not zp[i]\n",
    "            zp[n + j] = not zp[n + j]\n",
    "            pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "            shift -= coefficient\n",
    "\n",
    "    # term for sum(x_i*v_i)\n",
    "    for i in range(n):\n",
    "        coefficient = 0.5 * values[i]\n",
    "\n",
    "        xp = np.zeros(num_values, dtype=np.bool)\n",
    "        zp = np.zeros(num_values, dtype=np.bool)\n",
    "        zp[i] = not zp[i]\n",
    "        pauli_list.append([coefficient, Pauli(zp, xp)])\n",
    "        shift -= coefficient\n",
    "\n",
    "    return WeightedPauliOperator(paulis=pauli_list), shift\n",
    "\n",
    "\n",
    "def get_solution(x, values):\n",
    "    return x[:len(values)]\n",
    "\n",
    "\n",
    "def knapsack_value_weight(solution, values, weights):\n",
    "    value = np.sum(solution * values)\n",
    "    weight = np.sum(solution * weights)\n",
    "    return value, weight\n",
    "\n",
    "def sample_most_likely(state_vector):\n",
    "    if isinstance(state_vector, dict) or isinstance(state_vector, OrderedDict):\n",
    "        # get the binary string with the largest count\n",
    "        binary_string = sorted(state_vector.items(), key=lambda kv: kv[1])[-1][0]\n",
    "        x = np.asarray([int(y) for y in reversed(list(binary_string))])\n",
    "        return x\n",
    "    else:\n",
    "        n = int(np.log2(state_vector.shape[0]))\n",
    "        k = np.argmax(np.abs(state_vector))\n",
    "        x = np.zeros(n)\n",
    "        for i in range(n):\n",
    "            x[i] = k % 2\n",
    "            k >>= 1\n",
    "        return x\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking the solution with an exact eigensolver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "[1. 0. 0. 0.]\n",
      "result=[1. 0.]\n"
     ]
    }
   ],
   "source": [
    "# values = [680, 120, 590, 178]\n",
    "# weights = [13, 6, 15, 9]\n",
    "# w_max = 32\n",
    "\n",
    "values = [8, 2]\n",
    "weights = [3, 1]\n",
    "w_max = 3\n",
    "\n",
    "qubitOp, offset = get_knapsack_operator(values, weights, w_max)\n",
    "\n",
    "algo_input = EnergyInput(qubitOp)\n",
    "ee = ExactEigensolver(qubitOp, k=1)\n",
    "result = ee.run()\n",
    "\n",
    "most_lightly = result['eigvecs'][0]\n",
    "x = sample_most_likely(most_lightly)\n",
    "\n",
    "print(x)\n",
    "print('result=' + str(x[:len(values)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running with a statevector simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 1. 0. 1.]\n",
      "result usig statevector_simulator =[0. 1.]\n"
     ]
    }
   ],
   "source": [
    "seed = 10598\n",
    "\n",
    "spsa = SPSA(max_trials=10000)\n",
    "ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')\n",
    "vqe = VQE(qubitOp, ry, spsa)\n",
    "\n",
    "backend = BasicAer.get_backend('statevector_simulator')\n",
    "quantum_instance = QuantumInstance(backend, seed_simulator=seed, seed_transpiler=seed)\n",
    "\n",
    "result_statevector = vqe.run(quantum_instance)\n",
    "\n",
    "most_lightly_sv = result_statevector['eigvecs'][0]\n",
    "x_statevector = sample_most_likely(most_lightly_sv)\n",
    "\n",
    "print(x_statevector)\n",
    "print('result usig statevector_simulator =' + str(x_statevector[:len(values)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 1 0 0 1]\n",
      "result usig qasm_simulator =[0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "# run quantum algorithm with shots\n",
    "seed = 10598\n",
    "\n",
    "spsa = SPSA(max_trials=1000)\n",
    "ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')\n",
    "vqe = VQE(qubitOp, ry, spsa)\n",
    "\n",
    "backend = BasicAer.get_backend('qasm_simulator')\n",
    "quantum_instance = QuantumInstance(backend, shots=1024, seed_simulator=seed, seed_transpiler=seed)\n",
    "\n",
    "result_shots = vqe.run(quantum_instance)\n",
    "\n",
    "most_lightly_shots = result_shots['eigvecs'][0]\n",
    "x_shots = sample_most_likely(most_lightly_shots)\n",
    "\n",
    "print(x_shots)\n",
    "print('result usig qasm_simulator =' + str(x_shots[:len(values)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "math.log(8, 2)"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}