Possible error at TSP Quadratic Program calculation

Environment

Qiskit Optimization version: 0.6.1
Python version: 3.12.3
Operating system: macOS Sonoma 14.5

What is happening?

Dear qiskit-optimization team, I believe the implementation of the to_quadratic_program() function of the TSP application could be improved. The problem I'm facing is that I'm currently using custom graphs, and the function is apparently only compatible with geometric graphs. The random instance generator of TSP uses this kind of graphs, so the problem does not arise when using the provided instance generator. However, I want to use a different type of graph in my implementation.

How can we reproduce the issue?

My graph has the following edges:

[(0, 1), (0, 3), (1, 2), (1, 3), (2, 3), (2, 4), (3, 4)]

When trying to convert the TSP instance of my graph into a quadratic program I get the following error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/.local/share/virtualenvs/BA-PvaVEQ8H/lib/python3.12/site-packages/networkx/classes/reportviews.py:1088, in OutEdgeView.__getitem__(self, e)
   1087 try:
-> 1088     return self._adjdict[u][v]
   1089 except KeyError as ex:  # Customize msg to indicate exception origin
KeyError: 2

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
Cell In[48], line 1
----> 1 qp = tsp_custom.to_quadratic_program()

File ~/.local/share/virtualenvs/BA-PvaVEQ8H/lib/python3.12/site-packages/qiskit_optimization/applications/tsp.py:50, in Tsp.to_quadratic_program(self)
     48 x = {(i, k): mdl.binary_var(name=f"x_{i}_{k}") for i in range(n) for k in range(n)}
     49 # Original:
---> 50 tsp_func = mdl.sum(
     51     self._graph.edges[i, j]["weight"] * x[(i, k)] * x[(j, (k + 1) % n)]
     52     for i in range(n)
     53     for j in range(n)
     54     for k in range(n)
     55     if i != j
     56 )
     58 # New:
     59 # tsp_func = mdl.sum(
     60 #     self._graph.edges[i, j]["weight"] * x[(i, k)] * x[(j, (k + 1) % n)]
   (...)
     63 #     if i != j
     64 # )
     66 mdl.minimize(tsp_func)
File ~/.local/share/virtualenvs/BA-PvaVEQ8H/lib/python3.12/site-packages/docplex/mp/model.py:3342, in Model.sum(self, args)
   3329 def sum(self, args):
   3330     """ Creates an expression summing over an iterable over expressions or variables.
   3331     This method expects an iterable over any type of expression: quadrayic expression,
   3332     linear expression, variables, constants.
   (...)
   3340     :return: A Docplex expression.
   3341     """
-> 3342     return self._aggregator.sum(args)
File ~/.local/share/virtualenvs/BA-PvaVEQ8H/lib/python3.12/site-packages/docplex/mp/aggregator.py:198, in ModelAggregator.sum(self, sum_args)
    196 def sum(self, sum_args):
    197     if is_iterator(sum_args):
--> 198         sum_res = self._sum_with_iter(sum_args)
    199     elif is_numpy_ndarray(sum_args):
    200         sum_res = self._sum_with_iter(sum_args.flat)

File ~/.local/share/virtualenvs/BA-PvaVEQ8H/lib/python3.12/site-packages/docplex/mp/aggregator.py:221, in ModelAggregator._sum_with_iter(self, args)
    219 qcc = None
    220 number_validation_fn = checker.get_number_validation_fn()
--> 221 for item in args:
    222     if isinstance(item, LinearOperand):
    223         for lv, lk in item.iter_terms():
