はじめに

現在、量子アニーリングの学習を進めています。
最適化問題を扱うにあたっては数学の知識が不可欠ですが、数式を見るとどうしても身構えてしまいます…。
そこで「数式 → プログラム」に分解しながら、実装ベースで定式化を理解していくことにしました。
本記事では、巡回セールスマン問題（TSP）を題材に、QUBOの構造をシンプルに捉えることを目指して書いています。
※初学者の視点で書いているため、間違いやご指摘などありましたらコメントで教えていただけると嬉しいです！

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

複数の都市や観光地などいろんな場所を巡ろうとした時に、距離や時間などが最短になる組み合わせを見つける問題です。
以下の説明は ChatGPT からの引用です。

巡回セールスマン問題（TSP: Traveling Salesman Problem）は、複数の都市をすべて一度ずつ訪れ、元の都市に戻るときの「総移動距離（またはコスト）」を最小にする巡回ルートを求める問題です。

都市の数が増えると、可能な巡回ルートの数は急激に増加します（都市数が N ならパターンは (N-1)! 通り）。そのため、厳密解を求めるのが非常に難しくなる組合せ最適化問題のひとつです。

本記事では、完全グラフ（すべての都市間に移動可能） を前提とし、TSPを QUBO（Quadratic Unconstrained Binary Optimization）形式に落とし込む例を扱います。

TSPの最小定式化（完全グラフ・2制約のみ）

都市数：$N$
バイナリ変数：$x_{a,i}$ は「時刻 $a$ に都市 $i$ を配置するとき $1$、それ以外は $0$」

距離行列：$d_{ij} \ge 0$（都市 $i$ → $j$ の距離）

$\text{minimize}$（目的関数）:

$$
\sum_{a=0}^{N-1}\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}
d_{ij} \cdot x_{a,i} \cdot x_{(a+1) \bmod N,, j}
$$

$\text{subject to}$（制約）:

各時刻に1都市のみ配置：

$$
\sum_{i=0}^{N-1} x_{a,i} = 1 \quad (a = 0, \dots, N-1)
$$

各都市は1回だけ訪問：

$$
\sum_{a=0}^{N-1} x_{a,i} = 1 \quad (i = 0, \dots, N-1)
$$

バイナリ制約：

$$
x_{a,i} \in {0, 1}
$$

ざっくりと各項を咀嚼してみる

巡回セールスマン問題に限らず、多くの最適化問題では基本的に次の2つの構造を持っています。

目的項（目的関数）：解として「良い」状態（＝距離が短い、コストが小さい）を評価する項目
制約項：解として「ダメな」状態を罰則で排除するための項目

この 2 つを足し合わせて、最終的なスコア（＝エネルギー）を計算し、「一番良い状態＝最小スコアの組み合わせ」を探すのが QUBO モデルの基本です。

この記事では、これらの構造を具体的に理解するために、都市数を 4 つと仮定して、
目的項 → 制約項の順に、それぞれどういう意味なのかを簡単な例と一緒に見ていこうと思います。

図解してみよう

まず “状態” を行列で可視化し、バイナリ変数 $x_{a,i}$ が何を意味するかを掴みます。
経路の順番は時刻と表現していきます。

今回は4の都市で、それぞれの都市には1回だけ(重複なし)訪問する必要があります。
そのため、前提として4 回の移動（＝4 ステップ）が発生します。図の例だとDを起点とした時、Dを起点とするとD→A→B→Cの順で周り、最後にC→Dの都市に戻ります。
これを図示すると以下の表のようになります。

右の表では

0 ：『その都市を訪れていない』
1 ：『その都市を訪れている』

のようにt/fの形で表すことができています。行(時刻)で見ると1つの配列のように見ることができ、これがバイナリ変数というものになります。
例えば都市Aにいる時は、[1, 0, 0, 0]のようになりますが、[1, 1, 0, 0]や[1, 1, 1, 1]のように複数の都市にいる状態や、[0, 0, 0, 0]のようにどの都市にも存在しないような状態が表すことができます。
つまり一つの時刻に固定した時は2^4(0~1の組み合わせが2^都市数)=16通り表すことができます。
このような内容を一旦前提として以下の目的項について説明ができればと思います。

目的項（移動コスト）

$\sum$で囲まれている式であるため、プログラミング的な視点で言えばforループになっていると思います。
1つのステップ(時刻1→2への移動)を書き出すと以下のような形になるかと思います。

式の説明にあるように$d_{i,j}$は都市同士の距離を格納する2重配列のようなイメージです。またバイナリ変数$x$は、0か1の値を持つ配列として考えることができます。それぞれ$i$と$j$を0から3まで数値を加算していきそれぞれの配列から0か1をひっぱってきているだけになります。
ここまで書き下すと以下のようなイメージとなります。

