3-PyQUBO with OpenJij

Open in Colab

In this chapter, we explain how to convert the cost function to QUBO with PyQUBO, Simulated Annealing, and pass variables to OpenJij. Let us show you “Creek Coverage Problem” as an example.

We install pyqubo with the following command using pip

[ ]:
!pip install pyqubo

Formulation of QUBO with PyQUBO

PyQUBO is a intuitive library for fomulating QUBO & Ising models. In the previous chapters, we have shown the case withou PyQUBO. In the previous chapters, we had to formulate QUBO, then expand the expressions ourselves and put them into the Python script. However, we can eliminate that hassles with PyQUBO.

PyQUBO is a handy library that can help us reduce the computational and implementation errors in our QUBO and Ising model transformations.

Let us use PyQUBO as as example for the Creek Coverage Problem.

For more details of this problem, see also here (T-Wave: creek coverage problem (only Japanese)).

We introduce the formulation of the creek coverage problem in QUBO representation.

This problem is whethre the graph \(G=(V, E)\) can be covered by \(n\) creeks.

It can be expressed in QUBO as follows.

\[H = A\sum_v \left(1-\sum^n_{i=1} x_{v, i}\right)^2 + B \sum^n_{i=1}\left[ \frac{1}{2}\left(-1+\sum_{v \in V} x_{v,i}\right)\sum_{v \in V} x_{v, i} - \sum_{(u, v)\in E} x_{u,i} x_{v, i}\right]\]

First term is constraint where only one color is painted on each vertex. Second shows how close the split subgraph is to creek (complete graph). Both term must be zero. However, we treat the first term as a penalty term, and second as a cost(objective function).

Let us represent this QUBO using PyQUBO.

We give the Graph and the number of creek \(n\) as follows in this time.

[1]:
# set the number of vertex
N_VER = 8
# set the number of colors
N_COLOR = 4
# set the graph. define them which vertices are connected to each other
graph = [(0,1), (0,2), (1,2), (5,6), (2,3), (2,5), (3,4), (5,7), (7, 6)]

Formulation with PyQUBO

We import the required classes from PyQUBO.

[2]:
from pyqubo import Array, Constraint, solve_qubo

At First, we prepare variables for representing QUBO. We set an array of variables using Array. In this time, we need the number of (N_VER) x (N_COLOR), therefore we set shape argument as follows.

[3]:
x = Array.create('x', shape=(N_VER,N_COLOR), vartype='BINARY')
We could create binary variable ‘x’ with N_VER rows and N_COLOR columns.
Next, we set QUBO. PyQUBO allows us to follow the formula.
[4]:
# define first term (constraint)
H_A = Constraint(sum((1-sum(x[v,i] for i in range(1,N_COLOR)))**2 for v in range(N_VER)), label='HA')
# define seconde term (cost, objective function)
H_B = sum((-1+sum(x[v,i] for v in range (N_VER)))/2*(sum(x[v,i] for v in range (N_VER))) - sum(x[u,i]*x[v,i] for (u,v) in graph) for i in range (1,N_COLOR))
# set the entire Hamiltonian
Q = H_A+H_B

We can use Constraint function in the first term to make the script recognize that “this term is constraint”. The cost function is easily converted to QUBO (Python dictionary type) with Q.compile().to_qubo().

In OpenJij and D-Wave Ocean, QUBO is assumed to be represented by a Python dictionary type.

We can run it on each solver by .compile.

[5]:
# compile this model
model = Q.compile()
qubo, offset = model.to_qubo()

qubo is set to QUBO and offset is set to the constant that appears when it is converted to QUBO.

A Simulated Annealing solver solve_qubo(qubo) in PyQUBO is now deprecated. It is recommended to call D-Wave Ocean SDK dwave-neal directly.

[6]:
# use SA on neal
import neal
sampler = neal.SimulatedAnnealingSampler()
raw_solution = sampler.sample_qubo(qubo)
print(raw_solution)
  x[0][1] x[0][2] x[0][3] x[1][1] x[1][2] x[1][3] ... x[7][3] energy num_oc.
0       0       0       1       0       0       1 ...       0   -8.0       1
['BINARY', 1 rows, 1 samples, 24 variables]

.first.sample extracts the lowest energy of all derived solutions.

[7]:
raw_solution.first.sample
[7]:
{'x[0][1]': 0,
 'x[0][2]': 0,
 'x[0][3]': 1,
 'x[1][1]': 0,
 'x[1][2]': 0,
 'x[1][3]': 1,
 'x[2][1]': 0,
 'x[2][2]': 0,
 'x[2][3]': 1,
 'x[3][1]': 1,
 'x[3][2]': 0,
 'x[3][3]': 0,
 'x[4][1]': 1,
 'x[4][2]': 0,
 'x[4][3]': 0,
 'x[5][1]': 0,
 'x[5][2]': 1,
 'x[5][3]': 0,
 'x[6][1]': 0,
 'x[6][2]': 1,
 'x[6][3]': 0,
 'x[7][1]': 0,
 'x[7][2]': 1,
 'x[7][3]': 0}

Let’s look at the solutions obtained. We can see that is is stored in a dictionary type with the string as the key like ‘x[0][0]’: 1.

In this form, it is difficult to analyze from now on.

PyQUBO has the .decode_sample() function to convert it into a more manageable form.

[8]:
# decode (convert) result into a manageable form
decoded_sample = model.decode_sample(raw_solution.first.sample, vartype="BINARY")
# below is for visualization
# .array(variable, index) extracts the specific element of the index
x_solution = {}
for i in range(N_VER):
    x_solution[i] = {}
    for j in range(1,N_COLOR):
        x_solution[i][j] = decoded_sample.array('x', (i, j))
x_solution
[8]:
{0: {1: 0, 2: 0, 3: 1},
 1: {1: 0, 2: 0, 3: 1},
 2: {1: 0, 2: 0, 3: 1},
 3: {1: 1, 2: 0, 3: 0},
 4: {1: 1, 2: 0, 3: 0},
 5: {1: 0, 2: 1, 3: 0},
 6: {1: 0, 2: 1, 3: 0},
 7: {1: 0, 2: 1, 3: 0}}

We can see three groups, (0,1,2), (3, 4), (5,6,7).

This solution forms creeks on each of the graphs given this time.

.constraints(only_broken=True) shows how the penalty term is broken (when penalty is not equal to 0).

[9]:
print(decoded_sample.constraints(only_broken=True))
{}

In this time, We can see empty dictionary because constraint is satisfied.

decode function is very useful that can automatically check the constraints is satisfied or not.

Run with OpenJij

We just solved creek coverage problem in PyQUBO SA. Next, let’s use OpenJij.

Similarly, OpenJij can perform SA, but we use SQA, which is not implemented in PyQUBO.

[10]:
# import OpenJij
import openjij as oj

# solve this problem with SQA
sampler = oj.SQASampler()
# substitute into QUBO what we created using .to_qubo
response = sampler.sample_qubo(Q=qubo)

We can run with other algorithms and machines by replaceing sampler part.

Finally, we decode the result with OpenJij using PyQUBO decoder.

[11]:
# get the state of lowest energy
dict_solution = response.first.sample
# decode
decoded_sample = model.decode_sample(raw_solution.first.sample, vartype="BINARY")
# visualize
# .array(variable, index) extracts the specific element
x_solution = {}
for i in range(N_VER):
    x_solution[i] = {}
    for j in range(1,N_COLOR):
        x_solution[i][j] = decoded_sample.array('x', (i, j))
x_solution
[11]:
{0: {1: 0, 2: 0, 3: 1},
 1: {1: 0, 2: 0, 3: 1},
 2: {1: 0, 2: 0, 3: 1},
 3: {1: 1, 2: 0, 3: 0},
 4: {1: 1, 2: 0, 3: 0},
 5: {1: 0, 2: 1, 3: 0},
 6: {1: 0, 2: 1, 3: 0},
 7: {1: 0, 2: 1, 3: 0}}

Conclusion

We learned how to formulate it using PyQUBO and how it works with OpenJij.

Procedures are as follows.

  1. set up variables in pyqubo.Array

  2. formulate QUBO

  3. compile QUBO and convert it to a dictionary type

  4. solve optimization problems using solvers such as OpenJij that accepts dictionary type QUBOs.

  5. decode solution as a dictionary with the subscript as a key.

PyQUBO is useful and powerful tool for formulationg and evaluating constraints. When we use in conjunction with OpenJij, which provides a variety of solvers, it provides comfortable development experience.

Reference:PyQUBO official document https://pyqubo.readthedocs.io/en/latest/reference/array.html?highlight=arry%20create