Example NP problems published in A. Lucas, Front. Phys. (2014) - Hamiltonian cycle problem¶

This example code implements the Hamiltonian cycle problem introduced in the paper A. Lucas, "Ising formulations of many NP problems", Front. Phys. (2014) using Fixstars Amplify. Other NP-complete and NP-hard problems introduced in the same paper are also discussed below (the corresponding sections in the paper are shown in the brackets).

Graph partitioning problem (Sec. 2.2).
Maximum clique problem (Sec. 2.3)
Exact cover problem (Sec. 4.1)
Set packing problem (Sec. 4.2)
Minimum vertex cover problem (Sec. 4.3)
Satisfiability problem (SAT) (Sec. 4.4)
Minimum maximum matching problem (Sec. 4.5)
Graph coloring problem (Sec. 6.1)
Clique cover problem (Sec. 6.2)
Job sequencing problem with integer lengths (Sec. 6.3)
Hamiltonian cycle problem (Sec. 7.1)
Directed feedback vertex set problem (Sec. 8.3)
Minimum feedback edge set problem (Sec. 8.5)
Graph isomorphism problem (Sec. 9)

Hamiltonian cycle problem¶

Given a graph $G$, we call a closed path a Hamiltonian cycle if it passes through all vertices of $G$ once and returns to the origin. In general, when the size of the graph is large, it is difficult to determine in realistic time whether a Hamiltonian cycle exists in the graph.

Here, we use Fixstars Amplify to solve this Hamiltonian cycle problem. This problem corresponds to section 7.1 of A. Lucas, Front. Phys. (2014).

Problem definition¶

First, we create a graph $G$ to be considered in this sample program using NetworkX. The number of vertices is $N$.

In [ ]:

import networkx as nx
import numpy as np
    
N = 5  # Number of vertices of the graph
G = nx.Graph()
G.add_nodes_from(range(N))
elist = [(0, 1), (0, 2), (1, 2), (1, 3), (1, 4), (2, 4), (3, 4)]
    
pos = nx.circular_layout(G)  # Save the layout of the graph
G.add_edges_from(elist)
    
nx.draw_networkx(G, node_size=600, font_color="w", pos=pos)

Formulation¶

Decision variables¶

Let us consider $N\times N$ binary decision variables $q$, representing which vertex to pass and when. That is, a component $q_{k, i}$ of the binary decision variables corresponds to passing the vertex $i$ at $k$-th visit ($=1$) or not ($=0$). For example, when the binary variables are as follows, it corresponds to a closed path $0 \rightarrow 1 \rightarrow 3 \rightarrow 4 \rightarrow 2 \rightarrow 0$ in the above graph.

Order \ Index of vertex	0	1	2	3	4
1st	1	0	0	0	0
2nd	0	1	0	0	0
3rd	0	0	0	1	0
4th	0	0	0	0	1
5th	0	0	1	0	0

Objective function¶

Since the Hamiltonian cycle problem is to find one that satisfies the conditions, no objective function is considered.

Constraints¶

For $q$ to represent a Hamiltonian cycle, we need the following:

The $k$-th vertex must be a single vertex. We can rephrase this condition as there being exactly a single $1$ in each row of the binary variable table $q$.
Each vertex must be passed through exactly $1$ times. We can rewrite this condition as there being exactly a single $1$ in each column of the binary variable table $q$.
No transfers are allowed between vertices that are not connected. That is, when no edge connects between vertices $i$ and $j$, both $q_{k, i}$ and $q_{k+1, j}$ must not be $1$.

The above conditions 1-3 can be written in mathematical expressions as follows, respectively.

\begin{align*} \sum_{i=0}^{N-1} q_{k, i} = 1 & \quad \text{for} \quad k \in \{0, 1, \ldots, N-1\} \\ \sum_{k=0}^{N-1} q_{k, i} = 1 & \quad \text{for} \quad i \in \{0, 1, \ldots, N-1\} \\ q_{k, i}q_{k+1, j} = 0 & \quad \text{for} \quad k \in \{0, 1, \ldots, N-1\}, (i, j) \notin E. \end{align*}

Here, $E$ denotes the edge set of $G$.

Also, when the binary variables $q$ satisfy all conditions 1-3, $q$ corresponds to a Hamiltonian cycle of $G$.

Implementation¶

Using the problem and formulation described above, let us implement and solve the problem. First, use BinarySymbolGenerator in Fixstars Amplify SDK to create $N\times N$ binary decision variables $q$.

In [ ]:

from amplify import BinarySymbolGenerator
    
gen = BinarySymbolGenerator()
q = gen.array(N, N)

Next, we create the constraints corresponding to the conditions 1 and 2. Since these are conditions that there is only one $1$ in each row and column of $q$, we can write them using one_hot.

In [ ]:

from amplify.constraint import one_hot
    
row_constraints = [one_hot(q[k, :]) for k in range(N)]
col_constraints = [one_hot(q[:, i]) for i in range(N)]

We can print the abovementioned constraints and check that the one_hot condition is correctly imposed on each row and column.

In [ ]:

row_constraints

In [ ]:

col_constraints

Next, we create the constraint corresponding to the condition 3. Condition 3 is the condition that $q_{k, i}q_{k+1, j} = 0$ ($i$ and $j$ are two vertices not connected by an edge). Note that $q_{k+1, j}$ implies $q_{0, j}$ when $k=N-1$.

In [ ]:

from amplify.constraint import equal_to
    
edge_constraints = [
    equal_to(
        q[k, i] * q[(k + 1) % N, j], 0
    )  # %N is the treatment related to the above note.
    for k in range(N)
    for i in range(N)
    for j in range(N)
    if i != j and (i, j) not in G.edges and (j, i) not in G.edges
]

In [ ]:

edge_constraints[:6]

The necessary constraints are in place. Finally, we combine to create the QUBO model.

In [ ]:

from amplify import BinaryQuadraticModel
    
model = BinaryQuadraticModel(
    sum(row_constraints) + sum(col_constraints) + sum(edge_constraints)
)

Let us set the client and execute the solver with Fixstars Amplify Annealing Engine (AE).

In [ ]:

from amplify.client import FixstarsClient
from amplify import Solver
    
client = FixstarsClient()
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # If you use Amplify in a local environment or Google Colaboratory, enter your Amplify API token.
client.parameters.timeout = 1000
    
# Define and execute the solver
solver = Solver(client)
result = solver.solve(model)

Now, we check whether there is a Hamiltonian cycle in the graph. Since Solver automatically filters the solutions that satisfy the constraints, if the result is not empty, you know that a solution has been found that satisfies the constraints.

In [ ]:

if len(result) == 0:
    print("No Hamiltonian cycle has been found.")
else:
    print("A Hamiltonian cycle has been found.")

Finally, let us visualize the solution. Below is the found Hamiltonian cycle in orange color.

In [ ]:

# Setting default attributes for edges
for edge in G.edges.values():
    edge["color"] = "k"
    edge["width"] = 1.0
    
# Setting attributes for edges constituting the found Hamiltonian cycle
values = q.decode(result[0].values)
route = np.where(values == 1)[1]
for i, j in zip(route, np.roll(route, -1)):
    G.edges[i, j]["color"] = "C1"
    G.edges[i, j]["width"] = 2.0
    
# Visualize
edge_color = [edge["color"] for edge in G.edges.values()]
edge_width = [edge["width"] for edge in G.edges.values()]
nx.draw_networkx(
    G, node_size=600, font_color="w", pos=pos, edge_color=edge_color, width=edge_width
)