Stabilize float32 AMICATorchNG for the GPU fast path

[neuromechanist]

From the #63 GPU benchmark on an RTX 4090: float32 is 5-8x faster than float64 (CUDA float32 35 ms/it vs CPU float64 281 ms/it) but collapses to NaN on full-size data.

Observed

On the full sample EEG (30504 samples), float32 descends healthily for ~22 iterations (-3.51 -> -3.44) then a mixture component's responsibility mass underflows in float32's ~7-digit range; the exact-EM 0/0 yields non-finite mu/beta/alpha and the degenerate-fit contract (#50) stops with nan_ll (fails loudly, good). float32 DOES converge on smaller slices (4096 samples: -3.234, matching float64). So the math is right; it's precision headroom.

Goal

Make float32 stable on full-size recordings so the GPU fast path (5-8x) is usable in production, keeping float64 as the parity default.

Directions

• float32-specific flooring of mixture responsibilities / dalpha_n before the exact-EM divisions (the collapse point).
• Mixed precision: preprocessing + the covariance/sphere already run float64; consider keeping the LL accumulation / a few sensitive reductions in float64 while the bulk runs float32.
• Investigate whether a larger invsigmin/minrho or a responsibility floor tied to torch.finfo(float32).tiny prevents the collapse without hurting float64.

Acceptance

• float32 completes a full-data fit without NaN and reaches an LL within a documented tolerance of float64.
• float64 parity (Verify NG init-basin vs Newton bit-parity via Fortran load_* init-matching #24) unchanged. Re-run benchmarks/benchmark_gpu.py on a CUDA host and record the stable speedup.

[neuromechanist]

Investigation: float32 is not cheaply stabilizable (PR #71 documents this)

Reproduced: float32 diverges to NaN on the full 30504-sample data across every seed, Newton on AND off (crash iter 9-81). Converges only on small slices (4096 samples -> -3.234, matching float64).

Root cause: NOT a divide-by-zero. Instrumented the M-step denominators (dmu_d/dbeta_d/dalpha_n) -- all healthy (O(1e3)) right up to the crash; then dmu_d goes NaN because a parameter already diverged the prior iteration and the next E-step overflows. It's precision-driven instability of the natural-gradient EM at ~7 sig digits.

Cheap targeted fix tried -> rejected: responsibility logsumexp/softmax in float64 (density + matmuls left float32) did NOT stabilize it (still NaN, all seeds). So the instability is in the density and/or the 30504-sample sufficient-stat accumulation.

Path forward (deferred): a real fix needs mixed precision (float64 density/LL/accumulation + float32 matmuls). But those sensitive parts dominate runtime (~50%+), so the realistic net speedup is only ~1.5-2x, not the raw 10-19x -- low payoff for a moderate refactor. Meanwhile float64-CUDA (4.5x, bit-safe) is the shipped production GPU path.

Keeping this open as the mixed-precision refactor, but deprioritized given the payoff analysis. Documented as experimental in the dtype docstring + perf_findings.md.

[neuromechanist]

Reopening as deferred: the investigation (documented in #71 / .context/issue-63/perf_findings.md) concluded float32 needs a mixed-precision refactor whose realistic payoff is only ~1.5-2x. Keeping this open as a deprioritized future enhancement, not active work. float64-CUDA (4.5x, bit-safe) is the production GPU path.

[neuromechanist]

Parking per decision: float32 mixed-precision is deprioritized (realistic ~1.5-2x payoff, documented in perf_findings.md). float64-CUDA (4.5x) is the production GPU path. MPS acceleration is being researched separately as its own pathway.