はじめに

こんにちは！keisukeです．
今回は，QUBOモデル(特に，二次割当問題)におけるエネルギー関数を詳しく見ていきます．
「アニーリングマシンで問題を解く！」の1個前の段階のお話です．

TL;DR

• 二次割当問題とは
• 二次割当問題をQUBOモデルに落とし込む
• 二次割当問題のエネルギー関数を展開して，エネルギー関数の中身を詳しく見る．

二次割当問題(QAP)

二次割当問題，通称QAPでは，複数の工場と，工場と同数の地区が与えられます．
任意の工場間には物流量が，任意の地区間には距離が定義されています．
二次割当問題とは，「各工場間の物流量」と「工場が割り当てられた各地区間の距離」の積の総和が最小となるような，工場の地区への割り当てを求める問題です．

ではここからは，QAPを定式化していきましょう．
$F$ 個の工場の集合を $Fac = \{f_1, f_2, ..., f_F\}$，工場を割り当てるための $F$ 個の地区の集合 $Loc = \{l_1, l_2, ..., l_F\}$ と定義します．
任意の工場間には物流量が与えられており，工場 $f_i, f_j (1 \leq i, j \leq F )$ 間の物流量を $w(f_i, f_j)$ とします．
また，任意の地区間にも距離が与えられており，工場 $f_i$ が割り当てられた地区と工場 $f_j$ が割り当てられた地区間の距離を $d(f_i, f_j)$ とします．
したがって，「各工場間の物流量」と「工場が割り当てられた各地区間の距離」の積の総和 $C$ は，次の式(1)のように表せます．

$$C = \sum_{i=1}^{F-1} \sum_{j=1, i+1}^F w(f_i, f_j) d(f_i, f_j) \tag1$$

なお，QAPには以下の2つの制約があります．

1. 工場配置制約：任意の工場はただ一つの地区にのみ必ず配置する
2. 地区内工場制約：任意の地区にはただ一つの工場が割り当てられる

以上がQAPの問題定義および定式化(とその制約)です．

QUBOモデルに落とし込む

その前に

QAPをQUBOモデルに落とし込む際，以下のようなバイナリ変数を定義します．

$$x_{i,k} = \begin{cases} 1 & (工場 f_i を地区 l_k に配置するとき)\\ 0 & (上記以外) \tag2 \end{cases}$$

コスト関数

コスト関数は，「各工場間の物流量」と「工場が配置された各地区間の距離」の積の総和です．
エネルギー関数は次のように書けます．

$$\mathrm{H}_A = \sum_{i=1}^F \sum_{j=1}^F \sum_{k=1}^F \sum_{l=1}^F w(f_i, f_j) d(l_k, l_l) x_{i, k} x_{j, l} \tag3$$

$\mathrm{H_A}$ の最小値は $h_a$ をとります．
$h_a$ は問題に依存します．

制約項1：工場配置制約

これは，$i$ 番目の工場 $f_i$ はただ一つの地区にのみ配置するという制約です．
本制約をある時刻 $i$ 番目の工場 $f_i$ について定式化すると，次の式のようになります．

$$\sum_{k=1}^F x_{i,k} = 1~(1\leq i\leq F) \tag4$$

全ての工場 $f_i$ が上の式(4)を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_B}$ を導入すると，エネルギー関数は次のように書けます．

$$\mathrm{H}_B = \sum_{i=1}^F \Biggl( \sum_{k=1}^F x_{i, k} - 1 \Biggr)^2 \tag5$$

$\mathrm{H_B}$ の最小値は 0 をとります．

制約項2：地区内工場制約

これは，$k$ 番目の地区 $l_k$ にはただ一つの工場が配置されているという制約です．
本制約を $k$ 番目の地区 $l_k$ について定式化すると，次の式のようになります．

$$\sum_{i=1}^F x_{i,k} = 1~(1 \leq k \leq F) \tag6$$

全ての地区 $l_k$ が上の式(6)を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_C}$ を導入すると，エネルギー関数は次のように書けます．

$$\mathrm{H}_C = \sum_{k=1}^F \Biggl( \sum_{i=1}^F x_{i, k} - 1 \Biggr)^2 \tag7$$

$\mathrm{H_C}$ の最小値は 0 をとります．

最終的なエネルギー関数

以上より，最終的なエネルギー関数は，正のハイパーパラメータα を用いて次のように表せます．

$$\mathrm{H} = \mathrm{H}_A + \alpha (\mathrm{H}_B + \mathrm{H}_C) \tag8$$

