QUBO＋＋入門：行列操作（順列生成）

変数のベクトル

QUBO++は変数のベクトルをqbpp::Vectorを用いて扱います．
その変数の和（sum）や積（product）を計算するための関数として次の２つが用意されています．

• qbpp::sum(): 全要素の和を返します．
• qbpp::product()：全要素の積を返します．
  以下のプログラムでは、要素数4の変数ベクトルxを作り，それらの和と積を求め表示しています．

#include "qbpp.hpp"

int main() {
  qbpp::Vector<qbpp::Var> x = qbpp::var("x", 4);
  std::cout << "sum(x) = " << qbpp::sum(x) << std::endl;
  std::cout << "product(x) = " << qbpp::product(x) << std::endl;
}

理解を深めるため，コード中では型を明示してqbpp::Vector<qbpp::Var>と記述していますが，autoを使用することで簡略化することも可能です．

• 変数の生成：
  qbpp::varを用いて，大きさ4のベクトルを表す変数xを生成します．引数として行列のサイズを指定します．
• 和の計算：qbpp::sumを使って全変数の和を求めます．
• 積の計算：qbpp::productを使って全変数の積を求めます．
  このプログラムを実行すると次の結果が得られます．

sum(x) = x[0] +x[1] +x[2] +x[3]
product(x) = x[0]*x[1]*x[2]*x[3]

変数の行列

C++では，ネストしたstd::vector，具体的にはstd::vector<std::vector<int>>を使用して行列を表現することが一般的です．QUBO++も同様に，ネストしたqbpp::Vectorを用いて変数や式の行列を扱うことができます．
QUBO++では，このような行列に対して，和（sum）や積（product）を計算するための関数が専用関数が用意されています．
ただし，一次元ベクトル用のqbpp::sum()やqbpp::product()は用いず，行列に対してはバグを防ぐために、あえて別の関数名が使われています．

• qbpp::total_sum(): 全要素の和を返します．
• qbpp::vector_sum(): 低次元（行）のベクトルをその和で置き換えたベクトルを返します．
• qbpp::total_product(): 全要素の積を返します．
• qbpp::vector_product(): 低次元（行）のベクトルをその積で置き換えたベクトルを返します．
• qbpp::transpose(): 行列を転置（行と列の入れ替え）します．

以下は，$2 \times 3$の行列$x$におけるインデックスを表しています．

$$\begin{pmatrix} x[0][0] & x[0][1] & x[0][2] \\ x[1][0] & x[1][1] & x[1][2] \end{pmatrix}$$

この行列を使って，上述の関数の動作を説明します．

#include "qbpp.hpp"

int main() {
  qbpp::Vector<qbpp::Vector<qbpp::Var>> x = qbpp::var("x", 2, 3);
  qbpp::Expr total_sum = qbpp::total_sum(x).simplify_as_binary();
  qbpp::Vector<qbpp::Expr> vector_sum =
      qbpp::vector_sum(x).simplify_as_binary();
  std::cout << "total_sum = " << total_sum << std::endl;
  std::cout << "vector_sum = " << vector_sum << std::endl;

  qbpp::Expr total_product = qbpp::total_product(x).simplify_as_binary();
  qbpp::Vector<qbpp::Expr> vector_product =
      qbpp::vector_product(x).simplify_as_binary();
  std::cout << "total_product = " << total_product << std::endl;
  std::cout << "vector_product = " << vector_product << std::endl;

  qbpp::Vector<qbpp::Vector<qbpp::Expr>> transpose =
      qbpp::transpose(x).simplify_as_binary();
  std::cout << "tranpsose = " << transpose << std::endl;
  qbpp::Vector<qbpp::Expr> vector_sum_transpose =
      qbpp::vector_sum(transpose).simplify_as_binary();
  std::cout << "vector_sum(transpose) = " << vector_sum_transpose << std::endl;
}

ここでも理解を深めるために型を明示していますが，autoを用いて簡潔に記述することも可能です．

• 変数の生成：
  qbpp::varを用いて，大きさ$2 \times 3$の行列を表す変数xを生成します．引数として行列のサイズを指定します．
• 全要素の和と行ごとの和の計算：qbpp::total_sumを使って全変数の和を，qbpp::vector_sumを使って行ごとの和を計算し，それぞれを出力します．
• 全要素の積と行ごとの積の計算：qbpp::total_productを使って全変数の積を，qbpp::vector_productを使って行ごとの積を計算し，それぞれを出力します．
• 行列の転置：
  行列xを転置して新たな行列tを作成します．さらに転置後の行列tに対して行ごとの和を計算し出力します．これにより，元の行列xの列ごとの和が求まります．

以下は，このプログラムを実行した場合の出力です．

