Return the final Newton iterate as the Laplace latent mean, and fix the solver-3 log determinant

stan/math/mix/functor/laplace_base_rng.hpp

@@ -56,7 +56,7 @@ inline auto laplace_base_rng(LLFunc&& ll_fun, LLArgs&& ll_args,
       std::forward<CovarArgs>(covar_args));
   auto md_est = internal::laplace_marginal_density_est(
       ll_fun, std::forward<LLArgs>(ll_args), covariance_train, options, msgs);
-  Eigen::VectorXd mean_train = covariance_train * md_est.theta_grad;
+  Eigen::VectorXd mean_train = md_est.theta;
   Eigen::MatrixXd Sigma;
   if (md_est.solver_used == 1) {
     Eigen::MatrixXd V_dec

stan/math/mix/functor/laplace_marginal_density_estimator.hpp

@@ -900,13 +900,27 @@ struct LUSolver {
    * @brief Compute log determinant from LU factorization.
    *
    * @note This uses the diagonal of the combined LU matrix produced by Eigen
-   * (equivalently the diagonal of U). It does not account for the sign of the
-   * permutation; callers assume `det(B) > 0` in the Laplace correction.
+   * (equivalently the diagonal of U). B = I + Sigma * W is not symmetric, so
+   * with partial pivoting the diagonal of U can have negative entries also
+   * when `det(B) > 0`. The sign of `det(B)` is the sign of the permutation
+   * times the signs of the diagonal entries of U.
    *
-   * @return Sum of log of the LU diagonal entries.
+   * @return log(det(B)) = sum(log(abs(diag(U))))
+   * @throws std::domain_error if `det(B)` is not positive. Then
+   * `Sigma^-1 + W` is not positive definite at the final iterate, for example
+   * at a saddle point, and the Laplace approximation is not defined.
    */
   double compute_log_determinant() const {
-    return lu.matrixLU().diagonal().array().log().sum();
+    const auto& u_diag = lu.matrixLU().diagonal();
+    const bool odd_sign = (lu.permutationP().determinant() < 0)
+                          != ((u_diag.array() < 0.0).count() % 2 == 1);
+    if (odd_sign || (u_diag.array() == 0.0).any()) {
+      throw std::domain_error(
+          "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).");
+    }
+    return u_diag.array().abs().log().sum();
   }
 
   /**

test/unit/math/laplace/laplace_latent_solve_test.cpp

@@ -242,3 +242,70 @@ TEST(laplace_latent_tol_rng, solver_2_draw_uses_laplace_covariance) {
   Eigen::VectorXd draw_ref = stan::math::multi_normal_rng(mean, Sigma, rng_ref);
   EXPECT_MATRIX_NEAR(draw_ref, draw, 1e-8);
 }
+
+namespace {
+// Longitudinal part of a nonlinear joint model for tumor size:
+// 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).
+// The log likelihood is not concave in theta.
+struct tumor_size_likelihood {
+  template <typename Theta>
+  auto operator()(const Theta& theta, const Eigen::VectorXd& y,
+                  const Eigen::VectorXd& t, const Eigen::VectorXd& mu,
+                  double sigma, std::ostream* pstream) const {
+    using stan::math::exp;
+    using T = stan::value_type_t<Theta>;
+    const T bsld = mu(0) * exp(theta(0));
+    const T d = mu(1) * exp(theta(1));
+    const T g = mu(2) * exp(theta(2));
+    const T phi = stan::math::inv_logit(stan::math::logit(mu(3)) + theta(3));
+    Eigen::Matrix<T, Eigen::Dynamic, 1> m(y.size());
+    for (Eigen::Index i = 0; i < y.size(); ++i) {
+      m(i) = bsld * (phi * exp(-d * t(i)) + (1 - phi) * exp(g * t(i)));
+    }
+    return stan::math::normal_lpdf(y, m, stan::math::multiply(sigma, m));
+  }
+};
+
+struct diagonal_covariance {
+  Eigen::MatrixXd operator()(const Eigen::VectorXd& sd,
+                             std::ostream* pstream) const {
+    return sd.array().square().matrix().asDiagonal();
+  }
+};
+
+auto tumor_size_solve(double tolerance) {
+  Eigen::VectorXd y(12);
+  y << 201.04, 78.30, 48.01, 44.39, 52.91, 50.49, 48.63, 67.29, 52.21, 72.43,
+      102.43, 105.06;
+  Eigen::VectorXd t = Eigen::VectorXd::LinSpaced(12, 0.0, 693.0);
+  Eigen::VectorXd mu(4);
+  mu << 117.4, 9.121e-3, 1.913e-3, 0.3781;
+  Eigen::VectorXd sd(4);
+  sd << 0.4925, 0.9217, 1.120, 1.299;
+  const double sigma = 0.09826;
+  Eigen::VectorXd theta_0(4);
+  theta_0 << 0.3636, -0.3385, 1.9794, -1.9322;
+  return stan::math::laplace_latent_solve_tol(
+      tumor_size_likelihood{}, std::forward_as_tuple(y, t, mu, sigma), 4,
+      diagonal_covariance{}, std::make_tuple(sd),
+      std::make_tuple(theta_0, tolerance, 100, 3, 100, false), nullptr);
+}
+}  // namespace
+
+// The returned mean is the final Newton iterate.  It used to be
+// K * grad log p(y | theta), which differs from the iterate by
+// K * grad log p(theta | y) before exact convergence.  Here the posterior is
+// much narrower than the prior, and with tolerance 1e-6 the old value was
+// 0.53 posterior sd from the mode, while the final iterate is 0.0012 sd from
+// it.
+TEST(laplace_latent_solve, mean_is_final_iterate_tumor_size) {
+  auto [mean_ref, chol_ref] = tumor_size_solve(1e-14);
+  const Eigen::VectorXd post_sd
+      = (chol_ref * chol_ref.transpose()).diagonal().array().sqrt();
+  auto [mean_est, chol_est] = tumor_size_solve(1e-6);
+  const double max_dist
+      = ((mean_est - mean_ref).array().abs() / post_sd.array()).maxCoeff();
+  EXPECT_LT(max_dist, 0.01);
+}

test/unit/math/laplace/laplace_lu_log_determinant_test.cpp

@@ -0,0 +1,113 @@
+#include <stan/math.hpp>
+#include <stan/math/mix.hpp>
+#include <test/unit/util.hpp>
+
+#include <gtest/gtest.h>
+#include <cmath>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+
+namespace {
+struct normal_likelihood {
+  template <typename Theta>
+  auto operator()(const Theta& theta, const Eigen::VectorXd& y, double s,
+                  std::ostream* pstream) const {
+    return stan::math::normal_lpdf(y, theta, s);
+  }
+};
+
+// theta(0) has two Cauchy observations y, which make log p(theta | y)
+// bimodal in theta(0); theta(1) has one normal observation z.
+struct cauchy_normal_likelihood {
+  template <typename Theta>
+  auto operator()(const Theta& theta, const Eigen::VectorXd& y, double s,
+                  double z, double sz, std::ostream* pstream) const {
+    return stan::math::cauchy_lpdf(y, theta(0), s)
+           + stan::math::normal_lpdf(z, theta(1), sz);
+  }
+};
+
+struct fixed_covariance {
+  Eigen::MatrixXd operator()(const Eigen::MatrixXd& K,
+                             std::ostream* pstream) const {
+    return K;
+  }
+};
+}  // namespace
+
+// Solver 3 factorizes B = I + K W with partial pivoting.  Here the rows are
+// swapped and U has a negative diagonal entry although det(B) = 226 > 0.  The
+// log determinant used to be sum(log(diag(U))), which is NaN.  For a normal
+// likelihood the Laplace approximation is exact, so all solvers must return
+// log N(y | 0, K + s^2 I).
+TEST(laplace_marginal_tol, all_solvers_exact_with_negative_lu_pivot) {
+  Eigen::MatrixXd K(2, 2);
+  K << 1.0, 1.5, 1.5, 4.0;
+  Eigen::VectorXd y(2);
+  y << 0.3, -0.2;
+  const double s = std::sqrt(0.1);
+  const Eigen::MatrixXd B
+      = Eigen::MatrixXd::Identity(2, 2) + K / (s * s);  // W = I / s^2
+  Eigen::PartialPivLU<Eigen::MatrixXd> lu(B);
+  ASSERT_LT(lu.matrixLU().diagonal().minCoeff(), 0.0);
+  ASSERT_GT(B.determinant(), 0.0);
+  const double exact = stan::math::multi_normal_lpdf(
+      y, Eigen::VectorXd::Zero(2),
+      (K + s * s * Eigen::MatrixXd::Identity(2, 2)).eval());
+  Eigen::VectorXd theta_0 = Eigen::VectorXd::Zero(2);
+  for (int solver : {1, 2, 3}) {
+    SCOPED_TRACE("solver = " + std::to_string(solver));
+    const double lml = stan::math::laplace_marginal_tol<false>(
+        normal_likelihood{}, std::forward_as_tuple(y, s), 2, fixed_covariance{},
+        std::make_tuple(K),
+        std::make_tuple(theta_0, 1e-12, 1000, solver, 1000, false), nullptr);
+    EXPECT_NEAR(exact, lml, 1e-8);
+  }
+}
+
+// Started near the saddle point (0, 0) of log p(theta | y), the Newton
+// iteration of solver 3 converges to it.  There Sigma^-1 + W is not positive
+// definite and the Laplace approximation is not defined.  This must be
+// reported as such: laplace_marginal_tol() used to return NaN or a finite
+// value, and laplace_latent_solve_tol() threw from cholesky_decompose().
+// Started away from the saddle point, solver 3 finds a maximum.
+TEST(laplace_marginal_tol, solver_3_saddle_point_throws) {
+  Eigen::MatrixXd K(2, 2);
+  K << 4.0, 1.0, 1.0, 1.0;
+  Eigen::VectorXd y(2);
+  y << -3.0, 3.0;
+  const double s = 0.3;
+  const double z = 0.0;
+  const double sz = 0.1;
+  const std::string reason = "det(I + Sigma * W) is not positive";
+  Eigen::VectorXd theta_near(2);
+  theta_near << 0.01, 3.0;
+  auto ops_near = std::make_tuple(theta_near, 1e-10, 100, 3, 100, false);
+  try {
+    stan::math::laplace_marginal_tol<false>(
+        cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2,
+        fixed_covariance{}, std::make_tuple(K), ops_near, nullptr);
+    FAIL() << "laplace_marginal_tol: expected std::domain_error";
+  } catch (const std::domain_error& e) {
+    EXPECT_NE(std::string(e.what()).find(reason), std::string::npos)
+        << e.what();
+  }
+  try {
+    stan::math::laplace_latent_solve_tol(
+        cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2,
+        fixed_covariance{}, std::make_tuple(K), ops_near, nullptr);
+    FAIL() << "laplace_latent_solve_tol: expected std::domain_error";
+  } catch (const std::domain_error& e) {
+    EXPECT_NE(std::string(e.what()).find(reason), std::string::npos)
+        << e.what();
+  }
+  Eigen::VectorXd theta_away(2);
+  theta_away << 0.3, 3.0;
+  auto [mean, chol] = stan::math::laplace_latent_solve_tol(
+      cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2,
+      fixed_covariance{}, std::make_tuple(K),
+      std::make_tuple(theta_away, 1e-10, 100, 3, 100, false), nullptr);
+  EXPECT_NEAR(2.9385, mean(0), 1e-4);
+  EXPECT_NEAR(0.0097, mean(1), 1e-4);
+}