エネルギー関数 $\mathrm{H}$ は最小値 $h_a$ を取り，このとき基底解となります．
基底解が得られたときのスピンが最適解となります．
なお，エネルギー関数 $\mathrm{H}$ には定数項が出現しますが，本質的にはあってもなくても変わらないので，今回は定数項は無視します．

エネルギー関数を展開する

再掲ですが，QAPのエネルギー関数は以下の3つから構成されています．
$\begin{aligned} \mathrm{H}_A &= \sum_{i=1}^F \sum_{j=1}^F \sum_{k=1}^F \sum_{l=1}^F w(f_i, f_j) d(l_k, l_l) x_{i, k} x_{j, l}\\ \mathrm{H}_B &= \sum_{i=1}^F \Biggl( \sum_{k=1}^F x_{i, k} - 1 \Biggr)^2\\ \mathrm{H}_C& = \sum_{k=1}^F \Biggl( \sum_{i=1}^F x_{i, k} - 1 \Biggr)^2 \end{aligned}$

これらのエネルギー関数を何かしらのプログラミング言語で実装しようとしたとき，どうやって書けば良いか困ってしまうことはありませんか？
$\mathrm{H}_A$ はforループを4つ使って愚直に書けば良さそうですが，$\mathrm{H}_B, \mathrm{H}_C$ はシグマの中にシグマの2乗が含まれています．これは難しそうだ…

でも大丈夫！！！
ここから下を読めば，一発で分かります．
ということで早速，$\mathrm{H}_B, \mathrm{H}_C$ を展開していきましょう．

$\mathrm{H}_B$ の展開

$\mathrm{H}_B$ を展開します．
$\begin{aligned} \mathrm{H}_B &= \sum_{i=1}^F \Biggl( \sum_{k=1}^F x_{i, k} - 1 \Biggr)^2 \\ &= \sum_{i=1}^F \Biggl(2\sum_{k=1}^{F-1} \sum_{l=k+1}^{F} x_{i, k} x_{i, l} - \sum_{k=1}^{F} x_{i, k} + 1 \Biggr) \\ &= 2 \sum_{i=1}^F \sum_{k=1}^{F-1} \sum_{l=k+1}^{F} x_{i, k} x_{i, l} - \sum_{i=1}^{F} \sum_{k=1}^{F}x_{i, k} + F \end{aligned}$

なぜ1行目から2行目になるのか？
これは簡単な例で試してみると，そうなることが分かります．
ここでは，工場数 $F=4$ の場合で考えてみましょう．
(以下，2乗の部分だけを計算するので，$\mathrm{H}_B$ の一番外側のシグマ $\sum_{i=1}^{F}$ は除外します)
$\begin{aligned} \Biggl( \sum_{k=1}^{F=4} x_k - 1 \Biggr)^2 &= \biggl( (x_1 + x_2 + x_3 + x_4) - 1 \biggr)^2 \\ &= x_1^2 + x_2^2 + x_3^2 + x_4^2 \\ &~~~~~~~~~+2x_1x_2 + 2x_1x_3 + 2x_1x_4 \\ &~~~~~~~~~+2x_2x_3 + x_2x_4 + 2x_3x_4 \\ &~~~~~~~~~-2x_1 - 2x_2 - 2x_3 - 2x_4 + 1\\ &~~~~~(… 変数x_kはバイナリ変数なので x_k^2 = x_k となるから) \\ &= 2(x_1x_2 + x_1x_3 + x_1x_4 + x_2x_3 + x_2x_4 + x_3x_4) \\ &~~~~~~~-x_1 - x_2 - x_3 - x_4 + 1 \\ &= 2\sum_{k=1}^{4-1} \sum_{l=k+1}^{4} x_k x_l - \sum_{k=1}^{4} x_k + 1 \\ &= 2\sum_{k=1}^{F-1} \sum_{l=k+1}^{F} x_k x_l - \sum_{k=1}^{F} x_k + 1 \end{aligned}$
ほら！これで展開ができました．
お前ホンマかいな…という人は，神ツールである
Wolfram alphaで計算させた結果もご覧になってください．

… というわけで，エネルギー関数に2乗が含まれている場合も，今回の展開で最終的に得られた式をプログラムに落とし込めば，うまく行きそうですね！(というかうまく行きます笑)

$\mathrm{H}_C$ の展開も，扱う文字が違うだけで内容はほぼ全く同じなので，ぜひ練習として $\mathrm{H}_C$ を展開してみてください．

おわりに

今回は，QUBOモデルで頻繁に出てくる2乗式を詳しく見ました．
このように，実際に手を動かすことが大事となってくるので，皆さんもQUBOと格闘する時はどんどん手を動かしていきましょう！