diff --git a/doxygen/doxygen.cfg b/doxygen/doxygen.cfg
index f6e2a29f0dc..28942e98d4a 100644
--- a/doxygen/doxygen.cfg
+++ b/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. \
diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp
index d7a76dd4a92..d91f5e3df32 100644
--- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp
+++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp
@@ -10,6 +10,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>
@@ -75,7 +76,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,7 +96,7 @@ 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>;
@@ -102,10 +109,25 @@ using laplace_options_user_supplied = laplace_options<true>;
  */
 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));
+  return std::make_tuple(Eigen::VectorXd::Zero(theta_0_size).eval(),
+                         ops.tolerance, ops.max_num_steps, ops.solver,
+                         ops.line_search.max_iterations,
+                         static_cast<int>(ops.allow_fallthrough));
+}
+
+/**
+ * User function for generating laplace options tuple
+ * @tparam ThetaVec An Eigen vector type for user supplied initial theta
+ * @param theta_0 User supplied initial theta
+ * @return tuple representing laplace options exposed to user.
+ */
+template <typename ThetaVec, require_eigen_t<ThetaVec>* = nullptr>
+inline auto generate_laplace_options(ThetaVec&& theta_0) {
+  auto ops = laplace_options_default{};
+  return std::make_tuple(std::forward<ThetaVec>(theta_0), ops.tolerance,
+                         ops.max_num_steps, ops.solver,
+                         ops.line_search.max_iterations,
+                         static_cast<int>(ops.allow_fallthrough));
 }
 
 namespace internal {
@@ -135,41 +157,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);
diff --git a/stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp b/stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp
index 2d73059fc24..973aa074ac9 100644
--- a/stan/math/mix/prob/laplace_latent_bernoulli_logit_rng.hpp
+++ b/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
diff --git a/stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp b/stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp
index 3fdb669afc6..3cd8df1fbfd 100644
--- a/stan/math/mix/prob/laplace_latent_neg_binomial_2_log_rng.hpp
+++ b/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,6 +75,8 @@ 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
  */
@@ -79,15 +84,15 @@ 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
diff --git a/stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp b/stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp
index d2bbc33a157..b330d0fee5e 100644
--- a/stan/math/mix/prob/laplace_latent_poisson_log_rng.hpp
+++ b/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
diff --git a/stan/math/mix/prob/laplace_latent_rng.hpp b/stan/math/mix/prob/laplace_latent_rng.hpp
index d7aedc47969..9b09d850024 100644
--- a/stan/math/mix/prob/laplace_latent_rng.hpp
+++ b/stan/math/mix/prob/laplace_latent_rng.hpp
@@ -25,6 +25,8 @@ 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
@@ -32,15 +34,17 @@ namespace math {
 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
diff --git a/stan/math/mix/prob/laplace_marginal.hpp b/stan/math/mix/prob/laplace_marginal.hpp
index 7874eef753f..169eb256043 100644
--- a/stan/math/mix/prob/laplace_marginal.hpp
+++ b/stan/math/mix/prob/laplace_marginal.hpp
@@ -22,20 +22,25 @@ namespace math {
  * @param[in] L_f a function which returns the log likelihood.
  * @param[in] l_args A tuple of arguments to pass to the log likelihood.
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \laplace_options
  * \msg_arg
  */
 template <bool propto = false, typename LFun, typename LArgs, typename CovarFun,
           typename CovarArgs, typename OpsTuple>
 inline auto laplace_marginal_tol(LFun&& L_f, LArgs&& l_args,
+                                 int hessian_block_size,
                                  CovarFun&& covariance_function,
                                  CovarArgs&& covar_args, OpsTuple&& ops,
                                  std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_marginal_density(
       std::forward<LFun>(L_f), std::forward<LArgs>(l_args),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args),
-      internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops)), msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), msgs);
 }
 
 /**
@@ -53,17 +58,20 @@ inline auto laplace_marginal_tol(LFun&& L_f, LArgs&& l_args,
  * @param[in] L_f a function which returns the log likelihood.
  * @param[in] l_args A tuple of arguments to pass to the log likelihood
  * \laplace_common_args
+ * @param[in] hessian_block_size Block size for the Hessian approximation with
+ * respect to the latent gaussian variable theta.
  * \msg_arg
  */
 template <bool propto = false, typename LFun, typename LArgs, typename CovarFun,
           typename CovarArgs>
