ベスパリブ

プログラミングを主とした日記・備忘録です。ベスパ持ってないです。

Qiskit Global Summer School 2022 参加記

Qiskit Global Summer School 2022

Qiskit Global Summer School 2022 に参加した記録です。

EDT時間で 7/18 ~ 7/29 の2週間、毎日(土日除く)動画がアップされるのでそれを観ながら勉強しつつ、演習のテキストブックもこなすというもの。

相談等はDiscordでできます。日本人用のチャンネルもあり。

Exercise 1

Qiskitの基本操作と、計算複雑性(complexity)の話。

量子回路の計算量はゲートの数と量子回路の深さが関係している。量子回路の深さはテトリスのように計算できる

lab1_ex1 ~ lab1_ex6

特に言うことなし。Qiskitの量子ビットは右側が下位ビットであることに注意する。

lab1_ex7

テキストに qft_rotations が与えられているので、それの回路数を数える問題。

下図の赤枠内のゲートを数えればよい。

赤枠の量子ゲートの個数を数えればよい

Pゲートの個数は等差数列の和で求めた。

def lab1_ex7(n:int) -> int:
    #Here we want you to build a function calculating the number of gates needed by the fourier rotation for n qubits.
    
    numberOfGates=0 
    #
    # FILL YOUR CODE IN HERE
    # TRY TO FIND AN EXPLIZIT (NON RECURSIVE) FUNCTION
    #
    # 等差数列の和の公式を使う
    numberOfGates = (n-1)*n//2 + n
    return numberOfGates

lab1_ex8

量子回路の深さは浅いほどよい。  |00 \dots 0⟩ + |11 \dots 1⟩を深さ5で作るには、どうすればよいかという問題。

  • H(0)
  • CX(0,1)
  • CX(0,2), CX(1,3)
  • CX(0,4), CX(1,5), CX(2,6), CX(3,7)
  • ...

という風にやれば深さ5で作れる。

def lab1_ex8():
    qc = QuantumCircuit(16) #Same as above
    
    #Step 1: Preparing the first qubit in superposition
    qc.h(0)
    #
    #
    # FILL YOUR CODE IN HERE
    #
    N = 16
    c_end = 1
    is_end = False
    while 1:
        if is_end:
            break

        for i in range(0, c_end):
            if i+c_end == N:
                is_end = True
                break
            qc.cx(i, i+c_end)
        c_end *= 2
    return qc

Exercise 2

Qiskit Opflow ライブラリと、ハミルトニアンの話。

詰まりどころがあったりするので、このあたりからDiscordのピン留めされたメッセージは読まないといけないと思った。

grade_lab2_ex1 ~ grade_lab2_ex7

あんまり言うことなし。ライブラリの使い方の問題。

grade_lab2_ex8

# Define a parameter t for the time in the time evolution operator
t = Parameter('t')

# Follow the instructions above to define a time-evolution operator
###INSERT CODE BELOW THIS LINE

# time_evolution_operator = (H).exp_i() # Wrong answar but pass.
time_evolution_operator = (H*t).exp_i()

###DO NOT EDIT BELOW THIS LINE

print(time_evolution_operator)

初期の頃は t がなくてもパスした。

grade_lab2_ex9 ~ grade_lab2_ex10

これもライブラリの使い方の問題。MatrixEvolution を使ってハミルトニアンの時間発展をしてみようという問題。

grade_lab2_ex11

# Define the number of steps needed to reach the previously set total time-evolution
steps = int(evolution_time/time_step_value)

# Compose the operator for a Trotter step several times to generate the 
# operator for the full time-evolution
###INSERT CODE BELOW THIS LINE

total_time_evolution_operator = I
for i in range(steps):
    total_time_evolution_operator @= time_evolution_operator
pt_evo = PauliTrotterEvolution()
total_time_evolution_circuit = pt_evo.convert(total_time_evolution_operator)
total_time_evolution_circuit = total_time_evolution_circuit.bind_parameters({t: time_step_value})
    
###DO NOT EDIT BELOW THIS LINE

print(total_time_evolution_circuit)

