I have been looking at QUBO (quadratically unconstrained binary optimization) problems, and wanted to try to make a (probably not very efficent) branch & bound iterativley, and use this as a space to voice ideas. The QUBO problem is essentually just quadratic programming were the optimization variables are binary, and there are not any other constraints. Below is the general definition of the problem that I am looking at. Without loss of generality, we can assume that $Q = Q^T \succ 0$. That being said, even with a convex objective, this is a known NP-Hard Problem [1].
Here, I don’t want to spend a lot of time on the describing the B&B (Branch & Bound) algorithm, but the idea is that when we try to solve problems like above, tryically there is a combinatorial number of subproblems, that would have to be considered to verify a solution if we just enumerate all possible combinations. However, by instead exploring how the problems are interrelated, and removing entire branches of subproblems from consideration based on pruning or bounding rules, we can significantly reduce the number of subproblems that we need to consider. Here is a good reference to B&B algorithms, if you are unfamiliar with the concept then I would suggest you read this first [2]. I want to appologies as this is going to be rather code heavy, subsiquent B&B postings will be much lighter, as they will not need to repeat all of this boilerplate.
To start of with, I want to write the skeleton of the QUPO problem object. Here, this is effectivley a data class that has some simple assessory functions.
class QUBOProblem:
def __init__(self, Q, c):
self.Q = Q
self.c = c
def num_x(self):
return Q.shape[0]
def eval_obj(self,x):
return 0.5*x.T@self.Q@x + self.c.T@x
def is_feasible(self, x):
zeros = numpy.zeros(self.num_x()).reshape(-1,1)
ones = zeros + 1
near_zero = numpy.isclose(x, zeros)
near_ones = numpy.isclose(x, ones)
either = numpy.logical_or(near_zero, near_ones)
return numpy.all(either)
With the problem definition worked out, we need make a B&B node object to store all of our subproblem information and allows for easy tracking of the B&B tree. The basic utilities that it provides are solve the subproblem, and branch which will generate the next set of subproblems.
class BBNode:
def __init__(self, problem, fixed: dict = None, parent:BBNode = None):
if fixed is None:
fixed = {}
self.fixed = fixed
self.children = None
self.obj = None
self.sol = None
self.problem = problem
self.parent = parent
def num_fixed(self):
return len(self.fixed)
def solve(self, sub_problem_solver):
# hand node to the sub problem solver
self.sol = sub_problem_solver(self)
self.obj = self.problem.eval_obj(self.sol)
def branch(self, branch_stratagy):
branch_choice = branch_stratagy(self)
# Generate the branch to the Right & Left
zero_branch = copy.deepcopy(self.fixed)
one_branch = copy.deepcopy(self.fixed)
# Add fixed value
zero_branch[branch_choice] = 0
one_branch[branch_choice] = 1
# generate the children nodes
zero_node = BBNode(self.problem, zero_branch, parent=self)
one_node = BBNode(self.problem, one_branch, parent=self)
self.children = [zero_node, one_node]
return zero_node, one_node
With that done, we can write an initial interface for writing the overall solver that will operate on the branch and bound tree. Here this keeps track of all of the subproblems, in node_list, and the bounds and solutions. The iteratee function does a single round of selecting a node, solving the subproblem, and then branching (if nessisary). The solve function, calls iterate until there are no longer any subproblems to solve.
class QUBO_solver:
def __init__(self, problem, sub_problem_solver, branch_strat, node_strat, initial_sol = None):
if initial_sol is None:
initial_sol = numpy.zeros(problem.num_x()).reshape(-1,1)
self.problem = problem
self.initial_sol = initial_sol
self.initial_bound = self.problem.eval_obj(self.initial_sol)
self.sub_problem_solver = sub_problem_solver
self.root = BBNode(self.problem, fixed=None)
self.node_list = [self.root]
self.checked_nodes = 0
self.current_upper_bound = self.initial_bound
self.best_feasible_sol = self.initial_sol
self.branch_strat = branch_strat
self.node_strat = node_strat
def iterate(self):
# if there is nothing to do, early return
if len(self.node_list) == 0:
return
# select a node
node = self.node_strat(self.node_list)
self.node_list.remove(node)
# if the value of the parent of the node we selected is above the current
# feasible solution then we know the solution of the current node is above
if node.parent is not None:
if node.parent.obj > self.current_upper_bound:
return
# iterate the number of nodes we solved
self.checked_nodes += 1
# solve the sub_problem
node.solve(self.sub_problem_solver)
# check if we should branch or cull
node_obj = node.obj
# if the obj is above then we do not have to branch
if node_obj <= self.current_upper_bound:
feasible = self.problem.is_feasible(node.sol)
if feasible:
self.current_upper_bound = node.obj
self.best_feasible_sol = node.sol
if node.num_fixed() < self.problem.num_x() and not feasible:
z, o = node.branch(self.branch_strat)
self.node_list.append(z)
self.node_list.append(o)
def solve(self):
# some book keeping
start_time = time.time()
iter_count = 0
# while we have subproblems run an iteration
while len(self.node_list) != 0:
# do a single B&B iteration
self.iterate()
# display to screen some information about progress
if iter_count % 100 == 0:
curr_time = time.time() - start_time
print(f'{curr_time:4.3}',iter_count, float(self.current_upper_bound), len(self.node_list), self.checked_nodes)
iter_count += 1
print(f'{curr_time:4.3}',iter_count, float(self.current_upper_bound), len(self.node_list), self.checked_nodes)
Now, with B&B we have a choice on how to relax out problem to make it easy to solve the subproblems, as we are already taking absolutly no agency in how the B&B tree progresses other then bounding, we could make a very simple relaxation of the subproblems that would probably make a good number of people cringe. The idea of most relaxations is to make a simpler problem that is solvable that provide decent bounds for the B&B algorithm, or some other variation [3]. However here, we are relaxing the constraints on our optimization variables from $0$ or $1$ to any real number. This is a bit of an agressive relaxation (that is not good to be clear), but it allows us to solve everything using a single matrix operation (so it is fine for now). In either case, we can (and will) make better subproblem solvers that will have different relaxations; but for now we are doing the most simple thing.
\[\begin{align} \min_{x\in\mathbb{B}^n} \frac{1}{2}x^TQx + c^Tx \rightarrow \min_{x\in\mathbb{R}^n} \frac{1}{2}x^TQx + c^Tx \end{align}\]This leads us to running our first implimentation of the solver, without any of the bells and wistles but a version 0.0.-1 of the algorithm.
N = 20
numpy.random.seed(1696979520)
Q = numpy.random.randn(N,N)/20
Q = Q.T@Q + numpy.eye(N)
c = numpy.random.randn(N).reshape(-1,1)
p = QUBOProblem(Q, c)
qubo_solver = QUBO_solver(p, solve_eq_qp, random_branch_strat,random_node_strat)
qubo_solver.solve()
For this seed on my computer, it took $10.4$ seconds to solve and we see the answer is $-4.5582…$ Which is quite good considering how bad the relaxations are, and that node selection and branch variables are completely random. From the original $2^{20}$ subproblems, only $4412$ subproblems were required to be solved, which is a reduction of $238x$. This can be massivley improved on with, tighter subproblems, better branch variable and node selection.
Resources
[1] - QUBO
[2] - B&B Algorithm Overview
[3] - B&C Algorithm Overview
Code Reference
def random_node_strat(node_list):
return numpy.random.choice(node_list)
def random_branch_strat(node):
branchable = [i for i in range(node.problem.num_x()) if i not in node.fixed]
return numpy.random.choice(branchable)
def solve_eq_qp(node):
"""Solve equality constrained QP problem (the simple way)"""
# reduce the wordy-ness
p = node.problem
# number of equalities
num_eq = len(node.fixed)
num_var = p.num_x()
# if this is the root node we can skip
if num_eq == 0:
return numpy.linalg.solve(p.Q, -p.c)
# if it is fully constrained
if num_eq == num_var:
sol = numpy.zeros(num_var)
for var_index, value in node.fixed.items():
sol[var_index] = value
return sol.reshape(-1,1)
eq_block = numpy.zeros(num_var*num_eq).reshape(num_eq,-1)
for row_index,(var_index, value) in enumerate(node.fixed.items()):
eq_block[row_index,var_index] = 1
eq_rhs = numpy.array(list(node.fixed.values())).reshape(-1,1)
zeros = numpy.zeros(num_eq*num_eq).reshape(num_eq, num_eq)
LHS = numpy.block([[p.Q, eq_block.T],[eq_block, zeros]])
RHS = numpy.block([[-p.c],[eq_rhs]])
return numpy.linalg.solve(LHS, RHS)[:num_var]