Embedded Laplace: when the first Newton step is rejected, the log density at the initial guess is returned as the log marginal (wrong by a lot) #3469

I am working on fixing this, mostly just taking notes and documenting.

Description

internal::laplace_marginal_density_est (stan/math/mix/functor/laplace_marginal_density_estimator.hpp) can return a value that is not a Laplace approximation, without an error. In part 1 below, nothing is written to msgs. In part 2, the only output is the generic solver-fallback warnings (each printed once per process), which do not say that the returned value is wrong.

1. The initial guess is returned as the mode. run_newton_loop first evaluates the full Newton step. retry_evaluate halves the step while the objective or its directional derivative is not finite (or the evaluation throws), as long as the halved step stays above min_alpha (laplace_line_search_options: tau = 0.5, min_alpha = 1e-8). When every trial down to 2⁻²⁶ ≈ 1.5e-8 of the Newton step fails, the proposal is invalid and wolfe_line_search is not called. wolfe_status is then StepTooSmall and search_failed is set. The loop returns solver.build_result(state, ...) for state.prev(), the last accepted iterate. On the first iteration that is the initial guess θ = 0 (a = 0), so the result is the objective at θ = 0 minus half the log determinant there.

   This happens when, at θ = 0, the likelihood curvature W and the prior precision are both so small that the full Newton step is many orders of magnitude too long. In the example below (μ = −30, α = 1e5), W = e^μ ≈ 9.4e-14, and the prior precision along the vector of ones is 1/(5α²) = 2e-11. The Newton step from θ = 0, (K⁻¹ + W)⁻¹ ∇ log p(y | θ), therefore sets every θᵢ to (ȳ − e^μ) / (W + 1/(5α²)) ≈ 2.8 / 2.0e-11 = 1.39e11 (80-digit arithmetic gives the same). The mode is at 31.03, and exp(μ + θ) overflows for every trial step down to the last one, 2⁻²⁶ of the Newton step (θ ≈ 2.1e3).

2. +inf (or NaN) when solvers 1 and 2 fail and solver 3 factors a singular I + K W. When K is numerically singular and the entries of W^½ K W^½ at the current iterate exceed about 1/ε, the identity in B = I + W^½ K W^½ is lost and B is singular in floating point. In the example, K = α² J + 1e-8 I, whose jitter is lost from α = 1e5 on, so K is rank one there. At μ = 1, α = 1e9 (W ≈ e at θ = 0, entries of W^½ K W^½ about 2.7e18):

   • solver 1 fails (Cholesky failed after adding jitter up to 0.000100);
   • solver 2 fails (Cholesky of covariance failed at start);
   • with the default allow_fallthrough, solver 3 then factors the singular I + K W with LU.

   Solver 3's Newton proposal is not finite, so the loop again returns state.prev(). The three solvers share one Newton state. When the solvers fail already at θ = 0, as at the GP-model points in the notes (where θ = a = 0 is returned), that is the initial guess. LUSolver::compute_log_determinant then sums the logs of the LU diagonal:

   • a zero pivot gives log det = −∞, so the log marginal is obj - 0.5 * log_det = +inf;
   • a negative pivot gives NaN.

   Mathematically det(B) ≥ 1, so log det ≥ 0. Before B becomes singular, the result already loses accuracy, probably from rounding in B. It is 0.013 too high at μ = 1, α = 1e7 below, where the entries of W^½ K W^½ at the mode are 2.8e14.

Stan's Newton optimizer accepts a log density of +inf as an improvement (its step halving loops while (f1 < f0)); BFGS and L-BFGS treat a non-finite value as a failed evaluation. This was found with the Newton optimizer on a Poisson GP model using laplace_marginal_poisson_log. The optimizer overshot because of a separate bug in the scale of its finite-difference Hessian in stan (fix not yet submitted). It reached α = 4.7e110, ρ = 3.1e27, where the log marginal is +inf, and reported Log joint probability = inf in all its remaining iterations.

Both parts date from the Newton/Wolfe loop of #3250 (first in v5.3.0).

Example

Five Poisson counts y = (3, 1, 4, 1, 5) with log rates μ + θᵢ and θ ~ N(0, α² J + 1e-8 I), where J is the all-ones matrix. This is gp_exp_quad_cov with a length scale much longer than the span of the inputs, plus a 1e-8 jitter added to the diagonal with add_diag:

