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

Intel^R Core^TM i9-12900K
(E-Cores disabled)

OS

Ubuntu 22.04

Python 3.11

• amplify 1.0.0

• pyqubo 1.4.0

• dimod 0.12.14

Benchmark results

32 cities (1,024 binary variables)

Formulation	Model construction	QUBO construction	Total
Amplify	1.309 ms	0.998 ms	2.307 ms 🏆
PyQUBO	266.8 ms	24.63 ms	291.5 ms (x126.3)
dimod BQM (index)	22.53 ms	41.58 ms	64.12 ms (x27.8)
dimod BQM (symbol)	1028 ms	58.48 ms	1087 ms (x471.2)
dimod CQM	895.9 ms	N/A ms	895.9 ms (x388.4)

100 cities (10,000 binary variables)

Formulation	Model construction	QUBO construction	Total
Amplify	0.063 s	0.077 s	0.140 s 🏆
PyQUBO	8.519 s	1.360 s	9.879 s (x70.7)
dimod BQM (index)	0.706 s	1.419 s	2.126 s (x15.2)
dimod BQM (symbol)	30.83 s	1.935 s	32.77 s (x234.6)
dimod CQM	26.93 s	N/A	26.93 s (x192.8)

317 cities (100,489 binary variables)

Formulation	Model construction	QUBO construction	Total
Amplify	2.992 s	2.946 s	5.938 s 🏆
PyQUBO	279.7 s	69.23 s	349.0 s (x58.8)
dimod BQM (index)	31.05 s	52.05 s	83.10 s (x14.0)
dimod BQM (symbol)	991.5 s	70.83 s	1062 s (x178.9)
dimod CQM	855.3 s	N/A	855.3 s (x144.0)