#!pip install qulacs
#!pip install qulacsvis
#!pip install matplotlib
#!pip install numpy
#!pip install scipy

import matplotlib.pyplot as plt
import numpy as np
from qulacs import Observable, QuantumState, PauliOperator, QuantumCircuit, ParametricQuantumCircuit, DensityMatrix
from qulacs.gate import  Identity,X, Y, H, Z, S, Sdag, DepolarizingNoise, DephasingNoise, AmplitudeDampingNoise, TwoQubitDepolarizingNoise

擬確率法を実装してみる

シミュレーション・エラーの入れ方は外挿法のときと共通。

エラーなしの場合

１次元横磁場イジング模型（４サイト）の基底状態をVQEで求める。解く対象は外挿法のものと共通。 \begin{align} H=\sum_{i}(JZ_iZ_{i+1}+hX_i) \end{align}

#１次元横磁場イジング模型のハミルトニアンを定義
n_qubits = 4
transverse_Ising_hamiltonian = Observable(n_qubits)
J = -1.0
h = -1.0
for i in range(n_qubits):
    transverse_Ising_hamiltonian.add_operator(J, f"Z {i} Z {(i+1)%n_qubits}")
    transverse_Ising_hamiltonian.add_operator(h, f"X {i}")

## まずはエラーなしの場合
from qulacsvis import circuit_drawer
from scipy.optimize import minimize

#状態ベクトル、ansatz回路を定義
n_qubits = 4
state = QuantumState(n_qubits)
circuit_depth = 4
circuit = ParametricQuantumCircuit(n_qubits)

#ansatz回路を指定
def Circuit(n_qubits,circuit_depth):
    circuit = ParametricQuantumCircuit(n_qubits)
    for d in range(circuit_depth):
        for i in range(n_qubits):
            circuit.add_parametric_RY_gate(i, 0.0)
        for i in range(n_qubits):
            circuit.add_CZ_gate(i, (i+1)%n_qubits)
    return circuit
circuit = Circuit(n_qubits,circuit_depth)
params = np.zeros(circuit.get_parameter_count())

#ansatz回路を表示
circuit_drawer(circuit, "mpl")

#コスト関数を定義。ここではエネルギー期待値がコスト関数。
def get_cost(params):
    state = QuantumState(n_qubits)
    for i, p in enumerate(params):
        circuit.set_parameter(i,p)
    circuit.update_quantum_state(state)
    return transverse_Ising_hamiltonian.get_expectation_value(state)

#上で定義したansatz回路のもとで、コスト関数を最小化する。
params_initial = np.zeros_like(params)
minimized_costfunction = minimize(get_cost, params_initial, method='BFGS')
exact_energy = transverse_Ising_hamiltonian.solve_ground_state_eigenvalue_by_power_method(state, 50)
#コスト関数の最小値を表示
print("VQEで求めたコスト関数の最小値", minimized_costfunction.fun)
print("基底状態のエネルギー", np.real(exact_energy))

VQEで求めたコスト関数の最小値 -5.072893485692576
基底状態のエネルギー -5.203541892294178

#エラーがない場合のansatz回路のパラメータを取得
params_true = minimized_costfunction.x

エラーありの場合

外挿法のときと同様に、各ゲートに以下のエラーを作用させます。

Depolarizing noise: \(\mathcal{E}[\rho]=(1-p)\rho+\frac{p}{3} (X\rho X+Y\rho Y+Z\rho Z)\)

## 次にdepolarizing noiseありの場合（１次元横磁場イジング模型）

#密度演算子、ansatz回路を定義
density_state = DensityMatrix(n_qubits)
noisy_circuit = ParametricQuantumCircuit(n_qubits)
error_rate = 0.01

#ansatz回路を指定
def NoisyCircuit(error_rate):
    noisy_circuit = ParametricQuantumCircuit(n_qubits)
    for d in range(circuit_depth):
        for i in range(n_qubits):
            noisy_circuit.add_parametric_RY_gate(i, 0.0)
            noisy_circuit.add_gate(DepolarizingNoise(i, error_rate))
        for i in range(n_qubits):
            noisy_circuit.add_CZ_gate(i, (i+1)%n_qubits)
            noisy_circuit.add_gate(DepolarizingNoise(i, error_rate))
            noisy_circuit.add_gate(DepolarizingNoise((i+1)%n_qubits, error_rate))
    return noisy_circuit

noisy_circuit = NoisyCircuit(error_rate)

#コスト関数を定義。ここではエネルギー期待値がコスト関数。
def noisy_get_cost(params):
    density_state = DensityMatrix(n_qubits)
    for i, p in enumerate(params):
        noisy_circuit.set_parameter(i,p)
    noisy_circuit.update_quantum_state(density_state)
    return transverse_Ising_hamiltonian.get_expectation_value(density_state)

#エラーなし、エラーありの場合で期待値を計算
minimized_costfunction_noiseless = get_cost(params_true)
minimized_costfunction_noise = noisy_get_cost(params_true)

#コスト関数の最小値を表示
print("VQEで求めたコスト関数の最小値(エラーなし)：", minimized_costfunction_noiseless)
print("VQEで求めたコスト関数の最小値(エラーあり)：",minimized_costfunction_noise)

VQEで求めたコスト関数の最小値(エラーなし)： -5.072893485692576
VQEで求めたコスト関数の最小値(エラーあり)： -3.8479582178671157

エラーの影響で期待値が変わっていることが確認できます。次にサンプリングにより、期待値を求めます。

def get_cost_sampling(shots,params,error_rate):
    density_state = DensityMatrix(n_qubits)
    noisy_circuit = NoisyCircuit(error_rate)
    for i, p in enumerate(params):
        noisy_circuit.set_parameter(i,p)
    noisy_circuit.update_quantum_state(density_state)
    density_copy_ZZ = DensityMatrix(n_qubits)
    density_copy_X = DensityMatrix(n_qubits)
    density_copy_ZZ = density_state.copy()
    density_copy_X = density_state.copy()
    mask_list_ZZ = [0b1100,0b0110,0b0011,0b1001]
    mask_list_X = [0b1000,0b0100,0b0010,0b0001]
    for i in range(n_qubits):
        gate = H(i)
        gate.update_quantum_state(density_copy_X)
    noise_sampling_ZZ = density_copy_ZZ.sampling(shots)
    noise_sampling_X = density_copy_X.sampling(shots)
    estimated_ZZ_average = 0.0
    estimated_X_average = 0.0
    for s in noise_sampling_ZZ:
        for mask in mask_list_ZZ:
            bitcount = bin(s & mask).count("1")
            estimated_ZZ_average += (-1)**bitcount/shots
    for s in noise_sampling_X:
        for mask in mask_list_X:
            bitcount = bin(s & mask).count("1")
            estimated_X_average += (-1)**bitcount/shots
    data_sampling_average = h * estimated_X_average + J * estimated_ZZ_average
    return data_sampling_average

#エラー率を固定して、ヒストグラムを作成
error_rate = 0.02
n_test = 200
shots = 1000
no_mitigation_list = np.zeros(n_test)
no_mitigation_noisy_list = np.zeros(n_test)

for i in range(n_test):
    density_state = DensityMatrix(n_qubits)
    noiseless_circuit = NoisyCircuit(0.0)
    no_mitigation_list[i] = get_cost_sampling(shots,params_true,0.0)
    noisy_circuit = NoisyCircuit(error_rate)
    no_mitigation_noisy_list[i] = get_cost_sampling(shots,params_true,error_rate)
expec_value = noisy_get_cost(params_true)

plt.hist(no_mitigation_list, alpha=0.2, bins=20, label=f"Noiseless")
plt.hist(no_mitigation_noisy_list, alpha=0.2, bins=20, label=f"Noisy")
plt.axvline(expec_value, color="black")
plt.axvline(get_cost(params_true), color="red")
plt.legend(loc="upper left", fontsize=13)
plt.show()

擬確率法

