Merge pull request #29 from ICHEC/rt

Rt

lecture-11/lecture-11-note.ipynb

@@ -0,0 +1,244 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Quadratic unconstrained binary optimization (QUBO)\n",
+    "\n",
+    "## Example 1\n",
+    "\n",
+    "Consider a set of N positive integers (though one doesn't need to) $U = \\{n_1, n_2, n_3, \\dots, n_N\\}$, and imagine the goal is to construct two disjoint subsets of this, called A and B, such that $A\\cup B=U$ and $A\\cap B=\\empty$, such that the sum of the numbers in each subset is equal or as close to each other as possible.\n",
+    "\n",
+    "Formally, the sets A and B are to be constructed so that\n",
+    "\n",
+    "$$f(\\{n_i\\}) = (S_A - S_B)^2~~\\text{where}$$\n",
+    "$$S_A = \\sum_{i\\in A} n_i, \\quad S_B = \\sum_{i\\in B}n_i$$\n",
+    "\n",
+    "is minimum. This is a `cost function` that has to be minimized.\n",
+    "\n",
+    "What we do is, assign to each number $n_i$ a binary variable $x_i$, such that $x_i = 1$ if $n_i\\in A$ or otherwise $x_i = 0$ if $n_i \\in B$.\n",
+    "\n",
+    "Now we can write the partial sum as\n",
+    "$$\n",
+    "S_A = \\sum_{i\\in U} n_i x_i \\\\S_B = \\sum_{i\\in U}n_i (1-x_i)\n",
+    "$$\n",
+    "\n",
+    "The cost function can now be written as \n",
+    "\n",
+    "\\begin{align}\n",
+    "f(\\{n_i\\}) &= (S_A - S_B)^2 = S_A^2 + S_B^2 - 2S_A S_B\\\\\n",
+    "&= (\\sum_{i\\in U} n_i x_i)^2 + (\\sum_{i\\in U}n_i (1-x_i))^2 - 2 (\\sum_{i\\in U} n_i x_i)\\sum_{j\\in U}n_j (1-x_j)\\\\\n",
+    "&= \\sum_{ij\\in U}\\left[\n",
+    "    n_i n_j x_i x_j + n_i n_j (1-x_i)(1-x_j) - 2 n_i n_j x_i (1-x_j)\n",
+    "    \\right]\\\\\n",
+    "    &= \\sum_{ij\\in U} n_i n_j (x_i x_j + 1 - x_i - x_j + x_i x_j - 2x_i + 2 x_i x_j)\\\\\n",
+    "    &= 4\\sum_{ij\\in U} n_i n_j x_i x_j - 4\\sum_{ij\\in U}(n_i n_j x_i) + N^2\\\\\n",
+    "    &= 4\\sum_{ij\\in U} n_i n_j x_i x_j - 4S\\sum_{i\\in U}n_i x_i + N^2 \n",
+    "\\end{align}\n",
+    "\n",
+    "Where $S_U = \\sum_{i\\in U} n_i$ is the sum of all the numbers in the set U. Thus we need to minimize the `f` over the discrete space of vector **x** $ = \\{x_1, x_2, \\dots, x_N\\}$. Note that the cost function `f` is quadratic in $x_i$\n",
+    "\n",
+    "$$\n",
+    "f(x) = 4\\sum_{ij\\in U} n_i n_j x_i x_j - 4S\\sum_{i\\in U}n_i x_i + N^2 \n",
+    "$$\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The cost function can be expressed in matrix form as\n",
+    "\n",
+    "$$\n",
+    "f(x) = 4 x^{T}Qx + V^{T}x + N^2\n",
+    "$$\n",
+    "\n",
+    "where $Q_{ij} = 4n_i n_j$ and $V_i = 4S n_i$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from typing import Any\n",
+    "import numpy as np\n",
+    "\n",
+    "class qubo:\n",
+    "\n",
+    "    def __init__(self, ni: np.ndarray) -> None:\n",
+    "        self.ni = np.asarray(ni, dtype=int)\n",
+    "        self.Q = 4 * np.outer(self.ni, self.ni)\n",
+    "        self.sum = self.ni.sum()\n",
+    "        self.V = 4 * self.sum * self.ni\n",
+    "        self.N = self.ni.shape[0]\n",
+    "    def __call__(self, x: Any) -> float:\n",
+    "        f = x @ (Q @ x) + v.dot(x) + self.N**2\n",
+    "        return f\n",
+    "    \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x = np.arange(6)\n",
+    "\n",
+    "qfunc = qubo(x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\u001b[0;31mSignature:\u001b[0m       \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mouter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mCall signature:\u001b[0m  \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mouter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mType:\u001b[0m            _ArrayFunctionDispatcher\n",
+      "\u001b[0;31mString form:\u001b[0m     <function outer at 0x1194779c0>\n",
+      "\u001b[0;31mFile:\u001b[0m            ~/miniforge/envs/dwave/lib/python3.12/site-packages/numpy/core/numeric.py\n",
+      "\u001b[0;31mDocstring:\u001b[0m      \n",
+      "Compute the outer product of two vectors.\n",
+      "\n",
+      "Given two vectors `a` and `b` of length ``M`` and ``N``, repsectively,\n",
+      "the outer product [1]_ is::\n",
+      "\n",
+      "  [[a_0*b_0  a_0*b_1 ... a_0*b_{N-1} ]\n",
+      "   [a_1*b_0    .\n",
+      "   [ ...          .\n",
+      "   [a_{M-1}*b_0            a_{M-1}*b_{N-1} ]]\n",
+      "\n",
+      "Parameters\n",
+      "----------\n",
+      "a : (M,) array_like\n",
+      "    First input vector.  Input is flattened if\n",
+      "    not already 1-dimensional.\n",
+      "b : (N,) array_like\n",
+      "    Second input vector.  Input is flattened if\n",
+      "    not already 1-dimensional.\n",
+      "out : (M, N) ndarray, optional\n",
+      "    A location where the result is stored\n",
+      "\n",
+      "    .. versionadded:: 1.9.0\n",
+      "\n",
+      "Returns\n",
+      "-------\n",
+      "out : (M, N) ndarray\n",
+      "    ``out[i, j] = a[i] * b[j]``\n",
+      "\n",
+      "See also\n",
+      "--------\n",
+      "inner\n",
+      "einsum : ``einsum('i,j->ij', a.ravel(), b.ravel())`` is the equivalent.\n",
+      "ufunc.outer : A generalization to dimensions other than 1D and other\n",
+      "              operations. ``np.multiply.outer(a.ravel(), b.ravel())``\n",
+      "              is the equivalent.\n",
+      "tensordot : ``np.tensordot(a.ravel(), b.ravel(), axes=((), ()))``\n",
+      "            is the equivalent.\n",
+      "\n",
+      "References\n",
+      "----------\n",
+      ".. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*, 3rd\n",
+      "       ed., Baltimore, MD, Johns Hopkins University Press, 1996,\n",
+      "       pg. 8.\n",
+      "\n",
+      "Examples\n",
+      "--------\n",
+      "Make a (*very* coarse) grid for computing a Mandelbrot set:\n",
+      "\n",
+      ">>> rl = np.outer(np.ones((5,)), np.linspace(-2, 2, 5))\n",
+      ">>> rl\n",
+      "array([[-2., -1.,  0.,  1.,  2.],\n",
+      "       [-2., -1.,  0.,  1.,  2.],\n",
+      "       [-2., -1.,  0.,  1.,  2.],\n",
+      "       [-2., -1.,  0.,  1.,  2.],\n",
+      "       [-2., -1.,  0.,  1.,  2.]])\n",
+      ">>> im = np.outer(1j*np.linspace(2, -2, 5), np.ones((5,)))\n",
+      ">>> im\n",
+      "array([[0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j],\n",
+      "       [0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j],\n",
+      "       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
+      "       [0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j],\n",
+      "       [0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j]])\n",
+      ">>> grid = rl + im\n",
+      ">>> grid\n",
+      "array([[-2.+2.j, -1.+2.j,  0.+2.j,  1.+2.j,  2.+2.j],\n",
+      "       [-2.+1.j, -1.+1.j,  0.+1.j,  1.+1.j,  2.+1.j],\n",
+      "       [-2.+0.j, -1.+0.j,  0.+0.j,  1.+0.j,  2.+0.j],\n",
+      "       [-2.-1.j, -1.-1.j,  0.-1.j,  1.-1.j,  2.-1.j],\n",
+      "       [-2.-2.j, -1.-2.j,  0.-2.j,  1.-2.j,  2.-2.j]])\n",
+      "\n",
+      "An example using a \"vector\" of letters:\n",
+      "\n",
+      ">>> x = np.array(['a', 'b', 'c'], dtype=object)\n",
+      ">>> np.outer(x, [1, 2, 3])\n",
+      "array([['a', 'aa', 'aaa'],\n",
+      "       ['b', 'bb', 'bbb'],\n",
+      "       ['c', 'cc', 'ccc']], dtype=object)\n",
+      "\u001b[0;31mClass docstring:\u001b[0m\n",
+      "Class to wrap functions with checks for __array_function__ overrides.\n",
+      "\n",
+      "All arguments are required, and can only be passed by position.\n",
+      "\n",
+      "Parameters\n",
+      "----------\n",
+      "dispatcher : function or None\n",
+      "    The dispatcher function that returns a single sequence-like object\n",
+      "    of all arguments relevant.  It must have the same signature (except\n",
+      "    the default values) as the actual implementation.\n",
+      "    If ``None``, this is a ``like=`` dispatcher and the\n",
+      "    ``_ArrayFunctionDispatcher`` must be called with ``like`` as the\n",
+      "    first (additional and positional) argument.\n",
+      "implementation : function\n",
+      "    Function that implements the operation on NumPy arrays without\n",
+      "    overrides.  Arguments passed calling the ``_ArrayFunctionDispatcher``\n",
+      "    will be forwarded to this (and the ``dispatcher``) as if using\n",
+      "    ``*args, **kwargs``.\n",
+      "\n",
+      "Attributes\n",
+      "----------\n",
+      "_implementation : function\n",
+      "    The original implementation passed in."
+     ]
+    }
+   ],
+   "source": []
+  },
+  {
+   "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.12.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}