[avehtari]

Summary

Two independent fixes in the embedded Laplace approximation.

This PR was Claude assisted

1. The latent mean is the final Newton iterate (stan/math/mix/functor/laplace_base_rng.hpp). laplace_base_rng() computed the mean as covariance_train * md_est.theta_grad, that is K * grad log p(y | theta_hat). This equals theta_hat only at an exact mode. Before exact convergence the two differ by K * grad log p(theta | y). The stopping rule bounds the change of the objective, not this difference, and the factor K makes the difference large when the posterior is much narrower than the prior. In example A with tolerance 1e-6, the final iterate is 0.0012 posterior sd from the mode, but the returned mean was 0.53 sd from it, at a log density 7.1 lower. The mean is now md_est.theta. This affects laplace_latent_solve[_tol]() and the mean of the draws of laplace_latent[_tol]_rng() and laplace_latent[_tol]_{bernoulli_logit,poisson_log,neg_binomial_2_log}_rng().

2. The log determinant of solver 3 uses |U_ii| and checks the sign (LUSolver::compute_log_determinant() in stan/math/mix/functor/laplace_marginal_density_estimator.hpp). Solver 3 factorizes the non-symmetric B = I + K W with partial pivoting, so U can have negative diagonal entries when det(B) > 0. The old sum(log(diag(U))) was then NaN, and laplace_marginal*() returned NaN without a message, with solver 3 or after a fallthrough to solver 3, also for a Gaussian likelihood (example B). The log determinant is now sum(log(abs(diag(U)))), and the sign of det(B) is the permutation sign times the signs of the U_ii. Since det(B) = det(K) det(K^-1 + W) and det(K) > 0, a non-positive det(B) means that K^-1 + W is not positive definite at the final iterate (for example at a saddle point), where the Laplace approximation is not defined. The function then throws std::domain_error with a message that names this cause. (A positive det(B) does not prove positive definiteness: the check finds only an odd number of negative eigenvalues.) compute_log_determinant() also runs in laplace_latent_solve*() and laplace_latent*_rng(), so they throw the same message (example C).

Examples

• A: the longitudinal part of the joint model for one group: 12 observations y ~ normal(m, sigma * m) with m(t) = bsld * (phi * exp(-d * t) + (1 - phi) * exp(g * t)), and 4 latent variables on the log or logit scale of (bsld, d, g, phi) with diagonal K. The posterior sd is 5.6 to 9.4 times smaller than the prior sd. Call: laplace_latent_solve_tol() with solver 3. The reference mode is from optim() followed by Newton steps.
• B: y ~ normal(theta, s) with y = (0.3, -0.2), s = sqrt(0.1), K = [1 1.5; 1.5 4]. det(B) = 226, but the LU diagonal is (15, -15.07) because the rows are swapped. The exact log marginal is multi_normal_lpdf(y | 0, K + s^2 I).
• C: 2 latent variables, K = [4 1; 1 1]. theta[1] has two Cauchy observations y = (-3, 3), s = 0.3; theta[2] has one Gaussian observation z = 0, sz = 0.1. The log posterior has two maxima and a saddle point at (0, 0). From the start (0.01, 3), the Newton iteration of solver 3 converges to the saddle point.

example	call	develop	this PR
A, tol 1e-6	laplace_latent_solve_tol()	0.53 sd, 7.12 lower	0.0012 sd, 2.3e-6 lower
A, default tol	laplace_latent_solve_tol()	0.068 sd, 0.135 lower	0.00015 sd, 3.8e-8 lower
B	laplace_marginal_tol(), solver 3	NaN	-2.376754 (exact)
B	laplace_marginal_tol(), solvers 1, 2	-2.376754	-2.376754
C	laplace_marginal_tol()	NaN	new exception
C	laplace_latent_solve_tol()	cholesky_decompose exception	new exception

For A, "sd" is the largest distance from the reference mode in posterior sd units, and "lower" is the drop in log density from the mode. The new exception is det(I + Sigma * W) is not positive at the final iterate, .... The develop column is from CmdStan 2.40.0; A and B give the same values with develop 43bc568. The "this PR" column is from CmdStan develop with this branch.

R script (cmdstanr; the Stan programs are inline)

# pr_laplace_mean_lu_logdet.R
#
# Examples for a stan-dev/math PR (self-contained; the Stan programs are
# inline):
#
#   A. laplace_latent_solve_tol() returned K * grad log p(y | theta_hat)
#      instead of the final Newton iterate theta_hat.  Before exact
#      convergence the two differ by K * grad log p(theta | y), which is
#      large when the posterior is much narrower than the prior.
#   B. laplace_marginal_tol() with solver 3 returned NaN when the LU
#      factorization of B = I + K W has a negative pivot (also for a
#      Gaussian likelihood).
#   C. At a saddle point of log p(theta | y), the solver-3 log determinant
#      now reports that det(B) is not positive.
#
#   Rscript tests/pr_laplace_mean_lu_logdet.R
#   CMDSTAN_PATH=~/proj/cmdstan STAN_MATH=~/proj/stan_math_laplace_mean/ \
#     Rscript tests/pr_laplace_mean_lu_logdet.R

library(cmdstanr)
library(posterior)
if (nzchar(Sys.getenv("CMDSTAN_PATH"))) {
  set_cmdstan_path(Sys.getenv("CMDSTAN_PATH"))
}
cpp <- list()
if (nzchar(Sys.getenv("STAN_MATH"))) {
  cpp$MATH <- Sys.getenv("STAN_MATH")
}
compile <- function(code) {
  cmdstan_model(write_stan_file(code), cpp_options = cpp,
                force_recompile = length(cpp) > 0)
}
cat(sprintf("CmdStan %s at %s, Stan Math %s\n", cmdstan_version(),
            cmdstan_path(), if (is.null(cpp$MATH)) "of the CmdStan" else cpp$MATH))

# One fixed_param draw; returns the values of `vars` (NA if the draw failed)
# and the "Exception" lines of the output.
run_once <- function(mod, data, vars) {
  out <- capture.output(
    fit <- mod$sample(data = data, fixed_param = TRUE, chains = 1,
                      iter_sampling = 1, refresh = 0, show_messages = TRUE,
                      sig_figs = 12),
    type = "output")
  vals <- tryCatch(as.numeric(as_draws_matrix(fit$draws(vars))),
                   error = function(e) NA_real_)
  list(vals = vals,
       exception = unique(sub(".*Exception: (.*?) \\(in '.*", "\\1",
                              grep("Exception", out, value = TRUE),
                              perl = TRUE)))
}

