定式化の高速化

大規模な組合せ最適化問題を定式化したい場合、Python の for 文などを用いて定式化すると大きく時間がかかることがあります。 Amplify SDK は大規模な最適化問題を実用的な時間で定式化するための高速な定式化手段を提供しています。 例えば、PolyArray クラスの Numpy-like なメソッドを用いることで、大幅な高速化が可能です。

このページでは、巡回セールスマン問題のバイナリ変数による定式化を例に、高速な定式化を行ういくつかの方法について説明します。

巡回セールスマン問題の定式化

都市数が \(n\) である巡回セールスマン問題は、以下のように \((n+1)n\) 個のバイナリ変数を用いて定式化されます。ただし \(d\) は \(d_{i,j}\) が都市 \(i\) と都市 \(j\) の距離を表す距離行列であるとします。

参考

ここでは、定式化の意味は重要ではないため、変数や目的関数、制約条件が何を指すのかについては解説しません。巡回セールスマン問題の定式化の詳細については 巡回セールスマン問題 を参照してください。

\[\begin{split} \begin{align} \text{minimize} \quad & \sum_{0 \leq i, j, k < n} d_{i, j} q_{k, i} q_{k+1, j} & \\ \text{subject to} \quad & \sum_{0 \leq i < n} q_{k, i} = 1 \quad & \text{for} \quad k \in \{0, 1, \ldots, n - 1\}, \\ & \sum_{0 \leq k < n} q_{k, i} = 1 \quad & \text{for} \quad i \in \{0, 1, \ldots, n - 1\}, \\ & q_{0, i} = q_{n, i} \quad & \text{for} \quad i \in \{0, 1, \ldots, n - 1\}, \\ & q_{k, i} \in \{0, 1\} & \end{align} \end{split}\]

定式化に入る前に、距離行列 distance (= \(d\)) を作成しておきます。今回は、distance はランダムな対称行列として生成します。

import numpy as np

NUM_CITIES = 100

# ランダムな対称行列を生成
distance = np.zeros((NUM_CITIES, NUM_CITIES))
for i in range(NUM_CITIES):
    for j in range(i+1, NUM_CITIES):
        distance[i, j] = distance[j, i] = np.random.rand()

注釈

実用的には各都市の座標を生成して距離行列を求めると良いでしょう。例えば scipy を用いて次のような方法で高速に計算できます。

from scipy.spatial import distance

locations = np.random.random((NUM_CITIES, 2)) # 座標をランダムに生成
distance = distance.cdist(locations, locations, metric='euclidean')

ナイーブな定式化

for 文やリスト内包表記を用いると、以下のように定式化することが可能です。

from amplify import Poly, VariableGenerator, one_hot

# 決定変数を作成
gen = VariableGenerator()
q = gen.array("Binary", NUM_CITIES + 1, NUM_CITIES)

# q[NUM_CITIES, i] = q[0, i] となるように設定
for i in range(NUM_CITIES):
    q[NUM_CITIES, i] = q[0, i]

# 目的関数を計算
objective = 0
for i in range(NUM_CITIES):
    for j in range(NUM_CITIES):
        for k in range(NUM_CITIES):
            objective += distance[i, j] * q[k, i] * q[k + 1, j]

# 制約条件の構築
row_one_hot_constraints = sum(
    one_hot(sum(q[i, k] for k in range(NUM_CITIES))) for i in range(NUM_CITIES)
)
col_one_hot_constraints = sum(
    one_hot(sum(q[i, k] for i in range(NUM_CITIES))) for k in range(NUM_CITIES)
)

# モデルの作成
model = objective + (row_one_hot_constraints + col_one_hot_constraints)

定式化の時間は以下のようになります (実行環境により異なります)。

目的関数	制約条件	全体
1.55 s	13.3 ms	1.56 s

改善その 1: ビルトインの sum 関数を使わない

上記のコードの最初の行に

from amplify import sum

を追加し、ビルトインの sum() 関数の代わりに amplify.sum() 関数を使用します。amplify.sum() 関数は Amplify SDK の提供するクラスの和に特化した関数です。すると、定式化の時間は以下のようになります (実行環境により異なります)。

目的関数	制約条件	全体
1.55 s	3.95 ms	1.55 s

amplify の amplify.sum() 関数を import するだけで、制約条件の定式化が約 3 倍高速化されました。全体の時間からの差分は小さいですが、巡回セールスマン問題の場合、制約条件の個数 (\(2n\)) が目的関数の項数 (\(O\left(n^3\right)\)) と比べて小さいためです。定式化に時間がかかる部分に sum() 関数が使われている場合、この最適化が非常に効果的になることがあります。

注釈

上記では from amplify import sum とすることにより sum() 関数をオーバーライドしています。Amplify SDK の amplify.sum() 関数は Amplify SDK の提供するクラス以外に対しては自動的にビルトインの sum() 関数にフォールバックするため、置き換えが可能になっています。

