[jachymb]

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

Submission Checklist

• Run unit tests: ./runTests.py src/test/unit
• Run cpplint: make cpplint
• Declare copyright holder and open-source license: see below

Summary

stan::model::grad_hess_log_prob estimates the Hessian from gradients at $x \pm \varepsilon e_j$ and $x \pm 2\varepsilon e_j$, $\varepsilon = 10^{-3}$, and symmetrises it:

$$H_{jk} \approx \frac{1}{2\varepsilon} \sum_{i=1}^{4} c_i \left(\partial_k f(x + h_i e_j) + \partial_j f(x + h_i e_k)\right), \quad h = (-2, -1, 1, 2),\varepsilon, \quad c = \left(\tfrac{1}{12}, -\tfrac{2}{3}, \tfrac{2}{3}, -\tfrac{1}{12}\right).$$

The coefficients were multiplied by $\varepsilon/2$ instead of $1/(2\varepsilon)$, so the Hessian came out $\varepsilon^2 = 10^{-6}$ times too small. The only caller in Stan is optimization::newton_step, so its Newton step $-H^{-1}g$ was $10^6$ times too long. newton_step halves the step until the log density does not decrease, which took about 20 halvings (one gradient evaluation each) per iteration, and accepted the first step that did not lower it. On models where Stan Math wrongly returns a log density of $+\infty$ far from the mode (embedded Laplace with a numerically singular covariance, part 2 of stan-dev/math#3469), that first step can land there.

Open PR #3429 changes newton_step to compute its Hessian with stan::math::internal::finite_diff_hessian_auto, which also avoids this bug for the Newton optimizer. This PR fixes grad_hess_log_prob itself, which is public API, and it fixes Newton until #3429 is merged. The two PRs touch different files.

For discussion

Consider deprecating grad_hess_log_prob once #3429 is in. After #3429 nothing in Stan calls it. A GitHub code search finds no other caller apart from vendored copies of the Stan headers and two options of rstan's unexported doctor_cppcode, which rstan does not use. Its Hessian has been $10^6$ times too small since v2.13.0 (2016) without a report. The alternatives are stan::model::hessian (nested autodiff, with propto and jacobian fixed to true) and Stan Math's internal finite_diff_hessian_auto, which laplace_sample already uses; grad_hess_log_prob could also stay as a thin wrapper over the latter. Either way, until it is removed, this fix keeps it from returning a Hessian $10^6$ times too small.

Intended Effect

grad_hess_log_prob returns the Hessian, and the Newton optimizer needs about a tenth of the gradient evaluations. Against develop, algorithm=newton with iter=100 on 7 GLM-type models × 20 seeds and 5 embedded-Laplace models × 5 seeds: 0.06–0.15 times develop's mean evaluations per model (geometric mean 0.09), and all 165 runs end at the best log density found, against 137 on develop. Develop hits the iteration limit in 42 runs and ends at a log density of $+\infty$ in 12 of the 15 Laplace GP runs.

How to Verify

• New test ModelUtil.grad_hess_log_prob_values (src/test/unit/model/grad_hess_log_prob_test.cpp): log density $-x^2/2$ at $x = 1.5$. It fails on develop, which returns the Hessian $-10^{-6}$ instead of $-1$.
• Existing tests that reach the changed code pass unchanged: ModelUtil.grad_hess_log_prob and src/test/unit/services/optimize/newton_test.cpp, newton_jacobian_test.cpp.

Side Effects

Measured with an in-process A/B harness against develop (gradient evaluation counts, platform independent); I can share the harness. Newton iterates change. On the GLM-type models the runs end closer to the mode (median distance to the best log density $10^{-12}$ instead of $4 \cdot 10^{-8}$). grad_hess_log_prob has no other caller in Stan.

Documentation

None.

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):

Jáchym Barvínek jachymb@gmail.com

By submitting this pull request, the copyright holder is agreeing to 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/)

[avehtari]

Excellent work! Crazy this wasn't spotted bewfore

Claude suggested additional test:

The new test has one parameter, so it checks only the diagonal. A test with two parameters also checks the off-diagonal elements and the symmetrization. The existing test model src/test/test-models/good/optimization/rosenbrock.stan has a cross term, and its gradient is a cubic polynomial, so the 4-point differences are exact up to rounding. The test must be in a separate file, because valid.hpp and rosenbrock.hpp both define stan_model, for example src/test/unit/model/grad_hess_log_prob_2d_test.cpp:

#include <stan/model/grad_hess_log_prob.hpp>
#include <stan/io/empty_var_context.hpp>
#include <test/test-models/good/optimization/rosenbrock.hpp>
#include <gtest/gtest.h>

TEST(ModelUtil, grad_hess_log_prob_2d_values) {
  // log density -((1 - x)^2 + 100 (y - x^2)^2) at (x, y) = (0.5, 1):
  // gradient (151, -150), Hessian [[98, 200], [200, -200]].  The gradient is
  // a cubic polynomial, so the 4-point differences are exact up to rounding.
  // The off-diagonal elements check the symmetrized cross terms.
  stan::io::empty_var_context data_var_context;
  stan_model model(data_var_context, 0, static_cast<std::stringstream*>(0));
  std::vector<double> params_r{0.5, 1.0};
  std::vector<int> params_i;
  std::vector<double> gradient, hessian;
  double lp = stan::model::grad_hess_log_prob<true, true, stan_model>(
      model, params_r, params_i, gradient, hessian, 0);
  EXPECT_FLOAT_EQ(-56.5, lp);
  ASSERT_EQ(2, gradient.size());
  EXPECT_FLOAT_EQ(151, gradient[0]);
  EXPECT_FLOAT_EQ(-150, gradient[1]);
  ASSERT_EQ(4, hessian.size());
  EXPECT_NEAR(98, hessian[0], 1e-8);
  EXPECT_NEAR(200, hessian[1], 1e-8);
  EXPECT_NEAR(200, hessian[2], 1e-8);
  EXPECT_NEAR(-200, hessian[3], 1e-8);
}

I compiled this test with googletest 1.11.0 against the develop header and against your header, outside the stan test framework (the model header generated with stanc 2.40.0). With develop, all 4 Hessian elements fail (98e-6, 200e-6, 200e-6, -200e-6). With your header the test passes; the errors are 0.9e-11 to 2.0e-11 at -O0 and -O3, so the tolerance 1e-8 has a margin of about 500.

On the interaction with stan-dev/math#3469: newton_step() treats only std::domain_error as a rejected point, and accepts +inf and NaN because f1 < f0 is false for them. stan-dev/math#3438 makes the solver-3 log determinant throw std::domain_error when det(I + Σ W) is not positive, so with both PRs the Newton optimizer rejects those points instead of stopping at +inf. Part 1 of stan-dev/math#3469 (the value at the initial guess) is not covered by either PR.

[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.69% slower
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.17	0.16	1.02	2.07% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.88	0.98	-2.26% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.0	-0.14% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.1	1.0	0.46% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.05	8.16	0.99	-1.3% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.62	3.76	0.96	-4.05% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.9	32.98	1.0	-0.24% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.39	0.39	1.0	0.04% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.79	6.84	0.99	-0.79% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	170.03	167.75	1.01	1.34% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.23	3.21	1.01	0.61% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.6	0.99	-0.74% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.66	25.79	1.0	-0.48% slower
performance.compilation	371.1	386.55	0.96	-4.16% slower
Mean result: 0.9934290392355323

Jenkins Console Log
Jenkins Build Stages
Commit hash: 6a335dd02ea35aac291d2cb571d03092a8b16ad7

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.887
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