# What additional options are available in PEC?#

The application of probabilistic error cancellation (PEC) with Mitiq requires two main steps:

1. Building OperationRepresentation objects expressing ideal gates as linear combinations of noisy gates;

2. Estimating expectation values with PEC, making use of the representations obtained in the previous step.

Both steps can be achieved with Mitiq in different ways and with different options.

In the following code we use Qiskit as a frontend, but the workflow is the same for other frontends.

import numpy as np
import qiskit # Frontend library
from mitiq import pec  # Probabilistic error cancellation module


## Building OperationRepresentation objects#

Given the superoperator of an ideal gate $$\mathcal G$$, its quasi-probability representation is:

$\mathcal G = \sum_\alpha \eta_\alpha \mathcal O_\alpha$

where $$\{\eta_\alpha\}$$ are real coefficients and $$\{\mathcal O_\alpha\}$$ are the implementable noisy gates, i.e., those which can be actually applied by a noisy quantum computer.

For trace-preserving operations, the coefficients $$\{\eta_\alpha\}$$ form a quasi-probability distribution, i.e.:

$\sum_\alpha \eta_\alpha = 1, \quad \gamma = \| \eta \|_1= \sum_\alpha |\eta_\alpha| \ge 1.$

The value of $$\gamma$$ is related to the negative “volume” of the distribution and quantifies to the sampling cost of PEC (see What is the theory behind PEC?).

### Defining NoisyOperation objects#

The noisy operations $$\{\mathcal O_\alpha\}$$ in the equation above can correspond to single noisy gates. However, it is often useful to define a noisy operation as a sequence multiple noisy gates. To have this flexibility, we associate to each noisy operation a small QPROGRAM, i.e., a quantum circuit defined via a supported frontend. For example a basis of noisy operations that are useful to represent the Hadamard gate in the presence of depolarizing noise is:

basis_circuit_h = qiskit.QuantumCircuit(1)
basis_circuit_h.h(0)

basis_circuit_hx = qiskit.QuantumCircuit(1)
basis_circuit_hx.h(0)
basis_circuit_hx.x(0)

basis_circuit_hy = qiskit.QuantumCircuit(1)
basis_circuit_hy.h(0)
basis_circuit_hy.y(0)

basis_circuit_hz = qiskit.QuantumCircuit(1)
basis_circuit_hz.h(0)
basis_circuit_hz.z(0)

basis_circuits = [basis_circuit_h, basis_circuit_hx, basis_circuit_hy, basis_circuit_hz]

for c in basis_circuits:
print(c)

   ┌───┐
q: ┤ H ├
└───┘
┌───┐┌───┐
q: ┤ H ├┤ X ├
└───┘└───┘
┌───┐┌───┐
q: ┤ H ├┤ Y ├
└───┘└───┘
┌───┐┌───┐
q: ┤ H ├┤ Z ├
└───┘└───┘


Each element of basis_circuits describes “how to physically implement” a noisy operation $$\mathcal O_\alpha$$ on a noisy backend. To completely characterize a noisy operation we can also specify the actual (non-unitary) quantum channel associated to it. In Mitiq, this can be done using the NoisyOperation class.

For example, assuming that each of the previous basis circuits is affected by a final depolarizing channel, the following code cell generates the corresponding NoisyOperation objects.

from mitiq.pec.representations import local_depolarizing_kraus
from mitiq.pec.channels import kraus_to_super

# Compute depolarizing superoperator
BASE_NOISE = 0.2
depo_super = kraus_to_super(local_depolarizing_kraus(BASE_NOISE, num_qubits=1))

# Define the superoperator matrix of each noisy operation
super_matrices = [
depo_super @ kraus_to_super([qiskit.quantum_info.Operator(c).data])
for c in basis_circuits
]

# Define NoisyOperation objects combining circuits with channel matrices
noisy_operations = [
pec.NoisyOperation(circuit=c, channel_matrix=m)
for c, m in zip(basis_circuits, super_matrices)
]

