# はじめに
こんにちは！
普段，アニーリングマシン(厳密にはイジングマシンですが)を使って数理最適化を研究しているkeisukeと申します！
今回は初投稿となるのですが，自分が研究の過程で「これは美しい…！」と感嘆したものについて投稿したいと思います！


# TL;DR
- 巡回セールスマン問題とは？
- 巡回セールスマン問題をQUBOモデルに落とし込む
- QUBO行列はどういうものなのか，実際に覗いてみる
- QUBO行列は美しい！


# 巡回セールスマン問題(TSP)とは
巡回セールスマン問題，通称TSPとは，複数の都市が与えられており，全ての都市を一度だけ訪問して出発地点に戻ってくるような巡回路の中で，総移動距離が最小となるような都市の訪問順序を求める問題です．
TSPは組合せ最適化問題としては一番有名な問題なので，ご存知の方も多いのではないでしょうか．

ここからは，TSPを定式化していきます．
$N$ 個の都市の集合を $C'=\{c_1, c_2,...,c_N\}$ とします．
任意の都市 $c_i, c_j (1 \leq i,j \leq N)$ 間には距離が定義されており，その距離を $d(c_i,c_j)$ とします．また，$i$ 番目の都市 $c_i$，$j$ 番目の都市 $c_j$ に関するバイナリ変数を $x_{i,j}$ を以下のように定義します．
$$
x_{i,j} = \begin{cases}
    1 & (c_iの次にc_jを訪問するとき) \\
    0 & (上記以外) \tag1
\end{cases}
$$
以上より，全ての都市を一度だけ訪問して出発地点に戻って くるような巡回路の総距離 $C$ は，次の式のように表せます．
$$
C = \sum_{i=1}^{N} \sum_{j=1}^N d(c_i, c_j) x_{i, j} \tag2
$$

なお，TSPには以下の2つの制約があります．
1. 同時訪問制約：一度に訪れることのできる都市はただ一つのみ
2. 訪問回数制約：任意の都市にはただ一度しか訪問しない

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


# QUBOモデルに落とし込む
## その前に
今回は，QUBOがどういったモデルなのかは割愛させていただきます．
多分，僕がここで書くよりググった方が遥かに分かりやすいから…泣

TSPをQUBOモデルにマッピングする際，以下のようなバイナリ変数を定義します．
$$
x_{i,j} = \begin{cases}
    1 & (時刻 t に都市 c_i を訪問する)\\
    0 & (上記以外) \tag3
\end{cases}
$$

## コスト関数
コスト関数は，ある地点から出発して全ての都市を一度だけ訪問し，元の出発地点に戻ってくるような巡回路の総移動距離です．
エネルギー関数は次のように書けます．
$$
\mathrm{H}_A = \sum_{t=1}^N \sum_{i=1}^N \sum_{j=1}^N d(c_i, c_j) x_{t, i} x_{t+1, j} \tag4
$$
$\mathrm{H_A}$ の最小値は $h_a$ をとります．
$h_a$ は問題に依存します．


## 制約項1：同時訪問制約
これは，時刻 $t$ に訪問する都市はただ一つのみであるという制約です．
本制約をある時刻 $t$ について定式化すると，次の式のようになります．
$$
\sum_{i=1}^N x_{t,i} = 1~(1\leq t\leq N) \tag5
$$

全ての時刻 $t$ が上の式を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_B}$ を導入すると，エネルギー関数は次のように書けます．
$$
\mathrm{H}_B = \sum_{t=1}^N ( \sum_{i=1}^N x_{t, i} - 1)^2 \tag6
$$
$\mathrm{H_B}$ の最小値は 0 をとります．


## 制約項2：訪問回数制約
これは，$i$ 番目の都市 $c_i$ はただ一度しか訪問しないという制約です．
本制約を $i$ 番目の都市 $c_i$ について定式化すると，次の式のようになります．
$$
\sum_{t=1}^N x_{t,i} = 1~(1 \leq i \leq N) \tag7
$$

全ての都市 $c_i$ が上の式を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_C}$ を導入すると，エネルギー関数は次のように書けます．
$$
\mathrm{H}_C = \sum_{i=1}^N ( \sum_{t=1}^N x_{t, i} - 1)^2 \tag8
$$
$\mathrm{H_C}$ の最小値は 0 をとります．


## 最終的なエネルギー関数
以上より，最終的なエネルギー関数は，正のハイパーパラメータα を用いて次のように表せます．
$$
\mathrm{H} = \mathrm{H}_A + \alpha (\mathrm{H}_B + \mathrm{H}_C) \tag9
$$
エネルギー関数 $\mathrm{H}$ は最小値 $h_a$ を取り，このとき基底解となります．
基底解が得られたときのスピンが最適解となります．
なお，エネルギー関数 $\mathrm{H}$ には定数項が出現しますが，本質的にはあってもなくても変わりません．


# QUBO行列の中身を見ていく
さて，最終的なエネルギー関数 $\mathrm{H}$ をQUBO行列として取得します．
QUBO行列の中身はどのようになっているのか見ていきましょう．
今回は，都市数 $N=5$，全都市間の距離 $d(c_i, c_j)=3$，ハイパーパラメータ $\alpha=1$ の場合でやってみます．

