From 383a0b0659db3eb211dd22624c86c81e1cadc946 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Mon, 12 Oct 2026 11:25:47 +0200 Subject: [PATCH 1/3] Tests: ordered_logistic(_glm)_lpmf derivatives against closed forms --- .../prob/ordered_logistic_glm_lpmf_test.cpp | 10 + .../math/mix/prob/ordered_logistic_test.cpp | 22 ++ .../prob/ordered_logistic_test_helpers.hpp | 230 ++++++++++++++++++ 3 files changed, 262 insertions(+) create mode 100644 test/unit/math/mix/prob/ordered_logistic_test_helpers.hpp diff --git a/test/unit/math/mix/prob/ordered_logistic_glm_lpmf_test.cpp b/test/unit/math/mix/prob/ordered_logistic_glm_lpmf_test.cpp index 649d5cf2f46..f9c95f0df04 100644 --- a/test/unit/math/mix/prob/ordered_logistic_glm_lpmf_test.cpp +++ b/test/unit/math/mix/prob/ordered_logistic_glm_lpmf_test.cpp @@ -1,5 +1,6 @@ #include #include +#include TEST_F(AgradRev, mathMixScalFun_ordered_logistic_glm_lpmf) { auto f = [](const auto y) { @@ -85,3 +86,12 @@ TEST_F(AgradRev, << "cuts_ffv(" << j << ").val_.val_.adj() non-finite"; } } + +// Against closed forms at and near the cuts, in the tails and for close cuts. +TEST_F(AgradRev, mathMixScalFun_ordered_logistic_glm_lpmf_derivatives) { + auto lpmf + = [](const auto& y, const auto& x, const auto& beta, const auto& cuts) { + return stan::math::ordered_logistic_glm_lpmf(y, x, beta, cuts); + }; + ordered_logistic_test::expect_all(lpmf, true); +} diff --git a/test/unit/math/mix/prob/ordered_logistic_test.cpp b/test/unit/math/mix/prob/ordered_logistic_test.cpp index 28d4e5856ff..cd3a1c67fc7 100644 --- a/test/unit/math/mix/prob/ordered_logistic_test.cpp +++ b/test/unit/math/mix/prob/ordered_logistic_test.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -722,3 +723,24 @@ TEST_F(AgradRev, ProbDistributionsOrdLog_top_category_value_and_gradient) { EXPECT_FLOAT_EQ(c[0].adj(), 0.0); EXPECT_FLOAT_EQ(c[1].adj(), -0.574442516811659); } + +// Against closed forms at and near the cuts, in the tails and for close cuts. +TEST_F(AgradRev, ProbDistributionsOrdLog_derivatives) { + auto lpmf + = [](const auto& y, const auto& x, const auto& beta, const auto& cuts) { + return stan::math::ordered_logistic_lpmf( + y, stan::math::multiply(x, beta), cuts); + }; + ordered_logistic_test::expect_all(lpmf, false); +} + +// The same with one cut vector per location. +TEST_F(AgradRev, ProbDistributionsOrdLog_derivatives_cuts_array) { + auto lpmf + = [](const auto& y, const auto& x, const auto& beta, const auto& cuts) { + std::vector> cuts_arr(x.rows(), cuts); + return stan::math::ordered_logistic_lpmf( + y, stan::math::multiply(x, beta), cuts_arr); + }; + ordered_logistic_test::expect_all(lpmf, false); +} diff --git a/test/unit/math/mix/prob/ordered_logistic_test_helpers.hpp b/test/unit/math/mix/prob/ordered_logistic_test_helpers.hpp new file mode 100644 index 00000000000..a67510fe278 --- /dev/null +++ b/test/unit/math/mix/prob/ordered_logistic_test_helpers.hpp @@ -0,0 +1,230 @@ +#ifndef TEST_UNIT_MATH_MIX_PROB_ORDERED_LOGISTIC_TEST_HELPERS_HPP +#define TEST_UNIT_MATH_MIX_PROB_ORDERED_LOGISTIC_TEST_HELPERS_HPP + +#include +#include +#include +#include +#include + +/** + * Closed-form derivative checks in v = (beta, cuts). Class k with + * a = theta - c_{k-1}, b = theta - c_k has log p = log inv_logit(a) + * + log inv_logit(-b) + log(1 - exp(c_{k-1} - c_k)), a sum of univariate + * functions of linear forms in v. Tolerances are relative to the same sums + * over absolute values. + */ +namespace ordered_logistic_test { + +/** One univariate summand f(l . v): value and derivatives 1 to 3. */ +struct term { + Eigen::VectorXd l; + double f[4]; +}; + +template +std::vector terms(const T_y& y, const T_x& x, const Eigen::VectorXd& beta, + const Eigen::VectorXd& cuts) { + using stan::math::inv_logit; + const int P = beta.size(); + const int C = cuts.size(); + const int N = std::max(x.rows(), stan::math::size(y)); + std::vector out; + for (int n = 0; n < N; ++n) { + const int k = stan::math::size(y) == 1 ? stan::get(y, 0) : stan::get(y, n); + Eigen::VectorXd xn = x.row(x.rows() == 1 ? 0 : n).transpose(); + const double theta = xn.dot(beta); + for (int side = 0; side < 2; ++side) { + const int i = side == 0 ? k - 2 : k - 1; // index of the cut + if (i < 0 || i >= C) { + continue; + } + const double t = theta - cuts(i); + const double p = inv_logit(t), q = inv_logit(-t); + term tm{Eigen::VectorXd::Zero(P + C), {0, 0, -p * q, -p * q * (q - p)}}; + tm.l.head(P) = xn; + tm.l(P + i) = -1; + tm.f[0] = side == 0 ? stan::math::log_inv_logit(t) + : stan::math::log_inv_logit(-t); + tm.f[1] = side == 0 ? q : -p; + out.push_back(tm); + } + if (k > 1 && k <= C) { + const double d = cuts(k - 1) - cuts(k - 2); + const double r = 1 / std::expm1(d); + term tm{Eigen::VectorXd::Zero(P + C), + {stan::math::log1m_exp(-d), r, -r * (1 + r), + r * (1 + r) * (1 + 2 * r)}}; + tm.l(P + k - 1) = 1; + tm.l(P + k - 2) = -1; + out.push_back(tm); + } + } + return out; +} + +/** + * Checks value, gradient (var, fvar, fvar, fvar>), + * Hessian (fvar, fvar>, fvar>), third derivatives + * and the derivative of the value itself (fvar, fvar>) of + * `lpmf(y, x, beta, cuts)` in (beta, cuts). + */ +template +void expect_derivatives(const F& lpmf, const T_y& y, const T_x& x, + const Eigen::VectorXd& beta, + const Eigen::VectorXd& cuts) { + using stan::math::fvar; + using stan::math::var; + const int P = beta.size(); + const int M = P + cuts.size(); + Eigen::VectorXd v(M); + v << beta, cuts; + const double tiny = std::numeric_limits::min(); + double lp = 0, lp_tol = 0; + Eigen::VectorXd g = Eigen::VectorXd::Zero(M), g_tol = g; + Eigen::MatrixXd H = Eigen::MatrixXd::Zero(M, M), H_tol = H; + std::vector T(M, H), T_tol(M, H); + for (const auto& tm : terms(y, x, beta, cuts)) { + Eigen::VectorXd a = tm.l.cwiseAbs(); + lp += tm.f[0]; + lp_tol += std::fabs(tm.f[0]); + g += tm.f[1] * tm.l; + g_tol += std::fabs(tm.f[1]) * a; + H += tm.f[2] * tm.l * tm.l.transpose(); + H_tol += std::fabs(tm.f[2]) * a * a.transpose(); + for (int i = 0; i < M; ++i) { + T[i] += tm.f[3] * tm.l(i) * tm.l * tm.l.transpose(); + T_tol[i] += std::fabs(tm.f[3]) * a(i) * a * a.transpose(); + } + } + + auto f = [&](const auto& u) { + return lpmf(y, x, u.head(P).eval(), u.tail(M - P).eval()); + }; + double fx; + Eigen::VectorXd grad, grad_fwd; + Eigen::MatrixXd H_fv, H_ffd, H_ffv; + std::vector grad_H; + stan::math::hessian(f, v, fx, grad, H_fv); + stan::math::hessian(f, v, fx, grad_fwd, H_ffd); + stan::math::grad_hessian(f, v, fx, H_ffv, grad_H); + EXPECT_NEAR(fx, lp, 1e-13 * lp_tol + tiny); + for (int i = 0; i < M; ++i) { + const double gt = 1e-13 * g_tol(i) + tiny; + EXPECT_NEAR(grad(i), g(i), gt) << i; + EXPECT_NEAR(grad_fwd(i), g(i), gt) << i; + for (int j = 0; j < M; ++j) { + const double ht = 1e-13 * H_tol(i, j) + tiny; + EXPECT_NEAR(H_fv(i, j), H(i, j), ht) << i << ", " << j; + EXPECT_NEAR(H_ffd(i, j), H(i, j), ht) << i << ", " << j; + EXPECT_NEAR(H_ffv(i, j), H(i, j), ht) << i << ", " << j; + for (int k = 0; k < M; ++k) { + EXPECT_NEAR(grad_H[i](j, k), T[i](j, k), 1e-13 * T_tol[i](j, k) + tiny) + << i << ", " << j << ", " << k; + } + } + } + + // first order: var, fvar + stan::math::nested_rev_autodiff nested; + Eigen::Matrix v_v(M); + for (int k = 0; k < M; ++k) { + v_v(k) = v(k); + } + f(v_v).grad(); + for (int i = 0; i < M; ++i) { + const double gt = 1e-13 * g_tol(i) + tiny; + EXPECT_NEAR(v_v(i).adj(), g(i), gt) << i; + Eigen::Matrix, Eigen::Dynamic, 1> v_fd(M); + for (int k = 0; k < M; ++k) { + v_fd(k) = fvar(v(k), i == k); + } + EXPECT_NEAR(f(v_fd).d_, g(i), gt) << i; + } + + // derivative of the value itself + nested.set_zero_all_adjoints(); + Eigen::Matrix, Eigen::Dynamic, 1> v_fv(M); + for (int k = 0; k < M; ++k) { + v_fv(k) = fvar(v(k), 0); + } + f(v_fv).val_.grad(); + for (int i = 0; i < M; ++i) { + const double gt = 1e-13 * g_tol(i) + tiny; + EXPECT_NEAR(v_fv(i).val_.adj(), g(i), gt) << i; + Eigen::Matrix>, Eigen::Dynamic, 1> v_ff(M); + for (int k = 0; k < M; ++k) { + v_ff(k) = fvar>(fvar(v(k), i == k), 0); + } + EXPECT_NEAR(f(v_ff).val_.d_, g(i), gt) << i; + } +} + +/** + * Runs `expect_derivatives` at locations equal to a cut (where the third + * derivative is zero, so only orders 1 and 2 test the kink), in the tails, + * between nearly equal cuts, at locations where theta - c_k rounds to + * theta - c_{k-1}, with a single cut, for every class, vector and scalar y, + * and (if `row_x`) a row-vector x. + */ +template +void expect_all(const F& lpmf, bool row_x) { + // dyadic x and beta: x_n beta = 2 x_n0 - x_n1 hits the dyadic cuts exactly + Eigen::VectorXd beta(2); + beta << 2, -1; + Eigen::VectorXd cuts(3); + cuts << -1, 0.5, 2; + std::vector y{1, 2, 2, 3, 3, 4, 1, 4, 2, 3, 1, 4, 2, 3}; + Eigen::VectorXd theta(14); + theta << -1, -1, 0.5, 0.5, 2, 2, -40, 40, -40, 40, 750, -750, 0.25, 1.25; + Eigen::MatrixXd x(14, 2); + x.col(1) << 0.5, -1, 1, 0.25, -0.5, 2, 1.5, -0.75, 0.25, 0.125, -1, -0.25, + 0.75, 0.5; + x.col(0) = (theta + x.col(1)) / 2; + expect_derivatives(lpmf, y, x, beta, cuts); + for (int k = 1; k <= 4; ++k) { + expect_derivatives(lpmf, k, x, beta, cuts); + } + + // cuts about 1e-9 apart; theta at, between and far from them (theta - cut + // rounds at theta = 1000.2) + const double gap = std::ldexp(1.0, -30) + std::ldexp(1.0, -50); + Eigen::VectorXd close(3); + close << -1, 0.5, 0.5 + gap; + std::vector y_close{3, 3, 3, 3, 3, 2, 4}; + Eigen::MatrixXd x_close = Eigen::MatrixXd::Zero(7, 2); + x_close.col(0) << 0.25, 0.25 + gap / 4, 0.25 + gap / 2, 0.75, 500.1, 0.5, + 0.75; + expect_derivatives(lpmf, y_close, x_close, beta, close); + + // theta - c_k rounds to theta - c_{k-1} (theta = 1e17 and 1e14) + Eigen::MatrixXd x_far(1, 2); + x_far << 5e16, 0; + expect_derivatives(lpmf, 2, x_far, beta, Eigen::VectorXd::LinSpaced(2, 0, 1)); + x_far << 5e13, 0; + expect_derivatives(lpmf, 2, x_far, beta, + Eigen::VectorXd::LinSpaced(2, 0, 1e-3)); + + // one cut; theta = 0, 0.5, -40 + Eigen::VectorXd one_cut = Eigen::VectorXd::Zero(1); + Eigen::MatrixXd x_one(3, 2); + x_one << 0, 0, 0.5, 0.5, -20, 0; + expect_derivatives(lpmf, std::vector{1, 2, 1}, x_one, beta, one_cut); + for (int k = 1; k <= 2; ++k) { + expect_derivatives(lpmf, k, x_one, beta, one_cut); + } + + if (row_x) { + for (int n : {0, 2, 6, 7, 10}) { + Eigen::RowVectorXd x_row = x.row(n); + expect_derivatives(lpmf, y, x_row, beta, cuts); + for (int k = 1; k <= 4; ++k) { + expect_derivatives(lpmf, k, x_row, beta, cuts); + } + } + } +} + +} // namespace ordered_logistic_test + +#endif From f66c89a1a2030e56ca9e06cdc39fd61eb273aeb7 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Mon, 12 Oct 2026 11:25:48 +0200 Subject: [PATCH 2/3] ordered_logistic(_glm)_lpmf: one sign partition, cancellation-free partials --- .../kernels/ordered_logistic_glm_lpmf.hpp | 15 +++-- .../opencl/kernels/ordered_logistic_lpmf.hpp | 15 +++-- .../prim/prob/ordered_logistic_glm_lpmf.hpp | 43 +++++++------ stan/math/prim/prob/ordered_logistic_lpmf.hpp | 60 +++++++++---------- 4 files changed, 70 insertions(+), 63 deletions(-) diff --git a/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp b/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp index 751dc9be697..bac3ad89e4c 100644 --- a/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp +++ b/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp @@ -83,19 +83,22 @@ static constexpr const char* ordered_logistic_glm_kernel_code = STRINGIFY( logp -= log1p_exp(-cut2); } if (y != 1 && y != N_classes) { - logp += log1m_exp(cut1 - cut2); + logp += log1m_exp(cut_y2 - cut_y1); } if (need_location_derivative || need_cuts_derivative) { - double exp_cuts_diff = exp(cut_y2 - cut_y1); d1 = inv_logit(-cut2); - d1 -= exp_cuts_diff / (exp_cuts_diff - 1); - d2 = 1 / (1 - exp_cuts_diff); - d2 -= inv_logit(-cut1); - + d2 = inv_logit(cut1); if (need_location_derivative) { location_derivative[gid] = d1 - d2; } + if (need_cuts_derivative) { + // -1 / expm1(c_y - c_{y-1}), zero for the first and last class + const double cuts_diff = cut_y2 - cut_y1; + const double q = exp(cuts_diff) / expm1(cuts_diff); + d1 -= q; + d2 -= q; + } } } } diff --git a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp index 4098872e7e6..fd2e2ee6334 100644 --- a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp +++ b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp @@ -79,19 +79,22 @@ static constexpr const char* ordered_logistic_kernel_code = STRINGIFY( logp -= log1p_exp(-cut2); } if (y != 1 && y != N_classes) { - logp += log1m_exp(cut1 - cut2); + logp += log1m_exp(cut_y2 - cut_y1); } if (need_lambda_derivative || need_cuts_derivative) { - double exp_cuts_diff = exp(cut_y2 - cut_y1); d1 = inv_logit(-cut2); - d1 -= exp_cuts_diff / (exp_cuts_diff - 1); - d2 = 1 / (1 - exp_cuts_diff); - d2 -= inv_logit(-cut1); - + d2 = inv_logit(cut1); if (need_lambda_derivative) { lambda_derivative[gid] = d1 - d2; } + if (need_cuts_derivative) { + // -1 / expm1(c_y - c_{y-1}), zero for the first and last class + const double cuts_diff = cut_y2 - cut_y1; + const double q = exp(cuts_diff) / expm1(cuts_diff); + d1 -= q; + d2 -= q; + } } } } diff --git a/stan/math/prim/prob/ordered_logistic_glm_lpmf.hpp b/stan/math/prim/prob/ordered_logistic_glm_lpmf.hpp index e5855241666..df54557189c 100644 --- a/stan/math/prim/prob/ordered_logistic_glm_lpmf.hpp +++ b/stan/math/prim/prob/ordered_logistic_glm_lpmf.hpp @@ -131,11 +131,18 @@ inline return_type_t ordered_logistic_glm_lpmf( Array cut1 = as_array_or_scalar(location) - cuts_y1; + // exp(-|cut|) on the same sign partition as every select below + Array exp_m_abs_cut1 + = (cut1 > 0.0).select(-cut1, cut1).exp(); + Array exp_m_abs_cut2 + = (cut2 > 0.0).select(-cut2, cut2).exp(); + Array cuts_diff = cuts_y2 - cuts_y1; + // Not immediately evaluating next two expressions benefits performance auto m_log_1p_exp_cut1 - = (cut1 > 0.0).select(-cut1, 0) - (-cut1.abs()).exp().log1p(); + = (cut1 > 0.0).select(-cut1, 0) - exp_m_abs_cut1.log1p(); auto m_log_1p_exp_m_cut2 - = (cut2 <= 0.0).select(cut2, 0) - (-cut2.abs()).exp().log1p(); + = (cut2 <= 0.0).select(cut2, 0) - exp_m_abs_cut2.log1p(); T_partials_return logp(0); if constexpr (is_vector::value) { @@ -145,7 +152,7 @@ inline return_type_t ordered_logistic_glm_lpmf( .select(m_log_1p_exp_cut1, y_vec.cwiseEqual(N_classes).select( m_log_1p_exp_m_cut2, - m_log_1p_exp_m_cut2 + log1m_exp(cut1 - cut2).array() + m_log_1p_exp_m_cut2 + log1m_exp(cuts_diff).array() + m_log_1p_exp_cut1)) .sum(); } else { @@ -154,7 +161,7 @@ inline return_type_t ordered_logistic_glm_lpmf( } else if (y_seq[0] == N_classes) { logp = m_log_1p_exp_m_cut2.sum(); } else { - logp = (m_log_1p_exp_m_cut2 + log1m_exp(cut1 - cut2).array() + logp = (m_log_1p_exp_m_cut2 + log1m_exp(cuts_diff).array() + m_log_1p_exp_cut1) .sum(); } @@ -162,19 +169,15 @@ inline return_type_t ordered_logistic_glm_lpmf( auto ops_partials = make_partials_propagator(x_ref, beta_ref, cuts_ref); if constexpr (is_any_autodiff_v) { - Array exp_m_abs_cut1 = (-cut1.abs()).exp(); - Array exp_m_abs_cut2 = (-cut2.abs()).exp(); - Array exp_cuts_diff = exp(cuts_y2 - cuts_y1); - Array inv_logit_neg_cut2 = (cut2 > 0).select( - exp_m_abs_cut2 / (1 + exp_m_abs_cut2), 1 / (1 + exp_m_abs_cut2)); - Array inv_logit_neg_cut1 = (cut1 > 0).select( - exp_m_abs_cut1 / (1 + exp_m_abs_cut1), 1 / (1 + exp_m_abs_cut1)); - Array d1 - = inv_logit_neg_cut2 - exp_cuts_diff / (exp_cuts_diff - 1); - Array d2 - = 1 / (1 - exp_cuts_diff) - inv_logit_neg_cut1; + Array inv_logit_neg_cut2 + = (cut2 > 0).select(exp_m_abs_cut2, T_partials_return(1)) + / (1 + exp_m_abs_cut2); + Array inv_logit_cut1 + = (cut1 > 0).select(T_partials_return(1), exp_m_abs_cut1) + / (1 + exp_m_abs_cut1); if constexpr (is_any_autodiff_v) { - Matrix location_derivative = d1 - d2; + Matrix location_derivative + = inv_logit_neg_cut2 - inv_logit_cut1; if constexpr (is_autodiff_v) { if constexpr (T_x_rows == 1) { edge<0>(ops_partials).partials_ @@ -187,8 +190,7 @@ inline return_type_t ordered_logistic_glm_lpmf( if constexpr (is_autodiff_v) { if constexpr (T_x_rows == 1) { edge<1>(ops_partials).partials_ - = (location_derivative * x_val.replicate(N_instances, 1)) - .transpose(); + = (location_derivative.sum() * x_val).transpose(); } else { edge<1>(ops_partials).partials_ = (location_derivative * x_val).transpose(); @@ -196,6 +198,11 @@ inline return_type_t ordered_logistic_glm_lpmf( } } if constexpr (is_autodiff_v) { + // -1 / expm1(c_y - c_{y-1}), zero for the first and last class + Array q + = cuts_diff.exp() / cuts_diff.expm1(); + Array d1 = inv_logit_neg_cut2 - q; + Array d2 = inv_logit_cut1 - q; for (int i = 0; i < N_instances; i++) { int c = y_seq[i]; if (c != N_classes) { diff --git a/stan/math/prim/prob/ordered_logistic_lpmf.hpp b/stan/math/prim/prob/ordered_logistic_lpmf.hpp index db2b7d86d35..f077fa38d58 100644 --- a/stan/math/prim/prob/ordered_logistic_lpmf.hpp +++ b/stan/math/prim/prob/ordered_logistic_lpmf.hpp @@ -6,13 +6,12 @@ #include #include #include -#include #include -#include -#include +#include #include #include #include +#include #include #include #include @@ -30,25 +29,13 @@ namespace math { * will be the dot product of a vector of regression coefficients * and a vector of predictors for the outcome * - \f[ - \frac{\partial }{\partial \lambda} = - \begin{cases}\\ - -\mathrm{logit}^{-1}(\lambda - c_1) & \mbox{if } k = 1,\\ - -(((1-e^{c_{k-1}-c_{k-2}})^{-1} - \mathrm{logit}^{-1}(c_{k-2}-\lambda)) + - ((1-e^{c_{k-2}-c_{k-1}})^{-1} - \mathrm{logit}^{-1}(c_{k-1}-\lambda))) - & \mathrm{if } 1 < k < K, \mathrm{and}\\ - \mathrm{logit}^{-1}(c_{K-2}-\lambda) & \mathrm{if } k = K. - \end{cases} - \f] - \f[ \frac{\partial }{\partial \lambda} = \begin{cases} -\mathrm{logit}^{-1}(\lambda - c_1) & \text{if } k = 1,\\ - -(((1-e^{c_{k-1}-c_{k-2}})^{-1} - \mathrm{logit}^{-1}(c_{k-2}-\lambda)) + - ((1-e^{c_{k-2}-c_{k-1}})^{-1} - \mathrm{logit}^{-1}(c_{k-1}-\lambda))) + \mathrm{logit}^{-1}(c_{k-1} - \lambda) - \mathrm{logit}^{-1}(\lambda - c_k) & \text{if } 1 < k < K, \text{ and}\\ - \mathrm{logit}^{-1}(c_{K-2}-\lambda) & \text{if } k = K. + \mathrm{logit}^{-1}(c_{K-1} - \lambda) & \text{if } k = K. \end{cases} \f] * @@ -144,16 +131,23 @@ inline return_type_t ordered_logistic_lpmf(const T_y& y, Array cut2 = lambda_val - cuts_y2; Array cut1 = lambda_val - cuts_y1; + // exp(-|cut|) on the same sign partition as every select below + Array exp_m_abs_cut1 + = (cut1 > 0.0).select(-cut1, cut1).exp(); + Array exp_m_abs_cut2 + = (cut2 > 0.0).select(-cut2, cut2).exp(); + Array cuts_diff = cuts_y2 - cuts_y1; + // Not immediately evaluating next two expressions benefits performance auto m_log_1p_exp_cut1 - = (cut1 > 0.0).select(-cut1, 0) - (-cut1.abs()).exp().log1p(); + = (cut1 > 0.0).select(-cut1, 0) - exp_m_abs_cut1.log1p(); auto m_log_1p_exp_m_cut2 - = (cut2 <= 0.0).select(cut2, 0) - (-cut2.abs()).exp().log1p(); + = (cut2 <= 0.0).select(cut2, 0) - exp_m_abs_cut2.log1p(); if constexpr (is_vector::value) { Eigen::Map, Eigen::Dynamic, 1>> y_vec( y_seq.data(), y_seq.size()); - auto log1m_exp_cuts_diff = log1m_exp(cut1 - cut2); + auto log1m_exp_cuts_diff = log1m_exp(cuts_diff); logp = y_vec.cwiseEqual(1) .select(m_log_1p_exp_cut1, y_vec.cwiseEqual(K).select(m_log_1p_exp_m_cut2, @@ -167,7 +161,7 @@ inline return_type_t ordered_logistic_lpmf(const T_y& y, } else if (y_seq[0] == K) { logp = m_log_1p_exp_m_cut2.sum(); } else { - logp = (m_log_1p_exp_m_cut2 + log1m_exp(cut1 - cut2).array() + logp = (m_log_1p_exp_m_cut2 + log1m_exp(cuts_diff).array() + m_log_1p_exp_cut1) .sum(); } @@ -175,21 +169,21 @@ inline return_type_t ordered_logistic_lpmf(const T_y& y, auto ops_partials = make_partials_propagator(lambda_ref, c_ref); if constexpr (is_any_autodiff_v) { - Array exp_m_abs_cut1 = (-cut1.abs()).exp(); - Array exp_m_abs_cut2 = (-cut2.abs()).exp(); - Array exp_cuts_diff = exp(cuts_y2 - cuts_y1); - Array inv_logit_neg_cut2 = (cut2 > 0).select( - exp_m_abs_cut2 / (1 + exp_m_abs_cut2), 1 / (1 + exp_m_abs_cut2)); - Array inv_logit_neg_cut1 = (cut1 > 0).select( - exp_m_abs_cut1 / (1 + exp_m_abs_cut1), 1 / (1 + exp_m_abs_cut1)); - Array d1 - = inv_logit_neg_cut2 - exp_cuts_diff / (exp_cuts_diff - 1); - Array d2 - = 1 / (1 - exp_cuts_diff) - inv_logit_neg_cut1; + Array inv_logit_neg_cut2 + = (cut2 > 0).select(exp_m_abs_cut2, T_partials_return(1)) + / (1 + exp_m_abs_cut2); + Array inv_logit_cut1 + = (cut1 > 0).select(T_partials_return(1), exp_m_abs_cut1) + / (1 + exp_m_abs_cut1); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = d1 - d2; + partials<0>(ops_partials) = inv_logit_neg_cut2 - inv_logit_cut1; } if constexpr (is_autodiff_v) { + // -1 / expm1(c_y - c_{y-1}), zero for the first and last class + Array q + = cuts_diff.exp() / cuts_diff.expm1(); + Array d1 = inv_logit_neg_cut2 - q; + Array d2 = inv_logit_cut1 - q; for (int i = 0; i < N; i++) { int c = y_seq[i]; if (c != K) { From 73bcf000cd63fb2dd2f24c06be4b324c4f140591 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 19 Oct 2026 09:11:29 +0300 Subject: [PATCH 3/3] Test the OpenCL ordered_logistic kernels against references in the tails --- .../rev/ordered_logistic_glm_lpmf_test.cpp | 30 ++++++++++ .../opencl/rev/ordered_logistic_lpmf_test.cpp | 26 +++++++++ .../opencl/rev/ordered_logistic_tail_refs.hpp | 56 +++++++++++++++++++ 3 files changed, 112 insertions(+) create mode 100644 test/unit/math/opencl/rev/ordered_logistic_tail_refs.hpp diff --git a/test/unit/math/opencl/rev/ordered_logistic_glm_lpmf_test.cpp b/test/unit/math/opencl/rev/ordered_logistic_glm_lpmf_test.cpp index 9eb063f8a89..6b9a3fc30e0 100644 --- a/test/unit/math/opencl/rev/ordered_logistic_glm_lpmf_test.cpp +++ b/test/unit/math/opencl/rev/ordered_logistic_glm_lpmf_test.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include using Eigen::Array; @@ -222,4 +223,33 @@ TEST(ProbDistributionsOrderedLogisitcGLM, opencl_matches_cpu_big) { ordered_logistic_glm_lpmf_functor_propto, y, x, beta, cuts); } +TEST(ProbDistributionsOrderedLogisitcGLM, opencl_rev_tails_against_references) { + using ordered_logistic_tail_refs::expect_rel; + using stan::math::var_value; + // location = x * beta with x = 1, so d/dbeta = d/dlambda + for (const auto& r : ordered_logistic_tail_refs::refs()) { + vector y{r.y}; + Matrix x(1, 1); + x << 1.0; + Matrix beta(1); + beta << r.lambda; + Matrix cuts(2); + cuts << r.c1, r.c2; + matrix_cl y_cl(y); + matrix_cl x_cl(x); + var_value> beta_cl = stan::math::to_matrix_cl(beta); + var_value> cuts_cl = stan::math::to_matrix_cl(cuts); + var lp + = stan::math::ordered_logistic_glm_lpmf(y_cl, x_cl, beta_cl, cuts_cl); + lp.grad(); + Eigen::VectorXd d_beta = stan::math::from_matrix_cl(beta_cl.adj()); + Eigen::VectorXd d_cuts = stan::math::from_matrix_cl(cuts_cl.adj()); + expect_rel(r.val, lp.val(), "value", r.lambda); + expect_rel(r.d_lambda, d_beta(0), "d/dbeta", r.lambda, r.d_lambda_abs_tol); + expect_rel(r.d_c1, d_cuts(0), "d/dc1", r.lambda); + expect_rel(r.d_c2, d_cuts(1), "d/dc2", r.lambda); + stan::math::recover_memory(); + } +} + #endif diff --git a/test/unit/math/opencl/rev/ordered_logistic_lpmf_test.cpp b/test/unit/math/opencl/rev/ordered_logistic_lpmf_test.cpp index bf8fefd862d..a909dc8b835 100644 --- a/test/unit/math/opencl/rev/ordered_logistic_lpmf_test.cpp +++ b/test/unit/math/opencl/rev/ordered_logistic_lpmf_test.cpp @@ -1,5 +1,6 @@ #ifdef STAN_OPENCL #include +#include #include #include #include @@ -175,4 +176,29 @@ TEST(ProbDistributionsOrderedLogisitc, opencl_matches_cpu_big) { ordered_logistic_lpmf_functor_propto, y, lambda, cuts); } +TEST(ProbDistributionsOrderedLogisitc, opencl_rev_tails_against_references) { + using ordered_logistic_tail_refs::expect_rel; + using stan::math::var_value; + for (const auto& r : ordered_logistic_tail_refs::refs()) { + vector y{r.y}; + Eigen::VectorXd lambda(1); + lambda << r.lambda; + Eigen::VectorXd cuts(2); + cuts << r.c1, r.c2; + matrix_cl y_cl(y); + var_value> lambda_cl = stan::math::to_matrix_cl(lambda); + var_value> cuts_cl = stan::math::to_matrix_cl(cuts); + var lp = stan::math::ordered_logistic_lpmf(y_cl, lambda_cl, cuts_cl); + lp.grad(); + Eigen::VectorXd d_lambda = stan::math::from_matrix_cl(lambda_cl.adj()); + Eigen::VectorXd d_cuts = stan::math::from_matrix_cl(cuts_cl.adj()); + expect_rel(r.val, lp.val(), "value", r.lambda); + expect_rel(r.d_lambda, d_lambda(0), "d/dlambda", r.lambda, + r.d_lambda_abs_tol); + expect_rel(r.d_c1, d_cuts(0), "d/dc1", r.lambda); + expect_rel(r.d_c2, d_cuts(1), "d/dc2", r.lambda); + stan::math::recover_memory(); + } +} + #endif diff --git a/test/unit/math/opencl/rev/ordered_logistic_tail_refs.hpp b/test/unit/math/opencl/rev/ordered_logistic_tail_refs.hpp new file mode 100644 index 00000000000..c41814ab029 --- /dev/null +++ b/test/unit/math/opencl/rev/ordered_logistic_tail_refs.hpp @@ -0,0 +1,56 @@ +#ifndef TEST_UNIT_MATH_OPENCL_REV_ORDERED_LOGISTIC_TAIL_REFS_HPP +#define TEST_UNIT_MATH_OPENCL_REV_ORDERED_LOGISTIC_TAIL_REFS_HPP + +#include +#include +#include +#include + +namespace ordered_logistic_tail_refs { + +// Category y, location lambda and cut points (c1, c2); the log probability +// and its derivatives with respect to lambda, c1 and c2. mpmath, 60 digits, +// rounded to double. For close cut points the lambda derivative is a +// difference of two values near 1/2; d_lambda_abs_tol > 0 gives an absolute +// tolerance for it. +struct ref { + int y; + double lambda; + double c1; + double c2; + double val; + double d_lambda; + double d_c1; + double d_c2; + double d_lambda_abs_tol; +}; + +// clang-format off +// NOLINTBEGIN(whitespace/line_length) +inline std::vector refs() { + return { + {2, 1e17, 0, 1, -1e17, -1.0, -0.5819767068693265, 1.5819767068693265, 0}, + {2, -1e17, 0, 1, -1e17, 1.0, -1.5819767068693265, 0.5819767068693265, 0}, + {2, 0.3, 0, 1, -1.416216438741067, 0.09374525535650713, -1.0075341900576675, 0.9137889347011603, 0}, + {1, -40, 0, 1, -4.248354255291589e-18, -4.248354255291589e-18, 4.248354255291589e-18, 0, 0}, + {3, 41, 0, 1, -4.248354255291589e-18, 4.248354255291589e-18, 0, -4.248354255291589e-18, 0}, + {2, 0, 0, 9.313225746154785e-10, -22.18070977791825, 2.3283064365386963e-10, -1073741824.0, 1073741824.0, 1e-15}}; +} +// NOLINTEND +// clang-format on + +// Relative check, absolute below the smallest normal double, or absolute +// with tolerance abs_tol > 0. +inline void expect_rel(double expected, double actual, const char* what, + double lambda, double abs_tol = 0) { + const double tol = abs_tol > 0 + ? abs_tol + : 1e-13 + * std::fmax(std::fabs(expected), + std::numeric_limits::min()); + EXPECT_NEAR(expected, actual, tol) << what << " at lambda = " << lambda; +} + +} // namespace ordered_logistic_tail_refs + +#endif