Qiskit Runtimeを触ってみる

量子コンピュータ応用の一つと考えられている量子化学計算では、量子計算と古典計算をハイブリッドしたアルゴリズムが提案されています。このようなアルゴリズムを実際に実施するには、「クラウド上の量子コンピュータ」と「ユーザのローカルの古典コンピュータ」の間を何度も往復して通信する必要があります。これは多くの実行時間を要し、効率的ではありません。この量子コンピューティングのワークロードのボトルネックを解消するため、IBMでは「IBM Quantumシステムと連結された古典コンピュータ上」で、ユーザが量子古典ハイブリッドのプログラムを束ねて実行できるQiskit Runtimeを提供しています。IBMの記事[※]を参照すると、Qiskit Runtimeで実行された量子化学計算のプログラムでは最大で120倍程度のスピードアップが見られた、とあります。

本記事ではQiskit Runtimeの初歩的な利用について説明をしたいと思います。
後半では期待値計算の応用として水素化ナトリウム分子(NaH)のエネルギー計算[※2]を実機で行ってみたいと思います。以下、本記事はqiskitのtutorial[※3]を参考に執筆しました。

[※] https://www.ibm.com/blogs/think/jp-ja/qiskit-runtime-for-useful-quantum-computing/
[※2]"Quantum Chemistry as a Benchmark for Near-Term Quantum Computers"
A. J. McCaskey, Z. P. Parks, J. Jakowski, S. V. Moore, T. Morris, T. S. Humble, R. C. Pooser,
https://arxiv.org/abs/1905.01534
[※3]https://qiskit.org/ecosystem/ibm-runtime/index.html

Qiskit Runtime IBM Clientのインストール

Qiskit Runtimeにアクセスするためのパッケージをインストールします。
「qiskit-ibmq-provider」の非推奨化に伴い、「qiskit-ibm-runtime」への移行[※4]が推奨されています。

pip install qiskit-ibm-runtime

本記事の執筆当時(2023/5)、最新versionは0.10.0です。

[※4] https://qiskit.org/ecosystem/ibm-runtime/locale/ja_JP/migrate/migrate-setup.html

Qiskit Runtime へのアクセスの確認

qiskit_ibm_runtimeにあるQiskitRuntimServiceクラスを用いて、Qiskit Runtimeにアクセスしてみます。
以下のようにQiskitRuntimServiceを初期化します。

from qiskit_ibm_runtime import QiskitRuntimeService, Session, Estimator, Sampler
from qiskit.test.reference_circuits import ReferenceCircuits

hub = "ibm-q"
group = "open"
project = "main"
service = QiskitRuntimeService(channel="ibm_quantum", 
                               instance=f"{hub}/{group}/{project}")

もし「APIトークンが見つかりません」といったエラー出れば、以下のようにsave_accout()メソッドでアカウントを保存します。一度保存すれば、次回以降この操作は必要ありません。

QiskitRuntimeService.save_account(channel="ibm_quantum", 
                                  token="MY_IBM_QUANTUM_TOKEN", overwrite=True)

Qiskit Runtimeにアクセス可能か確認するため、Qiskit Runtimeのtutorialページ[※5]を参考に、「Hello, World!」を出力する以下のPrototypeプログラム[※6]を走らせてみます。

program_inputs = {'iterations': 1}
runtime_inputs = {}
options = {"backend": "ibmq_qasm_simulator"}
job = service.run(program_id='hello-world',
                  options=options,
                  inputs=runtime_inputs,
                )
# print(f"job id: {job.job_id()}")
result = job.result()
print(result)

コードに間違いがなければ

Hello, World!

と出力されます。これでQiskit Runtimeを利用する準備が出来ました。

[※5] https://qiskit.org/ecosystem/ibm-runtime/getting_started.html
Tutorialページにある「Test your setup」では、options = {"backend_name": "ibmq_qasm_simulator"}とありますが、これではエラーとなります。ここではoptions = {"backend": "ibmq_qasm_simulator"}
と変更しています。

[※6]特定の研究用途に合わせたプログラム(vqe, qaoaなど)が既にいくつか用意されています。これらは後述するQiskit RuntimeのPrimitiveで構成されています。PrototypeプログラムはIBM Quantumにサインインし、「Qiskit Runtime」のメニューから確認することが出来ます。

