1-An Introduction to OpenJij

Open in Colab

OpenJij is a heuristic optimization library of the Ising model and QUBO. It has a Python interface, therefore it can be easily written in Python, although the core part of the optimization is implemented with C++.

Let’s install OpenJij by pip.
(If you’re using Colab, you’ll need to install cmake first $ !pip install cmake.)
[30]:
# !pip install openjij

The Ising model

The Ising model is a model treated in physics and is written as follows.

\[H(\{\sigma_i\}) = \sum_{i > j} J_{ij}\sigma_i \sigma_j + \sum_{i=1}^N h_i \sigma_i\]
\[\sigma_i \in \{-1, 1\}, i=1,\cdots N\]
\(H(\{\sigma_i\})\) is called Hamiltonian, it is like an energy or a cost function.
\(\sigma_i\) is a binary variable it takes only two values \(\{1, -1\}\).

\(\sigma_i\) corresponds to the physical quantity called spin in physics, it is sometimes called a spin variable or spin simply. Spin seems a micro magnet. The value of this variable corresponds to the physical state (the orientation of the magnet). \(\sigma_i = -1\) means the magnet is up and 1 means magnet is down.

\(H\) depends on the combination of variables \(\{\sigma_i\} = \{\sigma_1, \sigma_2, \cdots, \sigma_N\}\).
\(J_{ij}, h_i\) represents the problem, and they are called interaction coefficients, longitudinal magnetic fields respectively.

OpenJij is a library which searches a combination of spin variables \(\{\sigma_i\}\) that minimizes \(H(\{\sigma_i\})\) when given \(J_{ij}\) and \(h_i\).

Let me show you a example.

Solve a problem with OpenJij

We assume that the number of variables is \(N=5\), and the longitudinal field and interaction are

\[h_i = -1~\text{for} ~\forall i, ~ J_{ij} = -1~\text{for} ~\forall i, j\]

respectively.

In this case, since all interaction coefficients are negative, the energy is lower when each spin variable takes the same value. Also, since all longitudinal magnetic fields are negative, the energy is lower when each spin takes 1. Therefore, we can predict the answer of this problem is \(\{\sigma_i\} = \{1, 1, 1, 1, 1\}\).

Now let’s solve the Ising model using OpenJij.

[1]:
import openjij as oj

# Create interaction coefficients and longitudinal magnetic fields that represent the problem.
# OpenJij accepts problems in dictionary type.

N = 5
h = {i: -1 for i in range(N)}
J = {(i, j): -1 for i in range(N) for j in range(i+1, N)}

print('h_i: ', h)
print('Jij: ', J)
h_i:  {0: -1, 1: -1, 2: -1, 3: -1, 4: -1}
Jij:  {(0, 1): -1, (0, 2): -1, (0, 3): -1, (0, 4): -1, (1, 2): -1, (1, 3): -1, (1, 4): -1, (2, 3): -1, (2, 4): -1, (3, 4): -1}
[2]:
# First, create an instance of Sampler to solve the problem.
sampler = oj.SASampler(num_reads=1)
# Solve the problem(h, J) by using the method of sampler.
response = sampler.sample_ising(h, J)

# The results (states) of calculation is in result.states.
print(response.states)

# If you want to see the result with subscript, use response.samples.
print([s for s in response.samples()])
[[1 1 1 1 1]]
[{0: 1, 1: 1, 2: 1, 3: 1, 4: 1}]

Explanation of OpenJij

Then, I will explain the code showed above. OpenJij has two interfaces currently. One of them used above is the same interface as D-Wave Ocean. Therefore, if you get used to using OpenJij, easy changing to Ocean.

Another interface will not be explained here, but it can be easily extended by using the structure of OpenJij (graph, updater, and algorithm) directly. In the present situation, however, it would be sufficient to be able to use the interface which treated in the cell above.

Sampler

In the above, we created an instance of Sampler after creating the problem in the dictionary type.

sampler = oj.SASampler(num_reads=1)

Here, this Sampler is the object of selecting what algorithm and hardware you use. When you want to try the other algorithm, change this Sampler.

