The Graph Coloring problem describes the task of coloring regions on a map, such that no neighboring regions have the same color.
In this post, I will demonstrate how Grover’s algorithm can be used to solve the Graph Coloring problem.
Defining the problem
The code below defines the constraints required for coloring a simple 2×2 square.
num_regions = 4
clause_list = [[0, 1], [0, 2], [1, 3], [2, 3]]
Next, let’s import the necessary libraries:
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, QuantumRegister
According to the four-color theorem, any 2-dimensional map can be colored using just four colors (which is a really cool!). Therefore, we can represent a single color using 2 qubits.
Below, I define a function that will transform a bit stream to a list of four colors: red, yellow, green, and blue.
We’ll use this function when we have found our solution at the end.
def decode_color(k):
colors = {"00": (1, 0, 0), "01": (1, 1, 0), "10": (0, 1, 0), "11": (0, 0, 1)}
return [colors[k[2 * i : 2 * i + 2]] for i in range(len(k) // 2)]
Defining constraints
Now, let’s define a function that applies a single constraint to two neighboring regions.
What we want is to take a set of input qubits from region 1 (a, b), and region 2 (c, d), and set an output qubit, also called a clause qubit, to 1 if the result is legal, and 0 if it is illegal.
def XOR(qc, a, b, c, d, output):
qc.cx(a, output)
qc.cx(b, output)
qc.cx(c, output)
qc.cx(d, output)
Here is how we can construct a quantum circuit for a single pair of neighboring regions:
in_qubits = QuantumRegister(4, name="input")
out_qubit = QuantumRegister(1, name="output")
qc = QuantumCircuit(in_qubits, out_qubit)
XOR(qc, in_qubits[0], in_qubits[1], in_qubits[2], in_qubits[3], out_qubit)
qc.draw(output="mpl")
plt.savefig("./img/1-constraint-single.png")
Here is the quantum circuit plotted:
Scaling up to all regions
We can extend this to apply to all four regions and plot the result. What we end up with is a set of four clause qubits that indicate the legality of each individual constraint.
# Create separate registers to name bits
var_qubits = QuantumRegister(2 * num_regions, name="v") # variable bits
clause_qubits = QuantumRegister(
len(clause_list), name="c"
) # bits to store clause-checks
# Create quantum circuit
qc = QuantumCircuit(var_qubits, clause_qubits)
# Use XOR gate to check each clause
for i, clause in enumerate(clause_list):
XOR(
qc,
2 * clause[0],
2 * clause[0] + 1,
2 * clause[1],
2 * clause[1] + 1,
clause_qubits[i],
)
qc.draw(output="mpl")
plt.savefig("./img/2-constraint-multiple.png")
The next step is to take all of the individual constraints and combine them into a single output qubit that indicates the legality of the whole system. We can do this using a toffoli gate (also known as a CCNOT gate).
# Create separate registers to name bits
var_qubits = QuantumRegister(2 * num_regions, name="v")
clause_qubits = QuantumRegister(len(clause_list), name="c")
output_qubit = QuantumRegister(1, name="out")
qc = QuantumCircuit(var_qubits, clause_qubits, output_qubit)
# Compute clauses
for i, clause in enumerate(clause_list):
XOR(
qc,
2 * clause[0],
2 * clause[0] + 1,
2 * clause[1],
2 * clause[1] + 1,
clause_qubits[i],
)
# Flip 'output' bit if all clauses are satisfied
qc.mct(clause_qubits, output_qubit)
qc.draw(output="mpl")
plt.savefig("./img/3-constraint-multiple-toffoli.png")
Creating a quantum oracle
The code below defines the oracle function required by Grover’s algorithm.
The oracle is essentially what we have already written, plus an uncomputation stage to reset the clause-checking bits to 0. I do not fully understand this step, but it likely has something to do with the reversible nature of quantum computers.
def map_oracle(qc, clause_list, clause_qubits, output_qubit):
# Compute clauses
for i, clause in enumerate(clause_list):
XOR(
qc,
2 * clause[0],
2 * clause[0] + 1,
2 * clause[1],
2 * clause[1] + 1,
clause_qubits[i],
)
# Flip 'output' bit if all clauses are satisfied
qc.mct(clause_qubits, output_qubit)
# Uncompute clauses to reset clause-checking bits to 0
for i, clause in enumerate(clause_list):
XOR(
qc,
2 * clause[0],
2 * clause[0] + 1,
2 * clause[1],
2 * clause[1] + 1,
clause_qubits[i],
)
var_qubits = QuantumRegister(2 * num_regions, name="v")
clause_qubits = QuantumRegister(len(clause_list), name="c")
output_qubit = QuantumRegister(1, name="out")
cbits = ClassicalRegister(2 * num_regions, name="cbits")
qc = QuantumCircuit(var_qubits, clause_qubits, output_qubit, cbits)
map_oracle(qc, clause_list, clause_qubits, output_qubit)
qc.draw(output="mpl", fold=-1)
plt.savefig("./img/4-oracle.png")
Here is a plot of the circuit:
Final circuit
The code below is a full implementation of the algorithm.
It starts by defining a “diffuser” function required the diffusion step of the algorithm.
The code then applies two iterations of the quantum oracle, which seems to be enough for this circuit but should be applied √N times, where N is the problem’s domain.
import numpy as np
from qiskit import QuantumCircuit, ClassicalRegister
def diffuser(nqubits):
qc = QuantumCircuit(nqubits)
# Apply transformation |s> -> |00..0> (H-gates)
for qubit in range(nqubits):
qc.h(qubit)
# Apply transformation |00..0> -> |11..1> (X-gates)
for qubit in range(nqubits):
qc.x(qubit)
# Do multi-controlled-Z gate
qc.h(nqubits - 1)
qc.mct(list(range(nqubits - 1)), nqubits - 1) # multi-controlled-toffoli
qc.h(nqubits - 1)
# Apply transformation |11..1> -> |00..0>
for qubit in range(nqubits):
qc.x(qubit)
# Apply transformation |00..0> -> |s>
for qubit in range(nqubits):
qc.h(qubit)
# We will return the diffuser as a gate
U_s = qc.to_gate()
U_s.name = "U$_s$"
return U_s
var_qubits = QuantumRegister(2 * num_regions, name="v")
clause_qubits = QuantumRegister(len(clause_list), name="c")
output_qubit = QuantumRegister(1, name="out")
cbits = ClassicalRegister(2 * num_regions, name="cbits")
qc = QuantumCircuit(var_qubits, clause_qubits, output_qubit, cbits)
# Initialize 'out0' in state |->
qc.initialize([1, -1] / np.sqrt(2), output_qubit)
# Initialize qubits in state |s>
qc.h(var_qubits)
qc.barrier() # for visual separation
for i in range(2):
# Apply our oracle
map_oracle(qc, clause_list, clause_qubits, output_qubit)
qc.barrier() # for visual separation
# Apply our diffuser
qc.append(diffuser(2 * num_regions), range(2 * num_regions))
# Measure the variable qubits
qc.measure(var_qubits, cbits)
qc.draw(output="mpl", fold=-1)
plt.savefig("./img/5-final-circuit.png")
Here is the plot:
Simulating the circuit
The code below simulates the final circuit for the Graph Coloring problem using the Aer simulator.
The code starts by setting a backend simulator, transpiles the final circuit to optimize it, and runs the simulation.
from qiskit import Aer, assemble, transpile
from qiskit.visualization import plot_histogram
# Simulate and plot results
aer_simulator = Aer.get_backend("aer_simulator")
transpiled_qc = transpile(qc, aer_simulator)
qobj = assemble(transpiled_qc)
result = aer_simulator.run(qobj).result()
counts = result.get_counts()
plot_histogram(counts, figsize=(50, 15))
plt.savefig("./img/6-results.png")
This successfully captures all legal 32 results out of 2^8=256 possibilities. We can also see that most illegal possibilities happen 0 times!
Plotting the results
The final step is to plot the result.
We have to get the counts, sort by the most common outcomes, and filter out events that are due to undesirable effects:
counts = {
k: v for k, v in sorted(counts.items(), key=lambda item: item[1], reverse=True)
}
counts = {
k: v
for k, v in counts.items()
if v > sum(counts.values()) / 2 ** (2 * num_regions)
}
print([decode_color(k) for k in counts.keys()])
Here is the first result of 32 results:
[[(1, 0, 0), (1, 1, 0), (1, 1, 0), (1, 0, 0)], ...]
As you can see from the image, the Graph Coloring problem has been successfully solved using Grover’s algorithm! It is also worth noting that this particular solution only required two colors.
In this case, the solution set includes three-color solutions
There are also 4!(4−4)!=24 examples of four-colorings.
Applying the algorithm to a real map
I tried applying the algorithm to the provinces of Canada. The regions were defined as follows:
provinces = {
0: "Newfoundland and Labrador",
1: "Prince Edward Island",
2: "Nova Scotia",
3: "New Brunswick",
4: "Quebec",
5: "Ontario",
6: "Manitoba",
7: "Saskatchewan",
8: "Alberta",
9: "British Columbia",
10: "Yukon",
11: "Northwest Territories",
12: "Nunavut",
}
clause_list = [
[0, 4],
[2, 3],
[3, 4],
[4, 5],
[5, 6],
[6, 7],
[6, 12],
[7, 8],
[7, 11],
[8, 9],
[8, 11],
[9, 10],
[9, 11],
[10, 11],
[11, 12]
]
Unfortunately, this was too much for the algorithm to simulate and it ran out of memory.
I then tried to split it into two sections along the Manitoba-Ontario border. It was able to finish running in this configuration, but the results were not accurate, probably due to not enough iterations of the quantum oracle. I tried increasing this but this caused the script to hang 🤷
My classical computer is fast enough that it can find solutions via random sampling, although this workaround would not work for more complicated maps.
Conclusion
This post has shown how we can use simple heuristics (Grover’s Algorithm) to solve a well defined problem using a simulated quantum computer. Experiments like this have the potential to revolutionize the way we solve complex problems as the industry matures.
Until then, the industry needs not only better and more accessible hardware, but also more learning resources, algorithms, quantum languages, quantum compilers, and a larger pool of programmers.
This post is based on Qiskit’s post on Grover’s algorithm which you should definitely check out too!
