xray linear prediction: lightly damped modes are dropped by the `decay >= 0` root filter under noise; proposal for a synthetic fixture with Cramer-Rao validation and per-GPU baselines

CONTRIBUTING、既存の synthetic_trace fixture、_linear_prediction_modes_from_roots から始め、現在の動作を確認するために pytest tests/xray を実行します。コードを変更する前に、提案された fixture、検証出力、root pairing、docs/xray の baseline の場所を maintainers と確認します。完了の条件は、合意されたスコープ、mode loss に対する regression coverage、承認済みの artifact および benchmark schema です。

Component: cuphoton.xray (linear prediction). cuPhoton main at 9e91835 (2026-09-14). Observations 1 and 2 are on the NumPy path (Python 3.12, Windows); the GPU environment (WSL2 Ubuntu 24.04, CUDA 13.2, CuPy 14.1.1, locked --extra gpu) passes all 44 xray linear-prediction tests and was used for the baseline in the proposal.

While building a synthetic validation for the linear-prediction fit I found a behavior that looks worth a look, and I would like to propose a fixture and test around it before writing code, per CONTRIBUTING.

What I ran: the repository's own synthetic_trace(96) (two damped modes, 2.4 rad/s with decay 0.09 and 0.9 rad/s with decay 0.03, plus a constant), white Gaussian noise added at 40, 30, 20, 15, 10 and 5 dB signal-to-noise (signal rms 0.660), 300 draws per level, linear_prediction_numpy(..., n_components=6), which recovers the clean trace exactly.

Observed:

1. At 40 dB, 42 of 300 draws lose the 0.9 rad/s mode entirely. In every one of those draws selected_model_order is 6 and decaying_root_count is 3: the strong mode and a zero-frequency root survive, the weak, lightly damped mode does not. n_components of 5, 7 or 8 gives 42, 45 and 45 failures respectively, so it is not a model-order choice. Instrumenting the root set on the failing draws pins it to the decay >= 0 filter in _linear_prediction_modes_from_roots: in every failing draw the 0.9 rad/s pair is present in the active root set at the right frequency (0.908 to 0.910 rad/s) but with a fitted decay of about -0.004 to -0.009 (|z| = 0.9991 to 0.9996), so noise has moved a root whose true decay is +0.03 just inside the unit circle and the filter discards the mode as growing. With this mode's decay scatter of roughly 0.03 at 40 dB, a sign flip about 14 percent of the time is what one would expect. Long-lived phonon modes are exactly this case. A side observation from the same instrumentation: the truncation sorted_roots[:singular_value_count] after a lexicographic complex sort can cut between the two members of a conjugate pair (the 24th root in the failing draws is an unpaired high-frequency spurious root); harmless here, but the conjugate-pair split in select_count assumes it never happens.

2. Where the fit succeeds, the estimator's scatter sits above the Cramer-Rao lower bound computed from the exact Fisher information of the nine-parameter model: about 3.8, 3.0, 2.9, 2.7, 2.2 and 1.7 times the bound for the strong mode's frequency from 40 dB down to 5 dB, and up to 6 times for its decay rate at 40 dB. Frequency bias stays below 0.003 rad/s down to 20 dB and grows to 0.02 rad/s at 5 dB. None of this is surprising for a linear-prediction estimator, but a table like this in the docs would let a user decide when the numbers can be trusted and when a refinement step is needed.

Proposal, matching the ROADMAP items on synthetic fixtures and performance baselines:

• A synthetic fixture generator alongside synthetic_trace: known modes, seedable noise, optional detector artifacts (structured non-stationary flat-field drift), returning ground truth with the trace, usable without external data.
• A validation that sweeps signal-to-noise, trace length, mode separation and model order, reports estimator standard deviation against the Cramer-Rao bound, and records the mode-loss rate; one numeric artifact (summary.json) and one figure; a test asserting the estimator stays within a stated factor of the bound in the easy regime.
• A small change or option for the filter: keep roots with decay above a negative tolerance tied to the trace length (a mode that decays by less than a few percent over the record cannot be distinguished from a slowly growing one by its root magnitude alone), or select by fitted amplitude with the growing roots reported rather than silently dropped; plus a pairing step that tolerates an unpaired spurious root; with a regression test built on the fixture.
• Separately, the P1/P2/Savgol benchmarks on two GPU architectures with the architecture, driver, CUDA, dtype and backend recorded, as a markdown table under docs/xray. First data point, RTX 4050 Laptop (Ada, cc 8.9) under WSL2, CUDA 13.2, CuPy 14.1.1, locked GPU environment, best of 3: P1 at 96 samples is about 15 times slower on the GPU than on the CPU at 16, 256 and 2048 traces (CPU 0.9 ms per trace, GPU 14 ms per trace, batched only 8 to 10 percent faster than serial GPU); P2 and Savitzky-Golay batch about 16 times over serial GPU but still trail the CPU at 16 traces. I would profile before drawing a conclusion, but it suggests the batched P1 path serializes the eigenvalue solve. A second card follows.

I have the fixture, the sweep, an lpv command, tests and a docs section ready on a branch (ruff clean; pytest tests/xray passes on Linux with the CUDA 13 extras) and will open the pull request against this issue. Which artifact schema and file locations would you prefer for the validation output, and would you rather have the baselines in the same PR or a separate one?