#include <stan/math/mix.hpp>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>

struct constant_cov {
  template <typename T>
  Eigen::Matrix<T, -1, -1> operator()(const T& alpha, int n,
                                      std::ostream*) const {
    return stan::math::add_diag(
        Eigen::Matrix<T, -1, -1>::Constant(n, n, alpha * alpha), 1e-8);
  }
};

int main() {
  const std::vector<int> y{3, 1, 4, 1, 5}, idx{1, 2, 3, 4, 5};
  std::cout << std::setprecision(10);
  for (double mu : {-30.0, 1.0}) {
    for (double alpha : {1e3, 1e5, 1e7, 1e9}) {
      const Eigen::VectorXd mean = Eigen::VectorXd::Constant(5, mu);
      std::stringstream msgs;
      std::cout << "mu " << mu << " alpha " << alpha << ": ";
      try {
        std::cout << stan::math::laplace_marginal_poisson_log_lpmf(
            y, idx, mean, 1, constant_cov{}, std::forward_as_tuple(alpha, 5),
            &msgs);
      } catch (const std::exception& e) {
        std::cout << "throws " << e.what();
      }
      std::cout << "\n" << msgs.str();
    }
  }
}

The program was built from the repository root, after make had built the bundled TBB into lib/tbb (any unit-test build does this), and run with LD_LIBRARY_PATH=lib/tbb ./example (WSL2, g++ 12.2):

g++ -std=c++17 -O0 -DNDEBUG -I . -I lib/eigen_5.0.1 -I lib/boost_1.87.0 -I lib/sundials_6.1.1/include -I lib/tbb_2020.3/include example.cpp -o example -L lib/tbb -ltbb -pthread

It prints:

mu -30 alpha 1000: -17.57039523
mu -30 alpha 100000: -429.759639
mu -30 alpha 10000000: -431.6906935
mu -30 alpha 1000000000: -436.2852909
mu 1 alpha 1000: -17.56991717
mu 1 alpha 100000: -22.17508644
mu 1 alpha 10000000: -26.76754568
mu 1 alpha 1000000000: inf
[laplace_marginal_density] WARNING: solver fallback
  iteration:  0
  failed:     solver 1 (Diagonal Hessian-root Cholesky)
  reason:     laplace_marginal_density: Cholesky failed after adding jitter up to 0.000100
  action:     trying solver 2 (Covariance-root Cholesky)
note: this warning message will only be displayed once.
[laplace_marginal_density] WARNING: solver fallback
  iteration:  0
  failed:     solver 2 (Covariance-root Cholesky)
  reason:     laplace_marginal_density: Cholesky of covariance failed at start
  action:     trying solver 3 (General LU solver)
note: this warning message will only be displayed once.

The μ = −30 rows run first, while the warnings have not yet been printed, and write nothing to msgs.

Expected Output

The expected result is the Laplace approximation at the mode, or else an error. The reference values below were computed with 80-digit arithmetic (mpmath: Newton for the mode, then ll(θ*) − ½ θᵀ K⁻¹ θ − ½ log det(I + W^½ K W^½)):

At μ = −30 and α ≥ 1e5, develop returns the value at the initial guess (last column) to all printed digits. The two other differences, 2.3e-6 at μ = −30, α = 1e3 and 8.8e-7 at μ = 1, α = 1e5, are much smaller and not part of this issue.

When the mode cannot be found or the log determinant cannot be computed, a std::domain_error would let Stan's samplers and optimizers reject the point. A value at a point that is not the mode is not an approximation of the log marginal, and neither is +inf or NaN.

Open PR #3438 makes LUSolver::compute_log_determinant throw when det(B) is not positive. With its two changed headers, the example throws at μ = 1, α = 1e9 (det(I + Sigma * W) is not positive at the final iterate), and every other line is unchanged; the μ = −30 rows still return the value at θ = 0. On the GP model where this was found, the NaN case (ρ = 1, μ = 1, α = 1e8; LU pivots of both signs) turns into a finite −866.131, which is again the value at θ = 0 (log likelihood −330.91 there). So #3438 removes the +inf, but not part 1.

Current Version:

v5.4.0 (develop at e4745c8)