Embedded Laplace: run the requested solver when allow_fallthrough is on

stan/math/mix/functor/laplace_marginal_density_estimator.hpp

@@ -57,7 +57,11 @@ struct laplace_options_base {
    */
   double tolerance{internal::laplace_default_tolerance};  // 2
   /* Maximum number of steps*/
-  int max_num_steps{internal::laplace_default_max_num_steps};          // 3
+  int max_num_steps{internal::laplace_default_max_num_steps};  // 3
+  /**
+   * If true, a failure of the requested solver falls through to the next
+   * one (1 to 2 to 3). The requested solver always runs first.
+   */
   int allow_fallthrough{internal::laplace_default_allow_fallthrough};  // 4
   laplace_line_search_options line_search{
       internal::laplace_default_max_steps_line_search};  // 5
@@ -1260,7 +1264,8 @@ inline auto laplace_marginal_density_est(
         "solver 2 (Covariance-root Cholesky)", e);
   }
   try {
-    if (options.solver == 2 || options.allow_fallthrough) {
+    if (options.solver == 2
+        || (options.allow_fallthrough && options.solver == 1)) {
       CholeskyKSolver solver(state, covariance);
       return run_newton_loop(solver, state, options, step_iter, ll_fun, ll_args,
                              covariance, update_fun, msgs);

test/unit/math/laplace/laplace_solver_fallback_test.cpp

@@ -34,6 +34,28 @@ double run(bool allow_fallthrough, std::ostream* msgs) {
       msgs);
 }
 
+// Strictly concave everywhere, so each of the three solvers converges on
+// its own and none of them has to fall back to the next one.
+struct ConcaveLikelihood {
+  template <typename Theta>
+  auto operator()(const Theta& theta, std::ostream* /*msgs*/) const {
+    return -0.5 * stan::math::square(theta(0) - 1.0)
+           - 0.5 * stan::math::square(theta(1) + 1.0);
+  }
+};
+
+// The solver that produced the estimate for the given likelihood.
+template <typename LL>
+int used_solver(const LL& ll, int solver, bool allow_fallthrough) {
+  Eigen::VectorXd theta0 = Eigen::VectorXd::Zero(2);
+  stan::math::laplace_options_user_supplied options(
+      theta0, 1e-10, 100, 1, solver, 100, allow_fallthrough);
+  return stan::math::internal::laplace_marginal_density_est(
+             ll, std::tuple<>{},
+             Eigen::MatrixXd(Eigen::MatrixXd::Identity(2, 2)), options, nullptr)
+      .solver_used;
+}
+
 }  // namespace
 
 TEST(LaplaceSolverFallback, FallbackIsLoggedToStream) {
@@ -72,3 +94,16 @@ TEST(LaplaceSolverFallback, DisabledFallthroughThrowsWithReason) {
         << e.what();
   }
 }
+
+// The requested solver runs whether or not the fallthrough is allowed; the
+// fallthrough only decides what happens after that solver has failed.
+TEST(LaplaceSolverFallback, RequestedSolverRunsWithFallthroughAllowed) {
+  for (int solver = 1; solver <= 3; ++solver) {
+    EXPECT_EQ(solver, used_solver(ConcaveLikelihood{}, solver, true))
+        << "solver " << solver;
+    EXPECT_EQ(solver, used_solver(ConcaveLikelihood{}, solver, false))
+        << "solver " << solver;
+  }
+  // Guard: after solver 1 fails, the fallthrough lands on solver 2.
+  EXPECT_EQ(2, used_solver(PartlyConvexLikelihood{}, 1, true));
+}