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

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:

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.
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?

@jimsweeney

Thanks for this detailed write-up @jimsweeney, and congrats on submitting the first cuPhoton issue! GPU perf baselines in a separate PR would be preferred, as that'll allow us to review the numerical behavior independently.

For the validation work, I’d suggest:

Reusable fixture and sweep code under src/cuphoton/xray/, with tests under tests/xray/.
A new docs/xray/LINEAR-PREDICTION-VALIDATION.md describing the model, assumptions, schema, and reproduction commands, linked from the XRay README.
A caller-selected output directory containing summary.json and the figure. Keep generated runs outside the checkout; attach a representative result to the PR.

We don’t currently have a dedicated schema for this validation. A small schema_version: 1 summary would be a good starting point, with:

config: ground-truth parameters and units, sampling, noise/SNR definition, seeds, trial counts, model order, and mode-matching criteria.
runtime: source revision, package versions, backend, dtype, and relevant hardware details. Please reuse cuphoton.core.runtime.runtime_metadata where applicable/if applicable.
results: one entry per sweep condition, with attempted/successful/failed trials, mode-loss rates, parameter bias, standard deviation, RMSE, and the corresponding Cramér–Rao bounds. Please distinguish variance bounds from standard-deviation bounds and label statistics computed only over successful recoveries.
artifacts: relative filenames for the derived figure and any additional numeric products.

It sounds like you've already prepped a command (lpv), which is presumably linear-prediction-validate, which is fine and fits the existing naming.

For the benchmark PR, a table under docs/xray/ sounds good. Please include reproduction commands, hardware/software details, warmup, repetitions, GPU synchronization, timing statistic, and whether transfers are included.

We can review the root-filter tolerance and conjugate-pair handling with the reproducer and proposed code before settling on a change to the default behavior.

As a side-bar, we have related fit-diagnostics, optional iterative-fitting, and GPU-batching work in progress. We’ve also encountered decay-sign sensitivity near the unit circle when comparing numerical paths. Your noise-driven reproducer and Cramér–Rao/mode-loss study would complement that work.

If you have PRs ready or close to ready, feel free to post when able. We'll figure out the overlap/details from there.

Thanks, that layout is what I have followed. Three pull requests:

Add synthetic fixture and Cramer-Rao validation for xray linear prediction (refs #1) #2: the fixture, the Cramer-Rao sweep and linear-prediction-validate with the schema_version 1 summary and --output-dir, documented in docs/xray/LINEAR-PREDICTION-VALIDATION.md;
Document how to run and interpret XRay benchmarks #3: the benchmark table under docs/xray/ with the measurement procedure, plus a profile of the P1 GPU path: cusolver.xgeev runs once per companion matrix even on the batched call, which is why batching gains only a few percent at 96 samples;
Add experimental per-trace and batched mode refinement #4: the per-trace and batched refinement, stacked on Add synthetic fixture and Cramer-Rao validation for xray linear prediction (refs #1) #2, posted because it overlaps the iterative-fitting and GPU-batching work you mentioned; take what is useful.

On the decay-sign sensitivity near the unit circle: that is the same mechanism the sweep exposes (a lightly damped mode whose fitted decay crosses zero is dropped), so the per-level loss rate in the summary should give a number for it on any fixture you care about.

@jimsweeney

Thanks @jimsweeney! I'll review (hopefully today, more likely tomorrow) and get back to you.

@jimsweeney

Hey @jimsweeney, I've reviewed the PRs, they're good, a few tweaks will be necessary; can you set me up to push to your fork?

There's a slight hiccup with timing, we've got another patch release coming up soon (around the end of the month) and we've made a bunch of cuphoton.xray changes that'll impact these PRs (and there's a small amount of overlap in what you've done and some of the stuff we've done). After this patch release we're intending to switch more of our development process over here to Github publicly so in-process PRs are visible and can bet tracked etc.

I'm going to investigate how we can include your work in this next release whilst still preserving you as the contributor. I'll review over the next few days and figure out a plan. If you could enable push access for me to your repo that'll at least let me convey initial agent feedback.

Thanks for the review, Trent. I've sent you a collaborator invite on jimsweeney/cuPhoton with push access, so you can put feedback straight on the branches.

The release plan sounds right to me. I'm fine with you folding the work into the patch release in whatever shape fits the xray changes you already have in flight, and I appreciate you keeping the contributor attribution. If it's easier to rebase my branches onto your updated xray code once it lands, tell me and I'll do that on my side.

Add synthetic damped-mode validation for XRay