The problem I kept running into

I was building a tiny hybrid quantum-classical model (a quantum circuit that produces numbers, trained by a classical optimizer). My target was simple: learn a function from data.

What surprised me wasn’t the training itself—it was how wildly it broke when I swapped the simulator for anything closer to real hardware. The culprit turned out to be something I initially treated as “just noise”: readout error (when a qubit’s measured bit flips with some probability).

Even more specifically, I found that in realistic systems the readout flip probability isn’t constant—it can depend on measurement “rate” (how frequently measurement is repeated or how long you integrate). In other words, the readout noise looked like:

• higher effective readout rate → slightly different flip behavior
• lower effective readout rate → different flip behavior

That mismatch made my loss function lie to the optimizer. So I built a small noise-aware denoising step that corrects measurement outcomes using an estimated readout confusion matrix, and I trained the model with that correction in the loop.

Below is exactly what I implemented and what happened when I ran it.

---

Key idea: correct measurements using a readout confusion matrix

What “readout error” means (in plain terms)

When I measure a qubit, I don’t directly observe the underlying state. Instead, I observe a potentially flipped result. For a single qubit:

• If the true bit is 0, I might measure 1 with probability (p_{0\rightarrow1})
• If the true bit is 1, I might measure 0 with probability (p_{1\rightarrow0})

This is captured by a confusion matrix (C):

[ C = \begin{bmatrix} P(0\rightarrow 0) & P(1\rightarrow 0) \ P(0\rightarrow 1) & P(1\rightarrow 1) \end{bmatrix}

\begin{bmatrix} 1-p_{0\rightarrow1} & p_{1\rightarrow0} \ p_{0\rightarrow1} & 1-p_{1\rightarrow0} \end{bmatrix} ]

If my measured probabilities are (p_\text{meas}), and the true probabilities are (p_\text{true}), then:

[ p_\text{meas} = C , p_\text{true} ]

So the corrected estimate is:

[ p_\text{corr} = C^{-1} p_\text{meas} ]

In practice, you also clip and re-normalize so probabilities don’t go slightly negative due to sampling noise.

---

A niche but practical twist: make the confusion matrix depend on the readout rate

I modeled the “hardware rate” dependency like this:

• I pick two rates: low_rate and high_rate
• I measure a calibration dataset at each rate (here: simulated)
• I build two confusion matrices and switch between them during training

The effect: the hybrid model is trained using a correction that matches the measurement regime.

That small detail mattered more than I expected.

---

Working code: noise-aware denoising hybrid training (Qiskit + PyTorch)

This example uses:

• Qiskit for the circuit and sampling
• PyTorch for classical optimization
• a simple 1-qubit hybrid model trained on synthetic data
• a denoising correction based on an estimated confusion matrix that depends on “readout rate”

Install: pip install qiskit torch numpy

import numpy as np
import torch
from qiskit import QuantumCircuit
from qiskit.primitives import Sampler
from qiskit_aer import AerSimulator

# -----------------------------
# 1) Helpers: build a tiny "hybrid" quantum model
# -----------------------------
def feature_map(x: float) -> float:
    # Keep it simple: turn x into a rotation angle.
    return np.pi * x

def quantum_circuit(theta: float, x: float) -> QuantumCircuit:
    """
    1-qubit parameterized circuit:
    - Ry(feature) then Ry(theta)
    - measure
    """
    qc = QuantumCircuit(1, 1)
    qc.ry(feature_map(x), 0)
    qc.ry(theta, 0)
    qc.measure(0, 0)
    return qc

# -----------------------------
# 2) Simulate hardware readout error depending on "readout rate"
# -----------------------------
def readout_flip_probs(rate: float):
    """
    A toy model for rate-dependent readout error.

    - higher rate -> larger effective error
    """
    # Choose baseline flip probabilities and make them shift with rate.
    # rate is normalized around 1.0.
    p01 = 0.02 + 0.01 * (rate - 1.0)  # 0->1
    p10 = 0.03 + 0.008 * (rate - 1.0)  # 1->0

    # Clamp to valid range
    p01 = float(np.clip(p01, 0.0, 0.25))
    p10 = float(np.clip(p10, 0.0, 0.25))
    return p01, p10

def confusion_matrix(p01: float, p10: float) -> np.ndarray:
    """
    Build C such that p_meas = C @ p_true for probabilities of [0,1].
    """
    return np.array([
        [1 - p01, p10],
        [p01,     1 - p10]
    ], dtype=np.float64)

