今回は2022年のインターン向けの技術指導の5回目になります。これまでは古典機械学習、最適化、量子組合せ最適化、量子機械学習を見てきました。最後は量子最適化を使った機械学習になります。
今回実行するのはクラスタリングと線形回帰計算です。これらのモデルはイジングモデルを使って実装をすることができます。最適化問題では、問題を直接マッピングするほか、モデルを作って間接的に問題を解くこともできます。

クラスタリング

クラスタリングは条件によって指定のグループに分ける問題で、今回はノード間の距離によって指定の数のグループに分ける問題を考えます。今回はイジングモデルを利用しますが、このモデルでは、ノード数*クラスタ数だけ量子ビットを必要とします。

QUBO

定式化は下記になります。今回はノードを複数のクラスタにコピーを作り、クラスタに所属するノードは1、所属しないノードは0とすることによって、クラスタリングを設定できます。項は２つあります。こちらはクラスタ内の距離の総和が最小になるコスト関数。

H_A = \sum d_{ij} q_i q_j\ \ \ \ (q_i \in \{0, 1\})

また、制約条件は一つのノードはすべてのクラスタ内で一つだけ1になるという条件になります。

H_B = \sum(\sum_{i}^n q_i - 1)^2

ここで、量子もつれを利用すると、二番目の制約をなくすこともできます。その場合には２つ目のハミルトニアンとして下記のものを利用します。これは|01>と|10>を交換する式になっています。

H_{XY} = (X_0 X_1 + Y_0 Y_1)/2 = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from blueqat import Circuit
from blueqat.utils import qaoa
from blueqat.pauli import X,Y,Z
from blueqat.pauli import qubo_bit as q

import networkx as nx
import matplotlib.pyplot as plt

n = 10
m = 30
seed = 15

options = {'node_size': 100}

G = nx.gnm_random_graph(n, m, seed = seed)
nx.draw(G, **options)

<Figure size 432x288 with 1 Axes>

hamiltonian = sum(q(e[0])*q(e[1]) + q(e[0]+n)*q(e[1]+n)  for e in G.edges)
step = 1

#mixer and init state
mixer = 0.0
init = Circuit()
for i in range(n):
    mixer += 0.5*X[i]*X[i+n] + 0.5*Y[i]*Y[i+n]
    init.h[i].cx[i,i+n].x[i]

result = qaoa(hamiltonian, step, init, mixer)
b = result.circuit.run(shots=10)
sample = b.most_common(1)[0][0]
print("sample:"+ str(sample))

nx.draw(G, **options, node_color=[int(s) for s in list(sample[:n])])

線形回帰計算

次に線形回帰問題です。こちらは論文で構築方法が書いてあります。

Reference: https://arxiv.org/abs/2008.02355
まず目的変数のベクトルY、説明変数の行列Xと重みベクトルのwを考えます。説明変数Xと重みwから予測されるYとの誤差を最小するようにモデルを考えます。

min_wE(x) = ||Xw - Y||^2

これを展開すると、

min_wE(x) = w^TX^TXw -2w^TX^TY + Y^TY

ここで、 $Y^TY$ は定数なので、最適化問題の定式化から除外

min_wE(x) = w^TX^TXw -2w^TX^TY

こちらを解きます。ここで問題になるのが、QUBOだと離散しか扱えないので、 $w = P \hat w$ をつかってエンコードします。

min_wE(x) = \hat w^TP^TX^TX P \hat w -2 \hat w^TP^TX^TY

結局、 $A=P^TX^TXP$ として、 $b=P^TX^TY$ とすると、

min_w E(x) = \hat w^T A \hat w + \hat w^T b

を最小化すればよいことになります。

import numpy as np

K = 2 # bit number for weight
d = 3 # Number of features

p = [2 ** (-i-1) for i in range(K)]
I = np.eye(d)
P = np.kron(I, p)

w = np.array([0.25, 0.75, 0.5])
X = np.random.rand(100, 3)

y = X @ w + np.random.normal(scale = 0.05, size = X.shape[0])

A = np.dot(np.dot(P.T, X.T), np.dot(X, P))
b = -2 * np.dot(np.dot(P.T, X.T), y)
QUBO = A + np.diag(b)
QUBO = np.triu(QUBO + QUBO.T) - np.eye(QUBO.shape[0]) * QUBO

from blueqat.utils import qaoa
from blueqat.pauli import from_qubo

step = 2
hamiltonian = from_qubo(QUBO)

result = qaoa(hamiltonian, step)
b = result.circuit.run(shots=100)
sample = b.most_common(1)[0][0]
print(sample)

w_qa = P @ [int(i) for i in list(sample)]

print("Predicted weight:", w_qa)
print("True weight:", w)

Predicted weight: [0.25 0.75 0.75]
True weight: [0.25 0.75 0.5 ]

たまに正解が出ます。

インターン講座5：QAOAでクラスタリングと線形回帰計算

Yuichiro Minato

クラスタリング

QUBO

線形回帰計算