Numerically stable derivatives of ordered_logistic_lpmf and ordered_logistic_glm_lpmf

[jachymb]

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

Summary

Follow-up to #3327 by @avehtari . At a location equal to a cut point the first-order partials are correct, but abs() has derivative 0 there, so the value's own derivative and the Hessian are wrong (cuts (0, 1), x = 1, y = 1, β = 0: −0.5 and −0.25 become 0). Two more defects are independent of the cut-point case. When two adjacent cut points are only Δ apart, the category between them has probability of about Δ, so the partial derivatives with respect to those cuts are of size 1/Δ. develop computes them as the difference of two nearly equal terms of that size, and the digits cancel: in the test case with Δ = 2⁻³⁰ ≈ 9.3e-10 and five observations in the narrow category, the cut partial is −5368703998.36 (about −5/Δ), but develop returns −5368703995.86, correct to 9 of its 16 digits. Forward mode and the OpenCL kernels lose the same digits. And a very large location overflows the value: ordered_logistic_lpmf(2, 1e17, [0, 1]') returns -inf instead of its finite value of about −1e17.

The new code uses one sign partition for the value and all partials, takes Δ from the cut points and uses these forms, whose cut partials add same-sign terms (σ is the inverse logit):

$$ \begin{aligned} &a = \theta - c_{k-1}, \quad b = \theta - c_k, \quad \Delta = c_k - c_{k-1}, \quad q = \frac{e^{-\Delta}}{\mathrm{expm1}(-\Delta)} \le 0, \\ &\log p = \log\sigma(a) + \log\sigma(-b) + \log(1 - e^{-\Delta}), \\ &\frac{\partial \log p}{\partial \theta} = \sigma(-a) - \sigma(b), \quad \frac{\partial \log p}{\partial c_{k-1}} = q - \sigma(-a), \quad \frac{\partial \log p}{\partial c_k} = \sigma(b) - q. \end{aligned} $$

The OpenCL kernels of both functions compute the same forms.

Includes a merge of #3418, which I recommend merging first. Until it is, the diff shown here includes its changes.

Tests

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

• New test/unit/math/mix/prob/ordered_logistic_test_helpers.hpp: value and first to third derivatives against closed forms at cut points, ±40, ±750, close cuts, far locations and one cut, used by ordered_logistic_glm_lpmf_test.cpp and two tests in ordered_logistic_test.cpp. They fail on develop.
• No new OpenCL test; the existing OpenCL tests compare with the CPU.

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.49–1.06 over all paths, Hessians 0.72–0.90; no path is above 1.06 after re-measurement. Arena memory per call and the tape nodes of an fvar<var> sweep are lower on every path. Values within 1 ulp of develop. OpenCL was not timed; the kernels gain one expm1 per element with autodiff cuts.

Release notes

Fixed ordered_logistic_lpmf and ordered_logistic_glm_lpmf for large locations, close cut points and derivatives at a cut point.

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

[jachymb]

Jenkins getting OOMkilled likely for the same reason #3418 did. Let's merge that first and then rerun.

[WardBrian]

We can also increase the memory available to that branch, but I don’t think we can change Jenkinsfiles from fork branches so we’d just have to do it in develop and then merge that in here

[jachymb]

@WardBrian I think after merging #3418 this should be able to rerun without a memory increase (because this includes the changes from #3418 )

[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.01	0.73% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	1.01	0.73% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.35	6.35	1.0	0.09% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.95	-5.28% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.3	8.31	1.0	-0.08% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.48	4.49	1.0	-0.23% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	168.33	200.53	0.84	-19.13% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.89	0.89	1.0	0.15% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.62	0.62	1.0	-0.35% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	43.16	43.26	1.0	-0.22% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.6	0.98	-1.98% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	21.0	20.96	1.0	0.18% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	0.99	-1.23% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.19	3.18	1.0	0.41% faster
performance.compilation	392.8	392.38	1.0	0.11% faster
Mean result: 0.9848641148146273

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

[jachymb]

Passed after merging #3418 as predicted :)