Getting Started

Installation

If you use pip, just type

pip install pyqubo

You can install from the source code like

git clone https://github.com/recruit-communications/pyqubo.git
cd pyqubo
python setup.py install

QUBO and Ising Model

If you want to solve a combinatorial optimization problem by quantum or classical annealing machines, you need to represent your problem in QUBO (Quadratic Unconstrained Binary Optimization) or Ising model. PyQUBO converts your problem into QUBO or Ising model format.

The objective function of QUBO is defined as:

\[\sum_{i \le j}q_{ij}x_{i}x_{j}\]

where \(x_{i}\) represents a binary variable which takes 0 or 1, and \(q_{ij}\) represents a quadratic coefficient. Note that \(q_{ii}x_{i}x_{i}=q_{ii}x_{i}\), since \(x_{i}^2=x_{i}\). Thus, the above expression includes linear terms of \(x_{i}\).

The objective function of Ising model is defined as:

\[\sum_{i}h_{i}s_{i} + \sum_{i<j}J_{ij}s_{i}s_{j}\]

where \(s_{i}\) represents spin variable which takes -1 or 1, \(h_{i}\) represents an external magnetic field and \(J_{ij}\) represents an interaction between spin \(i\) and \(j\).

Basic Usage

With PyQUBO, you can construct QUBOs with 3 steps:

  1. Define the hamiltonian.
>>> from pyqubo import Spin
>>> s1, s2, s3, s4 = Spin("s1"), Spin("s2"), Spin("s3"), Spin("s4")
>>> H = (4*s1 + 2*s2 + 7*s3 + s4)**2
  1. Compile the hamiltonian to get a model.
>>> model = H.compile()
  1. Call ‘to_qubo()’ to get QUBO coefficients.
>>> qubo, offset = model.to_qubo()
>>> pprint(qubo) # doctest: +SKIP
{('s1', 's1'): -160.0,
 ('s1', 's2'): 64.0,
 ('s1', 's3'): 224.0,
 ('s1', 's4'): 32.0,
 ('s2', 's2'): -96.0,
 ('s2', 's3'): 112.0,
 ('s2', 's4'): 16.0,
 ('s3', 's3'): -196.0,
 ('s3', 's4'): 56.0,
 ('s4', 's4'): -52.0}
>>> print(offset)
196.0

In this example, you want to solve Number Partitioning Problem with a set S = {4, 2, 7, 1}. The hamiltonian \(H\) is represented as

\[H = (4 s_{1} + 2 s_{2} + 7 s_{3} + s_{4})^2\]

where \(s_{i}\) is a \(i\) th spin variable which indicates a group the \(i\) th number should belong to. In PyQUBO, spin variables are internally converted to binary variables via the relationship \(x_{i} = (s_{i}+1)/2\). The QUBO coefficents and the offset returned from Model.to_qubo() represents the following objective function:

\[\begin{split}&-160 x_{1}x_{1} + 64 x_{1}x_{2} + 224 x_{1}x_{3} + 32 x_{1}x_{4} - 96 x_{2}x_{2}\\ &+ 112 x_{2}x_{3} + 16 x_{2}x_{4} - 196 x_{3}x_{3} + 56 x_{3}x_{4} - 52 x_{4}x_{4} + 196\end{split}\]
  1. Call ‘to_ising()’ to get Ising coefficients.

If you want to get the coefficient of the Ising model, just call Model.to_ising() method like below.

>>> linear, quadratic, offset = model.to_ising()
>>> pprint(linear) # doctest: +SKIP
{'s1': 0.0, 's2': 0.0, 's3': 0.0, 's4': 0.0}
>>> pprint(quadratic) # doctest: +SKIP
{('s1', 's2'): 16.0,
 ('s1', 's3'): 56.0,
 ('s1', 's4'): 8.0,
 ('s2', 's3'): 28.0,
 ('s2', 's4'): 4.0,
 ('s3', 's4'): 14.0}
>>> print(offset)
70.0

where linear represents external magnetic fields \(h\), quadratic represents interactions \(J\) and offset represents the constant value in the objective function below.

