Vectorized time-marginalized likelihood offsets by the BATCH max, not the per-sample max: lnL underflows to -inf above rho ~ 40 and mcsamplerAV collapses (device-independent; affects a delivered O4b result)

[oshaughnessy-junior]

Summary

DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop offsets its time integral by the
batch maximum instead of the per-extrinsic-sample maximum:

# RIFT/likelihood/factored_likelihood.py
lnLmax  = xpy.max(lnL_t)                      # lnL_t is (npts_extrinsic, npts_time) -- NO axis
L_t     = xpy.exp(lnL_t - lnLmax, out=lnL_t)
L       = simps(L_t, dx=deltaT, axis=-1)
lnL     = lnLmax + xpy.log(L, out=L)

Every row is shifted by the batch peak, so any extrinsic sample whose own peak lnL
sits more than ~745 nats below the loudest sample in the batch has exp() underflow to
0 across the whole time axis, L = 0, and log(L) = -inf. The likelihood returns
-inf for a sample at which it is perfectly finite.

Since a typical prior draw has lnL ~ 0 and the peak has lnL ~ rho^2/2, the defect
switches on once max lnL exceeds the float64 underflow budget (~745 nats), i.e. around
rho ~ 40, and then applies to the bulk of the prior.

This is not device-dependent. The path is selected by opts.gpu, and --force-xpy
keeps opts.gpu true even when cupy is absent, so --gpu --force-xpy on a machine with
no CUDA reproduces it exactly. cupy is not implicated. It is also not the CPU/GPU
Simpson divergence (#204), which is bounded to [ln(2/3), ln(4/3)] nats and cannot
produce -inf.

Present at tag 0.0.17.4 (line 2033), on rift_O4d (2055), and at head.

The one line, measured

Real S250114ax H1/L1 strain, SEOBNRv5PHM lmax=4, srate 4096, rho ~ 82,
--sampler-method AV --n-max 4e6 --n-eff 100, seeds 3001/3002/3003, one container
(rift_o4d_cc90-120_cuda128_20260717.sif), one RIFT tree (bd5c6fa), one host
(ldas-pcdev11), single-CPU. Backend proved in-process in every arm
(mcsamplerGPU.xpy_default is cupy asserted; no cupy / Override --gpu counted).

Finite likelihood values in the first 10 000 draws:

arm	backend	finite / 10000
no --gpu	numpy	9929
--gpu --force-xpy, cupy loaded, device 0	cupy	6
--gpu --force-xpy, cupy not loaded	numpy	11

Same tree, same seed. The flag, not the device.

Full budget (4e6), paired, only that one line differing between arms:

arm	chi-square ESS	Pareto k-hat	reported lnL
--gpu, as shipped	1.0 / 1.0 / 1.0	21.8 / 16.2 / 23.6	3365.57 / 3346.74 / 3381.48
--gpu, per-row offset	146.6 / 152.9 / 102.5	0.44 / 0.34 / 0.51	3370.95 / 3371.32 / 3370.83
CPU reference (no --gpu)	102.6 / 104.0 / 123.0	0.51 / 0.64 / 0.47	3371.48 / 3370.71 / 3371.11

The shipped --gpu arm scatters over 34.7 nats across three seeds while reporting
sigma = 1.03 — understating its own error by more than an order of magnitude (k-hat
16-24 says as much). With the per-row offset it agrees with the CPU reference to 0.07
nats, inside the 0.26-0.38 nat replicate spread.

Production blast radius

The O4b campaign test_o4b_eventlist_6 is 141 runs through one container
(Reviewed-GWTC5-RIFT-20250919.sif) with an identical flag set
(--gpu --force-xpy --force-gpu-only --sampler-method AV --n-eff 10 --n-max 4e6 --time-marginalization, distance sampled). Reading each run's own
sqrt(var)/res column out of all.net (that column is 1/sqrt(chi-square ESS) in
this code):

	median chi-square ESS	fraction of points that exhausted --n-max
S250114ax / rift-v5PHM-calmarg	1.73	98.2 %
S250114ax / rift-NRSur7dq4	1.83	100 %
the other 138 runs	77 (q05 59, q95 90)	0.31 % (q95 2.9 %)

Exactly two runs in the campaign have max lnL > 745, and they are these two.
Nothing in the campaign sits between max lnL 666 and 2905, so the campaign brackets
the onset rather than resolving it — but it contains no counterexample in either
direction.

--force-gpu-only is sys.exit(35) when cupy fails to load, so every row of the
delivered S250114ax all.net was produced with cupy on a real device: this is not the
silent-numpy-fallback failure mode.

For scale: the delivered S250114ax intrinsic lnL surface spans 23 nats across the
top 10 % of its 5327 points — the surface CIP fits — against the per-point scatter
measured above. (That scatter was measured at srate 4096 with --internal-use-lnL, not
at production's srate 1024, so treat it as an order of magnitude, not a transfer.)

Grep the class, do not fix one line

The correct idiom xpy.max(..., axis=-1, keepdims=True) is already used at five sites
(factored_likelihood.py:3127, factored_likelihood_freqresponse.py:457,
factored_likelihood_with_rotation.py:948, time_marginalization_quadrature.py:764,
study_stencil_lnL_sensitivity.py:396). Three sites lack it:

site	expression	status
factored_likelihood.py:2953	lnLmax = xpy.max(lnL_t)	confirmed live (this issue)
factored_likelihood.py:2181	lnLmax = np.max(lnL_t_accum), lnL_t_accum is (npts_extrinsic, npts)	same shape, same exposure — not exercised here
factored_likelihood.py:3132	m_c = xpy.max(lnL_t_c) feeding a global running_max for the calmarg sum	in-loop calmarg only — not exercised here

The two unexercised ones are the same construction, not the same diagnosis: each needs
its own reproduction before anything is claimed about it.

What I am not proposing

Per standing policy I am not submitting a patch. Changing this moves production
numbers on the --gpu path for every loud event, and a change of that class needs
end-to-end known-answer runs, not a diff. The one-line per-row offset above is a
diagnostic intervention run in a private snapshot to confirm the mechanism, and it
is reported as evidence, not as a proposal.

Two consequences worth deciding on separately:

1. Whether any delivered --gpu result above rho ~ 40 needs re-running.
2. Whether mcsamplerAV should refuse to report a result at k-hat > 10 / ESS < 2 rather
   than returning a finite lnL with a sigma that understates its error 17-fold. The
   [AV COLLAPSE] diagnostic already prints; it does not gate.

Evidence, harnesses and logs: ~/av_gpu_ownership on ldas-pcdev11/CIT shared home
(run_arm.sh carries the in-process backend assertion; Code/ and Code_fix/ are the
pinned snapshots). Records:
RIFT_roboto_paper:development/OPEN_gpu_adaptive_volume_collapse.md.

[oshaughnessy-junior]

Fix PRs, at RO'S's direction (this supersedes the "diagnose and report, no patch" line in the issue body):

target	PR	fix site
oshaughn/research-projects-RIT:rift_O4c	oshaughn#187 (cross-fork)	line 2033
oshaughnessy-junior/research-projects-RIT:rift_O4d	#234	line 2953

Both sources are branches on the oshaughnessy-junior fork. Both carry the same change — xpy.max(lnL_t, axis=-1, keepdims=True) plus lnLmax[...,0] on the add-back — a regression test, and CI wiring in the same commit.

axis=-1 alone is not the fix. Without keepdims=True the offset is (npts_extrinsic,) and broadcasts along the time axis instead — silently wrong on a non-square batch, and a ValueError only when the shapes happen not to conform. Both PRs' tests cover that mutation.

Verified independently of the agents that wrote the PRs, by reverting the fix inside a throwaway worktree of each PR head and re-running that PR's own test:

• factored_likelihood: offset the vectorized time marginalization per extrinsic sample, not per batch (#232) #234 / test_noloop_time_marg_row_offset.py — 4 passed 1 skipped patched, 2 failed reverted: array([1685.71347063, -inf]). Mutated to bare axis=-1: 3 failed.
• D1: AlternateIteration's extrinsic stage read the previous iteration's grid #187 / test_time_marginalization_perrow_offset.py — 5 passed patched, 1 failed reverted: lnL = [1491.28, 891.28, -inf, -inf, -inf].

Two notes for reviewers:

1. This is not a no-op and is not claimed to be. Per-row vs batch offsets change float64 rounding even where nothing underflows; measured agreement on well-resolved cases is ~1e-14 nats (≤ a few ulps). Above the underflow budget it changes the answer — that is the point.
2. rift_O4c had no .github/workflows and an empty # TODO: write some unit-test stage, so D1: AlternateIteration's extrinsic stage read the previous iteration's grid #187 also creates the job that runs its test. An unlisted test file never runs in this CI.

The two sibling sites named in the issue (factored_likelihood.py ~2181, and the in-loop calmarg running_max, which does not exist on rift_O4c) are untouched in both PRs.

[oshaughnessy-junior]

• Closed in O4d
• at time of writing, still open in O4c