Final Project: Pauli Twirling

Final project: Pauli twirling

Pauli twirling, also known as randomized compiling, is a common step applied before many noise mitigation techniques used in the field of quantum computing. With Pauli twirling it is possible to tailor general noise channels into a stochastic Pauli channels, which have a relatively simple model, while at the same time keeping the logical circuit unchanged.

In this final project, we will explore the effect of Pauli twirling on the different noise channels that we have covered in the course, implement a custom experiment class to be used in conjunction with other qiskit-experiments classes and simulate the effect of Pauli twirling on the noise channels.

Pauli twirling was first introduced as randomized compiling in J. Wallman and J. Emerson, Phys. Rev. A 94, 052325 (2015). It works by producing classical mixtures of the noise channel with gates applied on either side of the noise channel. The gates are often chosen to be the Pauli gates $\{I, X, Y, Z\}$. The twirling action $\mathcal{T}_i$ on a quantum channel $\Lambda(\rho)$ is described as \begin{align} \mathcal{T}_i\Lambda &=\mathcal{P}_i\Lambda\mathcal{P}_i^\dagger \end{align} where the super-operator $\mathcal{P}_i$ acts on an operator as \begin{align} \mathcal{P}_i\Lambda&= \sigma_i^\dagger\Lambda\sigma_i. \end{align}

The resulting channel is then produced by averaging over all the twirled channels over the chosen gates, in our case the Pauli gates \begin{align} \widetilde{\Lambda} &= \frac{1}{4}\sum_{i=0}^3{\mathcal{T}_i\Lambda} \\ &= \frac{\Lambda + \mathcal{P}_{x}\Lambda \mathcal{P}_{x}^\dagger + \mathcal{P}_{y}\Lambda \mathcal{P}_{y}^\dagger + \mathcal{P}_{z} \Lambda \mathcal{P}_{z}^\dagger}{4}. \end{align}

\begin{align} \widetilde{\Lambda}(\rho) &= \frac{1}{4}\sum_{i=0}^3{\mathcal{T}_i\Lambda(\rho)} \\ &= \frac{1}{4}\sum_{i=0}^3{\mathcal{P}_i(\Lambda(\mathcal{P}_i^\dagger(\rho)))} \\ &= \frac{\Lambda \rho + \mathcal{P}_{x}\Lambda \mathcal{P}_{x}^\dagger \rho + \mathcal{P}_{y}\Lambda \mathcal{P}_{y}^\dagger \rho + \mathcal{P}_{z} \Lambda \mathcal{P}_{z}^\dagger \rho}{4}. \end{align}

We start with a density matrix $\rho$ and apply a channel $\lambda$ on it: $$\Lambda \rho = \Lambda(\rho)$$ Next, we twirl only the channel: $$\widetilde{\Lambda} = \mathcal{T}_i\Lambda = \mathcal{P}_i\Lambda\mathcal{P}_i^\dagger$$ When we apply this new channel on $\rho$ again, we obtain $$\widetilde{\Lambda}(\rho) = \mathcal{P}_i\Lambda\mathcal{P}_i^\dagger\rho = \mathcal{P}_i(\Lambda(\mathcal{P}_i^\dagger(\rho))).$$ The right Pauli super-operator acts on the input of the twirled channel and not the channel itself.

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
import numpy as np

Task 1 (7p)

Algebraically compute the resulting channel after applying twirling on:

1) Bit-flip channel $$ \Lambda_d (\rho) = (1-p) I \rho I + p (X \rho X)$$ 2) Coherent rotation $$ \Lambda_r (\rho) = \cos^2{\left(\frac{\theta}{2}\right)} I \rho I + \sin^2{\left(\frac{\theta}{2}\right)} X \rho X + \frac{i}{2} (\sin{(\theta)} I \rho X - \sin{(\theta)} X \rho I)$$ 3) Amplitude damping channel $$ \Lambda_a (\rho) = A_0 \rho A_0^\dagger + A_1 \rho A_1^\dagger$$ where $A_0 = \frac{(1+\sqrt{1-p})}{2} I + \frac{(1-\sqrt{1-p})}{2} Z$ and $A_1 = \frac{\sqrt{p}}{2} (X + iY)$.

Write the resulting channels as Pauli channels.

Task 2 (8p)

Write custom experiment and analysis classes using e.g. the qiskit-experiments framework for extracting density matrices after twirling noise channels using state tomography. The experiment class should accept at least one circuit with a noise channel as input and the analysis class should return a density matrix as output.

There are many possible ways of implementing this, but to give you an idea where to start, we provide you with a template below. Your experiment class could inherit from BaseExperiment or TomographyExperiment and implement the __init__ and circuit methods to produce the circuits necessary for Pauli twirling. Your analysis class could inherit from BaseAnalysis and implement the _run_analysis method to average the results over different twirls obtained from the experiment. You can use the instructions on the qiskit-experiments Documentation for help. Qiskit-experiments is fully open source, so looking at the source code can help as well.

from qiskit import QuantumCircuit
from typing import Tuple, List, Optional, Sequence
import matplotlib
from qiskit.providers.backend import Backend
from qiskit_experiments.library.tomography import TomographyExperiment
from qiskit_experiments.framework import (
    BaseAnalysis,
    Options,
    ExperimentData,
    AnalysisResultData
)

class TwirledStateTomographyAnalysis(BaseAnalysis):
    """Analysis class for state tomography with twirling."""

    @classmethod
    def _default_options(cls) -> Options:
        """Set default analysis options."""

        options = super()._default_options()
        options.plot = False
        options.ax = None
        return options

    def _run_analysis(
        self,
        experiment_data: ExperimentData
    ) -> Tuple[List[AnalysisResultData], List["matplotlib.figure.Figure"]]:
        """Run the analysis."""

        data_objects = experiment_data.data()

        # Combine twirled results and perform state tomography
        result = AnalysisResultData()

        return result, []

class TwirledStateTomography(TomographyExperiment):
    """State tomography with twirling."""

    def __init__(self,
                 circuit: QuantumCircuit,
                #  add your own arguments here...
                 physical_qubits: Sequence[int] = None,
                 measurement_indices: Sequence[int] = None,
                 backend: Optional[Backend] = None):
        """Initialize the experiment."""
        if physical_qubits is None:
            physical_qubits = tuple(range(circuit.num_qubits))
        super().__init__(circuit=circuit,
                        backend=backend,
                        physical_qubits=physical_qubits,
                        measurement_indices=measurement_indices,
                        analysis=TwirledStateTomographyAnalysis(),
                        )

    def circuits(self) -> List[QuantumCircuit]:
        """Generate the list of circuits to be run."""
        circuits = []
        # Generate circuits and populate metadata here

        return circuits

    @classmethod
    def _default_experiment_options(cls) -> Options:
        """Set default experiment options here."""
        options = super()._default_experiment_options()
        return options

Task 3 (8p)

Initialize a qubit in a state with non-zero populations and coherences. Apply each of the noise channels below to that qubit and perform state tomography with and without twirling applied on the channel. For each noise channel, plot on the same figure $\rho_{00}$, $\rho_{11}$, $Re(\rho_{01})$, $Im(\rho_{01})$ for the twirled and non-twirled cases for different values of $p$ or $\theta$.

Noise channels:

Bit-flip
- Pauli channel with $(1-p, p, 0, 0)$.
Coherent rotation
- Apply the RX gate with a small angle.
Amplitude damping channel

You may use the noise models implemented in qiskit_aer.noise. You are also allowed to use the implementations of the noise channels from previous projects, however beware that you may encounter numerical instabilities.

p_values = np.linspace(0, 1, 10)
theta_values = np.linspace(-np.pi, np.pi, 10)

Task 4 (7p)

Finally, let us invert or mitigate the twirled noise channels to extract the noiseless density matrix. After twirling, the noise channel should be in the form of a Pauli channel, which is relatively easy to invert if the noise coefficients are known. The inverted channel is not guaranteed to be a physical channel, as the coefficients may be negative. The procedure to invert common noise channels is described in S. Mangini, et al. EPJ Quantum Technol. 9, 29 (2022).

The inverse of a general Pauli channel with coefficients $(p_{0}, p_{x}, p_{y}, p_{z})$ is found to be $$\begin{aligned} & \Lambda^{-1}_{\boldsymbol{p}}(\Lambda) = \beta _{0} \Lambda + \beta _{1} \sigma _{x} {\Lambda} \sigma _{x} + \beta _{2} \sigma _{y} {\Lambda} \sigma _{y} + \beta _{3} \sigma _{z} {\Lambda} \sigma _{z}\quad \text{with} \\ & \beta _{0} =\frac{1}{4} \biggl(1+\frac{1}{1-2(p_{x}+p_{y})}+ \frac{1}{1-2(p_{x}+p_{z})}+\frac{1}{1-2(p_{y}+p_{z})} \biggr), \\ & \beta _{1} = \frac{1}{4} \biggl(1-\frac{1}{1-2(p_{x}+p_{y})}- \frac{1}{1-2(p_{x}+p_{z})}+\frac{1}{1-2(p_{y}+p_{z})} \biggr), \\ & \beta _{2} = \frac{1}{4} \biggl(1-\frac{1}{1-2(p_{x}+p_{y})}+ \frac{1}{1-2(p_{x}+p_{z})}-\frac{1}{1-2(p_{y}+p_{z})} \biggr), \\ & \beta _{3} = \frac{1}{4} \biggl(1+\frac{1}{1-2(p_{x}+p_{y})}- \frac{1}{1-2(p_{x}+p_{z})}-\frac{1}{1-2(p_{y}+p_{z})} \biggr) . \end{aligned}$$

Using the coefficients of the twirled noise channels from Task 1, invert each of the density matrices you found in Task 3 (for each $p$ and $\theta$) and compare $\rho_{00}$, $\rho_{11}$, $Re(\rho_{01})$, $Im(\rho_{01})$ to the non-inverted ones in e.g. a plot.

Hint: The expressions simplify a lot if you substitute the coefficients from Task 1.