Merge pull request #3286 from stan-dev/fix/laplace-args

update laplace api to have hessian_block_size as a non-optional argument

Author: SteveBronder

doxygen/doxygen.cfg

@@ -2747,7 +2747,6 @@ ALIASES += laplace_options="\
   - theta_0 the initial guess for the Laplace approximation. \
   - tolerance controls the convergence criterion when finding the mode in the Laplace approximation. \
   - max_num_steps maximum number of steps before the Newton solver breaks and returns an error. \
-  - hessian_block_size Block size of Hessian of log likelihood w.r.t latent Gaussian variable theta. \
   - solver Type of Newton solver. Each corresponds to a distinct choice of B matrix (i.e. application SWM formula): \
     1. computes square-root of negative Hessian. \
     2. computes square-root of covariance matrix. \

stan/math/mix/functor/laplace_marginal_density_estimator.hpp

@@ -1,6 +1,7 @@
 #ifndef STAN_MATH_MIX_FUNCTOR_LAPLACE_MARGINAL_DENSITY_ESTIMATOR_HPP
 #define STAN_MATH_MIX_FUNCTOR_LAPLACE_MARGINAL_DENSITY_ESTIMATOR_HPP
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/generate_laplace_options.hpp>
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
 #include <stan/math/mix/functor/wolfe_line_search.hpp>
 #include <stan/math/rev/meta.hpp>
@@ -10,6 +11,7 @@
 #include <stan/math/mix/functor/barzilai_borwein_step_size.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/quad_form_diag.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/iter_tuple_nested.hpp>
 #include <unsupported/Eigen/MatrixFunctions>
 #include <cmath>