\[16 s_{1}s_{2} + 56 s_{1}s_{3} + 8 s_{1}s_{4} + 28 s_{2}s_{3} + 4 s_{2}s_{4} + 14 s_{3}s_{4} + 70\]

Variable: Binary and Spin

When you define a hamiltonian, you can use Binary or Spin as a variable.

Example: If you want to define a hamiltonian with binary variables \(x \in \{0, 1\}\), use Binary.

>>> from pyqubo import Binary
>>> x1, x2 = Binary('x1'), Binary('x2')
>>> exp = 2
>>> from pyqubo import Binary
>>> x1, x2 = Binary('x1'), Binary('x2')
>>> exp = 2
>>> from pyqubo import Binary
>>> x1, x2 = Binary('x1'), Binary('x2')
>>> exp = 2*x1*x2 + 3*x1
>>> pprint(exp.compile().to_qubo()) # doctest: +SKIP
({('x1', 'x1'): 3.0, ('x1', 'x2'): 2.0, ('x2', 'x2'): 0.0}, 0.0)

Example: If you want to define a hamiltonian with spin variables \(s \in \{-1, 1\}\), use Spin.

>>> from pyqubo import Spin
>>> s1, s2 = Spin('s1'), Spin('s2')
>>> exp = 2*s1*s2 + 3*s1
>>> pprint(exp.compile().to_qubo()) # doctest: +SKIP
({('s1', 's1'): 2.0, ('s1', 's2'): 8.0, ('s2', 's2'): -4.0}, -1.0)

Array of Variables

Array class represents a multi-dimensional array of Binary or Spin.

Example: You can access each element of the matrix with an index like:

>>> from pyqubo import Array
>>> x = Array.create('x', shape=(2, 3), vartype='BINARY')
>>> x[0, 1] + x[1, 2]
(Binary(x[0][1])+Binary(x[1][2]))

Example: You can use Array to represent multiple spins in the example of partitioning problem above.

>>> from pyqubo import Array
>>> numbers = [4, 2, 7, 1]
>>> s = Array.create('s', shape=4, vartype='SPIN')
>>> H = sum(n * s for s, n in zip(s, numbers))**2
>>> model = H.compile()
>>> qubo, offset = model.to_qubo()
>>> pprint(qubo) # doctest: +SKIP
{('s[0]', 's[0]'): -160.0,
 ('s[0]', 's[1]'): 64.0,
 ('s[0]', 's[2]'): 224.0,
 ('s[0]', 's[3]'): 32.0,
 ('s[1]', 's[1]'): -96.0,
 ('s[1]', 's[2]'): 112.0,
 ('s[1]', 's[3]'): 16.0,
 ('s[2]', 's[2]'): -196.0,
 ('s[2]', 's[3]'): 56.0,
 ('s[3]', 's[3]'): -52.0}

Placeholder

If you have a parameter that you will probably update, such as the strength of the constraints in your hamiltonian, using Placeholder will save your time. If you define the parameter by Placeholder, you can specify the value of the parameter after compile. This means that you don’t have to compile repeatedly for getting QUBOs with various parameter values. It takes longer time to execute a compile when the problem size is bigger. In that case, you can save your time by using Placeholder.

Example: If you have an objective function \(2a+b\), and a constraint \(a+b=1\) whose hamiltonian is \((a+b-1)^2\) where \(a,b\) is qbit variable, you need to find the penalty strength \(M\) such that the constraint is satisfied. Thus, you need to create QUBO with different values of \(M\). In this example, we create QUBO with \(M=5.0\) and \(M=6.0\).

In the first code, we don’t use placeholder. In this case, you need to compile the hamiltonian twice to get a QUBO with \(M=5.0\) and \(M=6.0\).

