accessing_higher_energy_states.ipynb

Kernel: Python 3

Qiskit Pulseで高エネルギー状態へアクセスする

ほとんどの量子アルゴリズム/アプリケーションでは、$|0\rangle$と$|1\rangle$によって張られた2次元空間で計算が実行されます。ただし、IBMのハードウェアでは、通常は使用されない、より高いエネルギー状態も存在します。このセクションでは、Qiskit Pulseを使ってこれらの状態を探索することにフォーカスを当てます。特に、$|2\rangle$ 状態を励起し、$|0\rangle$、$|1\rangle$、$|2\rangle$の状態を分類するための識別器を作成する方法を示します。

このノートブックを読む前に、前の章を読むことをお勧めします。また、Qiskit Pulseのスペック（Ref 1）も読むことをお勧めします。

物理学的背景

ここで、IBMの量子ハードウェアの多くの基礎となっている、トランズモンキュービットの物理学的な背景を説明します。このシステムには、ジョセフソン接合とコンデンサーで構成される超伝導回路が含まれています。超伝導回路に不慣れな方は、こちらのレビュー (Ref. 2)を参照してください。このシステムのハミルトニアンは以下で与えられます。

$$H = 4 E_C n^2 - E_J \cos(\phi),$$

ここで、$E_C, E_J$はコンデンサーのエネルギーとジョセフソンエネルギーを示し、$n$は減衰した電荷数演算子で、$\phi$はジャンクションのために減衰した磁束です。$\hbar=1$として扱います。

トランズモンキュービットは$\phi$が小さい領域で定義されるため、$E_J \cos(\phi)$をテイラー級数で展開できます（定数項を無視します）。

$$E_J \cos(\phi) \approx \frac{1}{2} E_J \phi^2 - \frac{1}{24} E_J \phi^4 + \mathcal{O}(\phi^6).$$

$\phi$の二次の項$\phi^2$は、標準の調和振動子を定義します。その他の追加の項はそれぞれ非調和性をもたらします。

$n \sim (a-a^\dagger), \phi \sim (a+a^\dagger)$の関係を使うと（$a^\dagger,a$は生成消滅演算子）、システムは以下のハミルトニアンを持つダフィング（Duffing）振動子に似ていることを示せます。

$$H = \omega a^\dagger a + \frac{\alpha}{2} a^\dagger a^\dagger a a,$$

$\omega$は、$0\rightarrow1$の励起周波数($\omega \equiv \omega^{0\rightarrow1}$)を与え、$\alpha$は$0\rightarrow1$の周波数と$1\rightarrow2$の周波数の間の非調和です。必要に応じて駆動の条件を追加できます。

標準の2次元部分空間へ特化したい場合は、$|\alpha|$ を十分に大きくとるか、高エネルギー状態を抑制する特別な制御テクニックを使います。

目次

0. はじめに

1. 0と1の状態の識別

   1. 0->1の周波数スイープ

   2. 0->1のラビ実験

   3. 0,1の識別器を構築する

2. 0,1,2の状態の識別

   1. 1->2の周波数の計算

      1. サイドバンド法を使った1->2の周波数スイープ

   2. 1->2のラビ実験

   3. 0,1,2の識別器を構築する

3. 参考文献

0. はじめに

まず、依存関係をインポートし、いくつかのデフォルトの変数を定義します。量子ビット0を実験に使います。公開されている単一量子ビットデバイスであるibmq_armonkで実験を行います。

import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import curve_fit
from scipy.signal import find_peaks

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import train_test_split

import qiskit.pulse as pulse
import qiskit.pulse.library as pulse_lib
from qiskit.compiler import assemble
from qiskit.pulse.library import SamplePulse

from qiskit.tools.monitor import job_monitor

import warnings
warnings.filterwarnings('ignore')
from qiskit.tools.jupyter import *
%matplotlib inline

from qiskit import IBMQ
IBMQ.load_account()
provider = IBMQ.get_provider(hub='ibm-q', group='open', project='main')
backend = provider.get_backend('ibmq_armonk')

backend_config = backend.configuration()
assert backend_config.open_pulse, "Backend doesn't support Pulse"

dt = backend_config.dt

backend_defaults = backend.defaults()

# 単位変換係数 -> すべてのバックエンドのプロパティーがSI単位系(Hz, sec, etc)で返される
GHz = 1.0e9 # Gigahertz
MHz = 1.0e6 # Megahertz
us = 1.0e-6 # Microseconds
ns = 1.0e-9 # Nanoseconds

qubit = 0 # 分析に使う量子ビット
default_qubit_freq = backend_defaults.qubit_freq_est[qubit] # デフォルトの量子ビット周波数単位はHz 
print(f"Qubit {qubit} has an estimated frequency of {default_qubit_freq/ GHz} GHz.")

#（各デバイスに固有の）データをスケーリング
scale_factor = 1e-14

# 実験のショット回数
NUM_SHOTS = 1024

### 必要なチャネルを収集する
drive_chan = pulse.DriveChannel(qubit)
meas_chan = pulse.MeasureChannel(qubit)
acq_chan = pulse.AcquireChannel(qubit)

Qubit 0 has an estimated frequency of 4.9744534273002605 GHz.

いくつか便利な関数を追加で定義します。

def get_job_data(job, average):
    """すでに実行されているジョブからデータを取得します。
    引数:
        job (Job): データが必要なジョブ
        average (bool): Trueの場合、データが平均であると想定してデータを取得。
                        Falseの場合、シングルショット用と想定してデータを取得。
    返し値:
        list: ジョブの結果データを含むリスト
    """
    job_results = job.result(timeout=120) # タイムアウトパラメーターは120秒にセット
    result_data = []
    for i in range(len(job_results.results)):
        if average: # 平均データを得る
            result_data.append(job_results.get_memory(i)[qubit]*scale_factor) 
        else: # シングルデータを得る
            result_data.append(job_results.get_memory(i)[:, qubit]*scale_factor)  
    return result_data

def get_closest_multiple_of_16(num):
    """16の倍数に最も近いものを計算します。
    パルスが使えるデバイスが16サンプルの倍数の期間が必要なためです。
    """
    return (int(num) - (int(num)%16))