改善その 2: for 文やリスト内包表記を使わない

目的関数は for 文を使って定式化していましたが、Python の for 文は見た目が悪い上に低速です。できるだけ numpy や Amplify SDK の numpy-like な多項式配列である PolyArray のメソッドを使用して定式化することを試みます。

まず、\(q\) の上 \(n\) 行からなる \(n \times n\) 行列を \(A\), 下 n 行からなる \(n\times n\) 行列を \(B\) とおきます。すると目的関数は

\[ \sum_{0 \leq i, j, k < n} d_{i, j} q_{k, i} q_{k+1, j} = \sum_{0 \leq i, j, k < n} d_{i, j} A_{k, i} B_{k, j} \]

と書けます。\(A\) と \(d\) に着目すると、

\[ (Ad)_{k, j} = \sum_{0 \leq k, j < n} A_{k, i} d_{i, j} \]

なので、目的関数 \(\displaystyle \sum_{0 \leq i, j, k < n} d_{i, j} A_{k, i} B_{k, j}\) は 「(\(A\) と \(d\) の行列積) と \(B\) の要素積」の総和であることができます。したがって、目的関数は Amplify を用いて以下のように書けます。

q1 = q[:-1]
q2 = q[1:]
objective = ((q1 @ distance) * q2).sum()

また、制約条件は、one_hot() 関数に axis キーワード引数を渡すことで

row_one_hot_constraints = one_hot(q1, axis=1)
col_one_hot_constraints = one_hot(q1, axis=0)

と書けます。全体としては以下のようなコードになります。

from amplify import VariableGenerator, one_hot

# 決定変数を作成
gen = VariableGenerator()
q = gen.array("Binary", NUM_CITIES + 1, NUM_CITIES)

# q[NUM_CITIES, i] = q[0, i] となるように設定
q[-1, :] = q[0, :]

# q のスライスを作成
q1 = q[:-1]
q2 = q[1:]

# 目的関数を計算
objective = ((q1 @ distance) * q2).sum()

# 制約条件の構築
row_one_hot_constraints = one_hot(q1, axis=1)
col_one_hot_constraints = one_hot(q1, axis=0)

model = objective + (row_one_hot_constraints + col_one_hot_constraints)

この改善により、かかる時間は以下のようになります (実行環境により異なります)。

目的関数	制約条件	全体
152.1 ms	0.870 ms	154.4 ms

改善その 3: 難しく考えない

改善その 2 では、目的関数を式変形して numpy のメソッドをうまく使える形に持っていくことを試みましたが、実際に式変形を考えるのは大変です。Amplify は、このような関数をアインシュタインの縮約記法を用いて機械的に定式化することができる einsum() 関数をサポートしています。

まず、「改善その 2」の場合と同様に、\(q\) の上 \(n\) 行からなる \(n \times n\) 行列を \(A\), 下 n 行からなる \(n\times n\) 行列を \(B\) とおくことにより、目的関数は

\[ \sum_{0 \leq i, j, k < n} d_{i, j} q_{k, i} q_{k+1, j} = \sum_{0 \leq i, j, k < n} d_{i, j} A_{k, i} B_{k, j} \]

と書くことができます。この式に登場する行列 \(d\), \(A\), \(B\) の添え字 \(i\), \(j\), \(k\) が動く範囲は、それぞれが対応する行または列の長さと一致していることに注意してください。このような場合、einsum() 関数を適用すると目的関数の計算を以下のように書くことができます。

# 目的関数を計算
objective = einsum("ij,ki,kj->", distance, q1, q2)

ここで、einsum() 関数の第 1 引数はアインシュタインの縮約記法による積の添え字の並びを表し、"ij,ki,kj->" は、\(d_{i, j} A_{k, i} B_{k, j}\) という積を表しています。そして、第 2 引数以降はそれぞれの添え字に対応する配列を与えます。

全体としては以下のようなコードになります。

from amplify import VariableGenerator, one_hot, einsum

# 決定変数を作成
gen = VariableGenerator()
q = gen.array("Binary", NUM_CITIES + 1, NUM_CITIES)

# q[NUM_CITIES, i] = q[0, i] となるように設定
q[-1, :] = q[0, :]

# q のスライスを作成
q1 = q[:-1]
q2 = q[1:]

# 目的関数を計算
objective = einsum("ij,ki,kj->", distance, q1, q2)

# 制約条件の構築
row_one_hot_constraints = one_hot(q1, axis=1)
col_one_hot_constraints = one_hot(q1, axis=0)

# モデルの作成
model = objective + (row_one_hot_constraints + col_one_hot_constraints)

この定式化にかかる時間は次の通りです (実行環境により異なります)。

目的関数	制約条件	全体
58.11 ms	0.870 ms	60.4 ms

以上により、目的関数と全体の計算時間は 25倍以上、制約条件の構築は 15倍以上 の高速化が達成されました。