[NumPy Parity] Neural Networks Project to fully train byte-identical weights on training

Overview

Embed a call-exact NumPy twin (mnist_mlp.py) next to the C# MNIST MLP demo (examples/NeuralNetwork.NumSharp/MnistMlp/) and make both sides produce byte-identical weights after a full 100-epoch training run — verified by .npy byte-equality (NumSharp's writer is already byte-identical to np.save, so the gate is literal file-hash equality). This extends the differential-fuzz philosophy (NumPy is the oracle, bit-exact or excused-loudly) from single ops to a complete composed workload: init → forward → loss → backward → Adam × 4,600 steps.

Problem

The demo trains and converges, but nothing proves NumSharp's training math is NumPy's. Per-op fuzz tiers pass individually, yet two documented excuse classes sit directly on the weight path, and the GEMM diverges at real sizes despite being fuzz-green at the corpus' K≤4 shapes.

Empirical probe results (this session; numpy==2.4.2 win-amd64 wheel, AVX2+FMA host; NumSharp replayed NumPy-saved .npy inputs and bit-compared):

Every other weight-affecting op is IEEE-exact by construction and already fuzz-green: add/sub/mul/div (incl. bias broadcast and the NDExpr fused Max(in0+in1,0) / in0*Greater(in1,0) kernels — plain exact ops), maximum, sqrt (correctly-rounded), astype f64→f32, max(axis=1) (order-independent), one-hot integer writes, view slicing. The softmax → CE backward is (softmax − labels) * (1/128) — exact (1/128 is a power of two). np.log affects only the printed loss, not weights — log parity is an optional stretch goal.

Why the GEMM diverges — and why no stride trick can fix it (NumPy 2.4.2 source finding)

numpy/_core/src/umath/matmul.c.src (@TYPE@_matmul, the "copy if not blasable, see gh-12365 & gh-23588" branch): for float32/float64 matrix@matrix NumPy copies any non-blasable operand into a contiguous temp and still calls cblas — the pure-C matmul_inner_noblas loop (acc += a*b, ascending k) is unreachable in the stock wheel for f32 mat@mat (only zero-dim or >BLAS_MAXSIZE reach it). Probes confirm:

• np.dot ≡ np.matmul bit-equal (0/16,384 diffs).
• BLAS result ≠ sequential mul+add chain ≠ sequential-FMA chain (per-element references, 8 probes each) → scipy-openblas sgemm uses an arch-specific multi-accumulator register scheme.
• Even BLAS-vs-BLAS differs by transpose variant: forcing a copy via a negative-stride column view changed 5,487/16,384 results (different transA/B kernel).

So byte-matching the stock wheel means matching that OpenBLAS binary's accumulation, not any portable algorithm.

Proposal

GEMM strategy (decide via spike, one seam: ParityMatmul)

NumSharp's existing GEBP FMA GEMM is untouched — the parity path is opt-in for the harness/gates.

Task checklist