The algorithm treated in OpenJij is a heuristic and stochastic algorithm, it may return the different solutions each time the problem is solved, so it may not always obtain the optimal solution. Therefore, we solve the problem many times and search for a good solution. So we call Sampler which means that we sample the solutions.

SASampler treated in the above cell is a sampler that samples the solution by using an algorithm called Simulated annealing. For example, we can choose other samples below.

  • SQASampler: An algorithm for simulating quantum annealing called Simulated Quantum Annealing (SQA) with a classical computer.

  • GPUSQASampler: It runs SQA on GPU. It can treat only a special structure called Chimera graph.

sample_ising(h, J)

When you want to solve the Ising model, you can use .sample_ising(h, J). As I will explain later, when you treat QUBO which equivalent to the Ising model optimization, use .sample_qubo(Q).

Response

.sample_ising(h, J) returns the Response class. This class contains results of solutions and energies that Sampler solved.

  • .states:

    • type: list (list (int))

    • This list contains solutions to the number of iteration. > In physics, an array of spins called a state. Since .states stores the solutions to the number of iteration, we defined as .states from the meaning that multiple states are stored.

  • .energies:

    • type: list (float)

    • This list contains the energy of each solution to the number of iteration.

The Response class inherits the SampleSet class of D-Wave’s dimod. More detailed information can be found at the following links SampleSet: dimod documantation.

Let’s look at the acutual code.

[3]:
# keys of h, J dictionaries can treat not only numbers.
h = {'a': -1, 'b': -1}
J = {('a', 'b'): -1, ('b', 'c'): 1}
# # Try solving 10 times by SA at a time. With the argument called num_reads.
sampler = oj.SASampler(num_reads=10)
response = sampler.sample_ising(h, J)
print(response.first.sample)
print(response.first.energy)
{'a': 1, 'b': 1, 'c': -1}
-4.0
[4]:
print(response.states)
[[ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]
 [ 1  1 -1]]

You can see that the response.states contains 10 solutions when you open it. Since this problem is easy, returned the same answer is [1, -1, 1] every iterations.

The parameters passed in the constructor such as num_reads can be overwritten by the method that executes sampling such as .sample_ising.

response = sampler.sample_ising(h, J, num_reads=2)
response.states
> [[1, 1, -1],[1, 1, -1]]
[41]:
# Next, let's look the energies.
response.energies
[41]:
array([-4., -4., -4., -4., -4., -4., -4., -4., -4., -4.])

Of course, the energies also take the same value every iterations. The solutions in response.states are list type, so you confuse the corresponding to the strings "a", "b", "c" when you set the problem. That is in response.indices.

[5]:
response.indices
[5]:
['a', 'b', 'c']

.first is useful if you only want to know the state with the lowest energy value.

[6]:
response.first
[6]:
Sample(sample={'a': 1, 'b': 1, 'c': -1}, energy=-4.0, num_occurrences=1)

Let’s try to solve QUBO

When you want to solve social problems, it is often easier to formulate as QUBO(Quadratic Unconstraited Binary Optimization, ) than the Ising model. Basically, it returns the same solution when using the Ising model.

Also, QUBO is often called as UBQP (Unconstraited binary quadratic optimization problems).

QUBO is written as follows.

\[H(\{q_i\}) = \sum_{i\geq j} Q_{ij}q_i q_j\]
\[q_i \in \{0, 1\}\]

Binary variables take 0 or 1, that is difference between QUBO and the Ising model.

There are other ways to define \(\sum\) and \(Q_{ij}\) (for example, make \(Q_{ij}\) a symmetric matrix), but in this time we have formulated it as above.

The Ising model and QUBO are equivalent for the reason why they be converted to each other. You can convert these equation: \(q_i=(\sigma_i + 1)/2\)

In the case of QUBO, you search a combination of \(\{q_i\}\) that minimizes \(H(\{q_i\})\) under the \(Q_{ij}\) which is given by a problem.

Also, since \(q_i^2 = q_i\) (because \(q_i\) takes only 0 or 1), you can write as fllows

\[H(\{q_i\}) = \sum_{i > j} Q_{ij}q_i q_j + \sum_i Q_{ii} q_i\]

