QUBO++入門：HUBOによる素因数分解

QUBO++のHUBO対応

QUBO++ が HUBO（高次の非制約二値最適化） に全面対応しました．
生成できる式の次数に上限はありません．なお，3 種類のソルバー（Easy Solver / Exhaustive Solver / ABS3 GPU Solver）は，次数削減を行うことなく 8 次までの HUBO を直接扱えます．

HUBO(High-order Unconstrained Binary Optimization)

QUBO（Quadratic Unconstrained Binary Optimization）が二値変数の二次式を対象とするのに対し，三次以上の積項を含むように拡張した問題を HUBO と呼びます．たとえば，二値変数$a$,$b$,$c$に対して，次の多項式$f$は3次項$6abc$を含んでいるため，QUBOではなくHUBOです．

$$\begin{align*} f(a,b,c) &=(a+b+c)(a+b+c-2)^2 \\ &= a+b+c -2ab-2ac-2bc+6abc \end{align*}$$

このような式$f$を最小する変数への二値割り当てを求めるのが HUBO 問題です．上の式では，$a+b+c=0$または$a+b+c=2$を満たすときに最小値$0$をとり，それ以外では1以上の値になります．したがって最適解は，$(a,b,c)=(0,0,0),(1,1,0),(1,0,1),(0,1,1)$の４通りです．

HUBOを解くQUBO++プログラム

先ほどの $f(a,b,c)$ を記号処理で生成して最適解を求める QUBO++ の最小実装例は次のとおりです．

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

int main() {
  auto a = qbpp::var("a");
  auto b = qbpp::var("b");
  auto c = qbpp::var("c");

  auto f = (a + b + c) * qbpp::sqr(a + b + c - 2);
  f.simplify_as_binary();

  std::cout << "f = " << f << std::endl;

  auto solver = qbpp::exhaustive_solver::ExhaustiveSolver(f);
  auto sol    = solver.search();
  std::cout << sol << std::endl;
}

3つの二値変数a，b，cを定義し，qbpp::Exprオブジェクトfに$f(a,b,c)$を構築しています．メンバ関数simplify_as_binary()は式fを展開・整理を行います．その後，fを引数としてExhaustive Solverのオブジェクトを生成し，メンバ関数search()を用いて最初に見つかった最適解を求めて，solに保存しています．

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

f = a +b +c -2*a*b -2*a*c -2*b*c +6*a*b*c
0:{{a,0},{b,0},{c,0}}

展開・整理した式が表示され，最初に見つかった最適解を表示しています．

すべての最適解を求める

すべての最適解を取得するには，次のようにsearch_optimal_solutions()を用います．

  auto sol    = solver.search_optimal_solutions();
  std::cout << sol << std::endl;

これを実行すると次のように全ての最適解が表示されます．

f = a +b +c -2*a*b -2*a*c -2*b*c +6*a*b*c
0:0:{{a,0},{b,0},{c,0}}
1:0:{{a,0},{b,1},{c,1}}
2:0:{{a,1},{b,0},{c,1}}
3:0:{{a,1},{b,1},{c,0}}

すべての割当を列挙

最適でないものも含めてすべての解を求めるには，次のようにsearch_all_solutions()を用います．

  auto sol    = solver.search_all_solutions();
  std::cout << sol << std::endl;

この場合，次のように全8通りの出力が得られます．

f = a +b +c -2*a*b -2*a*c -2*b*c +6*a*b*c
0:0:{{a,0},{b,0},{c,0}}
1:0:{{a,0},{b,1},{c,1}}
2:0:{{a,1},{b,0},{c,1}}
3:0:{{a,1},{b,1},{c,0}}
4:1:{{a,0},{b,0},{c,1}}
5:1:{{a,0},{b,1},{c,0}}
6:1:{{a,1},{b,0},{c,0}}
7:3:{{a,1},{b,1},{c,1}}

HUBOによる素因数分解

HUBO 問題の入門例として素因数分解を扱います．実務上は，HUBO への変換で式が巨大化するため計算効率は良くありませんが，定式化の流れが分かりやすい題材です．ごく簡単な例として，$15$の$p=3$，$q=5$への素因数分解をHUBO問題として定式化し，解を求めます．すなわち，$15$が与えられて，$3$と$5$を求めます．整数$p$と$q$を，それぞれ3個のの二値変数$p_0,p_1,p_2$と$q_0,q_1,q_2$を用いて次の式で表します．

$$\begin{align*} p &= 1+p_0+2p_1+4p_2\\ q &= 1+q_0+2q_1+4q_2 \end{align*}$$

これにより，$p$と$q$は1以上8以下の値を保持する整数変数として扱うことができます．次のHUBO式$f(p,q)$は6個の二値変数の4次式になり，最適解0は$pq=15$を満たすときに達成されます．

