Set Partitioning on Dirac

Device: Dirac-1

The set partitioning problem is an optimization problem which selects sets $S_i$ $S_{i}$ from a collection $S$ $S$ such that each member $x\in X=\bigcup_i S_i$ $x \in X = ⋃_{i} S_{i}$ is included in exactly one $S_i$ $S_{i}$ of the selected sets (see Operations Research Journal Vol. 17, No. 5). It has applications in logistics, design, manufacturing and scheduling, among others.

Formulation

The set partitioning problem is formulated by creating constraints which specify that only one out of all the sets $S_i$ $S_{i}$ of which a particular member $s$ $s$ is selected. This looks like

\sum_{i,x\in S_i} s_i = 1\quad \forall x\in X

\sum_{i, x \in S_{i}} s_{i} = 1 \forall x \in X

where $s_i\in\{0,1\}$ $s_{i} \in {0, 1}$ indicates set $S_i$ $S_{i}$ is selected. In addition to a constraint for each member, there is an objective function which measures the cost, weight or benefit of selecting certain sets from $S$ $S$ . The objective function could be any form, but we can solve linear or quadratic objective functions with Dirac. A quadratic objective function has coefficients $c_{ij},c_{i}$ $c_{ij}, c_{i}$ for quadratic and linear terms, respectively. With the variables $s_i$ $s_{i}$ , we have

\sum_j\sum_i c_{ij}s_is_j + \sum_i c_is_i.

\sum_{j} \sum_{i} c_{ij} s_{i} s_{j} + \sum_{i} c_{i} s_{i} .

Implementation

With constraint-based solvers, it would be sufficient to implement the objective and constraints directly in the modeling or matrix format required by the solver, but an addition step is required for unconstrained solving. Penalties can be created from the constraints $As=b$ $A s = b$ using

P(s)=s^T(A^TA)s-(2b^TA)^Ts+b^Tb.

P (s) = s^{T} (A^{T} A) s - (2 b^{T} A)^{T} s + b^{T} b .

When all constraints are met, $P(s)=0$ $P (s) = 0$ and when any constraint is violated, $P(s)>0$ $P (s) > 0$ . There is a difficulty in combining the objective function and penalties in that a scalar value of the objective function for a constraint violating solution could be less than 0 or at least less than the penalty, which results in a value of the total function less than 0. This will result in an optimizer finding an infeasible solution, unless a multiplier is applied to $P(s)$ $P (s)$ . Sufficiently large multipliers will guarantee that no infeasible solution will take on a value of the total function which is less than any value for a feasible solution.

Here we have a value which is known to be sufficient, $\alpha$ $α$ . We apply it to $P(s)$ $P (s)$ to get the function

Q(s)=\sum_j\sum_i c_{ij}s_is_j + \sum_i c_is_i + \alpha P(s)

Q (s) = \sum_{j} \sum_{i} c_{ij} s_{i} s_{j} + \sum_{i} c_{i} s_{i} + α P (s)

Demonstration

In [1]:

import numpy as np
from qci_client import QciClient

Creating the Data

$S$ $S$ is a collection of 9 different sets. The members of the sets are the letters A through F. $W$ $W$ are the weights of each subset $S_i$ $S_{i}$ .

In [2]:

np.random.seed(21)
S = [{"A", "B", "C"}, {"D", "E", "F"}, {"A", "F"}, {"B", "E"}, {"C", "D"}, {"A"},
{"B"}, {"C", "D", "E"}, {"B", "C"}]
W = [100 * np.random.random() for S_i in S]

$X$ $X$ is the union of all $S_i$ $S_{i}$ .

In [3]:

X = set()
for S_i in S:
X = X.union(S_i)
X

Out [3]:

{'A', 'B', 'C', 'D', 'E', 'F'}

Create the constraints by indicating if a member is in a subset with a 1 in the position for the variable $s_i$ $s_{i}$ for every $x$ $x$ .

In [4]:

A = []
for x in X:
row = [1 if x in S_i else 0 for S_i in S]
A.append(row)
A = np.array(A)
A

Out [4]:

