add md file

docs/source/examples/examples.md

@@ -29,6 +29,7 @@ ZNE and CDR with Cirq: 1D Ising Simulation <quantum_simulation_1d_ising.md>
 ZNE with PennyLane + Cirq: Energy of molecular Hydrogen <molecular_hydrogen_pennylane.md>
 ZNE with BQSKit compiled circuits <bqskit.md>
 ZNE on Stim backend with Cirq: Logical randomized benchmarking circuits <zne_logical_rb_cirq_stim.md>
+ZNE and PT with Cirq: ZNE with Noise Tailoring <pt_zne.md>
 LRE vs ZNE: comparing performance and overhead <lre-zne-comparison.md>
 PEC on a Braket simulator: Mirror circuits <pec_tutorial.md>
 PEC with Cirq: Learning representations <learning-depolarizing-noise.md>

docs/source/examples/pt_zne.md

@@ -0,0 +1,335 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.1
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+```{tags} cirq, zne, pt, basic
+```
+# ZNE with Noise Tailoring
+
+This tutorial helps us understand the cases when noise tailoring is helpful for improving the performance of 
+quantum error mitigation techniques. 
+
+We will
+
+- compare coherent noise to incoherent noise
+- Pauli twirl coherent noise to be incoherent
+- compare the performance of [Zero-Noise Extrapolation](../guide/zne.md) (ZNE)
+on its own and in combination with [Pauli Twirling](../guide/pt.md) (PT) when the circuit is subjected to coherent noise
+
+## Coherent noise vs. Incoherent noise
+
+**Coherent noise** is a reversible process as long as the noisy unitary transformation is known beforehand which is not always the case. These types of noise maintain the purity of the state. But in a quantum circuit, they are easily carried across the circuit such that they scale up quadratically. Dealing with coherent noise requires a large resource overhead.  
+
+**Incoherent noise** is a process that entangles the quantum system with its environment i.e. this type of noise is irreversible. The system and the environment are in a mixed state. It scales linearly in the small error limit. 
+The noise channel can be described using Pauli operators which makes it easy to analyze and simulate. Worst case error rate is directly proportional to the **average gate infidelity**. 
+
+```{note}
+If $\mathcal{F}$ is the fidelity defining the success of preparing an arbitrary pure state $\rho$, then
+$1-\mathcal{F}$ is the **average gate infidelity**.
+```
+
+For example, a depolarizing noise channel is a stochastic noise channel where a noiseless process is probabilistically mixed with orthogonal errors. If $\rho$ is a single qubit state, $p$ is the probabilistic error rate and $\epsilon(\rho)$ is the noise channel:
+
+```{math}
+\epsilon(\rho) = (1-p) \rho + p \frac{I}{2}
+```
+$\frac{I}{2}$ is the maximally mixed state which can be described using Paulis. 
+
+```{math}
+\frac{I}{2} = \frac{1}{4} (\rho + X \rho X + Y \rho Y + Z \rho Z)
+```
+
+Thus, the depolarizing channel can be redescribed using Paulis as shown below. 
+
+```{math}
+\epsilon(\rho) = (1-\frac{3p}{4}) \rho + \frac{p}{4} (X \rho X + Y \rho Y + Z \rho Z)
+
+```
+
+### Pauli Transfer Matrix
+
+Let $\Lambda(rho)$ be a $n-$qubit noise channel where $K_i$ are the Kraus operators.
+
+```{math}
+
+\Lambda(\rho) = \sum_{i=1} K_i \rho {K_i}^\dagger
+```
+If $P_i$ and $P_j$ are $n-$qubit Paulis, the following expression defines the entries of a Pauli Transfer Matrix (PTM). Here, 
+$i$ defines the rows while $j$ defines the columns of the PTM. All entries of the PTM are real and in the interval $[-1, 1]$.
+
+
+```{math}
+(R_{\Lambda})_{ij} = \frac{1}{2^n} \text{Tr} \{ P_i \Lambda(P_j)\}
+```
+A PTM allows us to distinguish between the two types of noise. The off-diagonal terms of the PTM are due to the effect of coherent noise while the diagonal terms are due to incoherent noise. PTM of a circuit is a product of the PTM of each layer in the circuit. Due to these, it is straightforward to see how coherent noise carries across different layers in the circuit and how incoherent errors are easier to deal with in the small error limit. 
+
+The known fault tolerant thresholds for stochastic noise are higher than coherent noise which makes the former a 'preferable' type of noise compared to the latter. If we want to not deal with coherent noise, Pauli twirling can be used to tailor coherent noise to incoherent noise. Same as the expression above, when a coherent noise channel is Pauli twirled, the noise channel can be described using Paulis after averaging over multiple Pauli twirled circuits. Refer to the [Pauli Twirling user guide](../guide/pt.md) for additional information. 
+
+## Using Pauli Twirling in Mitiq
+
+```{code-cell} ipython3
+# Utility functions for Pauli Twirling.
+
+import cirq
+import numpy as np
+import numpy.typing as npt
+from cirq.circuits import Circuit
+
+from mitiq.pec.channels import _circuit_to_choi, choi_to_super
+from mitiq.utils import matrix_to_vector, vector_to_matrix
+
+pauli_unitary_list = [
+    cirq.unitary((cirq.I)),
+    cirq.unitary((cirq.X)),
+    cirq.unitary((cirq.Y)),
+    cirq.unitary((cirq.Z)),
+]
+
+
+def _n_qubit_paulis(num_qubits: int) -> list[npt.NDArray[np.complex64]]:
+    """Get a list of n-qubit Pauli unitaries."""
+    if not num_qubits >= 1:
+        raise ValueError("Invalid number of qubits provided.")
+
+    # get the n-qubit paulis from the Pauli group
+    # disregard the n-qubit paulis with complex phase
+
+    n_qubit_paulis = pauli_unitary_list
+    for i in range(num_qubits - 1):
+        n_qubit_paulis = n_qubit_paulis
+        empty_pauli_list = []
+        for j in pauli_unitary_list:
+            for k in n_qubit_paulis:
+                new_pauli = np.kron(j, k)
+                empty_pauli_list.append(new_pauli)
+        n_qubit_paulis = empty_pauli_list
+    return n_qubit_paulis
+
+
+def _pauli_vectorized_list(num_qubits: int) -> list[npt.NDArray[np.complex64]]:
+    """Define a function to create a list of vectorized matrices.
+
+    If the density matrix of interest has more than n>1 qubits, the
+    Pauli group is used to generate n-fold tensor products before
+    vectorizing the unitaries.
+    """
+    n_qubit_paulis = _n_qubit_paulis(num_qubits)
+    output_pauli_vec_list = []
+    for i in n_qubit_paulis:
+        # the matrix_to_vector function stacks rows in vec form
+        # transpose is used here to instead stack the columns
+        matrix_trans = np.transpose(i)
+        output_pauli_vec_list.append(matrix_to_vector(matrix_trans))
+    return output_pauli_vec_list
+
+
+def ptm_matrix(circuit: Circuit, num_qubits: int) -> npt.NDArray[np.complex64]:
+    """Find the Pauli Transfer Matrix (PTM) of a circuit."""
+
+    superop = choi_to_super(_circuit_to_choi(circuit))
+
+    vec_pauli = _pauli_vectorized_list(num_qubits)
+    n_qubit_paulis = _n_qubit_paulis(num_qubits)
+    ptm_matrix = np.zeros([4**num_qubits, 4**num_qubits], dtype=complex)
+
+    for i in range(len(vec_pauli)):
+        superop_on_pauli_vec = np.matmul(superop, vec_pauli[i])
+        superop_on_pauli_matrix_transpose = vector_to_matrix(
+            superop_on_pauli_vec
+        )
+        superop_on_pauli_matrix = np.transpose(
+            superop_on_pauli_matrix_transpose
+        )
+
+        # ptm_matrix_row = []
+        for j in range(len(n_qubit_paulis)):
+            pauli_superop_pauli = np.matmul(
+                n_qubit_paulis[j], superop_on_pauli_matrix
+            )
+            ptm_matrix[j, i] = (0.5**num_qubits) * np.trace(
+                pauli_superop_pauli
+            )
+
+    return ptm_matrix
+```
+
+Let us consider a simple circuit of a CNOT gate. We are going to subject this circuit to two types of noise and compare
+their respective PTMs.
+
+```{code-cell} ipython3
+import numpy as np
+import seaborn as sns
+import matplotlib.pylab as plt
+
+from cirq import LineQubit, Circuit, CNOT, Ry, depolarize, X, Y, Z
+
+q0 = LineQubit(0)
+q1 = LineQubit(1)
+
+circuit = Circuit(CNOT(q0, q1))
+print(circuit)
+```
+```{code-cell} ipython3
+ptmcnot = ptm_matrix(circuit, 2)
+ax = sns.heatmap(ptmcnot.real, linewidth=0.5)
+print("Ideal CNOT PTM")
+plt.show()
+```
+```{code-cell} ipython3
+# PTM of a noisy CNOT gate: depolarizing noise
+noisy_circuit_incoherent = circuit.with_noise(depolarize(p=0.3))
+print(noisy_circuit_incoherent)
+
+ptmcnot = ptm_matrix(noisy_circuit_incoherent, 2)
+ax = sns.heatmap(ptmcnot.real, linewidth=0.5)
+print("\n Pauli Transfer Matrix of noisy CNOT (incoherent)")
+plt.show()
+```
+```{code-cell} ipython3
+# PTM of a noisy CNOT gate: Rz
+noisy_circuit_coherent = circuit.with_noise(Ry(rads=np.pi/12))
+print(noisy_circuit_coherent)
+
+ptmcnot = ptm_matrix(noisy_circuit_coherent, 2)
+ax = sns.heatmap(ptmcnot.real, linewidth=0.5)
+print("Pauli Transfer Matrix of noisy CNOT (coherent)")
+plt.show()
+```
+If we compare the PTM of the ideal CNOT gate to those when the gate was subjected to incoherent noise and coherent noise, 
+there are additional sources of errors to deal with when coherent noise is acting on the CNOT gate. These can be reduced or tailored to be close to how the incohrent noise PTM appears through Pauli Twirling.
+
+```{code-cell} ipython3
+from mitiq.pt import generate_pauli_twirl_variants
+
+# Generate twirled circuits
+NUM_TWIRLED_VARIANTS = 3
+twirled_circuits = generate_pauli_twirl_variants(
+    circuit, num_circuits=NUM_TWIRLED_VARIANTS)
+print("Example ideal twirled circuit", twirled_circuits[0], sep="\n")
+```
+Now, lets add coherent noise to the CNOT gate in each twirled circuit.
+```{code-cell} ipython3
+
+noisy_twirled_circuits = []
+
+for circ in twirled_circuits:
+    split_circuit = Circuit(circ[0], circ[1], Ry(rads=np.pi/12)(q0), Ry(rads=np.pi/12)(q1), circ[-1])
+    noisy_twirled_circuits.append(split_circuit)
+
+print("Example noisy twirled circuit", noisy_twirled_circuits[0], sep="\n")
+```
+
+The twirled PTM is averaged over each noisy twirled circuit such that the new PTM is close to that of the PTM of incoherent noise. We skip the step in this section as we require a very large number of twirled circuits to demonstrate the desired effect of averaging over multiple numpy arrays.
+
+## Noisy ZNE
+
+Lets define a larger circuit of CNOT, CZ and H gates. 
+
+```{code-cell} ipython3
+
+from cirq import LineQubit, Circuit, CZ, CNOT, H
+
+q0, q1, q2, q3  = LineQubit.range(4)
+circuit = Circuit(
+    H(q0),
+    CNOT.on(q0, q1),
+    CZ.on(q1, q2),
+    CNOT.on(q2, q3),
+)
+
+print(circuit)
+
+```
+
+We are going to add coherent noise to this circuit and then get the error-mitigated expectation value. For a detailed discussion on this, refer to the [ZNE user guide](../guide/zne-1-intro.md). 
+
+As we are using a simulator, we have to make sure the twirling gates are sandwiching the coherent noise channel in the circuit.
+For this, `get_noise_model` is used to add noise to CZ/CNOT gates. See [PT user guide](../guide/pt-1-intro.md) for more. 
+
+```{code-cell} ipython3
+from numpy import pi
+from cirq import CircuitOperation, CXPowGate, CZPowGate, DensityMatrixSimulator
+from cirq.devices.noise_model import GateSubstitutionNoiseModel
+
+def get_noise_model(noise_level: float) -> GateSubstitutionNoiseModel:
+    """Substitute each CZ and CNOT gate in the circuit
+    with the gate itself followed by an Rx rotation on the output qubits.
+    """
+    rads = pi / 2 * noise_level
+    def noisy_c_gate(op):
+        if isinstance(op.gate, (CZPowGate, CXPowGate)):
+            return CircuitOperation(
+                Circuit(
+                    op.gate.on(*op.qubits), 
+                    Ry(rads=rads).on_each(op.qubits),
+                ).freeze())
+        return op
+
+    return GateSubstitutionNoiseModel(noisy_c_gate)
+
+def execute(circuit: Circuit, noise_level: float):
+    """Returns Tr[ρ |0⟩⟨0|] where ρ is the state prepared by the circuit."""
+    return (
+        DensityMatrixSimulator(noise=get_noise_model(noise_level=noise_level))
+        .simulate(circuit)
+        .final_density_matrix[0, 0]
+        .real
+    )
+
+
+# Set the intensity of the noise
+NOISE_LEVEL = 7
+
+
+# Compute the expectation value of the |0><0| observable
+# in both the noiseless and the noisy setup
+ideal_value = execute(circuit, noise_level=0.0)
+noisy_value = execute(circuit, noise_level=NOISE_LEVEL)
+
+NUM_TWIRLED_VARIANTS = 300
+twirled_circuits = generate_pauli_twirl_variants(circuit, num_circuits=NUM_TWIRLED_VARIANTS)
+
+# Average results executed over twirled circuits
+from functools import partial
+from mitiq import Executor
+pt_vals = Executor(partial(execute, noise_level=NOISE_LEVEL)).evaluate(twirled_circuits)
+twirled_result = np.average(pt_vals)
+
+
+print(f"Error without twirling: {abs(ideal_value - noisy_value) :.3}")
+print(f"Error with twirling: {abs(ideal_value - twirled_result) :.3}")
+```
+
+## Combining Pauli Twirling with ZNE
+
+```{code-cell} ipython3
+from mitiq.zne import execute_with_zne
+
+executor=partial(execute, noise_level=NOISE_LEVEL)
+zne_pt_vals = []
+
+for i in twirled_circuits:
+    zne_pt_vals.append(execute_with_zne(i, executor))
+
+mitigated_result = np.average(zne_pt_vals)
+
+print(f"Error without twirling: {abs(ideal_value - noisy_value) :.3}")
+print(f"Error with ideal and twirling: {abs(ideal_value - twirled_result) :.3}")
+print(f"Error with ideal and ZNE + PT: {abs(ideal_value - mitigated_result) :.3}")
+
+```
+Accordingly, just PT or ZNE do not work that well compared to a combination of PT and ZNE.
+
+## Conclusion
+
+In this tutorial, we've shown how to use a noise tailoring method with Zero-Noise Extrapolation.
+If you're interested in finding out more about these techniques, check out their respective sections of the users guide: [ZNE](../guide/zne.md), [Pauli Twilring](../guide/pt.md).