# Formulation Benchmarks#

We have benchmarked the formulation of mathematical optimization models provided as a library in Python against the Amplify SDK. We measured the execution time to create a model and output it as QUBO using the traveling salesperson problem formulation as an example, assuming using the QUBO solver.

However, each library, including the Amplify SDK, covers different features depending on the formulation method. Here is a summary of the features and policies of each library, formulated as follows.

-

Amplify

PyQUBO

dimod
BQM (index)

dimod
BQM (symbol)

dimod
CQM*

Symbolic operation

β

β

β

β

β

Objective function

β

β

β

β

β

Constraint

β

β**

β***

β***

β*

Higher order polynomial

β

β

β

β

β

Coefficient matrix

β

β

β

β

β

Variable type

B/S/I/R

B/S

B/S

B/S

B/S/I/R

Supported machines

Various

Depends on user

D-Wave only

D-Wave only

D-Wave only

B: Binary, S: Ising Spin, I: Integer, R: Real

*: Model creation only, as QUBO output (conversion of constraints to the penalty functions) is not available.
**: Must define penalty function
***: Constraints are expressed by adding penalty functions to the objective function

Amplify
import amplify

def tsp_for_amplify(ncity: int, distances: np.ndarray, dmax: float):
q = amplify.VariableGenerator().array("Binary", ncity + 1, ncity)
q[-1, :] = q[0, :]

# Objective function
objective: amplify.Poly = amplify.einsum(
"ij,ki,kj->", distances, q[:-1], q[1:]
)

# Constraints
constraints: amplify.ConstraintList = amplify.one_hot(
q[:-1], axis=1
) + amplify.one_hot(q[:-1], axis=0)

return objective + dmax * constraints

class BenchTspAmplify:
def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
self.model = tsp_for_amplify_v1(ncity, distances, dmax)

def to_qubo(self):
self.model.to_unconstrained_poly()

PyQUBO
import pyqubo

def tsp_for_pyqubo(ncity: int, distances: np.ndarray, dmax: float):
# from https://github.com/recruit-communications/pyqubo/blob/master/notebooks/TSP.ipynb
# NOTE: https://github.com/recruit-communications/pyqubo/blob/master/benchmark/benchmark.py
# is not valid for TSP

x = pyqubo.Array.create("c", (ncity, ncity), "BINARY")

# Constraint not to visit more than two cities at the same time.
time_const = 0.0
for i in range(ncity):
# If you wrap the hamiltonian by Const(...), this part is recognized as constraint
time_const += pyqubo.Constraint(
(sum(x[i, j] for j in range(ncity)) - 1) ** 2, label=f"time{i}"
)

# Constraint not to visit the same city more than twice.
city_const = 0.0
for j in range(ncity):
city_const += pyqubo.Constraint(
(sum(x[i, j] for i in range(ncity)) - 1) ** 2, label=f"city{j}"
)

# distance of route
feed_dict = {}

distance = 0.0
for i in range(ncity):
for j in range(ncity):
for k in range(ncity):
# we set the constant distance
distance += distances[i, j] * x[k, i] * x[(k + 1) % ncity, j]

# Construct hamiltonian
A = pyqubo.Placeholder("A")
H = distance + A * (time_const + city_const)

feed_dict["A"] = dmax

# Compile model
return H.compile(), feed_dict

class BenchTspPyQubo:
def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
self.model, self._feed_dict = tsp_for_pyqubo(ncity, distances, dmax)

def to_qubo(self):
self.model.to_qubo(index_label=False, feed_dict=self._feed_dict)


dimod BQM (index)
import dimod

def tsp_for_dimod_bqm(ncity: int, distances: np.ndarray, dmax: float):
bqm = dimod.BinaryQuadraticModel(ncity * ncity, dimod.BINARY)

# Objective function
for n in range(ncity):
for i in range(ncity):
for j in range(ncity):
bqm.add_quadratic(
n * ncity + i,
((n + 1) % ncity) * ncity + j,
distances[i, j],
)

