IBM Qでアダマールテストをやってみよう

1. はじめに

アダマールテストは有名な量子アルゴリズムの一つです。

アダマールテストとそのマイナーチェンジによって、ユニタリー演算子の任意の状態での期待値を推定することができます。さらには、ユニタリー演算子の固有状態での期待値を求めることで、ユニタリー演算子の固有値の位相推定が可能となります。
この事実自体が量子計算を行う上で非常に強力です。加えて、アダマールテストの発展である位相推定アルゴリズムは素因数分解のアルゴリズム（Shorのアルゴリズム）に利用されているため、アダマールテストは応用でも力を発揮しています。

アダマールテストの仕組みについては、量子アルゴリズムの基本：アダマールテストの確認 が非常に詳しいです。本記事でも多大に参考にさせていただいています。
どのように期待値計算をするかはそちらを参照していただくことにして、本記事ではアダマールテストの紹介自体は簡単にとどめ、テストの実装と実行について主に書き連ねていくことにします。

2. アダマールテストの簡単な紹介

アダマールテストは、次の量子回路で表現され、

|0> ---H---*---H---M---
|phi> -----U-----------

状態$| \psi \rangle$に対する$U$の期待値の実部が

{Re \langle \psi |U| \psi \rangle) = p_{0} - p_{1}
}

と求められます。

また、アダマールテストのマイナーチェンジは、下の量子回路で与えられます。

|0> ---H---S---*---H---M---
|phi> ---------U-----------

制御ユニタリーゲートの手前で、補助量子ビットに位相ゲート$S$が付け加えられています。

状態$| \psi \rangle$に対する$U$の期待値の虚部が、

{Im \langle \psi | U | \psi \rangle = p_{1}^{\prime} - p_{0}^{\prime}
}

で計算できます。

実際に量子コンピュータを用いてアダマールテストを行う際には、多数の測定をして（サンプリングを行い）、確率分布{$ p_1,p_0 $}を推定することになります。サンプリング回数を増やすことにより、サンプリングで得られる確率分布が理想の確率分布に近づいていき、推定の精度を上げることができます。

3. アダマールテストを実装・実行する

それでは、アダマールテストを行ってみましょう。

IBM Q Experienceを用いて、Qiskitで実装したアダマールテストを実行してみます。IBM Q Experienceはクラウド上の量子コンピュータで、Qiskitは量子コンピュータソフトウェア開発のフレームワークです。
初めてIBM Q Experienceを利用する方は、初めにアカウントの登録とAPIの取得が必要です。3-2. QiskitとIBM Q Experienceの使い方が詳しいので、参考にしていただきたいです。

3.1. Qiskitのインストールおよび利用するモジュールのインポート

#google colablatory上では!が必要
!pip install qiskit

#必要なモジュールのインポート
from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute, Aer
from qiskit.qasm import pi
from qiskit.tools.visualization import plot_histogram, circuit_drawer
import numpy as np

## "MY_API_TOKEN" に、上で取得したtokenを入れる
IBMQ.save_account('MY_API_TOKEN')
# 自分のアカウント情報をloadする。（あらかじめ IBMQ.save_account を実行しておく必要がある. 複数のアカウントを使い分ける時はここで行う)
provider = IBMQ.load_account()

これで、IBM QとQiskitが利用できるようになりました。
それでは、アダマールテストを作製していきましょう。

3.2. アダマールテスト（量子回路）の作製

def hadamard_test_real_circuit():

  q = QuantumRegister(3)
  c = ClassicalRegister(3)
  qc = QuantumCircuit(q, c)

  qc.h(q[0])
  qc.ch(q[0],q[1])
  qc.cu3(pi/8,-pi/2,pi/2,q[0], q[2])
  qc.h(q[0])

  qc.measure(q[0], c[0])

  return qc

def hadamard_test_img_circuit():

  qc = QuantumCircuit(3, 3)

  qc.h(0)
  qc.s(0)
  qc.ch(0,1)
  qc.cu3(pi/8,-pi/2,pi/2, 0,2)
  qc.h(0)

  qc.measure(0, 0)

  return qc

まず、期待値の実部と虚部を求めるアダマールテストをそれぞれ定義しました。
期待値の実部を求めるアダマールテストではregisterでbitを定義してみました。registerを定義する方が、複雑なコードになった場合に書きやすくなると思います。