次に、駆動パルスと測定のためのいくつかのデフォルトパラメーターを含めます。命令スケジュールマップから（バックエンドデフォルトから）measureコマンドをプルして、新しいキャリブレーションでアップデートされるようにします。

# 駆動パルスのパラメーター (us = マイクロ秒)
drive_sigma_us = 0.075                     # ガウシアンの実際の幅を決めます
drive_samples_us = drive_sigma_us*8         # 切り捨てパラメーター
                                          # ガウシアンには自然な有限長がないためです。

drive_sigma = get_closest_multiple_of_16(drive_sigma_us * us /dt)       # ガウシアンの幅の単位はdt
drive_samples = get_closest_multiple_of_16(drive_samples_us * us /dt)   # 切り捨てパラメーターの単位はdt

# この量子ビットに必要な測定マップインデックスを見つける
meas_map_idx = None
for i, measure_group in enumerate(backend_config.meas_map):
    if qubit in measure_group:
        meas_map_idx = i
        break
assert meas_map_idx is not None, f"Couldn't find qubit {qubit} in the meas_map!"

# 命令スケジュールマップからデフォルトの測定パルスを取得
inst_sched_map = backend_defaults.instruction_schedule_map
measure = inst_sched_map.get('measure', qubits=backend_config.meas_map[meas_map_idx])

1. $|0\rangle$ と $|1\rangle$の状態の識別

このセクションでは、標準の$|0\rangle$と$|1\rangle$の状態の識別器を構築します。識別器のジョブは、meas_level=1の複素数データを取得し、標準の$|0\rangle$の$|1\rangle$の状態(meas_level=2)に分類することです。これは、前の章の多くと同じ作業です。この結果は、このNotebookがフォーカスしている高エネルギー状態に励起するために必要です。

1A. 0->1 周波数のスイープ

識別器の構築の最初のステップは、前の章でやったのと同じように、我々の量子ビット周波数をキャリブレーションすることです。

def create_ground_freq_sweep_program(freqs, drive_power):
    """基底状態を励起して周波数掃引を行うプログラムを作成します。
     ドライブパワーに応じて、これは0->1の周波数なのか、または0->2の周波数なのかを明らかにすることができます。
     引数：
         freqs (np.ndarray(dtype=float))：スイープする周波数のNumpy配列。
         drive_power (float)：ドライブ振幅の値。
     レイズ：
         ValueError：75を超える頻度を使用すると発生します。
         　　　　　　　現在、これを実行しようとすると、バックエンドでエラーが投げられます。
     戻り値：
         Qobj：基底状態の周波数掃引実験のプログラム。
    """
    if len(freqs) > 75:
        raise ValueError("You can only run 75 schedules at a time.")
    
    # スイープ情報を表示
    print(f"The frequency sweep will go from {freqs[0] / GHz} GHz to {freqs[-1]/ GHz} GHz \
using {len(freqs)} frequencies. The drive power is {drive_power}.")
    
    # 駆動パルスを定義
    ground_sweep_drive_pulse = pulse_lib.gaussian(duration=drive_samples,
                                                  sigma=drive_sigma,
                                                  amp=drive_power,
                                                  name='ground_sweep_drive_pulse')
    # スイープのための周波数を定義
    schedule = pulse.Schedule(name='Frequency sweep starting from ground state.')
    
    schedule |= pulse.Play(ground_sweep_drive_pulse, drive_chan)
    schedule |= measure << schedule.duration
    
    # define frequencies for the sweep
    schedule_freqs = [{drive_chan: freq} for freq in freqs]

    # プログラムを組み立てる
    # 注：それぞれが同じことを行うため、必要なスケジュールは1つだけです;
    # スケジュールごとに、ドライブをミックスダウンするLO周波数が変化します
    # これにより周波数掃引が可能になります
    ground_freq_sweep_program = assemble(schedule,
                                         backend=backend, 
                                         meas_level=1,
                                         meas_return='avg',
                                         shots=NUM_SHOTS,
                                         schedule_los=schedule_freqs)
    
    return ground_freq_sweep_program

# 75個の周波数で推定周波数の周りに40MHzを掃引します
num_freqs = 75
ground_sweep_freqs = default_qubit_freq + np.linspace(-20*MHz, 20*MHz, num_freqs)
ground_freq_sweep_program = create_ground_freq_sweep_program(ground_sweep_freqs, drive_power=0.3)

The frequency sweep will go from 4.954453427300261 GHz to 4.994453427300261 GHz using 75 frequencies. The drive power is 0.3.

ground_freq_sweep_job = backend.run(ground_freq_sweep_program)

print(ground_freq_sweep_job.job_id())
job_monitor(ground_freq_sweep_job)

5f9683d461672700130444cd
Job Status: job has successfully run

# ジョブのデータ（平均）を取得する
ground_freq_sweep_data = get_job_data(ground_freq_sweep_job, average=True)

データをローレンツ曲線に適合させ、キャリブレーションされた周波数を抽出します。

def fit_function(x_values, y_values, function, init_params):
    """Fit a function using scipy curve_fit."""
    fitparams, conv = curve_fit(function, x_values, y_values, init_params)
    y_fit = function(x_values, *fitparams)
    
    return fitparams, y_fit

# Hz単位でのフィッティングをします
(ground_sweep_fit_params, 
 ground_sweep_y_fit) = fit_function(ground_sweep_freqs,
                                   ground_freq_sweep_data, 
                                   lambda x, A, q_freq, B, C: (A / np.pi) * (B / ((x - q_freq)**2 + B**2)) + C,
                                   [7, 4.975*GHz, 1*GHz, 3*GHz] # フィッティングのための初期パラメーター
                                   )

# 注：シグナルの実数部のみをプロットしています
plt.scatter(ground_sweep_freqs/GHz, ground_freq_sweep_data, color='black')
plt.plot(ground_sweep_freqs/GHz, ground_sweep_y_fit, color='red')
plt.xlim([min(ground_sweep_freqs/GHz), max(ground_sweep_freqs/GHz)])
plt.xlabel("Frequency [GHz]", fontsize=15)
plt.ylabel("Measured Signal [a.u.]", fontsize=15)
plt.title("0->1 Frequency Sweep", fontsize=15)
plt.show()