Namely, the diagonal components of \(Q_{ij}\) is corresponds to the coefficients of the first-order term of \(q_i\).

Let’s solve QUBO with OpenJij

[44]:
# Q_ij (dict)
Q = {(0, 0): -1, (0, 1): -1, (1, 2): 1, (2, 2): 1}
sampler = oj.SASampler()
# use .sample_qubo to solve QUBO
response = sampler.sample_qubo(Q)
print(response.states)
[[1 1 0]]

In QUBO, the solutions are made by 0 or 1 because the variables take 0 or 1.

Let’s solve a little difficult problem

You want to solve a problem of QUBO which has \(Q_{ij}\) given by 50 random variables.

[15]:
N = 100
# Create Q_ij randomly
import random
Q = {(i, j): random.uniform(-1, 1) for i in range(N) for j in range(i+1, N)}

# Solve by OpenJij
sampler = oj.SASampler()
response = sampler.sample_qubo(Q, num_reads=100)
[16]:
# Look at some energies.
response.energies[:5]
[16]:
array([-192.6801102, -192.6801102, -192.6801102, -192.6801102,
       -192.6801102])

You will realize that the energies takes a different values than before. If you make \(Q_{ij}\) randomly, the problem becomes difficult generally. Therefore, SASampler does not necessarily give the same solution every time. Now let’s look at the histogram of energies.

[17]:
import matplotlib.pyplot as plt
plt.hist(response.energies, bins=15)
plt.xlabel('Energy', fontsize=15)
plt.ylabel('Frequency', fontsize=15)
plt.show()
../_images/en_001-Introduction_27_0.svg

The lower the energy, the better, but you can find that the high energy solutions appear sometimes. However, almost half of them seem to have the lowest energy.

Openjij solved the problem 100 times. Since the optimal solution is the lowest energy, the lowest energy solution among the solved (sampled) solutions should be close to the optimal solution. Therefore, we search for the solution from .states.

Caution: SA can not always give an optimal solution. Therefore, there is no guarantee that the solution with the lowest energy in this time is an optimal solution. It is an approximate solution.

[18]:
import numpy as np

min_samples = response.first
min_samples
[18]:
Sample(sample={0: 1, 1: 1, 2: 1, 3: 1, 4: 0, 5: 1, 6: 0, 7: 1, 8: 1, 9: 1, 10: 1, 11: 1, 12: 1, 13: 1, 14: 1, 15: 0, 16: 1, 17: 1, 18: 1, 19: 1, 20: 1, 21: 0, 22: 1, 23: 1, 24: 0, 25: 1, 26: 0, 27: 1, 28: 1, 29: 1, 30: 1, 31: 1, 32: 0, 33: 0, 34: 1, 35: 1, 36: 1, 37: 1, 38: 1, 39: 0, 40: 0, 41: 1, 42: 1, 43: 1, 44: 1, 45: 1, 46: 1, 47: 1, 48: 1, 49: 1, 50: 1, 51: 0, 52: 1, 53: 1, 54: 1, 55: 1, 56: 0, 57: 1, 58: 1, 59: 1, 60: 1, 61: 1, 62: 0, 63: 0, 64: 1, 65: 1, 66: 1, 67: 0, 68: 0, 69: 1, 70: 0, 71: 1, 72: 1, 73: 1, 74: 1, 75: 1, 76: 0, 77: 0, 78: 0, 79: 1, 80: 0, 81: 1, 82: 1, 83: 1, 84: 0, 85: 0, 86: 0, 87: 0, 88: 0, 89: 0, 90: 1, 91: 1, 92: 1, 93: 0, 94: 1, 95: 1, 96: 1, 97: 0, 98: 1, 99: 1}, energy=-192.6801102019457, num_occurrences=1)

Now you get the lowest energy solution. The solution which is contained in this min_states is the approximate solution in this time. Contained solutions have the same energy, but the arrangement of spins may differ in some cases (In the case of “degenerated state”).

So, after choosing the lowest energy state as above, choose the spin orientation from min_state as you like. Then, we solved the problem approximately.

Next section, “2-Evaluation”, we will explain the index(“time to solution” and “the residual energy” etc.) for evaluate the solution.

[ ]: