Formulation Benchmarks

We 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. We used 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 CQM1	PyQBPP double
Symbolic operation	✅	✅	❌	✅	✅	✅
Array operation and einsum	✅	❌5	❌	❌	❌	✅
Objective function	✅	✅	✅	✅	✅	✅
Constraint	✅	✅2	❌3	❌3	✅1	✅
Automatic penalty function	✅	❌	✅6	✅6	❌1	✅
Higher order polynomial	✅	✅	❌	❌	❌	✅
Coefficient matrix	✅	❌	❌	❌	❌	❌
Variable type	B/S/I/R	B/S	B/S	B/S	B/S/I/R	B4
Automatic variable encoding	✅	❌7	❌	❌	❌1	❌8
Supported machines	Various	Depends on user	D-Wave only	D-Wave only	D-Wave only	Various
Model file input and output	LP/QPLIB	❌	❌	❌	LP9	❌
Type hint	✅	❌	✅	✅	✅	❌

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

The dimod columns give the features of dimod alone. Other packages of the D-Wave Ocean SDK give some of these features.

1. Model creation only, as QUBO output (conversion of constraints to the penalty functions) is not available.

2. Must define penalty function

3. Constraints are expressed by adding penalty functions to the objective function

4. Declares binary variables only. An integer comes from onehot_to_int, and a spin from binary_to_spin

5. Array is a container of variables. It gives no element-wise operation between two arrays and no broadcast

6. Generates a penalty term with add_linear_equality_constraint and similar methods, for a linear equality or inequality constraint only

7. Selects the encoding of an integer with an explicit class. Gives no real variable

8. Integer variables only. Needs an explicit conversion with binarize

9. LP file only

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(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

PyQBPP

import pyqbpp.double as qbpp


def tsp_for_pyqbpp(ncity: int, distances: np.ndarray, dmax: float):
    variables = qbpp.var("x", shape=(ncity, ncity))

    # Constraints
    model = (
        qbpp.sum(qbpp.constrain(qbpp.vector_sum(variables, axis=1), equal=1))
        + qbpp.sum(qbpp.constrain(qbpp.vector_sum(variables, axis=0), equal=1))
    ) * dmax

    # Objective function
    distance = qbpp.array(distances.ravel().tolist(), shape=(ncity, ncity))
    successors = qbpp.concat([variables[1:], variables[:1]], axis=0)
    model += qbpp.einsum("jk,ij,ik->", distance, variables, successors)

    return model


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

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

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,
        BenchTspPyQbpp,
    ]:
        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 results

PyQBPP (double) is not measured above 10,000 binary variables.

Benchmark environment

CPU

12th Gen Intel(R) Core(TM) i9-12900K
(E-Cores disabled)

OS

Linux-6.8.0-137-generic-x86_64-with-glibc2.43

Python 3.12.12

• amplify 1.7.0

• amplify 1.0.5

• pyqubo 1.5.0

• dimod 0.12.22

• pyqbpp 2026.8.16

Formulation time at every size

Each point gives the total formulation time, so a lower point is faster.

32 cities (1,024 binary variables)

Formulation	Model construction	QUBO construction	Total
Amplify v1.7	0.24 ms	0.31 ms	0.55 ms 🏆
PyQUBO	142.60 ms	23.28 ms	165.88 ms (x301.3)
dimod BQM (index)	12.80 ms	29.18 ms	41.97 ms (x76.2)
dimod BQM (symbol)	611.76 ms	41.36 ms	653.12 ms (x1186.3)
dimod CQM	531.95 ms	N/A	531.95 ms (x966.2)
PyQBPP (double)	3.34 ms	4.94 ms	8.28 ms (x15.0)

100 cities (10,000 binary variables)

Formulation	Model construction	QUBO construction	Total
Amplify v1.7	3.27 ms	5.35 ms	8.61 ms 🏆
PyQUBO	4.820 s	1.430 s	6.250 s (x725.6)
dimod BQM (index)	431.44 ms	1.038 s	1.469 s (x170.6)
dimod BQM (symbol)	18.425 s	1.447 s	19.872 s (x2307.0)
dimod CQM	15.433 s	N/A	15.433 s (x1791.7)
PyQBPP (double)	57.64 ms	61.36 ms	119.00 ms (x13.8)

317 cities (100,489 binary variables)

Formulation	Model construction	QUBO construction	Total
Amplify v1.7	104.85 ms	180.43 ms	285.28 ms 🏆
PyQUBO	166.365 s	66.841 s	233.206 s (x817.5)
dimod BQM (index)	23.740 s	39.400 s	63.140 s (x221.3)
dimod BQM (symbol)	590.709 s	52.804 s	643.513 s (x2255.8)
dimod CQM	532.464 s	N/A	532.464 s (x1866.5)
PyQBPP (double)	N/A	N/A	N/A