Skip to content

Instantly share code, notes, and snippets.

@sndyuk

sndyuk/qubo_tutorial.ipynb

Last active February 16, 2023 18:15
Show Gist options
  • Save sndyuk/40cd214b44d14b264ae1798ff0e64964 to your computer and use it in GitHub Desktop.
Save sndyuk/40cd214b44d14b264ae1798ff0e64964 to your computer and use it in GitHub Desktop.
Download ZIP
QUBO Tutorial - Quadratic unconstrained binary optimization
Display the source blob
Display the rendered blob
Raw
qubo_tutorial.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "QUBO_Tutorial.ipynb",
"version": "0.3.2",
"provenance": [],
"collapsed_sections": [],
"include_colab_link": true
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/sndyuk/40cd214b44d14b264ae1798ff0e64964/qubo_tutorial.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"metadata": {
"id": "UyJ6ISvF9gwn",
"colab_type": "code",
"colab": {}
},
"source": [
"import numpy as np"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "-avWc61w9gwr",
"colab_type": "text"
},
"source": [
"# 1. QUBO expression\n",
"\n",
"- is written as a 2 dimentional matrix and binary variables.\n",
"- needs to be upper triangular matrix because it's an one of expression of quadratic form.\n",
"\n",
"Sample quadratic form with 3 variables(x, y, z) and the coeficients(10 for (x, x) => x^2, 20 for y^2, 30 for z^2, 2 for (x, y), 1 for (x, z)) .\n",
"\n",
"\\begin{equation}\n",
"A = 10x^2 + 20y^2 + 30z^2 + \\frac{2+2}{2}xy + \\frac{1+1}{2}xz\n",
"\\end{equation}\n",
"\n",
"Let's convert this to matrix.\n",
"\n",
"\\begin{equation}\n",
" A = \\left(\n",
" \\begin{array}{ccc}\n",
" 10 & 2 & 1 \\\\\n",
" 2 & 20 & 0 \\\\\n",
" 1 & 0 & 30\n",
" \\end{array}\n",
" \\right)\n",
" \\left(\n",
" \\begin{array}{ccc}\n",
" x \\\\\n",
" y \\\\\n",
" z\n",
" \\end{array}\n",
" \\right)\n",
"\\end{equation}\n",
"\n",
"Then fold lower off-diagonal to upper side to make it simpler since this is symmetric matrix.\n",
"\n",
"\\begin{equation}\n",
" A = \\left(\n",
" \\begin{array}{ccc}\n",
" 10 & 4 & 2 \\\\\n",
" 0 & 20 & 0 \\\\\n",
" 0 & 0 & 30\n",
" \\end{array}\n",
" \\right)\n",
" \\left(\n",
" \\begin{array}{ccc}\n",
" x \\\\\n",
" y \\\\\n",
" z\n",
" \\end{array}\n",
" \\right)\n",
"\\end{equation}\n"
]
},
{
"cell_type": "code",
"metadata": {
"id": "Lgvsag8Q9gws",
"colab_type": "code",
"outputId": "012459cc-a680-4ed2-b67a-4f1f3479829b",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
}
},
"source": [
"## Quadratic form matrix\n",
"Q = [\n",
" [10, 4, 2],\n",
" [0, 20, 0],\n",
" [0, 0, 30]\n",
"]\n",
"\n",
"# QUBO consists of the Quadratic form and its variables(= solution) which are all binary.\n",
"# - The 3x3 QUBO has 3 bits solution space(X).\n",
"# 4x4 QUBO has 4 bits.\n",
"# 5x5 QUBO has 5 bits and so on.\n",
"\n",
"# The solution is written as vector(or 1(col)xN(row) dimentional matrix).\n",
"# x: 0 or 1\n",
"X = [1, 1, 0]\n",
"\n",
"\n",
"# - The energy of the solution is evaluated by the logic.\n",
"# It's a summation of each rows and columns.\n",
"energy = .0\n",
"SOLUTION = np.array(X)\n",
"for i, c in enumerate(SOLUTION):\n",
" for j, d in enumerate(SOLUTION[i:]):\n",
" if c and d:\n",
" energy += Q[i][j + i]\n",
"print(energy)"
],
"execution_count": 5,
"outputs": [
{
"output_type": "stream",
"text": [
"34.0\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "iO_eZ6uG9gww",
"colab_type": "code",
"outputId": "2620a1e0-4529-4540-c378-6e3332e42330",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
}
},
"source": [
"# Also, it can be written as the mathematical formula.\n",
"# energy = \\inv{X} * Q * X\n",
"\n",
"_X = np.matrix(X).T # To 1 col, n row matrix.\n",
"_Xt = _X.T\n",
"_Q = np.matrix(Q)\n",
"\n",
"energy = _Xt * _Q * _X\n",
"# OR: (_Xt.dot(_Q)*_X).sum(axis=1)\n",
"\n",
"print(float(energy))"
],
"execution_count": 6,
"outputs": [
{
"output_type": "stream",
"text": [
"34.0\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "faKMpNcj9gw0",
"colab_type": "text"
},
"source": [
"# 2. Graph to QUBO"
]
},
{
"cell_type": "code",
"metadata": {
"id": "erzBquST9gw0",
"colab_type": "code",
"colab": {}
},
"source": [
"# An element of the QUBO can be a node of the Graph.\n",
"\n",
"# In this example, there are 3 nodes:\n",
"# node a have weight 10.\n",
"# node b have weight 0.\n",
"# node c have weight 5.\n",
"# Every weight of the nodes are written as diagonal matrix:\n",
"Q = [\n",
" [10, 0, 0],\n",
" [0, 0, 0],\n",
" [0, 0, 5]\n",
"]\n",
"\n",
"# There are some relations between each nodes:\n",
"# a-b: 400\n",
"# a-c: 3\n",
"# b-c: 200\n",
"# Every relation between each node are written as off-diagonal matrix:\n",
"Q = [\n",
" [10, 400, 3],\n",
" [0, 0, 200],\n",
" [0, 0, 5]\n",
"]\n",
"\n",
"# Important observation:\n",
"# If there is the relation below:\n",
"# c-a: 500 which means the edge has the specific direction(c->a).\n",
"# The matrix would be this but QUBO can't evaluate this kind of matrix because QUBO must be upper triangler matrix:\n",
"Q = [\n",
" [10, 400, 3],\n",
" [0, 0, 200],\n",
" [500, 0, 5]\n",
"]\n",
"# It needs to be this:\n",
"# a-c: +500 which means the edge doesn't have the specific direction.\n",
"Q = [\n",
" [10, 400, 3 + 500],\n",
" [0, 0, 200],\n",
" [0, 0, 5]\n",
"]\n",
"# we need to expand the solution space to handle the backward direction. E.g,\n",
"Q2 = [\n",
" # Forward direction(a->c)\n",
" [10, 400, 3, 0, 0, 0],\n",
" [0, 0, 200, 0, 0, 0],\n",
" [0, 0, 5, 0, 0, 0],\n",
" # Backward direction(c->a)\n",
" [0, 0, 0, 0, 0, 500],\n",
" [0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0],\n",
"]\n",
"\n",
"# The weight and the relation can be called 'cost'. And the summation of the cost can be called 'energy'."
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "em2Bpxjw9gw3",
"colab_type": "text"
},
"source": [
"# 3. QUBO expression part.2\n",
"\n",
"\\begin{equation}\n",
"\\sum_{ij} (\\sum_{k=1}^{9} x_{ijk} - 1)^2\n",
"\\end{equation}\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "KLkSkqUWO0D3",
"colab_type": "text"
},
"source": [
"### The formula to matrix\n",
"\n",
"\\begin{equation}\n",
"\\sum_{ij} (\\sum_{k=1}^{9} x_{ijk} - 1)^2\\\\\n",
"= \\sum_{ij} ((\\sum x_{ijk})^2 + (2 \\times (\\sum x_{ijk}) \\times -1) + (-1)^2)\\\\\n",
"\\end{equation}\n",
"\n",
"For example, calculate the first part of the formula$$(\\sum x_{ijk})^2$$ using the sample quadratic form.\n",
"\n",
"\n",
"Sample quadratic form:\n",
"\\begin{equation}\n",
" A = \\left(\n",
" \\begin{array}{ccc}\n",
" 10 & 0 & 0 \\\\\n",
" 0 & 20 & 0 \\\\\n",
" 0 & 0 & 30\n",
" \\end{array}\n",
" \\right)\n",
"\\end{equation}\n",
"\n",
"This can be\n",
"\n",
"\\begin{equation}\n",
"= (10x_{0} + 20x_{1} + 30x_{2})^2\\\\\n",
"= (10x_{0} + 20x_{1} + 30x_{2})(10x_{0} + 20x_{1} + 30x_{2})\\\\\n",
"= 100x_{0}^2 + 200x_{0}x_{1} + 300x_{0}x_{2} + ... + 900x_{2}^2\\\\\n",
"\\end{equation}\n",
"\n",
"when x is binary, x^2 and x are equivalent. \n",
"\n",
"\\begin{equation}\n",
"= 100x_{0} + 200x_{0}x_{1} + 300x_{0}x_{2} + ... + 900x_{2}\\\\\n",
"\\end{equation}\n",
"\n",
"To matrix:\n",
"\n",
"\\begin{equation}\n",
" A = \\left(\n",
" \\begin{array}{ccc}\n",
" 100 & 200 & 300 \\\\\n",
" 200 & 400 & 600 \\\\\n",
" 300 & 600 & 900\n",
" \\end{array}\n",
" \\right)\n",
"\\end{equation}\n",
"\n",
"To upper trianguler matrix:\n",
"\n",
"\\begin{equation}\n",
" A = \\left(\n",
" \\begin{array}{ccc}\n",
" 100 & 400 & 600 \\\\\n",
" 0 & 400 & 1200 \\\\\n",
" 0 & 0 & 900\n",
" \\end{array}\n",
" \\right)\n",
"\\end{equation}\n",
"\n",
"### To python"
]
},
{
"cell_type": "code",
"metadata": {
"id": "9mxe3yN-9gw4",
"colab_type": "code",
"colab": {}
},
"source": [
"# This part explain how we build QUBO from the math formula. Using SUDOKU example here: https://www.topcoder.com/challenges/30081256\n",
"def get_variable_id(N, i, j, k):\n",
" return (i * N + j) * N + k\n",
"\n",
"N = 9\n",
"solution_space = N * N * N"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "WA_szQg59gw7",
"colab_type": "code",
"outputId": "6af6603b-cf8e-470d-b09b-d9592e109bd6",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 153
}
},
"source": [
"# Example code 1. It's slow but descriptive.\n",
"Q = np.zeros((N * N * N, N * N * N))\n",
"constant = 0\n",
"\n",
"# \\sum_{ij}\n",
"for i in range(N):\n",
" for j in range(N):\n",
" # \\sum x_{ijk} - 1\n",
" Q_tmp = np.zeros((N * N * N, N * N * N))\n",
" constant_tmp = 0\n",
" for k in range(N):\n",
" ijk = get_variable_id(N, i, j, k)\n",
" Q_tmp[ijk][ijk] += 1\n",
" constant_tmp += -1\n",
"\n",
" # (\\sum x_{ijk} - 1)^2\n",
" # = (\\sum x_{ijk})^2 + (2 * x_{ijk} * -1) + -1^2\n",
" # part.1 part.2 part.3\n",
"\n",
" # part.1\n",
" diagonal = np.diag(Q_tmp)\n",
" A = np.outer(diagonal, diagonal)\n",
" A = (np.triu(A, k=1) + np.triu(A.T)) # Fold lower triangular to upper side.\n",
" \n",
" # part.2\n",
" # 2 * x_{ijk} * -1\n",
" B = (2 * np.diag(diagonal) * constant_tmp)\n",
" Q_tmp = A + B\n",
"\n",
" # part.3\n",
" # -1^2\n",
" constant_tmp = constant_tmp ** 2\n",
" \n",
" Q += Q_tmp\n",
" constant += constant_tmp\n",
"\n",
"print(Q)\n",
"print('constant:', constant)"
],
"execution_count": 9,
"outputs": [
{
"output_type": "stream",
"text": [
"[[-1. 2. 2. ... 0. 0. 0.]\n",
" [ 0. -1. 2. ... 0. 0. 0.]\n",
" [ 0. 0. -1. ... 0. 0. 0.]\n",
" ...\n",
" [ 0. 0. 0. ... -1. 2. 2.]\n",
" [ 0. 0. 0. ... 0. -1. 2.]\n",
" [ 0. 0. 0. ... 0. 0. -1.]]\n",
"constant: 81\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "tNc51KU-9gw_",
"colab_type": "code",
"outputId": "68c2a0fe-6513-4c3f-e842-f17be90e777c",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 153
}
},
"source": [
"# Example code 2. It's fast.\n",
"\n",
"Q = np.zeros((N * N * N, N * N * N))\n",
"constant = 0\n",
"\n",
"# \\sum_{ij}\n",
"for i in range(N):\n",
" for j in range(N):\n",
" # (\\sum x_{ijk} - 1)^2 for each ij\n",
" # = (\\sum x_{ijk})^2 + (2 * x_{ijk} * -1) + -1^2\n",
" # part.1 part.2 part.3\n",
"\n",
" for k1 in range(N):\n",
" # part.1\n",
" ijk1 = get_variable_id(N, i, j, k1)\n",
" for k2 in range(N):\n",
" ijk2 = get_variable_id(N, i, j, k2)\n",
" if ijk1 <= ijk2:\n",
" Q[ijk1][ijk2] += 1\n",
" else:\n",
" # Keep it as upper triangular matrix.\n",
" Q[ijk2][ijk1] += 1\n",
"\n",
" # part.2\n",
" # 2 * x_{ijk} * -1\n",
" Q[ijk1][ijk1] += -2\n",
"\n",
" # part.3\n",
" # -1^2\n",
" # constant only affects the total energy. It's not a part of QUBO matrix.\n",
" constant += 1\n",
"\n",
"print(Q)\n",
"print('constant:', constant)"
],
"execution_count": 10,
"outputs": [
{
"output_type": "stream",
"text": [
"[[-1. 2. 2. ... 0. 0. 0.]\n",
" [ 0. -1. 2. ... 0. 0. 0.]\n",
" [ 0. 0. -1. ... 0. 0. 0.]\n",
" ...\n",
" [ 0. 0. 0. ... -1. 2. 2.]\n",
" [ 0. 0. 0. ... 0. -1. 2.]\n",
" [ 0. 0. 0. ... 0. 0. -1.]]\n",
"constant: 81\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "Kj4-1mli9gxC",
"colab_type": "code",
"colab": {}
},
"source": [
""
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "xpEEdqA29gxF",
"colab_type": "code",
"colab": {}
},
"source": [
""
],
"execution_count": 0,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment