# 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#

Combining the noisy operations defined above, we obtain a NoisyBasis.

noisy_basis = pec.NoisyBasis(*noisy_operations)


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_basis).

from mitiq.pec.representations import find_optimal_representation

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

Optimal representation:
q_0: ───H─── = -0.091*(q_0: ───H───Z───)-0.091*(q_0: ───H───X───)-0.091*(q_0: ───H───Y───)+1.273*(q_0: ───H───)


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

# This is just a reordering of noisy_operations to match quasi_dist
reordered_noisy_operations = h_rep.noisy_operations

# Manual definition of the OperationRepresentation
basis_expansion = dict(zip(reordered_noisy_operations, quasi_dist))
manual_h_rep = pec.OperationRepresentation(
ideal=ideal_operation, basis_expansion=basis_expansion
)

# 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


### 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_0: ───H───Z───.
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.545454514545455

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

(-0.09090908909090867, -0.09090908909090911, -0.09090908909090867, 1.2727272472727287)

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

array([0.05882353, 0.05882353, 0.05882353, 0.82352941])


## 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 ├┤ X ├
└───┘└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ X ├┤ H ├
└───┘└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ Y ├┤ H ├┤ Z ├┤ H ├
└───┘└───┘└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ H ├
└───┘└───┘└───┘└───┘
┌───┐┌───┐┌───┐┌───┐┌───┐
q: ┤ H ├┤ H ├┤ H ├┤ H ├┤ Z ├
└───┘└───┘└───┘└───┘└───┘
The corresponding noisy expectation values are:
0.33615999999999974
0.6638399999999994
0.6310719999999995
0.7047999999999993
0.6638399999999994


### 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]

[-1.9176570667226647,
-3.786939157464224,
3.6000109483900684,
4.020599418806919,
-3.786939157464224]


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.2840201074783067


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.