擬確率法を実装して、エラーを抑制する。

エラーの逆写像は以下の通り。

\(\mathcal{E}^{-1}[\rho]=\frac{3-p}{3-4p}\rho-\frac{p}{3-4p} (X\rho X+Y\rho Y+Z\rho Z)\)

この逆写像そのものは量子回路上で実装できないので、確率的に実装する。

#逆写像を確率的に実装するゲートを定義
#パリティの値も取得する。恒等演算子のときは＋１、それ以外だと－１にする。
import random
def recovery_gate_depolarizing(i,error_rate):
    ran = random.uniform(0, 1)
    q1 = (3.0-error_rate)/(3.0-4.0*error_rate)
    q2 = error_rate/(3.0-4.0*error_rate)
    p1 = q1/(q1+3.0*q2)
    p2 = (q1+q2)/(q1+3.0*q2)
    p3 = (q1+2.0*q2)/(q1+3.0*q2)
    if 0 < ran < p1:
        gate = Identity(i)
        parity = 1.0
    elif p1 < ran < p2:
        gate = X(i)
        parity = -1.0
    elif p2 < ran < p3:
        gate = Y(i)
        parity = -1.0
    elif p3 < ran < 1.0:
        gate = Z(i)
        parity = -1.0
    outcome = [gate,parity]
    return outcome

エラーあり量子回路に対して擬確率法を実行する手順は以下の通り。 - 上で定義した確率的なゲートを各depolarizing noiseの直後で実行 - 確率的に選ばれたゲートが恒等演算子の場合は＋１、それ意外（X,Y,Z）の場合は－１としてパリティの値を定義。すべてのゲートについてパリティの積を計算し、測定結果にかける - この操作を何度も繰り返し、測定結果の平均値を求める - 求めた平均値に対して、\(C^{N_g}\)をかける。\(C\)は各ゲートのコスト（\(=\frac{3.0+2.0p}{3.0-4.0p}\)）、\(N_g\)は実行するゲート数。

#上で定義した確率的なゲートを脱分極エラーの直後に挿入
def NoisyCircuit_recovery(error_rate):
    noisy_circuit_PEC = ParametricQuantumCircuit(n_qubits)
    parity_total = 1.0
    for d in range(circuit_depth):
        for i in range(n_qubits):
            noisy_circuit_PEC.add_parametric_RY_gate(i, 0.0)
            noisy_circuit_PEC.add_gate(DepolarizingNoise(i, error_rate))
            gate_parity = recovery_gate_depolarizing(i,error_rate)
            noisy_circuit_PEC.add_gate(gate_parity[0])
            #ここで、各ゲートのパリティの合計を計算する
            parity_total = parity_total*gate_parity[1]
        for i in range(n_qubits):
            noisy_circuit_PEC.add_CZ_gate(i, (i+1)%n_qubits)
            noisy_circuit_PEC.add_gate(DepolarizingNoise(i, error_rate))
            gate_parity = recovery_gate_depolarizing(i,error_rate)
            noisy_circuit_PEC.add_gate(gate_parity[0])
            parity_total = parity_total*gate_parity[1]
            noisy_circuit_PEC.add_gate(DepolarizingNoise((i+1)%n_qubits, error_rate))
            gate_parity = recovery_gate_depolarizing((i+1)%n_qubits,error_rate)
            noisy_circuit_PEC.add_gate(gate_parity[0])
            #ここで、各ゲートのパリティの合計を計算する
            parity_total = parity_total*gate_parity[1]
    outcome = [noisy_circuit_PEC,parity_total]
    return outcome

noisy_circuit_PEC = NoisyCircuit_recovery(error_rate)[0]
parity_tot = NoisyCircuit_recovery(error_rate)[1]