-inline auto laplace_marginal(LFun&& L_f, LArgs&& l_args,
+inline auto laplace_marginal(LFun&& L_f, LArgs&& l_args, int hessian_block_size,
                              CovarFun&& covariance_function,
                              CovarArgs&& covar_args, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_marginal_density(
       std::forward<LFun>(L_f), std::forward<LArgs>(l_args),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), laplace_options_default{}, msgs);
+      std::forward<CovarArgs>(covar_args), options, msgs);
 }
 
 }  // namespace math
diff --git a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
index 357788dd983..f8d521f7392 100644
--- a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
@@ -49,6 +49,8 @@ struct bernoulli_logit_likelihood {
  * statistics.
  * @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
  * \msg_arg
  */
@@ -56,14 +58,16 @@ template <bool propto = false, typename Mean, typename CovarFun,
           typename CovarArgs, typename OpsTuple>
 inline auto laplace_marginal_tol_bernoulli_logit_lpmf(
     const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, OpsTuple&& ops,
-    std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, OpsTuple&& ops, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_marginal_density(
       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)), msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), msgs);
 }
 
 /**
@@ -81,19 +85,22 @@ inline auto laplace_marginal_tol_bernoulli_logit_lpmf(
  * statistics.
  * @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.
  * \msg_arg
  */
 template <bool propto = false, typename Mean, typename CovarFun,
           typename CovarArgs>
 inline auto laplace_marginal_bernoulli_logit_lpmf(
     const std::vector<int>& y, const std::vector<int>& n_samples, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args,
-    std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_marginal_density(
       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{}, msgs);
+      std::forward<CovarArgs>(covar_args), options, msgs);
 }
 
 }  // namespace math
diff --git a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
index a328061ad23..843f9d2e4dd 100644
--- a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
@@ -73,6 +73,8 @@ struct neg_binomial_2_log_likelihood {
  * @param[in] eta non-marginalized model parameters for the likelihood.
  * @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
  * \msg_arg
  */
@@ -80,14 +82,16 @@ template <bool propto = false, typename Eta, typename Mean, typename CovarFun,
           typename CovarArgs, typename OpsTuple>
 inline auto laplace_marginal_tol_neg_binomial_2_log_lpmf(
     const std::vector<int>& y, const std::vector<int>& y_index, const Eta& eta,
-    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
-    OpsTuple&& ops, std::ostream* msgs) {
+    Mean&& mean, int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, OpsTuple&& ops, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_marginal_density(
       neg_binomial_2_log_likelihood{},
       std::forward_as_tuple(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)), msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), msgs);
 }
 
 /**
@@ -106,19 +110,22 @@ inline auto laplace_marginal_tol_neg_binomial_2_log_lpmf(
  * @param[in] eta Parameter argument for likelihood function.
  * @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.
  * \msg_arg
  */
 template <bool propto = false, typename Eta, typename Mean, typename CovarFun,
           typename CovarArgs>
 inline auto laplace_marginal_neg_binomial_2_log_lpmf(
     const std::vector<int>& y, const std::vector<int>& y_index, const Eta& eta,
-    Mean&& mean, CovarFun&& covariance_function, CovarArgs&& covar_args,
-    std::ostream* msgs) {
+    Mean&& mean, int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_marginal_density(
       neg_binomial_2_log_likelihood{},
       std::forward_as_tuple(eta, y, y_index, std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), laplace_options_default{}, msgs);
+      std::forward<CovarArgs>(covar_args), options, msgs);
 }
 
 struct neg_binomial_2_log_likelihood_summary {
@@ -167,6 +174,8 @@ struct neg_binomial_2_log_likelihood_summary {
  * @param[in] eta non-marginalized model parameters for the likelihood.
  * @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
  * \msg_arg
  */
@@ -175,15 +184,17 @@ template <bool propto = false, typename Eta, typename Mean, typename CovarFun,
 inline auto laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
     const std::vector<int>& y, const std::vector<int>& n_per_group,
     const std::vector<int>& counts_per_group, const Eta& eta, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, OpsTuple&& ops,
-    std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, OpsTuple&& ops, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_marginal_density(
       neg_binomial_2_log_likelihood_summary{},
       std::forward_as_tuple(eta, y, n_per_group, counts_per_group,
                             std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args),
-      internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops)), msgs);
+      std::forward<CovarArgs>(covar_args), std::move(options), msgs);
 }
 
 /**
@@ -202,6 +213,8 @@ inline auto laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
  * @param[in] eta non-marginalized model parameters for the likelihood.
  * @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.
  * \msg_arg
  */
 template <bool propto = false, typename Eta, typename Mean, typename CovarFun,
