Symmetry-inferred `einsum` rejects an algebraically symmetric, all-finite float32 result with `max deviation = inf`; the check is `np.allclose(atol=1e-6)` on near-zero off-diagonals and the `inf` is a hard-coded sentinel

[jigonyoo]

Environment

	macOS arm64	Linux x86_64
Python	3.11.2	3.11.15
numpy	2.4.6	2.4.4
flopscope	0.12.1	0.12.1

The same input can pass on one platform and fail on the other (float32 accumulation order differs). Two inputs are supplied so the failure reproduces on both.

Minimal reproduction

The contraction R[a,b] = Σ_ij C[i,j] W[i,a] W[j,b] is exactly symmetric in (a,b) whenever C is symmetric. When C is tagged with flops.as_symmetric and the same tracked W object is passed in both slots, flopscope infers output symmetry and validates it.

import numpy as np, flopscope as flops, flopscope.numpy as fnp

# case A — 2x2, bisected on macOS arm64 (fails there; passes on Linux x86_64)
C_A = np.array([[622.4113159179688, 0.0050632040947675705],
                [0.0050632040947675705, 1357.942626953125]], dtype=np.float32)
W_A = np.array([[1183.504638671875, -151.5677947998047],
                [-16.548728942871094, -5003.23779296875]], dtype=np.float32)

# case B — 4x4, found on Linux x86_64 (fails there)
C_B = np.array([[ 2.561786 ,  4.6759605, -1.7735146,  1.6202407],
                [ 4.6759605, 29.628193 ,  6.8060284, 22.88294  ],
                [-1.7735146,  6.8060284, 10.365577 , 13.281559 ],
                [ 1.6202407, 22.88294  , 13.281559 , 42.211308 ]], dtype=np.float32)
W_B = np.array([[-136.78256 ,  -15.034283,   55.022297, -148.38255 ],
                [ 311.29367 , -206.92221 ,  -92.97633 ,   34.676193],
                [-241.32414 ,   96.26745 ,  107.79915 ,  -63.68914 ],
                [ 344.3328  ,    5.667633,   42.61409 , -187.45699 ]], dtype=np.float32)

def run(C, W, *, tagged):
    with flops.BudgetContext(flop_budget=10_000_000, quiet=True):
        cov = fnp.asarray(C)
        if tagged:
            cov = flops.as_symmetric(cov, symmetry=(0, 1))
        w = fnp.asarray(W)                          # same tracked object, both slots
        return np.asarray(fnp.einsum("ij,ia,jb->ab", cov, w, w))

Full reproducer with reporting (flopscope_symmetry_repro_v2.py) is at the bottom of this issue.

Actual result

Linux x86_64 / numpy 2.4.4 (pasted from the reproducer):

[A_2x2_macos_arm64]
  numpy einsum         finite=True  max|r|=3.40069e+10  max_asym=2.0625  allclose(r,r.T)=True
  flopscope untagged   finite=True  max|r|=3.40069e+10  max_asym=2.0625  allclose(r,r.T)=True
  flopscope tagged     finite=True  max|r|=3.40069e+10  max_asym=2.0625  allclose(r,r.T)=True

[B_4x4_linux_x86_64]
  numpy einsum         finite=True  max|r|=9.53528e+06  max_asym=0.25  allclose(r,r.T)=False
  flopscope untagged   finite=True  max|r|=9.53528e+06  max_asym=0.25  allclose(r,r.T)=False
  flopscope tagged     SymmetryError: Tensor not symmetric along axes (0, 1): max deviation = inf (tolerance: atol=1e-06, rtol=1e-05).

macOS arm64 / numpy 2.4.6 (case A, pasted from the original 2×2 reproducer):

numpy optimized finite= True max_abs= 34006855680.0 max_asymmetry= 8.0 allclose_transpose= False
flopscope untagged finite= True max_abs= 34006855680.0 max_asymmetry= 8.0 allclose_transpose= False
flopscope tagged SymmetryError Tensor not symmetric along axes (0, 1): max deviation = inf (tolerance: atol=1e-06, rtol=1e-05).
flopscope tagged, covariance x0.01 finite= True max_abs= 340068480.0 max_asymmetry= 0.0 allclose_transpose= True

Why this is not an overflow

Every input and every output is finite on both platforms. Scaling C by 0.01 makes the tagged call pass. The inf in the message is not measured.

Root cause (flopscope 0.12.1 source, line numbers verified against the installed package)

The validation is np.allclose with atol=1e-6, rtol=1e-5:

• flopscope/_symmetric.py:146 — def _check_generators(array, group, *, atol=1e-6, rtol=1e-5)
• flopscope/_symmetric.py:160 — if not np.allclose(array, array.transpose(perm), atol=atol, rtol=rtol):
• (also _symmetric.py:65, :334)

np.allclose passes element (a,b) iff |r_ab − r_ba| ≤ atol + rtol·|r_ba|. For an off-diagonal entry that is small relative to the matrix scale, rtol·|r_ba| is tiny and the test degenerates to atol = 1e-6, which float32 accumulation noise at that scale exceeds. In case B the offending element is:

result matrix scale: max|r| = 9.535e+06
fails at (a,b)=(0,2): |r_ab − r_ba| = 0.1167   allowed = 1e-6 + 1e-5·|r_ba| = 0.08097   (|r_ba| = 8097)

i.e. the entry is three orders of magnitude below the matrix scale, and the tolerance shrinks with it while the rounding noise does not.

The inf is a hard-coded sentinel raised on the failure branch rather than the measured deviation:

