Permutation matrix generation
Many combinatorial optimization problems are permutation-based in the sense that the objective is to find an optimal permutation. As a fundamental technique for formulating such optimization problems, a matrix of binary variables is used in their QUBO formulation.
Permutation matrix
Let $X=(x_{i,j})$ ($0\leq i,j\leq n-1$) is a matrix of $n\times n$ binary values. The matrix $X$ is called a permutation matrix if and only if every row and every column has exactly one entry equal to 1.
A permutation matrix represents a permutation of $n$ numbers $(0,1,\ldots,n-1)$, where $x_{i,j} = 1$ if and only if the $i$-th element is $j$.
QUBO formulation for permutation matrices
A binary variable matrix $X=(x_{i,j})$ ($0\leq i,j\leq n-1$) stores a permutation matrix if and only if the sum of each row and each column is 1. Thus, the following QUBO function takes the minimum value 0 if and only if $X$ stores a permutation matrix:
\[\begin{aligned} f(X) &= \sum_{i=0}^{n-1}\left(1-\sum_{j=0}^{n-1}x_{i,j}\right)^2+\sum_{j=0}^{n-1}\left(1-\sum_{i=0}^{n-1}x_{i,j}\right)^2 \end{aligned}\]PyQBPP program for generating permutation matrices
We can design a PyQBPP program based on the formula $f(X)$ above as follows:
import pyqbpp as qbpp
x = qbpp.var("x", 4, 4)
f = qbpp.expr()
for i in range(4):
s = qbpp.expr()
for j in range(4):
s += x[i][j]
f += qbpp.sqr(1 - s)
for j in range(4):
s = qbpp.expr()
for i in range(4):
s += x[i][j]
f += qbpp.sqr(1 - s)
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
result = solver.search({"best_energy_sols": 0})
for k, sol in enumerate(result.sols()):
row = [sol.get_vector(x[i]) for i in range(4)]
print(f"Solution {k} : {row}")
In this program, var("x", 4, 4) returns a nested Array of size $4\times 4$ named x. For an Expr object f, two double for-loops build the formula for $f(X)$. Using the Exhaustive Solver, all optimal solutions are computed and stored in sols. All solutions in sols are displayed one-by-one using sol.get_vector(). This program outputs all 24 permutations.
QUBO formulation using array functions and operations
Using vector_sum(), we can compute the row-wise and column-wise sums of a matrix x of binary variables:
vector_sum(x, 1): Computes the sum of each row ofxand returns an array of sizencontaining these sums.vector_sum(x, 0): Computes the sum of each column ofxand returns an array of sizencontaining these sums.
For these two arrays of size n, sqr() squares each element, and sum() computes the sum of all elements.
The following program implements a QUBO formulation using these array functions and operations:
import pyqbpp as qbpp
x = qbpp.var("x", 4, 4)
f = qbpp.sum(qbpp.sqr(qbpp.vector_sum(x, 1) - 1)) + qbpp.sum(qbpp.sqr(qbpp.vector_sum(x, 0) - 1))
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
result = solver.search({"best_energy_sols": 0})
for k, sol in enumerate(result.sols()):
perm = []
for i in range(4):
for j in range(4):
if sol(x[i][j]) == 1:
perm.append(j)
print(f"Solution {k}: {perm}")
Assignment problem and its QUBO formulation
Let $C = (c_{i,j})$ be a cost matrix of size $n \times n$. The assignment problem for $C$ is to find a permutation $p:\lbrace 0,1,\ldots, n-1\rbrace \rightarrow \lbrace 0,1,\ldots, n-1\rbrace$ that minimizes the total cost:
\[\begin{aligned} g(p) &= \sum_{i=0}^{n-1}c_{i,p(i)} \end{aligned}\]We can use a permutation matrix $X = (x_{i,j})$ of size $n \times n$ for a QUBO formulation of this problem by defining:
\[\begin{aligned} g(X) &= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}x_{i,j} \end{aligned}\]We combine the permutation constraint $f(X)$ and the total cost $g(X)$:
\[\begin{aligned} h(X) &= P\cdot f(X)+g(X) \end{aligned}\]Here, $P$ is a sufficiently large positive constant that prioritizes the permutation constraints.
PyQBPP program for the assignment problem
In this program, the cost matrix c is defined as a 2D Array using qbpp.Array() with a nested Python list. qbpp.Array() automatically converts nested lists into nested Array objects, so multi-dimensional arrays can be created concisely. The element-wise product c * x then computes $c_{i,j} \cdot x_{i,j}$ for all entries.
import pyqbpp as qbpp
c = qbpp.Array([[58, 73, 91, 44],
[62, 15, 87, 39],
[78, 56, 23, 94],
[11, 85, 68, 72]])
x = qbpp.var("x", 4, 4)
f = qbpp.sum(qbpp.vector_sum(x, 1) == 1) + qbpp.sum(qbpp.vector_sum(x, 0) == 1)
g = qbpp.sum(c * x)
h = 1000 * f + g
h.simplify_as_binary()
solver = qbpp.EasySolver(h)
sol = solver.search({"time_limit": 1.0})
print("sol =", sol)
result = []
for i in range(4):
for j in range(4):
if sol(x[i][j]) == 1:
result.append(j)
print("Result :", result)
for i in range(len(result)):
print(f"c[{i}][{result[i]}] = {c[i][result[i]]}")
We use the Easy Solver to find a solution of h. The time limit for searching is set to 1.0 seconds by passing a parameter dictionary to search(). The output of this program is as follows:
Result : [3, 1, 2, 0]
c[0][3] = 44
c[1][1] = 15
c[2][2] = 23
c[3][0] = 11
NOTE For an expression
fand an integerm,f == mreturns an expressionsqr(f - m), which takes the minimum value 0 if and only if the equalityf == mis satisfied.
置換行列の生成
多くの組合せ最適化問題は、最適な置換を見つけることが目的であるという意味で、置換に基づいています。 このような最適化問題を定式化するための基本的な手法として、QUBO定式化ではバイナリ変数の行列が使用されます。
置換行列
$X=(x_{i,j})$ ($0\leq i,j\leq n-1$) を $n\times n$ のバイナリ値の行列とします。 行列 $X$ が置換行列であるとは、すべての行とすべての列にちょうど1つの1のエントリがある場合に限ります。
置換行列は $n$ 個の数 $(0,1,\ldots,n-1)$ の置換を表し、$x_{i,j} = 1$ であることと $i$ 番目の要素が $j$ であることが同値です。
置換行列のQUBO定式化
バイナリ変数行列 $X=(x_{i,j})$ ($0\leq i,j\leq n-1$) が置換行列を格納しているのは、各行と各列の合計が1である場合に限ります。 したがって、以下のQUBO関数は $X$ が置換行列を格納している場合に限り最小値0をとります:
\[\begin{aligned} f(X) &= \sum_{i=0}^{n-1}\left(1-\sum_{j=0}^{n-1}x_{i,j}\right)^2+\sum_{j=0}^{n-1}\left(1-\sum_{i=0}^{n-1}x_{i,j}\right)^2 \end{aligned}\]置換行列を生成するPyQBPPプログラム
上記の式 $f(X)$ に基づいて、以下のようにPyQBPPプログラムを設計できます:
import pyqbpp as qbpp
x = qbpp.var("x", 4, 4)
f = qbpp.expr()
for i in range(4):
s = qbpp.expr()
for j in range(4):
s += x[i][j]
f += qbpp.sqr(1 - s)
for j in range(4):
s = qbpp.expr()
for i in range(4):
s += x[i][j]
f += qbpp.sqr(1 - s)
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
result = solver.search({"best_energy_sols": 0})
for k, sol in enumerate(result.sols()):
row = [sol.get_vector(x[i]) for i in range(4)]
print(f"Solution {k} : {row}")
このプログラムでは、var("x", 4, 4) が x という名前の $4\times 4$ サイズのネストされた Array を返します。 Expr オブジェクト f に対して、2つの二重forループが $f(X)$ の式を構築します。 Exhaustive Solverを使用して、すべての最適解が計算され sols に格納されます。 sols 内のすべての解は sol.get_vector() を使用して1つずつ表示されます。 このプログラムは24個すべての置換を出力します。
配列関数と演算を使用したQUBO定式化
vector_sum() を使用して、バイナリ変数の行列 x の行方向および列方向の合計を計算できます:
vector_sum(x, 1):xの各行の合計を計算し、これらの合計を含むサイズnの配列を返します。vector_sum(x, 0):xの各列の合計を計算し、これらの合計を含むサイズnの配列を返します。
これら2つのサイズ n の配列に対して、sqr() は各要素を二乗し、sum() はすべての要素の合計を計算します。
以下のプログラムは、これらの配列関数と演算を使用してQUBO定式化を実装しています:
import pyqbpp as qbpp
x = qbpp.var("x", 4, 4)
f = qbpp.sum(qbpp.sqr(qbpp.vector_sum(x, 1) - 1)) + qbpp.sum(qbpp.sqr(qbpp.vector_sum(x, 0) - 1))
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
result = solver.search({"best_energy_sols": 0})
for k, sol in enumerate(result.sols()):
perm = []
for i in range(4):
for j in range(4):
if sol(x[i][j]) == 1:
perm.append(j)
print(f"Solution {k}: {perm}")
割当問題とそのQUBO定式化
$C = (c_{i,j})$ をサイズ $n \times n$ のコスト行列とします。 $C$ に対する割当問題は、総コストを最小化する置換 $p:\lbrace 0,1,\ldots, n-1\rbrace \rightarrow \lbrace 0,1,\ldots, n-1\rbrace$ を見つける問題です:
\[\begin{aligned} g(p) &= \sum_{i=0}^{n-1}c_{i,p(i)} \end{aligned}\]この問題のQUBO定式化には、サイズ $n \times n$ の置換行列 $X = (x_{i,j})$ を使用し、以下のように定義します:
\[\begin{aligned} g(X) &= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}x_{i,j} \end{aligned}\]置換制約 $f(X)$ と総コスト $g(X)$ を組み合わせます:
\[\begin{aligned} h(X) &= P\cdot f(X)+g(X) \end{aligned}\]ここで、$P$ は置換制約を優先するための十分に大きな正の定数です。
割当問題のPyQBPPプログラム
このプログラムでは、コスト行列 c をネストされたPythonリストから qbpp.Array() で2次元の Array として定義しています。 qbpp.Array() はネストされたリストを自動的にネストされた Array オブジェクトに変換するため、多次元配列を簡潔に作成できます。 要素ごとの積 c * x は全要素について $c_{i,j} \cdot x_{i,j}$ を計算します。
import pyqbpp as qbpp
c = qbpp.Array([[58, 73, 91, 44],
[62, 15, 87, 39],
[78, 56, 23, 94],
[11, 85, 68, 72]])
x = qbpp.var("x", 4, 4)
f = qbpp.sum(qbpp.vector_sum(x, 1) == 1) + qbpp.sum(qbpp.vector_sum(x, 0) == 1)
g = qbpp.sum(c * x)
h = 1000 * f + g
h.simplify_as_binary()
solver = qbpp.EasySolver(h)
sol = solver.search({"time_limit": 1.0})
print("sol =", sol)
result = []
for i in range(4):
for j in range(4):
if sol(x[i][j]) == 1:
result.append(j)
print("Result :", result)
for i in range(len(result)):
print(f"c[{i}][{result[i]}] = {c[i][result[i]]}")
Easy Solverを使用して h の解を求めます。 search() にパラメータ辞書を渡すことで、探索の制限時間を1.0秒に設定しています。 このプログラムの出力は以下の通りです:
Result : [3, 1, 2, 0]
c[0][3] = 44
c[1][1] = 15
c[2][2] = 23
c[3][0] = 11
NOTE 式
fと整数mに対して、f == mは式sqr(f - m)を返します。 これは等式f == mが満たされる場合に限り最小値0をとります。