第3回: アニーリング基礎3. クラスタリング問題、2次元スピンの1次元化

• 12人を2グループに分けたい。それぞれの、血液型、出身地に基づいて、できる限り同じ種類の人が同じグループに集まるように分類する。
• これはクラスタリング問題であるが、アニーリングを用いて解いてみる。
• この問題ではスピンを$x_{ig}$といった2次元で示すことになるが、アニーリングで解くためには1次元化が必要。その方法をPyquboを使用した場合、直接的に表現した場合の双方で検証してみる。

前回: 第2回: アニーリング基礎2.ナップザック問題、Optunaによるハイパーパラメタ最適化
参考: 東北大学大関先生講義 第3回, ソースコード
Keyword ：クラスタリング

1. 基本準備/設定

1.1 量子アニーリング用モジュールのインストール

まず、初めに、以下のモジュールをインストールする。

• dwave-ocean-sdk: Dwaveアニーリングマシンのモジュール
• openjij: アニーリングシュミレーター
• pyqubo: Qubo行列作成サポートモジュール
• optuna: ハイパーパラメタの最適化

pip install dwave-ocean-sdk
pip install openjij
pip install pyqubo
pip install optuna

1.2 Dwaveの初期設定

続いて、Dwaveのtoken、enpoint(アクセスURL)、使用するマシンの指定を行う。

token = "DEV-73ac79954ef2135ada52dXXXXXXXXXXXXXXXXXX"
endpoint = "https://cloud.dwavesys.com/sapi/"
solver = "Advantage_system4.1"

更に、Dwaveマシンを指定するDwaveSamplerと、Qbitへの割り当て行うEmbeddingCompositeをインポート。続いて、DWaveSamplerを用いてDwaveマシンを指定する。尚、Dwaveの動かしかた概要は、Qiitaの記事を参考：https://qiita.com/y_daito/items/aad74517927e393564ee

from dwave.system import DWaveSampler, EmbeddingComposite
dw_sampler = DWaveSampler(solver=solver, token=token, endpoint = endpoint)

1.3 各種モジュールのインポート

openjijからSQASampler、pyquboからArray、その他としてnumpy、pandas、matplotlib.pyplotをインポートする。

from openjij import SQASampler
from pyqubo import Array, Placeholder
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

これで下準備は完了。

2. 定式化

2.1 問題の定義から相性の整理まで

今回グループ分け行う人数$N=12人$で、これをグループ数$G=2$に分割する。それぞれの血液型(A or B or AB or O)、出身地(関西 or 関東)を以下とする。

N = 12
G = 2
data = [[1, "A", "関東"],[2,"A", "関西"],[3, "AB", "関東"], [4, "O", "関西"], [5, "B", "関西"], [6, "AB", "関東"],[7, "AB", "関東"],[8,"O", "関東"],[9, "B", "関東"], [10, "A", "関西"], [11, "AB", "関西"], [12, "O", "関東"]]

12人のデータは文字データであり、そのまま加工できないので、数値化が必要。血液データはA、B、AB、Oの4種類あるため、Aであれば[1,0,0,0]、Bであれば[0,1,0,0]と表現する。また、出身地は関東であれば[1,0]、関西であれば[0,1]とする。尚、これはOne-hot表示や二値化と呼ばれる。まず、二値化を行う関数を定義する。

### 血液型の二値化する関数の定義
def blood_01(blood):
  blood_list = ["A","B","AB","O"]
  _output = [0]*len(blood_list)
  _output[blood_list.index(blood)]=1
  return(_output)

### 出身地を二値化する関数の定義
def local_01(local):
  local_list = ["関東","関西"]
  _output = [0]*len(local_list)
  _output[local_list.index(local)]=1
  return(_output)

続いて、Dateの血液型、出身地を二値化し、それぞれをリストで連結し、6個のデータで構成されるリストに変換する。

data_01 = [blood_01(_data[1])+local_01(_data[2]) for _data in data]
data_01

<Output>
[[1, 0, 0, 0, 1, 0],
 [1, 0, 0, 0, 0, 1],
 [0, 0, 1, 0, 1, 0],
 [0, 0, 0, 1, 0, 1],
 [0, 1, 0, 0, 0, 1],
 [0, 0, 1, 0, 1, 0],
 [0, 0, 1, 0, 1, 0],
 [0, 0, 0, 1, 1, 0],
 [0, 1, 0, 0, 1, 0],
 [1, 0, 0, 0, 0, 1],
 [0, 0, 1, 0, 0, 1],
 [0, 0, 0, 1, 1, 0]]

例えば、1行目の[1,0,0,0,1,0]は左から4データ[1,0,0,0]は血液型のA型を意味し、右の2データ[1,0]は出身地の関東を意味する。
今回似ている人同士をまとめてグループにする。従い、それぞれが似てるかどうかを表す指標が必要となる。夫々のデータを掛け合わせると、0、1、2のいずれかの数値になるが、これは同じ属性が何個あるかを示している。
例えば、1行目の[1,0,0,0,1,0]を二乗すると2になるが、いずれもA型で関東出身であることを意味する。12人全員の相性$C$(CompatibilityのC)を確認するべく、先ほどのdata_01の内積をとる。ただし、先ほどの例のように自己の相性は存在しない為、対角成分は全て0にする。

C = np.dot(np.array(data_01),np.array(data_01).T)#正方行列でないので、内積のもう一方は転置した。
for i in range(len(C)):
  C[i][i] = 0
C

<Output>
array([[0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1],
       [1, 0, 0, 1, 1, 0, 0, 0, 0, 2, 1, 0],
       [1, 0, 0, 0, 0, 2, 2, 1, 1, 0, 1, 1],
       [0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1],
       [0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0],
       [1, 0, 2, 0, 0, 0, 2, 1, 1, 0, 1, 1],
       [1, 0, 2, 0, 0, 2, 0, 1, 1, 0, 1, 1],
       [1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 2],
       [1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1],
       [1, 2, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0],
       [0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0],
       [1, 0, 1, 1, 0, 1, 1, 2, 1, 0, 0, 0]])

エネルギー式の作成に於いては、Cを用いた相性が最大化するような数式にすればよい。一方で、定式化はスピンの設定を行ってから後で考える。

2.2 コスト関数の作成、エネルギー式の作成、ハミルトニアンの作成

今回のグループ分けは$N=12$人を、$G=2$に分類する。表でまとめると以下のようになる。

$x_{ig}$	1さん(i=1)	2さん(i=2)	3さん(i=3)	...
Group 1 (g=1)		1		...
Group 2 (g=2)	1		1	...

この場合、スピンは$x_{ig}$となり、$i$さんが、$g$グループに所属する場合を1、しない場合を0として考える。

さて、コスト関数であるが、以下が考えられる。

• それぞれが所属できるグループは1件のみである。
• 各グループの人数は6名である。

まず、1点目から定式化する。上の表の列の合計は必ず1になることを意味する。

C_1(x) = \lambda_1\sum_{i=1}^N(\sum_{g=1}^Gx_{ig}-1)^2

続いて、2点目を定式化する。上の表の行の合計は必ず6になることを意味する。

C_2(x) = \lambda_2\sum_{g=1}^G(\sum_{i=1}^Nx_{ig}-6)^2

最後にエネルギー式の定式化を行う。今回、似た者同士が集まればよい。すると、各グループの相性値が最大化されれば良い。従い、エネルギー式は以下とできる。

E(x) = -\sum_{g=1}^G(\sum_{i=1}^N\sum_{j=1}^NC_{ij}x_{ig}x_{jg})

従い、ハミルトニアンは以下となる。

H(x) = E(x) + C_{1}(x) + C_{2}(x) \\
H(x) = -\sum_{g=1}^G\sum_{i=1}^N\sum_{j=1}^NC_{ij}x_{ig}x_{jg} + \lambda_1\sum_{i=1}^N(\sum_{g=1}^Gx_{ig}-1)^2 + \lambda_2\sum_{g=1}^G(\sum_{i=1}^Nx_{ig}-6)^2

さて、ここまでスピン$x_{ig}$と二つの添え字で表現されてきたことを、さも当たり前に扱ってきた。しかし、これは$x$が$i$と$g$の2次元になっている。スピンそれ自体は1行で記載されている必要がある。Pyquboではスピンが複数となった場合でもこれまでと同様に極めて簡単に処理できるようになっている。
実装ではPyquboで計算を行う。最後に、勉強がてらにPyquboを使用せずにQubo行列を作成した場合も作成してみる。

3. 実装

3.1 Pyquboを用いたQubo行列の作成と、λ1、λ2を仮置きしたDwaveでの計算実行

まず、スピンの設定を行う。スピンが多次元である場合、shape=(x,y,z,…)と設定することで難なくスピンの設定が完了する。

### ライブラリのインポート
from pyqubo import Array, Placeholder

### コスト関数の影響度の設定
lambda1, lambda2 = Placeholder("lambda1"), Placeholder("lambda2")

### スピンの設定
x = Array.create(name="x",shape=(N,G), vartype = "BINARY")
x

<Output>
Array([[Binary('x[0][0]'), Binary('x[0][1]')],
       [Binary('x[1][0]'), Binary('x[1][1]')],
       [Binary('x[2][0]'), Binary('x[2][1]')],
       [Binary('x[3][0]'), Binary('x[3][1]')],
       [Binary('x[4][0]'), Binary('x[4][1]')],
       [Binary('x[5][0]'), Binary('x[5][1]')],
       [Binary('x[6][0]'), Binary('x[6][1]')],
       [Binary('x[7][0]'), Binary('x[7][1]')],
       [Binary('x[8][0]'), Binary('x[8][1]')],
       [Binary('x[9][0]'), Binary('x[9][1]')],
       [Binary('x[10][0]'), Binary('x[10][1]')],
       [Binary('x[11][0]'), Binary('x[11][1]')]])

続いて、先ほど整理した以下のハミルトニアンを実装する。

\begin{equation}
H(x) = -\sum_{g=1}^G\sum_{i=1}^N\sum_{j=1}^NC_{ij}x_{ig}x_{jg} + \lambda_1\sum_{i=1}^N(\sum_{g=1}^Gx_{ig}-1)^2 + \lambda_2\sum_{g=1}^G(\sum_{i=1}^Nx_{ig}-6)^2
\end{equation}

###第1項の実装
E = 0
for g in range(G):
  for i in range(N):
    for j in range(N):
      E = E - C[i][j]*x[i,g]*x[j,g]

###第2項実装
C1 = 0
for i in range(N):
  _c = 0
  for g in range(G):
    _c = _c + x[i,g]
  C1 = C1 + lambda1*(_c -1)**2

###第3項実装
C2 = 0
for g in range(G):
  _c = 0
  for i in range(N):
    _c = _c + x[i,g]
  C2 = C2 + lambda2*(_c - 6)**2

###ハミルトニアンの実装
H = E + C1 + C2

続いて、ハミルトニアンからQubo行列を取り出し、結果を確認してみる。$\lambda_1,\lambda_2$は仮でそれぞれ100として計算してみる。

###Qubo行列の取り出し
model = H.compile()
Qubo,offset = model.to_qubo(feed_dict={"lambda1":100,"lambda2":100})

###Dwaveで計算
sampler_dw = EmbeddingComposite(dw_sampler)
solution = sampler_dw.sample_qubo(Qubo, num_reads=10)
solution.record

<Output>
rec.array([([1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1], -8444., 1, 0.        ),
           ([1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1], -8444., 1, 0.        ),
           ([1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1], -8438., 1, 0.        ),
           ([1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1], -8434., 1, 0.        ),
           ([1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1], -8240., 1, 0.        ),
           ([1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0], -8238., 1, 0.        ),
           ([1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1], -8236., 1, 0.        ),
           ([1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0], -8234., 1, 0.        ),
           ([1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1], -8238., 1, 0.04166667),
           ([1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1], -8050., 1, 0.04166667)],
          dtype=[('sample', 'i1', (24,)), ('energy', '<f8'), ('num_occurrences', '<i8'), ('chain_break_fraction', '<f8')])

3.2 算出結果の評価方法の準備と、評価

得られた結果は、以下の形式となっている。
[1さんがグループ1である、1さんがグループ2である、2さんがグループ1である、2さんがグループ2である、…]
従い、グループ1、2それぞれに誰が所属するかを抽出する関数を作成する。

def answer(solution, i =0):
###奇数番目の結果をG1のリストに、偶数番目の結果をG2のリストに格納
  g1, g2 =[], []
  for num, x in enumerate(solution.record[i][0]):
    if not(num+1)%2==0:
      g1.append(x)
    else:
      g2.append(x)
  result = [g1,g2]
  return result

###テスト
answer(solution,1)

<Output>
[[1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0], [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]]

続いて、G1、G2に夫々どういう特徴の人がグループ分けされたか可視化する関数を定義する。

def view(solution,i =1):
  _answer = answer(solution,i)
###G1の結果を取り出し、その結果からG1、G2を分類するとともに、元データから夫々の特徴をリストで並べる。
  g1_detail, g2_detail = [], []
  for num, x in enumerate(_answer[0]):
    _temp = [data[num][1],data[num][2]]
    if x == 1:
      g1_detail.append(_temp)
    elif x == 0:
      g2_detail.append(_temp)
  result = [g1_detail,g2_detail]

###関数はPrintのみで、計算結果は返さないとした。
  print(f"{g1_detail}\n{g2_detail}")

###テスト
view(solution,1)

<Output>
[['A', '関東'], ['AB', '関東'], ['B', '関西'], ['AB', '関東'], ['B', '関東'], ['AB', '関西']]
[['A', '関西'], ['O', '関西'], ['AB', '関東'], ['O', '関東'], ['A', '関西'], ['O', '関東']]

更に得られた結果が制約条件を満たすか堂か、更にグループ毎の相性合計を算出する関数を定義する。

def evaluation(solution,tar =1):
  _answer = np.array(answer(solution,tar))
  _answer_t = _answer.T

  ###夫々が所属するグループ数が1件となっているかの確認
  #夫々が所属するグループ数を抽出
  _result_x = np.sum(_answer,axis=0) #<Output>: array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

  #抽出したグループ数の内、グループ数が1のものを0、それ以外を1に返還
  _result_x_conclusion = np.where(_result_x==1 ,0 ,1) #<Output>: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

  #上記のArray合計が0であれば条件達成。達成か堂かをok or errorで返す。
  if np.sum(_result_x_conclusion) == 0:
    _distribution = 1 #okの意
  else:
    _distribution = 0 #errorの意

  ###グループの人数が均等になっているかの確認
  _result_y = np.sum(_answer,axis=1)
  _result_y_conclusion = np.where(_result_y==6 ,0 ,1)
  #上記のArray合計が0であれば条件達成。達成か堂かをok or errorで返す。
  if np.sum(_result_y_conclusion) == 0:
    _group = 1 #okの意
  else:
    _group = 0 #errorの意

  ###Groupの相性合計を計算
  #コスト関数をそのまま流用し、グループ1、グループ2の夫々で計算する。
  g1,g2 = 0,0
  for g in range(G):
    for i in range(N):
      for j in range(N):
        if g == 0:
          #xigは_answerを転置したものに置き換えている。
          g1 = g1 + C[i][j]*_answer_t[i,g]*_answer_t[j,g]
        elif g == 1:
          g2 = g2 + C[i][j]*_answer_t[i,g]*_answer_t[j,g]
  _evaluation =[str(tar)+":",_distribution, _group, g1,g2]
  return _evaluation

###テスト
evaluation(solution,1)

<Output>
['1:', 1, 1, 22, 20]

それでは、仮計算結果を作成した関数evaluationに投げてみる。

for i in range(10):
  print(evaluation(solution,i))

<Output>
['0:', 1, 1, 20, 18]
['1:', 1, 1, 22, 20]
['2:', 1, 1, 20, 18]
['3:', 1, 1, 20, 22]
['4:', 0, 1, 16, 18]
['5:', 0, 0, 12, 26]
['6:', 0, 0, 18, 10]
['7:', 0, 1, 22, 16]
['8:', 0, 1, 18, 18]
['9:', 0, 0, 22, 26]

条件を満たしていない計算結果がいくつか確認される。それらを確認してみる。

まず$i=2$を確認する。

view(solution,2)
answer(solution,2)

<Output>
[['A', '関東'], ['A', '関西'], ['AB', '関東'], ['O', '関西'], ['AB', '関東'], ['AB', '関西']]
[['B', '関西'], ['AB', '関東'], ['O', '関東'], ['B', '関東'], ['A', '関西'], ['O', '関東']]
[[1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0], [0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1]]

続いて、Optunaを用いて$\lambda_1$, $\lambda_2$を最適化することで、より良い結果が得られるか確認する。

3.3 Optunaによるλ1,λ2の最適化

まず始めに、アニーリングによる計算を実行する関数を定義する。関数ではDwave以外に、Jijシュミレーターも選択できるようにする。

def a_cal(machine = "Dwave", case =100):
### Dwaveかどうかの選択
  if machine == "Dwave":
    sampler_dw = EmbeddingComposite(dw_sampler)
    solution = sampler_dw.sample_qubo(Qubo, num_reads=case)
    
### Jijかどうかの選択
  elif machine =="Jij":
    sampler_jij = SQASampler()
    solution = sampler_jij.sample_qubo(Qubo,num_reads=case)
  return(solution)

続いて、Optunaを用いてハイパーパラメタの最適化を実施する。詳細は、第2回: アニーリング基礎2.ナップザック問題、Optunaによるハイパーパラメタ最適化を参照。

###ライブラリの呼び出し
import optuna