_, cal_qubit_freq, _, _ = ground_sweep_fit_params
print(f"We've updated our qubit frequency estimate from "
      f"{round(default_qubit_freq/GHz, 7)} GHz to {round(cal_qubit_freq/GHz, 7)} GHz.")

We've updated our qubit frequency estimate from 4.9744534 GHz to 4.974489 GHz.

1B. 0->1 のラビ実験

次に、$0\rightarrow1 ~ \pi$パルスの振幅を計算するラビ実験を実行します。$\pi$パルスは、$|0\rangle$から$|1\rangle$の状態へ移動させるパルス（ブロッホ球上での$\pi$回転）だということを思い出してください。

# 実験の構成
num_rabi_points = 50 # 実験の数（つまり、掃引の振幅）

# 反復する駆動パルスの振幅値：0から0.75まで等間隔に配置された50の振幅
drive_amp_min = 0
drive_amp_max = 0.75
drive_amps = np.linspace(drive_amp_min, drive_amp_max, num_rabi_points)

# スケジュールを作成
rabi_01_schedules = []
# 駆動振幅すべてにわたってループ
for ii, drive_amp in enumerate(drive_amps):
    # 駆動パルス
    rabi_01_pulse = pulse_lib.gaussian(duration=drive_samples, 
                                       amp=drive_amp, 
                                       sigma=drive_sigma, 
                                       name='rabi_01_pulse_%d' % ii)
    
    # スケジュールにコマンドを追加
    schedule = pulse.Schedule(name='Rabi Experiment at drive amp = %s' % drive_amp)
    schedule |= pulse.Play(rabi_01_pulse, drive_chan)
    schedule |= measure << schedule.duration # 測定をドライブパルスの後にシフト
    rabi_01_schedules.append(schedule)

# プログラムにスケジュールを組み込む
# 注：較正された周波数で駆動します。
rabi_01_expt_program = assemble(rabi_01_schedules,
                                backend=backend,
                                meas_level=1,
                                meas_return='avg',
                                shots=NUM_SHOTS,
                                schedule_los=[{drive_chan: cal_qubit_freq}]
                                               * num_rabi_points)

rabi_01_job = backend.run(rabi_01_expt_program)

print(rabi_01_job.job_id())
job_monitor(rabi_01_job)

5f96870a834560001306fd69
Job Status: job has successfully run

# ジョブのデータ（平均）を取得する
rabi_01_data = get_job_data(rabi_01_job, average=True)

def baseline_remove(values):
    """Center data around 0."""
    return np.array(values) - np.mean(values)

# 注：データの実数部のみがプロットされます
rabi_01_data = np.real(baseline_remove(rabi_01_data))
(rabi_01_fit_params, 
 rabi_01_y_fit) = fit_function(drive_amps,
                               rabi_01_data, 
                               lambda x, A, B, drive_01_period, phi: (A*np.cos(2*np.pi*x/drive_01_period - phi) + B),
                               [4, -4, 0.5, 0])

plt.scatter(drive_amps, rabi_01_data, color='black')
plt.plot(drive_amps, rabi_01_y_fit, color='red')

drive_01_period = rabi_01_fit_params[2] 
# piの振幅計算でphiを計算
pi_amp_01 = (drive_01_period/2/np.pi) *(np.pi+rabi_01_fit_params[3])

plt.axvline(pi_amp_01, color='red', linestyle='--')
plt.axvline(pi_amp_01+drive_01_period/2, color='red', linestyle='--')
plt.annotate("", xy=(pi_amp_01+drive_01_period/2, 0), xytext=(pi_amp_01,0), arrowprops=dict(arrowstyle="<->", color='red'))
plt.annotate("$\pi$", xy=(pi_amp_01-0.03, 0.1), color='red')

plt.xlabel("Drive amp [a.u.]", fontsize=15)
plt.ylabel("Measured signal [a.u.]", fontsize=15)
plt.title('0->1 Rabi Experiment', fontsize=15)
plt.show()

print(f"Pi Amplitude (0->1) = {pi_amp_01}")

Pi Amplitude (0->1) = 0.22320005035637547

この結果を使って、$0\rightarrow1$ $\pi$パルスを定義します。

pi_pulse_01 = pulse_lib.gaussian(duration=drive_samples,
                                 amp=pi_amp_01, 
                                 sigma=drive_sigma,
                                 name='pi_pulse_01')

1C. 0,1 の識別器を構築する

これで、キャリブレーションされた周波数と$\pi$パルスを得たので、$|0\rangle$と$1\rangle$の状態の識別器を構築できます。識別器は、IQ平面においてmeas_level=1のデータを取って、それを$|0\rangle$または$1\rangle$を判別することで機能します。

$|0\rangle$と$|1\rangle$の状態は、IQ平面上で重心として知られているコヒーレントな円形の"ブロブ"を形成します。重心の中心は、各状態の正確なノイズのないIQポイントを定義します。周囲の雲は、様々なノイズ源から生成されたデータの分散を示します。

$|0\rangle$と$|1\rangle$間を識別（判別）するために、機械学習のテクニック、線形判別分析を適用します。この方法は量子ビットの状態を判別する一般的なテクニックです。

最初のステップは、重心データを得ることです。そのために、2つのスケジュールを定義します（システムが$|0\rangle$の状態から始まることを思い出しましょう。）：

1. $|0\rangle$の状態を直接測定します（$|0\rangle$の重心を得ます）。

2. $\pi$パルスを適用して、測定します（$|1\rangle$の重心を得ます）。

# 2つのスケジュールを作る

# 基底状態のスケジュール
zero_schedule = pulse.Schedule(name="zero schedule")
zero_schedule |= measure

# 励起状態のスケジュール
one_schedule = pulse.Schedule(name="one schedule")
one_schedule |= pulse.Play(pi_pulse_01, drive_chan) 
one_schedule |= measure << one_schedule.duration

# スケジュールをプログラムにアセンブルする
IQ_01_program = assemble([zero_schedule, one_schedule],
                          backend=backend,
                          meas_level=1,
                          meas_return='single',
                          shots=NUM_SHOTS,
                          schedule_los=[{drive_chan: cal_qubit_freq}] * 2)

