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 (Lucas[2]). 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.])