print(f"{len(noisy_operations)} NoisyOperation objects defined.")

4 NoisyOperation objects defined.


Note: A NoisyOperation can also be instantiated with channel_matrix=None. In this case, however, the quasi-probability distribution must be known to the user and cannot be derived by Mitiq with the procedure shown in the next section.

### Finding an optimal OperationRepresentation#

Similar to what we did for basis_circuits, we also define the ideal_operation that we aim to represent in the form of a QPROGRAM. Assuming that we aim to represent the Hadamard gate, we have:

ideal_operation = qiskit.QuantumCircuit(1)
ideal_operation.h(0)
print(f"The ideal operation to expand in the noisy basis is:\n{ideal_operation}")

The ideal operation to expand in the noisy basis is:
┌───┐
q: ┤ H ├
└───┘


The Mitiq function find_optimal_representation() can be used to numerically obtain an OperationRepresentation of the ideal_operation in the basis of the noisy implementable gates (noisy_operations).

from mitiq.pec.representations import find_optimal_representation

h_rep = find_optimal_representation(ideal_operation, noisy_operations)
print(f"Optimal representation:\n{h_rep}")

Optimal representation:
q_0: ───H─── = 1.273*(   ┌───┐
q: ┤ H ├
└───┘)-0.091*(   ┌───┐┌───┐
q: ┤ H ├┤ X ├
└───┘└───┘)-0.091*(   ┌───┐┌───┐
q: ┤ H ├┤ Y ├
└───┘└───┘)-0.091*(   ┌───┐┌───┐
q: ┤ H ├┤ Z ├
└───┘└───┘)


The representation is optimal in the sense that, among all the possible representations, it minimizes the one-norm of the quasi-probability distribution. Behind the scenes, find_optimal_representation() solves the following optimization problem:

$\begin{split}\gamma^{\rm opt} = \min_{\substack{ \{ \eta_{\alpha} \} \\ \{ \mathcal O_{ \alpha} \}}} \left[ \sum_\alpha |\eta_{\alpha}| \right], \quad \text{ such that} \quad \mathcal G = \sum_\alpha \eta_\alpha \mathcal O_\alpha \, .\end{split}$

### Manually defining an OperationRepresentation#

Instead of solving the previous optimization problem, an OperationRepresentation can also be manually defined. This approach can be applied if the user already knows the quasi-probability distribution $${\eta_\alpha}$$.

# We assume to know the quasi-distribution
quasi_dist = h_rep.coeffs

# Manual definition of the OperationRepresentation
manual_h_rep = pec.OperationRepresentation(
ideal=ideal_operation,
noisy_operations=noisy_operations,
coeffs=quasi_dist,
)

# Test that the manual definition is equivalent to h_rep
assert manual_h_rep == h_rep


Note: For the particular case of depolarizing noise, Mitiq can directly create the OperationRepresentation of an arbitrary ideal_operation, as shown in the next cell.

from mitiq.pec.representations.depolarizing import represent_operation_with_local_depolarizing_noise

depolarizing_h_rep = represent_operation_with_local_depolarizing_noise(
ideal_operation,
noise_level=BASE_NOISE,
)

assert depolarizing_h_rep == h_rep


### Qubit-independent representations#

It is possible to define a qubit-independent OperationRepresentation by setting the option is_qubit_dependent to False.

In this case, a signle OperationRepresentation representing a gate acting on some arbitrary qubits can be used to mitigate the same gate even if acting on different qubits.

qreg = qiskit.QuantumRegister(2)
circuit = qiskit.QuantumCircuit(qreg)
circuit.h(0)
circuit.h(1)

# OperationRepresentation defined on different qubits
rep_qreg = qiskit.QuantumRegister(1, "rep_reg")
ideal_op = qiskit.QuantumCircuit(rep_qreg)
ideal_op.h(rep_qreg)
hxcircuit = qiskit.QuantumCircuit(rep_qreg)
hxcircuit.h(0)
hxcircuit.x(0)
hzcircuit = qiskit.QuantumCircuit(rep_qreg)
hzcircuit.h(0)
hzcircuit.z(0)
noisy_hxop = pec.NoisyOperation(hxcircuit)
noisy_hzop = pec.NoisyOperation(hzcircuit)

rep = pec.OperationRepresentation(
ideal=ideal_op,
noisy_operations=[noisy_hxop, noisy_hzop],
coeffs=[0.5, 0.5],
is_qubit_dependent=False,
)
print(rep)
print("Using the same rep on a circuit with H gates acting on different qubits:")
sampled_circuits, _, _ = pec.sample_circuit(circuit, representations=[rep])
print(*sampled_circuits)

rep_reg_0: ───H─── = 0.500*(         ┌───┐┌───┐
rep_reg: ┤ H ├┤ X ├
└───┘└───┘)+0.500*(         ┌───┐┌───┐
rep_reg: ┤ H ├┤ Z ├
└───┘└───┘)
Using the same rep on a circuit with H gates acting on different qubits:

      ┌───┐┌───┐
q0_0: ┤ H ├┤ Z ├
├───┤├───┤
q0_1: ┤ H ├┤ Z ├
└───┘└───┘


### Methods of the OperationRepresentation class#

The main idea of PEC is to estimate the average with respect to a quasi-probability distribution over noisy circuits with a probabilistic Monte-Carlo approach. This can be obtained rewriting $\mathcal G = \sum_\alpha \eta_\alpha \mathcal O_\alpha$ as:

$\mathcal G = \gamma \sum_\alpha p(\alpha) \textrm{sign}(\eta_\alpha) \mathcal O_\alpha \quad p(\alpha):= |\eta_\alpha| / \gamma,$

where $$p(\alpha)$$ is a (positive) well-defined probability distribution. If we take a single sample from $$p(\alpha)$$, we obtain a noisy operation $$\mathcal O_\alpha$$ that should be multiplied by the sign of the associated coefficient $$\eta_\alpha$$ and by $$\gamma$$.

The method OperationRepresentation.sample() can be used for this scope:

noisy_op, sign, coeff = h_rep.sample()
print(f"The sampled noisy operation is: {noisy_op}.")
print(f"The associated coefficient is {coeff:g}, whose sign is {sign}.")

The sampled noisy operation is:    ┌───┐┌───┐
q: ┤ H ├┤ X ├
└───┘└───┘.
The associated coefficient is -0.0909091, whose sign is -1.


Note: try re-executing the previous cell to get different samples.

Other useful methods of OperationRepresentation are shown in the next cells.

# One-norm "gamma" quantifying the mitigation cost
h_rep.norm

1.545454514545456

# Quasi-probability distribution
print(h_rep.coeffs)
assert sum([abs(eta) for eta in h_rep.coeffs]) == h_rep.norm

[1.272727247272726, -0.09090908909091, -0.09090908909090956, -0.09090908909091]

# Positive and normalized distribution p(alpha)=|eta_alpha|/gamma
h_rep.distribution

[0.8235294117647043,
0.05882352941176524,
0.058823529411764955,
0.05882352941176524]


## Estimating expectation values with PEC#

The second main step of PEC is to make use of the previously defined OperationRepresentations to estimate expectation values with respect to a quantum state prepared by a circuit of interest. In the previous section we defined the representation of the Hadamard gate. So, for simplicity, we consider a circuit that contains only Hadamard gates.

### Defining a circuit of interest and an Executor#

circuit = qiskit.QuantumCircuit(1)
for _ in range(4):
circuit.h(0)
print(circuit)

   ┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ H ├
└───┘└───┘└───┘└───┘


In this case, the list of OperationRepresentations that we need for PEC is simply:

representations = [h_rep]


In general representations will contain as many representations as the number of ideal gates involved in circuit.