Primitives

準備が出来ましたので、実際にQiskit Runtimeを利用した計算をしてみたいと思います。
その前に、Qiskit Runtimeのtutorialに沿うと、Primitives[※7]について説明がありますので、少し見てみたいと思います。説明を読むと、Primitiveとは「モジュール化されたアルゴリズムやアプリケーションを簡単に構築するためのコア関数」であり、「(エラー緩和された)擬確率を生成するSampler」と「期待値を計算するEstimator」が現時点で用意されていると書かれています。以前、「IBMの量子コンピュータを触ってみる」の記事では、計算結果は各量子ビットのPauli行列のZ成分($\sigma_z$)の測定値として得られていました。ですが、何か実践的なアルゴリズムやアプリを考えると、実際に欲しいのは量子回路の結果からサンプリングされた各状態の確率分布であったり、量子回路から得られる状態ベクトルを使った物理量（ハミルトニアンなど）の期待値$\langle \psi | \hat{O}| \psi \rangle$だったりするかと思います。状態ベクトルは$\sigma_z$の測定値から再構成できそうですが、用意されているPrimitiveを使えばそのような手間は不要になりそうです。

具体的にどういうことが出来るのか試してみます。

[※7]https://qiskit.org/ecosystem/ibm-runtime/locale/ja_JP/tutorials/how-to-getting-started-with-sampler.html

1. ローカル環境でのPrimitiveの利用: Samplerの場合

Tutorial[※7]に従うと、Primitivesはローカルな環境でも実行できます。まずPrimitivesのひとつのSamplerを試してみます。これは（エラー緩和された）擬確率を生成するPrimitiveとなっています。以下、5量子ビットのGHZ状態の量子回路をSamplerで実行したコード例を示します。これまで量子回路の計算にAerSimulator()を用いていた代わりに、Sampler()を用いて実行します。

from qiskit import QuantumCircuit
from qiskit.primitives import Sampler

sampler = Sampler() # samplerクラスのインスタンス作成

qc = QuantumCircuit(5,5)
qc.h(1)
qc.cx(1,2)
qc.cx(1,0)
qc.cx(2,3)
qc.cx(3,4)
qc.measure([0,1,2,3,4], [0,1,2,3,4])
#qc.draw('mpl')
job = sampler.run(qc)

print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")

result = job.result()
print(f">>> {result}")
print(f"  > Quasi-distribution: {result.quasi_dists[0]}")

実行すると、以下のような結果が得られました。

>>> Job ID: d23cc89c-efc0-45dd-aacf-ecce118939ec
>>> Job Status: JobStatus.RUNNING
>>> SamplerResult(quasi_dists=[{0: 0.4999999999999999, 31: 0.4999999999999999}], metadata=[{}])
  > Quasi-distribution: {0: 0.4999999999999999, 31: 0.4999999999999999}

'00000'と'11111'の状態が確率1/2で得られています。(Quasi-distributionのキー'0'と'31'は、2進数表記でそれぞれ'00000'と'11111'です。) これは状態ベクトルを用いた理想的な計算結果です。 以前のAerSimulator()ではshot数(測定数)に対する各状態の観測回数が返していましたが、Sampler()では擬確率がを返しています。

1-2. Qiskit RuntimeでのPrimitiveの利用: Samplerの場合

次にQiskit Runtimeを利用したSamplerを実行してみます。Runtimeを利用すると、量子コンピュータ実機でもPrimitiveを用いた計算が可能です。Primitvesを利用するには「qiskit.primitives」から「qiskit_ibm_runtime」のSamplerクラスをインポートします。また、QiskitRuntimServiceクラスを用いてbackendのリソースを指定します。以下ではまずテストとしてbackendに"ibmq_qasm_simulator"を指定してみます。

from qiskit import QuantumCircuit
from qiskit_ibm_runtime import Sampler

hub = "ibm-q"
group = "open"
project = "main"

service = QiskitRuntimeService(channel="ibm_quantum", 
                               instance=f"{hub}/{group}/{project}")
backend = service.backend("ibmq_qasm_simulator") # リソースを指定
sampler = Sampler(session=backend) # samplerクラスのインスタンス作成

qc = QuantumCircuit(5,5)
qc.h(1)
qc.cx(1,2)
qc.cx(1,0)
qc.cx(2,3)
qc.cx(3,4)
qc.measure([0,1,2,3,4], [0,1,2,3,4])
job = sampler.run(qc)

print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")

result = job.result()
print(f">>> {result}")
print(f"  > Quasi-distribution: {result.quasi_dists[0]}")

実行すると、以下のような結果が得られました。

>>> Job ID: cho3gv78rmtc64q9lk90
>>> Job Status: JobStatus.RUNNING
>>> SamplerResult(quasi_dists=[{31: 0.5145, 0: 0.4855}], metadata=[{'shots': 4000}])
  > Quasi-distribution: {31: 0.5145, 0: 0.4855}

次にbackendとして実機を指定してみます。利用可能なリソースは

service.backends()

で確認できます。

ここではbackendを"ibmq_lima"に変更して実行してみます。

backend = service.backend("ibmq_lima") 

実行すると、以下のような結果が得られました。

>>> Job ID: cho3odinajhpa643hbog
>>> Job Status: JobStatus.DONE
>>> SamplerResult(quasi_dists=[{0: 0.46569039220687825, 1: 0.008854774325396378, 2: 0.0028919657623803015, 3: 0.01167043376528802, 4: -0.0004984133076128505, 5: 0.0006935132529740888, 6: 0.0015530826195646443, 7: 0.005548204156321979, 8: 0.003944077891295497, 9: -0.0004316381730416175, 11: -0.0010122364931930874, 12: 0.0021789760382525444, 13: 0.0010441478075698147, 14: 0.0011117074039608719, 15: 0.00852907588769621, 16: 0.005076273076044268, 17: 0.00038706815552610606, 19: 8.628212759404266e-05, 20: -0.0005410317079811839, 21: 0.0005547686829574295, 22: -0.0007920126802441734, 23: -0.007783425412390393, 24: 0.007592845573009159, 25: 0.001037118563184794, 26: -0.00013721919819593693, 27: 0.00039797280488139, 28: 0.015434799262404554, 29: 0.003736589375049321, 30: 0.013414317750528686, 31: 0.44976759048390086}], metadata=[{'shots': 4000, 'readout_mitigation_overhead': 3.796154539914829, 'readout_mitigation_time': 0.1071654986590147}])
  > Quasi-distribution: {0: 0.46569039220687825, 1: 0.008854774325396378, 2: 0.0028919657623803015, 3: 0.01167043376528802, 4: -0.0004984133076128505, 5: 0.0006935132529740888, 6: 0.0015530826195646443, 7: 0.005548204156321979, 8: 0.003944077891295497, 9: -0.0004316381730416175, 11: -0.0010122364931930874, 12: 0.0021789760382525444, 13: 0.0010441478075698147, 14: 0.0011117074039608719, 15: 0.00852907588769621, 16: 0.005076273076044268, 17: 0.00038706815552610606, 19: 8.628212759404266e-05, 20: -0.0005410317079811839, 21: 0.0005547686829574295, 22: -0.0007920126802441734, 23: -0.007783425412390393, 24: 0.007592845573009159, 25: 0.001037118563184794, 26: -0.00013721919819593693, 27: 0.00039797280488139, 28: 0.015434799262404554, 29: 0.003736589375049321, 30: 0.013414317750528686, 31: 0.44976759048390086}

実機の場合、少し他の状態も得られています。（擬確率なので負値もありますが、全ての状態の確率和をとると1.0になっています。）

Tips: 結果の回収

Qiskit Runtimeを用いて実機にジョブを投入する場合、待ち行列のために計算終了まで数日要することがあります。そのような場合、以下のように計算終了後にJob IDから計算結果を回収することができます。

job_id = 'MY_JOB_ID'
retrieve_job = service.job(job_id)
result = retrieve_job.result()

計算jobが正しくsubmitされていれば、仮に結果が回収できず途中でjupyter notebookをシャットダウンしたとしても、後からIBM Quantumにサインインして「Jobs」情報からJob IDを確認できます。

Session