MatrixEvolution の代わりに PauliTrotterEvolution を使ってハミルトニアンの時間発展してみようという問題。時間発展する際に e^{−iHt}を作ることになるが、MatrixEvolution は正確に指数化するが、PauliTrotterEvolution はTrotter-Suzuki 分解を使って近似的に指数化するらしい。

grade_lab2_ex12

# Set number of qubits
num_qubits = 3
# Define time parameter
t = Parameter('t')
# Set total evolution time
evolution_time_t = 2
# Set size of time-step for Trotter evolution
time_step_value_t = 0.1
# Define the number of steps
steps_t = int(evolution_time_t/time_step_value_t)
# Create circuit
tight_binding_circuit = QuantumCircuit(num_qubits)
# Add initial state preparation
tight_binding_circuit.x(0)

# Define the Hamiltonian, the time-evolution operator, the Trotter step and the total evolution
###INSERT CODE BELOW THIS LINE

# ハミルトニアン作成
H = (I^X^X) + (X^X^I) + (I^Y^Y) + (Y^Y^I)
# H = (X^X^I) + (I^X^X) + (I^Y^Y) + (Y^Y^I) # トロッター化は近似的な手法なので、ハミルトニアンの項の順序が異なると結果も若干違ってくる

# 時間発展オペレーター作成
time_evo_op = (H*t).exp_i()
full_time_evo_op = I
for i in range(steps_t):
    full_time_evo_op @= time_evo_op

# PauliTrotterEvolutionを使って全体の時間発展
pt_evo = PauliTrotterEvolution()
exp_circuit = pt_evo.convert(full_time_evo_op)
bind_exp_circuit = exp_circuit.bind_parameters({t: time_step_value_t})

# 回路を合体
full_time_evolution_circuit = tight_binding_circuit.compose(bind_exp_circuit.to_circuit()) 

###DO NOT EDIT BELOW THIS LINE

print(full_time_evolution_circuit)

初期状態を |001⟩の状態で、で時間発展させるというもの。

最初、ハミルトニアンの項の順番が異なっていてパスできず詰まってしまった。Discordのピン留めされたメッセージに注意書きは確認しようと心に決めた。

Exercise 3

量子ノイズの話。

シミュレータではノイズのない完璧な回転操作をシミュレートしてくれるが、実機ではそうはいかない。シミュレータでもノイズを載せてシミュレートできるので、それをしてみようという内容。

ノイズは from qiskit.providers.aer.noise import NoiseModelNoiseModel ライブラリで作成できる。

grade_lab3_ex1 ~ grade_lab3_ex3

特に言うことなし。

ところで、回路のシミュレータの実行方法が2種類(実際はもっとある)あって、qiskit.execute を呼ぶ方法と、backend.run を呼ぶ方法があるのだけど、何が違うのかわかってないし、同じだとしたらなぜこのような実装になっているのかもわからない。

# 以下の2つは何が違う?
result1 = qiskit.execute(circuit, shots_backend, shots=num_shots).result()
result2 = shots_backend.run(circuit, shots=num_shots).result()

grade_lab3_ex4

scipycurve_fit を使って、カーブフィッティング(曲線あてはめ)をしてみようという問題。

# Define a Gaussian function for the fit
def gaussian(x, a, mean, sigma):
    return a * (1/(sigma*np.sqrt(2*np.pi)))*np.exp(-((x - mean)**2 / (2 * sigma**2)))

# Set initial estimates for the fit: the maximum height of the histogram, the theoretical
# average and the theoretical standard deviation
height_estimate = np.max(yh)
mu_estimate = 0.5
sigma_estimate = np.sqrt(mu_estimate*(1-mu_estimate)/num_shots)

# Define list to store the estimated values of the average (mu_sigma[0]) and standard deviation (mu_sigma[1])
mu_sigma = []

# Use the curve_fit function to fit the experimental data using the definition of the Gaussian function
# which will return the estimate of the parameters.
###INSERT CODE BELOW THIS LINE

popt, pcov = curve_fit(gaussian, x_01_h, yh, [height_estimate, mu_estimate, sigma_estimate])  # popt:推定値, pcov: 共分散

mu_sigma = [0]*2
mu_sigma[0] = mu_estimate
mu_sigma[1] = sigma_estimate

###DO NOT EDIT BELOW THIS LINE   

print("The mean is ", mu_sigma[0])
print("The standard deviation is ", mu_sigma[1])

# Plot experimental data and the fit
plt.scatter(x_01_h, yh, label = 'data', color = 'red')
plt.plot(x_01, gaussian(x_01, *popt), label = 'fit', linestyle = '--')
plt.title(f"Using {num_shots} shots to estimate probability")
plt.xlim((0, 1))
plt.xlabel(r'probability of $\vert 0 \rangle$')
plt.ylabel('counts')
plt.legend()
plt.show()

ガウス関数のカーブフィッティング

grade_lab3_ex5 ~ grade_lab3_ex7

特に言うことなし。

grade_lab3_ex8

測定ノイズ(Measurement Noise)の話。

ノイズのあるシミュレータにおいて、 |\psi⟩状態のものを測定して、実際に得られた状態が |\phi⟩だったときの確率を p(\psi|\phi)とする(逆じゃない?と思ったが、ノイズモデルの構築#読み出しエラー - Qiskitを参照すると合ってる。ただ、測定エラーの軽減では行列の形が転置している)と、 |0⟩状態と |1⟩状態を基底としたときの混同行列(confusion matrix)A と表せる。

測定して  |0⟩ を得られる確率を p0 |1⟩ を得られる確率を p1として、

  •  p_{ideal} = [p0_{ideal}, p1_{ideal} ] を理想的な確率ベクトル
  •  p_{noisy} = [p0_{noisy}, p1_{noisy} ]をノイズのある確率ベクトル

とすると、 と表せる。

# Define the confusion matrix from the probabilities found above
###INSERT CODE BELOW THIS LINE

confusion_matrix = np.array([np.array([0.0 for _ in range(2)]) for __ in range(2)])
confusion_matrix[0][0] = p0_0
confusion_matrix[0][1] = p0_1
confusion_matrix[1][0] = p1_0
confusion_matrix[1][1] = p1_1

###DO NOT EDIT BELOW THIS LINE


print("Confusion matrix:")
print(confusion_matrix)

grade_lab3_ex9

 p_{noisy} から  p_{ideal} を求めることもできる。

これにより、ノイズのある環境で得られた確率から、逆行列 A^{-1}を使って「緩和」して、理想的な確率を求めることができる。

# Invert the confusion matrix
inverse_confusion_matrix = np.linalg.inv(confusion_matrix)

# Mitigate the counts using the inverse of the confusion matrix
###INSERT CODE BELOW THIS LINE

p_vector_mitigated = inverse_confusion_matrix.dot(p_vector_noisy)
        
###DO NOT EDIT BELOW THIS LINE


print("Mitigated probability of |0>: ", p_vector_mitigated[0])
print("Mitigated probability of |1>: ", p_vector_mitigated[1])

grade_lab3_ex10

コヒーレントノイズ(Coherent Noise)の話。

コヒーレンス時間という言葉があって、それは「重ね合わせ状態を維持できる時間」という意味で理解しているのだけど、ここでコヒーレントノイズはゲートによる回転操作などの干渉によって発生する「干渉ノイズ」という意味っぽい。

ノイズモデルの作り方はやってくれているので、問題としてやること自体は難しくない。

# Set value of rotation
theta = np.pi
# Instantiate a quantum circuit
circuit = qiskit.QuantumCircuit(quantum_register, classical_register)

# Add a parametrized RX rotation and bind the value of the parameter. (By default, parameter binding is not an in-place operation)
# Then measure the qubit and calculate probability of seeing |0> after rx(np.pi) 
###INSERT CODE BELOW THIS LINE

circuit.rx(theta, 0)
circuit.measure(0,0)

###DO NOT EDIT BELOW THIS LINE

grade_lab3_ex11

インコヒーレントノイズ(Incoherent Noise)の話。

「非干渉ノイズ」とでも訳すのだろうか。何も操作をしない id ゲート(idle操作)でも、ノイズの影響を受けるらしい。

ちなみにデコヒーレンス(Decoherence)という言葉もあって、これは「量子重ね合わせ状態の破壊」という意味。

問題としては、「id ゲートにおける5%の脱分極エラー(depolarizing error)のモデルを作れ」というもの。

depolarizing って何?という話なのだが、量子技術序論 - 東京大学 のP.143, 11.5節を流し読むと、お気持ちは理解できる。曰く、

  • pの確率でX軸まわりの反転が起きるエラーを「ビットフリップエラー」
  • pの確率でZ軸まわりの反転が起きるエラーを「位相フリップエラー」
  • pの確率でY軸まわりの反転が起きるエラーを「ビット-位相フリップエラー」

という。脱分極エラー(depolarizing error)とは、確率pで「ビットフリップエラー」「位相フリップエラー」「ビット-位相フリップエラー」のいずれかが起きるエラーのこと。だと思う。

Qiskitのテキストブックだと量子ボリュームの測定で出てくるが、説明は詳しくない。

Qiskitライブラリとしてはqiskit.providers.aer.noise.depolarizing_errorで提供されている。

# Define number of shots
num_shots_inchoherent = 10000

# Create an empty noise model
depolarizing_noise_model = NoiseModel()

# Define a depolarizing error on the identity gate for qubit zero and add it to the noise model
###INSERT CODE BELOW THIS LINE

error = depolarizing_error(0.05, 1)
depolarizing_noise_model.add_quantum_error(error, ["id"], qubits=[0])

###DO NOT EDIT BELOW THIS LINE

Exercise 4

最後のExercise。Spin-1/2 Model という量子系をシミュレートする。

grade_lab4_ex1

Exercise 2 でやったopflowライブラリを使って XXX Heisenberg ハミルトニアンを作成する。

# Returns the XXX Heisenberg model for 3 spin-1/2 particles in a line
def ex1_compute_H_heis3():
    # FILL YOUR CODE IN HERE
    
    H = (I^X^X)+(X^X^I)+(I^Y^Y)+(Y^Y^I)+(I^Z^Z)+(Z^Z^I)
    
    # Return Hamiltonian
    return H

grade_lab4_ex2

ハミルトニアンがわかれば、その系の量子状態が時間とともにどのように変化するかはシュレディンガーの方程式を使って解くことができる。

シュレディンガー方程式

opflowを使って、問題文に言われたとおりに実装するだけ。

# Returns the time evolution operator U_heis3(t) for a given time t assuming an XXX Heisenberg Hamiltonian for 3 spins-1/2 particles in a line
def ex2_compute_U_heis3(t):
    # FILL YOUR CODE IN HERE

    H = ex1_compute_H_heis3()
    U = (H*t).exp_i()
    
    return U

grade_lab4_ex3

時間発展演算子 Uができたので、量子系の任意の状態がどのように変化するかをシミュレーションできるようになった。XXX Heisenberg spin model を古典的にシミュレートしてとりあえず終了。

fidelity(フィデリティ・忠実度) という単語が出てくるが、量子情報理論の基本:フィデリティ - Qiitaによると、2つの純粋状態の類似度のことらしい。計算は簡単で、内積の絶対値をとるだけ。 |⟨\psi|\phi⟩|

混合状態のフィデリティも求めることができるが、コツがあるっぽい。

state tomography(状態トモグラフィー)という単語が出てくるが、量子情報理論の基本:量子状態トモグラフィ - Qiitaを読んでなるほどした。

(問題文訳)状態トモグラフィーとは、量子回路を複数回実行し、異なる基底で測定することで、量子回路の終端での正確な量子状態(位相情報まで)を決定する方法。最終的な量子状態が期待される状態に近ければ近いほど、フィデリティは高くなる。

ここにきてトロッターの近似について理解が進む。1-6の説明文にと書いてあるが、これはHが行列で非可換なので一般に e^{H_1 + H_2} = e^{H_1} \cdot e^{H_2} は成り立たないから。

なので、という近似を使う。Nは大きくするほど近似の精度は良くなる。これをリー・トロッターの積公式という。ということがシュレーディンガー方程式の時間発展を量子回路で表現 - Qiitaに書いてあった。

grade_lab4_ex3の問題としては、YY ゲートを作成すること。これは1-6の解説にある論文(F. Tacchino, et al., Quantum Computers as Universal Quantum Simulators: State-of-the-Art and Perspectives, Adv. Quantum Technol. 3 3 (2020) free arXiv version)に書いてある。

YY(δ)の量子回路を作れば良い。(論文から引用)

# Build a subcircuit for YY(t) two-qubit gate
def ex3_compute_YY_gate(t):
    # FILL YOUR CODE IN HERE
    
    YY_qr = QuantumRegister(2)
    YY_qc = QuantumCircuit(YY_qr, name="YY")
    
    YY_qc.rx(np.pi/2,[0,1])
    YY_qc.cnot(0,1)
    YY_qc.rz(2 * t, 1)
    YY_qc.cnot(0,1)
    YY_qc.rx(-np.pi/2,[0,1])
    
    YY = YY_qc.to_instruction()

    return YY

grade_lab4_ex4

ZZ ゲートの作成をすればOK。これも論文に書いてある。

ZZ ゲートの式(論文から引用)

量子回路ではこう。(論文から引用)

# Build a subcircuit for ZZ(t) two-qubit gate
def ex4_compute_ZZ_gate(t):
    # FILL YOUR CODE IN HERE

    ZZ_qr = QuantumRegister(2)
    ZZ_qc = QuantumCircuit(ZZ_qr, name="ZZ")
    
    ZZ_qc.cnot(0,1)
    ZZ_qc.rz(2 * t, 1)
    ZZ_qc.cnot(0,1)
        
    ZZ = ZZ_qc.to_instruction()

    return ZZ

lab4-optimization-contest

ノイズのないシミュレータで実行したとき、フィデリティは0.0003 ± 0.0002だった。

シミュレータで実行したときのstate tomography fidelity

これを実機のibmq_manilaというproviderで実行するとフィデリティは 0.0441 ± 0.0017 だった。

実機で実行したときのstate tomography fidelity

lab4-optimization-contestでは、このフィデリティをなるべく大きくする追加の最適化コンテストである(上記の例でシミュレータよりも実機でやったほうがフィデリティが高いのはなぜ?)。

フィデリティはstate_fidelityライブラリで計算される。使い方は量子操作の概要#任意の初期化 - Qiskitに書いてある。

import math
desired_vector = [
    1 / math.sqrt(16) * complex(0, 1),
    1 / math.sqrt(8) * complex(1, 0),
    1 / math.sqrt(16) * complex(1, 1),
    0,
    0,
    1 / math.sqrt(8) * complex(1, 2),
    1 / math.sqrt(16) * complex(1, 0),
    0]

state_fidelity(desired_vector, desired_vector)
# 0.9999999999999996

フィデリティを内積だと考えると直交している場合でフィデリティ0で最小になると思うので、平行にすれば値が大きくなるイメージだろうか。正規化されているので最大値は1っぽい。

近似精度を高くするためにざっと以下の方法が考えられる。

  • トロッターのNを大きくする
  • トロッター化のためのサブシステムの分割方法を変えてみる
  • トロッター化以外の方法を試す

Discordの議論を読んで思いつく方法

  • 回路の深さを減らす

トロッターのNを大きくする(trotter_stepsを増やす)

元々は trotter_steps = 4 だったので、これを大きくしてみることを試す。ノイズがないシミュレータの場合はtrotter_stepsを増やせば増やすほどフィデリティは100%に近づく(0.9くらい行った)。ただ実機の場合は回路が長くなるとノイズもそれだけ載るはずなので、ただ増やせばいいというわけではないはず。

  • trotter_steps = 100

trotter_steps=100の結果

Discordの議論を見てみてると、8か20がいいらしい。4が最悪で、1でも高い値がでるらしい(なんで?)。

  • trotter_steps = 8

trotter_steps=8の結果

たしかに100と8は大差ないスコアが出た。

トロッター化のためのサブシステムを変えてみる

論文の画像を見ると、(a) 6-CNOT decomposition (b) 3-CNOT decomposition (c) 3-U decomposition がある。今は(a)の方法で実装しているので、(b)と(c)を試す。

(c) の U_{xy} は以下で実装する。

U_xy(δ)(論文から引用)

(b) 3-CNOT decompositionの実装は以下になった。

###
# (b) 3-CNOT decomposition を実装
###
def trot_qc_3_cnot_decomposition(t, to_gate=True):
    qr = QuantumRegister(2)
    qc = QuantumCircuit(qr, name='3-CNOT decomposition')
        
    qc.cnot(0,1)
    qc.rx(2*t - math.pi/2, 0)
    qc.h(0)
    qc.rz(2*t, 1)
    qc.cnot(0,1)
    qc.h(0)
    qc.append(RZGate(2*t).inverse(), [1]) # qc.rzdg(2*t, 1)
    qc.cnot(0,1)
    qc.rx(math.pi/2, 0)
    qc.append(RXGate(math.pi/2).inverse(), [1]) # qc.rxdg(math.pi/2, 1)

    if not to_gate: return qc
    
    return qc.to_instruction()

Trot_qc_3_cnot_decom = trot_qc_3_cnot_decomposition(t, to_gate=False)
Trot_qc_3_cnot_decom.draw()

3-CNOT decomposition 回路

(c) 3-U decompositionの実装は以下のようになった。

###
# (c) 3-U decomposition を実装
###
def compute_U_xy(t):
    H = (X^X) + (Y^Y)
    U = (H*t).exp_i()
    return U

def U_xy(t):
    target_time = np.pi
    trotter_steps = 100

    U = compute_U_xy(t)
        
    pt_evo = PauliTrotterEvolution()
    exp_circuit = pt_evo.convert(U)
    
    bind_exp_circuit = exp_circuit.bind_parameters({t: target_time/trotter_steps})
    return bind_exp_circuit
# Uxy = U_xy(t).to_circuit()

def trot_qc_3_U(t, to_gate=True):
    qr = QuantumRegister(2)
    qc = QuantumCircuit(qr, name="3-U decomposition")
    
    qc = qc.compose(U_xy(t/2).to_circuit())
    qc.rx(math.pi/2, [0])
    qc.rx(math.pi/2, [1])
    qc = qc.compose(U_xy(t/2).to_circuit())
    qc.append(RXGate(math.pi/2).inverse(), [0])
    qc.append(RXGate(math.pi/2).inverse(), [1])
    qc.ry(math.pi/2, [0])
    qc.ry(math.pi/2, [1])
    qc = qc.compose(U_xy(t/2).to_circuit())
    qc.append(RYGate(math.pi/2).inverse(), [0])
    qc.append(RYGate(math.pi/2).inverse(), [1])
    
    return qc

qc_test = trot_qc_3_U(t)
qc_test.decompose().draw()

3-U decomposition 回路

ibmq_manilaを再現するシミュレータの作成

実機はibmq_manilaを使うことになっている。毎回実機にジョブを投げると時間がかかってしょうがないので、ローカルシミュレータにibmq_manilaと同程度のエラーを再現するノイズを作りたい。

ibmq_manilaのエラー率の調べ方は、 IBM Quantumに行き、左上のメニューから「Compute Resources」をクリックすると各ibmqの資源情報が見れる。そこの検索バーに「manila」と入力し、ibmq_manilaをクリックすると詳細が出てくる。「Table view」にそれっぽい情報があるので、これを参考にノイズを作ればいいと思う。

ibmq_manilaのスペック

たとえば、「Prob meas0 prep1」の項目が0.0426とは、「状態|1⟩で準備されたものが測定したら|0⟩だったときの確率」だと思う。

色々試行錯誤した結果、以下のようなノイズモデルを作成した。

def make_noise_model():
    # Measurement Noise
    noise_model = NoiseModel()
    # Measurement miss-assignement probabilities
    p_meas_prep_list = [[0.0426, 0.0138], [0.0304, 0.0256], [0.0288, 0.0164], [0.0636, 0.0182], [0.058, 0.0102]]
    
    for i in range(5):
        p0given1 = p_meas_prep_list[i][0] # Probability of measuuring |0> given the state is |1>
        p1given0 = p_meas_prep_list[i][1] # Probability of measuring |1> given the state is |0>
        readout_error = ReadoutError([[1 - p1given0, p1given0], [p0given1, 1 - p0given1]]) # Define readout error
        noise_model.add_readout_error(readout_error, qubits=[i]) # Add error to noise model
    
    # Pauli error
    pauli_error_list = [2.012e-4, 2.142e-4, 2.298e-4, 4.420e-4, 5.142e-4]
    for i in range(5):
        error = depolarizing_error(pauli_error_list[i], 1)
        noise_model.add_quantum_error(error, ["x", "y", "z"], qubits=[i])

    # U error
    u_error_list = [2.012e-4, 2.142e-4, 2.298e-4, 4.420e-4, 5.142e-4]
    for i in range(5):
        error = depolarizing_error(pauli_error_list[i], 1)
        noise_model.add_quantum_error(error, ["u1", "u2", "u3"], qubits=[i])

    # CNOT error (mixed_unitary_error使うべき?)
    cnot_error_list = [5.333e-3, 8.854e-3, 1.023e-2, 9.058e-3, 9.058e-3]
    for i in range(5):
        error = depolarizing_error(cnot_error_list[i], 1)
        noise_model.add_quantum_error(error, ["cnot"], qubits=[i])
    
    # Coherent Noise
    # Construct a 1 qubit over-rotation of the RX gate
    # epsilon = np.pi/200 # over rotation amount
    epsilon = np.pi/50 # over rotation amount
    epsilon_rotation = RXGate(epsilon).to_matrix() # get matrix representation
    over_rotation = coherent_unitary_error(epsilon_rotation)
    for i in range(5):
        noise_model.add_quantum_error(over_rotation, ['rx'], qubits=[i])
    epsilon_rotation = RYGate(epsilon).to_matrix() # get matrix representation
    over_rotation = coherent_unitary_error(epsilon_rotation)
    for i in range(5):
        noise_model.add_quantum_error(over_rotation, ['ry'], qubits=[i])
    epsilon_rotation = RZGate(epsilon).to_matrix() # get matrix representation
    over_rotation = coherent_unitary_error(epsilon_rotation)
    for i in range(5):
        noise_model.add_quantum_error(over_rotation, ['rz'], qubits=[i])
    
    # Incoherent Noise
    id_error_list = [2.012e-4, 2.142e-4, 2.298e-4, 4.420e-4, 5.142e-4]
    for i in range(5):
        error = depolarizing_error(id_error_list[i], 1)
        noise_model.add_quantum_error(error, ["id"], qubits=[i])
    
    return noise_model
noise_model_test = make_noise_model()
print(noise_model_test)
NoiseModel:
  Basis gates: ['cx', 'id', 'rx', 'ry', 'rz', 'sx', 'u1', 'u2', 'u3', 'x', 'y', 'z']
  Instructions with noise: ['rx', 'x', 'measure', 'cnot', 'rz', 'u2', 'y', 'id', 'z', 'ry', 'u3', 'u1']
  Qubits with noise: [0, 1, 2, 3, 4]
  Specific qubit errors: [('x', (0,)), ('x', (1,)), ('x', (2,)), ('x', (3,)), ('x', (4,)), ('y', (0,)), ('y', (1,)), ('y', (2,)), ('y', (3,)), ('y', (4,)), ('z', (0,)), ('z', (1,)), ('z', (2,)), ('z', (3,)), ('z', (4,)), ('u1', (0,)), ('u1', (1,)), ('u1', (2,)), ('u1', (3,)), ('u1', (4,)), ('u2', (0,)), ('u2', (1,)), ('u2', (2,)), ('u2', (3,)), ('u2', (4,)), ('u3', (0,)), ('u3', (1,)), ('u3', (2,)), ('u3', (3,)), ('u3', (4,)), ('cnot', (0,)), ('cnot', (1,)), ('cnot', (2,)), ('cnot', (3,)), ('cnot', (4,)), ('rx', (0,)), ('rx', (1,)), ('rx', (2,)), ('rx', (3,)), ('rx', (4,)), ('ry', (0,)), ('ry', (1,)), ('ry', (2,)), ('ry', (3,)), ('ry', (4,)), ('rz', (0,)), ('rz', (1,)), ('rz', (2,)), ('rz', (3,)), ('rz', (4,)), ('id', (0,)), ('id', (1,)), ('id', (2,)), ('id', (3,)), ('id', (4,)), ('measure', (0,)), ('measure', (1,)), ('measure', (2,)), ('measure', (3,)), ('measure', (4,))]

ノイズ有りシミュレータが作成できたので、各実装の回路をシミュレータと実機で実行してみて、ノイズの再現度を比べてみる。

(a) 元々のXXYYZZ の実装を使って、trotter_steps=8 でのノイズ有りシミュレータの実行結果と、実機での実行結果は以下のようになった。

XXYYZZ回路, trotter_steps=8 での実行結果

  • ノイズ有りシミュレータの結果がstate tomography fidelity on qasm_simulator = 0.3768 ± 0.0016
  • 実機の結果がstate tomography fidelity on qasm_simulator(ibmq_manila) = 0.3695 ± 0.0025

ノイズはいい感じに再現できた。

次に、(b) 3-CNOT decomposition回路を使って、trotter_steps=8でのノイズ有りシミュレータの実行結果と、実機での実行結果は以下のようになった。

3-CNOT回路, trotter_steps=8 での実行結果

  • ノイズ有りシミュレータの結果がstate tomography fidelity on qasm_simulator = 0.5231 ± 0.0015
  • 実機の結果がstate tomography fidelity on qasm_simulator(ibmq_manila) = 0.5284 ± 0.0030

次に、(c) 3-U decomposition回路を使って、trotter_steps=8でのノイズ有りシミュレータの実行結果と、実機での実行結果は以下のようになった。

3-U回路, trotter_steps=8 での実行結果

  • ノイズ有りシミュレータの結果がstate tomography fidelity on qasm_simulator = 0.4117 ± 0.0025
  • 実機の結果がstate tomography fidelity on qasm_simulator(ibmq_manila) = 0.2214 ± 0.0020

(c) の結果が低すぎる。なんか間違えてない?

C-NOTゲートの数を減らしたほうがいいのかなと思っていたので、どうしたもんかいのーといった感じ。

(b)の実装が一番スコアが高い結果となった。

Discordの議論を見る

ここでDiscordの議論を見ると、強化学習を使った方法でやって高いスコアを出した人がいたほか、4 cnot decomposition と General Error Mitigation と呼ばれる緩和技術を使って0.89722のスコアを出してる人がいた。

私はconfusion matrixを作ったりしたが、それをResultに適用する方法がいまいちわからなかった。

無理やりResultの結果を書き換えることはできて、

# confution_matrixを使ってjobの結果を緩和する
def mitigate_job_result(job_result):
    mat = make_confusion_matrix()
    # mat = make_confusion_matrix2()
    inv_mat = np.linalg.inv(mat)

    mitigated_exp_results = []
    for result in job_result.results:
        # p_noisy作成
        p_noisy = {}
        for key, value in result.data.counts.items():
            p_noisy[key] = value/result.shots
        
        # "0x0", "0x1", "0x2",... の順にp_vector_noisyを入れていく
        p_vector_noisy = []
        key_list = list(p_noisy.keys())
        key_list.sort()
        for key in key_list:
            p_vector_noisy.append(p_noisy[key])

        # 緩和した確率を得る
        p_vector_mitigated = inv_mat.dot(p_vector_noisy)

        # 緩和した結果を格納する
        result_dict = result.to_dict()
        new_data = {}
        new_dict = {}
        for key, p in zip(key_list, p_vector_mitigated):
            new_dict[key] = int(result.shots*p)
        if sum(new_dict.values()) != result.shots:
            total = sum(new_dict.values())
            while total > result.shots:
                total -= 1
                new_dict["0x0"] -= 1
            while total < result.shots:
                total += 1
                new_dict["0x6"] += 1
        
        new_data["counts"] = new_dict
        result_dict["data"] = new_data
        mitigated_exp_results.append(result.from_dict(result_dict))

    # 緩和後の結果に書き換える
    for i in range(len(job_result.results)):
        job_result.results[i] = mitigated_exp_results[i]

    # print(job_result.results)
    return job_result

とやったが、フィデリティを計算する関数でエラーが出てしまった。

CompleteMeasFitterを使えば良いともあって試したが、このライブラリがよくわからなくて、チュートリアルどおりにやるとフィデリティが下がってしまった。

あれこれ試したが時間切れになってしまった。結局ベストスコアは3-CNOT回路でtrotter_steps=8 の 0.5284 ± 0.0030 の結果で終了した。

おわり

去年と同様に修了証明書としてバッジをもらえた。

www.credly.com