Note: If a gate is in circuit but its OperationRepresentation is not listed in representations, Mitiq can still apply PEC. However, any errors associated to that gate will not be mitigated. In practice, all the gates without OperationRepresentations are treated by Mitiq as if they were noiseless.

The executor must be defined by the user since it depends on the specific frontend and backend (see the Executors section). Here, for simplicity, we import the basic execute_with_noise() function from the Qiskit utilities of Mitiq.

from mitiq import Executor
from mitiq.interface.mitiq_qiskit import execute_with_noise, initialized_depolarizing_noise

def execute(circuit):
"""Returns Tr[ρ |0⟩⟨0|] where ρ is the state prepared by the circuit
executed with depolarizing noise.
"""
circuit_copy = circuit.copy()
noise_model = initialized_depolarizing_noise(BASE_NOISE)
projector_on_zero = np.array([[1, 0], [0, 0]])
return execute_with_noise(circuit_copy, projector_on_zero, noise_model)

# Cast the execute function into an Executor class to record the execution history
executor = Executor(execute)


### Options for estimating expectation values#

In the “How do I use PEC?” section, we have shown how to apply PEC with the minimum amount of options by simply calling the function execute_with_pec() with the basic arguments circuit, executor, and representations.

In the next code-cell, we show additional options that can be used:

pec_value, pec_data = pec.execute_with_pec(
circuit,
executor,
observable=None, # In this example the observable is implicit in the executor
representations=representations,
num_samples=5, # Number of PEC samples
random_state=0, # Seed for reproducibility of PEC sampling
full_output=True, # Return "pec_data" in addition to "pec_value"
)


Similar to other error mitigation modules, observable is an optional argument of execute_with_pec(). If observable is None the executor must return an expectation value, otherwise the executor must return a mitiq.QuantumResult from which the expectation value of the input observable can be computed. See the Executors section for more details.

Another option that can be used, instead of num_samples, is precision. Its default value is 0.03 and quantifies the desired estimation accuracy.

For a bounded observable $$\|A\|\le 1$$, precision approximates $$|\langle A \rangle_{ \rm ideal} - \langle A \rangle_{ \rm PEC}|$$ (up to constant factors and up to statistical fluctuations). In practice, precision is used by Mitiq to automatically determine num_samples, according to the formula: num_samples = $$(\gamma /$$ precision$$)^2$$, where $$\gamma$$ is the one-norm the circuit quasi-probability distribution. See “What is the theory behind PEC?” for more details on the sampling cost.

# Optional Executor re-initialization to clean the history
executor = Executor(execute)

pec_value, pec_data = pec.execute_with_pec(
circuit,
executor,
observable=None, # In this example, the observable is implicit in the executor.
representations=representations,
precision=0.5, # The estimation accuracy.
random_state=0, # Seed for reproducibility of probabilistic sampling of circuits.
full_output=True, # Return pec_data in addition to pec_value
)


Hint: The value of precision used above is very large, in order to reduce the execution time. Try re-executing the previous cell with smaller values of precision to improve the result.

### Analyzing the executed circuits#

As discussed in the Executors section, we can extract the execution history from executor. This is a way to see what Mitiq does behind the scenes which is independent from the error mitigation technique.

print(
f"During the previous PEC process, {len(executor.executed_circuits)} ",
"circuits have been executed."
)
print(f"The first 5 circuits are:")

for c in executor.executed_circuits[:5]:
print(c)

print(f"The corresponding noisy expectation values are:")
for c in executor.quantum_results[:5]:
print(c)

During the previous PEC process, 130  circuits have been executed.
The first 5 circuits are:
┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ H ├
└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ H ├
└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ H ├
└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ Y ├┤ H ├┤ X ├
└───┘└───┘└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ H ├
└───┘└───┘└───┘└───┘
The corresponding noisy expectation values are:
0.7047999999999993
0.7047999999999993
0.7047999999999993
0.6310719999999995
0.7047999999999993


### Analyzing data of the PEC process#

Beyond the executed circuits, one may be interested in analyzing additional data related to the specifc PEC technique. Setting full_output=True, this data is returned in pec_data as a dictionary.

print(pec_data.keys())

dict_keys(['num_samples', 'precision', 'pec_value', 'pec_error', 'unbiased_estimators', 'measured_expectation_values', 'sampled_circuits'])

print(pec_data["num_samples"], "noisy circuits have been sampled and executed.")

130 noisy circuits have been sampled and executed.


The unbiased raw results, whose average is equal to pec_value, are stored under the unbiased_estimators key. For example, the first 5 unbiased samples are:

pec_data["unbiased_estimators"][:5]

[4.020599418806928,
4.020599418806928,
4.020599418806928,
3.6000109483900764,
4.020599418806928]


The statistical error pec_error corresponds to pec_std / sqrt(num_samples), where pec_std is the standard deviation of the unbiased samples, i.e., the square root of the mean squared deviation of unbiased_estimators from pec_value.

pec_data["pec_error"]

0.28902565876334807


Instead of the error printed above, one could use more advanced statistical techniques to estimate the uncertainty of the PEC result. For example, given the raw samples contained in pec_data["unbiased_estimators"], one can apply a bootstrapping approach. Alternatively, a simpler but computationally more expensive approach is to perform multiple PEC estimations of the same expectation value and compute the standard deviation of the results.

## Applying learning-based PEC#

Mitiq also contains functions for learning the quasiprobability representations of an ideal gate from Clifford circuit simulations, with biased noise or depolarizing noise models. The resulting quasiprobability representations can be used to obtain an error-mitigated expectation value via the representations argument of pec.execute_with_pec(). The learning-based PEC workflow was inspired by the procedure described in Strikis et al. PRX Quantum (2021) .

The learning process is based on the execution of a set of training circuits on a noisy backend via a noisy executor and on a classical simulator via an ideal executor. The training circuits are near-Clifford approximations of the input circuit. During training, the noise strength parameter is used to calculate quasiprobability representations of the ideal gate with a depolarizing noise model. The representations are then input into pec.execute_with_pec() to obtain an error-mitigated expecation value from execution of the training circuit, for comparison with the ideal expecation value obtained from classical simulation of the training circuit. The optimizer used in the learning function is from scipy.optimize.minimize(). The default optimization method in the learning function is Nelder-Mead, as that appears to work best for this particular problem setup.

In addition to specifying the input operation, the circuit of interest, and the ideal and noisy executors, the user should specify the number of training circuits, the fraction of non-Clifford gates in the training circuits, an initial guess for noise strength, and in the case of biased noise, an initial guess for a noise bias. The user can also set options for the intermediate executions of pec.execute_with_pec() during the training process as a dictionary in pec_kwargs, specify the observable of which the expecation value is to be computed, and enter a dictionary of additional data and options including optimization method (supported by scipy.optimize.minimize()) and settings for the chosen optimization method.

Note

Using learn_depolarizing_noise_parameter and learn_biased_noise_parameters may require some tuning. One of the main challenges is setting a good value of num_samples in the PEC options pec_kwargs. Setting a small value of num_samples is typically necessary to obtain a reasonable execution time. On the other hand, using a number of PEC samples that is too small can result in a large statistical error, ultimately causing the optimization process to fail.

### Learning quasiprobability representations with a depolarizing noise model#

In cases where the noise of the backend can be approximated by a depolarizing noise model, pec.representations.learning.learn_depolarizing_noise_parameter() can be used to learn the noise strength epsilon associated to a set of input operations.

from qiskit_aer.noise import NoiseModel
from qiskit_aer.noise.errors.standard_errors import (
depolarizing_error,
)
from mitiq import Observable, PauliString
from mitiq.pec.representations.learning import (
learn_depolarizing_noise_parameter,
)

circuit = qiskit.QuantumCircuit(2)
circuit.rx(1.14 * np.pi, 0)
circuit.rz(0, 0)
circuit.cx(1, 0)
circuit.rx(1.71 * np.pi, 0)
circuit.rx(1.14 * np.pi, 1)

observable = Observable(PauliString("XZ"), PauliString("YY")).matrix()

# set up ideal simulator
def ideal_execute(circuit):
"""Simulate (training) circuits without noise"""
circuit_copy = circuit.copy()
noise_model = initialized_depolarizing_noise(0.0)
return execute_with_noise(circuit_copy, observable, noise_model)

epsilon = 0.05  # noise level for simulation in noisy executor

def noisy_execute(circuit):
"""Simulate circuit with a depolarizing noise model"""
circuit_copy = circuit.copy()
noise_model = NoiseModel()
return execute_with_noise(circuit_copy, observable, noise_model)

operations_to_learn = qiskit.QuantumCircuit(2)
operations_to_learn.cx(1, 0)

[success, epsilon_opt] = learn_depolarizing_noise_parameter(
operations_to_learn=[operations_to_learn],  # learn rep of Hadamard
circuit=circuit,
ideal_executor=Executor(ideal_execute),
noisy_executor=Executor(noisy_execute),
pec_kwargs={"num_samples": 500},
num_training_circuits=5,
fraction_non_clifford=0.2,

epsilon0=0.9 * epsilon,  # initial guess for epsilon
)


Upon completing the optimization loop, representations.learning.learn_depolarizing_noise_parameter() returns a flag indicating whether the optimizer exited successfully, in addition to the optimized noise strength, which can then be input into representations.depolarizing.represent_operation_with_local_depolarizing_noise() to calculate quasiprobability representations of the ideal gate in terms of noisy gates.

representations = represent_operation_with_local_depolarizing_noise(
operations_to_learn, epsilon_opt
)


### Learning quasiprobability representations with a biased noise model#

In cases where the noise of the backend can be approximated by a combined depolarizing and dephasing noise model with a bias factor, also referred to as a biased noise model, pec.representations.learning.learn_biased_noise_parameters can be used to learn the noise strength epsilon and noise bias eta associated to a set of input operations.

The single-qubit biased noise channel is given by:

$\mathcal{D}(\epsilon) = (1 - \epsilon)\mathbb{1} + \epsilon\Big(\frac{\eta}{\eta + 1} \mathcal{Z} + \frac{1}{3}\frac{1}{\eta + 1}(\mathcal{X} + \mathcal{Y} + \mathcal{Z})\Big)$

For multi-qubit operations, the noise channel used is the tensor product of the local single-qubit channels.

from mitiq.pec.representations.learning import (
learn_biased_noise_parameters,
)

epsilon = 0.05
eta = 0

# simulate biased noise occurs on the CNOT gates
def noisy_execute(circuit):
circuit_copy = circuit.copy()
noise_model = NoiseModel()
return execute_with_noise(circuit_copy, observable, noise_model)

operations_to_learn = qiskit.QuantumCircuit(2)
operations_to_learn.cx(1, 0)

[success, epsilon_opt, eta_opt] = learn_biased_noise_parameters(
operations_to_learn=[operations_to_learn],  # learn rep of CNOT
circuit=circuit,
ideal_executor=Executor(ideal_execute),
noisy_executor=Executor(noisy_execute),
pec_kwargs={"num_samples": 500, "random_state": 1},
num_training_circuits=5,
fraction_non_clifford=0.2,
training_random_state=np.random.RandomState(1),
epsilon0=1.01 * 0.05,  # initial guess for noise strength
eta0= eta + 0.01,  # initial guess for noise bias
)


Upon completing the optimization loop, representations.learning.learn_biased_noise_parameters() returns a flag indicating whether the optimizer exited successfully, in addition to the optimized noise strength and noise bias, which can then be input into represent_operation_with_local_biased_noise() to calculate quasiprobability representations of the ideal gate in terms of noisy gates.

from pec.representations import (
represent_operation_with_local_biased_noise,
)

representations = represent_operation_with_local_biased_noise(
operations_to_learn, epsilon_opt, eta_opt
)