• flopscope/_symmetric.py:58 — raise SymmetryError(axes=group, max_deviation=float("inf"))
• flopscope/_symmetric.py:325 — same
• flopscope/_pointwise.py:511, :536 — same

So the message tells the user "infinite deviation" when the actual deviation was 0.1167 (case B) or 8.0 (case A).

Expected behaviour / suggestion

For a contraction whose output symmetry follows algebraically from a tagged symmetric input and a repeated operand, one of:

1. Canonicalise the result (0.5·(R + Rᵀ)) instead of validating it — the symmetry is known a priori, not inferred from data; or
2. Scale the tolerance to the array, e.g. atol = eps_f32 · max|R| · k or rtol applied against max|R| rather than the per-element |r_ba|; or
3. If rejection is intentional, report the measured deviation (and the offending index) instead of the inf sentinel, and document that users must project explicitly.

Impact

• Any estimator that (a) tags a covariance as symmetric and (b) contracts it against a repeated weight operand — the natural way to write Wᵀ C W — can hit a hard exception on ordinary float32 inputs, with the failure depending on the grading machine's accumulation order. Under the Phase 2 grader that is the zero-prediction fallback for that MLP.
• In our own measurements the failure boundary sat far from He-initialised weights: on 45 genuine N(0, 2/n) MLPs at width 1024 the margin to the boundary was 53.8 standard deviations, so we saw no submission risk at the competition's shape. The issue is a correctness-of-reporting and tolerance-design bug, not a "my submission failed" report.
• The misleading inf cost us a day of investigation in the wrong direction (overflow) before the source was read.

Not a second bug

fnp.linalg.svd_cost does not exist in 0.12.1; fnp.linalg.svdvals_cost does. That was our error, not flopscope's.

Related: #256 (explicit-output einsum broadcast failure) — a different code path.

flopscope_symmetry_repro_v2.py — full reproducer

"""Minimal finite-value reproducer for flopscope 0.12.1 SymmetryError (two platforms).

The contraction  R[a,b] = sum_ij C[i,j] W[i,a] W[j,b]  is exactly symmetric in
(a, b) whenever C is symmetric.  When C is tagged with flops.as_symmetric and the
SAME tracked W object is passed twice, flopscope infers output symmetry and then
validates it with np.allclose(atol=1e-6, rtol=1e-5).  A few ulps of float32
accumulation-order noise on near-zero off-diagonal entries is enough to fail
that check, and the raised error reports "max deviation = inf" -- a hard-coded
sentinel, not a measured value.  No overflow is involved: every array is finite.

Because the noise depends on the platform's BLAS/accumulation order, the SAME
input may pass on one machine and fail on another.  Two inputs are included:
  case A  2x2, bisected on macOS arm64 / numpy 2.4.6  (fails there, passes on x86_64)
  case B  4x4, found on Linux x86_64 / numpy 2.4.4     (fails there)
"""
from __future__ import annotations
from importlib.metadata import version
import platform
import numpy as np
import flopscope as flops
import flopscope.numpy as fnp

CASES = {
    "A_2x2_macos_arm64": (
        np.array([[622.4113159179688, 0.0050632040947675705],
                  [0.0050632040947675705, 1357.942626953125]], dtype=np.float32),
        np.array([[1183.504638671875, -151.5677947998047],
                  [-16.548728942871094, -5003.23779296875]], dtype=np.float32),
    ),
    "B_4x4_linux_x86_64": (
        np.array([[ 2.561786 ,  4.6759605, -1.7735146,  1.6202407],
                  [ 4.6759605, 29.628193 ,  6.8060284, 22.88294  ],
                  [-1.7735146,  6.8060284, 10.365577 , 13.281559 ],
                  [ 1.6202407, 22.88294  , 13.281559 , 42.211308 ]], dtype=np.float32),
        np.array([[-136.78256 ,  -15.034283,   55.022297, -148.38255 ],
                  [ 311.29367 , -206.92221 ,  -92.97633 ,   34.676193],
                  [-241.32414 ,   96.26745 ,  107.79915 ,  -63.68914 ],
                  [ 344.3328  ,    5.667633,   42.61409 , -187.45699 ]], dtype=np.float32),
    ),
}

def run(C, W, *, tagged):
    with flops.BudgetContext(flop_budget=10_000_000, quiet=True):
        cov = fnp.asarray(C)
        if tagged:
            cov = flops.as_symmetric(cov, symmetry=(0, 1))
        w = fnp.asarray(W)                       # same tracked object in both slots
        return np.asarray(fnp.einsum("ij,ia,jb->ab", cov, w, w))

def describe(label, r):
    asym = float(np.max(np.abs(r - r.T)))
    print(f"  {label:<20} finite={bool(np.isfinite(r).all())}  max|r|={float(np.abs(r).max()):.6g}"
          f"  max_asym={asym:.6g}  allclose(r,r.T)={bool(np.allclose(r, r.T, atol=1e-6, rtol=1e-5))}")

def main():
    print("python", platform.python_version(), "|", platform.platform())
    print("numpy", np.__version__, "| flopscope", version("flopscope"))
    for name, (C, W) in CASES.items():
        print(f"\n[{name}]")
        describe("numpy einsum", np.einsum("ij,ia,jb->ab", C, W, W, optimize=True))
        describe("flopscope untagged", run(C, W, tagged=False))
        try:
            r = run(C, W, tagged=True)
            describe("flopscope tagged", r)
        except flops.SymmetryError as e:
            print(f"  {'flopscope tagged':<20} {type(e).__name__}: {e}")

if __name__ == "__main__":
    main()