# ── A. Returned mean ───────────────────────────────────────────────────────
tumor_code <- "
functions {
  real ll(vector theta, vector y, vector t, vector mu, real sigma) {
    real bsld = mu[1] * exp(theta[1]);
    real d = mu[2] * exp(theta[2]);
    real g = mu[3] * exp(theta[3]);
    real phi = inv_logit(logit(mu[4]) + theta[4]);
    vector[num_elements(y)] m
        = bsld * (phi * exp(-d * t) + (1 - phi) * exp(g * t));
    return normal_lpdf(y | m, sigma * m);
  }
  matrix k_fun(vector omega) {
    return diag_matrix(square(omega));
  }
}
data {
  int<lower=1> n;
  vector[n] y;
  vector[n] t;
  vector[4] mu;
  vector<lower=0>[4] omega;
  real<lower=0> sigma;
  vector[4] theta0;
  real<lower=0> tol;
}
generated quantities {
  vector[4] mode;
  {
    tuple(vector[4], matrix[4, 4]) lap = laplace_latent_solve_tol(
      ll, (y, t, mu, sigma), 4, k_fun, (omega,),
      (theta0, tol, 100, 3, 100, 0));
    mode = lap.1;
  }
}
"
y <- c(201.04, 78.30, 48.01, 44.39, 52.91, 50.49, 48.63, 67.29, 52.21,
       72.43, 102.43, 105.06)
t <- seq(0, 693, by = 63)
mu <- c(117.4, 9.121e-3, 1.913e-3, 0.3781)
omega <- c(0.4925, 0.9217, 1.120, 1.299)
sigma <- 0.09826
log_f <- function(th) {
  p <- plogis(qlogis(mu[4]) + th[4])
  m <- mu[1] * exp(th[1]) * (p * exp(-mu[2] * exp(th[2]) * t) +
                               (1 - p) * exp(mu[3] * exp(th[3]) * t))
  sum(dnorm(y, m, sigma * m, log = TRUE)) + sum(dnorm(th, 0, omega, log = TRUE))
}
mode_ref <- optim(rep(0, 4), function(z) -log_f(z), method = "BFGS",
                  control = list(reltol = 1e-15, maxit = 5000))$par
for (it in 1:5) {
  mode_ref <- mode_ref - solve(numDeriv::hessian(log_f, mode_ref),
                               numDeriv::grad(log_f, mode_ref))
}
post_sd <- sqrt(diag(solve(-numDeriv::hessian(log_f, mode_ref))))
cat("\nA. Returned mean.  Mode (optim + Newton):",
    sprintf("%.4f", mode_ref), "\n")
cat("   posterior sd:", sprintf("%.3f", post_sd), " prior sd:", omega, "\n")
mod_a <- compile(tumor_code)
for (tol in c(1e-6, 1.49012e-8, 1e-10)) {
  r <- run_once(mod_a, list(n = length(y), y = y, t = t, mu = mu,
                            omega = omega, sigma = sigma,
                            theta0 = c(0.3636, -0.3385, 1.9794, -1.9322),
                            tol = tol), "mode")
  cat(sprintf("   tol %-9g returned %s: %.2g posterior sd from the mode, log density %.3g lower\n",
              tol, paste(sprintf("%.4f", r$vals), collapse = " "),
              max(abs(r$vals - mode_ref) / post_sd),
              log_f(mode_ref) - log_f(r$vals)))
}

# ── B. Solver 3 with a negative LU pivot ───────────────────────────────────
normal_code <- "
functions {
  real ll(vector theta, vector y, real s) {
    return normal_lpdf(y | theta, s);
  }
  matrix k_fun(matrix K) {
    return K;
  }
}
data {
  vector[2] y;
  real<lower=0> s;
  matrix[2, 2] K;
  int<lower=1, upper=3> solver;
}
generated quantities {
  real lml = laplace_marginal_tol(ll, (y, s), 2, k_fun, (K,),
                                  (rep_vector(0.0, 2), 1e-10, 100, solver,
                                   100, 0));
  real lml_exact = multi_normal_lpdf(y | rep_vector(0.0, 2),
                                     K + diag_matrix(rep_vector(square(s), 2)));
}
"
K_b <- matrix(c(1, 1.5, 1.5, 4), 2)
s_b <- sqrt(0.1)
B <- diag(2) + K_b %*% diag(1 / s_b^2, 2)
cat(sprintf("\nB. Negative LU pivot.  B = I + K W: det %.4g, LU diagonal (partial pivoting) %s\n",
            det(B), paste(sprintf("%.4g", diag(as.matrix(Matrix::expand(Matrix::lu(B))$U))),
                          collapse = " ")))
mod_b <- compile(normal_code)
for (sv in 1:3) {
  r <- run_once(mod_b, list(y = c(0.3, -0.2), s = s_b, K = K_b, solver = sv),
                c("lml", "lml_exact"))
  cat(sprintf("   solver %d: laplace_marginal_tol %.6f, exact %.6f\n", sv,
              r$vals[1], r$vals[2]))
}

# ── C. Saddle point ────────────────────────────────────────────────────────
saddle_code <- "
functions {
  real ll(vector theta, vector y, real s, real z, real sz) {
    return cauchy_lpdf(y | theta[1], s) + normal_lpdf(z | theta[2], sz);
  }
  matrix k_fun(matrix K) {
    return K;
  }
}
data {
  vector[2] y;
  real<lower=0> s;
  real z;
  real<lower=0> sz;
  matrix[2, 2] K;
  vector[2] theta0;
  int<lower=0, upper=1> marginal;
}
generated quantities {
  real lml = not_a_number();
  vector[2] mode = rep_vector(not_a_number(), 2);
  if (marginal) {
    lml = laplace_marginal_tol(ll, (y, s, z, sz), 2, k_fun, (K,),
                               (theta0, 1e-10, 100, 3, 100, 0));
  } else {
    tuple(vector[2], matrix[2, 2]) lap = laplace_latent_solve_tol(
      ll, (y, s, z, sz), 2, k_fun, (K,), (theta0, 1e-10, 100, 3, 100, 0));
    mode = lap.1;
  }
}
"
mod_c <- compile(saddle_code)
data_c <- list(y = c(-3, 3), s = 0.3, z = 0, sz = 0.1,
               K = matrix(c(4, 1, 1, 1), 2), theta0 = c(0.01, 3))