###関数の定義
def objective(trial):

  ###ハイパーパラメタの設定
  _lambda1, _lambda2 = trial.suggest_uniform('lambda1', 100, 10000), trial.suggest_uniform('lambda2', 100, 10000)

  ###コンパイルとQuboの取り出し
  model = H.compile()
  Qubo, offset = model.to_qubo(feed_dict={"lambda1":_lambda1,"lambda2":_lambda2})

  ###a_calでアニーリング計算実施し、それをrateで正答率にする。
  _scoure = 0
  _solution = a_cal("Jij",30)
  for i in range(30):
    _temp = evaluation(_solution,i)
    _scoure = _scoure + _temp[1]*_temp[2]*_temp[3]*_temp[4]
    _scoure = -1*_scoure
  return _scoure

### Optunaでの最適化
study = optuna.create_study()
study.optimize(objective, n_trials=40)

### Optunaでの最適解を表示
best_lambda1, best_lambda2= [i for i in study.best_params.values()][0], [i for i in study.best_params.values()][1]
print(f"lambda1: {best_lambda1}, lambda2: {best_lambda2}")

<Output>
lambda1: 9350.162917124419, lambda2: 1314.048521265912

得られた最適化なパラメタを用いてハミルトニアンを再度コンパイルし、Dwaveを用いて計算してみる。

model = H.compile()
Qubo, offset = model.to_qubo(feed_dict={"lambda1":best_lambda1,"lambda2":best_lambda2})
solution = a_cal("Dwave")

続いて、計算結果を評価する。

for i in range(10):
  print(evaluation(solution,i))

<Output>
['0:', 1, 1, 16, 20]
['1:', 1, 1, 20, 18]
['2:', 1, 1, 18, 16]
['3:', 1, 1, 22, 20]
['4:', 1, 1, 22, 20]
['5:', 1, 1, 22, 20]
['6:', 1, 1, 18, 18]
['7:', 1, 1, 20, 18]
['8:', 1, 1, 18, 22]
['9:', 1, 1, 22, 18]

例えば、$i=2$を確認してみる。

view(solution,2)

<Output>
[['A', '関西'], ['O', '関西'], ['AB', '関東'], ['AB', '関東'], ['B', '関東'], ['O', '関東']]
[['A', '関東'], ['AB', '関東'], ['B', '関西'], ['O', '関東'], ['A', '関西'], ['AB', '関西']]

上記となった。
さて、Dwaveから取得する結果には、2列目にハミルトニアンの値が含まれる。
Optunaを用いてパラメタ設定を行う時、わざわざ最適化の目的関数を設定してきたが、ハミルトニアンの値を最小化すればよいはず。

上記仮定が正しいかどうか検証してみる。

###ライブラリの呼び出し
import optuna

###関数の定義
def objective(trial):

  ###ハイパーパラメタの設定
  _lambda1, _lambda2 = trial.suggest_uniform('lambda1', 100, 10000), trial.suggest_uniform('lambda2', 100, 10000)

  ###コンパイルとQuboの取り出し
  model = H.compile()
  Qubo, offset = model.to_qubo(feed_dict={"lambda1":_lambda1,"lambda2":_lambda2})

  ###a_calでアニーリング計算実施し、それをrateで正答率にする。
  _scoure = 0
  _solution = a_cal("Jij",30)
  for result in _solution.record:
    _scoure = _scoure + result[1]
    _scoure = _scoure /10000
  return _scoure

### Optunaでの最適化
study = optuna.create_study()
study.optimize(objective, n_trials=20)

### Optunaでの最適解を表示
test_lambda1, test_lambda2= [i for i in study.best_params.values()][0], [i for i in study.best_params.values()][1]
print(f"lambda1: {test_lambda1}, lambda2: {test_lambda2}")

<Output>
lambda1: 3859.6742293944476, lambda2: 1563.7173749286092

得られた最適化なパラメタを用いてハミルトニアンを再度コンパイルし、Dwaveを用いて計算してみる。

model = H.compile()
Qubo, offset = model.to_qubo(feed_dict={"lambda1":test_lambda1,"lambda2":test_lambda2})
solution = a_cal("Dwave")

続いて、計算結果を評価する。

for i in range(10):
  print(evaluation(solution,i))

<Output>
['0:', 1, 1, 20, 20]
['1:', 1, 1, 16, 18]
['2:', 1, 1, 18, 16]
['3:', 1, 1, 26, 22]
['4:', 1, 1, 20, 22]
['5:', 1, 1, 22, 18]
['6:', 1, 1, 16, 18]
['7:', 1, 1, 18, 16]
['8:', 1, 1, 20, 16]
['9:', 1, 1, 22, 26]

概ね同様の精度が確認された。まず始めに、ハミルトニアンを使用し、Workし無い場合、個別に設定して進めればよい。

4. Pyquboを使用せずに2次元スピンを1次元化する方法

Pyquboにより2次元スピンを1次元に変換できることを確認した。ここでは、理解促進の為に、Pyquboを使用しない場合の方法を整理する。

4.1 定式化 / ハミルトニアンの展開

既に整理済の以下のハミルトニアンを展開する。

H(x) = -\sum_{g=1}^G\sum_{i=1}^N\sum_{j=1}^NC_{ij}x_{ig}x_{jg}+
\lambda_1\sum_{i=1}^N(\sum_{g=1}^Gx_{ig}-1)^2+
\lambda_2\sum_{g=1}^G(\sum_{i=1}^Nx_{ig}-6)^2

まず、第2項を展開する。

H(x)_{2nd}=
\lambda_1\sum_{i=1}^N
(\sum_{g=1}^Gx_{ig}-1)(\sum_{l=1}^Gx_{il}-1)
\\
H(x)_{2nd}=
\lambda_1\sum_{i=1}^N\sum_{g=1}^G\sum_{l=1}^Gx_{ig}x_{il}-
2\lambda_1\sum_{i=1}^N\sum_{g=1}^Gx_{ig}x_{ig}+
\lambda_1

続いて、第3項を展開する。

H(x)_{3rd}=
\lambda_2\sum_{g=1}^G
(\sum_{i=1}^Nx_{ig}-6)(\sum_{j=1}^Nx_{jg}-6)
\\
H(x)_{3rd}=
\lambda_2\sum_{g=1}^G\sum_{i=1}^N\sum_{j=1}^Nx_{ig}x_{jg}-
12\lambda_2\sum_{g=1}^G\sum_{i=1}^Nx_{ig}x_{ig}+
36\lambda_2

全てを合計すると、

H(x)=-\sum_{g=1}^G\sum_{i=1}^N\sum_{j=1}^NC_{ij}x_{ig}x_{jg}+
\lambda_1\sum_{i=1}^N\sum_{g=1}^G\sum_{l=1}^Gx_{ig}x_{il}-
2\lambda_1\sum_{i=1}^N\sum_{g=1}^Gx_{ig}x_{ig}+
\lambda_2\sum_{g=1}^G\sum_{i=1}^N\sum_{j=1}^Nx_{ig}x_{jg}-
12\lambda_2\sum_{g=1}^G\sum_{i=1}^Nx_{ig}x_{ig}+
\lambda_1+
36\lambda_2

続いて、Qubo行列を作成していくが、$x_{ig}$の扱いを考える。$i$は12人、$g$は2グループ迄をとる。仮にグループが1グループだけとすると(問題が成立しないが…)、Qubo行列は$12×12$の行列となる。一方で、問題は2グループなので、Qubo行列は$12×12×2$の3階のテンソルとなる。しかし、Qubo行列は2次式であり、2階以上のテンソルは対応不可。そこで、無理矢理、2階のテンソル、つまり行列に変換しないといけない。そこで、$i×g=24$を行もしくは列とする、$24×24$の行列とすればよい。例えば、2行、13列目は、2人目でｸﾞﾙｰﾌﾟ1に属す人と、1人目でｸﾞﾙｰﾌﾟ2に属す人の相互関係を意味する。それでは、実際に実装してみる。

4.2 実装

早速、実装を行う。

### Qubo行列の空箱の作成
_qubo = np.zeros((N*G)**2).reshape(N*G,N*G)
### lambda1、lambda2の設定
_lambda1 = 100000
_lambda2 = 10000

### 第1項
for g in range(G):
  for i in range(N):
    for j in range(N):
      _qubo[i+N*g,j+N*g] =_qubo[i+N*g,j+N*g] - C[i][j]

### 第2項
for i in range(N):
  for g in range(G):
    for l in range(G):
      _qubo[i+N*g,i+N*l] = _qubo[i+N*g,i+N*l] + _lambda1

### 第3項
for i in range(N):
  for g in range(G):
    _qubo[i+N*g,i+N*g] = _qubo[i+N*g,i+N*g] - 2*_lambda1

### 第4項
for g in range(G):
  for i in range(N):
    for j in range(N):
      _qubo[i+N*g,j+N*g] = _qubo[i+N*g,j+N*g] + _lambda2

### 第5項
for g in range(G):
  for i in range(N):
    _qubo[i+N*g,i+N*g] = _qubo[i+N*g,i+N*g] - 12*_lambda2

それでは、実際に計算を行う。

