Set Partitioning on Dirac

Device: Dirac-1

Introduction

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. In more mathematical literature, the set partitioning problem is sometimes referred to as exact cover.

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} .

Importance

While the problem of set partitioning may first seem esoteric, the underlying concepts are actually quite natural and often arise in the real world, as do closely related problems. It also has a natural constraint structure that often appears in many problems (eg. one-hot constraints), enforcing that every member of each set is included once and only once. Unlike the two-way one-hot constraints seen in travelling salesperson or the quadratic assignment problem, the one-hot constraints here do not form a regular structure. However, they do overlap since each set will have one constraint for each element in the set. A natural generalization of set partitioning is set cover, where the goal is to find the minimum weight configuration where each element is included at least once, but may be duplicated. A review of problems related to set partitioning and their applications can be found here.

Applications

One perhaps surprising application to set partitioning is corporate tax structuring. The set partitioning constraints originate from the fact that subsidiaries can be combined for tax purposes, but each of subsidiary cannot be contained in more than one conglomerate. Since taxes must be paid on all of them, each has to be contained somewhere. Another application is airline crew scheduling. In this setting, the sets are different round trips that could be assigned to a type of crew member, and exactly one of which is needed for each leg of a flight. This model straightforwardly extends to similar applications such as bus driver scheduling. The formation of quarantine groups that minimize the spread of disease is another application, in this case the sets are the potential groupings, and the members of the sets are individuals being quarantined. It is clearly undesirable to place someone in multiple such groups, but every member has to go somewhere, hence the set partitioning constraint.

Implementation

In this tutorial, we take a slightly different approach, which requires some explanation. We formulate the set partitioning problem as a quadratic linearly constrained binary optimization problem (QLCBO), but we construct the constraints explicitly. Instead of directly using the built-in constraint features of the software, we are sending the problem to the solver as a quadratic unconstrained binary optimization problem. This tutorial explicitly shows what the constrained solver will do for you behind the scenes. 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 additional step is required for unconstrained solving that we have elected to use in this example. 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([[1, 0, 1, 0, 0, 1, 0, 0, 0],
       [0, 1, 1, 0, 0, 0, 0, 0, 0],
       [1, 0, 0, 1, 0, 0, 1, 0, 1],
       [0, 1, 0, 0, 1, 0, 0, 1, 0],
       [1, 0, 0, 0, 1, 0, 0, 1, 1],
       [0, 1, 0, 1, 0, 0, 0, 1, 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 [11]:

token = "token_here"
api_url = "https://api.qci-prod.com"
client = QciClient(api_token=token, url=api_url)

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 [12]:

file_id = client.upload_file(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 [15]:

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

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

In [16]:

response = client.process_job(job_body=body)

Out [ ]:

2024-04-29 16:29:09 - Dirac allocation balance = 0 s (unmetered)
2024-04-29 16:29:09 - Job submitted: job_id='662fbcc5e15a79bd9d02c4a7'
2024-04-29 16:29:10 - QUEUED
2024-04-29 16:29:12 - RUNNING
2024-04-29 16:30:48 - COMPLETED
2024-04-29 16:30:51 - 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 [21]:

response

Out [21]:

{'job_info': {'job_id': '662fbcc5e15a79bd9d02c4a7',
  'job_submission': {'problem_config': {'quadratic_unconstrained_binary_optimization': {'qubo_file_id': '662fbc9198263204a365428a'}},
   'device_config': {'dirac-1': {'num_samples': 5}}},
  'job_status': {'submitted_at_rfc3339nano': '2024-04-29T15:29:09.875Z',
   'queued_at_rfc3339nano': '2024-04-29T15:29:09.88Z',
   'running_at_rfc3339nano': '2024-04-29T15:29:10.369Z',
   'completed_at_rfc3339nano': '2024-04-29T15:30:47.916Z'},
  'job_result': {'file_id': '662fbd2798263204a3654295', 'device_usage_s': 9}},
 'status': 'COMPLETED',
 'results': {'counts': [3, 2],
  'energies': [-6, -6],
  'solutions': [[0, 0, 1, 0, 0, 0, 1, 1, 0], [0, 1, 0, 0, 0, 1, 0, 0, 1]]}}

In [24]:

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 = response["results"]
    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 [25]:

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_2 S_6 S_7 
Set coverage 100.0%
Partition found: True
QUBO value: -6
Times sample found 2
Selected Sets
S_1 S_5 S_8 
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 [26]:

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

Out [26]:

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 [27]:

alpha = np.max(objective)
alpha

Out [27]:

72.09663468312299

In [28]:

Q = objective + alpha * P
Q

Out [28]:

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 [29]:

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

In [31]:

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

In [35]:

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

In [36]:

response = client.process_job(job_body=body)

Out [ ]:

2024-04-29 16:36:42 - Dirac allocation balance = 0 s (unmetered)
2024-04-29 16:36:43 - Job submitted: job_id='662fbe8be15a79bd9d02c4a9'
2024-04-29 16:36:44 - RUNNING
2024-04-29 16:38:21 - COMPLETED
2024-04-29 16:38:24 - Dirac allocation balance = 0 s (unmetered)

In [37]:

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

Out [ ]:

QUBO value: -398.7963365771029
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 [38]:

check_result(get_results(penalty_only_response), objective)

Out [ ]:

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

Conclusion

In this tutorial, we have demonstrated how to solve the set partitioning problem on our Dirac-1 device. Unlike other tutorials, we have explicitly demonstrated how to construct the constraints in a quadratic setting. A useful next step is to continue learning about quadratic linearly constrained binary optimization and to attempt some of the tutorials linked there, which use software to construct the constraints automatically. An alternative path is to look into quadratic unconstrained binary optimization, where the implementation of constraints is not necessary. Of course, another option is to start using our device to solve some of your own optimization problems.