• RNG oracle tier — committed corpus pinning seed → standard_normal/normal doubles → astype(f32) streams in CI (P1 is proven live but only state-round-trip tested today; Randomizer.Tests.cs never compares against real NumPy output)
• exp f32 port — port simd_exp_f32 (FMA polynomial) from src/numpy/numpy/_core/src/umath/loops_exponent_log.dispatch.c.src (clone already in-repo); elementwise → lane-independent → identical bits at any vector width on FMA hardware. Narrow the blanket unary ~ULP ≤2 MisalignedRegistry branch to exclude exp/Single (house pattern: scoped branch + tightness test) and add corpus cases asserting exp/f32 bit-exact. Stretch: same for log f32 → byte-identical printed losses
• GEMM spike + implement — S-A vs S-B behind one ParityMatmul seam; new matmul_parity corpus tier with K ∈ {255, 256, 257, 512, 784} (the KC=256 block boundary the current K≤4 tier never crosses) × transposed-view operands
• Sum floors — corpus cases pinning sum(axis=1) contiguous and sum(axis=0) strided f32 bit-exact at reduction lengths ≤128 (already passing; pin so it can't regress)
• Demo determinism changes (both twins identically) — synthetic data via np.random MT19937 (replacing System.Random + local Box-Muller, which has no NumPy analog); np.power(grad, 2) → grad*grad; Adam bias correction Math.Pow(β, t) → cached running products β_pow *= β (removes CRT-pow cross-platform risk); real-MNIST normalize mirrored as x * (1f/255f)
• Scalar-expression mirroring audit — document every scalar rounding point. Example of the discipline: C# 1 - Beta1 computed in f32 is 0x3DCCCCCE, while the literal 0.1f is 0x3DCCCCCD — the twin must write np.float32(1) - np.float32(0.9), never the simplified constant
• mnist_mlp.py twin — at examples/NeuralNetwork.NumSharp/MnistMlp/parity/, file-by-file mirror (Program / MlpTrainer / FullyConnectedFused / SoftmaxCrossEntropy / Adam / MnistLoader-synthetic) with C# file/line cross-references in headers
• Staged byte-gate harness — --parity-dump <dir> flag on both twins dumping .npy at gate points; comparison = file byte-equality; test/oracle/verify_nn_parity.py manual gate (needs Python + SDK, like verify_npy_interop.py); optional nightly
• Run the ladder, document — results + host pin into the example's CLAUDE.md; any surviving divergence into the MisalignedRegistry-style ledger (excused loudly, never silent)

Stage ladder (each gate isolates one pillar; fail-fast localization)

Evidence

• Probe harness: NumPy 2.4.2 generated .npy inputs; NumSharp loaded them (byte-exact reader) and recomputed; per-element uint32 bit-compare. RNG/sums bit-exact; exp ≤2 ULP on 39% of elements; GEMM as tabled above.
• matmul.c.src routing read from the in-repo NumPy 2.4.2 clone (src/numpy/): noblas loop *(typ*)op += val1 * val2 ascending-k (lines ~282–320); the copy-then-BLAS branch (gh-23588) makes it unreachable for stock-wheel f32 mat@mat.
• Fuzz ledger (test/NumSharp.UnitTest/Fuzz/README.md): unary ~ULP excuse (563 hits) and reduction summation/two-pass precision excuse (401 hits) are the two permanent branches this work narrows; the host-pinned-parity precedent already exists there (cast kernels reproduce the MSVC NumPy answer — "a NumSharp regression pin rather than portable NumPy parity").
• np.random.cs NextGaussian() is NumPy's exact polar method with gauss-cache (matches legacy_gauss incl. the x·d cached / y·d returned convention) — the structural reason P1 probes pass.
• NumSharp .npy writer is byte-for-byte np.save output (the NpyOracle gate, 286 cases) — which is what makes "weights parity" checkable as plain file equality.

Scope / Non-goals

• Scope: the float32 MnistMlp path only (FullyConnectedFused, SoftmaxCrossEntropy, Adam, MlpTrainer). Deterministic batch order (already no shuffle).
• Host pin: parity target is numpy==2.4.2 wheels on FMA-capable x86-64, single-threaded (OPENBLAS_NUM_THREADS=1 on the twin; NumSharp multithreading stays default-off). Same pin class as the MSVC cast-kernel precedent. Cross-OS byte-parity only under S-B.
• Non-goals: the vanilla NeuralNet/FullyConnected/Sigmoid/BCE scaffolding; Half/Decimal networks; loss-print parity (needs the log port — stretch); shuffling/validation-split features; replacing NumSharp's native GEMM (parity path is opt-in).

Breaking changes

Performance

Parity paths are opt-in; default fast paths untouched. The exp port replaces the MathF.Exp scalar call inside the IL kernel — NumPy's polynomial is Vector256/FMA-vectorizable, so expect ≥ current throughput; gate with the benchmark/ harness before/after. S-A GEMM parity runs at OpenBLAS speed on both sides; S-B makes only the python twin slow.

Related issues

• NumPy: numpy/numpy gh-23588, gh-12365 (the copy-then-BLAS matmul routing that closes the noblas door)
• In-repo: fuzz ledger excuse branches this narrows (unary ~ULP, reduction summation); example CLAUDE.md perf-drift (fusion probe now 0.73×, training 12.8 s — stale doc, separate follow-up)