Sampling from QUBO with OMMX Adapter

Here, we explain how to convert a problem to QUBO and perform sampling using the Traveling Salesman Problem as an example.

Fig. 5 Illustration of a man in a suit

The Traveling Salesman Problem (TSP) is about finding a route for a salesman to visit multiple cities in sequence. Given the travel costs between cities, we seek to find the path that minimizes the total cost. Let’s consider the following city arrangement:

# From ulysses16.tsp in TSPLIB
ulysses16_points = [
    (38.24, 20.42),
    (39.57, 26.15),
    (40.56, 25.32),
    (36.26, 23.12),
    (33.48, 10.54),
    (37.56, 12.19),
    (38.42, 13.11),
    (37.52, 20.44),
    (41.23, 9.10),
    (41.17, 13.05),
    (36.08, -5.21),
    (38.47, 15.13),
    (38.15, 15.35),
    (37.51, 15.17),
    (35.49, 14.32),
    (39.36, 19.56),
]

Let’s plot the locations of the cities.

%matplotlib inline
from matplotlib import pyplot as plt

x_coords, y_coords = zip(*ulysses16_points)
plt.scatter(x_coords, y_coords)
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Ulysses16 Points')
plt.show()

Let’s consider distance as the cost. We’ll calculate the distance \(d(i, j)\) between city \(i\) and city \(j\).

def distance(x, y):
    return ((x[0] - y[0])**2 + (x[1] - y[1])**2)**0.5

# Number of cities
N = len(ulysses16_points)
# Distance between each pair of cities
d = [[distance(ulysses16_points[i], ulysses16_points[j]) for i in range(N)] for j in range(N)]

Using this, we can formulate TSP as follows. First, let’s represent whether we are at city \(i\) at time \(t\) with a binary variable \(x_{t, i}\). Then, we seek \(x_{t, i}\) that satisfies the following constraints. The distance traveled by the salesman is given by:

\[ \sum_{t=0}^{N-1} \sum_{i, j = 0}^{N-1} d(i, j) x_{t, i} x_{(t+1 \% N), j} \]

However, \(x_{t, i}\) cannot be chosen freely and must satisfy two constraints: at each time \(t\), the salesman can only be in one city, and each city must be visited exactly once:

\[ \sum_{i=0}^{N-1} x_{t, i} = 1, \quad \sum_{t=0}^{N-1} x_{t, i} = 1 \]

Combining these, TSP can be formulated as a constrained optimization problem:

\[\begin{split} \begin{align*} \min \quad & \sum_{t=0}^{N-1} \sum_{i, j = 0}^{N-1} d(i, j) x_{t, i} x_{(t+1 \% N), j} \\ \text{s.t.} \quad & \sum_{i=0}^{N-1} x_{t, i} = 1 \quad (\forall t = 0, \ldots, N-1) \\ \quad & \sum_{t=0}^{N-1} x_{t, i} = 1 \quad (\forall i = 0, \ldots, N-1) \end{align*} \end{split}\]

The corresponding ommx.v1.Instance can be created as follows:

from ommx.v1 import DecisionVariable, Instance

x = [[
        DecisionVariable.binary(
            i + N * t,  # Decision variable ID
            name="x",           # Name of the decision variable, used when extracting solutions
            subscripts=[t, i])  # Subscripts of the decision variable, used when extracting solutions
        for i in range(N)
    ]
    for t in range(N)
]

objective = sum(
    d[i][j] * x[t][i] * x[(t+1) % N][j]
    for i in range(N)
    for j in range(N)
    for t in range(N)
)
place_constraint = [
    (sum(x[t][i] for i in range(N)) == 1)
        .set_id(t)  # type: ignore
        .add_name("place")
        .add_subscripts([t])
    for t in range(N)
]
time_constraint = [
    (sum(x[t][i] for t in range(N)) == 1)
        .set_id(i + N)  # type: ignore
        .add_name("time")
        .add_subscripts([i])
    for i in range(N)
]

instance = Instance.from_components(
    decision_variables=[x[t][i] for i in range(N) for t in range(N)],
    objective=objective,
    constraints=place_constraint + time_constraint,
    sense=Instance.MINIMIZE
)

The variable names and subscripts added to DecisionVariable.binary during creation will be used later when interpreting the obtained samples.

Sampling with OpenJij

To sample the QUBO described by ommx.v1.Instance using OpenJij, use the ommx-openjij-adapter.

from ommx_openjij_adapter import OMMXOpenJijSAAdapter

sample_set = OMMXOpenJijSAAdapter.sample(instance, num_reads=16, uniform_penalty_weight=20.0)
sample_set.summary

	objective	feasible
