# What additional options are available in PEC?#

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

Building

`OperationRepresentation`

objects expressing ideal gates as linear combinations of noisy gates;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:

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

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:

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

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 `OperationRepresentation`

s 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 `OperationRepresentation`

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