いきなりまとめ

scikit-learnのk-meansでクラスタリング
クラスタリングをイジングモデルで表現する
PyQUBO+Openjijでシミュレーテッドアニーリングを用いたクラスタリング

参考資料

クラスタリング問題とは？

クラスタリングとは、教師なし機械学習一つで、与えられたデータの集合をラベルなしにいくつかのクラスタリング（部分集合）へ分割することです。

データの準備

今回は明らかに分割できそうな、二次元平面上の点の集合をクラスタリングをしたいと思います。
次のように人工データを生成します。


# ライブラリをインポート
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()

出力が以下のようになります。（作り方が作り方なので、バキバキに分割できそうですね。）

k-means(k平均法)での実装

クラスタリングの例として、scikit-learnを用いたクラスタリングをしてみましょう。アルゴリズム等の詳細は参考資料をご覧ください。
それでは、上で作った人工データに対してscikit learnのk-meansでクラスタリングをしてみましょう。

from sklearn import cluster
# クラスタリング自体は一行で書ける。
kmeans= cluster.KMeans(n_clusters=2).fit(df1)

# 結果を表示
# クラスターごとに色を変えます
for idx, label in enumerate(kmeans.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")

結果は以下のようになりました。

きれいにクラスタリングできていますね（人工データのでそれはそう。）

イジングモデルによるクラスタリングの表現

今回のクラスタリングのハミルトニアン(最適化したい関数)は以下のように書けます。

$$
H = -\sum_{i,j\in D} \frac{1}{2}d_{ij}(1 - \sigma_i \sigma_j)
$$
$D$はデータの集合, $\sigma_i$, $\sigma_j$はスピン変数（-1,1のどちらかを取る）、$d_{ij}$は点$i,j$同士の距離を表します。
上のハミルトニアンでは、スピン変数$\sigma_i$(-1もしくは1)がその点$i$が所属するクラスを表します。したがって、ハミルトニアンの各項は２つスピン変数が同じであれば0、異なっていればその２つの点同士その距離$d_{ij}$になります。
つまり、ハミルトニアンを最小にすることは、「異なるクラスに所属する点同士の距離の和を最大化する」ことに対応しています。
なお、イジングモデルでのクラスタリングの最適化関数は先のk-meansでの最適化関数と異なっています。

PyQUBO+Openjijでの実装

今回は、PyQUBOでハミルトニアンを作り、SAはOpenjijで計算します。
(なお、PyQUBOにもシミュレーテッドアニーリングのソルバーがあるので、そちらを使っても良いかもしれません。)
どちらも、pip installできるのでアニーリングアルゴリズムに興味がある方は、お手軽に量子アニーリングをシミュレーションをするSQA(simulated quantum annealing)が実行できて便利です。
今回は、pandasのDataFrameを渡すと、クラスタリングしたラベルを返す関数を実装します。

# ライブラリをインポート
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のシミュレーテッド量子アニーリング(SQASampler)で解を求めたいので、サンプラーのインスタンス化
    sampler = oj.SQASampler(iteration=10, step_num=100)
    # 実際に解を求める
    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,  sum(response.energies)

では、実際に解を求めてみましょう。

labels, energy =clustering_openjij(df1[["x", "y"]])

結果を確認してみましょう。

print(labels)
>> [0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1]

可視化をしましょう。

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")

k-meansと同様にクラスタリングができていますね。

まとめ（再掲）

今回は以下のことを行いました。

二次元平面上で人工データを生成する
scikit-learnのk-meansでクラスタリング
クラスタリングをイジングモデルで表現する
PyQUBO+Openjijでシミュレーテッドアニーリングを用いたクラスタリング

シミュレーテッドアニーリングでクラスタリング問題を解く