IQ_01_job = backend.run(IQ_01_program)

print(IQ_01_job.job_id())
job_monitor(IQ_01_job)

5f968a606167270013044513
Job Status: job has successfully run

# （単一の）ジョブデータを取得します；0と1に分割します
IQ_01_data = get_job_data(IQ_01_job, average=False)
zero_data = IQ_01_data[0]
one_data = IQ_01_data[1]

def IQ_01_plot(x_min, x_max, y_min, y_max):
    """Helper function for plotting IQ plane for |0>, |1>. Limits of plot given
    as arguments."""
    # 0のデータは青でプロット
    plt.scatter(np.real(zero_data), np.imag(zero_data), 
                    s=5, cmap='viridis', c='blue', alpha=0.5, label=r'$|0\rangle$')
    # 1のデータは赤でプロット
    plt.scatter(np.real(one_data), np.imag(one_data), 
                    s=5, cmap='viridis', c='red', alpha=0.5, label=r'$|1\rangle$')

    # 0状態と1状態の平均に大きなドットをプロットします。
    mean_zero = np.mean(zero_data) # 実部と虚部両方の平均を取ります。
    mean_one = np.mean(one_data)
    plt.scatter(np.real(mean_zero), np.imag(mean_zero), 
                s=200, cmap='viridis', c='black',alpha=1.0)
    plt.scatter(np.real(mean_one), np.imag(mean_one), 
                s=200, cmap='viridis', c='black',alpha=1.0)
    
    plt.xlim(x_min, x_max)
    plt.ylim(y_min,y_max)
    plt.legend()
    plt.ylabel('I [a.u.]', fontsize=15)
    plt.xlabel('Q [a.u.]', fontsize=15)
    plt.title("0-1 discrimination", fontsize=15)

以下のように、IQプロットを表示します。青の重心は$|0\rangle$状態で、赤の重心は$|1\rangle$状態です。（注：プロットが見えないときは、Notebookを再実行してください。）

x_min = -5
x_max = 15
y_min = -5
y_max = 10
IQ_01_plot(x_min, x_max, y_min, y_max)

さて、実際に識別器を構築する時が来ました。先に述べたように、線形判別分析（Linear Discriminant Analysis, LDA）と呼ばれる機械学習のテクニックを使います。LDAは、任意のデータセットをカテゴリーのセット（ここでは$|0\rangle$と$|1\rangle$）に分類するために、各カテゴリーの平均の間の距離を最大化し、各カテゴリーの分散を最小化します。より詳しくは、こちら (Ref. 3)をご覧ください。

LDAは、セパラトリックス(separatrix)と呼ばれるラインを生成します。与えられたデータポイントがどちら側のセパラトリックスにあるかに応じて、それがどのカテゴリーに属しているかを判別できます。我々の場合、セパラトリックスの片側が$|0\rangle$状態で、もう一方の側が$|1\rangle$の状態です。

我々は、最初の半分のデータを学習用に使い、残りの半分をテスト用に使います。LDAの実装のためにscikit.learnを使います：将来のリリースでは、この機能は、Qiskit-Ignisに直接実装されてリリースされる予定です（ここを参照）。

結果データを判別に適したフォーマットになるように再形成します。

def reshape_complex_vec(vec):
    """
    複素数ベクトルvecを取り込んで、実際のimagエントリーを含む2d配列を返します。
    これは学習に必要なデータです。
     Args：
         vec (list)：データの複素数ベクトル
     戻り値：
         list：(real(vec], imag(vec))で指定されたエントリー付きのベクトル
    """
    length = len(vec)
    vec_reshaped = np.zeros((length, 2))
    for i in range(len(vec)):
        vec_reshaped[i]=[np.real(vec[i]), np.imag(vec[i])]
    return vec_reshaped

# IQベクトルを作成します（実部と虚部で構成されています）
zero_data_reshaped = reshape_complex_vec(zero_data)
one_data_reshaped = reshape_complex_vec(one_data)  

IQ_01_data = np.concatenate((zero_data_reshaped, one_data_reshaped))
print(IQ_01_data.shape) # IQデータの形を確認します

(2048, 2)

次に、学習用データとテスト用データを分割します。期待される結果（基底状態のスケジュールは0の配列、励起状態のスケジュールは1の配列）を含む状態ベクトルを使ってテストします。

#（テスト用に）0と1でベクトルを構築する
state_01 = np.zeros(NUM_SHOTS) # shotsは実験の回数
state_01 = np.concatenate((state_01, np.ones(NUM_SHOTS)))
print(len(state_01))

# データをシャッフルしてトレーニングセットとテストセットに分割します
IQ_01_train, IQ_01_test, state_01_train, state_01_test = train_test_split(IQ_01_data, state_01, test_size=0.5)

2048

最後に、モデルを設定して、学習します。学習精度が表示されます。

# LDAをセットアップします
LDA_01 = LinearDiscriminantAnalysis()
LDA_01.fit(IQ_01_train, state_01_train)

LinearDiscriminantAnalysis()

# シンプルなデータでテストします
print(LDA_01.predict([[0,0], [10, 0]]))

[0. 0.]

# 精度を計算します
score_01 = LDA_01.score(IQ_01_test, state_01_test)
print(score_01)

0.939453125

最後のステップは、セパラトリックスをプロットすることです。

# セパラトリックスを表示データの上にプロットします
def separatrixPlot(lda, x_min, x_max, y_min, y_max, shots):
    nx, ny = shots, shots

    xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),
                         np.linspace(y_min, y_max, ny))
    Z = lda.predict_proba(np.c_[xx.ravel(), yy.ravel()])
    Z = Z[:, 1].reshape(xx.shape)

    plt.contour(xx, yy, Z, [0.5], linewidths=2., colors='black')

IQ_01_plot(x_min, x_max, y_min, y_max)
separatrixPlot(LDA_01, x_min, x_max, y_min, y_max, NUM_SHOTS)

セラパトリックスのどちらのサイドがどの重心（つまり状態）に対応しているか確認します。IQ平面上の点が与えられると、このモデルはセラパトリックスのどちらの側にそれが置かれているかチェックし、対応する状態を返します。