Qiskit RuntimeでPrimitive(Sampler)を用いた計算が出来ました。さて、Qiskit Runtimeの利点を思い出しますと、冒頭にも述べたように「ローカルの古典コンピュータ」ではなく「IBM Quantum systemと連結した古典コンピュータ」で量子-古典間の閉じた計算が可能なことでした。VQEのような何度も量子と古典を行き来しながら量子回路のパラメタを更新するようなプログラムの場合、これまではパラメタの更新の度にローカル環境とIBM Quantum systemとの間で通信が必要でしたが、Qiskit Runtimeによってこれが解消されるので大幅な計算時間の短縮が期待できます。但し、投げた一連のジョブが他のユーザのジョブと競合し、IBM Quantum system側でキュー待ちの遅延が度々おきてしまっては、せっかくのQiskit Runtimeの恩恵に預かることができません。一度システム上でプログラムが走り出すと、反復計算が終わるまで他のジョブに割り込まれること無く一気にジョブが流れて欲しいと思います。

そこでQiskit Runtimeではこれを可能にするSessionクラス[※8]が用意されています。Sessionを用いると、ジョブの集まりをグループ化し量子コンピュータのジョブスケジューラに優先的に実行させることができます。これによって、同じ量子コンピュータで実行される他のユーザのジョブによる遅延を防ぐことができます。Session内でPrimitiveを利用する場合、以下のようなコードになります。

from qiskit_ibm_runtime import Session

with Session(service, backend=backend) as session:
    sampler = Sampler(session=session)
    job = sampler.run(qc)
    print(f">>> Job ID: {job.job_id()}")

Sessionがアクティブな間は、同じSession内の後続のジョブが他のキューに入ったジョブより優先されます。
但し、キャリブレーションなどのシステムジョブは、Sessionジョブより優先される[※8]そうですので、実行時間(max_execution_time)には留意する必要があるかも知れません。（計算の途中でキャリブレーションが怒ると、その前後で計算結果に影響があらわれるかもしれません）

[※8] https://qiskit.org/ecosystem/ibm-runtime/sessions.html

1-3. Session内でQiskit RuntimeでのPrimitiveを実行: Samplerの場合

まとめとして、Sessionクラスを用いてQiskit RuntimeのSamplerを実行してみます。上述の内容から、Sessionを使うと複数の量子回路のSamplerジョブを一気に走らせることができそうです。回路はsampler.run([qc])としてリストで渡すことが出来ます。以下に、量子ビットの数が2〜5のGHZ状態(量子ビット数が2の場合はBelle状態と呼ばれます)を実行するコードを示します。まず、量子回路を定義する関数を以下のようにします。

def GHZ(num_qubits):
    from qiskit import QuantumCircuit
    
    qc = QuantumCircuit(num_qubits,num_qubits)
    qc.h(0)    
    for i in range(0,num_qubits-1):
        qc.cx(i,i+1)
    qc.barrier()
    
    mes = [i for i in range(num_qubits)]
    qc.measure(mes, mes)
    
    return qc

続いて、Session内でQiskit Runtimeを用いたSamplerを実行するコードを示します。

from qiskit_ibm_runtime import QiskitRuntimeService, Session, Sampler

#Qiskit Runtime サービスの初期化
hub = "ibm-q"
group = "open"
project = "main"
service = QiskitRuntimeService(channel="ibm_quantum", instance=f"{hub}/{group}/{project}") 

# Qiskit Runtimeを利用するbackendの指定
backend = service.backend("ibmq_lima")

# Samplerで計算したい量子回路の構築
qc = []
for i in range(2,6):
    qc.append(GHZ(i))  # 量子ビットが2〜5個の各GHZ状態を生成する量子回路をリストに格納

# Session内でQiskit Runtime Samplerを実行
with Session(service, backend=backend) as session:
    sampler = Sampler(session=session)
    job = sampler.run(qc)   # Session内で複数のGHZ状態を実行する
    job_id = job.job_id()
    
# 結果の回収(ジョブが終了すると結果をresultに格納します)
print(job_id)
retrieve_job = service.job(job_id)

if retrieve_job.status().name == 'QUEUED' or retrieve_job.status().name == 'RUNNING':
    print(retrieve_job.status())
    pass
else:
    print(retrieve_job.status())
    retrieve_job = service.job(job_id)
    result = retrieve_job.result()
    