また、ここで定義した量子状態（q-bit）ははじめに初期化されるので、初期化された状態{ $|0>,|0>$}でのユニタリー演算子の期待値をこの量子回路では求めることになります。

そして、テストに用いるユニタリー演算子$U$には、$U = H_1 Rx_{2}(1/8 \pi)$を考えることにします（なんでもよいです）。下付き添え字は、q-bitを指定しています。$0$は制御bitに割り当てられるので、1以降の数を用いています。

補足
コードを読んで、Rxが見当たらないと不思議に思った方もいらっしゃると思います。実は、$U3(\theta, -\pi/2, \pi/2)$が$Rx(\theta)$と等価であるため、$CRx$を$CU3$に差し替えた表現になっています。これは、$CRx$が標準ゲートとして存在しないという事情によります。$CRy$も同様で、$CU3$を用いて表現しなければいけません。

3.3. アダマールテストを実行するための関数を定義する。

def excute_circuit(qc, local=None, simulator=False, shots=4096):
  from qiskit.providers.ibmq import least_busy

  if local is None:
    #least busyだったbackendを選ぶ simulator = TrueならCloud SM
    backend = least_busy(provider.backends(simulator=simulator, operational=True))
    if simulator:
      print("Simulator: ", backend)
    else:
      print("Least busy backend: ", backend)

  elif local == 'local':
    #localのシミュレータをbackendに
    backend = Aer.get_backend("qasm_simulator")
    print("Local simulator: ", backend)

  #量子回路qcを指定したバックエンド(backend_sim)でshots回実行する。
  result = execute(qc, backend, shots=shots).result()

  return result

お次は、量子回路を実行するための関数です。バックエンドに量子コンピュータの実機もしくはシミュレータを利用できるようになっています。実機IBM Qは何台かあるのですが、least_busyで一番すいているものを取り出しています。また、シミュレータにはLocalなものとCloud上のものが選択できます。simulatorとlocalがバックエンド選択のための引数になっています。

3.4. アダマールテストを実行してみた。

def hadamard_test_real(shots=4096, local=None,simulator=False):
  hrc = hadamard_test_real_circuit()
  result_real = excute_circuit(hrc,local,simulator,shots)

  # #結果を出力する。
  # print(result_real.get_counts())
  c_str = str(result_real.get_counts()).strip('{''}').split(',')
  c_0 = int(c_str[1].split(':')[1])
  c_1 = int(c_str[0].split(':')[1])

  pr_0 = c_0/shots
  pr_1 = c_1/shots
  exp_r = pr_0 - pr_1
  print('estimation of real part of expectation value', exp_r)

  return exp_r

def hadamard_test_img(shots=4096, local=None, simulator=False):
  hic = hadamard_test_img_circuit()
  result_img = excute_circuit(hic, local, simulator, shots)

  # #結果を出力する。
  # print(result_img.get_counts())
  c_str = str(result_img.get_counts()).strip('{''}').split(',')
  c_0 = int(c_str[1].split(':')[1])
  c_1 = int(c_str[0].split(':')[1])

  pi_0 = c_0/shots
  pi_1 = c_1/shots

  exp_i = pi_1 - pi_0
  print('estimation of imaginary part of expectation value', exp_i)

  return exp_i

　if __name__ == '__main__':
  exp_r = hadamard_test_real()
  exp_i = hadamard_test_img()
  print('\n')

  print('estimation of expectation value', exp_r+exp_i*1j)

上二つの関数を組み合わせた期待値の推定値を出力するプログラムです。shotsは測定回数に相当します。デフォルトは、IBM Q実機で4096回測定をするようになっています。
ここまでを実行してみましょう。

Least busy backend:  ibmq_vigo
estimation of real part of expectation value 0.7392578125
Least busy backend:  ibmq_ourense
estimation of imaginary part of expectation value -0.068359375

estimation of expectation value (0.7392578125-0.068359375j)

このような表示が出てくればOKです。実機なので実行スピードは遅く結果の表示までが長いのですが、その分本当に量子コンピュータを使ってるんだなという実感が僕には湧きました。げに熱いです。