こうなると単純に2都市間の距離と、バイナリ変数から引っ張ってきた0か1を積算した値をi=0→3、j=0→3まで全て足し合わせただけとなります！

プログラム的なイメージだと以下の通りだと思います。

# dは距離の二重配列
for i in range(4): # 今の時刻
    for j in range(4): # 次の時刻
        # x_a[i]: 今の時刻の0/1
        # x_next[j]: 次の時刻0/1
        cost += d[i][j] * x_a[i] * x_next[j]

最終的に全てのステップ（時刻1->2、時刻2->3、時刻3→4、時刻4->0）の算出結果を加算したものをコストとして出します。

先ほどバイナリ変数の説明で示した以下のような

[0,0,0,0]　　　　　・・・どの都市にもいない
[1,0,1,0], [1, 1, 1, 1]　・・・複数の都市にいる

のようなパターンもあり、こちらも計算されています。
例えば全てのバイナリ変数が0の時は、そのパターンが一番コストが低い組み合わせとなってしまいます。
このようなパターンを以下の制約項でペナルティーを科す処理を追加していきます。

制約項①：各時刻に都市は1つ

上述した$x$のバイナリ変数の配列において、本来取り得ない1が複数入った状態を取り除くために制約項が設定されております。
制約項①は各時刻にて、バイナリの配列に1が複数個含まれていないかを確認するものとなります。
式を見てみると、バイナリ変数の加算を表しており、例えば[1, 0, 0, 0]の場合は、
$1 + 0 + 0 + 0=1$となり、式自体が0となるのでペナルティーが0となります。
仮に[1, 1, 1, 0]のような本来取り得ない状態となれば、
$1 + 1 + 1 + 0 = 3$となり、$(3 - 1)^2=4$となり大きなペナルティがつくことになります。
このような算出を全ての時刻において行った加算結果が第一項目のペナルティとなります。

制約項②：各都市は1回だけ

第2項は、同じ都市を再訪問していないかを確認する制約項となります。
バイナリ変数を[1, 0, 0, 0]のように表現していましたが、
複数の時刻も含めた場合2重配列として捉えることができるかと思います。
例えば都市をD→A→B→Cの順で巡ると仮定した場合で取りうる各バイナリ変数は、

時刻1 = [0, 0, 0, 1]
時刻2 = [1, 0, 0, 0]
時刻3 = [0, 1, 0, 0]
時刻4 = [0, 0, 1, 0]

のような形となり、時系列として表すと

[
    [0, 0, 0, 1], # 時刻1
    [1, 0, 0, 0], # 時刻2
    [0, 1, 0, 0], # 時刻3
    [0, 0, 1, 0]  # 時刻4
]

のような配列となります。
配列を縦に見る、pythonでいうスライスを使うと[:, i]のような形で見ると都市Aについては[0,1,0,0]が得られて、この場合では1回のみ訪れていることがわかります。
A~D都市に対してそれぞれの時刻で各都市に何回訪れたかが確認できると思います。
あとは制約項①と同様に一回のみ訪れていることをチェックするようなイメージです。

終わりに

最終的なコストは目的項と制約項を加算結果となります。
ほとんどが制約項に引っかかるバイナリ変数の組み合わせですが、全てのバイナリ変数の組み合わせと、その組み合わせから不要なものを取り除くためのペナルティーをつけると認識すると考えやすくなるのかなと思います。
※本来であれば制約項に係数をつけてペナルティの効き具合を調整しますが簡単のため省略いたしました。

プログラム的に考えると、今回の例では以下のようなイメージになるかと思います。

# xs: バイナリ変数
# d: 都市間の距離の配列
# K: ペナルティ係数 組み合わせの中の最大コストの値などを設定
def calc_cost(xs, d, K=100):
    N = len(xs)
    # 行ペナルティ・列ペナルティ
    penalty = sum((sum(row) - 1)**2 for row in xs)
    penalty += sum((sum(xs[a][i] for a in range(N)) - 1)**2 for i in range(N))

    # 距離コスト
    cost = 0
    for a in range(N):
        i = xs[a].index(1)
        j = xs[(a+1)%N].index(1)
        cost += d[i][j]

    return cost + K * penalty

それぞれのバイナリ変数の組み合わせが入った時にどうなるかを考えるとわかりやすくなるのかなと思います。
そもそも最終的にプログラムにしてソルバーなどに入れるものであるので、
最初からプログラムとして考えるのも良さそうかなと思います。
バイナリ変数を使用するような最適化問題(ナップザック問題等)も同様の考え方で理解が深まるのではと考えております。