結果を表示するコードを以下のように作ってみます。

# 結果の確認
res = []
for i in range(0,4):
    res.append(result.quasi_dists[i])

q_index = []
quasi_prop = []
for j in range(0,4):
    q_index.append([i for i in res[j]])
    quasi_prop.append([res[j][i] for i in res[j]])

    
import matplotlib.pyplot as plt 
fig = plt.figure()

ax1 = fig.add_subplot(2, 2, 1)
ax1.bar(q_index[0], quasi_prop[0],label = '#qubit_num: %s '%(2))
ax1.set_xlabel('Index of qubits',fontsize='10')
ax1.set_ylabel('Quasi Probabilities',fontsize='10')
ax1.legend()


ax2 = fig.add_subplot(2, 2, 2)
ax2.set_xlabel('Index of qubits',fontsize='10')
ax2.set_ylabel('Quasi Probabilities',fontsize='10')
ax2.bar(q_index[1], quasi_prop[1],label = '#qubit_num: %s '%(3))
ax2.legend()


ax3 = fig.add_subplot(2, 2, 3)
ax3.set_xlabel('Index of qubits',fontsize='10')
ax3.set_ylabel('Quasi Probabilities',fontsize='10')
ax3.bar(q_index[2], quasi_prop[2],label = '#qubit_num: %s '%(4))
ax3.legend()


ax4 = fig.add_subplot(2, 2, 4)
ax4.set_xlabel('Index of qubits',fontsize='10')
ax4.set_ylabel('Quasi Probabilities',fontsize='10')
ax4.bar(q_index[3], quasi_prop[3],label = '#qubit_num: %s '%(5))
ax4.legend()

fig.tight_layout()

以下のような結果が得られました。

Qiskit Runtimeを用いることにより、実機(ibmq_lima)によるBelle状態および3〜5量子ビットのGHZ状態の結果が、上記コードによる1ジョブで得られました。

Options

補足として、Optionsクラスについて説明します。Primitivesには、さまざまなカテゴリーにグループ化された数々のオプションが用意されています。例えば、エラーの抑止レベル(トランスパイルの最適化などに関する)、軽減レベル(エラー緩和のレベル)、ショット数を指定するには以下のようにします。

from qiskit_ibm_runtime import Options
options = Options()
options.optimization_level = 3  # default = 3
options.resilience_level = 1    # default = 1
options.execution.shots = 2048  # default = 4000

作成したoptionsは、samplerクラスのインスタンスを作成する際に渡すことができます。
これらのオプションは run() を使って計算を行うときに適用されます。

：
with Session(service, backend=backend) as session:
    sampler = Sampler(session=session, options = options)
    job = sampler.run(qc)   
    job_id = job.job_id()
:

エラー抑制のオプションについては[※9]、エラー緩和のオプションについては[※10]が参考になります。
エラー抑制・緩和、あわせてこちらが大変参考になります[※11]。
その他の詳しいオプションについては[※12]にOptionクラスのクラス変数に関する詳細情報があります。

[※9]https://qiskit.org/ecosystem/ibm-runtime/locale/ja_JP/how_to/error-suppression.html
[※10]https://qiskit.org/ecosystem/ibm-runtime/locale/ja_JP/how_to/error-mitigation.html
[※11]https://qiskit.org/ecosystem/ibm-runtime/locale/ja_JP/tutorials/Error-Suppression-and-Error-Mitigation.html
[※12]https://qiskit.org/ecosystem/ibm-runtime/stubs/qiskit_ibm_runtime.options.Options.html#qiskit_ibm_runtime.options.Options

2. ローカル環境でのPrimitiveの利用(その2): Estimatorの場合

Primitveには、上述の擬確率を返すSampler()の他に、期待値を返すEstimator()が用意されています。こちらは、ハミルトニアンなどの期待値をもとに変分最適化行うQEやQAOAなどで重宝されそうです。Esimatorも、Samplerと同様、Esimatorクラスを初期化して実行し、計算したい量子回路(circuits)を渡しますが、Samplerと違い、求めたい期待値の演算子(observables)も渡します。またSessionクラスもSamplerと同様に利用可能です。具体的には以下のようになります。

from qiskit.primitives import Estimator

