Default float32 in TensorNetworkDecoder returns non-physical posteriors (negative coset masses) on degenerate syndromes

[killianadrover]

Default float32 in TensorNetworkDecoder returns non-physical posteriors (negative coset masses) on degenerate syndromes

Summary

TensorNetworkDecoder is a maximum-likelihood decoder that performs an exact (untruncated)
tensor-network contraction, but its constructor defaults to single precision while the docstring
states double precision. In float32, the contraction loses significance on rare high-weight /
low-p syndromes and returns negative coset masses, so decode() produces posterior
probabilities outside [0, 1] (observed from -3.08 to +6.89) — passed through with no
validity check. In float64 the same inputs are always valid.

Environment

• cudaq-qec (cudaq-qec-cu12) 0.6.0
• Python 3.12 ; backend: cutensornet on CUDA (also reproduced with the CPU/torch contractor)
• Code: rotated surface code, circuit-level SI1000 noise, exact (sliced, untruncated) contraction.

Documentation vs. actual default

In libs/qec/python/cudaq_qec/plugins/decoders/tensor_network_decoder.py,
TensorNetworkDecoder.__init__ (0.6.0 wheel: default line 111 / docstring line 131;
on main: 113 / 133):

• signature: dtype: str = "float32"
• docstring: dtype (str, optional): ... Defaults to "float64".

Two more places corroborate that float32 is the de facto default and is known to be imprecise:

• your own test pins it — libs/qec/python/tests/test_tensor_network_decoder.py:65:
  assert decoder._dtype == "float32" (decoder built with no dtype argument);
• .../tensor_network_utils/nm_optimizer.py:38-43:
  # Coarse for fp32 because ``1.0 - 1e-12`` rounds back to ``1.0``.
_PRIOR_EPS_BY_DTYPE = {"float64": 1e-12, "float32": 1e-6}.

So a user relying on the documented default silently runs in float32.

Minimal reproduction (the contradiction itself)

import cudaq_qec as qec, numpy as np
H = np.array([[1, 1, 0], [0, 1, 1]], dtype=np.float64)
A = np.array([[1, 1, 1]], dtype=np.float64)   # logical that commutes with the checks (A·Hᵀ = 0)
dec = qec.get_decoder("tensor_network_decoder", H, logical_obs=A,
                      noise_model=[0.05, 0.05, 0.05],
                      contract_noise_model=True, device="cpu")
print(dec._dtype)   # -> 'float32', despite the docstring stating 'float64'

(This snippet reproduces only the documented-default contradiction — it always prints float32. The
non-physical posteriors below do not occur on typical syndromes; they require the deliberately
selected degenerate / high-weight tail at low p — see How to trigger below.)

Observed impact (rotated surface code, SI1000 circuit-level noise, seed 0)

All numbers below are from the canonical rounds = distance regime and are reproducible.

Scope — this is a rare-tail effect, deliberately surfaced, not a per-decode error. On typical
syndromes float32 and float64 agree. The failures appear only when the contraction is hardest, so
we intentionally select those syndromes: degenerate near-ties and/or high error weight, at low p
(large dynamic range). We sample a large pool and keep the highest-weight (most degenerate) tail; that
selection is what makes the bug visible. A random/typical decode will not show it.

d = 3 (bulk): no decision effect. Sweeping p ∈ [1e-3, 2e-2] (100k shots/p, 8 values):
float32 and float64 give identical logical error rates at every p (ΔLER = 0, 0 disagreements,
exact McNemar p = 1.0). Only the soft output differs, and only at low p (|ΔP|max ≈ 3e-2). At
d = 3 the float32 default corrupts the posterior but not the decision.

d = 5 (rounds = 5): float32 returns non-physical posteriors, with a statistically significant
decision effect at p = 1e-3. On the 1000 highest-weight (degeneracy-prone) syndromes per config,
using the raw coset masses Z0, Z1 stored in both precisions:

config	float32 negative masses	float64 negative masses	float32 P range
p=1e-3, r=3	142 / 1000	0 / 1000	[-3.08, 5.35]
p=1e-3, r=5	134 / 1000	0 / 1000	[-0.78, 4.38]
p=5e-3, r=3	7 / 1000	0 / 1000	[-0.82, 6.89]
p=5e-3, r=5	0 / 1000	0 / 1000	[0.004, 0.99]

• float64 never returns a negative mass or an out-of-[0,1] posterior on any of the 5000 shots.
• A coset mass is a sum over the (signed, Hadamard-basis) contraction; a negative result
  indicates loss of significance among signed contributions in float32 — it excludes
  underflow-to-zero (which would give 0, not a negative value), though the precise numerical
  pathway is not isolated. float64 returns a valid (positive) mass on the same inputs.
• Decision impact — on this same high-weight tail (p=1e-3, r=3, 1000 shots): discordant
  decisions are essentially one-directional — float32 wrong / float64 right on 10 shots, the
  reverse on 1 (exact McNemar p = 0.012). On the affected (negative-mass) shots, float32
  LER 0.049 vs float64 0.007 — also worse than the MWPM baseline 0.028 (float32 wrong
  while MWPM right on 5/142). The effect vanishes at p=5e-3, r=5 (0 invalid posteriors).

(For completeness, the same non-physical output is observed through the public decode() API at
d = 7; that run is preliminary — non-canonical rounds, small sample — and not relied upon here.)

Why this matters in practice

• It defeats the purpose of an exact decoder. The TN-MLD exists to provide accuracy; in float32
  it returns invalid posteriors (negative / outside [0,1]) on the hard tail — degrading the very
  precision it is meant to deliver.
• It can underperform the cheap baseline. On the affected shots the float32 MLD has a higher
  logical error rate than MWPM (0.049 vs 0.028 above): one pays for an exact maximum-likelihood
  decoder and, on those shots, gets a result worse than a fast minimum-weight matcher.
• It corrupts the TN decoder's use as a soft-output oracle. The TN-MLD is commonly used as a
  high-accuracy reference / oracle to train or benchmark ML / neural-network decoders, which consume
  its soft posteriors as ground-truth labels. In float32 it feeds non-physical labels (negative
  probabilities, values > 1, NaN) into that pipeline — poisoning training targets and skewing any
  benchmark that treats the TN output as exact. The silent float32 default undermines this use case
  precisely where the oracle is supposed to be most trustworthy (the rare, hard, degenerate shots).

How to trigger the non-physical posteriors (reproduction)

The minimal snippet above shows only the default-precision contradiction. To observe the negative
masses / out-of-[0,1] posteriors on the real decoder you must hit the hard tail:

1. sample many shots at low p (e.g. 1e-3);
2. keep the highest-weight / most-degenerate syndromes (the tail) — random/typical syndromes will
   not trigger it;
3. decode them in float32 vs float64 and compare the posteriors / coset masses.
   A standalone CPU re-derivation script (reads the saved per-shot masses; no GPU or cudaq needed) and the
   full per-shot data are available on request.

Expected behavior

For a maximum-likelihood decoder: the documented and actual default precision should agree, and
decode() should never return a probability outside [0, 1].

Suggested resolutions (increasing scope)

1. Fix the docstring (and the public API reference) to say "float32", or
2. change the default to "float64" to match the documented contract, or
3. add a validity guard in decode() (clamp / raise / surface the existing underflow flag) so a
   non-physical coset mass or out-of-[0,1] posterior is never returned silently.

Happy to share the full per-shot data and figures, or to open a PR for whichever direction you prefer.