[jachymb]

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

Summary

The derivative of inv_logit is computed as val * (1 - val), which subtracts a number close to 1: it has a relative error of 3.6e-8 at x = 20 and is exactly 0 for x ≥ 37 (true value 8.5e-17). The first-order gradient of a model is accurate in absolute terms, but everything that differentiates it is not: the second derivatives of log1p_exp, log1m_inv_logit and the −W of laplace_marginal_bernoulli_logit_lpmf are 0 there (and those of log_inv_logit for x ≤ −37), and so are the forward-mode derivatives of every function whose fvar overload goes through inv_logit (Phi_approx, lub_constrain, the logistic cdf family). The x partial of log_inv_logit_diff cancels the same way (0 at (40, 0) instead of 8.5e-18).

The rev, fwd and OpenCL overloads now use these forms (rev with the value's exponential):

$$ \begin{aligned} &\sigma'(x) = \sigma(x) \sigma(-x) = \frac{e}{(1 + e)^2}, \quad e = \exp(-|x|), \\ &\frac{\partial \log(\sigma(x) - \sigma(y))}{\partial x} = \sigma(-x) + \frac{1}{\mathrm{expm1}(x - y)}, \qquad \frac{\partial \log(\sigma(x) - \sigma(y))}{\partial y} = -\sigma(y) - \frac{1}{\mathrm{expm1}(x - y)}. \end{aligned} $$

The rev inv_logit is also rewritten in the current callback style, with one overload for var and var_value<Matrix> and one for Eigen::Matrix<var>. The latter computes the value and the derivative in one scalar pass over the elements and stores the derivative for the reverse pass.

It also carries a forward-mode tail test of lub_constrain, whose x tangent and second derivative go through the fvar inv_logit on develop. The lb tangent (1 - σ(x)) and the second derivative of the log Jacobian at x = 0 are separate defects of prim/constraint/lub_constrain.hpp, fixed in #3447 (commit 456e257), which also evaluates inv_logit only at non-positive arguments: #3447 does not depend on this PR, and with it merged the test here is a guard. lub_constrain's reverse-mode fix is also in #3447. Phi_approx is fixed in forward mode the same way; its forward-mode test is on bugfix/phi-approx-stable-derivatives, which is stacked on this branch.

Tests

References are mpmath values (50–60 digits, see test/unit/math/mix/fun/inv_logit_tail_refs.hpp); tolerances are 1e-13 relative, floored at the smallest normal double. Every new test fails on develop unless marked.

• test/unit/math/mix/fun/inv_logit_test.cpp, log1p_exp_test.cpp, log_inv_logit_diff_test.cpp: first and second derivatives for var, var_value, fvar<double>, fvar<fvar<double>> and fvar<var> at x from −800 to 800; ...EqualArgs of log_inv_logit_diff pins the behaviour change below.
• Forward-mode second derivatives through inv_logit: log1mInvLogitFwdDerivativeTails, logInvLogitFwdDerivativeTails (mix/fun/log1m_inv_logit_test.cpp, log_inv_logit_test.cpp) and lub_constrain_fwd_derivative_tails (mix/constraint/lub_constrain_1_test.cpp).
• OpenCL (test/unit/math/opencl/rev/inv_logit_test.cpp, log_inv_logit_diff_test.cpp): rev_derivative_large_args, rev_gradient_large_args against the same references, and rev_non_const_lvalue, which does not compile on develop.

Run locally and passing:

• WSL, g++ 12.2, O=0: at 549b4b5, before the one-pass rewrite of the Matrix<var> overload, the six mix test files above; after the rewrite, mix/fun/inv_logit_test.cpp, mix/prob/student_t_qf_test.cpp (the one other test that instantiates the rewritten overload), rev/prob/ordered_logistic_test.cpp, the expression tests of inv_logit, the header check of rev/fun/inv_logit.hpp, runChecks.py, cpplint and clang-format 10.
• The same source files, as part of the combined branch of all seven fixes (bugfix/numerically-stable-derivatives at a3f1aaa): MSYS2 g++ 15.2 -O1 prim, rev, fwd, mix and Laplace tests of the changed functions and their callers; the expression tests; the OpenCL tests on the Intel OpenCL CPU runtime and on an AMD gfx1035 GPU; clang-format 10 and cpplint.

Not run: make test-headers, make doxygen, the full unit suite. CI passed on the previous head (549b4b5).

Side Effects

Benchmarked against develop in one process, -O3 -DNDEBUG, on MinGW (MSYS2 g++ 15.2) and glibc (WSL g++ 12.2); the Matrix<var> row after the one-pass rewrite of that overload, the other rows on the earlier head, whose code on those paths is unchanged. Time, this branch / develop:

path	this branch / develop
inv_logit of Matrix<var>, n = 10 to 10⁴	0.94–1.11 on glibc, 1.00–1.07 on MinGW (identical code: 0.90–1.11 and 1.00–1.04)
inv_logit of a scalar var	1.02–1.11
lub_constrain fvar<var> Hessian, through forward-mode inv_logit	0.91–1.45 on MinGW, unstable between runs; 0.97–0.99 on glibc
every other path, Hessians included	0.62–1.11

Memory: inv_logit of a scalar var 32 → 40 arena bytes per call, containers +8 bytes per element; the fvar<var> Hessian of log1p_exp +4%. Values within 3 ulp of develop. OpenCL was not timed.

Behaviour changes:

• log_inv_logit_diff(x, y) at x = y: ∂x is +∞ (the limit from x > y) instead of −∞.
• OpenCL inv_logit of a non-const lvalue var_value<matrix_cl<double>> now compiles.

Release notes

Fixed the derivatives of inv_logit and log_inv_logit_diff, 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/gen_gp_data.stan	0.04	0.04	1.02	1.89% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.17	0.16	1.01	0.68% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.9	0.95	-4.78% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	0.97	-2.87% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.1	0.11	0.95	-4.97% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.01	8.04	1.0	-0.32% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.86	3.7	1.04	4.14% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.2	33.16	0.97	-2.98% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	0.99	-0.9% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.93	6.85	1.01	1.18% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	137.12	135.34	1.01	1.3% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.22	3.22	1.0	0.0% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.61	0.61	1.0	0.04% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.54	25.16	1.02	1.51% faster
performance.compilation	399.28	387.3	1.03	3.0% faster
Mean result: 0.998621302648084

Jenkins Console Log
Jenkins Build Stages
Commit hash: a827cdd3e0feedf2727b243d86681a6a8e0fd2ba

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.762
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]