:
    
with Session(service, backend=backend) as session:
    estimator = Estimator(session=session)
    job = estimator.run(circuits=qc,observables=op)
    print(f">>> Job ID: {job.job_id()}")

簡単な例を見てみます。

3量子ビットの系を考え、各量子ビットに$R_Y(\theta_i) (i=0,1,2)$を作用させた量子回路を考えます。
初期量子ビット$|0\rangle$に$R_Y(\theta_i)$を作用させると,

$$R_Y(\theta_i)|0\rangle = \cos\frac{\theta_i}{2}|0\rangle + \sin\frac{\theta_i}{2}|1\rangle,$$

が得られます。従って、考えている3量子ビット系の量子回路からは

$$\psi =\Big(\cos\frac{\theta_2}{2}|0\rangle + \sin\frac{\theta_2}{2}|1\rangle\Big) \otimes \Big(\cos\frac{\theta_1}{2}|0\rangle + \sin\frac{\theta_1}{2}|1\rangle\Big) \otimes \Big(\cos\frac{\theta_0}{2}|0\rangle + \sin\frac{\theta_0}{2}|1\rangle\Big), $$

が得られます。さて、この量子状態を使って、$ZII$の期待値を考えます。つまり
$$\langle \psi |ZII| \psi\rangle,$$

を計算します。これは3番目の量子ビットにパウリ行列のZ成分が作用しますから、

$$ \langle \psi |ZII| \psi\rangle　= \Big(\cos\frac{\theta_2}{2}\Big)^2-\Big(\sin\frac{\theta_2}{2}\Big)^2,$$

です。例えば$[\theta_0, \theta_1, \theta_2 ]= [0.2,0.3,0.4]$とすると、

$$ \langle \psi |ZII| \psi\rangle　= \Big(\cos(0.2)\Big)^2-\Big(\sin(0.2)\Big)^2=0.921060.. ,$$
となります。この結果をローカル環境のEstimator()を使って確認してみます。

from qiskit.primitives import Estimator
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter 
from qiskit.quantum_info import SparsePauliOp

# psiの作成
qc = QuantumCircuit(3)
for i in range(3):
    theta = Parameter('theta[%s]' % str(i))
    qc.ry(theta, i)
    
# オブザーバブルの作成
op = SparsePauliOp(['ZII'], [1.0])

# estimatorのrun
theta = [[0.2, 0.3, 0.4]]
estimator = Estimator()
job = estimator.run(circuits=qc,observables=op, parameter_values=theta)
result = job.result()
exp_val = result.values
print('exp_val: ', exp_val)

ここで、qiskit.circuitのParameter()を用いて、$R_Y$の角度パラメタthetaを定義しました。このthetaはesimater.run()にparameter_valuesとして渡すことが出来ます。このようにParameter()を用いれば、量子回路のパラメタ値を色々変えて調べる際に便利です。さて結果は

exp_val:  [0.92106099]

となりました。先ほどの計算と一致しているようです。

2-2. Session内でQiskit RuntimeでのPrimitiveを実行(その2): Estimatorの場合

では、Qiskit Runtimeを使って、実機に計算を投げてみます。要領はSamplerと同じです。qiskit.primitivesからqiskit_ibm_runtimeに変更し、バックエンドに実機を指定します。ここではSessionも使ってコードを書いてみます。

from qiskit_ibm_runtime import QiskitRuntimeService, Session, Estimator
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter 
from qiskit.quantum_info import SparsePauliOp

# Qiskit Runtime initialization
hub = "ibm-q"
group = "open"
project = "main"
service = QiskitRuntimeService(channel="ibm_quantum", instance=f"{hub}/{group}/{project}")

# psi
qc = QuantumCircuit(3)
for i in range(3):
    theta = Parameter('theta[%s]' % str(i))
    qc.ry(theta, i)
  
# observable
op = SparsePauliOp(['ZII'], [1.0])

# backend 
backend = service.backend("ibmq_quito")

# parameters
theta = [0.2, 0.3, 0.4]

with Session(service, backend=backend) as session:
    estimator = Estimator(session=session)
    job = estimator.run(circuits=qc,observables=op, parameter_values=theta)
    job_id= job.job_id()
    
