Noise Simulation in PyQuil

Real quantum hardware is noisy. Gates are imperfect, qubits lose coherence over time, and measurements introduce errors. To write useful NISQ algorithms, you need to understand how noise affects your circuits and how to simulate that noise before running on hardware. PyQuil provides a flexible noise simulation framework built around Kraus operators and noise models that you can apply to the QVM.

Prerequisites

You should be comfortable writing basic PyQuil programs and running them on the QVM. Familiarity with density matrices and quantum channels is helpful but not strictly required.

pip install pyquil numpy matplotlib

You will also need the QVM and Quilc compiler running locally (from Rigetti’s Forest SDK), or access to QCS.

Quantum Noise and Kraus Operators

Noise in quantum computing is modeled using quantum channels. A quantum channel transforms a density matrix according to a set of Kraus operators {K_i}, where the output state is:

rho' = sum_i K_i @ rho @ K_i^dagger

The Kraus operators must satisfy the completeness relation: sum_i K_i^dagger @ K_i = I. This ensures the channel is trace-preserving, meaning probabilities still sum to one.

PyQuil lets you define custom Kraus operators for any gate and attach them to a noise model that the QVM uses during simulation.

Building a Depolarizing Noise Model

A depolarizing channel is one of the simplest noise models. For a single qubit, it replaces the state with the maximally mixed state (I/2) with probability p, and leaves it unchanged with probability 1-p. The Kraus operators for single-qubit depolarizing noise are:

import numpy as np
from pyquil.noise import KrausModel, NoiseModel
from pyquil import Program, get_qc
from pyquil.gates import H, CNOT, MEASURE, I as I_gate

def depolarizing_kraus_operators(p):
    """
    Return Kraus operators for single-qubit depolarizing noise
    with error probability p.
    """
    # Identity (no error) with probability 1 - 3p/4
    K0 = np.sqrt(1 - 3 * p / 4) * np.eye(2)
    # Bit flip (X error) with probability p/4
    K1 = np.sqrt(p / 4) * np.array([[0, 1], [1, 0]])
    # Phase flip (Z error) with probability p/4
    K2 = np.sqrt(p / 4) * np.array([[1, 0], [0, -1]])
    # Bit-phase flip (Y error) with probability p/4
    K3 = np.sqrt(p / 4) * np.array([[0, -1j], [1j, 0]])
    return [K0, K1, K2, K3]

# Create Kraus operators with 5% depolarizing error
p_error = 0.05
kraus_ops = depolarizing_kraus_operators(p_error)

# Verify completeness: sum of K_i^dagger @ K_i should equal identity
completeness = sum(k.conj().T @ k for k in kraus_ops)
print("Completeness check (should be identity):")
print(np.round(completeness, 6))

Expected Output

