Using Laplace approximation for logistic normal multinomial model

Yingnan_Gao · May 22, 2026, 5:51am

Hello Stan community,

I posted a thread regarding some difficult geometries I encountered for fitting logistic normal multinomial model with explicit latent residuals (so that the likelihood can be evaluated for observed counts), on the same day when Stan v.2.39 was released with the embedded Laplace approximation as a new feature.

I tried to adapt the example in the user manual for my LNM model, but I am not sure if I understand the procedures completely. So, following the suggestion by @stemangiola , I will post my code here with my understanding of their roles, and I would like to know if I construct the model with Laplace approximation correctly.

The Laplace approximation integrates out the latent residuals, and to do that it requires two custom functions ll_function and cov_function:

functions{
  // likelihood of observation given latent residuals + parameters:
  real ll_function(
      vector u, // latent ILR residuals
      int M, // number of taxa for convenience
      array[] int y, // observed counts
      vector beta, // ILR mean
      matrix Helmert // canonical Helmert matrix constant
  ) {
    real lp = 0;
    vector[M] mu = Helmert*(beta+u); // transform to CLR composition

    lp += multinomial_logit_lpmf(
      y | mu
    );
    return lp;
  }
  
  matrix cov_function(
    matrix L,
    int M // for some reason I cannot make (L) a tuple to match the signature
  ){
    return L*L';
  }
}

My understanding is that ll_function should calculate the likelihood of data given the parameters and the latent residuals in LNM, and cov_function returns the covariance matrix of the latent residuals. That is, the latent parameters to be integrated out are assumed to follow a 0-mean multivariate normal distribution controlled by the covariance matrix returned by cov_function. This is definitely the case in LNM.

These two functions should be supplied to laplace_marginal as the first and the fourth argument. Their corresponding parameters should be supplied as the second and the fifth argument. The Hessian block size (in LNM it should be M-1) is supplied as the third argument. I implement this step in the partial sum function to enable CPU acceleration.

laplace_marginal should handle the Hessian correction and all other procedures in Laplace approximation automatically and internally.

functions {
  real partial_sum_2_lpmf(
    // Parallel
    array[] int idx_y,
    int start,
    int end,
    
    // General
    array[,] int y, //Sliced, M dimensions

    // Variance-Covariance
    array[] matrix L, // ILR, M-1 dimensions
   
    // Fixed effects
    matrix beta, // ILR-transformed, M-1 dimensions 
    int M, 
    
    // Helmert matrix
    matrix Helmert
    ){
      
      int N = end-start+1; // Number of observations subsetted to this chunk
      
      // mu
      matrix[M-1, N] mu = (rep_matrix(1,N,1) * beta)';
 
      // target log-probability
      real target_lp = 0;
      

      for(n in 1:N){
        target_lp += laplace_marginal(
          ll_function, //custom log-likelihood function with latent residuals
          (M,to_array_1d(y[idx_y[n],]),ysum[idx_y[n]],L[1],mu[,n],Helmert), //parameters without latent residuals
          M-1, // Hessian block size
          cov_function, //covariance function for latent residuals
          (to_matrix(L[1]),M) // cov_function just converts L to VCoV
          );
        }
      return target_lp;
    }
}

The rest of the Stan code is essentially the same as that in my previous post, but with the latent residuals removed from parameters and transformed parameters block.

I reran the LNM model with Laplace approximation on the simulated data (3 taxa and up to 1000 samples) and found the inflated rhat was resolved and I receive no warnings of divergent transition or maximum tree depth reached for runs with >2 samples.

So the new embedded Laplace approximation seem to address my problem completely, provided that I applied it properly in the codes I shared above.

avehtari · May 27, 2026, 9:18am

The release notes mention that 1-tuple is written as (L, ). The additional , is needed to make it tuple (this is same syntax as in Python)

Awesome! Thanks for sharing! (I didn’t have time to verify the code, but if the results match what you expect it’s likely to be fine)

Yingnan_Gao · June 4, 2026, 7:13am

Thanks for tuple solution!

I am benchmarking on larger datasets (e.g., with 10 taxa instead of 3), for some runs I encountered the following warnings:

WARNING(laplace_marginal_density): max number of iterations: 500 exceeded

