Example NP problems published in A. Lucas, Front. Phys. (2014)  Directed feedback vertex set problem¶
This example code implements the directed feedback vertex set problem introduced in the paper A. Lucas, "Ising formulations of many NP problems", Front. Phys. (2014) using Fixstars Amplify. Other NPcomplete and NPhard 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)
Directed feedback vertex set problem¶
Given a directed graph $G$, a subset $F$ of the vertices of $G$ such that any closed path in $G$ passes through at least $1$ vertices in $F$ is called a directed feedback vertex set.
In other words, $F$ is a directed feedback vertex set of $G$ if it is impossible to start from a vertex not contained in $F$ and return to the original vertex only through vertices not included in $F$.
For example, in the graph below, the subset of vertices shown in orange is one of the directed feedback vertex sets.
The directed feedback vertex set problem is the problem of finding the feedback vertex set for a directed graph $G$, which has the minimum number of elements. The formulation of this example program follows that of Sec. 8.3 of A. Lucas, Front. Phys. (2014).
Problem definition¶
First, we create a directed graph $G$ using NetworkX as the feedback vertex set problem we will solve in this example program.
import networkx as nx
N = 8 # Number of vertices of the graph G
G = nx.DiGraph()
G.add_nodes_from(range(N))
edge_list = [
(0, 1),
(1, 2),
(2, 3),
(3, 4),
(4, 5),
(5, 6),
(6, 7),
(7, 0),
(4, 2),
(7, 1),
(7, 5),
]
G.add_edges_from(edge_list)
node_labels = {
0: "a",
1: "b",
2: "c",
3: "d",
4: "e",
5: "f",
6: "g",
7: "h",
}
pos = nx.circular_layout(G) # Save the layout of the graph G
nx.draw_networkx(G, node_size=600, font_color="w", labels=node_labels, pos=pos)
It is easy to see that the graph created by removing the vertices $e$ and $f$ from $G$ has no closed path. That is, $\{e, f\}$ is a feedback vertex set of $G$.
Also, since the closed path $c\rightarrow d\rightarrow e\rightarrow c$ and the closed path $f\rightarrow g\rightarrow h\rightarrow f$ have no common part, the number of elements in the feedback vertex set of $G$ is at least $2$. Therefore, for this problem, the minimum number of elements in the feedback vertex set of $G$ is $2$.
Formulation¶
Formulation guidelines¶
In the following, let $N$ be the number of vertices in $G$.
First, we will associate $N$ binary variables $y$ with each vertex to indicate whether a feedback vertex set $F$ includes the vertex. The binary variable is $0$ if the vertex is included in $F$ and $1$ if not.
Next, the condition "the subgraph $H$ of $G$ consisting of vertices not included in $F$ has no closed path", which is a paraphrase of the main problem, is further paraphrased into the condition "if we number the vertices of $H$ appropriately, we can make every directed edge of $H$ goes from a vertex with a small number to one with a large number" (proof: $\Rightarrow$ is straightforward, $\Leftarrow$ is via topological sorting).
This numbering can be expressed by using a binary variable table $x$ of $N\times N$ and setting the binary variable in the $v$ row $i$ column to $1$ when the number of vertex $v$ is $i$.
For example, the problem created above has the following graph.
If each vertex of this graph is colored and numbered as follows, the two orange points $e$ and $f$ are the feedback vertex set, and the edge connecting the blue vertices goes out from the smallest numbered vertex to the largest.
The table below shows the binary variables $y$ and $x$ corresponding to this way of selecting a feedback vertex set and numbering the edges. Here, all rows of $x$ corresponding to vertices in the feedback vertex set are assumed to be $0$.
a  b  c  d  e  f  g  h  

$y$  1  1  1  1  0  0  1  1 
$$x$$  0  1  2  3  4  5  6  7 

a  0  0  1  0  0  0  0  0 
b  0  0  0  1  0  0  0  0 
c  0  0  0  0  1  0  0  0 
d  0  0  0  0  0  1  0  0 
e  0  0  0  0  0  0  0  0 
f  0  0  0  0  0  0  0  0 
g  1  0  0  0  0  0  0  0 
h  0  1  0  0  0  0  0  0 
Objective function¶
Since the number of elements in the feedback vertex set should be as small as possible, the objective function is:
$$ \displaystyle \sum_{v=0}^{N1} y_v, $$
Note that $y_v$ is $0$ if the feedback vertex set includes the vertex $v$ and $1$ otherwise.
Constraints¶
For $y$ and $x$ to represent the feedback vertex set, we need the following.

