[jachymb]

Depends on the PR #3443

AI use disclosure: This was developed with the help of Claude Code (Opus 5.5, model ID claude-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

Dense BFGS (algorithm=bfgs, BFGSUpdate_HInv) sets its inverse Hessian estimate to $H = (s^\top y / y^\top y) I$ at a reset and then only applies BFGS updates. In the directions that no step has explored since, $H$ keeps that initial scale. Started where the curvature is large, it underestimates the inverse curvature there by orders of magnitude: the steps in those directions are far too short, and the relative gradient test $g^\top H g / \max(|f|, 1)$, which uses $H$, can stop the run far from the optimum. For $f(x) = \sum_{i=1}^5 (e^{i x_i} - i x_i)$ from $x_i = 5$, develop stops at $f - f^* = 7.6 \cdot 10^3$ with CmdStan's default tolerances.

Before each update except at a reset, $H$ is now multiplied by

$$\tau = \frac{s^\top y}{y^\top H y} \quad \text{if } \tau > 1,$$

Oren–Luenberger sizing restricted to enlarging $H$ (in the spirit of Al-Baali, Comput. Optim. Appl. 9 (1998) 191–203). Enlarging all of $H$ can overestimate the inverse curvature along earlier steps; a later pair along such a direction has $y^\top H y > s^\top y$, is not sized, and the update corrects $H$ there.

Depends on the PR #3443, which rewrites update as a rank-two update; until that is merged, the diff here also shows its commit, and only the last commit belongs to this PR. The sizing comes after $u = H y$ of that form: with $y^\top u$ in place of $y^\top H y$, both $H$ and $u$ are multiplied by $\tau$, and $y^\top u = s^\top y$ is used in $c$.

Intended Effect

Dense BFGS reaches the optimum where develop stops early or crawls after a start in a region of large curvature. Against develop, algorithm=bfgs with CmdStan's defaults on 20 model configurations (GLMs, a hierarchical logistic regression, embedded-Laplace GLMMs and GPs; 30–50 seeds each, 900 runs): 795 instead of 651 runs end at the best log density found, with 0.82 times the gradient evaluations (geometric mean over the configurations; 0.38 on a linear regression, 0.44 on Poisson regressions). On 9 standard test functions with 40 starts each, the evaluations fall on 7 (Beale 541 → 117, Wood 61 → 46), and all 40 runs instead of 32–34 reach the optimum on the extended Rosenbrock and a 30-dimensional quadratic with condition number $10^4$.

How to Verify

New tests, both failing on develop:

• OptimizationBfgsUpdate.BFGSUpdate_HInv_enlarges_small_estimate (bfgs_update_test.cpp): after a reset on curvature $10^4$ along $e_1$ and a pair with curvature 1 along $e_2$, $H e_3 = e_3$ (develop: $10^{-4} e_3$); a later pair with curvature 100 along $e_3$ does not shrink $H$ along $e_1$.
• OptimizationBfgsMinimizerExpSum.dense_bfgs_large_initial_curvature (bfgs_minimizer_test.cpp): the example above with tol_rel_grad = 1e7; it reaches $f = 5$, develop stops at $f - 5 = 7.6 \cdot 10^3$.

Changed expectations, because dense BFGS takes a different path:

• OptimizationBfgsMinimizer.minimize (bfgs_minimizer_test.cpp) stops on the absolute change of $f$ (return code 20) instead of the relative gradient (31);
• ServicesOptimize.rosenbrock (services/optimize/bfgs_test.cpp, Rosenbrock from 0) writes 22 iterates instead of 20, calls the interrupt 21 times instead of 19, and ends at (1.0000112, 1.0000218) instead of (1, 1). The L-BFGS version of the same test ends at (0.99998, 0.99997) on develop.

Existing tests that reach the changed code pass otherwise unchanged: bfgs_update_test.cpp, bfgs_test.cpp, bfgs_exponential_boundary_test.cpp, services/optimize/bfgs_jacobian_test.cpp.

Side Effects

Measured with an in-process A/B harness, this branch including the rank-two commit against develop; I can share the harness. At the defaults the rank-two commit alone leaves every run unchanged. Dense BFGS iterates change; line-search failures: 106 of 900 on both at the defaults, 296 → 299 with tol_rel_obj = tol_rel_grad = 1 (up on the embedded-Laplace, Poisson, categorical and hierarchical logistic models, down on ordered logistic, logistic and negative binomial regression). Evaluations rise by more than 5% here:

setting	branch / develop evaluations
defaults: negative binomial and linear regression with $N = 5000$	1.17, 1.10 (all runs at the best log density on both)
tol_rel_obj = tol_rel_grad = 1: four embedded-Laplace models	1.09–1.38 (rank-two commit alone: 0.80–1.24; all at the best log density on both)
test function: Rosenbrock, 40 starts	1.08

L-BFGS is not affected.

Documentation

Doc comment of BFGSUpdate_HInv::update.

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/)

[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.21% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	1.0	-0.28% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.87	0.85	1.02	1.86% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.01	0.84% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.1	1.03	2.57% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.23	8.32	0.99	-1.09% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.64	3.71	0.98	-1.87% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	33.23	32.79	1.01	1.33% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.39	1.0	-0.23% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.82	6.84	1.0	-0.26% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	168.97	167.01	1.01	1.16% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.25	3.24	1.0	0.17% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.6	1.0	0.22% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.84	24.98	1.03	3.35% faster
performance.compilation	385.1	380.27	1.01	1.25% faster
Mean result: 1.0063677959399362

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:                                 1494.912
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]

Looks good. Couple (Claude assisted) comments

• In ServicesOptimize.rosenbrock (services/optimize/bfgs_test.cpp), I suggest replacing the two new EXPECT_FLOAT_EQ(1.0000112, ...) and EXPECT_FLOAT_EQ(1.0000218, ...) with EXPECT_NEAR(1, ..., 1e-4). The test should check that BFGS ends near the optimum (1, 1), which is what its message "optimal value should be (1, 1)" says. Pinning the exact end point of the current path makes the test fail at every later change of the path, even when the result stays correct.
• Optimizer: update the dense BFGS inverse Hessian in O(n²) instead of O(n³) #3443 now has a commit that keeps only the lower triangle of H (see my comment there). This PR changes the same lines, so it needs a rebase onto it. The sizing carries over: _Hk *= tau scales the stored triangle, _u *= tau is unchanged, and w = (c / 2) s - rho u takes the new c. The Rosenbrock expectations above may then change by rounding, so please run the tests again after the rebase.