シミュレーティド・アニーリングでクラスタリングのベンチマークを取得

この記事は、株式会社Jijでのインターンシップ中に作成した記事です。
[会社のホームページ j-ij.com](https://j-ij.com/)

# この記事の内容
- PyQUBOでクラスタリングのハミルトニアン→QUBO変換
- OpenJijのSQAsamplerでクラスタリングを実装
- クラスタリングのベンチマークを取得(SA/SQA)

# 参考資料
- [シュミレーテッドアニーリングでクラスタリング問題を解く](https://qiita.com/akiFQC/items/05671dc82fac7f818330)

# ライブラリのインストール
PyQUBO, OpenJijがない人は、インストールしましょう。
※OpenJijは、ver0.0.7とver0.0.8でベンチマーク関数の仕様変更がありました。古いバージョンの方は、アップデートしましょう。

```
pip install openjij  # ver0.0.8
pip install pyqubo
```

# データの準備
平均の異なるガウス分布により、明らかに分割できそうな点の集合を用意します。
>
```
# ライブラリをインポート
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
>
# 人工データの生成
data = []
for i in range(100): # 100個のデータを生成
    p = np.random.uniform(0, 1)
    cls =1 if p>0.5 else -1 #点の配置は平均の異なる２つのガウス分布のどちらかに従う。
    data.append(np.random.normal(0, 0.5, 2) + np.array([cls, cls]))
>
df1 = pd.DataFrame(data, columns=["x", "y"], index=range(len(data)))
>
# 可視化
df1.plot(kind='scatter', x="x", y="y")
plt.show()
```
＜[シュミレーテッドアニーリングでクラスタリング問題を解く](https://qiita.com/akiFQC/items/05671dc82fac7f818330)　より引用＞

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/496581/f56c526d-fddc-2a6f-187c-4e18bdfcba31.png)

# PyQUBO+Openjijでの実装
参考にした記事の手法で、クラスタリングができているか確認します。

```
# ライブラリをインポート
from pyqubo import Array
import openjij as oj
from scipy.spatial import distance_matrix

# 関数を定義
def clustering_openjij(df):

    # 距離行列(d_ij)(各サンプル同士の距離を要素に持つ行列)を計算
    d_ij = distance_matrix(df, df)

    # spinを定義
    spin = Array.create("spin", shape= len(df), vartype="SPIN")

    # ハミルトニアンを定義
    H = - 0.5* sum(
        [d_ij[i,j]* (1 - spin[i]* spin[j]) for i in range(len(df)) for j in range(len(df))]
    )
    # ハミルトニアンをコンパイル
    model = H.compile()
    # コンパイルしたモデルをqubo化
    qubo, offset = model.to_qubo()
    # openjijのシュミレーテッドアニーリング(SASampler)で解を求めたいので、サンプラーのインスタンス化
    step_num_list = [20, 40, 60, 80, 100]
    iteration = 10
    
    for step_num in step_num_list:
        sampler = oj.SASampler(beta_max=10.0, beta_min=0.01, step_num=step_num, iteration=iteration)
    
        # 実際に解を求める
        response = sampler.sample_qubo(qubo)
        # pyquboのdecode_solutionを使いたいのでラベルと解を合わせて辞書化する。
        # openjijのsamplerではいくつか解が返ってくるので、response.energiesが最小のものを最適解として採用する。
        raw_solution = dict(zip(response.indices, response.states[np.argmin(response.energies)]))
        # pyquboのdecoded_solutionでデコードされた解、壊れたスピン、エネルギーを得る
        decoded_solution, broken, energy= model.decode_solution(raw_solution, vartype="SPIN")
        # 最終的にほしい各点ごとのラベルのリスト
        labels = [int(decoded_solution["spin"][idx] ) for idx  in range(len(df))]
    return labels, response.states[np.argmin(response.energies)], response.energies
```

クラスタリング用のラベル、最適解のスピン配列、残留エネルギーを表示します。

```
labels, correct_state, energies =clustering_openjij(df1[["x", "y"]])

print('Optimal Labels :  ', labels)  # クラスタリングに必要な代表ラベル
print('Correct State :  ', correct_state)  # 代表ラベルを与えるスピン配列
print('energies :  ', energies)  # サンプリング10回分の各残留エネルギー
```

残留エネルギーの値が収束しているのが確認できるので、、解が最適解までたどり着いたとみなしてよいでしょう。
よって、最小エネルギーを取るスピン配列を取得し、プロットに必要なラベルを計算しています。
では、実際にプロットしてみましょう。

```
for idx, label in  enumerate(labels):
    if label:
        plt.scatter(df1.loc[idx]["x"], df1.loc[idx]["y"], color="b")  
    else:
        plt.scatter(df1.loc[idx]["x"], df1.loc[idx]["y"], color="r")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/496581/4683956b-c8f5-7520-2443-278bec4713b3.png)

確かに、クラスタリングできていますね！！

# ベンチマーク取得
OpenJijにはデフォルトでTTS(Time To Solution)、残留エネルギー、成功確率を評価してくれる関数が付いています。
これを使って、SAsampler のベンチマークを取得します。

```
# 最適解は、先ほど計算したので、、
correct_state = correct_state

# ステップ数とアニーリングの反復数を与えます
step_num_list = list(range(20, 101, 20))  #[20, 40, 60, 80, 100]
iteration = 10

# 距離行列(d_ij)(各サンプル同士の距離を要素に持つ行列)を計算
df = df1[["x", "y"]]
d_ij = distance_matrix(df, df)

# spinを定義
spin = Array.create("spin", shape= len(df), vartype="SPIN")

# ハミルトニアンを定義
H = - 0.5* sum(
    [d_ij[i,j]* (1 - spin[i]* spin[j]) for i in range(len(df)) for j in range(len(df))]
    )
# ハミルトニアンをコンパイル
model = H.compile()
# コンパイルしたモデルをqubo化
qubo, offset = model.to_qubo()

# benchmark 関数で TTS 残留エネルギー 成功確率を計算。
result = oj.solver_benchmark(
                      solver= lambda time, **args: oj.SASampler(step_num=time, iteration=iteration).sample_qubo(qubo), 
                      time_list=step_num_list, solutions=[correct_state], p_r=0.99
            )
```

TTS(Time To Solution)、残留エネルギー、成功確率はそれぞれ次のようになります。

```
fig, (axL,axC,axR) = plt.subplots(ncols=3, figsize=(15,3))
plt.subplots_adjust(wspace=0.4)

fontsize = 10
axL.plot(result['time'], result['tts'])
axL.set_xlabel('annealing time', fontsize=fontsize)
axL.set_ylabel('TTS', fontsize=fontsize)

axC.plot(result['time'], result['residual_energy'])
axC.set_xlabel('annealing time', fontsize=fontsize)
axC.set_ylabel('Residual energy', fontsize=fontsize)

axR.plot(result['time'], result['success_prob'])
axR.set_xlabel('annealing time', fontsize=fontsize)
axR.set_ylabel('Success probability', fontsize=fontsize)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/496581/c373b79d-e856-0890-be88-4a14bb5d2c10.png)

TTSと成功確率が増加傾向になりましたね。
残留エネルギーについては、ほとんど同じ値が10回分出力されていましたので、確かにこのグラフの通りになりそうです。
同様のやり方で、シミュレーティド・量子アニーリング(SQA)などのベンチマークも取得できますので、試しに取ってみます。

```
# benchmark 関数で TTS 残留エネルギー 成功確率を計算。
result = oj.solver_benchmark(
                      solver= lambda time, **args: oj.SQASampler(step_num=time, iteration=iteration).sample_qubo(qubo), 
                      time_list=step_num_list, solutions=[correct_state], p_r=0.99
```

```
fig, (axL,axC,axR) = plt.subplots(ncols=3, figsize=(15,3))
plt.subplots_adjust(wspace=0.4)

fontsize = 10
axL.plot(result['time'], result['tts'])
axL.set_xlabel('annealing time', fontsize=fontsize)
axL.set_ylabel('TTS', fontsize=fontsize)

axC.plot(result['time'], result['residual_energy'])
axC.set_xlabel('annealing time', fontsize=fontsize)
axC.set_ylabel('Residual energy', fontsize=fontsize)

axR.plot(result['time'], result['success_prob'])
axR.set_xlabel('annealing time', fontsize=fontsize)
axR.set_ylabel('Success probability', fontsize=fontsize)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/496581/2eb4936b-59e6-93e1-30b7-d5e86fa4277c.png)

このように、クラスタリングアルゴリズムについてSAやSQAでベンチマークを取得することができました！
ただし、今回のクラスタリングアルゴリズムに用いた最適化関数（ハミルトニアン）は、k-meansに用いられる最適化関数と異なっています。
実際に、従来のソルバーと比較する際には、注意が必要です。