Condition 1: Vertices not included in $F$ are numbered $1$. That is, the row $v$ of $x$ are all $0$ if the feedback vertex set includes $v$, and it yields a single $1$ otherwise.

Condition 2: For an edge $u\rightarrow v$ of $G$, if $u$ and $v$ are not both in the feedback vertex set, then the number of $u$ is less than the number of $v$. That is, $x_{u, j}$ and $x_{v, i}$ must not both be $1$ for natural numbers $i \leq j$ at this time (note: $x_{u, \underline{i}}$ and $x_{v, \underline{j}}$ may both be $1$).
The condition 1 is: $$ \sum_{i=0}^{N1} x_{v, i} = y_v \quad \text{for} \quad v \in \{0, 1, \ldots, N1\}. $$
Also, from condition 1, if the feedback vertex set includes either $u$ or $v$, then $x_{u, j}$ and $x_{v, i}$ cannot both be $1$, so the condition "$u$ and $v$ are not both included in the feedback vertex set" in the condition 2 is naturally considered. Therefore, condition 2 is satisfied by the following equation.
$$ x_{u, j} x_{v, i} = 0 \quad \text{for} \quad (u, v) \in E, \ 0 \leq i \leq j < N. $$
Conversely, when the binary variables $y$ and $x$ satisfy conditions 1 and 2, the set of vertices whose corresponding $y$ is $y=0$ is the feedback vertex set, so these can be given as the constraint conditions.
Implementation¶
Using the problem and formulation described above, let us implement and solve the problem. First, use
BinarySymbolGenerator
in Fixstars Amplify SDK to create the binary variables $y$ and $x$.
from amplify import VariableGenerator
gen = VariableGenerator()
y = gen.array("Binary", shape=(N,))
x = gen.array("Binary", shape=(N, N))
Next, we create the objective function $\displaystyle \sum_v y_v$.
cost = y.sum()
Then, let us construct the constraint corresponding to the condition 1. Condition 1 represents that each vertex not included in $F$ is numbered, and we can rephrase this condition as the sum of each row of $x$ equals each element of $y$, as mentioned earlier.
First, we create a onedimensional array representing the difference between the sum of each row of
$x$
and each element of $y$. Then, by specifying an empty tuple in the axis
parameter of the
equal_to
function, we can generate constraints at once such that each element of this
onedimensional array is equal to 0.
from amplify import equal_to
diff = x.sum(axis=1)  y
constraint1 = equal_to(diff, 0, axis=())
Next, we create the constraint corresponding to the condition 2. The condition 2 is the constraint $x_{u, j} x_{v, i} = 0 \bigl((u, v) \in E, \ 0 \leq i \leq j < N\bigr)$.
from amplify import sum as amplify_sum
constraint2 = amplify_sum(
equal_to(x[u, j] * x[v, i], 0)
for u, v in G.edges
for i in range(N)
for j in range(i, N)
)
Now, we can combine the objective function and the constraints into a QUBO model.
Since the constraints are given to the Ising machine as penalty functions on the objective function, we determine the weights of the constraints by estimating values that are approximately equal to or slightly greater than the possible values of the objective function. In this case, the constraint weights are set to $2$.
model = cost + (constraint1 + constraint2) * 2
Let us set the client and execute the solver with Fixstars Amplify Annealing Engine (AE). Since
Amplify
SDK 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.
from amplify import FixstarsClient, solve
from datetime import timedelta
client = FixstarsClient()
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment or Google Colaboratory, enter your Amplify API token.
client.parameters.timeout = timedelta(milliseconds=1000) # timeout is 1000 ms
# Solve the problem
result = solve(model, client)
if len(result) == 0:
print("No solution has been found.")
else:
print("A solution has been found.")
Lastly, we will visualize the solution.
import numpy as np
values = result.best.values
y_values = y.evaluate(values)
x_values = x.evaluate(values)
numbering = {v: "" for v in G.nodes}
numbering.update(dict(np.argwhere(x_values == 1)))
colors = ["C0" if v == 1 else "C1" for v in y_values]
nx.draw_networkx(
G, node_size=600, node_color=colors, font_color="w", labels=numbering, pos=pos
)