>>> from pyqubo import Binary
>>> a, b = Binary('a'), Binary('b')
>>> M = 5.0
>>> H = 2
>>> from pyqubo import Binary
>>> a, b = Binary('a'), Binary('b')
>>> M = 5.0
>>> H = 2
>>> from pyqubo import Binary
>>> a, b = Binary('a'), Binary('b')
>>> M = 5.0
>>> H = 2*a + b + M*(a+b-1)**2
>>> model = H.compile()
>>> qubo, offset = model.to_qubo() # QUBO with M=5.0
>>> M = 6.0
>>> H = 2*a + b + M*(a+b-1)**2
>>> model = H.compile()
>>> qubo, offset = model.to_qubo() # QUBO with M=6.0

If you don’t want to compile twice, define \(M\) by Placeholder.

>>> from pyqubo import Placeholder
>>> a, b = Binary('a'), Binary('b')
>>> M = Placeholder('M')
>>> H = 2*a + b + M*(a+b-1)**2
>>> model = H.compile()
>>> qubo, offset = model.to_qubo(feed_dict={'M': 5.0})

You get a QUBO with different value of M without compile

>>> qubo, offset = model.to_qubo(feed_dict={'M': 6.0})

The actual value of the placeholder \(M\) is specified in calling Model.to_qubo() as a value of the feed_dict.

Decode Solution

When you get a solution from the Ising solver, Model.decode_solution() decodes the solution and returns decoded_solution in dictionary form.

Example: You are solving a partitioning problem.

>>> from pyqubo import Array
>>> numbers = [4, 2, 7, 1]
>>> s = Array.create('s', 4, 'SPIN')
>>> H = sum(n * s_i for s_i, n in zip(s, numbers))**2
>>> model = H.compile()
>>> qubo, offset = model.to_qubo()

Let’s assume that you get a solution {'s[0]': 0, 's[1]': 0, 's[2]': 1, 's[3]': 0} from the solver.

>>> raw_solution = {'s[0]': 0, 's[1]': 0, 's[2]': 1, 's[3]': 0} # solution from the solver
>>> decoded_solution, broken, energy = model.decode_solution(raw_solution, vartype='BINARY')
>>> pprint(decoded_solution)
{'s': {0: 0, 1: 0, 2: 1, 3: 0}}
>>> broken
{}
>>> energy
0.0

You can see that decoded_solution has the decoded solution of spin vector where \(i\) th element of the vector is accessed via s[i]. broken represents broken constraint which will be explained in the following section. energy represents energy of the problem.

Validation of Constraints

When the hamiltonian has constraints, you can let the compiler recognize the hamiltonian of the constraint with Constraint. When you decode the solution, the model let you know which constraints are broken. You don’t have to write additional programs for validation of the constraints.

Example: If you have an objective function \(2a+b\), and a constraint \(a+b=1\) whose hamiltonian is \((a+b-1)^2\) where \(a,b\) is qbit variable, you need to put \((a+b-1)^2\) in Constraint to tell the compiler that this hamiltonian is constraint i.e. it should be zero when the solution is not broken.

>>> from pyqubo import Binary, Constraint
>>> a, b = Binary('a'), Binary('b')
>>> M = 5.0 # strength of the constraint
>>> H = 2*a + b + M * Constraint((a+b-1)**2, label='a+b=1')
>>> model = H.compile()

Let’s assume that you get a solution

>>> from pyqubo import Binary, Constraint
>>> a, b = Binary('a'), Binary('b')
>>> M = 5.0 # strength of the constraint
>>> H = 2*a + b + M * Constraint((a+b-1)**2, label='a+b=1')
>>> model = H.compile()

Let’s assume that you get a solution

>>> from pyqubo import Binary, Constraint
>>> a, b = Binary('a'), Binary('b')
>>> M = 5.0 # strength of the constraint
>>> H = 2*a + b + M * Constraint((a+b-1)**2, label='a+b=1')
>>> model = H.compile()

Let’s assume that you get a solution {'a': 1, 'b': 1} from the solver which breaks the constraint \(a+b=1\).

>>> raw_solution = {'a': 1, 'b': 1}
>>> decoded_solution, broken, energy = model.decode_solution(raw_solution, vartype='BINARY')
>>> pprint(broken)
{'a+b=1': {'penalty': 1.0, 'result': {'a': 1, 'b': 1}}}

broken object contains the information about the broken constraint. If no constraint is broken, broken is empty.