アニーリングを用いたクラスタリング

このチュートリアルでは、アニーリングの応用の一例としてPyQUBOとOpenjijを利用したクラスタリングを取り上げます。

クラスタリング

クラスタリングとは与えられたデータを$n$個のクラスターに分けるというタスクです($n$は外部から与えられているとします)。簡単のため、今回はクラスター数が2を考えましょう。

クラスタリングのハミルトニアン

今回は、以下のハミルトニアンを最小化することでクラスタリングを行います。

$$H = - \sum_{i, j} \frac{1}{2}d_{i,j} (1 - \sigma _i \sigma_j)$$

$i, j$はサンプルの番号、$d_{i,j}$は２つのサンプル間の距離、$\sigma_i=\{-1,1\}$は２つのクラスターのどちらかに属しているかを表すスピン変数です。
このハミルトニアンの和の各項は

• $\sigma_i = \sigma_j$のとき、0

• $\sigma_i \neq \sigma_j$のとき、$d_{i,j}$

となります。右辺のマイナスに注意すると、ハミルトニアン全体では「異なるクラスに属しているサンプル同士の距離を最大にする$\{\sigma _1, \sigma _2 \ldots \}$の組を選びなさい」という問題に帰着することがわかります。

必要なライブラリのインポート

!pip install jijmodeling-transpiler

Requirement already satisfied: jijmodeling-transpiler in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (0.6.5)
Requirement already satisfied: jijmodeling>=1.0.0 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (1.0.3)
Requirement already satisfied: numpy<1.25.0,>=1.17.0 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (1.24.4)
Requirement already satisfied: typeguard in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (4.1.5)
Requirement already satisfied: pydantic in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (2.3.0)
Requirement already satisfied: mip in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (1.15.0)
Requirement already satisfied: dimod in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling-transpiler) (0.12.12)
Requirement already satisfied: pandas in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling>=1.0.0->jijmodeling-transpiler) (2.1.0)
Requirement already satisfied: orjson<4.0.0,>=3.8.0 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from jijmodeling>=1.0.0->jijmodeling-transpiler) (3.9.7)
Requirement already satisfied: cffi==1.15.* in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from mip->jijmodeling-transpiler) (1.15.1)
Requirement already satisfied: pycparser in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from cffi==1.15.*->mip->jijmodeling-transpiler) (2.21)
Requirement already satisfied: annotated-types>=0.4.0 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from pydantic->jijmodeling-transpiler) (0.5.0)
Requirement already satisfied: pydantic-core==2.6.3 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from pydantic->jijmodeling-transpiler) (2.6.3)
Requirement already satisfied: typing-extensions>=4.6.1 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from pydantic->jijmodeling-transpiler) (4.7.1)
Requirement already satisfied: importlib-metadata>=3.6 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from typeguard->jijmodeling-transpiler) (6.8.0)
Requirement already satisfied: zipp>=0.5 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from importlib-metadata>=3.6->typeguard->jijmodeling-transpiler) (3.16.2)
Requirement already satisfied: python-dateutil>=2.8.2 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from pandas->jijmodeling>=1.0.0->jijmodeling-transpiler) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from pandas->jijmodeling>=1.0.0->jijmodeling-transpiler) (2023.3.post1)
Requirement already satisfied: tzdata>=2022.1 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from pandas->jijmodeling>=1.0.0->jijmodeling-transpiler) (2023.3)
Requirement already satisfied: six>=1.5 in /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages (from python-dateutil>=2.8.2->pandas->jijmodeling>=1.0.0->jijmodeling-transpiler) (1.16.0)

# ライブラリのインポート
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
from scipy.spatial import distance_matrix

import openjij as oj
import jijmodeling as jm
from jijmodeling.transpiler.pyqubo.to_pyqubo import to_pyqubo

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[2], line 9
      7 import openjij as oj
      8 import jijmodeling as jm
----> 9 from jijmodeling.transpiler.pyqubo.to_pyqubo import to_pyqubo

ModuleNotFoundError: No module named 'jijmodeling.transpiler'

JijModelingとOpenJijを用いたクラスタリング

最初に、JijModelingを用いた上記のハミルトニアンを定式化します。JijModelingでは、スピン変数$\sigma_i$を扱えないので、バイナリ変数$x_i$で書けるように、$\sigma_i = 2x_i - 1$という関係を用いて書き直します。

problem = jm.Problem("clustering")
d = jm.Placeholder("d", dim=2)
N = d.shape[0].set_latex('N')
x = jm.Binary("x", shape=(N))
i = jm.Element("i", (0, N))
j = jm.Element("j", (0, N))
problem += (
    -1 / 2 * jm.Sum([i, j], d[i, j] * (1 - (2 * x[i] - 1) * (2 * x[j] - 1)))
)
problem

$$\begin{alignat*}{4}\text{Problem} & \text{: clustering} \\\min & \quad \left( -0.5 \right) \cdot \sum_{ i = 0 }^{ N - 1 } \sum_{ j = 0 }^{ N - 1 } d_{i,j} \cdot \left( 1 - \left( 2 \cdot x_{i} - 1 \right) \cdot \left( 2 \cdot x_{j} - 1 \right) \right) \\& x_{i_{0}} \in \{0, 1\}\end{alignat*}$$

人工データの生成

今回は人工的に二次元平面上の明らかに線形分離可能なデータを生成しましょう。

data = []
label = []
for i in range(100):
    # [0, 1]の乱数を生成
    p = np.random.uniform(0, 1)
    # ある条件を満たすときをクラス1、満たさないときを-1
    cls =1 if p>0.5 else -1
    # ある座標を中心とする正規分布に従う乱数を作成
    data.append(np.random.normal(0, 0.5, 2) + np.array([cls, cls]))
    label.append(cls)
# DataFrameとして整形
df1 = pd.DataFrame(data, columns=["x", "y"], index=range(len(data)))
df1["label"] = label

# データセットの可視化
df1.plot(kind='scatter', x="x", y="y")

<AxesSubplot: xlabel='x', ylabel='y'>

instance_data = {'d':distance_matrix(df1, df1)}

OpenJijを使ってクラスタリング問題を解く

数理モデルとデータができたので、早速、Openjijを使って問題を解いていきましょう。

pyq_obj, pyq_cache = to_pyqubo(problem, instance_data, {})
qubo, constant = pyq_obj.compile().to_qubo()
sampler = oj.SASampler()
response = sampler.sample_qubo(qubo)
result = pyq_cache.decode(response)

可視化をしてみましょう。

for idx in range(0, len(instance_data['d'])):
    if idx in result.record.solution["x"][0][0][0]:
        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")
plt.xlabel('x')
plt.ylabel('y')

Text(0, 0.5, 'y')

赤と青のクラスに分類できていることがわかります。