2. $|0\rangle$, $|1\rangle$, $|2\rangle$ の状態の識別

$0, 1$の識別器をキャリブレーションしたので、高エネルギー状態の励起に移ります。特に、$|2\rangle$の状態の励起にフォーカスし、$|0\rangle$と$|1\rangle$と$|2\rangle$の状態をそれぞれのIQデータポイントから判別する識別器を構築することに焦点を当てます。さらに高い状態（$|3\rangle$、$|4\rangle$など）の手順も同様ですが、明示的にテストはしていません。

高い状態の識別器を構築する手順は以下の通りです：

1. $1\rightarrow2$周波数を計算します。

2. $1\rightarrow2$のための$\pi$パルスの振幅を得るためにラビ実験を行います。そのためには、まず、$0\rightarrow1$ $\pi$パルスを適用して、$|0\rangle$から$|1\rangle$の状態にします。次に、上記で得た$1\rightarrow2$周波数において、駆動振幅のスイープを行います。

3. 3つのスケジュールを構成します： a. 0スケジュール：基底状態を測定するだけです。 b. 1スケジュール：$0\rightarrow1$ $\pi$パルスを適用し、測定します。 c. 2スケジュール：$0\rightarrow1$ $\pi$パルスを適用し、次に$1\rightarrow2$ $\pi$パルスを適用し測定します。

4. 各スケジュールのデータを学習用データとテスト用データのセットに分け、判別用のLDAモデルを構築します。

2A. 1->2 周波数の計算

キャリブレーションの最初のステップは、$1\rightarrow2$ の状態に移行するために必要な周波数を計算することです。これを行うには2つの方法があります：

1. 基底状態から周波数をスイープし、非常に高い電力をかけます。印加電力が十分に高い場合には、2つのピークが観測されます。1つはセクション 1で見つかった $0\rightarrow1$周波数で、もう一つは、$0\rightarrow2$周波数です。$1\rightarrow2$周波数は２つの差を取ることで得られます。残念ながら、ibmq_armonkでは、最大駆動電力$1.0$はこの遷移を起こすのに十分ではありません。代わりに、２番目の方法を使います。

2. $0\rightarrow1$ $\pi$パルスを適用して、$|1\rangle$状態を励起します。その後、$|1\rangle$状態のさらに上の励起に対して、周波数スイープを実行します。$0\rightarrow1$周波数より低いところで、$1\rightarrow2$周波数に対応した単一ピークが観測されるはずです。

サイドバンド法を使用した1->2 の周波数スイープ

上記の２番目の方法に従いましょう。$0\rightarrow 1$ $\pi$パルスを駆動するために、ローカル共振（local oscilattor, LO）周波数が必要です。これは、キャリブレーションされた$0\rightarrow1$周波数cal_qubit_freq（セクション1のラビ$\pi$パルスの構築を参照）によって与えられます。ただし、$1\rightarrow2$周波数の範囲をスイープするために、LO周波数を変化させる必要があります。残念ながら、Pulseのスペックでは、各スケジュールごとに、一つのLO周波数が必要です。

これを解決するには、LO周波数をcal_qubit_freqにセットし、をfreq-cal_qubit_freqで$1\rightarrow2$パルスの上にサイン関数を乗算します。ここでfreqは目的のスキャン周波数です。知られているように、正弦波サイドバンドを適用すると、プログラムのアセンブル時に手動で設定せずにLO周波数を変更可能です。

def apply_sideband(pulse, freq):
    """freq周波数でこのパルスに正弦波サイドバンドを適用します。
     引数：
         pulse (SamplePulse)：対象のパルス。
         freq (float)：スイープを適用するLO周波数。
     戻り値：
         SamplePulse：サイドバンドが適用されたパルス（freqとcal_qubit_freqの差で振動します）。
    """
    # 時間は0からdt*drive_samplesで、2*pi*f*tの形の正弦波引数になります
    t_samples = np.linspace(0, dt*drive_samples, drive_samples)
    sine_pulse = np.sin(2*np.pi*(freq-cal_qubit_freq)*t_samples) # no amp for the sine
    
    # サイドバンドが適用されたサンプルパルスを作成
    # 注：sq_pulse.samplesを実数にし、要素ごとに乗算する必要があります
    sideband_pulse = SamplePulse(np.multiply(np.real(pulse.samples), sine_pulse), name='sideband_pulse')
    
    return sideband_pulse    

プログラムをアセンブルするためのロジックをメソッドにラップして、プログラムを実行します。

def create_excited_freq_sweep_program(freqs, drive_power):
    """|1>状態を励起することにより、周波数掃引を行うプログラムを作成します。
     これにより、1-> 2の周波数を取得できます。
     較正された量子ビット周波数を使用して、piパルスを介して|0>から|1>の状態になります。
     |1>から|2>への周波数掃引を行うには、正弦係数を掃引駆動パルスに追加することにより、サイドバンド法を使用します。
     引数：
         freqs (np.ndarray(dtype=float))：掃引周波数のNumpy配列。
         drive_power (float)：駆動振幅の値。
     レイズ：
         ValueError：75を超える頻度を使用するとスローされます； 現在、75個を超える周波数を試行すると、
         　　　　　　　バックエンドでエラーがスローされます。
     戻り値：
         Qobj：周波数掃引実験用のプログラム。
    """
    if len(freqs) > 75:
        raise ValueError("You can only run 75 schedules at a time.")
        
    print(f"The frequency sweep will go from {freqs[0] / GHz} GHz to {freqs[-1]/ GHz} GHz \
using {len(freqs)} frequencies. The drive power is {drive_power}.")

    base_12_pulse = pulse_lib.gaussian(duration=drive_samples,
                                        sigma=drive_sigma,
                                        amp=drive_power,
                                        name='base_12_pulse')
    schedules = []
    for jj, freq in enumerate(freqs):
        
        # ガウシアンパルスにサイドバンドを追加
        freq_sweep_12_pulse = apply_sideband(base_12_pulse, freq)
        
        # スケジュールのコマンドを追加
        schedule = pulse.Schedule(name="Frequency = {}".format(freq))

        # 0->1のパルス、掃引パルスの周波数、測定を追加
        schedule |= pulse.Play(pi_pulse_01, drive_chan)
        schedule |= pulse.Play(freq_sweep_12_pulse, drive_chan) << schedule.duration 
        schedule |= measure << schedule.duration # 駆動パルスの後に測定をシフト

        schedules.append(schedule)

    num_freqs = len(freqs)
    
    # スケジュールを表示します
    display(schedules[-1].draw(channels=[drive_chan, meas_chan], label=True, scale=1.0))
    
    # 周波数掃引プログラムを組み込みます
    # 注：LOは各スケジュールでのcal_qubit_freqです；サイドバンドによって組み込みます
    excited_freq_sweep_program = assemble(schedules,
                                          backend=backend, 
                                          meas_level=1,
                                          meas_return='avg',
                                          shots=NUM_SHOTS,
                                          schedule_los=[{drive_chan: cal_qubit_freq}]
                                                         * num_freqs)
    
    return excited_freq_sweep_program

# 0->1周波数より下で1->2の周波数を見つけるために400 MHzを掃引します
num_freqs = 75
excited_sweep_freqs = cal_qubit_freq + np.linspace(-400*MHz, 30*MHz, num_freqs)
excited_freq_sweep_program = create_excited_freq_sweep_program(excited_sweep_freqs, drive_power=0.3)

# 確認のためにスケジュールの一例をプロットします

The frequency sweep will go from 4.574489007795566 GHz to 5.0044890077955655 GHz using 75 frequencies. The drive power is 0.3.

excited_freq_sweep_job = backend.run(excited_freq_sweep_program)

print(excited_freq_sweep_job.job_id())
job_monitor(excited_freq_sweep_job)

5f968da3183f64001367e5c7
Job Status: job has successfully run

# (平均の)ジョブデータを取得します
excited_freq_sweep_data = get_job_data(excited_freq_sweep_job, average=True)

# 注：シグナルの実部だけをプロットします
plt.scatter(excited_sweep_freqs/GHz, excited_freq_sweep_data, color='black')
plt.xlim([min(excited_sweep_freqs/GHz)+0.01, max(excited_sweep_freqs/GHz)]) # ignore min point (is off)
plt.xlabel("Frequency [GHz]", fontsize=15)
plt.ylabel("Measured Signal [a.u.]", fontsize=15)
plt.title("1->2 Frequency Sweep (first pass)", fontsize=15)
plt.show()

最小値が$4.64$ GHz近辺に見られます。いくつかの偽の最大値がありますが、それらは、$1\rightarrow2$周波数には大きすぎます。最小値が$1\rightarrow2$周波数に対応します。

相対最小関数を使って、この点の値を正確に計算します。これで、$1\rightarrow2$周波数の推定値が得られます。

# output_dataに相対的最小周波数を表示します；高さは下限（絶対値）を示します
def rel_maxima(freqs, output_data, height): 
    """output_dataに相対的な最小周波数を出力します（ピークを確認できます）；
    高さは上限（絶対値）を示します。
    高さを正しく設定しないと、ピークが無視されます。
    引数：
         freqs (list)：周波数リスト
         output_data (list)：結果のシグナルのリスト
         height (float)：ピークの上限（絶対値）
     戻り値：
         list：相対的な最小周波数を含むリスト
    """
    peaks, _ = find_peaks(output_data, height)
    print("Freq. dips: ", freqs[peaks])
    return freqs[peaks]

maxima = rel_maxima(excited_sweep_freqs, np.real(excited_freq_sweep_data), 10)
approx_12_freq = maxima

Freq. dips:  [4.62678631e+09]

上記で得られた推定値を使って、より正確な掃引を行います（つまり、大幅に狭い範囲で掃引を行います）。これによって、$1\rightarrow2$周波数のより正確な値を得ることができます。上下$20$ MHzをスイープします。

# 狭い範囲での掃引
num_freqs = 75
refined_excited_sweep_freqs = approx_12_freq + np.linspace(-20*MHz, 20*MHz, num_freqs)
refined_excited_freq_sweep_program = create_excited_freq_sweep_program(refined_excited_sweep_freqs, drive_power=0.3)

The frequency sweep will go from 4.6067863050928635 GHz to 4.6467863050928635 GHz using 75 frequencies. The drive power is 0.3.

refined_excited_freq_sweep_job = backend.run(refined_excited_freq_sweep_program)

print(refined_excited_freq_sweep_job.job_id())
job_monitor(refined_excited_freq_sweep_job)

5f9690fc7b787b00147e8cbe
Job Status: job has successfully run

# より正確な（平均）データを取得する
refined_excited_freq_sweep_data = get_job_data(refined_excited_freq_sweep_job, average=True)

標準ローレンツ曲線を用いて、このより正確な信号をプロットしてフィットします。

# Hzの単位でフィッティングする 
(refined_excited_sweep_fit_params, 
 refined_excited_sweep_y_fit) = fit_function(refined_excited_sweep_freqs,
                                     refined_excited_freq_sweep_data, 
                                     lambda x, A, q_freq, B, C: (A / np.pi) * (B / ((x - q_freq)**2 + B**2)) + C,
                                     [-12, 4.625*GHz, 0.05*GHz, 3*GHz] # フィッティングのための初期パラメーター
                                     )

# 注：シグナルの実数部のみをプロットしています
plt.scatter(refined_excited_sweep_freqs/GHz, refined_excited_freq_sweep_data, color='black')
plt.plot(refined_excited_sweep_freqs/GHz, refined_excited_sweep_y_fit, color='red')
plt.xlim([min(refined_excited_sweep_freqs/GHz), max(refined_excited_sweep_freqs/GHz)])
plt.xlabel("Frequency [GHz]", fontsize=15)
plt.ylabel("Measured Signal [a.u.]", fontsize=15)
plt.title("1->2 Frequency Sweep (refined pass)", fontsize=15)
plt.show()

_, qubit_12_freq, _, _ = refined_excited_sweep_fit_params
print(f"Our updated estimate for the 1->2 transition frequency is "
      f"{round(qubit_12_freq/GHz, 7)} GHz.")

Our updated estimate for the 1->2 transition frequency is 4.6263002 GHz.

2B. 1->2 ラビ実験

