tutorials/tutorial_02_solving_qubos.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7dca7602",
   "metadata": {},
   "source": [
    "*Copyright (C) 2022 Intel Corporation*<br>\n",
    "*SPDX-License-Identifier: BSD-3-Clause*<br>\n",
    "*See: https://spdx.org/licenses/*\n",
    "\n",
    "---\n",
    "\n",
    "# Quadratic Unconstrained Binary Optimization (QUBO) with Lava\n",
    "\n",
    "In previous work, we found that our first generation Loihi 1 chip can solve QUBOs more than 17x faster and more than 670x more energy efficient than the CPU-based solver QBSolv [1], which \"executes in two-thirds of the time of the best previously known implementation\" [2]. Lava translates this solver to Loihi 2 to exploit this chip's higher speed, efficiency, and new hardware features designed for optimization.\n",
    "\n",
    "This tutorial ramps up users who want to use Lava's solver for quadratic unconstraint binary optimization (QUBO) problems. It starts with an intuitive description of the algorithm behind the solver, before it exemplifies how to encode and solve a QUBO workload with 700 variables on CPU and Loihi2. An outlook will finally give a glimpse into the future feature releases of the Lava solver. \n",
    "\n",
    "An accompanying tutorial will soon be released to provide a more technical deep dive into Lava's QUBO solver."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33f5ee96",
   "metadata": {},
   "source": [
    "## Recommended tutorials before starting\n",
    "\n",
    "- [Installing Lava](https://github.com/lava-nc/lava/blob/main/tutorials/in_depth/tutorial01_installing_lava.ipynb \"Tutorial on Installing Lava\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c190a5a6",
   "metadata": {},
   "source": [
    "## Set up the environment"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7617dfdb",
   "metadata": {},
   "source": [
    "To solve QUBOs in Lava, we import the corresponding modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0864507",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interface for QUBO problems\n",
    "from lava.lib.optimization.problems.problems import QUBO\n",
    "# Generic optimization solver\n",
    "from lava.lib.optimization.solvers.generic.solver import OptimizationSolver, solve, SolverConfig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba620b77",
   "metadata": {},
   "source": [
    "In addition, we import auxiliary modules to generate the workloads and run the solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f255560d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import networkx as ntx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73d93682",
   "metadata": {},
   "source": [
    "If Loihi 2 hardware is available, we can take advantage of the large speed and energy efficiency of this chip to solve QUBOs. To access the chip, we must configure the following environment variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3338d93b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lava.utils.system import Loihi2\n",
    "Loihi2.preferred_partition = \"oheogulch\"\n",
    "loihi2_is_available = Loihi2.is_loihi2_available\n",
    "\n",
    "if loihi2_is_available:\n",
    "    # Enable SLURM, the workload manager used to distribute Loihi2 resources to users\n",
    "    os.environ['SLURM'] = '1'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b10596f",
   "metadata": {},
   "source": [
    "## An intuitive description of Lava's QUBO solver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab5535e0",
   "metadata": {},
   "source": [
    "A quadratic unconstrained binary optimization (QUBO) task is an NP-hard combinatorial optimization problem with a plethory of commercial applications, as exemplified later. Loihi offers a range of features useful for their fast and energy efficient solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d933e4ba",
   "metadata": {},
   "source": [
    "#### The QUBO formalism"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8e1426d",
   "metadata": {},
   "source": [
    "The goal of a QUBO is to minimize the cost function\n",
    "\n",
    "&emsp;&emsp; $\\underset{x}{\\text{min}}\\; \\mathbf{x}^T \\mathbf{Q} \\mathbf{x} $ ,\n",
    "\n",
    "which is subject to no constraints. This equation is defined by the QUBO matrix $\\mathbf{Q}\\in\\mathbb{R}^{n\\times n}$ and the vector of binary variables, $\\mathbf{x} \\in \\{0,1\\}^n$ ."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e94a9fa7",
   "metadata": {},
   "source": [
    "#### Lava's QUBO solver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56046a73",
   "metadata": {},
   "source": [
    "To solve QUBOs, Lava currently implements a spiking version of a Boltzmann machine that was inspired by the work of Jonke et al. [2].\n",
    "\n",
    "To give a rough intuition of the solver, it represents binary variables as neurons that either spike (variable value=1) or are silent (value=0).\n",
    "\n",
    "<br>\n",
    "\n",
    "<img src=\"https://user-images.githubusercontent.com/86950058/192372805-b974e39b-fae8-4a68-810b-f52a9363987b.png\" width=\"500\" align=\"center\"/>\n",
    "    \n",
    "Each neuron has an internal state representing the probability that the neuron should spike. When the state exceeds a threshold, the neuron spikes and the variable is thus switched to 1. Once a spike occurs at neuron $i$, the synapses forward the spikes according to the off-diagonal QUBO weights $Q_{ij}$ to other neurons $j$. Negative weights increase the chance that the connected neuron will spike, positive weights decrease the chance. <br>\n",
    "On-diagonal weights $Q_{ii}$ add a bias to the neurons state variable, which continuously increases the chance that the neuron will spike.\n",
    "\n",
    "<br>\n",
    "\n",
    "<img src=\"https://user-images.githubusercontent.com/86950058/192372894-a6f757e5-a6bf-47bb-8254-08de9188d029.png\" width=\"800\"  align=\"center\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38684903",
   "metadata": {},
   "source": [
    "## Example application: finding the Maximum Independent Set of a graph"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b3dcc34",
   "metadata": {},
   "source": [
    "The following example shows how to formulate and solve a problem as QUBO, using the NP-hard maximum independent set (MIS) problem from graph theory. The goal of an MIS task is to find the largest subset of vertices in a graph that are mutually unconnected. In the following graph, the purple nodes form such a set:\n",
    "\n",
    "<br>\n",
    "\n",
    "<img src=\"https://user-images.githubusercontent.com/86950058/192372990-ec4e5926-463c-4b30-810d-08a896446d8a.png\" width=\"250\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fa5e8f6",
   "metadata": {},
   "source": [
    "Hardware accelerators for this problem class are subject to substantial research efforts [3, 4] due to the large industrial need, as recently summarized by Wurtz et al. [5] and exemplified in the following figure."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e096240f",
   "metadata": {},
   "source": [
    "<img src=\"https://user-images.githubusercontent.com/86950058/197763251-c6e5929d-587e-400a-8392-f398d5176f7e.png\" width=\"800\"/>\n",
    "\n",
    "*a) For wireless communication such as 5G or WiFi, access points and sensors are grouped into small networks. To avoid interference, MIS can help to divide the networks into almost interference-free groups that may use the same frequencies [6]. <br>\n",
    "b) To find the largest error-correcting code, each binary vector $u_i$ is considered as the vertex of a graph. Two vertices $u_i$ and $u_j$ are joined by an edge if both vectors may end up in the same state if acted upon by certain errors, such as deletion or transposition of bits. The MIS describes the largest code that is robust to such errors [7]. <br>\n",
    "c) In the design of semiconductor chips, a via connects two adjacent net segments in adjacent metal layers. Redundant vias are inserted to tolerate single via failures and thus improving the manufacturing yield [8]. If a specific via can be replaced by a redundant double via, it is modeled as a vertex. If two vias cannot be simultaneously replaced by double vias due to design rule violations, their vertices are connected by an edge. The MIS determines the most reliable chip design with the maximum number of redundant vias.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4d8f8a4",
   "metadata": {},
   "source": [
    "The following code determines the MIS for a graph with 1000 nodes. For this, it translates the MIS problem into a QUBO workload, and then solves it in Lava."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5ff96c0",
   "metadata": {},
   "source": [
    "#### Define the QUBO workload"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a497c70",
   "metadata": {},
   "source": [
    "Lava provides an easy interface to encode QUBO problems, by providing the QUBO matrix $\\mathbf{Q}$.\n",
    "\n",
    "Please note:\n",
    "As shown below, the QUBO solver can run on CPU and Loihi 2 backend. At the moment, the Loihi 2 model is algorithmically more advanced and more performant. For illustration purposes, the following cell thus instantiates a less complex workload for the CPU backend. In the future, the algorithms for CPU and Loihi 2 will be aligned again. At that point, both backends will be able to solve workloads to the same optimality, while Loihi 2 will remain faster and more energy efficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "411835c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import utility functions to create and analyze MIS workloads\n",
    "from lava.lib.optimization.utils.generators.mis import MISProblem\n",
    "\n",
    "# Create an undirected graph with num_vertices vertices and a \n",
    "# probability of connection_prob that any two vertices are randomly connected\n",
    "if loihi2_is_available:\n",
    "    mis = MISProblem(num_vertices=1000, connection_prob=0.15, seed=0)\n",
    "else:\n",
    "    mis = MISProblem(num_vertices=15, connection_prob=0.9, seed=0)\n",
    "\n",
    "# Translate the MIS problem for this graph into a QUBO matrix\n",
    "q = mis.get_qubo_matrix(w_diag=1, w_off=8)\n",
    "\n",
    "# Create the qubo problem\n",
    "qubo_problem = QUBO(q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c570e079",
   "metadata": {},
   "source": [
    "#### Solve the QUBO on a CPU or Loihi 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5b888df",
   "metadata": {},
   "source": [
    "The QUBO matrix is provided to Lava's generic _OptimizationSolver_ for constraint optimization problems. A call to its _solve_ method then tackles the workload.\n",
    "\n",
    "As stopping condition, you can provide the solver with a maximum number of time steps (_timeout_) and/or a target cost.\n",
    "\n",
    "Lava will print an obtained solution whenever it found a local minimum. Then it continues its search until either stopping condition is fullfilled.\n",
    "\n",
    "The following cell will run the solver either on a CPU backend, or on a Loihi 2 chip if available. As mentioned above, Loihi 2 has a more advanced algorithm, and thus requires a different set of hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a672adae",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = OptimizationSolver(qubo_problem)\n",
    "\n",
    "# Provide hyperparameters for the solver\n",
    "# Guidance on the hyperparameter search will be provided in the deep dive tutorial\n",
    "np.random.seed(0)\n",
    "if loihi2_is_available:\n",
    "    target_cost = int(-50)\n",
    "    timeout = 200000\n",
    "    hyperparameters = {\n",
    "        \"neuron_model\": 'nebm-sa-refract',\n",
    "        'annealing_schedule': 'linear',\n",
    "        \"max_temperature\": int(15),\n",
    "       \"min_temperature\": int(0),\n",
    "       \"steps_per_temperature\": int(20),\n",
    "       \"delta_temperature\": int(1),\n",
    "       \"refract_scaling\": int(4),\n",
    "       \"refract\": int(0),\n",
    "       \"refract_counter\": 0,\n",
    "       \"init_value\": np.random.randint(0, 2, size=(qubo_problem.num_variables)),}\n",
    "else:\n",
    "    target_cost = int(-2)\n",
    "    timeout = 100\n",
    "    hyperparameters = {\n",
    "        \"neuron_model\": 'nebm',\n",
    "        \"temperature\": int(1),}\n",
    "\n",
    "# Define the backend\n",
    "if loihi2_is_available:\n",
    "    backend = 'Loihi2'\n",
    "else:\n",
    "    backend = 'CPU'\n",
    "\n",
    "# Solve the QUBO using Lava's OptimizationSolver on CPU\n",
    "# Change \"backend='Loihi2'\" if your system has physical access to this chip\n",
    "solver_report = solver.solve(\n",
    "    config=SolverConfig(\n",
    "        timeout=timeout,\n",
    "        hyperparameters=hyperparameters,\n",
    "        target_cost=target_cost,\n",
    "        backend=backend\n",
    "    )\n",
    ")\n",
    "\n",
    "solution_loihi = solver_report.best_state\n",
    "print(f'\\nSolution of the provided QUBO: {np.where(solution_loihi == 1.)[0]}.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9112d43a",
   "metadata": {},
   "source": [
    "To compare the solution with the optimal solution, we can compare it with the result of Python's Networkx package, which solves the underlying graph theoretical problem directly. For Loihi 2, the code is commented out as Networkx takes very long to find a solution to the optimization workload."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "440201d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not loihi2_is_available:\n",
    "    # Find the optimal solution to the MIS problem\n",
    "    solution_opt = mis.find_maximum_independent_set()\n",
    "\n",
    "    # Calculate the QUBO cost of the optimal solution\n",
    "    cost_opt = qubo_problem.evaluate_cost(solution=solution_opt)\n",
    "\n",
    "    # Calculate the QUBO cost of Lava's solution\n",
    "    cost_lava = qubo_problem.evaluate_cost(solution=solution_loihi)\n",
    "\n",
    "    print(f'QUBO cost of solution: {cost_lava} (Lava) vs. {cost_opt} (optimal)\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c296054",
   "metadata": {},
   "source": [
    "## Future features "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d08828c",
   "metadata": {},
   "source": [
    "Future releases will provide the following advancements and many more:\n",
    "- Compiler support for larger problems and faster problem compilation.\n",
    "- Performance improvements using advanced hyperparameter search, improved noise annealing schedules, additional algorithms, and parallel runs of solver instances."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccef5697",
   "metadata": {},
   "source": [
    "## How to learn more?\n",
    "\n",
    "We will soon release an accompanying deep dive tutorial to explain how...\n",
    "- the QUBO solver works under the hood.\n",
    "- optimization tasks can be encoded as QUBO.\n",
    "- hyperparameters affect the performance of the solver.\n",
    "- to benchmark solution accuracy, run time, and energy consumption between Loihi2 and CPU-based solvers.\n",
    "\n",
    "Watch this [space](https://github.com/lava-nc/lava-optimization) to learn about upcoming developments to the QUBO solver and the optimization toolbox in Lava in general. \n",
    "\n",
    "If you want to find out more about the implementation of the QUBO solver, have a look at the [Lava documentation](https://lava-nc.org/ \"Lava Documentation\") or dive into the [source code](https://github.com/lava-nc/lava-optimization/tree/main/src/lava/lib/optimization/solvers/generic/solver.py\n",
    "\"Source code of the generic constraint optimization solver\").\n",
    "\n",
    "To receive regular updates on the latest developments and releases of the Lava Software Framework please subscribe to the [INRC newsletter](http://eepurl.com/hJCyhb \"INRC Newsletter\")."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1a075a2",
   "metadata": {},
   "source": [
    "\n",
    "## References\n",
    "\n",
    "[1] Benchmarks ran with [QBSolv](https://github.com/dwavesystems/qbsolv) on an  Intel Core i9-7920X CPU @ 2.90GHz with 132GB DRAM. <br>\n",
    "Performance varies by use, configuration and other factors. Learn more at www.Intel.com/PerformanceIndex. <br>\n",
    "Performance results are based on testing as of dates shown in configurations and may not reflect all publicly available updates.  See backup for configuration details.  No product or component can be absolutely secure. \n",
    "\n",
    "[2] D-Wave Systems Inc: Booth, Reinhardt, Roy, _Partitioning Optimization Problems for Hybrid Classical/Quantum Execution_. [Technical report](https://docs.ocean.dwavesys.com/projects/qbsolv/en/latest/_downloads/bd15a2d8f32e587e9e5997ce9d5512cc/qbsolv_techReport.pdf \"DWave's technical report\"), 2017. <br>\n",
    "\n",
    "[3] Mallick, A., Bashar, M.K., Truesdell, D.S., Calhoun, B.H., Joshi, S., Shukla, N. Using synchronized oscillators to compute the maximum independent set. Nature Communications, 11, 2020. <br>\n",
    "\n",
    "[4] Ebadi, S., Keesling, A., Cain, M., Wang, T.T., Bluvstein, D., Semeghini, G., Omran, A., Liu, J.-G., Samajdar, R., Luo, X.-Z., Nash, B., Gao, X., Barak, B., Farhi, E., Sachdev, S., Gemelke, N., Zhou, L., Choi, S., Pichler, H., Wang, S.-T., Greiner, M., Vuletic, V.,. Lukin, M.D. Quantum optimization of maximum independent set using Rydberg atom arrays. Science, 376, 6598, 2022. <br>\n",
    "\n",
    "[5] Wurtz, J., Lopes, P., Gemelke, N., Keesling, A., Wang, S. Industry applications of neutral-atom quantum computing solving independent set problems. arXiv:2205.08500v1, 2022. <br>\n",
    "\n",
    "[6] Zhou, J., Wang, L., Wang, W. & Zhou, Q. Efficient graph-based resource allocation scheme using maximal independent set for randomly-deployed small star networks. Sensors 17, 2553 (2017). <br>\n",
    "\n",
    "[7] Butenko, S., Pardalos, P., Sergienko, I., Shylo, V. & Stetsyuk, P. Finding maximum independent sets in graphs arising from coding theory. Proceedings of the 17th ACM Symposium on Applied Computing, 542–546, 2002. <br>\n",
    "\n",
    "[8] Lee, K. Y. & Wang, T. C. Post-routing redundant via insertion for yield/reliability improvement. In Proceedings of the 2006 Asia and South Pacific Design Automation Conference, 303–308, 2006. <br>"
   ]
  }
 ],
 "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.6.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "0285f4d976e7d5c77fbc79f3712c5561c504e187cba953573dd269218b32f868"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}