@@ -209,14 +222,15 @@ template <bool propto = false, typename Eta, typename Mean, typename CovarFun,
 inline auto laplace_marginal_neg_binomial_2_log_summary_lpmf(
     const std::vector<int>& y, const std::vector<int>& n_per_group,
     const std::vector<int>& counts_per_group, const Eta& eta, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args,
-    std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_marginal_density(
       neg_binomial_2_log_likelihood_summary{},
       std::forward_as_tuple(eta, y, n_per_group, counts_per_group,
                             std::forward<Mean>(mean)),
       std::forward<CovarFun>(covariance_function),
-      std::forward<CovarArgs>(covar_args), laplace_options_default{}, msgs);
+      std::forward<CovarArgs>(covar_args), options, msgs);
 }
 
 }  // namespace math
diff --git a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
index 5f13c8eee12..8078ce72dd8 100644
--- a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
@@ -65,6 +65,8 @@ struct poisson_log_likelihood {
  * @param[in] y_index group to which each observation belongs
  * @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
  * \msg_arg
  */
@@ -72,13 +74,16 @@ template <bool propto = false, typename Mean, typename CovarFun,
           typename CovarArgs, typename OpsTuple>
 inline auto laplace_marginal_tol_poisson_log_lpmf(
     const std::vector<int>& y, const std::vector<int>& y_index, Mean&& mean,
-    CovarFun&& covariance_function, CovarArgs&& covar_args, OpsTuple&& ops,
-    std::ostream* msgs) {
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, OpsTuple&& ops, std::ostream* msgs) {
+  auto options
+      = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
+  options.hessian_block_size = hessian_block_size;
   return laplace_marginal_density(
       poisson_log_likelihood{},
       std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
       covariance_function, std::forward<CovarArgs>(covar_args),
-      internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops)), msgs);
+      std::move(options), msgs);
 }
 
 /**
@@ -95,21 +100,21 @@ inline auto laplace_marginal_tol_poisson_log_lpmf(
  * @param[in] y_index group to which each observation belongs
  * @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.
  * \msg_arg
  */
 template <bool propto = false, typename CovarFun, typename CovarArgs,
           typename Mean>
-inline auto laplace_marginal_poisson_log_lpmf(const std::vector<int>& y,
-                                              const std::vector<int>& y_index,
-                                              Mean&& mean,
-                                              CovarFun&& covariance_function,
-                                              CovarArgs&& covar_args,
-                                              std::ostream* msgs) {
+inline auto laplace_marginal_poisson_log_lpmf(
+    const std::vector<int>& y, const std::vector<int>& y_index, Mean&& mean,
+    int hessian_block_size, CovarFun&& covariance_function,
+    CovarArgs&& covar_args, std::ostream* msgs) {
+  auto options = laplace_options_default{hessian_block_size};
   return laplace_marginal_density(
       poisson_log_likelihood{},
       std::forward_as_tuple(y, y_index, std::forward<Mean>(mean)),
-      covariance_function, std::forward<CovarArgs>(covar_args),
-      laplace_options_default{}, msgs);
+      covariance_function, std::forward<CovarArgs>(covar_args), options, msgs);
 }
 
 }  // namespace math