cat("\nC. Saddle point at (0, 0); start (0.01, 3), solver 3\n")
r <- run_once(mod_c, c(data_c, list(marginal = 1L)), "lml")
cat(sprintf("   laplace_marginal_tol: %s\n",
            if (length(r$exception)) paste("exception:", r$exception)
            else sprintf("returned %.6g", r$vals)))
r <- run_once(mod_c, c(data_c, list(marginal = 0L)), "mode")
cat(sprintf("   laplace_latent_solve_tol: %s\n",
            if (length(r$exception)) paste("exception:", r$exception)
            else sprintf("returned %s", paste(sprintf("%.4g", r$vals), collapse = " "))))

Tests

New tests:

• laplace_latent_solve.mean_is_final_iterate_tumor_size (test/unit/math/laplace/laplace_latent_solve_test.cpp), example A: with tolerance 1e-6, the returned mean must be within 0.01 posterior sd of a solve with tolerance 1e-14.
• laplace_marginal_tol.all_solvers_exact_with_negative_lu_pivot (new file test/unit/math/laplace/laplace_lu_log_determinant_test.cpp), example B: asserts that det(B) > 0 and that the LU factor of B has a negative diagonal entry; solvers 1, 2 and 3 must return the exact log marginal.
• laplace_marginal_tol.solver_3_saddle_point_throws (same new file), example C: from (0.01, 3), laplace_marginal_tol() and laplace_latent_solve_tol() must throw std::domain_error with the new message; from (0.3, 3), solver 3 must return the maximum (2.9385, 0.0097).

All three tests fail without the changes in stan/ (as in the develop column above). With the changes, they pass, and so do the other test files in test/unit/math/laplace/. I stopped the slow laplace_disease_map_test suite in laplace_marginal_lpdf_test.cpp after its first case, which passed. clang-format and cpplint pass on the changed files.

Side Effects

• laplace_latent_solve*() and laplace_latent*_rng(): the old and new means agree up to the convergence error. Results change only where the iteration stopped before exact convergence, and they move toward the mode.
• laplace_marginal*() with solver 3, or after a fallthrough to solver 3: where a negative pivot gave NaN, it now returns the log marginal. Where det(B) <= 0, it now throws std::domain_error. Before, it returned NaN, or a finite value when the permutation was odd and all pivots were positive. In the model block, NaN and an exception both reject the proposal; in generated quantities, the exception replaces a silent NaN or a wrong value.
• Solvers 1 and 2 are not affected by the log-determinant change: they use the Cholesky factor of a symmetric positive definite matrix.
• The Newton iteration can still converge to a saddle point (example C). This PR only makes the error name the cause.

Release notes

Corrected behavior and improved stability of laplace_latent_solve(), laplace_latent_rng(), laplace_marginal() and their variants.

Checklist

• Copyright holder: Aalto University

  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.06	5.94% faster
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.16	1.0	0.22% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.85	0.86	0.99	-0.98% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	1.01	0.5% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.1	0.11	0.99	-0.65% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	7.96	8.36	0.95	-5.02% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.7	3.71	1.0	-0.07% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.97	33.27	0.99	-0.9% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	0.98	-1.99% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.82	6.79	1.0	0.45% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	135.45	165.53	0.82	-22.2% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.22	3.23	1.0	-0.26% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.6	1.0	-0.01% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.18	25.83	0.97	-2.57% slower
performance.compilation	406.66	378.62	1.07	6.89% faster
Mean result: 0.9897691373393953

Jenkins Console Log
Jenkins Build Stages
Commit hash: 29e32c4215464fd5550540ce9a448016d22be56f

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

An update after #3469 (embedded Laplace returns the value at the initial guess) and stan-dev/stan#3442 (the Newton optimizer's Hessian is 1e-6 too small).

I ran the example of #3469 with develop (e4745c8) and with this PR (3c5e80e), -O0 -DNDEBUG. Develop prints the values in the issue. This PR prints the same values, except at μ = 1, α = 1e9, where develop gives +inf and this PR throws det(I + Sigma * W) is not positive at the final iterate. On a grid of the same example (μ ∈ {−5, 0, 1, 2, 5}, α = 10^7 to 10^12, 105 points), develop gives +inf at 87 points, this PR throws at the same 87 points, and the other 18 values are identical.

This PR does not fix part 1 of #3469: the μ = −30, α ≥ 1e5 rows still return the value at θ = 0. #3469 also reports a GP-model point where develop gives NaN and this PR gives a finite value, which is again the value at θ = 0. The log determinant is correct there (det(B) > 0 with LU pivots of both signs); the wrong value comes from part 1. Before this PR, the NaN made the samplers and optimizers reject such a point. So it would be good to fix part 1 soon after this PR.

stan::optimization::newton_step() treats only std::domain_error as a rejected point, and it accepts +inf and NaN, because f1 < f0 is false for both. With this PR, the points where det(I + Σ W) is not positive throw std::domain_error, so the Newton optimizer rejects them. stan-dev/stan#3442 makes it less likely that the optimizer goes there.

[jachymb]

Just a heads up, #3473 fixes a somewhat related issue and should probably me merged only after this one. (no conflict)

