Quadratic Unconstrained Binary Optimization (QUBO)

Many \(\mathcal{NP}\)-hard discrete optimization problems that naturally arise in application fields such as finance, energy, healthcare, and machine learning, can be mapped to quadratically unconstrained binary optimization (QUBO) problems (Kochenberger et al.1), or equivalently to Ising Hamiltonians (Lucas2). Examples of such problems include max-cut, graph colouring, partitioning, and maximum independent set. A QUBO model is a mixed-integer quadratic program (MIQP) where the decision variables are restricted to take binary values and there are no explicit constraints.

Since solving optimization problems typically requires significant computing resources, research in the development of novel computing architectures designed to solve QUBO problems has flourished in recent years (Johnson et al.3, Matsubara et al.4, Farhi et al.5). Such devices are fast, general-purpose, and power efficient. However, QUBO is not nearly as expressive as MIP, so expressing an optimization problem in QUBO form can have significant drawbacks. For example, constraints have to be modeled using penalties variables, which can greatly expand the solution space and significantly increase the difficulty of the problem. While the lack of expressiveness of QUBO will probably be a significant barrier for solving general optimization problems, there are some problems that map quite nicely to QUBO form. These problems may see a significant performance boost as special-purpose approaches improve over time.

Problem Specification

We are given a set \(I\) of \(n\) items, weights \(q_i \in \mathbb{R}\) for each single item \(i \in I\), and weights \(q_{ij} \in \mathbb{R}\) for each pair of distinct items \(i,j \in I,~ i \neq j\).

The objective is to find a subset \(S^* \subseteq I\) of items such that the sum of weights of all single items in \(S^*\) and all pairs of distinct items in \(S^*\) is minimal, i.e.,

\[S^* = \arg \min_{S \subseteq I} \sum_{i \in S} q_i + \sum_{i,j \in S,~ i \neq j} q_{ij}\]

We arrange the weights in an upper triangular matrix \(Q \in \mathbb{R}^{n \times n}\) where entry \(q_{ij}\) with \(i < j\) is the weight for item pair \(i,j\), and entry \(q_{ii} = q_i\) is the weight for single item \(i\).

Note that the input matrix does not necessarily need to be in upper triangular format. We accept matrices \(Q'\) that are populated in an arbitrary way and accumulate symmetric entries, i.e., \(q_{ij} = q'_{ij} + q'_{ji}\) for all item pairs \(i,j\) with \(i < j\).

Background: Optimization Model

This Mod is implemented by formulating the QUBO problem as a Binary Quadratic Program (BQP). To do so, we define a binary decision vector \(x \in \{0,1\}^n\) with variables

\[\begin{split}x_i = \begin{cases} 1 & \text{if item}\,i \in S^*\\ 0 & \text{otherwise.} \\ \end{cases}\end{split}\]

The BQP is then formulated as:

\[\begin{split}\begin{align} \min \quad & x' Q x = \sum_{i \in I} \sum_{j \in I} q_{ij} x_i x_j & \\ \mbox{s.t.} \quad & x \in \{0, 1\}^n & \end{align}\end{split}\]

Note that weights \(q_i = q_{ii}\) for single items \(i \in I\) are correctly considered in the objective function since \(x_i x_i = x_i\) holds for binary variables.

The input data consisting of the item (pair) weights is defined as a matrix (see the description), either as a NumPy array ndarray or as a SciPy sparse matrix spmatrix.

Code

The example below solves a QUBO problem instance based on 3 items with single-item weights \(q_1 = 0,~ q_2 = -3,~ q_3 = 2\), and item-pair weights \(q_{12} = -1,~ q_{13} = -2,~ q_{23} = 3\), resulting in the following item weight matrix:

\[\begin{split}Q = \begin{pmatrix} 0 & -1 & -2\\ 0 & -3 & 3\\ 0 & 0 & 2 \end{pmatrix}\end{split}\]

We use a NumPy array to represent matrix \(Q\) (and alternatively we show the definition as a SciPy sparse matrix in a comment).

import numpy as np
import scipy.sparse as sp
from gurobi_optimods.qubo import solve_qubo

Q = np.array([[0, -1, -2], [0, -3, 3], [0, 0, 2]])

# weights = [-3, 2, -1, -2, 3]
# row = [1, 2, 0, 0, 1]
# col = [1, 2, 1, 2, 2]
# Q = sp.coo_array((weights, (row, col)), shape=(3, 3))

result = solve_qubo(Q)

Solution

The returned result is a data class containing the objective value and the solution itself as a NumPy ndarray.

>>> result
QuboResult(solution=array([1., 1., 0.]), objective_value=-4.0)
>>> result.objective_value
-4.0
>>> result.solution
array([1., 1., 0.])

1

Gary Kochenberger, Jin-Kao Hao, Fred Glover, Mark Lewis, Zhipeng Lü, Haibo Wang, and Yang Wang. The unconstrained binary quadratic programming problem: a survey. Journal of Combinatorial Optimization, 28:58–81, 2014.

2

Andrew Lucas. Ising formulations of many np problems. Frontiers in physics, 2:5, 2014.

3

Mark W Johnson, Mohammad HS Amin, Suzanne Gildert, Trevor Lanting, Firas Hamze, Neil Dickson, Richard Harris, Andrew J Berkley, Jan Johansson, Paul Bunyk, and others. Quantum annealing with manufactured spins. Nature, 473(7346):194–198, 2011.

4

Satoshi Matsubara, Hirotaka Tamura, Motomu Takatsu, Danny Yoo, Behraz Vatankhahghadim, Hironobu Yamasaki, Toshiyuki Miyazawa, Sanroku Tsukamoto, Yasuhiro Watanabe, Kazuya Takemoto, and others. Ising-model optimizer with parallel-trial bit-sieve engine. In Complex, Intelligent, and Software Intensive Systems: Proceedings of the 11th International Conference on Complex, Intelligent, and Software Intensive Systems (CISIS-2017), 432–438. Springer, 2018.

5

Edward Farhi, Jeffrey Goldstone, and Sam Gutmann. A quantum approximate optimization algorithm. arXiv preprint arXiv:1411.4028, 2014.