3.5. Qulacsを用いて期待値を計算してみる（not equal 推定）。

アダマールテストはあくまで期待値推定を行うアルゴリズムですので、実際の期待値からはズレています。実際の期待値を計算することで、そのズレを見積もってみましょう。

今考えていたユニタリー演算子は、適当に基底(z-basisなど）を考え4×4行列として表現できるでので手計算で求めてもいいのですが、ちょっと面倒です。
そこで、量子シミュレータQulacsでユニタリー演算子（量子ゲート）$U$を用意して、初期化された状態に対する$U$の期待値を計算してみます。

## Google Colaboratoryの場合 ・ Qulacsがインストールされていないlocal環境の場合のみ実行してください
!pip install qulacs

まず、Qulacsをインストールして

import numpy as np
from qulacs import QuantumState
from qulacs.gate import X, RX, H, merge
from qulacs.state import inner_product
n = 2
state = QuantumState(n)
state_bra = QuantumState(n)
state_ket = QuantumState(n)
# state_bra.set_zero_state()
# state_ket.set_zero_state()

#H
h_gate = H(0)
#RX
angle = np.pi / 8
rx_gate = RX(1, angle)

# ゲートを合成して新たなゲートを生成
# 第一引数が先に作用する
h_and_rx_gate = merge(h_gate, rx_gate)

# ゲート行列を取得
matrix = h_and_rx_gate.get_matrix()
print(matrix)

#状態を変化させる
h_and_rx_gate.update_quantum_state(state_ket)
print(state_ket.get_vector())

# 内積値の計算
value = inner_product(state_bra, state_ket)
print('calculation of expectation value', value)

以上のように期待値を計算します。
結果は、

estimation of expectation value (0.7392578125-0.068359375j)
calculation of expectation value (0.6935199226610737+0j)

と出るはずです（下段参照）。
先のアダマールテストも上段に表示しています。

実際の期待値から実部が0.05 程度、虚部が0.07 程度、アダマールテストで得られた結果はズレていることがわかります。このズレが妥当なのか少し検討してみます。

3.6. アダマールテストのエラーの原因を探る。

アダマールテストでは、総測定回数$shots$に対して、$0$ or $1$がそれぞれいくら出たかで、確率分布を推定します。それゆえ、その精度は$1/shots$で決まることになります。$shots = 4096$だとすると、$1/4096$ ~ $0.0002$程度の精度で推定できることになります。実際に得られたズレよりかなり小さいです。この原因の一つは、実機を使うことによるエラー（うまい表現ができませんが）にあると考えられます。

そこで、アダマールテストを量子シミュレータで計算することで、実際に、実機によるエラーが原因かどうか確かめてみることにしましょう。
'local'=localとして、期待値推定のプログラムを実行してみます。

Local simulator:  qasm_simulator
estimation of real part of expectation value -0.6796875
Local simulator:  qasm_simulator
estimation of imaginary part of expectation value -0.0322265625

estimation of expectation value (-0.6796875-0.0322265625j)
calculation of expectation value (0.6935199226610737+0j)

初期化された状態に対する$U$の期待値の量子シミュレータによる推定結果です。最下段に期待値を計算した結果も表示しています。実部が0.01 程度、虚部が0.03 程度、シミュレーションで得られた結果はズレていることがわかります。実機よりズレは小さくなりましたが、予測していた精度$1/4096$ ~ $0.0002$よりも100倍程度大きいです。ぐぬぬ…。
アダマールテストの精度に関する考え方が間違っていたのでしょうか。それとも、まだ考慮できていないエラーの源があるのでしょうか。少し考えてみることにします（うまく着地できなかったです…）

4. まとめ

量子アルゴリズムの基礎としてよく取り上げられるアダマールテストを実機とシミュレータを用いて実行しました。このような基礎的な量子アルゴリズムの実装・実行は初めて実機を触る経験としてちょうどいいと思いますので、皆さんぜひIBM Qを使ってアダマールテストで遊んでみてください。
最後の節では、エラーの原因探求をしようと思いましたが、予想とはずれて着地点を見失ってしまいました。申し訳ありません…。アウトプットをする習慣を続けたいので、ひとまずまとまったところだけで記事を投稿したいと思います。後程、勉強して追記します。