こちらの記事で紹介した Initializer と量子振幅推定(QAE)を使って、
積分計算を行ってみました。
計算式としては、Qiskit tutorial を少し平易にしたものになっています。
最後にこの計算がモンテカルロ法の高速化に寄与することについても触れたいと思います。
具体的には次の定積分の計算をします。
$$
\int_{-1}^{3} f(x)g(x) dx \
f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x - \mu)^2}{2 \sigma^2}} , g(x) = { \frac{(7x + 9)\pi}{256} }^2
$$
$\sigma^2 = 0.5, \mu = 1.0$とします。( $ g(x) $ に含まれる定数は後で値を合わせるために恣意的に設定しています)
Qiskitのライブラリで正規分布を表現する
Qiskit の NormalDistribution を使って、$f(x)$ を表現します。量子ビット数は3にします。
import qiskit
from qiskit_finance.circuit.library import NormalDistribution
from qiskit import QuantumCircuit, execute, Aer, IBMQ
from qiskit.visualization import plot_histogram
import numpy as np
# f(x)のパラメータを設定
num_qubits = 3
high = 3
low = -1
mu = 1.0
sigma = 0.5
normal = NormalDistribution(num_qubits, mu=mu, sigma=sigma, bounds=(low, high))
$f(x)$ が量子ビットで表現されたことを確かめるために 10,000 回観測して、確率振幅を推定します。
qc = QuantumCircuit(num_qubits, num_qubits)
qc.append(normal, list(range(num_qubits)))
for i in range(num_qubits):
qc.measure(i, i)
backend = Aer.get_backend('qasm_simulator')
results = execute(qc, backend=backend, shots=10000).result()
answer = results.get_counts()
print("測定結果：", answer)
# 確率振幅を求める
key = sorted(answer)
value = np.array([answer[k] for k in key], dtype='float64')
value /= sum(value)
print("確率振幅：", value)
測定結果： {'101': 1566, '001': 380, '111': 61, '000': 63, '011': 2915, '010': 1597, '110': 413, '100': 3005}
確率振幅： [0.0063 0.038 0.1597 0.2915 0.3005 0.1566 0.0413 0.0061]
次に、$f(x)$ を実際に計算した値を求めて、両者をグラフで表記します。すると、ほぼ一致していることが分かります。
%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
%matplotlib inline
# f(x)を直接計算する
fx = [np.exp( -(low + k*(high-low)/(2**num_qubits - 1) - mu)**2 / 2 / sigma ) / np.sqrt(2*np.pi*sigma) for k in range(2**num_qubits)]
key = sorted(answer)
plt.plot(key, value)
plt.plot(fx / sum(fx))
plt.xlabel('$x$')
plt.ylabel('$p(x)$')
f(x)g(x) を表現する
$g(x)$は2次関数となっていますが、
量子回路による実装を簡単にするために、$x$ が小さいときに $x^2 \sim \sin{x}$ と近似できることを利用します。
この近似によって、$g(x)$との積を回転ゲートで表現できるようになります。
この演算を行う量子回路を $U_A$ として次のように実装します。同時に逆演算を行う $U_A^\dagger$ も定義しておきます。
b = np.pi / 8
alpha = 0.5
def U_A(qc, n_encode, b, alpha):
qc.append(normal, range(1, n_encode+1))
qc.ry(alpha * b / 2**(n_encode - 1), 0)
for i in range(n_encode):
qc.cry(b / 2**(n_encode - i - 1), i+1, 0)
def U_Adg(qc, n_encode, b, alpha):
for i in reversed(range(n_encode)):
qc.cry(-b / 2**(n_encode - i - 1), i+1, 0)
qc.ry(-alpha * b / 2**(n_encode - 1), 0)
qc.append(normal.inverse(), range(1, n_encode+1))
qc = QuantumCircuit(num_qubits + 1, 1)
U_A(qc, num_qubits, b, alpha)
qc.measure(0, 0)
qc.draw('mpl')
では、この量子回路が正しい実装になっているかどうかを検証します。
例えば、
$q_1 = q_2 = q_3 = 0$ のとき、$q_0$ は $\pi/64$だけ回転されます。
このとき、$ x = -1 $ になるので、$g(-1) = (\pi/128)^2 \sim \sin^2{\pi/128}$ になります。
従って、$q_0$を回転させる計算と、$g(-1)$との積を求める計算が一致しています。
$q_1, q_2, q_3$ がとりうる値に対して、$q_0$の回転角と$g(x)$の値を表にまとめると以下のようになり、実装がうまくできていることが分かります。
| $q_1$ | $q_2$ | $q_3$ | $x$ | $q_0$ の回転角 | $g(x)$ |
|---|---|---|---|---|---|
| 0 | 0 | 0 | -1 | $\pi/64$ | $(\pi/128)^2$ |
| 0 | 0 | 1 | -3/7 | $3\pi/64$ | $(3\pi/128)^2$ |
| 0 | 1 | 0 | 1/7 | $5\pi/64$ | $(5\pi/128)^2$ |
| 0 | 1 | 1 | 5/7 | $7\pi/64$ | $(7\pi/128)^2$ |
| 1 | 0 | 0 | 9/7 | $9\pi/64$ | $(9\pi/128)^2$ |
| 1 | 0 | 1 | 13/7 | $11\pi/64$ | $(11\pi/128)^2$ |
| 1 | 1 | 0 | 17/7 | $13\pi/64$ | $(13\pi/128)^2$ |
| 1 | 1 | 1 | 3 | $15\pi/64$ | $(15\pi/128)^2$ |
量子振幅増幅を行う
この記事では解説を省きますが、量子振幅推定を行うために振幅増幅オラクル $Q$ を実装する必要があります。
その準備として、上で実装した $U_A$ に加えて、$ U_{S_0}, U_{S_x} $ を実装します。
$U_{S_0}$ は$q_0$以外の全ビットが0のときに符号反転を行う演算子で、$U_{S_x}$ は$q_0$が1のときに符号反転を行う演算子になります。
以上を組み合わせて $Q$ を実装することができます。
$$
Q = U_A U_{S_0} U_A^\dagger U_{S_x}
$$
def U_S0(qc, n_encode, b, alpha):
for i in range(1 + n_encode):
qc.x(i)
qc.h(0)
qc.mct(list(range(1, n_encode+1)), 0)
qc.h(0)
for i in range(1 + n_encode):
qc.x(i)
def U_Sx(qc):
qc.z(0)
# 振幅増幅オラクル
def Q(qc, n_encode):
U_Sx(qc)
U_Adg(qc, n_encode, b, alpha)
U_S0(qc, n_encode, b, alpha)
U_A(qc, n_encode, b, alpha)
振幅増幅オラクルの確認
振幅増幅オラクルが想定通りに実装されていることを確認するために、
振幅増幅を行った後、$q_0$の状態を測定します。
振幅増幅を2回行って$q_0$の状態を測定すると、1が測定される確率が高くなっており、振幅増幅されていることが分かります。
qc = QuantumCircuit(num_qubits + 1, 1)
U_A(qc, num_qubits, b, alpha)
Q(qc, num_qubits)
Q(qc, num_qubits)
qc.measure(0, 0)
results = execute(qc, backend=backend, shots=10000).result()
answer = results.get_counts()
print(answer)
{'1': 7210, '0': 2790}
MLQAE の実装
量子振幅推定の手法として、量子位相推定を必要としない振幅推定法(MLQAE)を採用します。
MLQAE の解説も詳しくは記しませんが、次の手順で振幅推定を行う手法です。
- 振幅増幅を0,1,2,4,8…回行った後(これをExponentially Incremental Sequenceという)にそれぞれ $q_0$ を測定し、状態1が測定される確率を記録
- 最尤推定法によって、記録された確率を満たす可能性が最も高い確率振幅を算出
尚、次の実装では確率振幅の算出に optuna を用いた探索を行っています。
import optuna
# Exponentially incremental sequence
def EIS(n_encode, n_iteration, shots):
m = [0] + [2**i for i in range(n_iteration)]
h = []
for i in range(n_iteration + 1):
qc = QuantumCircuit(1 + n_encode, 1)
U_A(qc, n_encode, b, alpha)
for j in range(m[i]):
Q(qc, n_encode)
qc.measure(0, 0)
results = execute(qc, backend=backend, shots=shots).result()
res = results.get_counts()
for key in res.keys():
if key[0] == '1':
h.append(res[key])
if len(res.keys()) == 1: # "1"が一度も観測されなかったときリストに0を追加する
h.append(0)
return m, h
# 尤度関数の最大化
def logML(m, h, shots):
def ML_loss(trial):
L = []
eps = 1e-30
for i in range(len(m)):
theta = trial.suggest_uniform('theta', 0, 2 * np.pi)
A = h[i] * np.log(np.sin((2. * m[i] + 1) * theta)**2 + eps)
B = (shots - h[i]) * np.log(np.cos((2. * m[i] + 1) * theta)**2 + eps)
Lk = A + B
L.append(Lk)
return -np.sum(L) #-np.prod(L)
study = optuna.create_study()
study.optimize(ML_loss, n_trials=50)
res = study.best_params['theta']
ans = (np.sin(res))**2
return ans
MLQAEによる振幅推定の実行
ようやく実装が出来たので、MLQAEで振幅推定を行います。
複数回実施して求められた確率振幅の平均と分散を表示させます。
n_iteration = 2
shots = 100
ans_list = []
for i in range(100):
m, h = EIS(num_qubits, n_iteration, shots)
ans = logML(m, h, shots)
ans_list.append(ans)
ave = np.mean(np.array(ans_list))
var = np.std(np.array(ans_list))**2
print("Mean:", ave)
print("Variance:", var)
Mean: 0.041630808767425066
Variance: 4.585873114994599e-05
振幅推定結果の評価
上の表に記載された $ g(x) $ と $f(x)$ を掛けたものを足し合わせることで簡易的に
$ \int_{-1}^{3} f(x)g(x) dx $ を計算します。
- MLQAEの結果：0.04168
- 簡易計算の結果：0.04207
- $\sin^2{x} \sim x^2$ の近似を行わない計算結果：0.04142
これより、MLQAEによって定積分の値が求まっていると言えます。
（わずか3量子ビットを使った計算のため厳密な評価は行いません。）
# 近似あり
amp = sum([((i*2 + 1)/128*np.pi)**2 *j for i, j in zip(range(8), value)])
print("近似あり：", amp)
# 近似なし
amp = sum([(np.sin((i*2 + 1)/128*np.pi))**2 *j for i, j in zip(range(8), fx / sum(fx))])
print("近似なし：", amp)
近似あり： 0.04232749481232815
近似なし： 0.04141884447243714
モンテカルロ法への適用可能性について
最後に、QAE がモンテカルロ・シミュレーションの計算を高速化するのに有用だと期待されている理由を見てみます。
モンテカルロ・シミュレーションの簡単な実装例として、
Wikipedia に記載されているオルンシュタイン・ウーレンベック過程の数式と、
それを解くために使われるオイラー・丸山法の実装を確認します。
$$
d Y_t = \theta \cdot (\mu - Y_t) dt + \sigma d W_t \
Y_0 = Y_{init}
$$
ここで、乱数を使う $ \sigma d W_t $ の部分に着目すると、
オイラー・丸山法では正規乱数から値を取得して多数の時系列グラフを作成し、それらの平均を求めています。
def dW(delta_t):
"""Sample a random number at each call."""
return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))
for _ in range(num_sims):
for i in range(1, ts.size):
t = t_init + (i - 1) * dt
y = ys[i - 1]
ys[i] = y + mu(y, t) * dt + sigma(y, t) * dW(dt)
plt.plot(ts, ys)
この部分の計算は（全く同じではありませんが）上記のQAEを使った積分計算に類似していることから、
将来的に QAE がモンテカルロ・シミュレーションの高速化に発展していくことが推察されるかと思います。