Binary Integer Linear Programming (BILP)

Given a vector \(c\) of size \(N\), a vector \(b\) of size \(m\) and a matrix \(S\) of size \(m \times N\), the problem consists in finding a binary vector \(x\) (i.e. vector composed of 0 and 1) of size \(N\) that maximizes the dot product \(c*x\), such as \(S * x = b\).

Solving this problem using the simulated quantum annealing method requires our BILP class and SQAQPU with \(N\) spins (one spin per binary value of \(x\)).

import numpy as np
from qat.opt import BILP
from qat.qpus import SQAQPU
from qat.core import Variable
from qat.core.spins import integer_to_spins

# Specify the problem; the solution is x = [0, 1, 1, 0]
c = np.array([0, 1, 1, 1])
S = np.array([[1, 0, 1, 1], [0, 1, 0, 1]])
b = np.array([1, 1])

# Impose constraints for the right encoding
B = 1
A = 10 * B
bilp_problem = BILP(c, S, b, A, B)

# Extract parameters for SQA
problem_parameters_dict = bilp_problem.get_sqa_initial_parameters()
temp_max = problem_parameters_dict["temp_max"]
temp_min = problem_parameters_dict["temp_min"]
gamma_max = problem_parameters_dict["gamma_max"]
gamma_min = problem_parameters_dict["gamma_min"]
nbshots = problem_parameters_dict['nbshots']
n_steps = problem_parameters_dict['n_steps']
n_trotters = problem_parameters_dict['n_trotters']

# Create a temperature and a gamma schedule
tmax = 1.0
t = Variable("t", float)
temp_t = temp_min * (t / tmax) + temp_max * (1 - t / tmax)
gamma_t = gamma_min * (t / tmax) + gamma_max * (1 - t / tmax)

# Create a job and send it to a QPU
problem_job = bilp_problem.to_job('sqa', gamma_t=gamma_t, tmax=tmax, nbshots=nbshots)
sqa_qpu = SQAQPU(temp_t=temp_t, n_steps=n_steps, n_trotters=n_trotters, seed=9999)
problem_result = sqa_qpu.submit(problem_job)

# Sample analysis
Q = bilp_problem.q_matrix
offset = bilp_problem.offset_q
energies = []
best_state = 0
for sample in problem_result: # raw_data is a list of Samples - one per shot
    x = np.array([int(b) for b in sample.state.bitstring])
    energy = -x.T @ Q @ x - offset
    energies.append(energy)
    if energy <= min(energies):
        best_state= sample.state.int

# Present best configuration
solution_configuration = integer_to_spins(best_state, Q.shape[0])
print("Solution configuration: \n" + str(solution_configuration) + "\n")
x = [1 if spin == 1 else 0 for spin in solution_configuration]
print("x is respectively: \n" + str(x) + "\n")

# Calculate c * x
print("Largest value of c * x:\n" + str(np.dot(c,x)))

Solution configuration: 
[ 1. -1. -1.  1.]

x is respectively: 
[1, 0, 0, 1]

Largest value of c * x:
1

This example is also detailed in this notebook on solving BILP problems with SQA.