# Constraint on each row
for n in range(ncity):
left = [(n * ncity + i, 1) for i in range(ncity)]
bqm.add_linear_equality_constraint(left, dmax, -1)

# Constraint on each column
for i in range(ncity):
left = [(n * ncity + i, 1) for n in range(ncity)]
bqm.add_linear_equality_constraint(left, dmax, -1)

return bqm

class BenchTspDimodBQM:
def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
self.model = tsp_for_dimod_bqm(ncity, distances, dmax)

def to_qubo(self):
self.model.to_qubo()

dimod BQM (symbol math)
import dimod

def tsp_for_dimod_bqm_sym(
ncity: int, distances: np.ndarray, dmax: float
) -> dimod.BinaryQuadraticModel:
bqm = dimod.BinaryQuadraticModel(ncity * ncity, dimod.BINARY)
vars = [
[dimod.Binary(f"{n},{i}") for i in range(ncity)] for n in range(ncity)
]

# Objective function
for n in range(ncity):
for i in range(ncity):
for j in range(ncity):
bqm += distances[i, j] * vars[n][i] * vars[(n + 1) % ncity][j]

# Constraint on each row
for n in range(ncity):
bqm += dmax * (sum(vars[n][i] for i in range(ncity)) - 1) ** 2

# Constraint on each column
for i in range(ncity):
bqm += dmax * (sum(vars[n][i] for n in range(ncity)) - 1) ** 2

return bqm  # type: ignore

class BenchTspDimodBQMSym:
def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
self.model = tsp_for_dimod_bqm_sym(ncity, distances, dmax)

def to_qubo(self):
self.model.to_qubo()

dimod CQM
import dimod

def tsp_for_dimod_cqm(ncity: int, distances: np.ndarray, dmax: float):
cqm = dimod.ConstrainedQuadraticModel()
vars = [
[dimod.Binary(f"{n},{i}") for i in range(ncity)] for n in range(ncity)
]

# Objective function
obj = 0.0
for n in range(ncity):
for i in range(ncity):
for j in range(ncity):
obj += distances[i, j] * vars[n][i] * vars[(n + 1) % ncity][j]
cqm.set_objective(obj)

# Constraint on each row
for n in range(ncity):
cqm.add_constraint(sum(vars[n]) == 1)

# Constraint on each column
for i in range(ncity):
cqm.add_constraint(sum(vars[n][i] for n in range(ncity)) == 1)

return cqm

class BenchTspDimodCQM:
def create_model(self, ncity: int, distances: np.ndarray, dmax: float):
self.model = tsp_for_dimod_cqm(ncity, distances, dmax)

def to_qubo(self):
pass

Benchmark code
import time

def make_distance(ncity: int) -> tuple[np.ndarray, float]:
rng = np.random.default_rng(12345)
x = rng.random(ncity)
y = rng.random(ncity)
distances = (
(x[:, np.newaxis] - x[np.newaxis, :]) ** 2
+ (y[:, np.newaxis] - y[np.newaxis, :]) ** 2
) ** 0.5
dmax: float = np.max(distances)  # type: ignore
return distances, dmax

for ncity in [32, 100, 317]:
distances, dmax = make_distance(ncity)

for bench_class in [
BenchTspAmplify,
BenchTspPyQubo,
BenchTspDimodBQM,
BenchTspDimodBQMSym,
BenchTspDimodCQM,
]:
bench = bench_class()
start = time.time()
bench.create_model(ncity, distances, dmax)
end = time.time()
t1 = end - start
start = time.time()
bench.to_qubo()
end = time.time()
t2 = end - start
print(f"{t1} {t2}")


Benchmark environment

CPU

IntelR CoreTM i9-12900K
(E-Cores disabled)

OS

Ubuntu 22.04

Python 3.11
• amplify 1.0.0

• pyqubo 1.4.0

• dimod 0.12.14