Fix full Wiener density underflow

[martonaronvarga]

Summary

Fixes #3329.

Adds an automatic log-scaled fallback for the full Wiener density when
natural-scale integration produces a zero, subnormal, or non-finite result.
The public api is unchanged.

Tests

• Prim Wiener full tests: 23/23 passed
• Rev Wiener full tests: 5/5 passed
• Standalone header compilation passed
• cpplint passed
• Benchmarks found no measurable slowdown in ordinary support regions

Side Effects

None expected.

Release notes

Fixed underflow in the full Wiener log density and its gradients for extreme
valid parameter values.

Checklist

• Copyright holder: Marton A. Varga

• the basic tests are passing

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[martonaronvarga]

@avehtari This is the PR fixing underflows mentioned in #3329. Unfortunately, I don't have permissions to assign reviewers, so please assign yourself for the time being.

[avehtari]

This doesn't compile for forward or mixed mode.** CI job "mix fun and prob tests 5" fails on all three platforms with 88 errors from three lines:

• line 852: std::isfinite(density) with var or fvar;
• line 854: wiener7_log_params{y_value, …}, which narrows T_partials_return to double;
• line 863: xmin/xmax passed as const Eigen::VectorXd&.

[martonaronvarga]

Sorry, forgot about mixed modes. Fixed in bd47838.

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.04	0.04	1.02	1.85% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	0.99	-0.97% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.86	1.01	0.7% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.0	-0.22% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.01	0.85% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	7.97	8.07	0.99	-1.32% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.74	3.74	1.0	-0.15% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.27	33.12	0.97	-2.62% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	1.01	1.34% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.95	6.95	1.0	-0.0% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	135.29	129.82	1.04	4.05% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.21	3.22	1.0	-0.24% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.61	0.99	-1.37% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.75	25.88	0.99	-0.5% slower
performance.compilation	407.17	399.05	1.02	1.99% faster
Mean result: 1.0025071235394414

Jenkins Console Log
Jenkins Build Stages
Commit hash: 44a50912ff3f2e71023ea29a6f91adf7c6b1cf87

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           52 bits physical, 57 bits virtual
CPU(s):                                  192
On-line CPU(s) list:                     0-191
Thread(s) per core:                      2
Core(s) per socket:                      48
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              25
Model:                                   17
Model name:                              AMD EPYC 9474F 48-Core Processor
Stepping:                                1
Frequency boost:                         enabled
CPU MHz:                                 1497.776
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.39
Virtualization:                          AMD-V
L1d cache:                               3 MiB
L1i cache:                               3 MiB
L2 cache:                                96 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-47,96-143
NUMA node1 CPU(s):                       48-95,144-191
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Not affected
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Enhanced / Automatic IBRS; IBPB conditional; STIBP always-on; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Mitigation; Clear CPU buffers
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good amd_lbr_v2 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba perfmon_v2 ibrs ibpb stibp ibrs_enhanced vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a avx512f avx512dq rdseed adx smap avx512ifma clflushopt clwb avx512cd sha_ni avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk avx512_bf16 clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin cppc arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif x2avic v_spec_ctrl vnmi avx512vbmi umip pku ospke avx512_vbmi2 gfni vaes vpclmulqdq avx512_vnni avx512_bitalg avx512_vpopcntdq la57 rdpid overflow_recov succor smca fsrm flush_l1d debug_swap

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[avehtari]

Thanks for working on this. While testing this PR, I (with help from Claude) tried an alternative: keep the existing integrator and scale the integrand by a per-observation log shift, which is a simpler change. In addition, the tests found several other defects in develop, so I propose a separate PR with all fixes. I did test also adding the additional fixes on top of this PR, but building on top the log shift approach did produce better results. Comparison with this PR (head bd47838; for double and var it gives the same values as c8b33c4, where I measured):