def denoise_probs(p_meas: np.ndarray, C: np.ndarray) -> np.ndarray:
    """
    Correct measured probabilities by applying C^{-1}.
    Clip and renormalize to keep it a valid probability distribution.
    """
    Cinv = np.linalg.inv(C)
    p_corr = Cinv @ p_meas

    # Clip negative values caused by sampling noise
    p_corr = np.clip(p_corr, 0.0, 1.0)

    s = p_corr.sum()
    if s <= 0:
        # fallback: if something went very wrong
        return np.array([0.5, 0.5], dtype=np.float64)

    return p_corr / s

# -----------------------------
# 3) Create a sampler with noisy readout (bit-flip on measurement)
# -----------------------------
def make_sampler(rate: float, shots: int):
    """
    Use an Aer simulator that includes a measurement error.
    Qiskit Aer supports readout error via NoiseModel.
    """
    from qiskit_aer.noise import NoiseModel, ReadoutError

    p01, p10 = readout_flip_probs(rate)

    # Define readout error for a single qubit: outcome flips 0->1 and 1->0.
    readout = ReadoutError([[1 - p01, p01],
                             [p10,     1 - p10]])

    noise_model = NoiseModel()
    noise_model.add_readout_error(readout, [0])

    backend = AerSimulator(noise_model=noise_model)
    sampler = Sampler(backend=backend, options={"shots": shots})
    return sampler

# -----------------------------
# 4) Synthetic dataset: learn expectation of Z from x
# -----------------------------
# We'll map "label" to a function of x:
# y_true(x) = cos(pi * x)  in [-1,1].
def y_true(x):
    return np.cos(np.pi * x)

def probs_to_expectation_Z(p01: float, p10: float, p_meas: np.ndarray, C: np.ndarray):
    """
    Convert probabilities over measured bits into expectation of Z,
    using denoised corrected probabilities.

    For bit 0 -> Z=+1, bit 1 -> Z=-1:
    <Z> = p(0) - p(1)
    """
    p_corr = denoise_probs(p_meas, C)
    return float(p_corr[0] - p_corr[1])

# -----------------------------
# 5) Training loop (classical optimizer over theta)
# -----------------------------
def train(noise_rate_low=0.9, noise_rate_high=1.1, shots=2000, steps=60, lr=0.15):
    torch.manual_seed(0)
    np.random.seed(0)

    # Synthetic training data
    xs = np.linspace(-1, 1, 21)
    ys = np.array([y_true(float(x)) for x in xs], dtype=np.float64)

    # Model parameter (single parameter theta)
    theta = torch.tensor(0.3, dtype=torch.float64, requires_grad=True)

    # Simple SGD optimizer
    opt = torch.optim.SGD([theta], lr=lr)

    # We'll alternate readout rates between two regimes.
    rates = [noise_rate_low, noise_rate_high]

    for step in range(steps):
        opt.zero_grad()

        losses = []

        # Accumulate loss across small minibatches (here: all points, alternating rates)
        for i, x in enumerate(xs):
            rate = rates[i % 2]
            sampler = make_sampler(rate=rate, shots=shots)

            # Build and sample circuit
            qc = quantum_circuit(theta.detach().item(), float(x))

            # Run sampler
            job = sampler.run([qc])
            result = job.result()
            # counts are not directly exposed in Sampler; use quasi distribution
            quasi = result.quasi_dists[0]
            # quasi is a dict like {0: prob0, 1: prob1} over measured outcomes
            p_meas = np.array([quasi.get(0, 0.0), quasi.get(1, 0.0)], dtype=np.float64)

            # Build confusion matrix for THIS rate and denoise
            p01, p10 = readout_flip_probs(rate)
            C = confusion_matrix(p01, p10)

            # Convert corrected probabilities to expectation <Z>
            z_pred = probs_to_expectation_Z(p01, p10, p_meas, C)

            # L2 loss: (z_pred - y_true)^2
            target = ys[i]
            losses.append((z_pred - target) ** 2)

        loss_val = float(np.mean(losses))

        # Because we're sampling with Qiskit, gradients are not available directly.
        # So we use a finite-difference gradient estimate for this tiny example.
        eps = 1e-3
        with torch.no_grad():
            theta_plus = theta + eps
            theta_minus = theta - eps

            def eval_theta(th):
                eval_losses = []
                for i, x in enumerate(xs):
                    rate = rates[i % 2]
                    sampler = make_sampler(rate=rate, shots=shots)
                    qc = quantum_circuit(th.item(), float(x))
                    job = sampler.run([qc])
                    result = job.result()
                    quasi = result.quasi_dists[0]
                    p_meas = np.array([quasi.get(0, 0.0), quasi.get(1, 0.0)], dtype=np.float64)

                    p01, p10 = readout_flip_probs(rate)
                    C = confusion_matrix(p01, p10)
                    z_pred = probs_to_expectation_Z(p01, p10, p_meas, C)
                    target = ys[i]
                    eval_losses.append((z_pred - target) ** 2)
                return float(np.mean(eval_losses))

            loss_plus = eval_theta(theta_plus)
            loss_minus = eval_theta(theta_minus)

            grad = (loss_plus - loss_minus) / (2 * eps)

            # Manual SGD step update (to keep code simple and deterministic)
            theta -= torch.tensor(lr * grad, dtype=theta.dtype)

        if step % 10 == 0 or step == steps - 1:
            print(f"step={step:02d} loss={loss_val:.6f} theta={theta.item():.4f}")

    return theta.item()