```
import numpy as np
import plotly.express as px

class TSP:
    def __init__(self):
        self.NUM_CITY = 5 # 都市数
        self.NUM_SPIN = self.NUM_CITY * self.NUM_CITY # QUBO行列に必要なスピン数
        self.distance_matrix = np.array([[0,3,3,3,3],[3,0,3,3,3],[3,3,0,3,3],[3,3,3,0,3],[3,3,3,3,0]]) # 都市間の距離行列
        self.ALPHA = 1 # ハイパーパラメータ α

        # 各スピンに番号を割り振る
        self.spin_index = np.arange(self.NUM_CITY * self.NUM_CITY).reshape(self.NUM_CITY, self.NUM_CITY)

        # コスト関数用のQUBO行列
        self.cost_weight = np.zeros((self.NUM_SPIN, self.NUM_SPIN))
        self.cost_bias = np.zeros((self.NUM_SPIN))
        self.cost_const = np.zeros((1))

        # 制約項のQUBO行列
        self.penalty_weight = np.zeros((self.NUM_SPIN, self.NUM_SPIN))
        self.penalty_bias = np.zeros((self.NUM_SPIN))
        self.penalty_const = np.zeros((1))

    # QUBO行列を作る
    def make_qubo_matrix(self):
        # ========== ここから式(4) ==========
        for t in range(self.NUM_CITY):
            for u in range(self.NUM_CITY):
                for v in range(self.NUM_CITY):
                    if t < self.NUM_CITY-1:
                        idx_i = self.spin_index[t, u]
                        idx_j = self.spin_index[t+1, v]
                    elif t == self.NUM_CITY-1:
                        idx_i = self.spin_index[t, u]
                        idx_j = self.spin_index[0, v]
                    coef  = self.distance_matrix[u, v]
                    if coef == 0:
                        continue
                    self.cost_weight[idx_i, idx_j] += coef
                    self.cost_weight[idx_j, idx_i] = self.cost_weight[idx_i, idx_j]

        # ========== ここから式(6) ==========
        # 二次の項
        for t in range(self.NUM_CITY):
            for u in range(self.NUM_CITY-1):
                for v in range(u+1, self.NUM_CITY):
                    idx_i = self.spin_index[t, u]
                    idx_j = self.spin_index[t, v]
                    coef = 2
                    self.penalty_weight[idx_i, idx_j] += self.ALPHA * coef
                    self.penalty_weight[idx_j, idx_i] = self.penalty_weight[idx_i, idx_j]

        # 一次の項
        for t in range(self.NUM_CITY):
            for u in range(self.NUM_CITY):
                idx = self.spin_index[t, u]
                coef = -1
                self.penalty_bias[idx] += self.ALPHA * coef

        # 定数項
        self.penalty_const[0] += self.ALPHA * self.NUM_CITY

        # ========== ここから式(8) ==========
        # 二次の項
        for u in range(self.NUM_CITY):
            for t in range(self.NUM_CITY-1):
                for tt in range(t+1, self.NUM_CITY):
                    idx_i = self.spin_index[t, u]
                    idx_j = self.spin_index[tt, u]
                    coef = 2
                    self.penalty_weight[idx_i, idx_j] += self.ALPHA * coef
                    self.penalty_weight[idx_j, idx_i] = self.penalty_weight[idx_i, idx_j]

        # 一次の項
        for u in range(self.NUM_CITY):
            for t in range(self.NUM_CITY):
                idx = self.spin_index[t, u]
                coef = -1
                self.penalty_bias[idx] += self.ALPHA * coef

        # 定数項
        self.penalty_const[0] += self.ALPHA * self.NUM_CITY

        # 式(9)のHのQUBO行列を作る
        qubo_matrix = np.zeros((self.NUM_SPIN, self.NUM_SPIN))
        qubo_matrix = self.cost_weight + self.penalty_weight
        for i in range(self.NUM_SPIN):
            qubo_matrix[i, i] = self.cost_bias[i] + self.penalty_bias[i]
        fig = px.imshow(qubo_matrix)
        fig.show()
```

上記のコードを実行すると，下記のような結果が得られます．
(ちなみにplotlyというpythonライブラリ，ものすごく便利なのでぜひご活用ください！)
各セル(マス？)は各スピン間の相互作用に対応しています．
今回は都市数が5ということで，QUBOモデルに必要なスピン数は $5 \times 5=25$ となります．
ですので，QUBO行列は $25 \times 25$ の行列となるわけです．
なお，下図では色と各値が対応しており，その関係は以下のようになっています．
- 紺色：-2
- 紫色：0
- 黄色：2
- オレンジ色：3

![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-0061-7526-97%2F17edf766f022ff?alt=media&token=090f003c-eabf-4de1-aa01-c10887d820e0)

QUBO行列は対称行列となっていることが分かりますね．
これは当然といえば当然です．
なぜなら，スピン $i$ から見たスピン $j$ との間の相互作用と，スピン $j$ から見たスピン $i$ との間の相互作用は同じになるからです(逆に違ったらおかしいですよね)．
それにしてもQUBO，う，美しい…！！！


# おわりに
いかがでしたでしょうか？
QUBOは対称行列となっており，それを可視化すると美しいものになるということについて投稿させていただきました．
皆さんも，TSPに限らずいろんなQUBOモデルを可視化してみて，その美しさを味わってみてください！

QUBO行列の中身を見てみよう

巡回セールスマン問題(TSP)とは

# はじめに
こんにちは！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](https://www.wolframalpha.com/input/?i=%28x_1+%2B+x_2+%2B+x_3+%2B+x_4+-+1%29%5E2&lang=ja)で計算させた結果もご覧になってください．

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

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


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

Keisuke

QUBO行列の中身を見てみよう

QUBOモデルを数式からしっかり理解しよう