From 1f33f4a93d29c6547dfb8ed1a08c5edf120c5329 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 19 Oct 2026 08:56:01 +0300 Subject: [PATCH 1/2] closed-form derivatives for the built-in embedded Laplace likelihoods --- stan/math/mix/functor/laplace_likelihood.hpp | 122 +++++++- .../laplace_marginal_bernoulli_logit_lpmf.hpp | 66 +++++ ...place_marginal_neg_binomial_2_log_lpmf.hpp | 154 ++++++++++ .../laplace_marginal_poisson_log_lpmf.hpp | 63 ++++ .../laplace_analytic_derivatives_test.cpp | 268 ++++++++++++++++++ 5 files changed, 669 insertions(+), 4 deletions(-) create mode 100644 test/unit/math/laplace/laplace_analytic_derivatives_test.cpp diff --git a/stan/math/mix/functor/laplace_likelihood.hpp b/stan/math/mix/functor/laplace_likelihood.hpp index ec5be6b36d1..75b7d23f468 100644 --- a/stan/math/mix/functor/laplace_likelihood.hpp +++ b/stan/math/mix/functor/laplace_likelihood.hpp @@ -10,6 +10,7 @@ #include #include +#include #include namespace stan { @@ -21,7 +22,46 @@ namespace math { */ namespace laplace_likelihood { +/** + * Specialize this trait with `value = true` for a likelihood functor `F` if + * the Hessian of its log likelihood with respect to `theta` is diagonal and + * its derivatives with respect to `theta` are known in closed form. The + * specialization must have three static member functions, which take + * `theta` and the arguments of `F::operator()` without the message stream, + * all with arithmetic scalar types, and return an `Eigen::VectorXd`: + * - `theta_grad(theta, args...)`: the gradient; + * - `neg_hessian_diag(theta, args...)`: the diagonal of the negative + * Hessian; + * - `third_diff(theta, args...)`: the diagonal of the third derivative. + * When all arguments are arithmetic, `theta_grad()`, `diagonal_hessian()`, + * `block_hessian()` and `third_diff()` below then use these functions + * instead of nested autodiff. + * @tparam F Type of the likelihood functor + */ +template +struct analytic_derivatives : std::false_type {}; + namespace internal { + +/** + * True if the closed-form derivatives of `F` (see `analytic_derivatives`) + * are used for these argument types: `F` has them, and `theta` and all + * arguments have arithmetic scalar types. + */ +template +inline constexpr bool use_analytic_derivatives_v + = analytic_derivatives>::value&& + std::is_arithmetic>>::value&& + std::conjunction_v< + std::is_arithmetic>>...>; + +template +using require_analytic_derivatives_t = require_t< + std::bool_constant>>; + +template +using require_not_analytic_derivatives_t = require_t< + std::bool_constant>>; /** * @tparam F A functor with `opertor()(Args&&...)` returning a scalar * @tparam Theta A class assignable to an Eigen vector type @@ -53,7 +93,8 @@ inline auto log_likelihood(F&& f, Theta&& theta, Stream* msgs, Args&&... args) { * @param args Variadic arguments for the likelihood function. */ template * = nullptr> + require_eigen_vector_vt* = nullptr, + require_not_analytic_derivatives_t* = nullptr> inline auto theta_grad(F&& f, Theta&& theta, Stream* msgs, Args&&... args) { using Eigen::Dynamic; using Eigen::Matrix; @@ -104,7 +145,8 @@ inline void ll_arg_grad(F&& f, Theta&& theta, Stream* msgs, Args&&... args) { * @param args Variadic arguments for the likelihood function. */ template * = nullptr> + require_eigen_vector_vt* = nullptr, + require_not_analytic_derivatives_t* = nullptr> inline auto diagonal_hessian(F&& f, Theta&& theta, Stream* msgs, Args&&... args) { using Eigen::Dynamic; @@ -134,7 +176,8 @@ inline auto diagonal_hessian(F&& f, Theta&& theta, Stream* msgs, * @param args Variadic arguments for the likelihood function. */ template * = nullptr> + require_eigen_vector_vt* = nullptr, + require_not_analytic_derivatives_t* = nullptr> inline auto block_hessian(F&& f, Theta&& theta, const Eigen::Index hessian_block_size, Stream* msgs, Args&&... args) { @@ -224,7 +267,8 @@ inline auto diff(F&& f, Theta&& theta, const Eigen::Index hessian_block_size, * @param args Variadic arguments for likelihood function. */ template * = nullptr> + require_eigen_vector_t* = nullptr, + require_not_analytic_derivatives_t* = nullptr> inline Eigen::VectorXd third_diff(F&& f, Theta&& theta, Stream&& msgs, Args&&... args) { nested_rev_autodiff nested; @@ -355,6 +399,76 @@ inline auto diff_eta_implicit(F&& f, V_t&& v, Theta&& theta, Stream* msgs, grad(f_sum.d_.vi_); } +/** + * Gradient of `f` with respect to `theta`, from the closed-form derivatives + * of `F` (see `analytic_derivatives`). + */ +template * = nullptr, + require_analytic_derivatives_t* = nullptr> +inline Eigen::VectorXd theta_grad(F&& f, Theta&& theta, Stream* msgs, + Args&&... args) { + return analytic_derivatives>::theta_grad(theta, args...); +} + +/** + * Diagonal of the negative Hessian of `f` with respect to `theta`, from the + * closed-form derivatives of `F` (see `analytic_derivatives`). + */ +template * = nullptr, + require_analytic_derivatives_t* = nullptr> +inline Eigen::VectorXd diagonal_hessian(F&& f, Theta&& theta, Stream* msgs, + Args&&... args) { + return analytic_derivatives>::neg_hessian_diag(theta, + args...); +} + +/** + * Negative Hessian of `f` with respect to `theta` as a sparse block-diagonal + * matrix, from the closed-form derivatives of `F` (see + * `analytic_derivatives`). The Hessian of such an `F` is diagonal; the + * matrix has the sparsity pattern of the autodiff version + * (`hessian_block_diag()`): full blocks of size `hessian_block_size`, with + * zero off-diagonal entries. + */ +template * = nullptr, + require_analytic_derivatives_t* = nullptr> +inline Eigen::SparseMatrix block_hessian( + F&& f, Theta&& theta, const Eigen::Index hessian_block_size, Stream* msgs, + Args&&... args) { + const Eigen::VectorXd neg_hessian_diag + = analytic_derivatives>::neg_hessian_diag(theta, args...); + const Eigen::Index theta_size = theta.size(); + Eigen::SparseMatrix neg_hessian(theta_size, theta_size); + neg_hessian.reserve( + Eigen::VectorXi::Constant(theta_size, hessian_block_size)); + const Eigen::Index n_blocks = theta_size / hessian_block_size; + for (Eigen::Index j = 0; j < n_blocks; ++j) { + const Eigen::Index start = j * hessian_block_size; + for (Eigen::Index c = start; c < start + hessian_block_size; ++c) { + for (Eigen::Index r = start; r < start + hessian_block_size; ++r) { + neg_hessian.insert(r, c) = r == c ? neg_hessian_diag(r) : 0.0; + } + } + } + neg_hessian.makeCompressed(); + return neg_hessian; +} + +/** + * Diagonal of the third derivative of `f` with respect to `theta`, from the + * closed-form derivatives of `F` (see `analytic_derivatives`). + */ +template * = nullptr, + require_analytic_derivatives_t* = nullptr> +inline Eigen::VectorXd third_diff(F&& f, Theta&& theta, Stream&& msgs, + Args&&... args) { + return analytic_derivatives>::third_diff(theta, args...); +} + } // namespace internal /** 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 79e1598486d..78e809988c0 100644 --- a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp +++ b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp @@ -3,6 +3,8 @@ #include #include +#include +#include #include #include #include @@ -58,6 +60,70 @@ struct bernoulli_logit_likelihood { } }; +namespace laplace_likelihood { +/** + * Closed-form derivatives of `bernoulli_logit_likelihood` with respect to + * `theta` (see `analytic_derivatives`). With `p_g = inv_logit(theta_g + + * mean_g)` and `q_g = 1 - p_g = inv_logit(-theta_g - mean_g)` (computed + * without cancellation), the gradient is `counts_g - n_g p_g`, the negative + * Hessian is `n_g p_g q_g`, and the third derivative is + * `-n_g p_g q_g (q_g - p_g)`. + */ +template <> +struct analytic_derivatives : std::true_type { + // sum of the counts and number of observations of each group + template + static void group_sums(Eigen::Index theta_size, const YVec& y, + const std::vector& y_index, + Eigen::VectorXd& counts, Eigen::VectorXd& n) { + counts = Eigen::VectorXd::Zero(theta_size); + n = Eigen::VectorXd::Zero(theta_size); + for (size_t i = 0; i < y_index.size(); i++) { + counts(y_index[i] - 1) += y[i]; + n(y_index[i] - 1) += 1; + } + } + + template + static void probabilities(const Theta& theta, const YVec& y, + const std::vector& y_index, const Mean& mean, + Eigen::VectorXd& counts, Eigen::VectorXd& n, + Eigen::VectorXd& p, Eigen::VectorXd& q) { + group_sums(theta.size(), y, y_index, counts, n); + Eigen::VectorXd theta_offset = add(theta, mean); + p = inv_logit(theta_offset); + q = inv_logit(-theta_offset); + } + + template + static Eigen::VectorXd theta_grad(const Theta& theta, const YVec& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, n, p, q; + probabilities(theta, y, y_index, mean, counts, n, p, q); + return counts - n.cwiseProduct(p); + } + + template + static Eigen::VectorXd neg_hessian_diag(const Theta& theta, const YVec& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, n, p, q; + probabilities(theta, y, y_index, mean, counts, n, p, q); + return n.cwiseProduct(p).cwiseProduct(q); + } + + template + static Eigen::VectorXd third_diff(const Theta& theta, const YVec& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, n, p, q; + probabilities(theta, y, y_index, mean, counts, n, p, q); + return -n.cwiseProduct(p).cwiseProduct(q).cwiseProduct(q - p); + } +}; +} // namespace laplace_likelihood + /** * Wrapper function around the laplace_marginal function for * a logistic Bernoulli likelihood. Returns the marginal density 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 699a54b914a..bc56d3444e3 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 @@ -2,6 +2,9 @@ #define STAN_MATH_MIX_PROB_LAPLACE_MARGINAL_NEG_BINOMIAL_2_LOG_LPMF_HPP #include +#include +#include +#include #include #include @@ -193,6 +196,157 @@ struct neg_binomial_2_log_likelihood_summary { } }; +namespace laplace_likelihood { +/** + * Closed-form derivatives of `neg_binomial_2_log_likelihood` with respect + * to `theta` (see `analytic_derivatives`). With the dispersion `eta_g` of + * group `g` (a scalar `eta` or one entry per group), + * `s_g = counts_g + eta_g n_g`, `p_g = inv_logit(theta_g + mean_g - + * log(eta_g))` and `q_g = 1 - p_g` (computed without cancellation), the + * gradient is `counts_g - s_g p_g`, the negative Hessian is `s_g p_g q_g`, + * and the third derivative is `-s_g p_g q_g (q_g - p_g)`. + */ +template <> +struct analytic_derivatives : std::true_type { + // gradient, negative Hessian and third derivative from the sums of the + // counts `counts`, `s = counts + eta n`, and x = theta + mean - log(eta) + static Eigen::VectorXd grad(const Eigen::VectorXd& counts, + const Eigen::VectorXd& s, + const Eigen::VectorXd& x) { + return counts - s.cwiseProduct(inv_logit(x)); + } + static Eigen::VectorXd neg_hess(const Eigen::VectorXd& s, + const Eigen::VectorXd& x) { + return s.cwiseProduct(inv_logit(x)).cwiseProduct(inv_logit(-x)); + } + static Eigen::VectorXd third(const Eigen::VectorXd& s, + const Eigen::VectorXd& x) { + const Eigen::VectorXd p = inv_logit(x); + const Eigen::VectorXd q = inv_logit(-x); + return -s.cwiseProduct(p).cwiseProduct(q).cwiseProduct(q - p); + } + + template + static void parts(const Theta& theta, const Eta& eta, + const std::vector& y, const std::vector& y_index, + const Mean& mean, Eigen::VectorXd& counts, + Eigen::VectorXd& s, Eigen::VectorXd& x) { + const Eigen::Index theta_size = theta.size(); + counts = Eigen::VectorXd::Zero(theta_size); + Eigen::VectorXd n = Eigen::VectorXd::Zero(theta_size); + for (size_t i = 0; i < y.size(); i++) { + counts(y_index[i] - 1) += y[i]; + n(y_index[i] - 1) += 1; + } + Eigen::VectorXd eta_g(theta_size); + if constexpr (is_stan_scalar::value) { + eta_g.setConstant(eta); + } else { + check_size_match("neg_binomial_2_log_likelihood", "eta", eta.size(), + "theta", theta_size); + for (Eigen::Index g = 0; g < theta_size; ++g) { + eta_g(g) = eta[g]; + } + } + s = counts + eta_g.cwiseProduct(n); + x = add(theta, mean) - eta_g.array().log().matrix(); + } + + template + static Eigen::VectorXd theta_grad(const Theta& theta, const Eta& eta, + const std::vector& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, s, x; + parts(theta, eta, y, y_index, mean, counts, s, x); + return grad(counts, s, x); + } + + template + static Eigen::VectorXd neg_hessian_diag(const Theta& theta, const Eta& eta, + const std::vector& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, s, x; + parts(theta, eta, y, y_index, mean, counts, s, x); + return neg_hess(s, x); + } + + template + static Eigen::VectorXd third_diff(const Theta& theta, const Eta& eta, + const std::vector& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, s, x; + parts(theta, eta, y, y_index, mean, counts, s, x); + return third(s, x); + } +}; + +/** + * Closed-form derivatives of `neg_binomial_2_log_likelihood_summary` with + * respect to `theta`; the same formulas as for + * `neg_binomial_2_log_likelihood`, with a scalar `eta`. + */ +template <> +struct analytic_derivatives + : std::true_type { + template + static void parts(const Theta& theta, const Eta& eta, + const std::vector& n_per_group, + const std::vector& counts_per_group, const Mean& mean, + Eigen::VectorXd& counts, Eigen::VectorXd& s, + Eigen::VectorXd& x) { + const Eigen::Index theta_size = theta.size(); + counts = Eigen::Map(counts_per_group.data(), + counts_per_group.size()) + .cast(); + const Eigen::VectorXd n = Eigen::Map( + n_per_group.data(), n_per_group.size()) + .cast(); + check_size_match("neg_binomial_2_log_likelihood_summary", + "counts_per_group", counts.size(), "theta", theta_size); + check_size_match("neg_binomial_2_log_likelihood_summary", "n_per_group", + n.size(), "theta", theta_size); + s = counts + eta * n; + x = add(theta, mean) - Eigen::VectorXd::Constant(theta_size, std::log(eta)); + } + + template + static Eigen::VectorXd theta_grad(const Theta& theta, const Eta& eta, + const std::vector& y, + const std::vector& n_per_group, + const std::vector& counts_per_group, + const Mean& mean) { + Eigen::VectorXd counts, s, x; + parts(theta, eta, n_per_group, counts_per_group, mean, counts, s, x); + return analytic_derivatives::grad(counts, s, + x); + } + + template + static Eigen::VectorXd neg_hessian_diag( + const Theta& theta, const Eta& eta, const std::vector& y, + const std::vector& n_per_group, + const std::vector& counts_per_group, const Mean& mean) { + Eigen::VectorXd counts, s, x; + parts(theta, eta, n_per_group, counts_per_group, mean, counts, s, x); + return analytic_derivatives::neg_hess(s, x); + } + + template + static Eigen::VectorXd third_diff(const Theta& theta, const Eta& eta, + const std::vector& y, + const std::vector& n_per_group, + const std::vector& counts_per_group, + const Mean& mean) { + Eigen::VectorXd counts, s, x; + parts(theta, eta, n_per_group, counts_per_group, mean, counts, s, x); + return analytic_derivatives::third(s, x); + } +}; +} // namespace laplace_likelihood + /** * Wrapper function around the laplace_marginal function for * a negative binomial likelihood. Uses the 2nd parameterization. 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 c9cb99fa629..7bc279a9291 100644 --- a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp +++ b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp @@ -2,6 +2,7 @@ #define STAN_MATH_MIX_PROB_LAPLACE_MARGINAL_POISSON_LOG_LPMF_HPP #include +#include #include #include #include @@ -53,6 +54,68 @@ struct poisson_log_likelihood { } }; +namespace laplace_likelihood { +/** + * Closed-form derivatives of `poisson_log_likelihood` with respect to + * `theta` (see `analytic_derivatives`). With `m_g = n_g exp(theta_g + + * mean_g)`, where `n_g` is the number of observations of group `g`, the + * gradient is `counts_g - m_g`, the negative Hessian is `m_g`, and the third + * derivative is `-m_g`. + */ +template <> +struct analytic_derivatives : std::true_type { + // sum of the counts and number of observations of each group + template + static void group_sums(Eigen::Index theta_size, const YVec& y, + const std::vector& y_index, + Eigen::VectorXd& counts, Eigen::VectorXd& n) { + counts = Eigen::VectorXd::Zero(theta_size); + n = Eigen::VectorXd::Zero(theta_size); + for (size_t i = 0; i < y_index.size(); i++) { + counts(y_index[i] - 1) += y[i]; + n(y_index[i] - 1) += 1; + } + } + + template + static void rates(const Theta& theta, const YVec& y, + const std::vector& y_index, const Mean& mean, + Eigen::VectorXd& counts, Eigen::VectorXd& m) { + Eigen::VectorXd n; + group_sums(theta.size(), y, y_index, counts, n); + Eigen::VectorXd theta_offset = add(theta, mean); + m = n.cwiseProduct(theta_offset.array().exp().matrix()); + } + + template + static Eigen::VectorXd theta_grad(const Theta& theta, const YVec& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, m; + rates(theta, y, y_index, mean, counts, m); + return counts - m; + } + + template + static Eigen::VectorXd neg_hessian_diag(const Theta& theta, const YVec& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, m; + rates(theta, y, y_index, mean, counts, m); + return m; + } + + template + static Eigen::VectorXd third_diff(const Theta& theta, const YVec& y, + const std::vector& y_index, + const Mean& mean) { + Eigen::VectorXd counts, m; + rates(theta, y, y_index, mean, counts, m); + return -m; + } +}; +} // namespace laplace_likelihood + /** * Wrapper function around the laplace_marginal function for * a log poisson likelihood. Returns the marginal density diff --git a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp new file mode 100644 index 00000000000..5ed15dcca23 --- /dev/null +++ b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp @@ -0,0 +1,268 @@ +#include +#include +#include +#include +#include +#include +#include + +namespace { + +using stan::math::bernoulli_logit_likelihood; +using stan::math::neg_binomial_2_log_likelihood; +using stan::math::neg_binomial_2_log_likelihood_summary; +using stan::math::poisson_log_likelihood; +using stan::math::laplace_likelihood::analytic_derivatives; +namespace lli = stan::math::laplace_likelihood::internal; + +// The likelihood `LL` behind a lambda: no `analytic_derivatives` +// specialization, so the derivative functions use nested autodiff. +template +auto autodiff_version() { + return [](auto&&... args) { return LL{}(args...); }; +} + +void expect_near_rel(const Eigen::VectorXd& expected, + const Eigen::VectorXd& actual, double tol, + const char* what) { + ASSERT_EQ(expected.size(), actual.size()) << what; + for (Eigen::Index i = 0; i < expected.size(); ++i) { + EXPECT_TRUE(std::isfinite(actual(i))) << what << " [" << i << "]"; + EXPECT_NEAR(expected(i), actual(i), + tol * std::fmax(1.0, std::fabs(expected(i)))) + << what << " [" << i << "]"; + } +} + +void expect_same_sparse(const Eigen::SparseMatrix& expected, + const Eigen::SparseMatrix& actual, double tol, + const char* what) { + ASSERT_EQ(expected.rows(), actual.rows()) << what; + ASSERT_EQ(expected.cols(), actual.cols()) << what; + EXPECT_EQ(expected.nonZeros(), actual.nonZeros()) << what; + const Eigen::MatrixXd e = expected; + const Eigen::MatrixXd a = actual; + for (Eigen::Index c = 0; c < e.cols(); ++c) { + for (Eigen::Index r = 0; r < e.rows(); ++r) { + EXPECT_NEAR(e(r, c), a(r, c), tol * std::fmax(1.0, std::fabs(e(r, c)))) + << what << " (" << r << ", " << c << ")"; + } + } +} + +// Each derivative function of laplace_likelihood::internal, for the +// analytic version (`LL{}`) and for the autodiff version, with +// hessian_block_size 1 and 2. +template +void expect_analytic_equals_autodiff(const Eigen::VectorXd& theta, double tol, + const Args&... args) { + static_assert(lli::use_analytic_derivatives_v, + "the analytic derivatives must be selected"); + auto ad = autodiff_version(); + std::ostream* msgs = nullptr; + expect_near_rel(lli::theta_grad(ad, theta, msgs, args...), + lli::theta_grad(LL{}, theta, msgs, args...), tol, + "theta_grad"); + expect_near_rel(lli::diagonal_hessian(ad, theta, msgs, args...), + lli::diagonal_hessian(LL{}, theta, msgs, args...), tol, + "diagonal_hessian"); + for (int block_size : {1, 2}) { + expect_same_sparse( + lli::block_hessian(ad, theta, block_size, msgs, args...), + lli::block_hessian(LL{}, theta, block_size, msgs, args...), tol, + "block_hessian"); + } + expect_near_rel(lli::third_diff(ad, theta, msgs, args...), + lli::third_diff(LL{}, theta, msgs, args...), tol, + "third_diff"); +} + +// 4 groups, 1 to 3 observations per group +const std::vector y_counts{0, 3, 7, 2, 1, 4, 12, 9}; +const std::vector y_binary{0, 1, 1, 0, 1, 1, 0, 1}; +const std::vector y_index{1, 1, 1, 2, 3, 3, 4, 4}; + +// a likelihood with the formulas of poisson_log_likelihood that counts the +// calls of its analytic derivatives +struct counted_poisson_likelihood { + template + auto operator()(Args&&... args) const { + return poisson_log_likelihood{}(std::forward(args)...); + } +}; +int n_analytic_calls = 0; + +} // namespace + +namespace stan { +namespace math { +namespace laplace_likelihood { +template <> +struct analytic_derivatives : std::true_type { + template + static Eigen::VectorXd theta_grad(const Args&... args) { + ++n_analytic_calls; + return analytic_derivatives::theta_grad(args...); + } + template + static Eigen::VectorXd neg_hessian_diag(const Args&... args) { + ++n_analytic_calls; + return analytic_derivatives::neg_hessian_diag( + args...); + } + template + static Eigen::VectorXd third_diff(const Args&... args) { + ++n_analytic_calls; + return analytic_derivatives::third_diff(args...); + } +}; +} // namespace laplace_likelihood +} // namespace math +} // namespace stan + +namespace { + +struct se_covariance { + template + auto operator()(const std::vector& x, const T1& sigma, const T2& rho, + std::ostream*) const { + return stan::math::add_diag(stan::math::gp_exp_quad_cov(x, sigma, rho), + 1e-8); + } +}; +const std::vector x_gp{0.0, 1.3, 2.1, 3.7}; + +} // namespace + +TEST(laplace_analytic_derivatives, trait) { + EXPECT_TRUE(analytic_derivatives::value); + EXPECT_TRUE(analytic_derivatives::value); + EXPECT_TRUE(analytic_derivatives::value); + EXPECT_TRUE( + analytic_derivatives::value); + auto ad = autodiff_version(); + EXPECT_FALSE(analytic_derivatives::value); + // a var argument selects the autodiff version + EXPECT_FALSE((lli::use_analytic_derivatives_v< + neg_binomial_2_log_likelihood, Eigen::VectorXd, stan::math::var, + std::vector, std::vector, Eigen::VectorXd>)); +} + +TEST(laplace_analytic_derivatives, used_by_laplace_marginal) { + const Eigen::VectorXd mean = Eigen::VectorXd::Constant(4, 0.5); + std::ostream* msgs = nullptr; + n_analytic_calls = 0; + double lp = stan::math::laplace_marginal( + counted_poisson_likelihood{}, + std::forward_as_tuple(y_counts, y_index, mean), 1, se_covariance{}, + std::forward_as_tuple(x_gp, 1.0, 1.5), msgs); + EXPECT_TRUE(std::isfinite(lp)); + EXPECT_GT(n_analytic_calls, 0); +} + +TEST(laplace_analytic_derivatives, equal_to_autodiff) { + Eigen::VectorXd theta(4); + theta << 0.3, -1.2, 1.7, 0.05; + const Eigen::VectorXd mean = Eigen::VectorXd::Constant(4, 0.4); + const double tol = 1e-12; + for (double shift : {0.0, -30.0, 30.0}) { + const Eigen::VectorXd th = theta.array() + shift; + expect_analytic_equals_autodiff(th, tol, y_counts, + y_index, mean); + expect_analytic_equals_autodiff( + th, tol, y_binary, y_index, mean); + expect_analytic_equals_autodiff( + th, tol, 2.5, y_counts, y_index, mean); + Eigen::VectorXd eta(4); + eta << 0.7, 2.5, 11.0, 1.3; + expect_analytic_equals_autodiff( + th, tol, eta, y_counts, y_index, mean); + const std::vector n_per_group{3, 1, 2, 2}; + const std::vector counts_per_group{10, 2, 5, 21}; + expect_analytic_equals_autodiff( + th, tol, 2.5, y_counts, n_per_group, counts_per_group, mean); + } +} + +TEST(laplace_analytic_derivatives, finite_at_large_theta) { + // exp(theta) overflows here; the inv_logit forms stay finite + Eigen::VectorXd theta(4); + theta << -800, 800, -750, 760; + const Eigen::VectorXd mean = Eigen::VectorXd::Zero(4); + expect_analytic_equals_autodiff( + theta, 1e-12, y_binary, y_index, mean); + expect_analytic_equals_autodiff( + theta, 1e-12, 2.5, y_counts, y_index, mean); +} + +TEST(laplace_analytic_derivatives, end_to_end_equal_to_autodiff) { + using stan::math::var; + const Eigen::VectorXd mean = Eigen::VectorXd::Constant(4, 0.5); + std::ostream* msgs = nullptr; + const double tol = 1e-7; + for (int solver : {1, 2, 3}) { + for (int block_size : {1, 2}) { + auto ops = std::make_tuple(Eigen::VectorXd::Zero(4).eval(), 1e-10, 200, + solver, 1000, false); + // value and gradient w.r.t. sigma, rho (and eta) of `lpmf` + auto run = [&](auto&& lpmf) { + var sigma = 1.0, rho = 1.5, eta = 2.5; + var lp = lpmf(sigma, rho, eta); + stan::math::grad(lp.vi_); + std::vector out{lp.val(), sigma.adj(), rho.adj(), eta.adj()}; + stan::math::recover_memory(); + return out; + }; + auto compare = [&](const std::vector& a, + const std::vector& b, const char* what) { + for (size_t i = 0; i < a.size(); ++i) { + EXPECT_NEAR(a[i], b[i], tol * std::fmax(1.0, std::fabs(a[i]))) + << what << " solver " << solver << " block size " << block_size + << " output " << i; + } + }; + compare(run([&](auto s, auto r, auto) { + return stan::math::laplace_marginal_tol( + autodiff_version(), + std::forward_as_tuple(y_counts, y_index, mean), block_size, + se_covariance{}, std::forward_as_tuple(x_gp, s, r), ops, + msgs); + }), + run([&](auto s, auto r, auto) { + return stan::math::laplace_marginal_tol_poisson_log_lpmf( + y_counts, y_index, mean, block_size, se_covariance{}, + std::forward_as_tuple(x_gp, s, r), ops, msgs); + }), + "poisson"); + compare(run([&](auto s, auto r, auto) { + return stan::math::laplace_marginal_tol( + autodiff_version(), + std::forward_as_tuple(y_binary, y_index, mean), block_size, + se_covariance{}, std::forward_as_tuple(x_gp, s, r), ops, + msgs); + }), + run([&](auto s, auto r, auto) { + return stan::math::laplace_marginal_tol_bernoulli_logit_lpmf( + y_binary, y_index, mean, block_size, se_covariance{}, + std::forward_as_tuple(x_gp, s, r), ops, msgs); + }), + "bernoulli"); + // eta as a var: the Newton iterations use the analytic derivatives, + // the gradient w.r.t. eta the autodiff ones + compare(run([&](auto s, auto r, auto e) { + return stan::math::laplace_marginal_tol( + autodiff_version(), + std::forward_as_tuple(e, y_counts, y_index, mean), + block_size, se_covariance{}, + std::forward_as_tuple(x_gp, s, r), ops, msgs); + }), + run([&](auto s, auto r, auto e) { + return stan::math::laplace_marginal_tol_neg_binomial_2_log_lpmf( + y_counts, y_index, e, mean, block_size, se_covariance{}, + std::forward_as_tuple(x_gp, s, r), ops, msgs); + }), + "neg_binomial_2 (eta var)"); + } + } +} From d7d750235b4ff79ef996a714b56ee2e8769bcc92 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 19 Oct 2026 16:54:49 +0300 Subject: [PATCH 2/2] fix cpplint whitespace/braces in the analytic-derivatives test --- test/unit/math/laplace/laplace_analytic_derivatives_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp index 5ed15dcca23..1e7c81558aa 100644 --- a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp +++ b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp @@ -88,7 +88,7 @@ const std::vector y_index{1, 1, 1, 2, 3, 3, 4, 4}; struct counted_poisson_likelihood { template auto operator()(Args&&... args) const { - return poisson_log_likelihood{}(std::forward(args)...); + return poisson_log_likelihood()(std::forward(args)...); } }; int n_analytic_calls = 0;