# check results
retrieve_job = service.job(job_id)
if retrieve_job.status().name == 'QUEUED' or retrieve_job.status().name == 'RUNNING':
    print(retrieve_job.status())
    pass
else:
    print(retrieve_job.status())
    retrieve_job = service.job(job_id)
    result = retrieve_job.result()
    exp_val = result.values
    print('exp_val: ', exp_val)
    

ここではバックエンドにibmq_quitoを用いました。結果は、

exp_val:  [0.90040915]

となりました。厳密な解析計算の結果と値はズレていますが、かといって大きく外れるわけでもなく大体合っているようです。

水素化ナトリウム分子のエネルギー計算

これまでの手順を踏まえると、論文https://arxiv.org/abs/1905.01534 にある図4の結果を計算できそうな気がしてきます。論文の詳細は省き、図4の計算について簡単に説明します。図4は、図2にあるUCC circuit ansatzと書かれている$R_Z$のパラメタ$\theta$を一つ含む量子回路から生成される量子状態を用いて、水素化ナトリウム(NaH)分子のハミルトニアンの期待値を$\theta$を変化させながらプロットしたものです。NaH分子の基底状態は、エネルギー期待値が最小になる$\theta$で実現します。いわゆるVQEでは適切なオプティマイザでパラメタ$\theta$を最適化しますが、ここでは単に$\theta$を振って、エネルギー期待値がどのように変化するのか見てみます。順を追ってコードを書いてみます。

3-1. UCC circuit ansatz

論文の図2の回路を作ります。

# UCC
import numpy as np
theta = Parameter('theta')

def UCC():
    qc = QuantumCircuit(4)
    qc.rx(np.pi/2,0)
    qc.h(1)
    qc.h(2)
    qc.h(3)
    qc.cx(0,1)
    qc.cx(1,2)
    qc.cx(2,3)

    qc.rz(theta,3)

    qc.cx(2,3)
    qc.cx(1,2)
    qc.cx(0,1)
    qc.rx(-np.pi/2,0)
    qc.h(1)
    qc.h(2)
    qc.h(3)
    
    return qc

確認します。

qc = UCC()
qc.draw('mpl')

以下のような量子回路ができました。

3-2. 水素化ナトリウム分子のハミルトニアン

次にオブザーバブルのハミルトニアンを用意します。これはqiskitにあるqiskit_natureを用いて用意します。

必要なパッケージをインストールします。

pip install qiskit-nature[pyscf]

問題なくインストールできれば、qiskit-natureパッケージを用いて、以下のようなコードで4量子ビット用のNaH分子(原子間距離は論文の1.914388Åに合わせました)のハミルトニアンが得られます。ここではqiskit-natureの詳細については省きますが、qiskit-natureについては[※13]が参考になります。

[※13] https://qiskit.org/ecosystem/nature/locale/ja_JP/tutorials/index.html

from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer
from qiskit_nature.second_q.mappers import JordanWignerMapper, QubitMapper

from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
from qiskit_nature.second_q.algorithms import GroundStateEigensolver

driver = PySCFDriver(
    atom="Na 0 0 0; H 0 0 1.914388",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)

problem = driver.run()

# 2電子2軌道(4qubitの計算)にモデル化
active_space_trafo = as_transformer = ActiveSpaceTransformer(2, 2,active_orbitals=[5,6])

model = active_space_trafo.transform(problem)
hamiltonian = model.hamiltonian
nuclear_repulsion_energy = hamiltonian.nuclear_repulsion_energy #核間反発エネルギー

solver = GroundStateEigensolver(JordanWignerMapper(),NumPyMinimumEigensolver())
extracted_transformer_energy = solver.solve(model).extracted_transformer_energy #modelによる補正エネルギー

second_q_op = hamiltonian.second_q_op() #第２量子化されたハミルトニアン
mapper = JordanWignerMapper() # Jordan-Wigner変換
qubit_jw_op = mapper.map(second_q_op) # qubitハミルトニアン

NaH分子は電子が全部で12個あります。したがって、下から順に軌道を詰めていくと6つの軌道が占有されます。
全ての電子配置を考えると厳密解が得られますが、ここでは4量子ビットの計算をしますので、2軌道に限定した
ようなモデルを考えています。具体的には5番目と6番目(軌道の通し番号を"0"から開始しています)の2軌道のみ
考えます。