$$\begin{align*} f(p,q) &= (pq-15)^2 \\ &= ((1+p_0+2p_1+4p_2)(1+q_0+2q_1+4q_2)-15)^2 \end{align*}$$

素因数分解を行うQUBO++プログラム

下のQUBO++プログラムは，$f(p,q)$のHUBO式を生成し，Exhaustive Solverで最適解を求めます．次のQUBO++では，まず，pとqを1以上8以下の値をとる整数変数として定義してます．念のため，pとqが二値変数を用いてどのように表現されているかを表示します．そして，等式p * q == 15を満たす場合に最適解0となるHUBO式をfとします．この等式は自動的に，qbpp::sqr(p * q - 15)に変換され，fに代入されます．fを展開・整理し，表示しています．

Exhaustive Solverにより，fを全探索し，得られたすべての最適解をsolsに代入し，各解を表示しています．

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

int main() {
  auto p = 1 <= qbpp::var_int("p") <= 8;
  auto q = 1 <= qbpp::var_int("q") <= 8;
  std::cout << "p = " << p << ", q = " << q << std::endl;

  auto f = p * q == 15;
  f.simplify_as_binary();
  std::cout << "f = " << f << std::endl;

  auto solver = qbpp::exhaustive_solver::ExhaustiveSolver(f);
  auto sols = solver.search_optimal_solutions();
  for (const auto& sol : sols) {
    std::cout << "p = " << sol(p) << ", q = " << sol(q) << ", sol = " << sol
              << std::endl;
  }
}

このプログラムを実行すると次の出力が得られます．この出力から，fが４次式であることがわかります．また，$(p,q)=(3,5),(5,3)$の2通りの解が得られています．

p = 1 +p[0] +2*p[1] +4*p[2], q = 1 +q[0] +2*q[1] +4*q[2]
f = 196 -27*p[0] -52*p[1] -96*p[2] -27*q[0] -52*q[1] -96*q[2] +4*p[0]*p[1] +8*p[0]*p[2] -21*p[0]*q[0] -36*p[0]*q[1] -48*p[0]*q[2] +16*p[1]*p[2] -36*p[1]*q[0] -56*p[1]*q[1] -48*p[1]*q[2] -48*p[2]*q[0] -48*p[2]*q[1] +96*p[2]*q[2] +4*q[0]*q[1] +8*q[0]*q[2] +16*q[1]*q[2] +12*p[0]*p[1]*q[0] +32*p[0]*p[1]*q[1] +96*p[0]*p[1]*q[2] +24*p[0]*p[2]*q[0] +64*p[0]*p[2]*q[1] +192*p[0]*p[2]*q[2] +12*p[0]*q[0]*q[1] +24*p[0]*q[0]*q[2] +48*p[0]*q[1]*q[2] +48*p[1]*p[2]*q[0] +128*p[1]*p[2]*q[1] +384*p[1]*p[2]*q[2] +32*p[1]*q[0]*q[1] +64*p[1]*q[0]*q[2] +128*p[1]*q[1]*q[2] +96*p[2]*q[0]*q[1] +192*p[2]*q[0]*q[2] +384*p[2]*q[1]*q[2] +16*p[0]*p[1]*q[0]*q[1] +32*p[0]*p[1]*q[0]*q[2] +64*p[0]*p[1]*q[1]*q[2] +32*p[0]*p[2]*q[0]*q[1] +64*p[0]*p[2]*q[0]*q[2] +128*p[0]*p[2]*q[1]*q[2] +64*p[1]*p[2]*q[0]*q[1] +128*p[1]*p[2]*q[0]*q[2] +256*p[1]*p[2]*q[1]*q[2]
p = 5, q = 3, sol = 0:{{p[0],0},{p[1],0},{p[2],1},{q[0],0},{q[1],1},{q[2],0}}
p = 3, q = 5, sol = 0:{{p[0],0},{p[1],1},{p[2],0},{q[0],0},{q[1],0},{q[2],1}}

さらなる高次のHUBO式

QUBO++ は 8 次項までの HUBO を次数低減なしで直接扱えます．技術的な合理性は薄いものの，8 次の例として「4 因子の積が所与値になるように整数変数を同時に求める」定式化を示します．

まず，先の素因数分解を 4 個の整数変数に拡張します．

$$\begin{align*} p &= 1+p_0+2p_1+4p_2\\ q &= 1+q_0+2q_1+4q_2\\ r &= 1+r_0+2r_1+4r_2\\ s &= 1+s_0+2s_1+4s_2 \end{align*}$$

これらが満たすべき制約を$pqrs=210 (=2\cdot 3\cdot 5\cdot 7)$と設定します．$(p,q,r,s)=(2,3,5,7)$や$(1,5,6,7)$などの整数の組が解となります．この解をもとめるためのHUBO式$f(p,q,r,s)$は，次の8次式となります．