total_sum = x[0][0] +x[0][1] +x[0][2] +x[1][0] +x[1][1] +x[1][2]
vector_sum = {[0],x[0][0] +x[0][1] +x[0][2]},{[1],x[1][0] +x[1][1] +x[1][2]}
total_product = x[0][0]*x[0][1]*x[0][2]*x[1][0]*x[1][1]*x[1][2]
vector_product = {[0],x[0][0]*x[0][1]*x[0][2]},{[1],x[1][0]*x[1][1]*x[1][2]}
tranpsose = {[0][0],x[0][0]},{[0][1],x[1][0]},{[1][0],x[0][1]},{[1][1],x[1][1]},{[2][0],x[0][2]},{[2][1],x[1][2]}
vector_sum(transpose) = {[0],x[0][0] +x[1][0]},{[1],x[0][1] +x[1][1]},{[2],x[0][2] +x[1][2]}

one-hotベクトルによる順列表現

組合せ最適化問題には，あらゆる順列の中から最適なものを選択するタイプの問題が多く存在します．そのような組合せ最適化問題をQUBO問題として扱うための準備として，QUBO++を用いて順列を生成するQUBO式を記述します．

QUBOでは，整数$0, 1, 2, \ldots, n-1$の順列を表現するために，$n \times n$の二値変数の行列を用います．ここで，各行はone-hotベクトルと仮定します．この1が存在する列番号によって，$0$から$n-1$までの整数を表現します．さらに，もし各列もone-hotベクトルであるならば，各行で表される整数はすべて異なります．したがって，その整数を並べることで順列となります．

以下の$4 \times 4$の行列は，各行および各列がone-hotベクトルとなっています．この行列が表す順列は，$[2, 0, 3, 1]$です．

$$\begin{pmatrix} 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ \end{pmatrix}$$

$n \times n$の二値変数の行列$x$において，各行および各列がone-hotベクトルとなる条件は，それぞれの行・列の要素の合計が1になることです．この条件は，次のQUBO式で表すことができます．

$$\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$$

この式は，QUBO++の等価演算子==を用いて次のように表すことができます．

$$\sum_{i=0}^{n-1}\left(\sum_{j=0}^{n-1}x_{i,j}==1\right)+\sum_{j=0}^{n-1}\left(\sum_{i=0}^{n-1}x_{i,j}==1\right)$$

順列生成QUBO＋＋プログラム

$n \times n$の変数行列xを生成し、各行および各列がone-hotベクトルのとき、エネルギー値が最小となるQUBO式fを作成します。以下は、$4 \times 4$の行列の場合のコードです。

#include "qbpp.hpp"
#include "qbpp_exhaustive_solver.hpp"

int main() {
  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);
  f.simplify_as_binary();
  std::cout << "f = " << f << std::endl;

  auto solver = qbpp::exhaustive_solver::ExhaustiveSolver(f);
  auto sol = solver.search();
  auto result = sol.get(x);
  auto onehot = qbpp::onehot_to_int(result);

  std::cout << sol << std::endl;
  std::cout << str(result, "result") << std::endl;
  std::cout << str(onehot, "onehot") << std::endl;
}

式fは1行で生成されています．

• qbpp::sum(x)は，行列xの各行をその和で置き換えたベクトルを返します．
• qbpp::sum(x) == 1は，そのベクトルの各要素が1である場合に最小値0を取るQUBO式を返します．
• qbpp::sum(qbpp::vector_sum(x) == 1)は，そのQUBO式の和を返します．
• qbpp::sum(qbpp::vector_sum(qbpp::transpose(x)) == 1)は，行列xの転置行列に対して同じ処理を行います．
  以上により，式fは各行および各列がone-hotベクトルとなる場合にエネルギーが最小値0を取るQUBO式となります．

Exhaustive Solverを用いて最適解solを求めます．そして，メンバ関数get()を使用して，solからxの値を取得し，それをresultとします．xが行列であるため，resultも同様に行列形式になります．さらに，関数qbpp::onehot_to_int()を用いて，resultの各行をone-hotベクトルから整数に変換し，onehotとして格納します．最後に，std::coutを用いてsol，result，およびonehotを出力します．

このプログラムを実行すると，下の出力が得られます．

