Qiskitで関数の勾配(微分)を計算する

量子コンピューターで関数の勾配、つまり微分を計算はよく登場します。
特に機械学習の分野ではモデルの学習時に関数の勾配情報を用いる手法が多数存在するため、量子機械学習の分野が特に注目を集めている今、勾配計算の重要性や高くなっているはずです。

量子コンピューターで関数の勾配情報を求める手法はいくつかありますが、その中でもよく目にするパラメータシフト法についての概要をQiskitを使った実装と併せて解説していきます。

スポンサーリンク

パラメータシフト法とは

パラメータシフト法とはパラメータ$\theta$の目的関数$\langle H(\theta)\rangle$の勾配を以下の計算式で求めることが出来るという計算手法です。

$$\frac{d\langle H(\theta)\rangle}{d\theta}=\frac{1}{2}(\langle H(\theta+\frac{\pi}{2})\rangle-\langle H(\theta-\frac{\pi}{2})\rangle)$$

目的関数はVQEのアルゴリズムにおいてはエネルギー期待値とみなせます。

この手法は上の式を見ると単純な差分法によく似ていますが、パラメータを大きくシフトさせて差分を取っているため差分法よりも誤差に堅牢であるという特徴があります。

$$\frac{d\langle H(\theta)\rangle}{d\theta}\sim\frac{1}{2h}(\langle H(\theta+h)\rangle-\langle H(\theta-h)\rangle)$$

勾配の導出

では上の式がどのように導出されるかを確認していきましょう。例としてVQEにおいてエネルギー期待値のパラメータ微分を求めることを考えます。

仮定として試行状態を作る量子回路を1量子ビットの回転ゲートの形で表されるとします。

この時エネルギー期待値は次のように表されます。

$$\langle H(\theta)\rangle = \langle\psi|H|\psi\rangle = \langle0|e^{i\theta X/2}He^{-i\theta X/2}|0\rangle$$

よってこれを$\theta$で微分すると下のように書けます。

$$\frac{d\langle H(\theta)\rangle}{d\theta}=\frac{i}{2}(\langle0|Xe^{i\theta X/2}He^{-i\theta X/2}|0\rangle-\langle0|e^{i\theta X/2}He^{-i\theta X/2}X|0\rangle)$$

この物理量はそのまま観測することが難しく、補助ビットを用いることでこの微分表式を直接評価できる手法が提示されています。しかし使用できる量子ビットに限りあるNISQデバイスにおいては補助ビットの数を極力減らしたいものです。

そこで、この式をパウリ演算子の性質を用いて変形してみます。

$$e^{-i\theta X/2} = \cos \frac{\theta}{2}I – i\sin \frac{\theta}{2}X$$

この式を利用することで最初に載せた式を導くことが出来ます。

$$\frac{d\langle H(\theta)\rangle}{d\theta}=\frac{1}{2}(\langle H(\theta+\frac{\pi}{2})\rangle-\langle H(\theta-\frac{\pi}{2})\rangle)$$

Qiskitで実装する

ではこの手法をQiskitで実装していきます。例題として2021年に開催されたQHACKという量子機械学習ハッカソンで出題された問題を利用してみます。

問題はこちらから見ることが出来ます。概要を説明すると、与えられたパラメータ化された量子回路、ハミルトニアンにおいて、ハミルトニアンの期待値の勾配を計算せよというものです。

与えられた量子回路は3量子ビットの回路で下図です。ハミルトニアンは$H = YIZ$です。

エネルギー期待値の勾配を計算するにあたって、今回は2つのアプローチで実装をしていきます。

  • Qiskitモジュールを利用しないで一から自分で実装する
  • QiskitのGradient Frameworkを利用して実装する

結果の正当性を確かめるために、配布されている正解データを使って検証していきます。入力パラメータは6つあるのでリストで渡します。

[1, 0.5, -0.765, 0.1, 0, -0.654]

出力として各パラメータにおける勾配を求めればよいので、出力も入力と同次元の一次元リストで返します。期待される出力は以下の値です。

[0, 0, 0, 0, -0.4553474723, 0]

1から自分で実装する

まずは必要なモジュールをインポートします。(後半で使用するモジュールも含んでいます。)