File ~/.local/share/virtualenvs/BA-PvaVEQ8H/lib/python3.12/site-packages/qiskit_optimization/applications/tsp.py:51, in <genexpr>(.0)
     48 x = {(i, k): mdl.binary_var(name=f"x_{i}_{k}") for i in range(n) for k in range(n)}
     49 # Original:
     50 tsp_func = mdl.sum(
---> 51     self._graph.edges[i, j]["weight"] * x[(i, k)] * x[(j, (k + 1) % n)]
     52     for i in range(n)
     53     for j in range(n)
     54     for k in range(n)
     55     if i != j
     56 )
     58 # New:
     59 # tsp_func = mdl.sum(
     60 #     self._graph.edges[i, j]["weight"] * x[(i, k)] * x[(j, (k + 1) % n)]
   (...)
     63 #     if i != j
     64 # )
     66 mdl.minimize(tsp_func)

File ~/.local/share/virtualenvs/BA-PvaVEQ8H/lib/python3.12/site-packages/networkx/classes/reportviews.py:1090, in OutEdgeView.__getitem__(self, e)
   1088     return self._adjdict[u][v]
   1089 except KeyError as ex:  # Customize msg to indicate exception origin
-> 1090     raise KeyError(f"The edge {e} is not in the graph.")

KeyError: 'The edge (0, 2) is not in the graph.'

What should happen?

I would expect the to_quadratic_program() function to be compatible with any graph, as is the case with the MaxCut application.

Any suggestions?

I tried to adjust the source code and modify it to match the MaxCut implementation of the function. In the following code snippet you can see the original implementation of the problematic code and my alternative (under # New).

    def to_quadratic_program(self) -> QuadraticProgram:
        """Convert a traveling salesman problem instance into a
        :class:`~qiskit_optimization.problems.QuadraticProgram`

        Returns:
            The :class:`~qiskit_optimization.problems.QuadraticProgram` created
            from the traveling salesman problem instance.
        """
        mdl = Model(name="TSP")
        n = self._graph.number_of_nodes()
        x = {(i, k): mdl.binary_var(name=f"x_{i}_{k}") for i in range(n) for k in range(n)}
        # Original:
        tsp_func = mdl.sum(
            self._graph.edges[i, j]["weight"] * x[(i, k)] * x[(j, (k + 1) % n)]
            for i in range(n)
            for j in range(n)
            for k in range(n)
            if i != j
        )

        # New:
        # tsp_func = mdl.sum(
        #     self._graph.edges[i, j]["weight"] * x[(i, k)] * x[(j, (k + 1) % n)]
        #     for i, j in self._graph.edges
        #     for k in range(n)
        #     if i != j
        # )

        mdl.minimize(tsp_func)
        for i in range(n):
            mdl.add_constraint(mdl.sum(x[(i, k)] for k in range(n)) == 1)
        for k in range(n):
            mdl.add_constraint(mdl.sum(x[(i, k)] for i in range(n)) == 1)
        op = from_docplex_mp(mdl)
        return op

If you compare the original implementation with that of MaxCut, you can see that TSP iterates over all possible edges, while I think it should be iterating over the existing edges only.

    def to_quadratic_program(self) -> QuadraticProgram:
        """Convert a Max-cut problem instance into a
        :class:`~qiskit_optimization.problems.QuadraticProgram`

        Returns:
            The :class:`~qiskit_optimization.problems.QuadraticProgram` created
            from the Max-cut problem instance.
        """
        mdl = Model(name="Max-cut")
        x = {i: mdl.binary_var(name=f"x_{i}") for i in range(self._graph.number_of_nodes())}
        for w, v in self._graph.edges:
            self._graph.edges[w, v].setdefault("weight", 1)
        objective = mdl.sum(
            self._graph.edges[i, j]["weight"] * x[i] * (1 - x[j])
            + self._graph.edges[i, j]["weight"] * x[j] * (1 - x[i])
            for i, j in self._graph.edges
        )
        mdl.maximize(objective)
        op = from_docplex_mp(mdl)
        return op

My implementation works with any graph, but the problem instance / Ising hamiltonian I get is different from that of the original implementation. I greatly appreciate your support. Best regards, Federico

qiskit-community / qiskit-optimization