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)
sols = solver.search_optimal_solutions()
for k, sol in enumerate(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 Vector 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 vector 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 a vector of sizencontaining these sums.vector_sum(x, 0): Computes the sum of each column ofxand returns a vector of sizencontaining these sums.
For these two vectors of size n, sqr() squares each element, and sum() computes the sum of all elements.
The following program implements a QUBO formulation using these vector 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)
sols = solver.search_optimal_solutions()
for k, sol in enumerate(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 Vector using qbpp.Vector() with a nested Python list. qbpp.Vector() automatically converts nested lists into nested Vector 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.Vector([[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)
solver.time_limit(1.0)
sol = solver.search()
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 calling the time_limit() method. 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)
sols = solver.search_optimal_solutions()
for k, sol in enumerate(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$ サイズのネストされた Vector を返します。 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)
sols = solver.search_optimal_solutions()
for k, sol in enumerate(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.Vector() で2次元の Vector として定義しています。 qbpp.Vector() はネストされたリストを自動的にネストされた Vector オブジェクトに変換するため、多次元配列を簡潔に作成できます。 要素ごとの積 c * x は全要素について $c_{i,j} \cdot x_{i,j}$ を計算します。
import pyqbpp as qbpp
c = qbpp.Vector([[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)
solver.time_limit(1.0)
sol = solver.search()
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 の解を求めます。 time_limit() メソッドを呼び出すことで、探索の制限時間を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をとります。