@@ -30,7 +32,7 @@ namespace math {
  */
 struct laplace_options_base {
   /* Size of the blocks in block diagonal hessian*/
-  int hessian_block_size{1};  // 0
+  int hessian_block_size{internal::laplace_default_hessian_block_size};  // 0
   /**
    * Which linear solver to use inside the Newton step.
    *
@@ -46,19 +48,20 @@ struct laplace_options_base {
    *    `Sigma = K_root * K_root^T` and form `B = I + K_root^T * W * K_root`.
    * 3. General LU: form `B = I + Sigma * W` and factorize with LU.
    */
-  int solver{1};  // 1
+  int solver{internal::laplace_default_solver};  // 1
   /**
    * Iterations end when the absolute change in the optimization objective
    * is less than this tolerance.
    *
    * Note: the objective used for convergence is the one optimized by the
    * Newton/Wolfe loop (not the final Laplace-corrected log marginal density).
    */
-  double tolerance{1.49012e-08};  // 2
+  double tolerance{internal::laplace_default_tolerance};  // 2
   /* Maximum number of steps*/
-  int max_num_steps{500};                   // 3
-  int allow_fallthrough{true};              // 4
-  laplace_line_search_options line_search;  // 5
+  int max_num_steps{internal::laplace_default_max_num_steps};          // 3
+  int allow_fallthrough{internal::laplace_default_allow_fallthrough};  // 4
+  laplace_line_search_options line_search{
+      internal::laplace_default_max_steps_line_search};  // 5
   laplace_options_base() = default;
   laplace_options_base(int hessian_block_size_, int solver_, double tolerance_,
                        int max_num_steps_, bool allow_fallthrough_,
@@ -75,7 +78,13 @@ template <bool HasInitTheta>
 struct laplace_options;
 
 template <>
-struct laplace_options<false> : public laplace_options_base {};
+struct laplace_options<false> : public laplace_options_base {
+  laplace_options() = default;
+
+  explicit laplace_options(int hessian_block_size_) {
+    hessian_block_size = hessian_block_size_;
+  }
+};
 
 template <>
 struct laplace_options<true> : public laplace_options_base {
@@ -89,25 +98,12 @@ struct laplace_options<true> : public laplace_options_base {
       : laplace_options_base(hessian_block_size_, solver_, tolerance_,
                              max_num_steps_, allow_fallthrough_,
                              max_steps_line_search_),
-        theta_0(std::forward<ThetaVec>(theta_0_)) {}
+        theta_0(value_of(std::forward<ThetaVec>(theta_0_))) {}
 };
 
 using laplace_options_default = laplace_options<false>;
 using laplace_options_user_supplied = laplace_options<true>;
 
-/**
- * User function for generating laplace options tuple
- * @param theta_0_size Size of user supplied initial theta
- * @return tuple representing laplace options exposed to user.
- */
-inline auto generate_laplace_options(int theta_0_size) {
-  auto ops = laplace_options_default{};
-  return std::make_tuple(
-      Eigen::VectorXd::Zero(theta_0_size).eval(),  // 0 -> 6
-      ops.tolerance, ops.max_num_steps, ops.hessian_block_size, ops.solver,
-      ops.line_search.max_iterations, static_cast<int>(ops.allow_fallthrough));
-}
-
 namespace internal {
 
 template <typename Options>
@@ -135,41 +131,36 @@ inline constexpr auto tuple_to_laplace_options(Options&& ops) {
           "the laplace approximation.");
     }
     if constexpr (!stan::is_inner_tuple_type_v<3, Ops, int>) {
-      static_assert(
-          sizeof(std::decay_t<Ops>*) == 0,
-          "ERROR:(laplace_marginal_lpdf) The fifth laplace argument is "
-          "expected to be an int representing the hessian block size.");
-    }
-    if constexpr (!stan::is_inner_tuple_type_v<4, Ops, int>) {
       static_assert(
           sizeof(std::decay_t<Ops>*) == 0,
           "ERROR:(laplace_marginal_lpdf) The fourth laplace argument is "
           "expected to be an int representing the solver.");
     }
-    if constexpr (!stan::is_inner_tuple_type_v<5, Ops, int>) {
+    if constexpr (!stan::is_inner_tuple_type_v<4, Ops, int>) {
       static_assert(
           sizeof(std::decay_t<Ops>*) == 0,
-          "ERROR:(laplace_marginal_lpdf) The sixth laplace argument is "
+          "ERROR:(laplace_marginal_lpdf) The fifth laplace argument is "
           "expected to be an int representing the max steps for the laplace "
           "approximaton's wolfe line search.");
     }
     constexpr bool is_fallthrough
         = stan::is_inner_tuple_type_v<
-              6, Ops, int> || stan::is_inner_tuple_type_v<6, Ops, bool>;
+              5, Ops, int> || stan::is_inner_tuple_type_v<5, Ops, bool>;
     if constexpr (!is_fallthrough) {
       static_assert(
           sizeof(std::decay_t<Ops>*) == 0,
-          "ERROR:(laplace_marginal_lpdf) The seventh laplace argument is "
+          "ERROR:(laplace_marginal_lpdf) The sixth laplace argument is "
           "expected to be an int representing allow fallthrough (0/1).");
     }
+    auto defaults = laplace_options_default{};
     return laplace_options_user_supplied{
         value_of(std::get<0>(std::forward<Ops>(ops))),
         std::get<1>(ops),
         std::get<2>(ops),
+        defaults.hessian_block_size,
         std::get<3>(ops),
         std::get<4>(ops),
-        std::get<5>(ops),
-        (std::get<6>(ops) > 0) ? true : false,
+        (std::get<5>(ops) > 0) ? true : false,
     };
   } else {
     return std::forward<Ops>(ops);

stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp

@@ -24,6 +24,8 @@ namespace math {
  * @param[in] n_samples Vector of number of trials.
  * @param[in] mean the mean of the latent normal variable.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \laplace_options
  * \rng_arg
  * \msg_arg
@@ -32,15 +34,16 @@ template <typename Mean, typename CovarFun, typename CovarArgs,
           typename OpsTuple, typename RNG>
 inline Eigen::VectorXd laplace_latent_tol_bernoulli_logit_rng(
     const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, OpsTuple&& ops,
-    RNG& rng, std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, OpsTuple&& ops, RNG& rng, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_base_rng(
       bernoulli_logit_likelihood{},
       std::forward_as_tuple(to_vector(y), n_samples, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args),
-      internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops)), rng,
-      msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), rng, msgs);
 }
 
 /**
@@ -59,20 +62,22 @@ inline Eigen::VectorXd laplace_latent_tol_bernoulli_logit_rng(
  * @param[in] n_samples Vector of number of trials.
  * @param[in] mean the mean of the latent normal variable.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \rng_arg
  * \msg_arg
  */
 template <typename Mean, typename CovarFun, typename CovarArgs, typename RNG>
 inline Eigen::VectorXd laplace_latent_bernoulli_logit_rng(
     const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, RNG& rng,
-    std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, RNG& rng, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_base_rng(
       bernoulli_logit_likelihood{},
       std::forward_as_tuple(to_vector(y), n_samples, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), laplace_options_default{}, rng,
-      msgs);
+      std::forward<CovarArgs>(covar_args), options, rng, msgs);
 }
 
 }  // namespace math

stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp

@@ -31,6 +31,8 @@ namespace math {
  * @param[in] eta Overdisperison parameter.
  * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \laplace_options
  * \rng_arg
  * \msg_arg
@@ -39,16 +41,17 @@ template <typename Eta, typename Mean, typename CovarFun, typename CovarArgs,
           typename OpsTuple, typename RNG>
 inline Eigen::VectorXd laplace_latent_tol_neg_binomial_2_log_rng(
     const std::vector<int>& y, const std::vector<int>& y_index, Eta&& eta,
-    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
-    OpsTuple&& ops, RNG& rng, std::ostream* msgs) {
+    Mean&& mean, int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, OpsTuple&& ops, RNG& rng, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_base_rng(
       neg_binomial_2_log_likelihood{},
       std::forward_as_tuple(std::forward<Eta>(eta), y, y_index,
                             std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args),
-      internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops)), rng,
-      msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), rng, msgs);
 }
 
 /**
@@ -72,22 +75,24 @@ inline Eigen::VectorXd laplace_latent_tol_neg_binomial_2_log_rng(
  * @param[in] eta Overdisperison parameter.
  * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \rng_arg
  * \msg_arg
  */
 template <typename Eta, typename Mean, typename CovarFun, typename CovarArgs,
           typename RNG>
 inline Eigen::VectorXd laplace_latent_neg_binomial_2_log_rng(
     const std::vector<int>& y, const std::vector<int>& y_index, Eta&& eta,
-    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
-    RNG& rng, std::ostream* msgs) {
+    Mean&& mean, int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, RNG& rng, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_base_rng(
       neg_binomial_2_log_likelihood{},
       std::forward_as_tuple(std::forward<Eta>(eta), y, y_index,
                             std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), laplace_options_default{}, rng,
-      msgs);
+      std::forward<CovarArgs>(covar_args), options, rng, msgs);
 }
 
 }  // namespace math

stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp

@@ -26,6 +26,8 @@ namespace math {
  * @param[in] y_index Index indicating which group each observation belongs to.
  * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \laplace_options
  * \rng_arg
  * \msg_arg
@@ -34,15 +36,16 @@ template <typename Mean, typename CovarFun, typename CovarArgs,
           typename OpsTuple, typename RNG>
 inline Eigen::VectorXd laplace_latent_tol_poisson_log_rng(
     const std::vector<int>& y, const std::vector<int>& y_index, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, OpsTuple&& ops,
-    RNG& rng, std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, OpsTuple&& ops, RNG& rng, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_base_rng(
       poisson_log_likelihood{},
       std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args),
-      internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops)), rng,
-      msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), rng, msgs);
 }
 
 /**
@@ -62,20 +65,22 @@ inline Eigen::VectorXd laplace_latent_tol_poisson_log_rng(
  * @param[in] y_index Index indicating which group each observation belongs to.
  * @param[in] mean The mean of the latent normal variable.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \rng_arg
  * \msg_arg
  */
 template <typename CovarFun, typename CovarArgs, typename RNG, typename Mean>
 inline Eigen::VectorXd laplace_latent_poisson_log_rng(
     const std::vector<int>& y, const std::vector<int>& y_index, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, RNG& rng,
-    std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, RNG& rng, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_base_rng(
       poisson_log_likelihood{},
       std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), laplace_options_default{}, rng,
-      msgs);
+      std::forward<CovarArgs>(covar_args), options, rng, msgs);
 }
 
 }  // namespace math

stan/math/mix/prob/laplace_latent_rng.hpp

@@ -25,22 +25,26 @@ namespace math {
  * @param[in] L_f Function that returns log likelihood.
  * @param[in] ll_args Arguments for likelihood function.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \laplace_options
  * \rng_arg
  * \msg_arg
  */
 template <typename LLFunc, typename LLArgs, typename CovarFun,
           typename CovarArgs, typename OpsTuple, typename RNG>
 inline auto laplace_latent_tol_rng(LLFunc&& L_f, LLArgs&& ll_args,
+                                   int hessian_block_size,
                                    CovarFun&& covariance_function,
                                    CovarArgs&& covar_args, OpsTuple&& ops,
                                    RNG& rng, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_base_rng(
       std::forward<LLFunc>(L_f), std::forward<LLArgs>(ll_args),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args),
-      internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops)), rng,
-      msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), rng, msgs);
 }
 
 /**
@@ -58,20 +62,23 @@ inline auto laplace_latent_tol_rng(LLFunc&& L_f, LLArgs&& ll_args,
  * @param[in] L_f Function that returns log likelihood.
  * @param[in] ll_args Arguments for likelihood function.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \rng_arg
  * \msg_arg
  */
 template <typename LLFunc, typename LLArgs, typename CovarFun,
           typename CovarArgs, typename RNG>
 inline auto laplace_latent_rng(LLFunc&& L_f, LLArgs&& ll_args,
+                               int hessian_block_size,
                                CovarFun&& covariance_function,
                                CovarArgs&& covar_args, RNG& rng,
                                std::ostream* msgs) {
-  return laplace_base_rng(std::forward<LLFunc>(L_f),
-                          std::forward<LLArgs>(ll_args),
-                          std::forward<CovarFun>(covariance_function),
-                          std::forward<CovarArgs>(covar_args),
-                          laplace_options_default{}, rng, msgs);
+  auto options = laplace_options_default{hessian_block_size};
+  return laplace_base_rng(
+      std::forward<LLFunc>(L_f), std::forward<LLArgs>(ll_args),
+      std::forward<CovarFun>(covariance_function),
+      std::forward<CovarArgs>(covar_args), options, rng, msgs);
 }
 
 }  // namespace math