<jupyter_start><jupyter_text>Experiment with adding noise originating from a IBMQ machine to the pennylane's builtin QNG optimizer and to a simple manual method of computing the Fisher information metric tensor.More info on the algorithm:https://pennylane.ai/qml/demos/tutorial_quantum_natural_gradient.htmlMore info on the noise model:https://qiskit.org/documentation/stubs/qiskit.providers.aer.noise.NoiseModel.htmlqiskit.providers.aer.noise.NoiseModelhttps://qiskit.org/documentation/stubs/qiskit.providers.aer.QasmSimulator.htmlqiskit.providers.aer.QasmSimulator<jupyter_code>import numpy as np
import qiskit
import pennylane as qml
from pennylane import expval, var

qiskit.IBMQ.load_account()

# Build noise model from backend properties
provider = qiskit.IBMQ.get_provider(group='open')
ibmq_backend = provider.get_backend('ibmq_burlington')
device_properties = ibmq_backend.properties()

noise_model = qiskit.providers.aer.noise.NoiseModel.from_backend(device_properties)

# Get coupling map from backend
coupling_map = ibmq_backend.configuration().coupling_map

# Get basis gates from noise model
basis_gates = noise_model.basis_gates

# Provision the the default device with noise
dev = qml.device('qiskit.aer', wires=3, noise_model=noise_model, 
                 basis_gates=basis_gates, coupling_map=coupling_map, backend='qasm_simulator')


@qml.qnode(dev)
def circuit(params):
    # |psi_0>: state preparation
    qml.RY(np.pi / 4, wires=0)
    qml.RY(np.pi / 3, wires=1)
    qml.RY(np.pi / 7, wires=2)

    # V0(theta0, theta1): Parametrized layer 0
    qml.RZ(params[0], wires=0) 
    qml.RZ(params[1], wires=1)

    # W1: non-parametrized gates
    qml.CNOT(wires=[0, 1])
    qml.CNOT(wires=[1, 2])

    # V_1(theta2, theta3): Parametrized layer 1
    qml.RY(params[2], wires=1)
    qml.RX(params[3], wires=2)

    # W2: non-parametrized gates
    qml.CNOT(wires=[0, 1])
    qml.CNOT(wires=[1, 2])

    return qml.expval(qml.PauliY(0))


params = np.array([0.432, -0.123, 0.543, 0.233])
g0 = np.zeros([2, 2])


def layer0_subcircuit(params):
    """This function contains all gates that
    precede parametrized layer 0"""
    qml.RY(np.pi / 4, wires=0)
    qml.RY(np.pi / 3, wires=1)
    qml.RY(np.pi / 7, wires=2)

@qml.qnode(dev)
def layer0_diag(params):
    layer0_subcircuit(params)
    return var(qml.PauliZ(0)), var(qml.PauliZ(1))


# calculate the diagonal terms
varK0, varK1 = layer0_diag(params)
g0[0, 0] = varK0 / 4
g0[1, 1] = varK1 / 4

@qml.qnode(dev)
def layer0_off_diag_single(params):
    layer0_subcircuit(params)
    return expval(qml.PauliZ(0)), expval(qml.PauliZ(1))


@qml.qnode(dev)
def layer0_off_diag_double(params):
    layer0_subcircuit(params)
    ZZ = np.kron(np.diag([1, -1]), np.diag([1, -1]))
    return expval(qml.Hermitian(ZZ, wires=[0, 1]))


# calculate the off-diagonal terms
exK0, exK1 = layer0_off_diag_single(params)
exK0K1 = layer0_off_diag_double(params)

g0[0, 1] = (exK0K1 - exK0 * exK1) / 4
g0[1, 0] = (exK0K1 - exK0 * exK1) / 4

#########################################
g1 = np.zeros([2, 2])


def layer1_subcircuit(params):
    """This function contains all gates that
    precede parametrized layer 1"""
    # |psi_0>: state preparation
    qml.RY(np.pi / 4, wires=0)
    qml.RY(np.pi / 3, wires=1)
    qml.RY(np.pi / 7, wires=2)

    # V0(theta0, theta1): Parametrized layer 0
    qml.RZ(params[0], wires=0)
    qml.RZ(params[1], wires=1)

    # W1: non-parametrized gates
    qml.CNOT(wires=[0, 1])
    qml.CNOT(wires=[1, 2])