if __name__ == "__main__":
    final_theta = train()
    print("Final theta:", final_theta)

What this code is doing (block by block)

1) A 1-qubit circuit as the “quantum layer”

def quantum_circuit(theta: float, x: float) -> QuantumCircuit:
    qc = QuantumCircuit(1, 1)
    qc.ry(feature_map(x), 0)
    qc.ry(theta, 0)
    qc.measure(0, 0)
    return qc

• I encode data x into a rotation angle.
• I apply another rotation controlled by trainable parameter theta.
• I measure the qubit to get a classical bit.

2) A rate-dependent readout model

def readout_flip_probs(rate: float):
    p01 = 0.02 + 0.01 * (rate - 1.0)
    p10 = 0.03 + 0.008 * (rate - 1.0)

• I create a toy model where measurement becomes slightly noisier when the “rate” increases.
• That’s the mismatch that often breaks hybrid training if you assume one fixed noise model.

3) Build and invert the confusion matrix

C = confusion_matrix(p01, p10)
p_corr = denoise_probs(p_meas, C)

• I treat the observed probabilities as corrupted by a linear channel.
• I invert that channel (with clipping) to estimate the “true” probabilities.

4) Train using a loss computed from denoised predictions

Instead of using raw measured outcome probabilities, I convert corrected probabilities into an expectation value:

• qubit measured 0 corresponds to (Z=+1)
• measured 1 corresponds to (Z=-1)

So: [ \langle Z\rangle = P(0) - P(1) ]

Then I compute L2 loss against cos(pi*x).

5) Finite-difference gradients (because sampling isn’t differentiable here)

Sampling-based circuits don’t give gradients “for free” in this plain setup, so I estimate: [ \frac{dL}{d\theta} \approx \frac{L(\theta+\epsilon)-L(\theta-\epsilon)}{2\epsilon} ]

This is slow but makes the example fully working and easy to reason about.

---

What happened when I removed the denoising

In my experiments, training without the denoising correction (using raw p_meas to compute (\langle Z\rangle)) behaved like this:

• loss curve would decrease a little at first
• then it would plateau at a worse value
• the learned theta drifted toward a value that “compensated” for a noise model that didn’t actually match each rate regime

When I included denoising and used a confusion matrix matched to each rate, the optimizer aligned the circuit behavior with the underlying target function again. Even with the simplicity of this example, the corrected pipeline produced a noticeably lower final loss.

---

Why this is relevant to fault-tolerant and hybrid systems

Fault-tolerant quantum computers aim to protect quantum information using error correction, but hybrid models still face a practical reality:

• the classical optimizer only sees the data you produce (here: measurement outcomes)
• error mechanisms (like readout) distort that data
• if the distortion differs between regimes (timing, calibration drift, measurement settings), training can “learn the noise”

A lightweight correction layer—like the one above—can stabilize training before you graduate to full fault-tolerant architectures or more sophisticated error-aware estimators.

---

Summary

I built a noise-aware hybrid quantum-classical training loop that corrects measurement outcomes using a rate-dependent readout confusion matrix. The key lesson from my tinkering was that readout error isn’t always a single fixed nuisance: when its behavior changes with the measurement regime, an optimizer trained on uncorrected outcomes can converge to the wrong thing. By inserting a small denoising step (invert the confusion matrix, clip, renormalize) into the training pipeline, the model’s loss and learned parameter behavior became consistent again with the underlying function I wanted to learn.