sample_id
15	79.406924	True
7	87.993033	True
3	91.841910	True
0	95.679199	True
10	96.609239	True
2	97.993181	True
12	101.666141	True
14	112.381833	True
13	92.036314	False
8	94.225558	False
9	97.529814	False
11	101.878419	False
5	102.115437	False
6	102.849302	False
1	111.974673	False
4	139.667391	False

OMMXOpenJijSAAdapter.sample returns ommx.v1.SampleSet, which stores the evaluated objective function values and constraint violations in addition to the decision variable values of samples. The SampleSet.summary property is used to display summary information. feasible indicates the feasibility to the original problem before conversion to QUBO. This is calculated using the information stored in removed_constraints of the qubo instance.

To view the feasibility for each constraint, use the summary_with_constraints property.

sample_set.summary_with_constraints

	objective	feasible	place[0]	place[1]	place[2]	place[3]	place[4]	place[5]	place[6]	place[7]	...	time[6]	time[7]	time[8]	time[9]	time[10]	time[11]	time[12]	time[13]	time[14]	time[15]
sample_id
15	79.406924	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
7	87.993033	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
3	91.841910	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
0	95.679199	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
10	96.609239	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
2	97.993181	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
12	101.666141	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
14	112.381833	True	True	True	True	True	True	True	True	True	...	True	True	True	True	True	True	True	True	True	True
13	92.036314	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
8	94.225558	False	True	True	True	True	True	True	True	False	...	True	True	True	True	False	True	True	True	True	True
9	97.529814	False	False	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
11	101.878419	False	True	False	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
5	102.115437	False	True	True	True	False	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
6	102.849302	False	True	True	True	True	True	True	True	False	...	True	True	True	True	False	True	True	True	True	True
1	111.974673	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True
4	139.667391	False	True	True	True	True	True	True	True	True	...	True	True	True	True	False	True	True	True	True	True

16 rows × 34 columns

For more detailed information, you can use the SampleSet.decision_variables and SampleSet.constraints properties.

sample_set.decision_variables.head(2)

	kind	lower	upper	name	subscripts	description	substituted_value	0	3	4	...	8	9	11	12	13	15	1	2	10	14
id
0	binary	0.0	1.0	x	[0, 0]	<NA>	<NA>	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	1.0	1.0	1.0	1.0
16	binary	0.0	1.0	x	[1, 0]	<NA>	<NA>	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

2 rows × 23 columns

sample_set.constraints.head(2)

	equality	used_ids	name	subscripts	description	removed_reason	value.0	value.1	value.2	value.3	...	feasible.14	feasible.6	feasible.15	feasible.12	feasible.8	feasible.0	feasible.5	feasible.13	feasible.1	feasible.4
id
0	=0	{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...	place	[0]	<NA>	uniform_penalty_method	0.0	0.0	0.0	0.0	...	True	True	True	True	True	True	True	True	True	True
1	=0	{16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 2...	place	[1]	<NA>	uniform_penalty_method	0.0	0.0	0.0	0.0	...	True	True	True	True	True	True	True	True	True	True

2 rows × 38 columns

To obtain the samples, use the SampleSet.extract_decision_variables method. This interprets the samples using the name and subscripts registered when creating ommx.v1.DecisionVariables. For example, to get the value of the decision variable named x with sample_id=1, use the following to obtain it in the form of dict[subscripts, value].

sample_id = 1
x = sample_set.extract_decision_variables("x", sample_id)
t = 2
i = 3
x[(t, i)]

1.0

Since we obtained a sample for \(x_{t, i}\), we convert this into a TSP path. This depends on the formulation used, so you need to write the processing yourself.

def sample_to_path(sample: dict[tuple[int, ...], float]) -> list[int]:
    path = []
    for t in range(N):
        for i in range(N):
            if sample[(t, i)] == 1:
                path.append(i)
    return path

Let’s display this. First, we obtain the IDs of samples that are feasible for the original problem.

feasible_ids = sample_set.summary.query("feasible == True").index
feasible_ids

Index([15, 7, 3, 0, 10, 2, 12, 14], dtype='int64', name='sample_id')

Let’s display the optimized paths for these samples.

fig, axie = plt.subplots(3, 3, figsize=(12, 12))

for i, ax in enumerate(axie.flatten()):
    if i >= len(feasible_ids):
        break
    s = feasible_ids[i]
    x = sample_set.extract_decision_variables("x", s)
    path = sample_to_path(x)
    xs = [ulysses16_points[i][0] for i in path] + [ulysses16_points[path[0]][0]]
    ys = [ulysses16_points[i][1] for i in path] + [ulysses16_points[path[0]][1]]
    ax.plot(xs, ys, marker='o')
    ax.set_title(f"Sample {s}, objective={sample_set.objectives[s]:.2f}")

plt.tight_layout()
plt.show()