順列行列生成
多くの組み合わせ最適化問題は、最適置換を見つけることが目的であるという意味で、置換ベースです。 このような最適化問題を定式化するための基本的な手法として、その QUBO 定式化では、2 値変数の行列が使用されます。
順列行列
$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$ である場合に限る。 例えば、上記の置換行列は置換 $(1,2,0,3)$ を表す。
置換行列の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}\]置換行列を生成するQUBO++プログラム
上記の式 $f(X)$ に基づいて、次のように QUBO++ プログラムを設計できる:
#include "qbpp.hpp"
#include "qbpp_exhaustive_solver.hpp"
int main() {
auto x = qbpp::var("x", 4, 4);
auto f = qbpp::expr();
for (size_t i = 0; i < 4; i++) {
auto s = qbpp::expr();
for (size_t j = 0; j < 4; j++) {
s += x[i][j];
}
f += qbpp::sqr(1 - s);
}
for (size_t j = 0; j < 4; j++) {
auto s = qbpp::expr();
for (size_t i = 0; i < 4; i++) {
s += x[i][j];
}
f += qbpp::sqr(1 - s);
}
f.simplify_as_binary();
auto solver = qbpp::exhaustive_solver::ExhaustiveSolver(f);
auto sols = solver.search_optimal_solutions();
for (size_t k = 0; k < sols.size(); k++) {
const auto& sol = sols[k];
std::cout << "Solution " << k << " : " << sol(x) << std::endl;
}
}
このプログラムでは、qbpp::var("x",4,4) qbpp::Var オブジェクト
を返すx。
オブジェクト qbpp::Expr オブジェクトfに対して、2重のforループが
$f(X)$の計算式を追加する。
網羅的ソルバーを用いて、全ての最適解が計算され
に格納されるsols。
の全ての解が sols は一つずつ表示される。
ここで、 sol(x) は値の行列を返す x を返します。このプログラムは、以下の24の順列をすべて出力します： solの値の行列を返す。
このプログラムは全24の順列を以下のように出力する：
Solution 0 : {{0,0,0,1},{0,0,1,0},{0,1,0,0},{1,0,0,0}}
Solution 1 : {{0,0,0,1},{0,0,1,0},{1,0,0,0},{0,1,0,0}}
Solution 2 : {{0,0,0,1},{0,1,0,0},{0,0,1,0},{1,0,0,0}}
Solution 3 : {{0,0,0,1},{0,1,0,0},{1,0,0,0},{0,0,1,0}}
Solution 4 : {{0,0,0,1},{1,0,0,0},{0,0,1,0},{0,1,0,0}}
Solution 5 : {{0,0,0,1},{1,0,0,0},{0,1,0,0},{0,0,1,0}}
Solution 6 : {{0,0,1,0},{0,0,0,1},{0,1,0,0},{1,0,0,0}}
Solution 7 : {{0,0,1,0},{0,0,0,1},{1,0,0,0},{0,1,0,0}}
Solution 8 : {{0,0,1,0},{0,1,0,0},{0,0,0,1},{1,0,0,0}}
Solution 9 : {{0,0,1,0},{0,1,0,0},{1,0,0,0},{0,0,0,1}}
Solution 10 : {{0,0,1,0},{1,0,0,0},{0,0,0,1},{0,1,0,0}}
Solution 11 : {{0,0,1,0},{1,0,0,0},{0,1,0,0},{0,0,0,1}}
Solution 12 : {{0,1,0,0},{0,0,0,1},{0,0,1,0},{1,0,0,0}}
Solution 13 : {{0,1,0,0},{0,0,0,1},{1,0,0,0},{0,0,1,0}}
Solution 14 : {{0,1,0,0},{0,0,1,0},{0,0,0,1},{1,0,0,0}}
Solution 15 : {{0,1,0,0},{0,0,1,0},{1,0,0,0},{0,0,0,1}}
Solution 16 : {{0,1,0,0},{1,0,0,0},{0,0,0,1},{0,0,1,0}}
Solution 17 : {{0,1,0,0},{1,0,0,0},{0,0,1,0},{0,0,0,1}}
Solution 18 : {{1,0,0,0},{0,0,0,1},{0,0,1,0},{0,1,0,0}}
Solution 19 : {{1,0,0,0},{0,0,0,1},{0,1,0,0},{0,0,1,0}}
Solution 20 : {{1,0,0,0},{0,0,1,0},{0,0,0,1},{0,1,0,0}}
Solution 21 : {{1,0,0,0},{0,0,1,0},{0,1,0,0},{0,0,0,1}}
Solution 22 : {{1,0,0,0},{0,1,0,0},{0,0,0,1},{0,0,1,0}}
Solution 23 : {{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}
注記 二値変数の行列は、
qbpp::Vectorクラスを用いて実装されます。 例えば、qbpp::var("x",4,4)はqbpp::Vector<qbpp::Vector<qbpp::Var>>オブジェクトを返します。 各qbpp::Varオブジェクトはx[i][j]で表され、x[i][j]の値はsolの値は、以下のいずれかの方法で取得できます。sol(x[i][j])またはx[i][j](sol).
ベクトル関数と演算を用いた順列行列のQUBO定式化
以下の2つの関数は、行列 x 以下の2つの関数がサポートされています:
qbpp::vector_sum(x): 各行の和を計算し、それらの和を含むサイズ $n$ のベクトルを返す。qbpp::transpose(x): $n\times n$ サイズの転置行列を返す。 これらの関数を用いて、行ごとの和と列ごとの和を含むベクトルは次の式で計算できる。qbpp::vector_sum(qbpp::transpose(x)): 列ごとの和のベクトル（サイズ $n$）を返す。
スカラー-ベクトル演算を用いて各要素から1を差し引きます:
qbpp::vector_sum(x) - 1: 各行ごとの和から1を引きます。qbpp::vector_sum(qbpp::transpose(x)) - 1: 各列ごとの和から1を差し引きます。 これらの2つの$n$次元のベクトルに対して、qbpp::sqr()各要素を二乗し、qbpp::sum()全要素の和を計算します。
以下のQUBO++プログラムは、これらのベクトル関数と演算を用いたQUBO形式を実装する:
#include "qbpp.hpp"
#include "qbpp_exhaustive_solver.hpp"
int main() {
auto x = qbpp::var("x", 4, 4);
auto f = qbpp::sum(qbpp::sqr(qbpp::vector_sum(x) - 1)) +
qbpp::sum(qbpp::sqr(qbpp::vector_sum(qbpp::transpose(x)) - 1));
f.simplify_as_binary();
auto solver = qbpp::exhaustive_solver::ExhaustiveSolver(f);
auto sols = solver.search_optimal_solutions();
for (size_t k = 0; k < sols.size(); k++) {
const auto& sol = sols[k];
const auto& val_x = qbpp::onehot_to_int(x(sol));
std::cout << "Solution " << k << ": " << val_x << std::endl;
}
}
このプログラムでは、x(sol) は割り当て値の行列を x を返す。 solに割り当てられた値の行列を返します。
これはサイズ $n\times n$ の整数行列です。
各ワンホットベクトqbpp::onehot_to_int()ルを対応する整数に変換するため、
qbpp::onehot_to_int(x(sol))は置換を整数ベクトルとして返します。
このプログラムは、以下の通り全ての順列を整数ベクトルとして出力する:
Solution 0: {3,2,1,0}
Solution 1: {3,2,0,1}
Solution 2: {3,1,2,0}
Solution 3: {3,1,0,2}
Solution 4: {3,0,2,1}
Solution 5: {3,0,1,2}
Solution 6: {2,3,1,0}
Solution 7: {2,3,0,1}
Solution 8: {2,1,3,0}
Solution 9: {2,1,0,3}
Solution 10: {2,0,3,1}
Solution 11: {2,0,1,3}
Solution 12: {1,3,2,0}
Solution 13: {1,3,0,2}
Solution 14: {1,2,3,0}
Solution 15: {1,2,0,3}
Solution 16: {1,0,3,2}
Solution 17: {1,0,2,3}
Solution 18: {0,3,2,1}
Solution 19: {0,3,1,2}
Solution 20: {0,2,3,1}
Solution 21: {0,2,1,3}
Solution 22: {0,1,3,2}
Solution 23: {0,1,2,3}
割当問題とその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}\]明らかに、$X$が置換$p$を表すときのみ、$g(p) = g(X)$が成立する。
置換行列 $f(X)$ と総コスト $g(X)$ の QUBO 形式を組み合わせることで、割り当て問題の QUBO 形式を得る：
\[\begin{aligned} h(X) &= P\cdot f(x)+g(x) \\ &=P\left(\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\right)+\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}x_{i,j} \end{aligned}\]ここで、$P$ は $f(X)$ に符号化された置換制約を優先させる十分に大きな正の定数である。
割当問題のためのQUBO++プログラム
割り当て問題のためのQUBO++プログラムを設計する準備が整った。
このプログラムでは、固定行列$C$（サイズ$4\times4$）がネストされた qbpp::Vectorとして与えられる。
$f(X)$ および $g(X)$ の式は、ベクトル関数と演算を用いて定義される。
ここで、 qbpp::vector_sum(x) == 1 は等式が成立する場合に最小値0を取るQUBO式を返す。
実際、これは qbpp::sqr(qbpp::vector_sum(x) - 1)と同じQUBO式を返す。
また、 c * x は、 c と xの要素ごとの積を計算して得られる行列を返す。
したがって qbpp::sum(c * x) は g(X).
#include "qbpp.hpp"
#include "qbpp_easy_solver.hpp"
int main() {
qbpp::Vector<qbpp::Vector<uint32_t>> c = {
{58, 73, 91, 44}, {62, 15, 87, 39}, {78, 56, 23, 94}, {11, 85, 68, 72}};
auto x = qbpp::var("x", 4, 4);
auto f = qbpp::sum(qbpp::vector_sum(x) == 1) +
qbpp::sum(qbpp::vector_sum(qbpp::transpose(x)) == 1);
auto g = qbpp::sum(c * x);
auto h = 1000 * f + g;
h.simplify_as_binary();
auto solver = qbpp::easy_solver::EasySolver(h);
solver.time_limit(1.0);
auto sol = solver.search();
std::cout << "sol = " << sol << std::endl;
auto result = qbpp::onehot_to_int(x(sol));
std::cout << "Result : " << result << std::endl;
for (size_t i = 0; i < result.size(); ++i) {
std::cout << "c[" << i << "][" << result[i] << "] = " << c[i][result[i]]
<< std::endl;
}
}
Easy Solverを使用して解を求める hの解を求める。
Easy Solverオブジェクト solver に対して hの検索時間制限は、 time_limit() メンバー関数を呼び出すことで設定されます。
結果の順列は resultに格納され、選択された c[i][j] 値が順番に出力される。
このプログラムの出力は以下の通りである:
sol = 93:{{x[0][0],0},{x[0][1],0},{x[0][2],0},{x[0][3],1},{x[1][0],0},{x[1][1],1},{x[1][2],0},{x[1][3],0},{x[2][0],0},{x[2][1],0},{x[2][2],1},{x[2][3],0},{x[3][0],1},{x[3][1],0},{x[3][2],0},{x[3][3],0}}
Result : {3,1,2,0}
c[0][3] = 44
c[1][1] = 15
c[2][2] = 23
c[3][0] = 11
注記 式
fと整数m,f == mは式qbpp::sqr(f - m)を返します。 この値は、等式f == mが成立する場合に限り最小値 0 を取る。