[jachymb]

Fixes #3469. The diff is only optically large: most of it is a block of code moved into a helper and re-indentation, so hiding whitespace in the diff view helps. Stacked on #3465 (its commit shows in the diff until it merges); suggested merge order: #3438, #3465, this PR.

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

Summary

laplace_marginal_density_est could return the Laplace approximation evaluated at the initial guess as the log marginal, or +inf / NaN, without a message (#3469). Three changes:

1. wolfe_line_search no longer accepts the zero step. The fallback best starts at step 0, which satisfies Armijo trivially, so a failed search was reported as an accepted step and the Newton loop "converged" at the starting point. The fallback now requires best.alpha() > 0.
2. run_newton_loop searches shorter steps when no step improves the objective. From θ = 0, W can be many orders of magnitude smaller than at the mode, so the Newton step $p$ overshoots by more than $1/\alpha_{\mathrm{min}}$ (min_alpha = 1e-8) and every trial step lowers the objective or overflows. The loop now retries the search with $p$ scaled by $\alpha_{\mathrm{min}}$. For a concave objective $f$, steps that short gain at most
   $$f(a + t p) - f(a) \le t \nabla f(a)^\top p \le \alpha_{\mathrm{min}} \nabla f(a)^\top p, \qquad 0 < t \le \alpha_{\mathrm{min}},$$
   so the retry stops when this bound is below tolerance or below the rounding error of $f$, or is not finite, or min_alpha ≥ 1.
3. A singular $I + \Sigma W$ throws std::domain_error. Mathematically $\det(I + W^{1/2} \Sigma W^{1/2}) \ge 1$, but when $\Sigma$ is numerically rank deficient and $\Sigma W$ swamps the identity, the matrix is singular in floating point: in the example at α = 1e9 the jitter is below the ulp of α², so all rows of $I + \Sigma W$ are equal. Then the Newton system has no solution and the LU factorization of solver 3 has zero pivots, so its log determinant is −∞ and develop returned +inf (NaN with a negative pivot). Solvers 1 and 2 already throw when their Cholesky fails in this situation; solver 3 now throws at the first factorization with a zero pivot, instead of iterating on a meaningless step and reporting the result at the end. A check after the solver chain also throws for any other non-finite log marginal (for example a non-finite log likelihood at the final iterate); it sits after the chain so that the chain cannot catch it and fall through into the next solver. Return the final Newton iterate as the Laplace latent mean, and fix the solver-3 log determinant #3438 makes the LU log determinant throw for a non-positive determinant at the end of the loop; the two merge cleanly and complement each other.

Example from the issue (y = (3, 1, 4, 1, 5), K = α² J + 1e-8 I); the reference values are Laplace approximations at the mode in 80-digit arithmetic (mpmath):

μ	α	develop	this PR	reference
−30	3e4	−429.757515545	−20.9711152603	−20.971115053
−30	1e5	−429.759638992	−22.1750801284	−22.1750873716
−30	1e7	−431.690693525	−26.7458217424	−26.7802575096
−30	1e9	−436.285290904	throws	−31.3854276956
1	1e9	+inf	throws	−31.3854276956

The develop values at μ = −30 are the approximation evaluated at θ = 0. At α = 1e9, $I + \Sigma W$ is singular in double precision at the mode (the jitter is below the ulp of α²), so throwing is the correct outcome. The remaining 0.034 error at α = 1e7 is not addressed here. On the N = 100 Poisson GP model where this was found (ρ = 3e27, μ = −51.6), α = 1e4 gives −311.555 instead of −21782.9, the value at θ = 0; from α = 1e8 on, where develop returned that value, +inf or NaN, the result now throws.

The line search of run_newton_loop now lives in internal::line_search_fails, so the retry is a short loop, and the solver chain is wrapped in a lambda so that the check runs after it. An alternative to the lambda, if preferred: split laplace_marginal_density_est into an _impl holding the chain and a thin wrapper holding the check.

Tests

Three new tests, all failing on develop:

• laplace_marginal_poisson_log_lpmf_newton.step_overshoots_mode compares μ = −30, α = 3e4 and 1e5 with the mpmath references (tolerance 1e-4, since the loop stops within its tolerance of the mode, where the log-determinant term still moves the value by up to 1e-5). Develop is off by about 408 in both. α = 1e5 needs change 2; α = 3e4 needs changes 1 and 2 together.
• laplace_marginal_poisson_log_lpmf_newton.singular_log_determinant expects std::domain_error at μ = 1, α = 1e9 with solver 3 and no fallthrough, where all rows of $I + \Sigma W$ are equal and the factorization has zero pivots. Develop returns +inf.
• WolfeLineSearch.NoImprovingStepFailsInsteadOfZeroStep scripts an update where every trial step lowers the objective and expects accept_ == false with status Fail. Develop returns Armijo with accept_ == true, the zero step.

All 18 files in test/unit/math/laplace pass at O=3; the header checks of both changed headers, cpplint and clang-format pass.

Side Effects

Well-conditioned problems are unaffected: where the line search succeeds, values and cost equal develop's. In ill-conditioned cases the number of likelihood evaluations can go up as well as down: up where the first Newton step overshoots and develop stopped at the starting point, since the retry searches on to the mode; down where $I + \Sigma W$ is singular, since solver 3 now fails at the first factorization.

Counts are likelihood calls (WSL, g++, -O0); otherwise the change adds scalar comparisons only, and memory is unchanged. Rows where develop's line search failed or the branch throws (all others are equal):

μ = −30	develop	this PR
α = 3e4	118	177
α = 1e5	57	187
α = 1e9	57 (value at θ = 0)	121 (throws)
α = 1e100	58 (+inf)	4 (throws)

Behaviour changes:

• a failed line search no longer counts as convergence (same returned iterate, one factorization pass fewer);
• a singular $I + \Sigma W$ at any iteration and any non-finite log marginal throw instead of being returned, also in the laplace_latent_* functions (before Return the final Newton iterate as the Laplace latent mean, and fix the solver-3 log determinant #3438 also for a negative LU pivot of a non-singular B);
• with max_steps_line_search = 0, μ = −30 and α = 1e5 or 1e7 throw instead of returning the value at θ = 0: the shorter step found by the retry is accepted blindly and lands where the factorization is singular.

Release notes

Embedded Laplace now finds the mode when the first Newton step overshoots, and throws instead of returning +inf or NaN.

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.01	0.87% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	0.97	-3.08% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.87	0.86	1.01	1.08% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	0.97	-3.13% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.1	0.11	0.98	-1.95% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.03	8.16	0.98	-1.6% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.67	3.66	1.0	0.23% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.76	32.74	1.0	0.06% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.39	1.0	-0.38% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.8	6.87	0.99	-0.94% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	171.31	166.38	1.03	2.88% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.24	3.23	1.0	0.17% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.61	0.99	-1.34% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.77	25.78	1.0	-0.03% slower
performance.compilation	366.85	387.35	0.95	-5.59% slower
Mean result: 0.991960860535368

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