@qml.qnode(dev)
def layer1_diag(params):
    layer1_subcircuit(params)
    return var(qml.PauliY(1)), var(qml.PauliX(2))

varK0, varK1 = layer1_diag(params)
g1[0, 0] = varK0 / 4
g1[1, 1] = varK1 / 4

@qml.qnode(dev)
def layer1_off_diag_single(params):
    layer1_subcircuit(params)
    return expval(qml.PauliY(1)), expval(qml.PauliX(2))


@qml.qnode(dev)
def layer1_off_diag_double(params):
    layer1_subcircuit(params)
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    YX = np.kron(Y, X)
    return expval(qml.Hermitian(YX, wires=[1, 2]))


# calculate the off-diagonal terms
exK0, exK1 = layer1_off_diag_single(params)
exK0K1 = layer1_off_diag_double(params)

g1[0, 1] = (exK0K1 - exK0 * exK1) / 4
g1[1, 0] = g1[0, 1]

from scipy.linalg import block_diag

g = block_diag(g0, g1)
print(np.round(g, 8))

#print(np.round(circuit.metric_tensor([params]), 8))

print(circuit.metric_tensor([params], diag_approx=True))
steps = 200
init_params = np.array([0.432, -0.123, 0.543, 0.233])

gd_cost = []
opt = qml.GradientDescentOptimizer(0.01)

print("Starting GD Optimizer run...")
theta = init_params
for _ in range(steps):
    theta = opt.step(circuit, theta)
    gd_cost.append(circuit(theta))
print("Done.")
    
qng_cost = []
opt = qml.QNGOptimizer(0.01, diag_approx=True)

print("Starting builtin QNG Optimizer run...")
theta = init_params
for _ in range(steps):
    theta = opt.step(circuit, theta)
    qng_cost.append(circuit(theta))
print("Done.")
    
from matplotlib import pyplot as plt

plt.style.use("seaborn")
plt.plot(gd_cost, "b", label="Vanilla gradient descent")
plt.plot(qng_cost, "g", label="Quantum natural gradient descent")

plt.ylabel("Cost function value")
plt.xlabel("Optimization steps")
plt.legend()
plt.show()<jupyter_output>Starting GD Optimizer run...
Done.
Starting builtin QNG Optimizer run...
Done.<jupyter_text>Now let's use measurement error mitigation.<jupyter_code>with pennylane_extra.qiskit_measurement_error_mitigation(shots=1024):
    gd_cost = []
    opt = qml.GradientDescentOptimizer(0.01)

    print("Starting GD Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit, theta)
        gd_cost.append(circuit(theta))
    print("Done.")
    
    qng_cost_mitigated = []
    opt = qml.QNGOptimizer(0.01, diag_approx=True)

    print("Starting builtin QNG Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit, theta)
        qng_cost_mitigated.append(circuit(theta))
    print("Done.")


plt.style.use("seaborn")
plt.plot(gd_cost, "b", label="Vanilla gradient descent")
plt.plot(qng_cost_mitigated, "g", label="Quantum natural gradient descent (with error mitigation)")

plt.ylabel("Cost function value")
plt.xlabel("Optimization steps")
plt.legend()
plt.show()
plt.plot(qng_cost, "b", label="QNGD")
plt.plot(qng_cost_mitigated, "g", label="QNGD with error mitigation")

plt.ylabel("Cost function value")
plt.xlabel("Optimization steps")
plt.legend()
plt.show()
cost_difference = [(a - b) for a, b in zip(qng_cost, qng_cost_mitigated)]

print('Average cost difference:', sum(cost_difference)/len(cost_difference))

plt.plot(cost_difference, "b", label="Cost difference (the more the better)")

plt.ylabel("Cost function value")
plt.xlabel("Optimization steps")
plt.ylim(-0.25, 0.25)
plt.legend()
plt.show()<jupyter_output>Average error difference: 0.059208984375<jupyter_text>Now let's compare these results with a new one using larger number of shots made during mitigation process<jupyter_code>with pennylane_extra.qiskit_measurement_error_mitigation(shots=4096):
    gd_cost = []
    opt = qml.GradientDescentOptimizer(0.01)

    print("Starting GD Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit, theta)
        gd_cost.append(circuit(theta))
    print("Done.")
    
    qng_cost_mitigated_4096 = []
    opt = qml.QNGOptimizer(0.01, diag_approx=True)

    print("Starting builtin QNG Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit, theta)
        qng_cost_mitigated_4096.append(circuit(theta))
    print("Done.")


plt.style.use("seaborn")
plt.plot(gd_cost, "b", label="Vanilla gradient descent")
plt.plot(qng_cost_mitigated_4096, "g", label="Quantum natural gradient descent (4096 shots)")

plt.ylabel("Cost function value")
plt.xlabel("Optimization steps")
plt.legend()
plt.show()
plt.figure(figsize=(16,12))

plt.plot(qng_cost, "r", label="QNGD")
plt.plot(qng_cost_mitigated, "b", label="QNGD with error mitigation (1024 mitigation shots)")
plt.plot(qng_cost_mitigated_4096, "g", label="QNGD with error mitigation (4096 mitigation shots)")

plt.ylabel("Cost function value")
plt.xlabel("Optimization steps")
plt.legend()
plt.show()
cost_difference = [(a - b) for a, b in zip(qng_cost_mitigated, qng_cost_mitigated_4096)]

print('Average cost difference:', sum(cost_difference)/len(cost_difference))

plt.plot(cost_difference, "b", label="Cost difference (the more the better)")
plt.ylabel("Cost function value")
plt.xlabel("Optimization steps")
plt.ylim(-0.25, 0.25)
plt.legend()
plt.show()<jupyter_output>Average error difference: -0.001484375<jupyter_text>Those were the sample configurations. Now we can measure the mean values for those experiments performed multiple times and compare them instead<jupyter_code>import pennylane_extra

n = 10
gd_cost_1024 = []
qng_cost_mitigated_1024 = []


with pennylane_extra.qiskit_measurement_error_mitigation(shots=1024):
    for i in range(n):
        opt = qml.GradientDescentOptimizer(0.01)

        print("Starting GD Optimizer run...")
        theta = init_params
        for _ in range(steps):
            theta = opt.step(circuit, theta)
            gd_cost_1024.append(circuit(theta))
        print("Done.")
        
        opt = qml.QNGOptimizer(0.01, diag_approx=True)

        print("Starting builtin QNG Optimizer run...")
        theta = init_params
        for _ in range(steps):
            theta = opt.step(circuit, theta)
            qng_cost_mitigated_1024.append(circuit(theta))
        print("Done.")
qng_cost_mitigated_4096 = []

with pennylane_extra.qiskit_measurement_error_mitigation(shots=4096):
    for i in range(n):
        opt = qml.QNGOptimizer(0.01, diag_approx=True)

        print("Starting builtin QNG Optimizer run...")
        theta = init_params
        for _ in range(steps):
            theta = opt.step(circuit, theta)
            qng_cost_mitigated_4096.append(circuit(theta))
        print("Done.")
qng_cost_no_mitigation = []
gd_cost_no_mitigation = []

for i in range(n):
    opt = qml.QNGOptimizer(0.01, diag_approx=True)

    print("Starting builtin QNG Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit, theta)
        qng_cost_no_mitigation.append(circuit(theta))
    print("Done.")
    
    opt = qml.GradientDescentOptimizer(0.01)

    print("Starting GD Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit, theta)
        gd_cost_no_mitigation.append(circuit(theta))
    print("Done.")
all_gd_cost_1024 = np.array(gd_cost_1024).reshape(n,-1)
mean_gd_cost_1024 = np.mean(all_gd_cost_1024, axis=0)

all_qng_cost_no_mitigation = np.array(qng_cost_no_mitigation).reshape(n,-1)
mean_qng_cost_no_mitigation = np.mean(all_qng_cost_no_mitigation, axis=0)

all_gd_cost_no_mitigation = np.array(gd_cost_no_mitigation).reshape(n,-1)
mean_gd_cost_no_mitigation = np.mean(all_gd_cost_no_mitigation, axis=0)

all_qng_cost_mitigated_1024 = np.array(qng_cost_mitigated_1024).reshape(n,-1)
mean_qng_cost_mitigated_1024 = np.mean(all_qng_cost_mitigated_1024, axis=0)

all_qng_cost_mitigated_4096 = np.array(qng_cost_mitigated_4096).reshape(n,-1)
mean_qng_cost_mitigated_4096 = np.mean(all_qng_cost_mitigated_4096, axis=0)
from matplotlib import pyplot as plt

def plot_opt_results(title, filename, *results):
    f = plt.figure(figsize=(16,12))
    plt.style.use("seaborn")
    for result in results:
        plt.plot(result[0], result[1], label=result[2])
    plt.title(title, fontsize=22)
    plt.ylabel("Cost function value", fontsize=22)
    plt.xlabel("Optimization steps", fontsize=22)
    plt.legend(prop={'size': 18})
    f.savefig(filename, bbox_inches='tight')
    plt.show()
plot_opt_results("QNG vs GD with noise", "qng_vs_gd_noise_mitigation.pdf",
                 (mean_qng_cost_mitigated_1024, "g", "Quantum natural gradient descent with error mitigation (1024 shots)"),
                 (mean_qng_cost_no_mitigation, "b", "Quantum natural gradient descent (no error mitigation)"),
                 (mean_gd_cost_1024, "black", "Vanilla gradient descent with error mitigation (1024 shots)"),
                 (mean_gd_cost_no_mitigation, "r", "Vanilla gradient descent (no error mitigation)"))
plot_opt_results("QNG (shots impact comparison)", "qng_vs_gd_noise_mitigation_shots.pdf",
                 (mean_qng_cost_mitigated_1024, "b", "Quantum natural gradient descent with error mitigation (1024 shots)"),
                 (mean_qng_cost_mitigated_4096, "r", "Quantum natural gradient descent with error mitigation (4096 shots)"))<jupyter_output><empty_output><jupyter_text>Below we also present optimization results without any noise (on an ideal simulator).<jupyter_code># Provision the the default device with no noise
dev_no_noise = qml.device("default.qubit", wires=3)


@qml.qnode(dev_no_noise)
def circuit_no_noise(params):
    # |psi_0>: state preparation
    qml.RY(np.pi / 4, wires=0)
    qml.RY(np.pi / 3, wires=1)
    qml.RY(np.pi / 7, wires=2)

    # V0(theta0, theta1): Parametrized layer 0
    qml.RZ(params[0], wires=0) 
    qml.RZ(params[1], wires=1)

    # W1: non-parametrized gates
    qml.CNOT(wires=[0, 1])
    qml.CNOT(wires=[1, 2])

    # V_1(theta2, theta3): Parametrized layer 1
    qml.RY(params[2], wires=1)
    qml.RX(params[3], wires=2)

    # W2: non-parametrized gates
    qml.CNOT(wires=[0, 1])
    qml.CNOT(wires=[1, 2])

    return qml.expval(qml.PauliY(0))
# same as above but as a separate cell
n = 10
steps = 200
init_params = np.array([0.432, -0.123, 0.543, 0.233])
qng_cost_no_noise = []
gd_cost_no_noise = []

for i in range(n):
    opt = qml.QNGOptimizer(0.01, diag_approx=True)

    print("Starting builtin QNG Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit_no_noise, theta)
        qng_cost_no_noise.append(circuit_no_noise(theta))
    print("Done.")
    
    opt = qml.GradientDescentOptimizer(0.01)

    print("Starting GD Optimizer run...")
    theta = init_params
    for _ in range(steps):
        theta = opt.step(circuit_no_noise, theta)
        gd_cost_no_noise.append(circuit_no_noise(theta))
    print("Done.")
all_qng_cost_no_noise = np.array(qng_cost_no_noise).reshape(n,-1)
mean_qng_cost_no_noise = np.mean(all_qng_cost_no_noise, axis=0)

all_gd_cost_no_noise = np.array(gd_cost_no_noise).reshape(n,-1)
mean_gd_cost_no_noise = np.mean(all_gd_cost_no_noise, axis=0)

plot_opt_results("QNG vs GD averaged without noise", "qng_vs_gd_no_noise_avg.pdf",
                 (mean_qng_cost_no_noise, "b", "Quantum natural gradient descent (no error mitigation)"),
                (mean_gd_cost_no_noise, "r", "Vanilla gradient descent (no error mitigation)"))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13358444877917536, 'fraction_numerical': 0.03571668796663752, 'mean_word_length': 4.2835820895522385, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 21, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '623245', 'license': 'Apache-2.0', 'n_tokens_mistral': 6080, 'n_tokens_neox': 5458, 'n_words': 1109, 'path': 'examples/qng_with_noise_demo.ipynb', 'repo_name': 'qkitchen/pennylane-extra', 'size': '14867'}