sampler_dw = EmbeddingComposite(dw_sampler)
solution = sampler_dw.sample_qubo(_qubo, num_reads=20)

続いて、計算結果を評価する。

for i in range(10):
  print(evaluation(solution,i))

<Output>
['0:', 0, 1, 18, 20]
['1:', 0, 1, 20, 18]
['2:', 0, 1, 22, 18]
['3:', 0, 1, 20, 20]
['4:', 0, 1, 18, 18]
['5:', 0, 1, 18, 16]
['6:', 0, 1, 18, 18]
['7:', 0, 1, 22, 18]
['8:', 0, 1, 18, 20]
['9:', 0, 1, 18, 20]

計算は完了したが、Pyquboを使用して計算した場合と比較して精度が極めて悪い。そこで、QUBO行列を小さくして、直接法とPyquboで比較してみる。

4.3 Pyqubo有無での比較検証

まず設定として、$N=3$、$G=2$、$\lambda1=3$、$\lambda2=5$、さらに$C$を以下の通り設定して進める。尚、定式化した数式はそのまま使用する。

_N, _G, _lambda1, _lambda2 = 3, 2, 3, 5
_C = [[0,2,1],[1,0,1],[0,2,0]]

まずは直接法の場合。

### Qubo行列の空箱の作成
_qubo = np.zeros((_N*_G)**2).reshape(_N*_G,_N*_G)

### 第1項
for g in range(_G):
  for i in range(_N):
    for j in range(_N):
      _qubo[i+_N*g,j+_N*g] =_qubo[i+_N*g,j+_N*g] - _C[i][j]

### 第2項
for i in range(_N):
  for g in range(_G):
    for l in range(_G):
      _qubo[i+_N*g,i+_N*l] = _qubo[i+_N*g,i+_N*l] + _lambda1

### 第3項
for i in range(_N):
  for g in range(_G):
    _qubo[i+_N*g,i+_N*g] = _qubo[i+_N*g,i+_N*g] - 2*_lambda1

### 第4項
for g in range(_G):
  for i in range(_N):
    for j in range(_N):
      _qubo[i+_N*g,j+_N*g] = _qubo[i+_N*g,j+_N*g] + _lambda2

### 第5項
for g in range(_G):
  for i in range(_N):
    _qubo[i+_N*g,i+_N*g] = _qubo[i+_N*g,i+_N*g] - 12*_lambda2

続いて、Pyquboの場合。

###スピンの設定
_x = Array.create(name="x",shape=(_N,_G), vartype = "BINARY")

###第1項の実装
_E = 0
for g in range(_G):
  for i in range(_N):
    for j in range(_N):
      _E = _E - _C[i][j]*_x[i,g]*_x[j,g]

###第2項実装
_C1 = 0
for i in range(_N):
  _c = 0
  for g in range(_G):
    _c = _c + _x[i,g]
  _C1 = _C1 + _lambda1*(_c -1)**2

###第3項実装
_C2 = 0
for g in range(_G):
  _c = 0
  for i in range(_N):
    _c = _c + _x[i,g]
  _C2 = _C2 + _lambda2*(_c - 6)**2

###ハミルトニアンの実装
_H = _E + _C1 + _C2

###QUBO行列の抜出
model = _H.compile()
Qubo,offset = model.to_qubo()

###Qubo行列で得られる座標を文字から数値へ変更する関数の定義
def change(test):
  _a =test[2:].replace("]","")
  _i = int(_a[:_a.find("[")])
  _g = int(_a[_a.find("[")+1:])
  _d = _i+_g*_N
  return _d

###座標データと数値の取り出し
_dim_list =[]
for _xy in Qubo:
  __dim_list = [change(_xy[0]),change(_xy[1])]
  _dim_list.append(__dim_list)

_value_list = []
for _value in Qubo.values():
  _value_list.append(_value)

###QUBO行列で表示
p_qubo = np.zeros((_N*_G)**2).reshape(_N*_G,_N*_G)

for k in range(len(_value_list)):
  p_qubo[_dim_list[k][0],_dim_list[k][1]]=_value_list[k]

出揃ったので、比較してみると以下。

p_qubo

<Output>
array([[-58.,   0.,   0.,   0.,   0.,   0.],
       [  7., -58.,   0.,   0.,   0.,   0.],
       [  9.,   7., -58.,   0.,   0.,   0.],
       [  6.,   0.,   0., -58.,   0.,   0.],
       [  0.,   6.,   0.,   7., -58.,   0.],
       [  0.,   0.,   6.,   9.,   7., -58.]])