array([[0, 1, 1, 0, 0, 0, 0, 0, 0],
       [1, 0, 0, 0, 1, 0, 0, 1, 1],
       [0, 1, 0, 0, 1, 0, 0, 1, 0],
       [1, 0, 0, 1, 0, 0, 1, 0, 1],
       [0, 1, 0, 1, 0, 0, 0, 1, 0],
       [1, 0, 1, 0, 0, 1, 0, 0, 0]])

In [5]:

b = np.ones((A.shape[0],))
b

Out [5]:

array([1., 1., 1., 1., 1., 1.])

In [6]:

J = A.T@A
h = -2 * b.T@A
n = h.shape[0]
h = h.reshape((n, 1))
offset = b.T@b
J, h, offset

Out [6]:

(array([[3, 0, 1, 1, 1, 1, 1, 1, 2],
        [0, 3, 1, 1, 1, 0, 0, 2, 0],
        [1, 1, 2, 0, 0, 1, 0, 0, 0],
        [1, 1, 0, 2, 0, 0, 1, 1, 1],
        [1, 1, 0, 0, 2, 0, 0, 2, 1],
        [1, 0, 1, 0, 0, 1, 0, 0, 0],
        [1, 0, 0, 1, 0, 0, 1, 0, 1],
        [1, 2, 0, 1, 2, 0, 0, 3, 1],
        [2, 0, 0, 1, 1, 0, 1, 1, 2]]),
 array([[-6.],
        [-6.],
        [-4.],
        [-4.],
        [-4.],
        [-2.],
        [-2.],
        [-6.],
        [-4.]]),
 6.0)

Solving the problem without an objective function should reveal if there are multiple solutions to the exact cover problem.

First, create a connection to the REST API using QciClient.

In [7]:

client = QciClient()

The next line creates a QUBO from the quadratic and linear portions of the penalty function by adding all the linear terms in the diagonal of the quadratic operator. This file is uploaded to the API.

In [8]:

P = J + np.diag(h.T[0])
P

Out [8]:

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

In [9]:

qubo_file = {"file_name": "penalty-only-sp-qubo", "file_config": {"qubo": {"data":P}}}

In [10]:

file_id = client.upload_file(qubo_file)["file_id"]

Using the file ID returned by the upload method, build the job body requesting the job to run on eqc1 (aka Dirac-1), returning 5 samples.

In [11]:

body = client.build_job_body("sample-qubo", qubo_file_id=file_id, job_params={
"sampler_type": "dirac-1", "nsamples": 5
})

The job type sample-qubo converts the QUBO into an Ising Hamiltonian before sending to Dirac-1 for sampling.

In [12]:

response = client.process_job(job_type="sample-qubo", job_body=body)

Out [ ]:

Dirac allocation balance = 0 s (unmetered)
Job submitted job_id='65f261379e39850da64a40fb'-: 2024/03/13 20:30:15
RUNNING: 2024/03/13 20:30:17
COMPLETED: 2024/03/13 20:31:22
Dirac allocation balance = 0 s (unmetered)

The iteration over the samples tests if the sample selected all the members of $X$ $X$ exactly once. 100% set coverage indicates that all members were included and a partition was found if no member is included in multiple sets.

In [13]:

def check_result(results, objective=None):
samples = results["solutions"]
energies = results["energies"]
counts = results["counts"]
for j, sample in enumerate(samples):
print("QUBO value:", energies[j])
print("Times sample found", counts[j])
sample = np.array(sample)
if objective is not None:
print("Objective Function value:", sample.T@objective@sample)
print("Selected Sets")
testX = set()
members = []
for i in range(len(S)):
if sample[i] == 1:
print(f"S_{i}", end=" ")
testX = testX.union(S[i])
members.extend(S[i])
print()
print(f"Set coverage {100*len(testX)/len(X)}%")
print(f"Partition found: {len(testX)==len(members)}")
def get_results(response):
if "results" in response and response["results"] is not None:
results_file_config = response["results"]["file_config"]
# results file config has one key, named by the job type detail
assert len(results_file_config) == 1, "Unknown results format"
results = list(results_file_config.values())[0]
else:
if "job_info" in response and "job_result" in response["job_info"]:
details = response["job_info"]["job_result"]
else:
details = None
raise RuntimeError(f"Execution failed. See details: {details}")
return results