• wiener_full_lpdf.hpp: +162/−35 (log shift, gradients for small sw and st0, two develop fixes), against +426 in this PR.
• Integration: the existing hcubature, with a fix for two bugs in it (+8/−14), against a second adaptive integrator of about 355 lines. With the two bugs, hcubature runs to the 6000-evaluation limit once the first box fails, and can count a split box twice.
• Speed of the 7-parameter lpdf with gradients (300 ordinary points, default precision, one core): 14 times faster with sw and st0, 11 times with st0 only, 0.95 times with sw only; about 1.2 times faster at the wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 point and scan (precision 1e-8). Most of the gain in the ordinary region comes from the hcubature fix.
• wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 point and scan (y − t0 = 3e-6 to 3e-5, sw = 1e-8, st0 = 1e-3): both give the log density to 5e-11. d/dsw: error at most 5e-10, against −1e8 at 3 of 16 scan points (develop's integration limits) and errors of 1e-5 to 4e-5 at 7 points (develop's end-point formula). Your fallback is accurate at the other 6 points.
• Extreme grid (128 points around the wiener_lpdf full model returns -inf for valid-support inputs in variability branches #3329 locations): finite at 108 points, −inf at 20 points (true log density −4.5e4 to −3.2e7), no exceptions. Errors below 1e-7 at 104 points; at the other 4 (a narrow 2-D peak) errors of 6.4e-4, with the same log-density error in this PR. This PR: finite at 112 points, but 24 more have log-density errors (0.13 to 0.47 at 8 points, 1e-5 to 4e-5 at 16) and 12 have d/dsw errors (develop's integration limits and series tolerance); exceptions at 16 points. The log-density errors come from the absolute tolerance log(1e-12) that the fallback passes to wiener5_density and wiener5_grad_*. With a tolerance relative to the factor in front of the series, params.log_error + fmin(0.0, wiener5_compute_log_error_term(reaction_time, params.a, params.v, relative_start, params.sv)) in wiener7_log_node_value, all 24 errors are below 2e-12 and the gradients there below 2e-9 (tested on this PR's head).
• The gradients for small sw and st0 use the idea of your fallback: differentiation under the integral. The end-point formulas multiply the integration error by 1/sw (1/st0); in develop the d/dsw error at sw = 1e-8 is up to 1.1 at the default precision.
• Other develop defects fixed: 5-parameter lpdf series truncated where the density is small (value errors up to 1e-2 at ordinary points); cancellation in P − F in the 4-parameter lccdf (−inf at 8 of 300 ordinary points) and in the 7-parameter lccdf (−inf below about 1e-16 of P; now finite down to about 1e-308); sign errors in the 4-parameter lcdf gradients (also in the 7-parameter lcdf and lccdf gradients, errors up to 3e-3); 7-parameter lccdf missing the probability mass of trials with t0 + τ > y (log-ccdf errors up to 1.4); d/dsw with an absolute series tolerance when st0 = 0 (relative errors 2e3).

EDIT: new PR is #3442

[martonaronvarga]

Thanks for the detailed comparison. I integrated the directly applicable fallback improvement in 33fbfab: each fallback node now makes the Wiener-5 series tolerance relative to its prefactor. I added a regression where the previous fallback returned -38.6717878 instead of the 50-digit reference -38.699512995; the corrected node matches the reference value and active scores to 1e-10. Local checks pass: mixed 4/4, prim 24/24, rev 5/5.\n\nI did not duplicate the broader hcubature, small-sw/st0, lcdf, and lccdf changes from #3442 here. Those form a coherent alternative implementation and are better reviewed in their existing PR. Given its smaller log-shift design and broader accuracy/performance results, #3442 appears to be the preferred merge candidate; this PR can remain available as an independently tested fallback comparison and then be closed as superseded if maintainers agree.

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.04	0.04	0.99	-0.66% slower
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.17	0.16	1.03	2.98% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.85	1.0	0.22% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.0	0.17% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	0.99	-1.43% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.02	8.1	0.99	-1.07% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.78	3.7	1.02	1.98% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	33.17	32.75	1.01	1.26% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.39	0.97	-3.45% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.94	6.97	1.0	-0.46% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	137.19	131.22	1.05	4.35% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.2	3.22	0.99	-0.55% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.61	0.6	1.01	0.87% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.86	25.7	1.01	0.62% faster
performance.compilation	389.11	385.77	1.01	0.86% faster
Mean result: 1.004145656666831

Jenkins Console Log
Jenkins Build Stages
Commit hash: 44a50912ff3f2e71023ea29a6f91adf7c6b1cf87

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           52 bits physical, 57 bits virtual
CPU(s):                                  192
On-line CPU(s) list:                     0-191
Thread(s) per core:                      2
Core(s) per socket:                      48
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              25
Model:                                   17
Model name:                              AMD EPYC 9474F 48-Core Processor
Stepping:                                1
Frequency boost:                         enabled
CPU MHz:                                 1493.824
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.39
Virtualization:                          AMD-V
L1d cache:                               3 MiB
L1i cache:                               3 MiB
L2 cache:                                96 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-47,96-143
NUMA node1 CPU(s):                       48-95,144-191
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Not affected
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Enhanced / Automatic IBRS; IBPB conditional; STIBP always-on; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Mitigation; Clear CPU buffers
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good amd_lbr_v2 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba perfmon_v2 ibrs ibpb stibp ibrs_enhanced vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a avx512f avx512dq rdseed adx smap avx512ifma clflushopt clwb avx512cd sha_ni avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk avx512_bf16 clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin cppc arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif x2avic v_spec_ctrl vnmi avx512vbmi umip pku ospke avx512_vbmi2 gfni vaes vpclmulqdq avx512_vnni avx512_bitalg avx512_vpopcntdq la57 rdpid overflow_recov succor smca fsrm flush_l1d debug_swap

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin