[jachymb]

AI use disclosure: This was developed with the help of Claude Code (Claude Opus 5.5, Claude Fable 5.1).

Summary

The reverse-mode partials of $y = l + (u - l),\sigma(x)$ (σ the inverse logit) use 1 - inv_logit(x): a relative error of 3.6e-8 at x = 20 and exactly 0 for x ≥ 37 (true ∂y/∂x 4.3e-16 there for l = −2, u = 3). This is a first-order error, visible in optimization with jacobian = false. $\sigma(-x)$ now comes from the exponential that $\sigma(x)$ already uses ($e^{-|x|}$ for scalars, the stored $e^x$ for matrices and OpenCL):

$$ \frac{\partial y}{\partial x} = (u - l),\sigma(x),\sigma(-x), \qquad \frac{\partial y}{\partial l} = \sigma(-x) = \frac{1}{1 + e^{x}}, \qquad \frac{\partial y}{\partial u} = \sigma(x). $$

Forward mode has the same defect through the fvar overload of inv_logit; see also #3446 .

Tests

References are closed forms; tolerances are 1e-13 relative, floored at the smallest normal double. Both fail on develop.

• New test/unit/math/rev/constraint/lub_constrain_test.cpp: all reverse-mode overloads (scalar, Matrix<var>, var_value, with and without lp) at x from −800 to 800 against closed forms.
• OpenCL: rev_derivatives_tails in test/unit/math/opencl/rev/constraint/lub_constrain_test.cpp.

Side Effects

Benchmarked against develop in one process, -O3 -DNDEBUG, on MinGW (MSYS2 g++ 15.2) and glibc (WSL g++ 12.2). Time, this branch / develop: 0.87–1.09 over all paths; after re-measurement the scalar var path is 1.06–1.09 and Matrix<var> with lp 1.00–1.07, on both platforms. The fvar<var> Hessian goes through forward-mode inv_logit and is unchanged. Memory: scalar +8 arena bytes per call (64 → 72), matrices +16 bytes per element. OpenCL was not timed.

Release notes

Fixed the reverse-mode derivatives of lub_constrain, inaccurate from x ≈ 20 and zero beyond 37.

Checklist

• Copyright holder: Jáchym Barvínek jachymb@gmail.com

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• 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/gp_regr.stan	0.22	0.23	0.99	-1.45% slower
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.97	-3.2% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.93	-7.91% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.61	0.62	0.99	-0.63% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.18	3.2	1.0	-0.48% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	43.14	43.08	1.0	0.14% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.61	0.96	-3.96% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.45	4.48	0.99	-0.83% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	21.02	22.67	0.93	-7.84% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.89	0.93	0.96	-4.26% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	0.49% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.3	8.31	1.0	-0.11% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.34	6.56	0.97	-3.45% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	170.46	169.46	1.01	0.59% faster
performance.compilation	385.01	388.87	0.99	-1.0% slower
Mean result: 0.9785669609363813

Jenkins Console Log
Jenkins Build Stages
Commit hash: 017f091b63e5f89b01de871d2a7ead2d78127331

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:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1496.698
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.56
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
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:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
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; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
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 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 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 ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

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]

Looks good. I fixed two things (Claude assisted) and there are two comments about PR text

Commit 8cae6d4: no division in the matrix reverse passes

(1.0 + exp_x).inverse() is a lazy Eigen expression, so each use of σ(-x) in a reverse pass divided again. The matrix overloads now do what your scalar overloads already do. They store exp(-|x|) instead of exp(x), compute σ(x) = select(x >= 0, 1, e) / (1 + e) with one division and without the INFTY guard, and take σ(-x) = select(σ(x) >= 1/2, e σ(x), 1 - σ(x)) in the reverse pass. 1 - σ(x) is exact to rounding for σ(x) <= 1/2. The lp overloads reuse the same array: log1p_exp(-|x|) becomes log1p(exp(-|x|)), which saves one exp per element. Arena memory does not change, and scalar and matrix inputs now get the same σ values.

Time relative to your head commit (Triton, Xeon E5-2680 v3, g++ 14.3, -O3 -DNDEBUG, one core, 5 interleaved rounds, median; inputs uniform in [-6, 6]; one call is value, sum, grad and recover_memory):

path	n = 10	100	1000	10000
Matrix<var>, scalar bounds	0.93	0.86	0.91	1.03
the same with lp	0.86	0.79	0.81	0.90
var_value, scalar bounds	1.12	0.87	0.82	0.83
the same with lp	0.94	0.69	0.67	0.67
Matrix<var>, var bounds, lp	0.90	0.80	0.82	0.95
Matrix<var>, vector lb, lp	0.90	0.84	0.87	0.89
Matrix<var>, vector lb and ub, lp	0.89	0.85	0.86	0.92

Relative to develop, the head commit took 1.01 to 1.17 times as long on these paths, and the new commit 0.70 to 1.14. One or two of the five rounds were up to 24 % slow on the shared node, so treat differences below about 5 % as noise. The var_value result at n = 10 (1.12) had a spread of 5 %, so it is probably real, about 30 ns per call.

Commit 456e257: forward-mode derivatives in prim/constraint/lub_constrain.hpp

Forward mode goes through the prim scalar overloads, which this PR did not change. On develop and on your head commit I measured two defects against mpmath (lb = -2, ub = 3):

• Autodiff of (ub - lb) * inv_logit(x) + lb gives the lb tangent as 1 - σ(x): relative error 3.6e-8 at x = 20 and 1.0 from x = 37. The x tangent and the second derivative in x go through the fvar inv_logit, which has the same cancellation until Accurate derivatives of inv_logit and log_inv_logit_diff at large arguments #3446 is merged.
• The log Jacobian uses -abs(x) - 2 log1p_exp(-abs(x)). The fwd abs returns fvar(0, 0) at exactly 0, so the second derivative there is 0 instead of -0.5, in fvar<fvar<double>> and in fvar<var>. x = 0 is where init = 0 starts.

For autodiff types the commit evaluates inv_logit only at a non-positive argument: ub - diff * inv_logit(-x) for x ≥ 0 and lb + diff * inv_logit(x) otherwise, with the branch on value_of_rec(x). That fixes the lb and x tangents and the second derivative without #3446. Arithmetic types keep diff * inv_logit(x) + lb, so double results do not change bitwise (prim/constraint/lub_transform_test compares them exactly). The log Jacobian branches on the sign in the same way instead of calling abs; for doubles that gives the same bits as before.

The new test test/unit/math/mix/constraint/lub_constrain_fwd_test.cpp checks the value, the tangents in lb, ub and x, the second derivative in x, and the first and second derivatives of the log Jacobian at 13 points from -700 to 700, in fvar<double>, fvar<fvar<double>> and fvar<var>. The references are mpmath values at 50 digits; the tolerance is 1e-13 relative, floored at the smallest normal double. With the prim header of your head commit, 42 assertions in both test cases fail: 40 at x ≥ 20 (the lb and x tangents and the second derivatives), and the second derivative of the log Jacobian at x = 0 in both nested modes. I put the test in a new file so that it does not conflict with the lub_constrain_1_test.cpp change in #3446. The two PRs change different files. With both merged into develop e4745c8, all 17 CPU suites that either PR changes or adds pass, including the forward-mode lub_constrain test of #3446.

Suggestions for the PR text and the tests:

1. The PR text says "1e-13 relative, floored at the smallest normal double", but expect_rel in the rev test uses 1e-13 * std::fabs(expected) and the OpenCL test has no floor either. An expected value of 0 (x = ±800, and dlogJ at x = 0) therefore needs an exact 0. These assertions pass today. I would add the floor, as the new forward-mode test does, or change the sentence.

2. The PR text says +16 bytes per element for matrices. Develop already keeps inv_logit_x in the arena and the PR adds one array, so it is 8 bytes. The forward-mode part of the PR text ("Forward mode has the same defect ... see also Accurate derivatives of inv_logit and log_inv_logit_diff at large arguments #3446") can now point to the second commit.

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.23	0.22	1.01	1.09% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.98	-2.24% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.33	6.59	0.96	-4.0% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.0	0.4% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.29	8.26	1.0	0.44% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.5	4.44	1.01	1.42% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	219.18	219.02	1.0	0.08% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.9	0.94	0.95	-4.96% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.63	0.61	1.02	2.44% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	43.39	43.19	1.0	0.44% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.6	0.98	-1.96% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	21.01	22.64	0.93	-7.72% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.01	0.74% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.18	3.17	1.0	0.25% faster
performance.compilation	374.53	378.29	0.99	-1.01% slower
Mean result: 0.9910349884894918

Jenkins Console Log
Jenkins Build Stages
Commit hash: a2236819d562837e97a9bfe97dc9ebd6d8613729

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:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1492.577
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.85
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
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:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
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; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
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 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 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 ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

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