これで、$1\rightarrow2$周波数の良い推定が得られたので、$1\rightarrow2$遷移のための$\pi$パルス振幅を得るためのラビ実験を行います。そのために、$0\rightarrow1$ $\pi$ パルスを適用してから、$1\rightarrow2$周波数において駆動振幅をスイープします（サイドバンド法を使います）。

# 実験の構成
num_rabi_points = 75 # 実験数（つまり掃引する振幅）

# 駆動振幅の繰り返し値：0から1.0の間で均等に配置された75個の振幅
drive_amp_min = 0
drive_amp_max = 1.0
drive_amps = np.linspace(drive_amp_min, drive_amp_max, num_rabi_points)

# スケジュールの作成
rabi_12_schedules = []

# すべての駆動振幅をループします
for ii, drive_amp in enumerate(drive_amps):
    
    base_12_pulse = pulse_lib.gaussian(duration=drive_samples,
                                       sigma=drive_sigma,
                                       amp=drive_amp,
                                       name='base_12_pulse')
    # 1->2の周波数においてサイドバンドを適用
    rabi_12_pulse = apply_sideband(base_12_pulse, qubit_12_freq)
    
    # スケジュールにコマンドを追加
    schedule = pulse.Schedule(name='Rabi Experiment at drive amp = %s' % drive_amp)
    schedule |= pulse.Play(pi_pulse_01, drive_chan) # 0->1
    schedule |= pulse.Play(rabi_12_pulse, drive_chan) << schedule.duration # 1->2のラビパルス
    schedule |= measure << schedule.duration # 駆動パルスの後に測定をシフト
    
    rabi_12_schedules.append(schedule)

# プログラムにスケジュールを組み込みます
# 注：LO周波数はcal_qubit_freqであり、0->1のpiパルスを作ります；
# サイドバンドを使って、1->2のパルス用に変更されます
rabi_12_expt_program = assemble(rabi_12_schedules,
                                backend=backend,
                                meas_level=1,
                                meas_return='avg',
                                shots=NUM_SHOTS,
                                schedule_los=[{drive_chan: cal_qubit_freq}]
                                               * num_rabi_points)

rabi_12_job = backend.run(rabi_12_expt_program)

print(rabi_12_job.job_id())
job_monitor(rabi_12_job)

5f96940d00cb39001375c22e
Job Status: job has successfully run

# ジョブデータ（平均）を取得します
rabi_12_data = get_job_data(rabi_12_job, average=True)

We plot and fit our data as before.

# 注：信号の実部のみプロットします。
rabi_12_data = np.real(baseline_remove(rabi_12_data))
(rabi_12_fit_params, 
 rabi_12_y_fit) = fit_function(drive_amps,
                            rabi_12_data, 
                            lambda x, A, B, drive_12_period, phi: (A*np.cos(2*np.pi*x/drive_12_period - phi) + B),
                            [3, 0.5, 0.9, 0])

plt.scatter(drive_amps, rabi_12_data, color='black')
plt.plot(drive_amps, rabi_12_y_fit, color='red')

drive_12_period = rabi_12_fit_params[2]
# piパルス用の振幅のためにphiを考慮します
pi_amp_12 = (drive_12_period/2/np.pi) *(np.pi+rabi_12_fit_params[3])

plt.axvline(pi_amp_12, color='red', linestyle='--')
plt.axvline(pi_amp_12+drive_12_period/2, color='red', linestyle='--')
plt.annotate("", xy=(pi_amp_12+drive_12_period/2, 0), xytext=(pi_amp_12,0), arrowprops=dict(arrowstyle="<->", color='red'))
plt.annotate("$\pi$", xy=(pi_amp_12-0.03, 0.1), color='red')

plt.xlabel("Drive amp [a.u.]", fontsize=15)
plt.ylabel("Measured signal [a.u.]", fontsize=15)
plt.title('Rabi Experiment (1->2)', fontsize=20)
plt.show()

print(f"Our updated estimate for the 1->2 transition frequency is "
      f"{round(qubit_12_freq/GHz, 7)} GHz.")
print(f"Pi Amplitude (1->2) = {pi_amp_12}")

Our updated estimate for the 1->2 transition frequency is 4.6263002 GHz.
Pi Amplitude (1->2) = 0.37256049920143336

この情報を使って、$1\rightarrow2$ $\pi$パルスを定義できます。（必ず、$1\rightarrow2$周波数でサイドバンドを追加してください。）

pi_pulse_12 = pulse_lib.gaussian(duration=drive_samples,
                                 amp=pi_amp_12, 
                                 sigma=drive_sigma,
                                 name='pi_pulse_12')
# このパルスがサイドバンドであることを再確認してください
pi_pulse_12 = apply_sideband(pi_pulse_12, qubit_12_freq)

2C. 0,1,2の識別器を構築する

とうとう、 $|0\rangle$と$|1\rangle$と$|2\rangle$状態の識別器を構築できます。手順はセクション1と同様ですが、$|2\rangle$状態のためにスケジュールを追加します。

3つのスケジュールがあります。（再度、私たちのシステムが$|0\rangle$から開始することを思い出してください）：

1. $|0\rangle$状態を直接測定します。（$|0\rangle$の重心を得ます。）

2. $0\rightarrow1$ $\pi$パルスを適用し、測定します。（$|1\rangle$の重心を得ます。）

3. $0\rightarrow1$ $\pi$パルスを適用した後、$1\rightarrow2$ $\pi$パルスを適用しそして測定します。（$|2\rangle$の重心を得ます。）

# 3つのスケジュールを作ります

# 基底状態のスケジュール
zero_schedule = pulse.Schedule(name="zero schedule")
zero_schedule |= measure

# 励起状態のスケジュール
one_schedule = pulse.Schedule(name="one schedule")
one_schedule |= pulse.Play(pi_pulse_01, drive_chan)
one_schedule |= measure << one_schedule.duration

# 励起状態のスケジュール
two_schedule = pulse.Schedule(name="two schedule")
two_schedule |= pulse.Play(pi_pulse_01, drive_chan)
two_schedule |= pulse.Play(pi_pulse_12, drive_chan) << two_schedule.duration
two_schedule |= measure << two_schedule.duration

プログラムを構築し、IQ平面上に重心をプロットします。

