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

Benchmark results#