Implementation of QNG on circuits that have been mitigated with ZNE

Hello pennylane help team, I am having difficulty applying QNG to a circuit that has been mitigated using ZNE, can you help me to solve this problem, I am very grateful for that.

import pennylane as qml
from pennylane import numpy as np
import matplotlib.pyplot as plt
from pennylane.transforms import mitigate_with_zne

# Define the molecule H2
symbols = ["H", "H"]
coordinates = np.array([0.0, 0.0, 0, 0.0, 0.0, 1.4])  # Bond length = 0.7 Å

# Load the electronic Hamiltonian for H2
H, qubits = qml.qchem.molecular_hamiltonian(
    symbols, coordinates, charge=0, basis="sto-3g"
)

# Describe noise model
fcond = qml.noise.wires_in(range(qubits))
noise = qml.noise.partial_wires(qml.DepolarizingChannel, 0.05)
noise_model = qml.NoiseModel({fcond: noise})

# Define ideal and noisy devices
dev_ideal = qml.device("default.mixed", wires=qubits)
dev_noisy = qml.add_noise(dev_ideal, noise_model=noise_model)

# define Ansatz
hf_state = np.array([1, 1, 0, 0], requires_grad=False)

def ansatz(params, wires=[0, 1, 2, 3]):
    for i in wires:
        qml.RZ(params[3 * i], wires=i)
        qml.RY(params[3 * i + 1], wires=i)
        qml.RZ(params[3 * i + 2], wires=i)
    qml.CNOT(wires=[2, 3])
    qml.CNOT(wires=[2, 0])
    qml.CNOT(wires=[3, 1])

# Initialize parameters
np.random.seed(0)
init_params = np.random.uniform(0, 2 * np.pi, size=12, requires_grad=True)
step_size = 0.5

# Define cost functions
@qml.qnode(dev_noisy, diff_method="parameter-shift")  
def cost_noisy(params):
    ansatz(params, wires=range(qubits))
    return qml.expval(H)

# define mitigate with zne
scale_factors = [1, 2, 3]

# cost mitigated func
cost_mitigated = mitigate_with_zne(cost_noisy,
    scale_factors=scale_factors,
    folding=qml.transforms.fold_global,
    extrapolate=qml.transforms.richardson_extrapolate,
)


opt_mitigated = qml.QNGOptimizer(step_size)

# Optimization loop for mitigated case
params_mitigated = init_params.copy()
energies_mitigated = []
for n in range(50):
    params_mitigated, energy = opt_mitigated.step_and_cost(cost_mitigated, params_mitigated)
    energies_mitigated.append(energy)
    print(f"Noisy Iteration {n + 1}: Energy = {energy:.8f} Ha")

This is the error from the code above

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-140-7a8f4d29cb27> in <cell line: 0>()
      5 energies_mitigated = []
      6 for n in range(50):
----> 7     params_mitigated, energy = opt_mitigated.step_and_cost(cost_mitigated, params_mitigated)
      8     energies_mitigated.append(energy)
      9     print(f"Noisy Iteration {n + 1}: Energy = {energy:.8f} Ha")

14 frames
/usr/local/lib/python3.11/dist-packages/pennylane/_grad.py in new_f(*args, **kwargs)
    259         output = f(*args, **kwargs)
    260         if output.__class__.__module__.split(".")[0] not in {"autograd", "pennylane", "numpy"}:
--> 261             raise ValueError(
    262                 f"autograd can only differentiate with respect to arrays, not {type(output)}. Ensure the output class is an autograd array."
    263             )

ValueError: autograd can only differentiate with respect to arrays, not <class 'tuple'>. Ensure the output class is an autograd array.

Hi @Dwi_Cahyo_Mariyanto ,

Thank you for posting your question here. I’m not sure why this is happening but I can indeed replicate your issue.

Let me check and get back to you.

1 Like

Hi @Dwi_Cahyo_Mariyanto,

I don’t know exactly why this is happening but this looks related to the metric tensor function.

Do you need to use QNG? Or are you ok with using AdamOptimizer for example? I’m asking because I don’t know how long it can take to sort out the issue with QNG.

1 Like

Hi @CatalinaAlbornoz

Thank you for helping me, yes I need QNG in my research, maybe can you recommend me to use another mitigation where I can still use QNG?

Hi @Dwi_Cahyo_Mariyanto ,

Changing the mitigation strategy can be a bit complicated but it might be the way to go in this case indeed.

Here’s why it currently doesn’t work:
In order to mitigate you need to run the circuit with different levels of noise. In order to do this you “fold” the circuit, creating versions of it that include more gates, and thus the noise depends directly on the scaling factor. If you draw the mitigated qnode (e.g. see the code below), you’ll see that the mitigated qnode is in fact a tuple of 3 circuits.

import pennylane as qml
from pennylane import numpy as np
import matplotlib.pyplot as plt
from pennylane.transforms import mitigate_with_zne

# Define the molecule H2
symbols = ["H", "H"]
coordinates = np.array([0.0, 0.0, 0, 0.0, 0.0, 1.4])  # Bond length = 0.7 Å

# Load the electronic Hamiltonian for H2
H, qubits = qml.qchem.molecular_hamiltonian(
    symbols, coordinates, charge=0, basis="sto-3g"
)