#上で定義した回路からサンプリングする
def get_cost_sampling_PEC(shots,params):
    estimated_ZZ_average = 0.0
    estimated_X_average = 0.0
    for j in range(shots):
        state = QuantumState(n_qubits)
        outcome = NoisyCircuit_recovery(error_rate)
        cost = (3.0+2.0* error_rate)/(3.0-4.0*error_rate)
        noisy_circuit_PEC = outcome[0]
        for i, p in enumerate(params):
            noisy_circuit_PEC.set_parameter(i,p)
        noisy_circuit_PEC.update_quantum_state(state)
        state_copy_ZZ = QuantumState(n_qubits)
        state_copy_X = QuantumState(n_qubits)
        state_copy_ZZ = state.copy()
        state_copy_X = state.copy()
        mask_list_ZZ = [0b1100,0b0110,0b0011,0b1001]
        mask_list_X = [0b1000,0b0100,0b0010,0b0001]
        for i in range(n_qubits):
            gate = H(i)
            gate.update_quantum_state(state_copy_X)
        noise_sampling_ZZ = state_copy_ZZ.sampling(1)
        noise_sampling_X = state_copy_X.sampling(1)
        for s in noise_sampling_ZZ:
            for mask in mask_list_ZZ:
                bitcount = bin(s & mask).count("1")
                #outcome[1]はパリティの値。パリティの値を作用させる。
                estimated_ZZ_average += (-1)**bitcount/shots * outcome[1]
        for s in noise_sampling_X:
            for mask in mask_list_X:
                bitcount = bin(s & mask).count("1")
                estimated_X_average += (-1)**bitcount/shots * outcome[1]
    #最後に測定結果を足し合わせる際に、全体にコストを掛ける。
    data_sampling_average = (h * estimated_X_average + J * estimated_ZZ_average) * cost**((3.0*n_qubits) * circuit_depth)
    return data_sampling_average

#擬確率法でエラーの影響を削減できていることを確認
shots = 10000

print("VQEで求めたコスト関数の最小値(エラーなし)：", minimized_costfunction_noiseless)
print("VQEで求めたコスト関数の最小値(エラーあり)：",minimized_costfunction_noise)
print("エラーを擬確率で削減した結果（統計誤差を含む）", get_cost_sampling_PEC(shots,params_true))

VQEで求めたコスト関数の最小値(エラーなし)： -5.072893485692576
VQEで求めたコスト関数の最小値(エラーあり)： -3.8479582178671157
エラーを擬確率で削減した結果（統計誤差を含む） -5.374169942495543

ヒストグラムを作成

（※計算時間が長い）

#エラー率を固定して、ヒストグラムを作成
error_rate = 0.02
n_test = 200
shots = 1000
PEC_list = np.zeros(n_test)

for i in range(n_test):
    state = QuantumState(n_qubits)
    noisy_circuit_PEC = ParametricQuantumCircuit(n_qubits)
    PEC_list[i] = get_cost_sampling_PEC(shots,params_true)

plt.hist(no_mitigation_noisy_list, alpha=0.2, bins=30, label=f"no mitigation")
plt.hist(PEC_list, alpha=0.2, bins=30, label=f"Mitigated")
plt.axvline(expec_value, color="black")
plt.axvline(get_cost(params_true), color="red")
plt.legend(loc="upper left", fontsize=13)
plt.show()

擬確率法により、確かにエラーが抑制できていることが分かります。バイアスが減ったぶん、分散が増えている様子が見えているので、サンプリングコストを計算してみます。

#得られたヒストグラムから、それぞれの標準偏差を求めてみる。標準偏差の比から測定コストが求まる。
std_no_mitigation = np.std(no_mitigation_noisy_list)
std_PEC = np.std(PEC_list)

cost_per_gate = (3.0+2.0* error_rate)/(3.0-4.0*error_rate)
total_cost = cost_per_gate**((3.0*n_qubits) * circuit_depth)

#計算して得られた標準偏差の増加分と、コストの見積もりの比較
print("計算して得られた標準偏差の増加分：",std_PEC/std_no_mitigation)
print("理論式から求めた測定コスト：",total_cost)

計算して得られた標準偏差の増加分： 9.703427076809103
理論式から求めた測定コスト： 6.911226777901118

測定コストの理論式は、少し小さめの値になっていることが分かりました。