[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.83% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	1.05	4.99% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.33	6.35	1.0	-0.45% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.0	-0.27% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.3	8.29	1.0	0.12% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.45	4.46	1.0	-0.11% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	217.99	217.51	1.0	0.22% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.9	0.91	0.98	-1.94% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.61	0.61	1.0	-0.02% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	43.03	42.45	1.01	1.35% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.59	0.99	-0.52% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.95	21.08	0.99	-0.6% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	-0.32% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.18	3.2	0.99	-0.54% slower
performance.compilation	380.45	359.52	1.06	5.5% faster
Mean result: 1.0059462127355647

Jenkins Console Log
Jenkins Build Stages
Commit hash: f53903f0bee49582a667973d404c4e486f943759

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

[stevieb-bot]

Action: Not Ready for Human Review

Commit checked: 3c5e80ef4ce8

Summary

Thank you for this PR and for the detailed examples and tests.

The PR says that laplace_latent_solve*() and the laplace_latent*_rng() functions now use the final Newton iterate as the mean, and that the solver-3 log determinant now uses |U_ii| and throws when det(B) is not positive.
It gives three examples (A, B, and C), each with a new unit test, and an R script that produced the table.

I read the changed code at 3c5e80e and ran the PR's tests in a container without network access, at the PR head and at the merge base e4745c8 with the PR's test files copied in.
I ran the R script from the description at both commits with CmdStan develop (e4510f7726), building the models against the Stan Math of each commit.
I also built and ran the example program from #3469 at both commits, because your comment gives the output of that program as evidence.

The output of each run is in this gist.
Its README.md gives the command and commit behind each file.

• Claims and Evidence: 17 Confirmed, 6 Partially Confirmed, 2 Not Verified.
• Description vs Code: 10 Confirmed, 1 Partially Confirmed.

Claims and Evidence

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: The new tests fail on develop and pass with the PR

Quote:

All three tests fail without the changes in stan/ (as in the develop column above).
With the changes, they pass

(PR description, Tests)

STATUS: Confirmed

All three new tests fail at the merge base with the PR's test files copied in, and pass at the head.
At the base, mean_is_final_iterate_tumor_size reports a distance of 0.527581 posterior sd, which matches the 0.53 sd of the develop column.
all_solvers_exact_with_negative_lu_pivot fails only for solver 3, where lml is -nan and the exact value is -2.3767541633169298.
solver_3_saddle_point_throws fails because laplace_marginal_tol() does not throw.
The test returns at that FAIL(), so the base run does not reach the laplace_latent_solve_tol() part of the test.

Evidence:

The three new tests in the PR.
Code: laplace_latent_solve_test.cpp, laplace_lu_log_determinant_test.cpp.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/laplace/laplace_latent_solve_test.cpp test/unit/math/laplace/laplace_lu_log_determinant_test.cpp at the head (8 and 2 tests pass) and at the base with the two test files copied in (1 of 8 and 2 of 2 fail).
The exit status of both runs is 0 because of --run-all.
The gtest summary lines give the result.
Output in the gist: new_tests-head.out, new_tests-base.out.

---

FACT: The other Laplace test files pass

Quote:

With the changes, they pass, and so do the other test files in test/unit/math/laplace/.

(PR description, Tests)

STATUS: Confirmed

All 19 test files in test/unit/math/laplace/ build and run at the head, including the full laplace_disease_map_test suite.
The gtest summaries add up to 254 passed, 144 skipped, and 0 failed.
The test files skip these cases themselves, with the message Skipping test for hessian of size N with hessian block size of B, where N is not a multiple of B (2 or 3).
The run took 74 minutes.

Evidence:

The existing test files in test/unit/math/laplace/.
Code: test/unit/math/laplace/.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/laplace at the head.
Output in the gist: laplace_dir-head.out.

---

FACT: Solver 3 returned NaN on develop for a Gaussian likelihood

Quote:

The old sum(log(diag(U))) was then NaN, and laplace_marginal*() returned NaN without a message, with solver 3 or after a fallthrough to solver 3, also for a Gaussian likelihood (example B).

(PR description, Summary)

STATUS: Confirmed

At the base, laplace_marginal_tol() with solver 3 returns -nan for example B, a Gaussian likelihood.
The base code computes lu.matrixLU().diagonal().array().log().sum() and writes nothing to msgs (base L909).
After solvers 1 and 2 fail and allow_fallthrough is set, the same LUSolver runs (base L1282-L1286), so the fallthrough path computes the same log determinant.

Evidence:

Example B as the test all_solvers_exact_with_negative_lu_pivot.
Code: all_solvers_exact_with_negative_lu_pivot.

What I read: The base log determinant at L909, and the solver-3 branch that runs directly or after a fallthrough at base L1282-L1286.

What I ran: The base run of the test prints lml evaluates to -nan under the trace solver = 3.
Output in the gist: new_tests-base.out.

---

FACT: Example B, solver 3: NaN on develop, the exact value with the PR

Quote:

| B | laplace_marginal_tol(), solver 3 | NaN | -2.376754 (exact) |

(PR description, table)

STATUS: Confirmed

The test compares the result with multi_normal_lpdf(y | 0, K + s^2 I), which is the exact log marginal for a Gaussian likelihood.
At the base, solver 3 gives -nan and the exact value is -2.3767541633169298.
At the head, the test passes with a tolerance of 1e-8.
The test uses tolerance 1e-12 and 1000 steps, where the R script uses 1e-10 and 100.
The R script prints solver 3: laplace_marginal_tol NaN, exact -2.376754 at the base and solver 3: laplace_marginal_tol -2.376754, exact -2.376754 at the head.

Evidence:

Example B as the test all_solvers_exact_with_negative_lu_pivot, and section B of the R script in the PR description.
Code: R script in the PR description, all_solvers_exact_with_negative_lu_pivot.

What I ran: The same test runs as in the fact check "The new tests fail on develop and pass with the PR", and the R script, at the head and at the base.
Output in the gist: new_tests-head.out, new_tests-base.out, r_script-head.out, r_script-base.out.

---

FACT: Example B, solvers 1 and 2: the exact value before and after

Quote:

| B | laplace_marginal_tol(), solvers 1, 2 | -2.376754 | -2.376754 |

(PR description, table)

STATUS: Confirmed

The same test loops over solvers 1, 2, and 3.
At the base, the only failure is under the trace solver = 3, so solvers 1 and 2 match the exact value within 1e-8.
At the head, all three pass.
The R script prints -2.376754 for solvers 1 and 2 at both commits.

Evidence:

Example B as the test all_solvers_exact_with_negative_lu_pivot, and section B of the R script in the PR description.
Code: R script in the PR description, all_solvers_exact_with_negative_lu_pivot.

What I ran: The same test runs as in the fact check "The new tests fail on develop and pass with the PR", and the R script, at the head and at the base.
Output in the gist: new_tests-head.out, new_tests-base.out, r_script-head.out, r_script-base.out.

---

FACT: The latent solve and the RNGs throw the same message at a saddle point

Quote:

compute_log_determinant() also runs in laplace_latent_solve*() and laplace_latent*_rng(), so they throw the same message (example C).

(PR description, Summary)

STATUS: Confirmed

laplace_base_rng() calls laplace_marginal_density_est() (L57-L58), whose Newton loop ends with compute_log_determinant() (L1043, L1054).
At the head, the test solver_3_saddle_point_throws passes, so laplace_latent_solve_tol() throws std::domain_error with the text det(I + Sigma * W) is not positive from the start (0.01, 3).

Evidence:

Example C as the test solver_3_saddle_point_throws.
Code: solver_3_saddle_point_throws.

What I read: The call from laplace_base_rng() into the estimator at laplace_base_rng.hpp L57-L58, and the calls to compute_log_determinant() at the end of the Newton loop.

What I ran: The head run of the test passes.
Output in the gist: new_tests-head.out.

---

FACT: The #3469 example gives the same values except at mu = 1, alpha = 1e9

Quote:

I ran the example of #3469 with develop (e4745c8) and with this PR (3c5e80e), -O0 -DNDEBUG.
Develop prints the values in the issue.
This PR prints the same values, except at μ = 1, α = 1e9, where develop gives +inf and this PR throws det(I + Sigma * W) is not positive at the final iterate.

(comment by @avehtari)

STATUS: Confirmed

I built the example program from #3469 with the command in the issue (-O0 -DNDEBUG) at both commits.
The base prints the eight values from the issue, including inf at mu = 1, alpha = 1e9.
The head prints the same seven other values and throws laplace_marginal_density: det(I + Sigma * W) is not positive at the final iterate, ... at mu = 1, alpha = 1e9.

Evidence:

The example program in #3469, which the comment cites.
Code: #3469 example.

What I ran: g++ -std=c++17 -O0 -DNDEBUG ... issue3469_example.cpp -o example -L lib/tbb -ltbb -pthread && LD_LIBRARY_PATH=lib/tbb ./example at the head and the base.
The source is the program from the issue, copied as written.
Output in the gist: issue3469_example-head.out, issue3469_example-base.out.

---

FACT: newton_step() rejects only std::domain_error and accepts +inf and NaN

Quote:

stan::optimization::newton_step() treats only std::domain_error as a rejected point, and it accepts +inf and NaN, because f1 < f0 is false for both.

(comment by @avehtari)

STATUS: Confirmed

In stan::optimization::newton_step() on stan-dev/stan develop at 33ebe6fc9e, the step halving loop is while (f1 < f0), and the only catch is for std::domain_error, which sets f1 = -1e100 (L53-L67).
With f1 equal to +inf or NaN, f1 < f0 is false, and the loop accepts the new point.

Evidence:

The code of newton_step() in stan-dev/stan.
Code: newton.hpp.

What I read: newton.hpp L53-L67 on stan-dev/stan develop.

---

FACT: With the PR, the Newton optimizer rejects points where det(I + Sigma W) is not positive

Quote:

With this PR, the points where det(I + Σ W) is not positive throw std::domain_error, so the Newton optimizer rejects them.

(comment by @avehtari)

STATUS: Confirmed

The new check throws std::domain_error when det(B) is not positive (L913-L923), and newton_step() catches std::domain_error, sets f1 = -1e100, and halves the step (newton.hpp L53-L67).
In the #3469 example, the head throws this error at mu = 1, alpha = 1e9.

Evidence:

The new check in the PR, the code of newton_step() in stan-dev/stan, and the #3469 example.
Code: compute_log_determinant(), newton.hpp.

What I read: The new check at L913-L923 and the catch in newton_step().

What I ran: The #3469 example at the head.
Output in the gist: issue3469_example-head.out.

---

FACT: Example A: the posterior sd is 5.6 to 9.4 times smaller than the prior sd

Quote:

The posterior sd is 5.6 to 9.4 times smaller than the prior sd.

(PR description, Examples)

STATUS: Confirmed

At both commits, the R script prints posterior sd: 0.088 0.146 0.119 0.150 prior sd: 0.4925 0.9217 1.12 1.299.
The prior sd over the posterior sd ranges from 5.6 for the first latent variable to 9.4 for the third.

Evidence:

The R script in the PR description (section A).
Code: R script in the PR description.

What I ran: The R script at the head and the base, as in the fact check above.
Output in the gist: r_script-head.out, r_script-base.out.

---

FACT: Table, example A with tolerance 1e-6

Quote:

| A, tol 1e-6 | laplace_latent_solve_tol() | 0.53 sd, 7.12 lower | 0.0012 sd, 2.3e-6 lower |

(PR description, table)

STATUS: Confirmed

The R script prints 0.53 sd and 7.12 lower at the base, and 0.0012 sd and 2.28e-06 lower at the head.
The table gives 0.53 sd, 7.12 lower, 0.0012 sd, and 2.3e-6 lower.

Evidence:

The R script in the PR description (section A).
Code: R script in the PR description.

What I ran: The R script at the head and the base.
Output in the gist: r_script-head.out, r_script-base.out.

---

FACT: Table, example A with the default tolerance

Quote:

| A, default tol | laplace_latent_solve_tol() | 0.068 sd, 0.135 lower | 0.00015 sd, 3.8e-8 lower |

(PR description, table)

STATUS: Confirmed

The R script runs the default tolerance as 1.49012e-08.
It prints 0.068 sd and 0.135 lower at the base, and 0.00015 sd and 3.83e-08 lower at the head.
The table gives 0.068 sd, 0.135 lower, 0.00015 sd, and 3.8e-8 lower.

Evidence:

The R script in the PR description (section A).
Code: R script in the PR description.

What I ran: The R script at the head and the base.
Output in the gist: r_script-head.out, r_script-base.out.

---

FACT: Table, example C with laplace_marginal_tol(): NaN on develop, the new exception with the PR

Quote:

| C | laplace_marginal_tol() | NaN | new exception |

(PR description, table)

STATUS: Confirmed

From the start (0.01, 3) with solver 3, the R script prints laplace_marginal_tol: returned NaN at the base.
At the head it prints laplace_marginal_tol: exception: laplace_marginal_density: det(I + Sigma * W) is not positive at the final iterate, ....
The test solver_3_saddle_point_throws agrees.
It fails at the base because the call does not throw, and it passes at the head.

Evidence:

The R script in the PR description (section C) and the test solver_3_saddle_point_throws.
Code: R script in the PR description, solver_3_saddle_point_throws.

What I ran: The R script at the head and the base, and the test runs in the fact check "The new tests fail on develop and pass with the PR".
Output in the gist: r_script-head.out, r_script-base.out, new_tests-head.out, new_tests-base.out.

---

FACT: Table, example C with laplace_latent_solve_tol(): cholesky_decompose exception on develop

Quote:

| C | laplace_latent_solve_tol() | cholesky_decompose exception | new exception |

(PR description, table)

STATUS: Confirmed

At the base, the R script prints laplace_latent_solve_tol: exception: cholesky_decompose: Matrix m is not positive definite.
At the head it prints the new det(I + Sigma * W) is not positive at the final iterate, ... exception.
The test solver_3_saddle_point_throws passes at the head.
At the base it stops at its first FAIL(), so for the base only the R script shows the cholesky_decompose exception.

Evidence:

The R script in the PR description (section C) and the test solver_3_saddle_point_throws.
Code: R script in the PR description, solver_3_saddle_point_throws.

What I ran: The R script at the head and the base, and the head test run in the fact check "The new tests fail on develop and pass with the PR".
Output in the gist: r_script-head.out, r_script-base.out, new_tests-head.out.

---

FACT: Where the table values come from

Quote:

The develop column is from CmdStan 2.40.0; A and B give the same values with develop 43bc568.
The "this PR" column is from CmdStan develop with this branch.

(PR description, table notes)

STATUS: Confirmed

My base run uses CmdStan develop at e4510f7726, which reports version 2.40.0, with Stan Math at the merge base e4745c8.
It prints every value in the develop column for all three examples.
My head run uses the same CmdStan with this branch at 3c5e80e and prints every value in the "this PR" column.
The merge base contains 43bc568.
I did not build the CmdStan 2.40.0 release or Stan Math at 43bc568 itself (see Not Checked).

Evidence:

The R script in the PR description.
Code: R script in the PR description.

What I ran: The R script at the head and the base.
Output in the gist: r_script-head.out, r_script-base.out.

---

FACT: Checklist: the basic tests are passing

Quote:

the basic tests are passing

(PR description, Checklist)

STATUS: Confirmed

All 34 CI checks at the reviewed commit pass, including the GitHub Actions unit test groups and the Jenkins run, which reports "This commit looks good".
In the container, all 19 test files in test/unit/math/laplace/ pass at the head, including the full laplace_disease_map_test suite that the Tests section says was stopped after its first case.
I did not run the full unit suite, the header checks, the dependency checks, or the doxygen build myself.
The cpplint item is a formatting check and is not part of this fact check.

Evidence:

The repository's unit tests, and the CI checks on the PR at the reviewed commit.
Code: CI checks, test/unit/math/laplace/.

What I read: The CI checks listed by gh pr checks 3438 for the reviewed commit, all with status pass.

What I ran: python3 runTests.py -j8 --run-all test/unit/math/laplace at the head.
Output in the gist: laplace_dir-head.out.

---

FACT: Checklist: the new changes are tested

Quote:

the new changes are tested

(PR description, Checklist)

STATUS: Confirmed

Each of the two source changes has a new test that fails at the merge base and passes at the head.
mean_is_final_iterate_tumor_size covers the new mean, and all_solvers_exact_with_negative_lu_pivot and solver_3_saddle_point_throws cover the new log determinant and the new exception.

Evidence:

The three new tests in the PR.
Code: laplace_latent_solve_test.cpp, laplace_lu_log_determinant_test.cpp.

What I ran: The same test runs as in the fact check "The new tests fail on develop and pass with the PR", at the head and at the base.
Output in the gist: new_tests-head.out, new_tests-base.out.

---

Partially Confirmed Facts

FACT: Example C converges to the saddle point from (0.01, 3)

Quote:

From the start (0.01, 3), the Newton iteration of solver 3 converges to the saddle point.

(PR description, Examples)

STATUS: Partially Confirmed

At the head, solver 3 throws the det(I + Sigma * W) is not positive error from (0.01, 3), and from (0.3, 3) it returns (2.9385, 0.0097) within 1e-4.
This shows that det(B) is not positive at the final iterate from (0.01, 3).
Neither the test nor the R script prints the final iterate, so the code does not show that the iteration ends at (0, 0).

Evidence:

Example C as the test solver_3_saddle_point_throws.
The final iterate from (0.01, 3) is not printed or checked.
Code: solver_3_saddle_point_throws.

What I ran: The head run of the test passes.
Output in the gist: new_tests-head.out.

Request: If you have code that shows the iteration from (0.01, 3) ends at (0, 0), for example a run that prints the final iterate, could you link it?
Otherwise the sentence could say that det(B) is not positive at the final iterate, which is what the test shows.

---

FACT: Old and new means agree up to the convergence error and move toward the mode

Quote:

laplace_latent_solve*() and laplace_latent*_rng(): the old and new means agree up to the convergence error.
Results change only where the iteration stopped before exact convergence, and they move toward the mode.

(PR description, Side Effects)

STATUS: Partially Confirmed

In example A, the base mean is 0.53 posterior sd from a tolerance-1e-14 solve, and the head mean is within 0.01 sd.
The R script measures example A against a reference mode at three tolerances.
At each tolerance the distance falls, from 0.53 to 0.0012, from 0.068 to 0.00015, and from 1.4e-08 to 9e-11 posterior sd.
The existing test all_solvers_match_closed_form_normal compares the mean with the exact posterior mean within 1e-8 and passes at the base and at the head (L164-L171).
The two Cauchy tests compare the mean only with another solve at the same commit, so they do not show that the old and new means agree.
The sentence is general.
These runs cover one example where the mean moved toward the mode and one where the old and new means agree.

Evidence:

Example A as the test mean_is_final_iterate_tumor_size and in section A of the R script, and the existing tests in laplace_latent_solve_test.cpp.
Code: R script in the PR description, laplace_latent_solve_test.cpp.

What I ran: The same test runs as in the fact check "The new tests fail on develop and pass with the PR", and the R script, at the head and at the base.
Output in the gist: new_tests-head.out, new_tests-base.out, r_script-head.out, r_script-base.out.

Request: Could you say that this is shown for example A and the closed-form normal test, or link a broader comparison if you ran one?

---

FACT: Where det(B) <= 0, develop returned NaN or a finite value

Quote:

where a negative pivot gave NaN, it now returns the log marginal.
Where det(B) <= 0, it now throws std::domain_error.
Before, it returned NaN, or a finite value when the permutation was odd and all pivots were positive.

(PR description, Side Effects)

STATUS: Partially Confirmed

At the head, the saddle point of example C throws, and so does the #3469 example at mu = 1, alpha = 1e9.
The finite case follows from the base code, which sums the logs of the pivots without checking the sign of the permutation (base L909).
Your comment on #3469 already reports +inf at mu = 1, alpha = 1e9.
The base run also prints inf there, and the head throws the new error.
The Side Effects list of earlier results leaves out +inf.

Evidence:

Examples B and C as tests, and the example program of #3469 that your comment cites.
Code: laplace_lu_log_determinant_test.cpp, #3469 example.

What I read: The base log determinant at L909 and the new check at L913-L923.

What I ran: The new tests at the head and the base, and the #3469 example program, built with the command in the issue and run at the head and the base.
Output in the gist: new_tests-head.out, new_tests-base.out, issue3469_example-head.out, issue3469_example-base.out.

Request: Could you add +inf to the list of earlier results, as for a zero pivot in the #3469 example?

---

FACT: NaN and an exception both reject the proposal in the model block

Quote:

In the model block, NaN and an exception both reject the proposal; in generated quantities, the exception replaces a silent NaN or a wrong value.

(PR description, Side Effects)

STATUS: Partially Confirmed

The generated-quantities part is shown by section C of the R script, which calls laplace_marginal_tol() in generated quantities.
At the base it prints laplace_marginal_tol: returned NaN without an exception, and at the head the chain output has Chain 1 Exception: laplace_marginal_density: det(I + Sigma * W) is not positive ....
The model-block part follows from the stan-dev/stan code on develop.
A std::domain_error in the log density sets the potential to +inf (base_hamiltonian.hpp L52-L58), and NUTS sets a NaN Hamiltonian to +inf (base_nuts.hpp L258-L260), so both reject the proposal.
In generated quantities, a std::domain_error is logged and the values are written as NaN (mcmc_writer.hpp L108-L132).
The runs do not show the case where the exception replaces a wrong value.

Evidence:

Section C of the R script in the PR description, and the stan-dev/stan code that handles the NaN and the exception.
The runs do not cover the wrong-value case.
Code: R script in the PR description.

What I read: base_hamiltonian.hpp L52-L58, base_nuts.hpp L258-L260, and mcmc_writer.hpp L108-L132 on stan-dev/stan develop at 33ebe6fc9e.

What I ran: The R script at the head and the base.
Output in the gist: r_script-head.out, r_script-base.out.

Request: Could you link a run for the wrong-value case, where the permutation is odd and all pivots are positive?
There develop returns a finite value and this PR throws.

---

FACT: The mu = -30, alpha >= 1e5 rows still return the value at theta = 0

Quote:

This PR does not fix part 1 of #3469: the μ = −30, α ≥ 1e5 rows still return the value at θ = 0.

(comment by @avehtari)

STATUS: Partially Confirmed

In the #3469 example, the head prints the same values as the base for mu = -30, alpha = 1e5, 1e7, and 1e9 (-429.759639, -431.6906935, and -436.2852909).
The claim that these are the values at theta = 0 rests on the last column of the table in #3469, which the issue says was computed with mpmath.
That code is not posted.

Evidence:

The example program in #3469.
The code for the reference values at theta = 0 is not posted.
Code: #3469 example.

What I ran: The #3469 example at the head and the base.
Output in the gist: issue3469_example-head.out, issue3469_example-base.out.

Request: Could you link the code for the value at theta = 0, or say that the run shows these rows are unchanged from develop?

---

FACT: Before this PR, the NaN made the samplers and optimizers reject such a point

Quote:

Before this PR, the NaN made the samplers and optimizers reject such a point.

(comment by @avehtari)

STATUS: Partially Confirmed

The stan-dev/stan code on develop agrees for NUTS and BFGS but not for the Newton optimizer.
NUTS sets a NaN Hamiltonian to +inf and so rejects the point (base_nuts.hpp L258-L260).
The BFGS function evaluation returns error code 2 for a non-finite value, which counts as a failed evaluation (bfgs.hpp L318-L326).
newton_step() accepts a NaN, because f1 < f0 is false (newton.hpp L53-L67), as the next paragraph of your comment also says.

Evidence:

The stan-dev/stan code for NUTS, BFGS, and the Newton optimizer.
Code: base_nuts.hpp, bfgs.hpp, newton.hpp.

What I read: The three permalinks above, on stan-dev/stan develop at 33ebe6fc9e.
The first two show a rejection of a NaN, and the third shows that newton_step() accepts it.

Request: Could you say which samplers and optimizers reject a NaN, since the next paragraph says the Newton optimizer accepts it?

---

Not Verified Facts

FACT: A grid of 105 points: +inf on develop and an exception with the PR at the same 87 points

Quote:

On a grid of the same example (μ ∈ {−5, 0, 1, 2, 5}, α = 10^7 to 10^12, 105 points), develop gives +inf at 87 points, this PR throws at the same 87 points, and the other 18 values are identical.

(comment by @avehtari)

STATUS: Not Verified

The comment does not include or link the code for the grid.
The program in #3469 covers 8 points, which the fact check "The #3469 example gives the same values except at mu = 1, alpha = 1e9" runs.

Evidence:

There is no code for the grid.

Request: Could you post the program for the grid of 105 points?

---

FACT: A GP-model point gives NaN on develop and a finite value at theta = 0 with the PR

Quote:

#3469 also reports a GP-model point where develop gives NaN and this PR gives a finite value, which is again the value at θ = 0.
The log determinant is correct there (det(B) > 0 with LU pivots of both signs); the wrong value comes from part 1.

(comment by @avehtari)

STATUS: Not Verified

#3469 gives the point (rho = 1, mu = 1, alpha = 1e8) and the value -866.131.
Neither the comment nor the issue includes the GP model or the data.

Evidence:

There is no code for the GP model or the data.

Request: Could you post the GP model and the data?

---

Description vs Code

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: The latent mean is now the final Newton iterate

Quote:

laplace_base_rng() computed the mean as covariance_train * md_est.theta_grad, that is K * grad log p(y | theta_hat).
[...]
The mean is now md_est.theta.

(PR description, Summary)

STATUS: Confirmed

At the base, laplace_base_rng() sets mean_train = covariance_train * md_est.theta_grad (base L59).
At the head, it sets mean_train = md_est.theta (L59).

---

FACT: The stopping rule bounds the change of the objective

Quote:

The stopping rule bounds the change of the objective, not this difference

(PR description, Summary)

STATUS: Confirmed

The Newton loop stops when std::abs(state.curr().obj() - state.prev().obj()) < options.tolerance (L1030-L1034) or when the line search fails.

---

FACT: The change affects the latent solve and the latent RNGs

Quote:

This affects laplace_latent_solve[_tol]() and the mean of the draws of laplace_latent[_tol]_rng() and laplace_latent[_tol]_{bernoulli_logit,poisson_log,neg_binomial_2_log}_rng().

(PR description, Summary)

STATUS: Confirmed

laplace_latent_solve and laplace_latent_solve_tol call laplace_base_rng<true>() (L47, L81).
The eight RNGs call laplace_base_rng() in laplace_latent_rng.hpp (L44, L78), laplace_latent_bernoulli_logit_rng.hpp (L42, L76), laplace_latent_poisson_log_rng.hpp (L44, L79), and laplace_latent_neg_binomial_2_log_rng.hpp (L50, L92).
laplace_base_rng() passes mean_train to multi_normal_rng() or returns it (L80-L82).

---

FACT: The solver-3 log determinant uses |U_ii| and checks the sign

Quote:

The log determinant is now sum(log(abs(diag(U)))), and the sign of det(B) is the permutation sign times the signs of the U_ii.

(PR description, Summary)

STATUS: Confirmed

LUSolver::compute_log_determinant() counts the negative entries of the LU diagonal, compares their parity with the sign of the permutation, and returns u_diag.array().abs().log().sum() (L913-L923).

---

FACT: A non-positive det(B) throws std::domain_error that names the cause

Quote:

The function then throws std::domain_error with a message that names this cause.

(PR description, Summary)

STATUS: Confirmed

When the sign is negative or a diagonal entry is zero, the function throws std::domain_error with the message laplace_marginal_density: det(I + Sigma * W) is not positive at the final iterate, so the negative Hessian of the log density is not positive definite there (for example at a saddle point). (L917-L921).

---

FACT: Solvers 1 and 2 are not affected by the log-determinant change

Quote:

Solvers 1 and 2 are not affected by the log-determinant change: they use the Cholesky factor of a symmetric positive definite matrix.

(PR description, Side Effects)

STATUS: Confirmed

In laplace_marginal_density_estimator.hpp, the diff changes only LUSolver::compute_log_determinant().
The three Cholesky solvers return 2.0 * llt_B.matrixLLT().diagonal().array().log().sum() (L573-L574, L703-L704, L814-L815), as at the base.

---

FACT: The test for example A

Quote:

laplace_latent_solve.mean_is_final_iterate_tumor_size (test/unit/math/laplace/laplace_latent_solve_test.cpp), example A: with tolerance 1e-6, the returned mean must be within 0.01 posterior sd of a solve with tolerance 1e-14.

(PR description, Tests)

STATUS: Confirmed

mean_is_final_iterate_tumor_size solves example A with tolerance 1e-14 and 1e-6, takes the posterior sd from the Cholesky factor of the 1e-14 solve, and checks that the largest distance in posterior sd is below 0.01 (L303-L311).

---

FACT: The test for example B

Quote:

laplace_marginal_tol.all_solvers_exact_with_negative_lu_pivot (new file test/unit/math/laplace/laplace_lu_log_determinant_test.cpp), example B: asserts that det(B) > 0 and that the LU factor of B has a negative diagonal entry; solvers 1, 2 and 3 must return the exact log marginal.

(PR description, Tests)

STATUS: Confirmed

all_solvers_exact_with_negative_lu_pivot asserts that the LU diagonal has a negative entry and that det(B) > 0, then checks solvers 1, 2, and 3 against multi_normal_lpdf(y | 0, K + s^2 I) within 1e-8 (L44-L67).

---

FACT: The test for example C

Quote:

laplace_marginal_tol.solver_3_saddle_point_throws (same new file), example C: from (0.01, 3), laplace_marginal_tol() and laplace_latent_solve_tol() must throw std::domain_error with the new message; from (0.3, 3), solver 3 must return the maximum (2.9385, 0.0097).

(PR description, Tests)

STATUS: Confirmed

solver_3_saddle_point_throws checks that laplace_marginal_tol() and laplace_latent_solve_tol() throw std::domain_error containing det(I + Sigma * W) is not positive from (0.01, 3), and that solver 3 from (0.3, 3) returns (2.9385, 0.0097) within 1e-4 (L75-L113).

---

FACT: Checklist: the changes are documented in the doxygen

Quote:

the code is written in idiomatic C++ and changes are documented in the doxygen

(PR description, Checklist)

STATUS: Confirmed

The doxygen of LUSolver::compute_log_determinant() now describes the negative pivots, the sign of det(B), the return value sum(log(abs(diag(U)))), and the std::domain_error (L899-L912).
The doxygen of laplace_base_rng() already says that it returns the posterior mean of the Laplace approximation, and the one-line change keeps that true (L13-L42).
Whether the code is idiomatic C++ is a judgment for the code review, which this fact check does not do.

---

Partially Confirmed Facts

FACT: The PR only makes the error at a saddle point name the cause

Quote:

The Newton iteration can still converge to a saddle point (example C).
This PR only makes the error name the cause.

(PR description, Side Effects)

STATUS: Partially Confirmed

The source diff changes two places, the mean in laplace_base_rng() (L59) and LUSolver::compute_log_determinant() (L913-L923).
run_newton_loop() and the solver steps are the same as at the base, so the iteration can still end at a saddle point.
For laplace_latent_solve_tol() at the saddle point of example C, the error changes from cholesky_decompose to the new message.
For laplace_marginal_tol(), the R script shows that the base returns NaN without an error and the head throws, so the PR adds an error there.

Request: Could the sentence say that laplace_marginal_tol() returned NaN on develop at the saddle point and now throws, while laplace_latent_solve_tol() threw before and now names the cause?

---

Not Checked

• The CmdStan 2.40.0 release and Stan Math at 43bc568, which the description names as the builds behind the develop column.
  The runs here build the R script only against the PR head and the merge base, with CmdStan develop at e4510f7726.
• The copyright line of the checklist, which is a statement about licensing.