$$\begin{align*} f(p,q,r,s) &= (pqrs-210)^2 \\ &= ((1+p_0+2p_1+4p_2)(1+q_0+2q_1+4q_2)(1+r_0+2r_1+4r_2)(1+s_0+2s_1+4s_2)-210)^2 \end{align*}$$

以下の QUBO++ プログラムはこの式を生成し、Exhaustive Solver で値0の最適解を列挙します．

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

int main() {
  auto p = 1 <= qbpp::var_int("p") <= 8;
  auto q = 1 <= qbpp::var_int("q") <= 8;
  auto r = 1 <= qbpp::var_int("r") <= 8;
  auto s = 1 <= qbpp::var_int("s") <= 8;
  std::cout << "p = " << p << ", q = " << q << ", r = " << r << ", s = " << s
            << std::endl;
  auto f = p * q * r * s == 2 * 3 * 5 * 7;
  f.simplify_as_binary();
  std::cout << "f = " << f << std::endl;
  auto solver = qbpp::exhaustive_solver::ExhaustiveSolver(f);
  auto sols = solver.search_optimal_solutions();
  for (const auto& sol : sols) {
    std::cout << "p = " << sol(p) << ", q = " << sol(q) << ", r = " << sol(r)
              << ", s = " << sol(s) << ", sol = " << sol << std::endl;
  }
}

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

p = 1 +p[0] +2*p[1] +4*p[2], q = 1 +q[0] +2*q[1] +4*q[2], r = 1 +r[0] +2*r[1] +4*r[2], s = 1 +s[0] +2*s[1] +4*s[2]
f = 43681 -417*p[0] -832*p[1] -1656*p[2] -417*q[0] -832*q[1] -1656*q[2] -417*r[0] -832*r[1] -1656*r[2] -417*s[0] -832*s[1] -1656*s[2] +4*p[0]*p[1] +8*p[0]*p[2] -411*p[0]*q[0] -816*p[0]*q[1] -1608*p[0]*q[2] -411*p[0]*r[0] -816*p[0]*r[1] -1608*p[0]*r[2] -411*p[0]*s[0] -816*p[0]*s[1] -1608*p[0]*s[2] +16*p[1]*p[2]  
省略
+16384*p[1]*p[2]*q[0]*q[2]*r[1]*r[2]*s[0]*s[2] +32768*p[1]*p[2]*q[0]*q[2]*r[1]*r[2]*s[1]*s[2] +4096*p[1]*p[2]*q[1]*q[2]*r[0]*r[1]*s[0]*s[1] +8192*p[1]*p[2]*q[1]*q[2]*r[0]*r[1]*s[0]*s[2] +16384*p[1]*p[2]*q[1]*q[2]*r[0]*r[1]*s[1]*s[2] +8192*p[1]*p[2]*q[1]*q[2]*r[0]*r[2]*s[0]*s[1] +16384*p[1]*p[2]*q[1]*q[2]*r[0]*r[2]*s[0]*s[2] +32768*p[1]*p[2]*q[1]*q[2]*r[0]*r[2]*s[1]*s[2] +16384*p[1]*p[2]*q[1]*q[2]*r[1]*r[2]*s[0]*s[1] +32768*p[1]*p[2]*q[1]*q[2]*r[1]*r[2]*s[0]*s[2] +65536*p[1]*p[2]*q[1]*q[2]*r[1]*r[2]*s[1]*s[2]
p = 1, q = 5, r = 7, s = 6, sol = 0:{{p[0],0},{p[1],0},{p[2],0},{q[0],0},{q[1],0},{q[2],1},{r[0],0},{r[1],1},{r[2],1},{s[0],1},{s[1],0},{s[2],1}}
p = 1, q = 5, r = 6, s = 7, sol = 0:{{p[0],0},{p[1],0},{p[2],0},{q[0],0},{q[1],0},{q[2],1},{r[0],1},{r[1],0},{r[2],1},{s[0],0},{s[1],1},{s[2],1}}
省略
p = 6, q = 7, r = 1, s = 5, sol = 0:{{p[0],1},{p[1],0},{p[2],1},{q[0],0},{q[1],1},{q[2],1},{r[0],0},{r[1],0},{r[2],0},{s[0],0},{s[1],0},{s[2],1}}
p = 6, q = 7, r = 5, s = 1, sol = 0:{{p[0],1},{p[1],0},{p[2],1},{q[0],0},{q[1],1},{q[2],1},{r[0],0},{r[1],0},{r[2],1},{s[0],0},{s[1],0},{s[2],0}}

出力より，8 次項を含む HUBO が正しく生成され，$0$の値を取る最適解の組が列挙されていることが分かります．

まとめ

QUBO++を用いてHUBO問題を解く簡単な例として，素因数分解を取り上げました．
なお，QUBO++は，http://qbpp.cs.hiroshima-u.ac.jp からダウンロード可能です．