I think this could be addressed by lifting the iteration limit to a higher number:

transformed data {
  tuple(vector[N], real, int, int, int, int, int) laplace_ops =
    generate_laplace_options(N);

  laplace_ops.3 = 1000; // maximum number of steps for optimizer.
}

and switch to laplace_marginal_tol instead of the default laplace_marginal?

      for(n in 1:N){
        target_lp += laplace_marginal_tol(
          ll_function, //custom log-likelihood function with latent residuals
          (M,to_array_1d(y[idx_y[n],]),ysum[idx_y[n]],L[1],mu[,n],Helmert), //parameters without latent residuals
          M-1, // Hessian block size
          cov_function, //covariance function for latent residuals
          (to_matrix(L[1]),M), // cov_function just converts L to VCoV
          laplace_ops
          );
        }
      return target_lp;

Is there any rule of thumb for setting this limit based on the size of data? Or it has to be case-specific and rely on manual curation?

avehtari · June 4, 2026, 7:26am

No rules, as it’s not just the size of the data but shape of log likelihood, too. You can also try setting theta_init if you have some useful information (has helped a lot in my experiments) and making tol bigger (the default is 1.5e-8, which is often more than needed)

stevebronder · June 4, 2026, 5:19pm

From the docs I’d also recommend trying one of the more robust solvers inside of laplace. The default solver 1 is a fast diagonal solver, while solver 3 is much more robust but a low slower.

Yingnan_Gao · June 5, 2026, 2:06am

Thanks for the replies!

I added the options following the example codes in the User manual:

functions {
  real partial_sum_2_lpmf(
    ... // other arguments
    // Laplace approximation options
    data tuple(vector, real, int, int, int, int) laplace_ops
    ){
      ... // codes for LL calculation with laplace_marginal_tol()
    }
}

data {
  // options for Laplace approximation
  real<lower=0> laplace_tol; // tolerance for optimizer
  int<lower=0> laplace_max_iter; // maximum number of steps for optimizer
  int<lower=1, upper=3> laplace_solver; // Newton solver type being used
}

transformed data{
  // For laplace approximation
  tuple(vector[M-1], real, int, int, int, int) laplace_ops = generate_laplace_options(M-1);
  laplace_ops.2 = laplace_tol;      // tolerance for optimizer
  laplace_ops.3 = laplace_max_iter; // maximum number of steps for optimizer
  laplace_ops.4 = laplace_solver;   // solver type being used
}

However, when I tried to compile the Stan program, I got errors with the message:

Compiling Stan program...
In file included from /var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:3:
In file included from stan/lib/stan_math/stan/math/mix.hpp:10:
In file included from stan/lib/stan_math/stan/math/mix/functor.hpp:11:
In file included from stan/lib/stan_math/stan/math/mix/functor/laplace_base_rng.hpp:5:
In file included from stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density.hpp:5:
stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density_estimator.hpp:157:30: error: no matching function for call to 'forward'
  157 |         value_of(std::get<0>(std::forward<Ops>(ops))),
      |                              ^~~~~~~~~~~~~~~~~
stan/lib/stan_math/stan/math/mix/prob/laplace_marginal.hpp:38:19: note: in instantiation of function template specialization 'stan::math::internal::tuple_to_laplace_options<const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int> &>' requested here
   38 |       = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
      |                   ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1401:23: note: in instantiation of function template specialization 'stan::math::laplace_marginal_tol<false, glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::ll_function_functor__, std::tuple<const int &, std::vector<int> &&, const int &, const Eigen::Matrix<double, -1, -1> &, Eigen::Block<Eigen::Matrix<double, -1, -1>, -1, 1, true> &&, const Eigen::Map<Eigen::Matrix<double, -1, -1>> &>, glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::cov_function_functor__, std::tuple<const Eigen::Matrix<double, -1, -1> &, const int &>, const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int> &>' requested here
 1401 |           stan::math::laplace_marginal_tol(ll_function_functor__(),
      |                       ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/Rtmp
mH8MbS/model-5fad2d23e8fa.hpp:531:12: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf<false, std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<double, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<double, -1, -1>>, int, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  531 |     return partial_sum_2_lpmf<propto__>(idx_y, (start + 1), (end + 1),
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/lib/stan_math/stan/math/prim/functor/reduce_sum.hpp:216:10: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf_rsfunctor__<false>::operator()<std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<double, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<double, -1, -1>>, int, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  216 |   return ReduceFunction()(std::forward<Vec>(vmapped), 0, vmapped.size() - 1,
      |          ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:2676:12: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob_impl<false, false, Eigen::Matrix<double, -1, 1>, Eigen::Matrix<int, -1, 1>, nullptr, nullptr, nullptr>' requested here
 2676 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/src/stan/model/model_base_crtp.hpp:93:50: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob<false, false, double>' requested here
   93 |     return static_cast<const M*>(this)->template log_prob<false, false, double>(
      |                                                  ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1467:3: note: in instantiation of member function 'stan::model::model_base_crtp<gl
m_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model>::log_prob' requested here
 1467 |   ~glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model()
      |   ^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__utility/forward.h:25:1: note: candidate function template not viable: 1st argument ('const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>') would lose const qualifier
   25 | forward(_LIBCPP_LIFETIMEBOUND __libcpp_remove_reference_t<_Tp>& __t) _NOEXCEPT {
      | ^                             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__utility/forward.h:31:1: note: candidate function template not viable: 1st argument ('const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>') would lose const qualifier
   31 | forward(_LIBCPP_LIFETIMEBOUND __libcpp_remove_reference_t<_Tp>&& __t) _NOEXCEPT {
      | ^                             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:3:
In file included from stan/lib/stan_math/stan/math/mix.hpp:33:
In file included from stan/lib/stan_math/stan/math/mix/prob.hpp:8:
stan/lib/stan_math/stan/math/mix/prob/laplace_marginal.hpp:38:9: error: no matching function for call to 'tuple_to_laplace_options'
   38 |       = internal::tuple_to_laplace_options(std::forward<OpsTuple>(ops));
      |         ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1401:23: note: in instantiation of function template specialization 'stan::math::laplace_marginal_tol<false, glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::ll_function_functor__, std::tuple<const int &, std::vector<int> &&, const int &, const Eigen::Matrix<stan::math::var_value<double>, -1, -1> &, Eigen::Block<Eigen::Matrix<stan::math::var_value<double>, -1, -1>, -1, 1, true> &&, const Eigen::Map<Eigen::Matrix<double, -1, -1>> &>, glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::cov_function_functor__, std::tuple<const Eigen::Matrix<stan::math::var_value<double>, -1, -1> &, const int &>, const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int> &>' requested here
 1401 |           stan::math::laplace_marginal_tol(ll_function_functor__(),
      |                       ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:531:12: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf<false, std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, in
t, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  531 |     return partial_sum_2_lpmf<propto__>(idx_y, (start + 1), (end + 1),
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/lib/stan_math/stan/math/prim/functor/reduce_sum.hpp:216:10: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf_rsfunctor__<false>::operator()<std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, int, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  216 |   return ReduceFunction()(std::forward<Vec>(vmapped), 0, vmapped.size() - 1,
      |          ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:2676:12: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob_impl<false, false, Eigen::Matrix<stan::math::var_value<double>, -1, 1>, Eigen::Matrix<int, -1, 1>, nullptr, nullptr, nullptr>' requested here
 2676 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/src/stan/model/model_base_crtp.hpp:98:50: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob<false, false, stan::math::var_value<doubl
e>>' requested here
   98 |     return static_cast<const M*>(this)->template log_prob<false, false>(theta,
      |                                                  ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1467:3: note: in instantiation of member function 'stan::model::model_base_crtp<glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model>::log_prob' requested here
 1467 |   ~glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model()
      |   ^
stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density_estimator.hpp:110:23: note: candidate template ignored: substitution failure [with Options = const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int> &]
  110 | inline constexpr auto tuple_to_laplace_options(Options&& ops) {
      |                       ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1401:11: error: no matching function for call to 'laplace_marginal_tol'
 1401 |           stan::math::laplace_marginal_tol(ll_function_functor__(),
      |           ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:531:12: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf<true, std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<double, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<double, -1, -1>>, int, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  531 |     return partial_sum_2_lpmf<propto__>(idx_y, (start + 1), (end + 1),
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/lib/stan_math/stan/math/prim/functor/reduce_sum.hpp:216:10: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf_rsfunctor__<true>::operator()<std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<double, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<double, -1, -1>>, int, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  216 |   return ReduceFunction()(std::forward<Vec>(vmapped), 0, vmapped.size() - 1,
      |          ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:2676:12: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_app
roximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob_impl<true, false, Eigen::Matrix<double, -1, 1>, Eigen::Matrix<int, -1, 1>, nullptr, nullptr, nullptr>' requested here
 2676 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/src/stan/model/model_base_crtp.hpp:115:50: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob<true, false, double>' requested here
  115 |     return static_cast<const M*>(this)->template log_prob<true, false>(theta,
      |                                                  ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1467:3: note: in instantiation of member function 'stan::model::model_base_crtp<glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model>::log_prob_propto' requested here
 1467 |   ~glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model()
      |   ^
stan/lib/stan_math/stan/math/mix/prob/laplace_marginal.hpp:32:13: 
note: candidate template ignored: substitution failure [with propto = false, LFun = ll_function_functor__, LArgs = tuple<const int &, vector<int, allocator<int>> &&, const int &, const Matrix<double, -1, -1, 0, -1, -1> &, Block<Matrix<double, -1, -1, 0, -1, -1>, -1, 1, true> &&, const Map<Matrix<double, -1, -1, 0, -1, -1>, 0, Stride<0, 0>> &>, CovarFun = cov_function_functor__, CovarArgs = tuple<const Matrix<double, -1, -1, 0, -1, -1> &, const int &>, OpsTuple = const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int> &]
   32 | inline auto laplace_marginal_tol(LFun&& L_f, LArgs&& l_args,
      |             ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1401:11: error: no matching function for call to 'laplace_marginal_tol'
 1401 |           stan::math::laplace_marginal_tol(ll_function_functor__(),
      |           ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:531:12: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf<true, std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, int, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  531 |     return partial_sum_2_lpmf<propto__>(idx_y, (start + 1), (end + 1),
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/lib/stan_math/stan/math/prim/functor/reduce_sum.hpp:216:10: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::partial_sum_2_lpmf_rsfunctor__<true>::operator()<std::vector<int>, int, unsigned long, int, std::vector<std::vector<int>>, std::vector<int>, std::vector<std::vector<double>>, std::vector<Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, Eigen::Transpose<const Eigen::Matrix<stan::math::var_value<double>, -1, -1>>, int, Eigen::Map<Eigen::Matrix<double, -1, -1>>, std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int>, nullptr>' requested here
  216 |   return ReduceFunction()(std::forward<Vec>(vmapped), 0, vmapped.size() - 1,
      |          ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:2676:12: note: in instantiation of fu
nction template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob_impl<true, false, Eigen::Matrix<stan::math::var_value<double>, -1, 1>, Eigen::Matrix<int, -1, 1>, nullptr, nullptr, nullptr>' requested here
 2676 |     return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
      |            ^
/Users/a1239028/.cmdstan/cmdstan-2.39.0/stan/src/stan/model/model_base_crtp.hpp:120:50: note: in instantiation of function template specialization 'glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model::log_prob<true, false, stan::math::var_value<double>>' requested here
  120 |     return 
static_cast<const M*>(this)->template log_prob<true, false>(theta,
      |                                                  ^
/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.hpp:1467:3: note: in instantiation of member function 'stan::model::model_base_crtp<glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model_namespace::glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model>::log_prob_propto' requested here
 1467 |   ~glm_logistic_normal_multinomial_isometric_single_laplace_approximation_model()
      |   ^
stan/lib/stan_math/stan/math/mix/prob/laplace_marginal.hpp:32:13: note: candidate template ignored: substitution failure [with propto = false, LFun = ll_function_functor__, LArgs = tuple<const int &, vector<int, allocator<int>> &&, const int &, const Matrix<var_value<double, void>, -1, -1, 0, -1, -1> &, Block<Matrix<var_value<double, void>, -1, -1, 0, -1, -1>, -1, 1, true> &&, const Map<Matrix<double, -1, -1, 0, -1, -1>, 0, Stride<0, 0>> &>, CovarFun = cov_function_functor__, CovarArgs = tuple<const Matrix<var_value<double, void>, -1, -1, 0, -1, -1> &, const int &>, OpsTuple = const std::tuple<Eigen::Matrix<double, -1, 1>, double, int, int, int, int> &]
   32 | inline auto laplace_marginal_tol(LFun&& L_f, LArgs&& l_args,
      |             ^
4 errors generated.
make: *** [/var/folders/11/bb8zxz8j0fx1lv3_s0mbhkr00000gp/T/RtmpmH8MbS/model-5fad2d23e8fa.o] Error 1

Scanning through it, it seems to be a tuple signature mismatch with stan_math function expectations, but I cannot figure out where I did things wrongly.

Also, the documentation for control parameters have a weird chunk of code that reads:

tuple(vector[theta_size], real, int, int, int, int, int) laplace_ops = generate_laplace_options(theta_size);
// Modify solver type
laplace_ops.5 = 2;
// Turn off fallthrough
laplace_ops.7 = 0;

There appears to be one extra int element in the returned tuple that differs from the signature in all other variants in the section.

Yingnan_Gao · June 5, 2026, 6:09am

Some update here:
I cannot locate the problem so I feed the error message into ChatGPT and tried a few work arounds. Moving laplace_marginal_tol to the model block from the custom partial sum function does not help, and the same error persists.

Eventually, ChatGPT suggest that I manually edit the file cmdstan-2.39.0/stan/lib/stan_math/stan/math/mix/functor/laplace_marginal_density_estimator.hpp with the substitution of replacing:

std::forward<Ops>(ops)

with

std::forward<Options>(ops)

There are two occurrence of this substitution, one in line 157 and one in line 166.

ChatGPT’s remark on the root issue is something like:

I now think this is an actual Stan Math forwarding bug inside tuple_to_laplace_options().

Your dynamically supplied tuple exposed it because of const-reference binding semantics in generated C++.

I have very limited knowledge on modern C++ coding, so I cannot assess if ChatGPT’s remark is accurate. But after I substitute the two lines as described above, I can compile the program without error (with laplace_marginal_tol inside model block). I still need to verify if putting it back to the custom partial sum function works, but the progress so far seems promising.

Tagging previous respondents: @stevebronder @avehtari

Maybe @Bob_Carpenter will also be interested in this issue?

stevebronder · June 5, 2026, 6:51am

That is a real bug! Thank you for catching that I’ll put in a fix for this tmrw

stevebronder · June 5, 2026, 4:13pm

github.com/stan-dev/math

Fix typo in tuple_to_laplace_options forwarding function (#3330)

develop ← fix/tuple_to_laplace_options-forwarding

opened 04:12PM - 05 Jun 26 UTC

SteveBronder

+2 -2

## Summary Fixes `tuple_to_laplace_options` typo so that we use `std::forward…<Options>` instead of `std::forward<Opts>` ## Tests None for one liner ## Checklist - [x] Copyright holder: Simons Foundation The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses: - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause) - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/) - [x] the basic tests are passing - unit tests pass (to run, use: `./runTests.py test/unit`) - header checks pass, (`make test-headers`) - dependencies checks pass, (`make test-math-dependencies`) - docs build, (`make doxygen`) - code passes the built in [C++ standards](https://github.com/stan-dev/stan/wiki/Coding-Style-and-Idioms) checks (`make cpplint`) - [x] the code is written in idiomatic C++ and changes are documented in the doxygen - [x] the new changes are tested

Bob_Carpenter · June 5, 2026, 6:30pm

@stevebronder and @avehtari know much more about this than me!

Topic		Replies	Views
Integrated Laplace approximation roadmap Developers	49	5530	January 13, 2025
CmdStan & Stan 2.37 release candidate General	20	1508	October 7, 2025
Algebraic sovler problems (Laplace approximation in Stan) Modeling	65	7130	April 13, 2018
Embedded Laplace - model.hpp gets C++ compiler novel Developers	5	214	June 13, 2025
Release of CmdStan 2.39.0 Announcements	15	310	June 5, 2026

Using Laplace approximation for logistic normal multinomial model

Related topics