Completeness check (should be identity):
[[1.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

Applying Noise to a Bell State Circuit

Now let’s build a noise model and attach it to the QVM. We will create a Bell state circuit and compare ideal versus noisy results.

from pyquil.noise import KrausModel, NoiseModel

# Define the Kraus model for H gate on qubit 0
h_kraus = KrausModel(
    gate="H",
    params=[],
    targets=(0,),
    kraus_ops=depolarizing_kraus_operators(0.05),
    fidelity=0.95,
)

# Define the Kraus model for CNOT gate on qubits 0, 1
# For a two-qubit gate, we need 4x4 Kraus operators
def two_qubit_depolarizing_kraus(p):
    """Two-qubit depolarizing channel Kraus operators."""
    d = 4  # dimension of two-qubit Hilbert space
    K0 = np.sqrt(1 - p) * np.eye(d)
    # Generate the remaining 15 Pauli tensor products
    paulis_1q = [
        np.eye(2),
        np.array([[0, 1], [1, 0]]),
        np.array([[0, -1j], [1j, 0]]),
        np.array([[1, 0], [0, -1]]),
    ]
    kraus_list = [K0]
    for i in range(4):
        for j in range(4):
            if i == 0 and j == 0:
                continue  # skip identity (already included)
            K = np.sqrt(p / 15) * np.kron(paulis_1q[i], paulis_1q[j])
            kraus_list.append(K)
    return kraus_list

cnot_kraus = KrausModel(
    gate="CNOT",
    params=[],
    targets=(0, 1),
    kraus_ops=two_qubit_depolarizing_kraus(0.03),
    fidelity=0.97,
)

# Build the noise model
noise_model = NoiseModel(
    gates=[h_kraus, cnot_kraus],
    assignment_probs={0: [0.99, 0.01], 1: [0.01, 0.99]},
    # assignment_probs models readout error:
    # P(measure 0 | state 0) = 0.99, P(measure 1 | state 0) = 0.01
)

Now run the Bell state circuit with and without noise.

from collections import Counter

# Build Bell state program
p = Program()
ro = p.declare('ro', 'BIT', 2)
p += H(0)
p += CNOT(0, 1)
p += MEASURE(0, ro[0])
p += MEASURE(1, ro[1])
p.wrap_in_numshots_loop(10000)

# Ideal simulation (no noise)
qc_ideal = get_qc('2q-qvm')
result_ideal = qc_ideal.run(qc_ideal.compile(p))
counts_ideal = Counter(map(tuple, result_ideal.readout_data['ro'].tolist()))
print("Ideal counts:", dict(counts_ideal))

# Noisy simulation
qc_noisy = get_qc('2q-qvm', noise_model=noise_model)
result_noisy = qc_noisy.run(qc_noisy.compile(p))
counts_noisy = Counter(map(tuple, result_noisy.readout_data['ro'].tolist()))
print("Noisy counts:", dict(counts_noisy))

Expected Output

Ideal counts: {(0, 0): 5021, (1, 1): 4979}
Noisy counts: {(0, 0): 4712, (1, 1): 4680, (0, 1): 298, (1, 0): 310}

In the ideal case, you only see correlated outcomes |00> and |11>. With noise, some fraction of shots produce the “wrong” outcomes |01> and |10>, which would never appear in a perfect Bell state. This is the signature of decoherence and gate errors breaking entanglement.

Adding T1/T2 Decoherence

T1 (relaxation) and T2 (dephasing) are the most important decoherence timescales for real qubits. T1 describes how quickly an excited qubit decays to the ground state. T2 describes how quickly phase information is lost. PyQuil provides built-in support for adding these effects.

from pyquil.noise import add_decoherence_noise

# Build a simple program
p = Program()
ro = p.declare('ro', 'BIT', 1)
p += H(0)
p += MEASURE(0, ro[0])
p.wrap_in_numshots_loop(10000)

# Add T1/T2 decoherence noise
# T1 = 30 microseconds, T2 = 15 microseconds
# gate_time = 60 nanoseconds (typical single-qubit gate time)
noisy_program = add_decoherence_noise(
    p,
    T1=30e-6,
    T2=15e-6,
    gate_time_1q=60e-9,
    gate_time_2q=150e-9,
)

print("Program with decoherence noise:")
print(noisy_program)

The add_decoherence_noise function inserts Kraus operators after each gate that model amplitude damping (T1) and phase damping (T2) for the specified gate durations. This gives you a physically motivated noise model without manually constructing Kraus operators.

Comparing Noise Levels

A useful exercise is sweeping the noise parameter and observing how the fidelity of your Bell state degrades.

import matplotlib.pyplot as plt

error_rates = np.linspace(0, 0.3, 20)
fidelities = []

for p_err in error_rates:
    h_km = KrausModel(
        gate="H", params=[], targets=(0,),
        kraus_ops=depolarizing_kraus_operators(p_err),
        fidelity=1.0 - p_err,
    )
    cnot_km = KrausModel(
        gate="CNOT", params=[], targets=(0, 1),
        kraus_ops=two_qubit_depolarizing_kraus(p_err),
        fidelity=1.0 - p_err,
    )
    nm = NoiseModel(gates=[h_km, cnot_km], assignment_probs={})

    prog = Program()
    ro_reg = prog.declare('ro', 'BIT', 2)
    prog += H(0)
    prog += CNOT(0, 1)
    prog += MEASURE(0, ro_reg[0])
    prog += MEASURE(1, ro_reg[1])
    prog.wrap_in_numshots_loop(5000)

    qc = get_qc('2q-qvm', noise_model=nm)
    result = qc.run(qc.compile(prog))
    counts = Counter(map(tuple, result.readout_data['ro'].tolist()))

    # Bell state fidelity: fraction of correlated outcomes
    correlated = counts.get((0, 0), 0) + counts.get((1, 1), 0)
    fidelity = correlated / 5000
    fidelities.append(fidelity)

plt.plot(error_rates, fidelities, 'o-')
plt.xlabel('Depolarizing Error Rate')
plt.ylabel('Bell State Fidelity')
plt.title('Bell State Fidelity vs. Noise')
plt.grid(True)
plt.savefig('noise_sweep.png')
plt.show()

You should see a smooth decline from fidelity near 1.0 at zero noise down to roughly 0.5 at high noise levels, where the output is essentially random.

Practical Implications for NISQ Algorithms

Understanding noise simulation matters for several reasons when working with near-term quantum devices:

1. Circuit depth limits. Every gate adds noise. Deeper circuits accumulate more errors, so NISQ algorithms must keep circuit depth as short as possible. Simulating noise helps you estimate the maximum useful depth for a given hardware error rate.

2. Error mitigation strategies. Techniques like zero-noise extrapolation (ZNE) and probabilistic error cancellation (PEC) require understanding the noise model. You can use PyQuil’s noise simulation to test these mitigation techniques before deploying on hardware.

3. Algorithm selection. Some algorithms are more noise-resilient than others. Variational algorithms (VQE, QAOA) tend to be more robust because their classical optimizer can partially compensate for noise, while algorithms requiring deep coherent circuits (like QPE) degrade quickly.

4. Hardware benchmarking. By matching your noise model parameters to actual device calibration data (T1, T2, gate fidelities, readout errors), you can predict how your program will perform on a specific QPU before using your allocation.

Summary

PyQuil’s noise simulation tools let you model realistic quantum errors using Kraus operators and noise models. You can construct custom depolarizing channels, add physically motivated T1/T2 decoherence, include readout errors, and sweep noise parameters to understand how your algorithms degrade. This is essential for developing practical NISQ applications, where understanding and mitigating noise is just as important as designing the quantum circuit itself.