3-3. Qiskit Runtimeを用いた各θのエネルギー計算

量子回路とオブザーバブルが準備できました。あとは$\theta$を振ってエネルギー(ハミルトニアンの期待値)計算を行います。まずはローカル環境でのシミュレータ(statevector simulator)で動作を確認します。

#ローカルでrun
from qiskit.primitives import Estimator

# quantum circuit
qc = []

# observable
op = []

# parameters
theta0 = []

# 量子回路、オブザーバブルを各パラメタ毎にリスト化
num = 100
for i in range(1,num+1):
    qc.append(UCC())
    op.append(qubit_jw_op)
    theta0.append([-np.pi + 2*np.pi*(i-1)/(num-1)])
    
estimator = Estimator()
job = estimator.run(circuits=qc,observables=op, parameter_values=theta0)
result = job.result()

total_energy = result.values+extracted_transformer_energy+ nuclear_repulsion_energy  

# 結果のプロット
import matplotlib.pyplot as plt

plt.plot(theta0,total_energy)

結果は以下のようになりました。

論文の図4に近い結果が得られています。では、これまでの手順をふまえ、Qiskit Runtimeを用いた計算を行います。

from qiskit_ibm_runtime import QiskitRuntimeService, Session, Estimator
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter 
from qiskit.quantum_info import SparsePauliOp

# Qiskit Runtime initialization
hub = "ibm-q"
group = "open"
project = "main"
service = QiskitRuntimeService(channel="ibm_quantum", instance=f"{hub}/{group}/{project}")

# quantum circuit
qc = []

# observable
op = []

# backend 
backend = service.backend("ibmq_qasm_simulator")

# parameters
theta = []

# 量子回路、オブザーバブルを各パラメタ毎にリスト化
num = 21
for i in range(1, num+1):
    qc.append(UCC())
    op.append(qubit_jw_op)
    theta.append([-np.pi + 2*np.pi*(i-1)/(num-1)])
    

# Sessionによる実行
with Session(service, backend=backend) as session:
    estimator = Estimator(session=session)
    job = estimator.run(circuits=qc,observables=op, parameter_values=theta)
    job_id= job.job_id()
    
# check results
retrieve_job = service.job(job_id)
if retrieve_job.status().name == 'QUEUED' or retrieve_job.status().name == 'RUNNING':
    print(retrieve_job.status())
    pass
else:
    print(retrieve_job.status())
    retrieve_job = service.job(job_id)
    result = retrieve_job.result()
    exp_val = result.values
    
    total_energy_qasm = result.values + extracted_transformer_energy + nuclear_repulsion_energy

実機"ibmq_lima"でも実行してみます。これは先ほどのコードで

:
# backend 
backend = service.backend("ibmq_lima")
:

total_energy_lima = result.values + extracted_transformer_energy + nuclear_repulsion_energy

と変更するだけです。結果は以下のようになりました。

実機でもシミュレータの結果と非常に良く合っていて驚きます。Estimatorでは、デフォルトでoptions.resilience_level=1のreadout_error_mitigationが適用されていますので、このように良い一致が見られていると考えられます。前述のOptionの設定を色々変えて(options.resilience_level=0にしてみる等)計算してみると、面白いかもしれません。

まとめ

今回、Qiskit Runtimeについて説明してみました。Qiskit Runtimeを用いるとIBM Quantum上で高速な量子古典ハイブリッドの計算が可能になります。またQiskit Runtimeに用意されているSessionを用いるとジョブを束ねて計算することが出来ます。コードを書く際にはPrimitvesという便利なクラスが用意されており、これらとRuntimeを上手く使うことで効率的な開発が進められます。本記事では触れることは出来ませんでしたが、Primitiveを使ってRuntime用のプログラムを書き、IBM Quantumへuploadし実行させることも可能です。

Qiskit Runtimeの威力が発揮されるのは、やはりVQEのような何度も反復計算を行うようなアルゴリズムと思います。Qiskit Runtimeのさらなる応用例については[※14]が大変参考になります。

[※14]https://qiskit.org/ecosystem/ibm-runtime/tutorials.html