#
Example NP problems published in A. Lucas, *Front. Phys.* (2014) - Minimal maximal matching
problem¶

This example code implements the **minimal maximal matching 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)
**Minimal maximal 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)

## Minimal maximal matching problem¶

For a graph $G$, $D$ is called a **maximal matching** if the subset $D$ of edges in $G$
satisfies the following.

- Edges contained in $D$ are not adjacent to each other.
- Edges not included in $D$ are always adjacent to one of the edges in $D$.

For example, the orange edges in the figure below are maximal matching. Note that the orange edges are not connected and that if even one of the black edges were painted orange, the orange edges would be connected and would not be a maximal matching.

The minimal maximal matching problem is the problem of finding, for a given graph, the maximal matching of that graph that has the smallest number of elements.

This example program uses Fixstars Amplify to find the minimum maximal matching. The formulation follows that of Sec. 4.5 of A. Lucas, Front. Phys. (2014).

## Problem definition¶

First, as a problem, create the following graph $G$ using NetworkX.

```
from matplotlib import pyplot as plt
import networkx as nx
import numpy as np
N = 6 # Number of vertices of the graph
G = nx.Graph()
G.add_nodes_from(range(N))
elist = [
(0, 1),
(0, 5),
(1, 2),
(1, 5),
(2, 3),
(2, 4),
(3, 4),
(4, 5),
] # Edges connecting two vertices
G.add_edges_from(elist)
pos = nx.circular_layout(G)
nx.draw_networkx(G, node_size=600, font_color="w", pos=pos)
```

## Formulation¶

Let $N$ be the number of vertices and $M$ the number of edges in $G$.

### Decision variables¶

Let each of $M$ binary variables $q$ correspond to each edge of $G$ to indicate whether the maximal matching $D$ contains the edge. $q=1$ if included in $D$, $0$ if not.

For example, for the following maximal matching, the binary variable $q$ would be as in the table below.

Edge $(u, v)$ | $$(0, 1)$$ | $$(0, 5)$$ | $$(1, 2)$$ | $$(1, 5)$$ | $$(2, 3)$$ | $$(2, 4)$$ | $$(3, 4)$$ | $$(4, 5)$$ |
---|---|---|---|---|---|---|---|---|

$$q$$ | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 |

### Objective function¶

Since the number of elements in $D$ should be as small as possible, we minimize $ \displaystyle \sum_{i = 0}^{M - 1} q_i$.

### Constraints¶

As explained earlier, $D$ being a maximal matching means that the following constraints are satisfied ($D$ is the subset of edges in $G$).

- Condition 1: Edges contained in $D$ are not adjacent to each other.
- Condition 2: Edges not included in $D$ are always adjacent to one of the edges in $D$.

Let us rephrase these conditions and express them in terms of $q$.

First, we can rephrase condition 1 as "no two adjacent edges are both contained in $D$":

$$ q_{v, u} q_{v, w} = 0 \quad \text{for} \quad (v, u), (v, w) \in E $$

Here, the element of the binary variable array $q$ corresponding to the edge $(u, v)$ is written as $q_{u, v}$. Also, $E$ denotes the edge set of $G$.

Next, we can rephrase condition 2 as "every edge of $G$ is necessarily adjacent to one of the edges of $D$". We can further rephrase this as "for any edge $(u, v)$ of $G$, either $u$ or $v$ is an endpoint of one of the edges of $D$". We can determine whether a vertex $v$ is an endpoint of any edge of $D$ by whether the sum of the corresponding binary variables is $1$ or $0$ for all edges out of $v$. So the condition 2 is:

$$ (1 - \sum_{(v, x) \in E} q_{v, x}) (1 - \sum_{(u, y) \in E} q_{u, y}) = 0 \quad \text{for} \quad (u, v)\in E. $$

## Implementation¶

Based on the problem and formulation described above, let us implement and solve the problem. First,
create $M$ binary variables $q$ using `BinarySymbolGenerator`

in Fixstars Amplify SDK.

```
from amplify import BinarySymbolGenerator
M = len(G.edges)
gen = BinarySymbolGenerator()
q = gen.array(M)
```

We can create the objective function according to the previous formulation. The objective function equals the number of elements in the maximal matching $D$ and is written as $\displaystyle \sum_{i = 0}^{M - 1} q_i$.

```
cost = q.sum()
```

In preparation for constructing the constraints, for each vertex $v$ of $G$, we create a function that returns a list of the indices of the edges coming out of $v$.

```
def incident_edge_indices(v):
return [i for i, e in enumerate(G.edges) if v in e]
```

Now, let us construct the constraint corresponding to condition 1. Condition 1 is that no two edges contained in the maximal matching $D$ are adjacent, i.e., both adjacent $2$ edges are not included in $D$ together, and is written as $q_{v, u} q_{v, w} = 0 \bigl((v, u), (v, w) \in E\bigr)$.

```
from itertools import combinations
from amplify.constraint import equal_to
constraint1 = [
equal_to(q[i] * q[j], 0)
for v in G.nodes
for i, j in combinations(incident_edge_indices(v), 2)
]
```

Then, we construct the constraint corresponding to condition 2. Condition 2 is that all edges are adjacent to one of the edges of $D$ and written as:

$$ \displaystyle(1 - \sum_{(v, x) \in E} q_{v, x}) (1 - \sum_{(u, y) \in E} q_{u, y}) = 0 \bigl((u, v)\in E\bigr). $$

```
constraint2 = [
equal_to(
(1 - sum([q[i] for i in incident_edge_indices(u)]))
* (1 - sum([q[i] for i in incident_edge_indices(v)])),
0,
)
for (u, v) in G.edges
]
```

The objective function and constraints constructed above are then combined and converted to a logical model.

Although not necessary in this case, when both the objective function and the constraints are present, it is usually a good idea to multiply some weights by the constraints because the constraints are given to the Ising machine as a penalty function for the objective function. The basic idea is to estimate and determine a weight value equal to or slightly larger than the possible values of the objective function.

```
from amplify import BinaryQuadraticModel
model = BinaryQuadraticModel(cost + sum(constraint1) + sum(constraint2))
```

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

automatically filters the solutions that satisfy the constraints, if the
`result`

is not empty, you know that there is a solution that satisfies the constraints.

```
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)
if len(result) == 0:
print("No solution has been found.")
else:
print("A solution has been found.")
```

Finally, let us visualize the result. You can try to solve the same problem for different graph shapes.

```
values = q.decode(result[0].values)
colors = ["k" if i == 0 else "C1" for i in values]
width = [1.0 if i == 0 else 2.0 for i in values]
nx.draw_networkx(
G, node_size=600, font_color="w", edge_color=colors, width=width, pos=pos
)
```