I reviwed this with help from Claude

1. Speed of inv_logit of Matrix<var>. I see the same slowdown on glibc. Setup: Triton, Xeon E5-2680 v3, g++ 14.3, -O3 -DNDEBUG, no -march, one core

n	this PR	store d = e/(1+e)^2 instead of e	scalar loop below
10	1.22	1.15	1.05
100	1.22	1.14	1.10
1000	1.24	1.21	1.19
10000	1.33	1.33	1.19

The var_value<VectorXd> path is at 0.99 to 1.04 and the scalar var path at 1.01, so only Matrix<var> is slower. The PR text attributes the glibc result to Eigen's vectorized exp. Develop also computes the values with vectorized Eigen code, logistic() in prim/fun/inv_logit.hpp. My guess is that the extra arena array causes the difference, because none of the variants that store one reaches the develop time. The fastest variant computes the value and d in one loop and does not copy x:

auto x_arena = to_arena(std::forward<T>(x));
const Eigen::Index n = x_arena.size();
arena_t<Eigen::ArrayXd> d(n);
arena_t<T> ret(n);
for (Eigen::Index i = 0; i < n; ++i) {
  const double xi = x_arena.coeff(i).val();
  const double e = std::exp(-std::fabs(xi));
  ret.coeffRef(i) = xi < 0 ? (xi < LOG_EPSILON ? e : e / (1.0 + e)) : 1.0 / (1.0 + e);
  d.coeffRef(i) = e / ((1.0 + e) * (1.0 + e));
}
reverse_pass_callback([x_arena, ret, d]() mutable {
  x_arena.adj().array() += ret.adj().array() * d;
});

For var_value, a reverse pass without a division ran at 0.92 to 0.94 of develop for n ≥ 100. It keeps e and computes v * (v >= 0.5).select(e * v, 1.0 - v) with v = res.val().array(). For v <= 1/2, 1 - v is exact to rounding, and the worst adjoint error was 5.0e-16. So there is a possibility for 20% speed-up

2. PR text. It lists PhiApproxFwdDerivativeTails in a new file mix/fun/Phi_approx_fwd_test.cpp. Commit 32e48bd moved that test to the Phi_approx branch, and the file is not in the diff.

3. Forward-mode lub_constrain. The PR text says that this PR fixes lub_constrain in forward mode through inv_logit. That holds for d/dx. Two more forward-mode defects were in prim/constraint/lub_constrain.hpp: the fvar tangent with respect to lb was computed as 1 - σ(x) (relative error 3.6e-8 at x = 20, 1.0 from x = 37), and at exactly x = 0 the nested second derivative of the log Jacobian was 0 instead of -0.5, because the fwd abs returns fvar(0, 0) there. I pushed a fix for both to Accurate lub_constrain derivatives for large positive arguments #3447 (commit 456e257, details there). It evaluates inv_logit only at non-positive arguments, so it does not depend on this PR. With both PRs merged into develop e4745c8, your lub_constrain_fwd_derivative_tails test and the other mix tests that this PR changes pass, as do the lub_constrain suites of Accurate lub_constrain derivatives for large positive arguments #3447.

[jachymb]

@avehtari Do I understand your review correctly that it's just commentary without a request for changes?

[avehtari]

Sorry it was not clear that they are not just comments, but suggestions for changes

1. suggests 20% speed-up which you could add, but I left it for you to decide whether it is worth it.
2. suggests to edit the PR text
3. suggests to edit the PR text

They are not show-stopper requests

[jachymb]

@avehtari I see. Sometimes the LLM comments are challenging to parse properly, but I'm very guilty of producing that myself too 😅

Anyway, thanks for the review. I checked that the case 1 is valid and then adopted it. Adjusted the PR text based on notes in 2, 3.

[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.23	1.0	0.3% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	1.0	0.11% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.34	6.32	1.0	0.37% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.99	-1.04% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.28	8.28	1.0	0.02% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.47	4.46	1.0	0.18% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	218.41	217.9	1.0	0.23% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.9	0.9	1.01	0.6% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.62	0.61	1.01	1.13% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	43.58	42.76	1.02	1.87% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.59	0.99	-0.68% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.98	20.99	1.0	-0.04% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	-0.48% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.19	3.18	1.01	0.5% faster
performance.compilation	361.19	372.2	0.97	-3.05% slower
Mean result: 1.0001265973234728

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:                                 1496.184
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