_qubo

<Output>
array([[-58.,   3.,   4.,   3.,   0.,   0.],
       [  4., -58.,   4.,   0.,   3.,   0.],
       [  5.,   3., -58.,   0.,   0.,   3.],
       [  3.,   0.,   0., -58.,   3.,   4.],
       [  0.,   3.,   0.,   4., -58.,   4.],
       [  0.,   0.,   3.,   5.,   3., -58.]])

直接的に作成したQubo行列は、上三角行列と下三角行列に数値が別れてしまっている。例えば、$[1,2]$の$3$と、$[2,1]$の$4$であるが、これは1さんと2さんの相性、もしくは2さんと1さんの相性を示すが、いずれも同じもののはず。ならば、使用するQubo行列を上下いずれかの三角行列に絞って、$[1,2]=0$、$[2,1]=3+4=7$とすべき。従い、直接的にQubo行列を作成する場合、最後に上下三角行列の統一が必要となる。

尚、この解釈に至るまで、Qiita:巡回セールスマン問題を量子アニーリングを用いて解くを参考にした。つまり、以下のような計算式をQubo計算時に追加すればよい。

### 上下三角行列の統一
for y in range(len(_qubo)):
  for x in range(y,len(_qubo)):
    if not x ==y:
      _qubo[x,y] = _qubo[x,y]+_qubo[y,x]
      _qubo[y,x] = 0
_qubo

<Output>
array([[-58.,   0.,   0.,   0.,   0.,   0.],
       [  7., -58.,   0.,   0.,   0.,   0.],
       [  9.,   7., -58.,   0.,   0.,   0.],
       [  6.,   0.,   0., -58.,   0.,   0.],
       [  0.,   6.,   0.,   7., -58.,   0.],
       [  0.,   0.,   6.,   9.,   7., -58.]])

実際の問題で、再度計算してみる。

4.3 直接法の再計算

### Qubo行列の空箱の作成
_qubo = np.zeros((N*G)**2).reshape(N*G,N*G)
### lambda1、lambda2の設定
lambda1 = best_lambda1
lambda2 = best_lambda2

### 第1項
for g in range(G):
  for i in range(N):
    for j in range(N):
      _qubo[i+N*g,j+N*g] =_qubo[i+N*g,j+N*g] - C[i][j]

### 第2項
for i in range(N):
  for g in range(G):
    for l in range(G):
      _qubo[i+N*g,i+N*l] = _qubo[i+N*g,i+N*l] + lambda1

### 第3項
for i in range(N):
  for g in range(G):
    _qubo[i+N*g,i+N*g] = _qubo[i+N*g,i+N*g] - 2*lambda1

### 第4項
for g in range(G):
  for i in range(N):
    for j in range(N):
      _qubo[i+N*g,j+N*g] = _qubo[i+N*g,j+N*g] + lambda2

### 第5項
for g in range(G):
  for i in range(N):
    _qubo[i+N*g,i+N*g] = _qubo[i+N*g,i+N*g] - 12*lambda2

### 上下三角行列の統一
for y in range(len(_qubo)):
  for x in range(y,len(_qubo)):
    if not x ==y:
      _qubo[y,x] = _qubo[y,x]+_qubo[x,y]
      _qubo[x,y] = 0

### Dwaveでの計算
sampler_dw = EmbeddingComposite(dw_sampler)
solution = sampler_dw.sample_qubo(_qubo, num_reads=20)
 
### 計算結果の評価
for i in range(10):
  print(evaluation(solution,i))

<Output>
['0:', 0, 1, 22, 20]
['1:', 0, 1, 16, 22]
['2:', 0, 1, 18, 20]
['3:', 0, 1, 16, 20]
['4:', 0, 1, 20, 22]
['5:', 0, 1, 24, 18]
['6:', 0, 1, 22, 16]
['7:', 0, 1, 20, 26]
['8:', 0, 1, 16, 24]
['9:', 0, 1, 20, 18]

Qubo行列は同じ内容のはずだが、結果が伴わない。その後、色々と試行錯誤を行ったが、判明せず。直接的な表現に関しては、今後の課題としたい。ただし、いずれにせよ、Pyqubo用いて2次元スピンを容易に解決できることを確認。とりあえずは、Pyquboを用いて進めていく。

おわりに

2次元スピンを使用する問題として、代表的なものが巡回セールスマン問題がある。 次回はこちらを解いてみるとする。

今後の課題 : 巡回セールスマン問題
積み残し: 補助スピン$y_{i}$の活用, TTSを用いた探索