In [14]:

results = get_results(response)
check_result(results)
# save the response
penalty_only_response = response

Out [ ]:

QUBO value: -6
Times sample found 3
Selected Sets
S_1 S_5 S_8 
Set coverage 100.0%
Partition found: True
QUBO value: -6
Times sample found 2
Selected Sets
S_2 S_6 S_7 
Set coverage 100.0%
Partition found: True

The objective function specified is a linear function. It is built into the diagonal of a matrix, as the linear portion of the penalty function.

In [15]:

objective = np.diag([W[i] for i in range(len(S))])
objective

Out [15]:

array([[ 4.87248808,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        , 28.91096598,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        , 72.09663468,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  2.16162499,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        , 20.59227653,
         0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         5.07732567,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        , 30.2271894 ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        , 66.39102946,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        , 30.81143932]])

In [16]:

alpha = np.max(objective)
alpha

Out [16]:

72.09663468312299

In [17]:

Q = objective + alpha * P
Q

Out [17]:

array([[-211.41741597,    0.        ,   72.09663468,   72.09663468,
          72.09663468,   72.09663468,   72.09663468,   72.09663468,
         144.19326937],
       [   0.        , -187.37893807,   72.09663468,   72.09663468,
          72.09663468,    0.        ,    0.        ,  144.19326937,
           0.        ],
       [  72.09663468,   72.09663468,  -72.09663468,    0.        ,
           0.        ,   72.09663468,    0.        ,    0.        ,
           0.        ],
       [  72.09663468,   72.09663468,    0.        , -142.03164437,
           0.        ,    0.        ,   72.09663468,   72.09663468,
          72.09663468],
       [  72.09663468,   72.09663468,    0.        ,    0.        ,
        -123.60099284,    0.        ,    0.        ,  144.19326937,
          72.09663468],
       [  72.09663468,    0.        ,   72.09663468,    0.        ,
           0.        ,  -67.01930901,    0.        ,    0.        ,
           0.        ],
       [  72.09663468,    0.        ,    0.        ,   72.09663468,
           0.        ,    0.        ,  -41.86944529,    0.        ,
          72.09663468],
       [  72.09663468,  144.19326937,    0.        ,   72.09663468,
         144.19326937,    0.        ,    0.        , -149.89887459,
          72.09663468],
       [ 144.19326937,    0.        ,    0.        ,   72.09663468,
          72.09663468,    0.        ,   72.09663468,   72.09663468,
        -113.38183004]])

In [18]:

qubo_file = {"file_name": "full-sp-qubo", "file_config": {"qubo": {"data": Q}}}

In [19]:

file_id = client.upload_file(qubo_file)["file_id"]

In [20]:

body = client.build_job_body("sample-qubo", qubo_file_id=file_id, job_params={
"sampler_type": "dirac-1", "nsamples": 5
})

In [21]:

response = client.process_job(job_type="sample-qubo", job_body=body)

Out [ ]:

Dirac allocation balance = 0 s (unmetered)
Job submitted job_id='65f2617c9e39850da64a40fd'-: 2024/03/13 20:31:24
RUNNING: 2024/03/13 20:31:26
COMPLETED: 2024/03/13 20:32:31
Dirac allocation balance = 0 s (unmetered)

In [22]:

result = get_results(response)
check_result(result, objective)

Out [ ]:

QUBO value: -398.79632568359375
Times sample found 5
Objective Function value: 33.78345405989442
Selected Sets
S_0 S_1 
Set coverage 100.0%
Partition found: True

Check the samples from the saved response with the objective function included to manually validate the minimization. Note: The QUBO value does not change because it represents the penalty-only QUBO and not including the objective function.

In [23]:

check_result(get_results(penalty_only_response), objective)

Out [ ]:

QUBO value: -6
Times sample found 3
Objective Function value: 64.79973097220724
Selected Sets
S_1 S_5 S_8 
Set coverage 100.0%
Partition found: True
QUBO value: -6
Times sample found 2
Objective Function value: 168.71485354205467
Selected Sets
S_2 S_6 S_7 
Set coverage 100.0%
Partition found: True