from qiskit import Aer
from qiskit.circuit import QuantumCircuit, ParameterVector
from qiskit.aqua import QuantumInstance
from qiskit.aqua.operators import I,X, Y, Z, StateFn, CircuitStateFn
from qiskit.aqua.operators.expectations import PauliExpectation, AerPauliExpectation
from qiskit.aqua.operators.converters import CircuitSampler
from qiskit.aqua.operators.gradients import Gradient
import numpy as np
import matplotlib as plt
%matplotlib inline

ハミルトニアンの期待値を求める関数を定義します。

# インスタンスを定義
backend = Aer.get_backend('qasm_simulator') 
q_instance = QuantumInstance(backend, shots=1024)

# ハミルトニアンの期待値
def exp_val(params):
    n = 3
    qc = QuantumCircuit(n)
    param_list = ParameterVector('Parameter', 2*n)

    for i in range(len(param_list)//n):
        qc.rx(param_list[3*i], 0)
        qc.ry(param_list[3*i+1], 1)
        qc.rz(param_list[3*i+2], 2)

        qc.cnot(0, 1)
        qc.cnot(1, 2)
        qc.cnot(2, 0)

    param_dict = dict(zip(param_list.params, params))
    qc.assign_parameters(param_dict, inplace=True)

    op = Z ^ I ^ Y # ハミルトニアンを定義
    psi = CircuitStateFn(qc) # 状態ベクトルを定義
    measurable_expression = StateFn(op, is_measurement=True).compose(psi) 

    # 期待値を計算
    expectation = AerPauliExpectation().convert(measurable_expression)

    sampler = CircuitSampler(q_instance).convert(expectation) 
    return sampler.eval().real

続いて勾配を計算します。パラメーターを$\pm \pi/2\だけシフトさせたハミルトニアンの期待値を計算してその差分を計算します。これを各パラメータにおいて計算します。

gradient = np.zeros_like(params) # 結果格納用のリスト
for i in range(len(params)):
    shifted = params.copy()
    shifted[i] += np.pi/2
    forward = exp_val(shifted)

    shifted[i] -= np.pi
    backward = exp_val(shifted)

    gradient[i] = 0.5 * (forward - backward)

結果を確認してみます。期待された結果を出力していて正しい結果を求められていることが分かりました。

print(np.round(gradient, 10))
# > [ 0.          0.          0.          0.         -0.45534747  0.        ]

QiskitのGradient Frameworkを利用して実装する

続いてはQiskitが提供しているGradient Frameworkを利用して先ほどよりも簡単に勾配を計算していきます。

公式Document: https://qiskit.org/documentation/tutorials/operators/02_gradients_framework.html

実装は簡単で、Gradientモジュールにハミルトニアンと量子状態を掛け合わせた演算子を渡すだけです。

def parameter_shift_gradient(params):
    n = 3
    qc = QuantumCircuit(n)
    param_list = ParameterVector('param_list', 2*n)

    for i in range(len(param_list)//n):
        qc.rx(param_list[3*i], 0)
        qc.ry(param_list[3*i+1], 1)
        qc.rz(param_list[3*i+2], 2)

        qc.cnot(0, 1)
        qc.cnot(1, 2)
        qc.cnot(2, 0)

    param_dict = dict(zip(param_list.params, params))

    obs = Z ^ I ^ Y
    psi = CircuitStateFn(primitive=qc, coeff=1.)
    operator = ~StateFn(obs) @ psi

    state_grad = Gradient(grad_method="param_shift").convert(operator=operator, params=param_list)
    result = state_grad.assign_parameters(param_dict).eval()
    return result

勾配を計算すると上と同じく期待された結果を得ることが出来ました。

print(np.round(np.array(parameter_shift_gradient(params)).real, 10))
# > [ 0.          0.          0.          0.         -0.45534747  0.        ]

まとめ

パラメータシフト法を用いることで量子回路上で関数の勾配を計算することが出来ました。この手法は期待値を評価する量子回路をそのまま用いて評価が可能という点で特徴的です。

さらにこの手法はより一般の場合についても拡張が可能なことが知られています。興味のある方はこちらもご覧ください。

G.E.Crooks. Gradients of parameterized quantum gates using the parameter-shift rule and gate
decomposition

今回掲載したソースコードは私のGitHubでも公開しているので全体をご覧になりたい方はこちらにアクセスしてください。

qiskit-implementation/parameter-shift.ipynb at main · daiki0623/qiskit-implementation
Contribute to daiki0623/qiskit-implementation development by creating an account on GitHub.

コメント

タイトルとURLをコピーしました