# プログラムにスケジュールを組み込みます
IQ_012_program = assemble([zero_schedule, one_schedule, two_schedule],
                           backend=backend,
                           meas_level=1,
                           meas_return='single',
                           shots=NUM_SHOTS,
                           schedule_los=[{drive_chan: cal_qubit_freq}] * 3)

IQ_012_job = backend.run(IQ_012_program)

print(IQ_012_job.job_id())
job_monitor(IQ_012_job)

5f969656d3b8890012128a9e
Job Status: job has successfully run

# （単一の）ジョブデータを取得します；0,1,2に分割します
IQ_012_data = get_job_data(IQ_012_job, average=False)
zero_data = IQ_012_data[0]
one_data = IQ_012_data[1]
two_data = IQ_012_data[2]

def IQ_012_plot(x_min, x_max, y_min, y_max):
    """0、1、2のIQ平面をプロットするための補助関数。引数としてプロットの制限を与えます。
    """
    # 0のデータは青でプロット
    plt.scatter(np.real(zero_data), np.imag(zero_data), 
                    s=5, cmap='viridis', c='blue', alpha=0.5, label=r'$|0\rangle$')
    # 1のデータは赤でプロット
    plt.scatter(np.real(one_data), np.imag(one_data), 
                    s=5, cmap='viridis', c='red', alpha=0.5, label=r'$|1\rangle$')
    # 2のデータは緑でプロット
    plt.scatter(np.real(two_data), np.imag(two_data), 
                    s=5, cmap='viridis', c='green', alpha=0.5, label=r'$|2\rangle$')

    # 0、1、2の状態の結果の平均を大きなドットでプロット
    mean_zero = np.mean(zero_data) # 実部と虚部それぞれの平均をとる
    mean_one = np.mean(one_data)
    mean_two = np.mean(two_data)
    plt.scatter(np.real(mean_zero), np.imag(mean_zero), 
                s=200, cmap='viridis', c='black',alpha=1.0)
    plt.scatter(np.real(mean_one), np.imag(mean_one), 
                s=200, cmap='viridis', c='black',alpha=1.0)
    plt.scatter(np.real(mean_two), np.imag(mean_two), 
                s=200, cmap='viridis', c='black',alpha=1.0)
    
    plt.xlim(x_min, x_max)
    plt.ylim(y_min,y_max)
    plt.legend()
    plt.ylabel('I [a.u.]', fontsize=15)
    plt.xlabel('Q [a.u.]', fontsize=15)
    plt.title("0-1-2 discrimination", fontsize=15)

x_min = -10
x_max = 20
y_min = -25
y_max = 10
IQ_012_plot(x_min, x_max, y_min, y_max)

今回は、$|2\rangle$状態に対応した3個目の重心が観測されます。（注：プロットが見えない場合は、Notebookを再実行してください。）

このデータで、識別器を構築します。再びscikit.learn を使って線形判別分析(LDA)を使います。

LDAのためにデータを形成することから始めます。

# IQベクトルを作成します（実部と虚部で構成されています）
zero_data_reshaped = reshape_complex_vec(zero_data)
one_data_reshaped = reshape_complex_vec(one_data)  
two_data_reshaped = reshape_complex_vec(two_data)  

IQ_012_data = np.concatenate((zero_data_reshaped, one_data_reshaped, two_data_reshaped))
print(IQ_012_data.shape) # IQデータの形を確認します

(3072, 2)

次に、学習用データとテスト用データを分割します（前回と同じように半分ずつです）。テスト用データは、0スケジュールの場合、0の配列が含まれたベクトルで、1スケジュールの場合、1の配列が含まれたベクトルで、2スケジュールの場合2の配列が含まれたベクトルです。

# （テスト用に）0と1と2の値が含まれたベクトルを構築します
state_012 = np.zeros(NUM_SHOTS) # 実験のショット数
state_012 = np.concatenate((state_012, np.ones(NUM_SHOTS)))
state_012 = np.concatenate((state_012, 2*np.ones(NUM_SHOTS)))
print(len(state_012))

# データをシャッフルして学習用セットとテスト用セットに分割します
IQ_012_train, IQ_012_test, state_012_train, state_012_test = train_test_split(IQ_012_data, state_012, test_size=0.5)

3072

最後に、モデルを設定して学習します。学習の精度が出力されます。

# LDAを設定します
LDA_012 = LinearDiscriminantAnalysis()
LDA_012.fit(IQ_012_train, state_012_train)

LinearDiscriminantAnalysis()

# シンプルなデータでテストします
print(LDA_012.predict([[0, 0], [-10, 0], [-15, -5]]))

[0. 1. 1.]

# 精度を計算します
score_012 = LDA_012.score(IQ_012_test, state_012_test)
print(score_012)

0.826171875

最後のステップは、セパラトリックスのプロットです。

IQ_012_plot(x_min, x_max, y_min, y_max)
separatrixPlot(LDA_012, x_min, x_max, y_min, y_max, NUM_SHOTS)

3つの重心を得たので、セラパトリックスは線ではなく、2つの線の組み合わせを含む曲線になります。$|0\rangle$、$|1\rangle$と$|2\rangle$の状態を区別するために、私たちのモデルは、IQ上の点がセラパトリックスのどの側にあるかどこにあるかチェックし、それに応じて分類します。

3. 参考文献

1. D. C. McKay, T. Alexander, L. Bello, M. J. Biercuk, L. Bishop, J. Chen, J. M. Chow, A. D. C ́orcoles, D. Egger, S. Filipp, J. Gomez, M. Hush, A. Javadi-Abhari, D. Moreda, P. Nation, B. Paulovicks, E. Winston, C. J. Wood, J. Wootton, and J. M. Gambetta, “Qiskit backend specifications for OpenQASM and OpenPulse experiments,” 2018, https://arxiv.org/abs/1809.03452.

2. Krantz, P. et al. “A Quantum Engineer’s Guide to Superconducting Qubits.” Applied Physics Reviews 6.2 (2019): 021318, https://arxiv.org/abs/1904.06560.

3. Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011, https://scikit-learn.org/stable/modules/lda_qda.html#id4.

import qiskit.tools.jupyter
%qiskit_version_table