diff --git a/test/unit/math/laplace/aki_ex_test.cpp b/test/unit/math/laplace/aki_ex_test.cpp
index 284ee8c48a3..60047348fd8 100644
--- a/test/unit/math/laplace/aki_ex_test.cpp
+++ b/test/unit/math/laplace/aki_ex_test.cpp
@@ -85,7 +85,7 @@ TEST(WriteArrayBodySimple, ExceededIteration) {
     double ll_laplace_val{0};
     try {
       ll_laplace_val = stan::math::laplace_marginal(
-          poisson_re_log_ll_functor(), std::forward_as_tuple(y, mu),
+          poisson_re_log_ll_functor(), std::forward_as_tuple(y, mu), 1,
           cov_fun_functor(), std::tuple<double, int>(sigmaz, 1), pstream);
     } catch (const std::domain_error& e) {
       // Log bad values to CSV files
@@ -135,7 +135,7 @@ TEST(WriteArrayBodySimple, ExecutesBodyWithHardcodedData) {
       try {
         ll_laplace_val = stan::math::laplace_marginal(
             poisson_re_log_ll_functor(),
-            std::forward_as_tuple(y[i - 1], mu[i - 1]), cov_fun_functor(),
+            std::forward_as_tuple(y[i - 1], mu[i - 1]), 1, cov_fun_functor(),
             std::tuple<double, int>(sigmaz, 1), pstream);
       } catch (const std::domain_error& e) {
         ADD_FAILURE() << "LAPLACE FAILURE: y and mu for i = " << i << ": ("
@@ -164,7 +164,7 @@ TEST(WriteArrayBodySimple, ExecutesBodyWithHardcodedData) {
       }
     }
     auto ll_laplace_all = stan::math::laplace_marginal(
-        poisson_re_log_ll_functor(), std::forward_as_tuple(y, mu),
+        poisson_re_log_ll_functor(), std::forward_as_tuple(y, mu), 1,
         cov_fun_functor(), std::tuple<double, int>(sigmaz, N), pstream);
     stan::test::relative_tolerance sum_rel_tol(3e-2);
     expect_near_rel("sum laplace vs integrated sum", ll_laplace,
diff --git a/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp b/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
index 4442bcf7066..c37eb4fb392 100644
--- a/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
@@ -76,7 +76,7 @@ TEST(laplace_bernoulli_logit_rng, two_dim_diag) {
   boost::random::mt19937 rng;
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_bernoulli_logit_rng(
-      sums, n_samples, mean, diagonal_kernel_functor{},
+      sums, n_samples, mean, 1, diagonal_kernel_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   // Compute exact mean and covariance
diff --git a/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp
index 89ea619819f..58317a5ebc2 100644
--- a/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_bernoulli_logit_lpmf_test.cpp
@@ -44,7 +44,7 @@ TEST_P(laplace_marginal_bernoulli_logit_lpmf, phi_dim500) {
   Eigen::Matrix<double, Eigen::Dynamic, 1> phi_dbl{{1.6, 1}};
   using stan::math::test::sqr_exp_kernel_functor;
   double target = laplace_marginal_bernoulli_logit_lpmf(
-      y, n_samples, 0, sqr_exp_kernel_functor{},
+      y, n_samples, 0, hessian_block_size, sqr_exp_kernel_functor{},
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), nullptr);
   // Benchmark against gpstuff.
   constexpr double tol = 8e-4;
@@ -56,10 +56,10 @@ TEST_P(laplace_marginal_bernoulli_logit_lpmf, phi_dim500) {
   auto f = [&](auto&& alpha, auto&& rho) {
     try {
       return laplace_marginal_tol_bernoulli_logit_lpmf(
-          y, n_samples, mean, sqr_exp_kernel_functor{},
+          y, n_samples, mean, hessian_block_size, sqr_exp_kernel_functor{},
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
diff --git a/test/unit/math/laplace/laplace_marginal_lpdf_moto_test.cpp b/test/unit/math/laplace/laplace_marginal_lpdf_moto_test.cpp
index 7c43c1da675..f00974cf3c4 100644
--- a/test/unit/math/laplace/laplace_marginal_lpdf_moto_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_lpdf_moto_test.cpp
@@ -148,12 +148,12 @@ TEST_P(laplace_motorcyle_gp_test, gp_motorcycle_val) {
   LAPLACE_SKIP_IF_INVALID_TEST_COMBO(hessian_block_size, dim_theta);
 
   laplace_marginal_tol<false>(
-      normal_likelihood{}, std::forward_as_tuple(y, n_obs),
+      normal_likelihood{}, std::forward_as_tuple(y, n_obs), hessian_block_size,
       covariance_motorcycle_functor{},
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1), phi_dbl(2), phi_dbl(3),
                             n_obs),
-      std::make_tuple(theta0, tolerance, max_num_steps, hessian_block_size,
-                      solver_num, max_steps_line_search, true),
+      std::make_tuple(theta0, tolerance, max_num_steps, solver_num,
+                      max_steps_line_search, true),
       &output_stream);
 }
 
@@ -182,11 +182,11 @@ TEST_P(laplace_motorcyle_gp_test, gp_motorcycle_ad) {
     try {
       return laplace_marginal_tol<false>(
           normal_likelihood{}, std::forward_as_tuple(y, n_obs),
-          covariance_motorcycle_functor{},
+          hessian_block_size, covariance_motorcycle_functor{},
           std::forward_as_tuple(x, phi_01_v(0), phi_01_v(1), phi_rest_v(0),
                                 phi_rest_v(1), n_obs),
-          std::make_tuple(theta0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -256,11 +256,11 @@ TEST_P(laplace_motorcyle_gp_test, gp_motorcycle2_val) {
   LAPLACE_SKIP_IF_INVALID_TEST_COMBO(hessian_block_size, dim_theta);
   laplace_marginal_tol<false>(
       normal_likelihood2{}, std::forward_as_tuple(y, n_obs, sigma_global),
-      covariance_motorcycle_functor{},
+      hessian_block_size, covariance_motorcycle_functor{},
       std::forward_as_tuple(x, length_scale_f, length_scale_g, sigma_f, sigma_g,
                             n_obs),
-      std::make_tuple(theta0, tolerance, max_num_steps, hessian_block_size,
-                      solver_num, max_steps_line_search, true),
+      std::make_tuple(theta0, tolerance, max_num_steps, solver_num,
+                      max_steps_line_search, true),
       &output_stream);
 }
 
@@ -289,11 +289,11 @@ TEST_P(laplace_motorcyle_gp_test, gp_motorcycle2_ad) {
     try {
       return laplace_marginal_tol<false>(
           normal_likelihood2{}, std::forward_as_tuple(y, n_obs, sigma_global_v),
-          covariance_motorcycle_functor{},
+          hessian_block_size, covariance_motorcycle_functor{},
           std::forward_as_tuple(x, length_scale_v(0), length_scale_v(1),
                                 sigma_v(0), sigma_v(1), n_obs),
-          std::make_tuple(theta0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
diff --git a/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp b/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp
index 1b1f6549adf..f494662b5a0 100644
--- a/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_lpdf_test.cpp
@@ -55,7 +55,7 @@ TEST_P(laplace_marginal_lpdf, poisson_log_phi_dim_2) {
 
   double target = laplace_marginal<false>(
       poisson_log_likelihood2{}, std::forward_as_tuple(sums),
-      stan::math::test::squared_kernel_functor{},
+      hessian_block_size, stan::math::test::squared_kernel_functor{},
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), &output_stream);
 
   // TODO(Charles): benchmark target against gpstuff.
@@ -72,10 +72,10 @@ TEST_P(laplace_marginal_lpdf, poisson_log_phi_dim_2) {
 
     target = laplace_marginal_tol<false>(
         poisson_log_likelihood2{}, std::forward_as_tuple(sums),
-        stan::math::test::squared_kernel_functor{},
+        hessian_block_size, stan::math::test::squared_kernel_functor{},
         std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)),
-        std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                        solver, max_steps_line_search, true),
+        std::make_tuple(theta_0, tolerance, max_num_steps, solver,
+                        max_steps_line_search, true),
         &output_stream);
     EXPECT_NEAR(-2.53056, value_of(target), tol);
   }
@@ -91,10 +91,10 @@ TEST_P(laplace_marginal_lpdf, poisson_log_phi_dim_2) {
     try {
       return laplace_marginal_tol<false>(
           poisson_log_likelihood2{}, std::forward_as_tuple(sums),
-          stan::math::test::squared_kernel_functor{},
+          hessian_block_size, stan::math::test::squared_kernel_functor{},
           std::forward_as_tuple(x_v, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -134,7 +134,7 @@ TEST_P(laplace_disease_map_test, laplace_marginal) {
   {
     double marginal_density = laplace_marginal<false>(
         poisson_log_exposure_likelihood{}, std::forward_as_tuple(ye, y),
-        stan::math::test::sqr_exp_kernel_functor{},
+        hessian_block_size, stan::math::test::sqr_exp_kernel_functor{},
         std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), &output_stream);
 
     constexpr double tol = 6e-4;
@@ -149,10 +149,10 @@ TEST_P(laplace_disease_map_test, laplace_marginal) {
     try {
       return laplace_marginal_tol<false>(
           poisson_log_exposure_likelihood{}, std::forward_as_tuple(ye, y),
-          stan::math::test::sqr_exp_kernel_functor{},
+          hessian_block_size, stan::math::test::sqr_exp_kernel_functor{},
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -203,7 +203,7 @@ TEST_P(laplace_marginal_lpdf, bernoulli_logit_phi_dim500) {
   LAPLACE_SKIP_IF_INVALID_TEST_COMBO(hessian_block_size, dim_theta);
   double target = laplace_marginal<false>(
       bernoulli_logit_likelihood{}, std::forward_as_tuple(y),
-      stan::math::test::sqr_exp_kernel_functor{},
+      hessian_block_size, stan::math::test::sqr_exp_kernel_functor{},
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)), &output_stream);
 
   constexpr double tol = 3e-4;
@@ -218,10 +218,10 @@ TEST_P(laplace_marginal_lpdf, bernoulli_logit_phi_dim500) {
     try {
       return laplace_marginal_tol<false>(
           bernoulli_logit_likelihood{}, std::forward_as_tuple(y),
-          stan::math::test::sqr_exp_kernel_functor{},
+          hessian_block_size, stan::math::test::sqr_exp_kernel_functor{},
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
diff --git a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp
index 8572cf102b3..cf2c078a2f3 100644
--- a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_lpmf_test.cpp
@@ -46,10 +46,11 @@ TEST_P(laplace_marginal_neg_binomial_log_lpmf, phi_dim_2) {
   auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
     try {
       return laplace_marginal_tol_neg_binomial_2_log_lpmf(
-          y, y_index, eta, 0, stan::math::test::squared_kernel_functor{},
+          y, y_index, eta, 0, hessian_block_size,
+          stan::math::test::squared_kernel_functor{},
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -86,20 +87,22 @@ TEST_P(laplace_disease_map_test, laplace_marginal_neg_binomial_2_log_lpmf) {
   constexpr int max_num_steps = 1000;
   auto smoke = [&](auto&& alpha, auto&& rho, auto&& eta_arg) {
     return laplace_marginal_tol_neg_binomial_2_log_lpmf(
-        y, y_index, eta_arg, mean, stan::math::test::sqr_exp_kernel_functor{},
+        y, y_index, eta_arg, mean, hessian_block_size,
+        stan::math::test::sqr_exp_kernel_functor{},
         std::forward_as_tuple(x, alpha, rho),
-        std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                        solver_num, max_steps_line_search, true),
+        std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                        max_steps_line_search, true),
         &output_stream);
   };
   smoke(phi_dbl[0], phi_dbl[1], eta);
   auto f = [&](auto&& alpha, auto&& rho, auto&& eta_arg) {
     try {
       return laplace_marginal_tol_neg_binomial_2_log_lpmf(
-          y, y_index, eta_arg, mean, stan::math::test::sqr_exp_kernel_functor{},
+          y, y_index, eta_arg, mean, hessian_block_size,
+          stan::math::test::sqr_exp_kernel_functor{},
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
diff --git a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp
index 5c3ad8214f8..fba6e8ac345 100644
--- a/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_neg_binomial_log_summary_lpmf_test.cpp
@@ -52,11 +52,11 @@ TEST_P(laplace_marginal_neg_binomial_log_summary_lpmf, phi_dim_2) {
   auto f = [&](auto&& alpha, auto&& rho, auto&& eta) {
     try {
       return laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
-          y, n_per_group, counts_per_group, eta, 0,
+          y, n_per_group, counts_per_group, eta, 0, hessian_block_size,
           stan::math::test::squared_kernel_functor{},
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -100,22 +100,22 @@ TEST_P(laplace_disease_map_test,
   constexpr int max_num_steps = 1000;
   auto smoke = [&](auto&& alpha, auto&& rho, auto&& eta_arg) {
     return laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
-        y, n_per_group, counts_per_group, eta_arg, 0,
+        y, n_per_group, counts_per_group, eta_arg, 0, hessian_block_size,
         stan::math::test::sqr_exp_kernel_functor{},
         std::forward_as_tuple(x, alpha, rho),
-        std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                        solver_num, max_steps_line_search, true),
+        std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                        max_steps_line_search, true),
         &output_stream);
   };
   smoke(phi_dbl[0], phi_dbl[1], eta);
   auto f = [&](auto&& alpha, auto&& rho, auto&& eta_arg) {
     try {
       return laplace_marginal_tol_neg_binomial_2_log_summary_lpmf(
-          y, n_per_group, counts_per_group, eta_arg, mean,
+          y, n_per_group, counts_per_group, eta_arg, mean, hessian_block_size,
           stan::math::test::sqr_exp_kernel_functor{},
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
diff --git a/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp b/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp
index 0b40f33b8da..6b9d86db888 100644
--- a/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp
+++ b/test/unit/math/laplace/laplace_marginal_poisson_log_lpmf_test.cpp
@@ -59,9 +59,10 @@ TEST_P(laplace_marginal_poisson_log_lpmf, phi_dim_2) {
   auto f = [&](auto&& alpha, auto&& rho) {
     try {
       return laplace_marginal_tol_poisson_log_lpmf(
-          y, y_index, 0, sq_kernel, std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, allow_fallthrough),
+          y, y_index, 0, hessian_block_size, sq_kernel,
+          std::forward_as_tuple(x, alpha, rho),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, allow_fallthrough),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -124,9 +125,10 @@ TEST_P(laplace_marginal_poisson_log_lpmf, log_phi_dim_2) {
   auto f = [&](auto&& alpha, auto&& rho) {
     try {
       return laplace_marginal_tol_poisson_log_lpmf(
-          y, y_index, log(ye), sq_kernel, std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          y, y_index, log(ye), hessian_block_size, sq_kernel,
+          std::forward_as_tuple(x, alpha, rho),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -174,10 +176,10 @@ TEST_P(laplace_marginal_poisson_log_lpmf, mean_argument) {
   constexpr double tolerance = 1e-12;
   constexpr int max_num_steps = 500;
   double marginal_density = laplace_marginal_tol_poisson_log_lpmf(
-      y, y_index, mu, diag_covariance{},
+      y, y_index, mu, hessian_block_size, diag_covariance{},
       std::tuple<double, int>(sigmaz, dim_theta),
-      std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                      solver_num, max_steps_line_search, true),
+      std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                      max_steps_line_search, true),
       &output_stream);
 
   EXPECT_FLOAT_EQ(-6.7098737, marginal_density);
@@ -200,10 +202,11 @@ TEST_P(laplace_disease_map_test, laplace_marginal_poisson_log_lpmf) {
   constexpr int max_num_steps = 500;
 
   double marginal_density = laplace_marginal_tol_poisson_log_lpmf(
-      y, y_index, log(ye), stan::math::test::sqr_exp_kernel_functor(),
+      y, y_index, log(ye), hessian_block_size,
+      stan::math::test::sqr_exp_kernel_functor(),
       std::forward_as_tuple(x, phi_dbl(0), phi_dbl(1)),
-      std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                      solver_num, max_steps_line_search, true),
+      std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                      max_steps_line_search, true),
       &output_stream);
 
   double tol = 6e-4;
@@ -212,10 +215,11 @@ TEST_P(laplace_disease_map_test, laplace_marginal_poisson_log_lpmf) {
   auto f = [&](auto&& alpha, auto&& rho) {
     try {
       return laplace_marginal_tol_poisson_log_lpmf(
-          y, y_index, log(ye), stan::math::test::sqr_exp_kernel_functor(),
+          y, y_index, log(ye), hessian_block_size,
+          stan::math::test::sqr_exp_kernel_functor(),
           std::forward_as_tuple(x, alpha, rho),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
diff --git a/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp b/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
index 4ace1fc9bdd..3e248b8874e 100644
--- a/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
@@ -103,7 +103,7 @@ TEST(laplace_latent_neg_binomial_2_log_rng, count_two_dim_diag) {
 
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_neg_binomial_2_log_rng(
-      y, y_index, eta, 0, diagonal_kernel_nb_functor{},
+      y, y_index, eta, 0, 1, diagonal_kernel_nb_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   double tol = 1e-3;
@@ -116,7 +116,7 @@ TEST(laplace_latent_neg_binomial_2_log_rng, count_two_dim_diag) {
   for (int i = 0; i < n_sim; i++) {
     rng.seed(2025 + i);
     Eigen::MatrixXd theta_pred = laplace_latent_neg_binomial_2_log_rng(
-        y, y_index, eta, 0, diagonal_kernel_nb_functor{},
+        y, y_index, eta, 0, 1, diagonal_kernel_nb_functor{},
         std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
     theta_dim0(i) = theta_pred(0);
diff --git a/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp b/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
index c2d358a1586..6c5a13670f0 100644
--- a/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
@@ -23,7 +23,7 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_likelihood) {
   boost::random::mt19937 rng;
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_poisson_log_rng(
-      y, y_index, 0, stan::math::test::diagonal_kernel_functor{},
+      y, y_index, 0, 1, stan::math::test::diagonal_kernel_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   // double tol = 1e-3;
@@ -36,7 +36,7 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_likelihood) {
   for (int i = 0; i < n_sim; i++) {
     rng.seed(2025 + i);
     Eigen::MatrixXd theta_pred = laplace_latent_poisson_log_rng(
-        y, y_index, 0, stan::math::test::diagonal_kernel_functor{},
+        y, y_index, 0, 1, stan::math::test::diagonal_kernel_functor{},
         std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
     theta_dim0(i) = theta_pred(0);
@@ -73,7 +73,7 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_exp_likelihood) {
 
   rng.seed(1954);
   Eigen::MatrixXd theta_pred_exp = laplace_latent_poisson_log_rng(
-      y, y_index, log(ye), stan::math::test::diagonal_kernel_functor{},
+      y, y_index, log(ye), 1, stan::math::test::diagonal_kernel_functor{},
       std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
   EXPECT_NEAR(theta_benchmark(0), theta_pred_exp(0), tol);
@@ -84,7 +84,7 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_exp_likelihood) {
   for (int i = 0; i < n_sim; i++) {
     rng.seed(2025 + i);
     Eigen::MatrixXd theta_pred = laplace_latent_poisson_log_rng(
-        y, y_index, log(ye), stan::math::test::diagonal_kernel_functor{},
+        y, y_index, log(ye), 1, stan::math::test::diagonal_kernel_functor{},
         std::forward_as_tuple(phi(0), phi(1)), rng, nullptr);
 
     theta_dim0(i) = theta_pred(0);
diff --git a/test/unit/math/laplace/laplace_types_test.cpp b/test/unit/math/laplace/laplace_types_test.cpp
index 3bbac39b1c8..c1599aedd8e 100644
--- a/test/unit/math/laplace/laplace_types_test.cpp
+++ b/test/unit/math/laplace/laplace_types_test.cpp
@@ -94,9 +94,9 @@ TEST(laplace, theta_0_as_expression_issue_3196) {
       std::tuple<const std::vector<int>&, Eigen::Matrix<double, -1, 1>>(
           y, stan::math::add(stan::math::add(offset, alpha),
                              stan::math::multiply(X, beta))),
-      cov_fun(), std::tuple<double, int>(sigmaz, N),
+      hessian_block_size, cov_fun(), std::tuple<double, int>(sigmaz, N),
       std::tuple{stan::math::rep_vector(0.0, N), tolerance, max_num_steps,
-                 hessian_block_size, solver_num, max_steps_line_search, 1},
+                 solver_num, max_steps_line_search, 1},
       nullptr));
   auto arena_init = stan::math::to_arena(stan::math::rep_vector(0.0, N));
   EXPECT_NO_THROW(stan::math::laplace_marginal_tol<false>(
@@ -104,9 +104,9 @@ TEST(laplace, theta_0_as_expression_issue_3196) {
       std::tuple<const std::vector<int>&, Eigen::Matrix<double, -1, 1>>(
           y, stan::math::add(stan::math::add(offset, alpha),
                              stan::math::multiply(X, beta))),
-      cov_fun(), std::tuple(sigmaz, N),
-      std::make_tuple(arena_init, tolerance, max_num_steps, hessian_block_size,
-                      solver_num, max_steps_line_search, true),
+      hessian_block_size, cov_fun(), std::tuple(sigmaz, N),
+      std::make_tuple(arena_init, tolerance, max_num_steps, solver_num,
+                      max_steps_line_search, true),
       nullptr));
 }
 
@@ -156,10 +156,10 @@ TEST_P(laplace_types, poisson_log_phi_dim_2_tuple_extended) {
       return laplace_marginal_tol<false>(
           poisson_log_likelihood_tuple_expanded{},
           std::forward_as_tuple(sums, eta1_tuple, eta2, eta3),
-          stan::math::test::squared_kernel_functor{},
+          hessian_block_size, stan::math::test::squared_kernel_functor{},
           std::forward_as_tuple(x, std::make_tuple(phi_dbl(0), phi_dbl(1))),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -222,10 +222,10 @@ TEST_P(laplace_types, poisson_log_phi_dim_2_tuple) {
     try {
       return laplace_marginal_tol<false>(
           poisson_log_likelihood2{}, std::forward_as_tuple(sums),
-          stan::math::test::squared_kernel_functor{},
+          hessian_block_size, stan::math::test::squared_kernel_functor{},
           std::forward_as_tuple(x_v, std::make_tuple(alpha, rho)),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -244,10 +244,10 @@ TEST_P(laplace_types, poisson_log_phi_dim_2_tuple) {
       return laplace_marginal_tol<false>(
           poisson_log_likelihood_tuple{},
           std::forward_as_tuple(sums, std::make_tuple(eta1, eta2)),
-          stan::math::test::squared_kernel_functor{},
+          hessian_block_size, stan::math::test::squared_kernel_functor{},
           std::forward_as_tuple(x, std::make_tuple(alpha_rho(0), alpha_rho(1))),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;
@@ -326,11 +326,11 @@ TEST_P(laplace_types, poisson_log_phi_dim_2_array_tuple) {
       alpha_tuple.push_back(std::make_tuple(alpha_rho(0), alpha_rho(1)));
       return laplace_marginal_tol<false>(
           poisson_log_likelihood_array_tuple{},
-          std::forward_as_tuple(sums, eta_tuple),
+          std::forward_as_tuple(sums, eta_tuple), hessian_block_size,
           stan::math::test::squared_kernel_functor{},
           std::forward_as_tuple(x, alpha_tuple),
-          std::make_tuple(theta_0, tolerance, max_num_steps, hessian_block_size,
-                          solver_num, max_steps_line_search, true),
+          std::make_tuple(theta_0, tolerance, max_num_steps, solver_num,
+                          max_steps_line_search, true),
           &output_stream);
     } catch (const std::exception& e) {
       std::stringstream fail_msg;