f = 8 -2*x[0][0] -2*x[0][1] -2*x[0][2] -2*x[0][3] -2*x[1][0] -2*x[1][1] -2*x[1][2] -2*x[1][3] -2*x[2][0] -2*x[2][1] -2*x[2][2] -2*x[2][3] -2*x[3][0] -2*x[3][1] -2*x[3][2] -2*x[3][3] +2*x[0][0]*x[0][1] +2*x[0][0]*x[0][2] +2*x[0][0]*x[0][3] +2*x[0][0]*x[1][0] +2*x[0][0]*x[2][0] +2*x[0][0]*x[3][0] +2*x[0][1]*x[0][2] +2*x[0][1]*x[0][3] +2*x[0][1]*x[1][1] +2*x[0][1]*x[2][1] +2*x[0][1]*x[3][1] +2*x[0][2]*x[0][3] +2*x[0][2]*x[1][2] +2*x[0][2]*x[2][2] +2*x[0][2]*x[3][2] +2*x[0][3]*x[1][3] +2*x[0][3]*x[2][3] +2*x[0][3]*x[3][3] +2*x[1][0]*x[1][1] +2*x[1][0]*x[1][2] +2*x[1][0]*x[1][3] +2*x[1][0]*x[2][0] +2*x[1][0]*x[3][0] +2*x[1][1]*x[1][2] +2*x[1][1]*x[1][3] +2*x[1][1]*x[2][1] +2*x[1][1]*x[3][1] +2*x[1][2]*x[1][3] +2*x[1][2]*x[2][2] +2*x[1][2]*x[3][2] +2*x[1][3]*x[2][3] +2*x[1][3]*x[3][3] +2*x[2][0]*x[2][1] +2*x[2][0]*x[2][2] +2*x[2][0]*x[2][3] +2*x[2][0]*x[3][0] +2*x[2][1]*x[2][2] +2*x[2][1]*x[2][3] +2*x[2][1]*x[3][1] +2*x[2][2]*x[2][3] +2*x[2][2]*x[3][2] +2*x[2][3]*x[3][3] +2*x[3][0]*x[3][1] +2*x[3][0]*x[3][2] +2*x[3][0]*x[3][3] +2*x[3][1]*x[3][2] +2*x[3][1]*x[3][3] +2*x[3][2]*x[3][3]
0:{{x[0][0],0},{x[0][1],0},{x[0][2],0},{x[0][3],1},{x[1][0],1},{x[1][1],0},{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],0},{x[3][1],1},{x[3][2],0},{x[3][3],0}}
{result[0][0],0},{result[0][1],0},{result[0][2],0},{result[0][3],1},{result[1][0],1},{result[1][1],0},{result[1][2],0},{result[1][3],0},{result[2][0],0},{result[2][1],0},{result[2][2],1},{result[2][3],0},{result[3][0],0},{result[3][1],1},{result[3][2],0},{result[3][3],0}
{onehot[0],3},{onehot[1],0},{onehot[2],2},{onehot[3],1}

この出力では，順列[2, 1, 3, 0]が得られています．

全ての順列を求める

Exhaustive Solverは，メンバ関数search_optimal_solutions()を使用することで，全ての最適解を列挙することができます．以下は，この機能を利用して全ての順列を表示するQUBO++プログラムです．

#include "qbpp.hpp"
#include "qbpp_exhaustive_solver.hpp"

int main() {
  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);
  f.simplify_as_binary();

  auto solver = qbpp::exhaustive_solver::ExhaustiveSolver(f);
  auto sols = solver.search_optimal_solutions();
  for (const auto& sol : sols) {
    std::cout << qbpp::onehot_to_int(sol.get(x)) << std::endl;
  }
}

search_optimal_solutions()を使用して，全ての最適解をsolsに保存します．それをforループで取り出し，表示しています．

プログラムを実行すると，次のような出力が得られ，全24通りの順列が表示されます．

{[0],3},{[1],2},{[2],1},{[3],0}
{[0],3},{[1],2},{[2],0},{[3],1}
{[0],3},{[1],1},{[2],2},{[3],0}
{[0],3},{[1],1},{[2],0},{[3],2}
{[0],3},{[1],0},{[2],2},{[3],1}
{[0],3},{[1],0},{[2],1},{[3],2}
{[0],2},{[1],3},{[2],1},{[3],0}
{[0],2},{[1],3},{[2],0},{[3],1}
{[0],2},{[1],1},{[2],3},{[3],0}
{[0],2},{[1],1},{[2],0},{[3],3}
{[0],2},{[1],0},{[2],3},{[3],1}
{[0],2},{[1],0},{[2],1},{[3],3}
{[0],1},{[1],3},{[2],2},{[3],0}
{[0],1},{[1],3},{[2],0},{[3],2}
{[0],1},{[1],2},{[2],3},{[3],0}
{[0],1},{[1],2},{[2],0},{[3],3}
{[0],1},{[1],0},{[2],3},{[3],2}
{[0],1},{[1],0},{[2],2},{[3],3}
{[0],0},{[1],3},{[2],2},{[3],1}
{[0],0},{[1],3},{[2],1},{[3],2}
{[0],0},{[1],2},{[2],3},{[3],1}
{[0],0},{[1],2},{[2],1},{[3],3}
{[0],0},{[1],1},{[2],3},{[3],2}
{[0],0},{[1],1},{[2],2},{[3],3}

まとめ

• 順列の表現方法：$n \times n$ の二値行列を用いることで$0, 1, \ldots, n-1$ の順列を表現できます．
• 順列を求めるQUBO式：$n \times n$ の二値変数行列において，各行および各列が one-hot ベクトルとなるように制約を設定することで，最適解が順列を表すQUBO 式を構築できます．

QUBO++情報

• QUBO＋＋ツールとABS2 QUBO Solverの詳細情報：QUBO++に関する公式サイトです．ツールの特徴や使用方法，ABS2 QUBO Solverとの連携について詳しく説明されています．
• 【研究成果】組合せ最適化問題を解くためのプログラミング用ツールQUBO＋＋を無償公開します：広島大学からの公式リリース情報です．QUBO++の開発背景や研究成果について記載されています．
• 利用条件：QUBO++は非商用利用や評価研究目的に限り無償で利用可能です．商用利用や定常的な業務での使用を希望する場合は，ライセンス契約が必要となります．