fix integrate_1d_gauss_kronrod gradient speed

[avehtari]

Making this PR has been assisted by Claude. I have checked and edited all new code, comments, and PR text

Summary

The gradient computation for integrate_1d_gauss_kronrod calls the Boost adaptive Gauss Kronrod again, but Boost computes each parameter's gradient integral with the same Boost adaptive Gauss–Kronrod call. Boost's bisection criterion is relative to the leaf estimate, so a gradient component that is tiny relative to f and dominated by round-off (e.g. the score of a saturated observation deep in a likelihood tail) never converges and bisects to max_depth (2¹⁵ × 21 evaluations) on every call, although its contribution to ∂ log I is negligible. This can be fixed by shifting the integral integrand as
g_i(x) = ∂f/∂θ_i (x) + c f(x)
and returning ∫ g_i − c I, where I is the value integral already computed and c != 0 is a constant
such that ∂f/∂θ_i ≈ −c·f does not occur for common integrands; c = 1 is a bad choice because log-density gradients of logit-type likelihoods saturate at ±1 in the tails, which turns the shifted integrand back into noise. The inverse golden ratio is far from every integer and is not a structural constant of any link.

As this is needed only for integrate_1d_gauss_kronrod, integrate_1d_adjoint gets an option whether the shift is used. The plain integrate_1d and integrate_1d_double_exponential are unchanged: their L1-relative termination handles ordinary gradient integrands, and enabling the shift there costs 7–32 % more gradient evaluations on 10 of 35 suite integrands with no gain.

Tests

Measured with the fix (same model, same point):

	gradient	vs integrate_1d	max |Δ grad| vs integrate_1d
GK, unpatched	7582 ms	134×	2.9e-11
GK, patched	156 ms	2.75×	2.9e-11

49× faster gradients, results unchanged to 1e-11 (patched vs unpatched
GK: 1.0e-11). Per subject the counter shows 668 850 → 2 478 gradient
evaluations. The residual 2.75× is GK's per-integral cost on this
integrand (value-only GQ is 2.1× DE), not a gradient effect.

Effect on every integrand in the existing GK reverse-mode test suite

Integrand evaluations (value alone; all gradient components together)
and wall for value + gradients (mean of 200, -O3), unpatched →
patched. Same arguments and tolerances as the tests. #var = number
of var inputs (parameters + var endpoints).

integrand (test)	#var	value	gradients unpatched → patched	wall µs
f1 exp(x)+θ, [0.2, 0.7]	1	21	21 → 21	1.3 → 1.2
f2 exp(θ cos 2πx)+θ, var endpoints	3	21	23 → 23	2.9 → 2.9
f3 exp(x)+θ₀^2.5+2θ₁³+2θ₂, var endpoints	5	21	65 → 65	11 → 7.6
f3, params only	3	21	63 → 63	11 → 5.7
f10 1/(1+x⁴/1.7), (−∞, ∞), no params	0	147	0 → 0	6.3 → 2.7
f11 beta(5, 3), [0, 1]	2	21	294 → 336	29 → 33
f11 beta(3, 5), [0, 1]	2	21	294 → 336	20 → 23
f12 N(5.7, 1), (−∞, ∞)	2	357	672 → 672	44 → 42
f13 x+θ₀+θ₁, endpoints = params	4	21	44 → 44	1.2 → 1.2
Cauchy pdf, (−∞, ∞)	2	231	546 → 462	38 → 33
Exponential pdf, [0, ∞)	1	105	231 → 147	11 → 7.6
Normal pdf, (−∞, ∞)	2	441	1008 → 924	57 → 53
Student-t pdf, (−∞, ∞)	3	651	3129 → 2205	323 → 241
Uniform pdf, endpoints = params	4	21	44 → 44	2.5 → 2.3
noisy gradient (new test), [−1, 1]	1	21	throws → 21	— → 2.6
odd, zero integral (new test)	1	21	(n/a) → 21	— → 0.8
negative integral (new test)	1	21	(n/a) → 21	— → 0.8

Unpatched, the noisy-gradient case does not merely hit max_depth: it
then fails Stan's post-check (error estimate 2.4e-13 exceeds max(rel_tol × L1, abs_tol) with L1 ≈ 0) and throws, i.e. in a
model block every proposal is rejected. Elsewhere the change is
neutral or a mild gain (heavy-tailed pdfs 15–35% fewer gradient
evaluations); the one regression is f11 (beta kernel: the gradient
integrands carry a log x factor that is easier than f on its own),
+14 % gradient evaluations at unchanged wall. All 17 pre-existing
tests pass at their tolerances.

A gradient component that is analytically zero but autodiffs to
round-off noise reproduces it without any likelihood machinery
(cos² + sin² = 1). Added to test/unit/math/rev/functor/integrate_1d_gauss_kronrod_test.cpp
as GradientShift_noisy_gradient (evaluation-count bound + accuracy;
fails on the unpatched code, passes on the patched) together with
GradientShift_guards (exactly-zero integral → shift disabled; negative
integral → shift active), and the existing 17 tests pass unchanged.

Side Effects

1. Cost when a gradient integrand is easier than f. `f = spike(x)
   • θ`: the gradient integrand is the constant 1 and took 21
     evaluations; shifted it re-resolves the spike, 189. Bound: gradient
     cost ≤ (p + 1) × value cost, against 2^15 × 21 per component before.
2. Absolute, not relative, gradient accuracy. The shifted gradient
   inherits the value integral's absolute error (≲ tol × L1(f)), so a
   component with |∂I/∂θ| ≲ tol × L1(f) loses relative precision.
   Measured (f = spike + ε θ x, spike width 0.01, tol = 1e-6):
   unshifted exact; shifted absolute error 3e-19 … 3e-18, i.e. relative
   6e-15 at ε = 1e-3 down to 5e-10 at ε = 1e-9 — round-off level in
   practice, because leaves that pass the value's test are usually far
   below tol. For HMC this is the right accuracy: value and gradient
   enter the Hamiltonian at the same absolute scale.

Release notes

integrate_1d_gauss_kronrod gradient computation was modified to have much improved worst performance timing

Checklist

• Copyright holder: (fill in copyright holder information)

Aki Vehtari

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

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

• the new changes are tested

[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.0	0.13% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	0.99	-1.44% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.87	0.99	-0.94% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	0.97	-3.14% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.1	1.02	2.2% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.34	8.37	1.0	-0.4% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.78	3.77	1.0	0.29% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	33.11	33.36	0.99	-0.76% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	1.0	-0.14% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.61	6.62	1.0	-0.12% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	135.01	134.67	1.0	0.26% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.2	3.24	0.99	-1.31% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.59	0.59	1.0	-0.43% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	24.95	25.03	1.0	-0.32% slower
performance.compilation	403.38	385.73	1.05	4.37% faster
Mean result: 0.9990894014826632

Jenkins Console Log
Jenkins Build Stages
Commit hash: 7e5ee47a0634cc5f5f59c91e7e28b9f4519385fb

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:                                 1495.986
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.04
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 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

[SteveBronder]

The inverse golden ratio is far from every integer and is not a structural constant of any link.

Magical constants like this are a weird agent smell. Give me a minute to read into this more

[avehtari]

Magical constants like this are a weird agent smell. Give me a minute to read into this more

The justification why the constant should not be 1/2 or 1 is clear ,as they match saturated gradient for some models and link functions. There are other constants that could be used. This really as an ad hoc fix, but unfortunately the better fix would require changing Boost code to add absolute tolerance threshold.