# Describe noise model
fcond = qml.noise.wires_in(range(qubits))
noise = qml.noise.partial_wires(qml.DepolarizingChannel, 0.05)
noise_model = qml.NoiseModel({fcond: noise})

# Define ideal and noisy devices
dev_ideal = qml.device("default.mixed", wires=qubits)
dev_noisy = qml.add_noise(dev_ideal, noise_model=noise_model)

# define Ansatz
hf_state = np.array([1, 1, 0, 0], requires_grad=False)

def ansatz(params, wires=[0, 1, 2, 3]):
    for i in wires:
        qml.RZ(params[3 * i], wires=i)
        qml.RY(params[3 * i + 1], wires=i)
        qml.RZ(params[3 * i + 2], wires=i)
    qml.CNOT(wires=[2, 3])
    qml.CNOT(wires=[2, 0])
    qml.CNOT(wires=[3, 1])

# Initialize parameters
np.random.seed(0)
init_params = np.random.uniform(0, 2 * np.pi, size=12, requires_grad=True)
step_size = 0.5

# Define cost functions
@qml.qnode(dev_noisy, diff_method="parameter-shift")  
def cost_noisy(params):
    ansatz(params, wires=range(qubits))
    return qml.expval(H)

# define mitigate with zne
scale_factors = [1, 2, 3]

# cost mitigated func
cost_mitigated = mitigate_with_zne(cost_noisy,
    scale_factors=scale_factors,
    folding=qml.transforms.fold_global,
    extrapolate=qml.transforms.richardson_extrapolate,
)

# Draw the mitigated qnode
print(qml.draw(cost_mitigated)(init_params))

Output:

0: ──RZ(3.45)──RY(4.49)──RZ(3.79)────╭X────┤ ╭<𝓗>
1: ──RZ(3.42)──RY(2.66)──RZ(4.06)────│──╭X─┤ ├<𝓗>
2: ──RZ(2.75)──RY(5.60)──RZ(6.05)─╭●─╰●─│──┤ ├<𝓗>
3: ──RZ(2.41)──RY(4.97)──RZ(3.32)─╰X────╰●─┤ ╰<𝓗>

0: ──RZ(3.45)──RY(4.49)──RZ(3.79)────╭X────────╭X†───────────────────────────────────────────────
1: ──RZ(3.42)──RY(2.66)──RZ(4.06)────│──╭X─╭X†─│─────────────────────────────────────────────────
2: ──RZ(2.75)──RY(5.60)──RZ(6.05)─╭●─╰●─│──│───╰●──╭●───RZ(6.05)†──RY(5.60)†──RY(5.60)───RZ(6.05)
3: ──RZ(2.41)──RY(4.97)──RZ(3.32)─╰X────╰●─╰●──────╰X†──RZ(3.32)†──RY(4.97)†──RZ(2.41)†──RZ(2.41)

─────────────────────────╭X────┤ ╭<𝓗>
─────────────────────────│──╭X─┤ ├<𝓗>
──────────────────────╭●─╰●─│──┤ ├<𝓗>
───RY(4.97)──RZ(3.32)─╰X────╰●─┤ ╰<𝓗>

0: ──RZ(3.45)──RY(4.49)──RZ(3.79)────╭X────────╭X†──RZ(3.79)†──RY(4.49)†──RZ(3.45)†──RZ(3.45)─
1: ──RZ(3.42)──RY(2.66)──RZ(4.06)────│──╭X─╭X†─│────RZ(4.06)†──RY(2.66)†──RZ(3.42)†──RZ(3.42)─
2: ──RZ(2.75)──RY(5.60)──RZ(6.05)─╭●─╰●─│──│───╰●──╭●──────────RZ(6.05)†──RY(5.60)†──RZ(2.75)†
3: ──RZ(2.41)──RY(4.97)──RZ(3.32)─╰X────╰●─╰●──────╰X†─────────RZ(3.32)†──RY(4.97)†──RZ(2.41)†

───RY(4.49)──RZ(3.79)──────────────╭X────┤ ╭<𝓗>
───RY(2.66)──RZ(4.06)──────────────│──╭X─┤ ├<𝓗>
───RZ(2.75)──RY(5.60)──RZ(6.05)─╭●─╰●─│──┤ ├<𝓗>
───RZ(2.41)──RY(4.97)──RZ(3.32)─╰X────╰●─┤ ╰<𝓗>

Some optimizers can handle this, but for QNGO the issue is that you need to calculate the metric tensor for a single qnode, hence why it doesn’t work in this situation.

One (possible) solution:
Note that this is just an idea, so I don’t know if it will work. You could try implementing the mitigation yourself (manually instead of using PennyLane’s functions) and using QNGO for some new cost function that you’d need to define.
I don’t know if this is practical nor if it’s the best solution, but it’s something you could explore.

1 Like

Hi @Dwi_Cahyo_Mariyanto,

My colleague Josh suggested another strategy that could (maybe) work.
You can define a gradient transform (e.g., by inheriting from the param-shift transform) that returns quantum_gradient(tape) * metric_tensor(tape) .Then, in your QNode you can set diff_method=custom_qngrad , and use standard GradientDescent optimization. This is not straightforward, but it’s an option in case you want to try it.

1 Like