From 7bf50285979375a1733818aa83ffbc3bee52c4db Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 13 Sep 2026 14:53:55 +0300 Subject: [PATCH 01/10] Stabilize logistic location gradients --- stan/math/opencl/prim/logistic_lpdf.hpp | 8 +------- stan/math/prim/prob/logistic_lpdf.hpp | 6 ++---- test/unit/math/opencl/rev/logistic_lpdf_test.cpp | 14 ++++++++++++++ test/unit/math/rev/prob/logistic_lpdf_test.cpp | 14 ++++++++++++++ 4 files changed, 31 insertions(+), 11 deletions(-) create mode 100644 test/unit/math/rev/prob/logistic_lpdf_test.cpp diff --git a/stan/math/opencl/prim/logistic_lpdf.hpp b/stan/math/opencl/prim/logistic_lpdf.hpp index cfcfb2b5ba7..f6a9b1dd23d 100644 --- a/stan/math/opencl/prim/logistic_lpdf.hpp +++ b/stan/math/opencl/prim/logistic_lpdf.hpp @@ -83,14 +83,8 @@ inline return_type_t logistic_lpdf( auto y_deriv = elt_multiply( elt_divide(2.0, 1.0 + exp(y_minus_mu_div_sigma)) - 1.0, inv_sigma); - auto exp_mu_div_sigma = exp(elt_multiply(mu_val, inv_sigma)); auto mu_deriv = elt_multiply( - 1.0 - - 2.0 - * elt_divide( - exp_mu_div_sigma, - exp_mu_div_sigma + exp(elt_multiply(y_val, inv_sigma))), - inv_sigma); + tanh(0.5 * y_minus_mu_div_sigma), inv_sigma); auto sigma_deriv = elt_multiply(-elt_multiply(y_deriv, y_minus_mu) - 1.0, inv_sigma); diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp index 181ce1c5650..ce4bb340016 100644 --- a/stan/math/prim/prob/logistic_lpdf.hpp +++ b/stan/math/prim/prob/logistic_lpdf.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -82,11 +83,8 @@ inline return_type_t logistic_lpdf(const T_y& y, } } if constexpr (is_autodiff_v) { - const auto& exp_mu_div_sigma = to_ref(exp(mu_val * inv_sigma)); edge<1>(ops_partials).partials_ - = (1 - - 2 * exp_mu_div_sigma / (exp_mu_div_sigma + exp(y_val * inv_sigma))) - * inv_sigma; + = tanh(0.5 * y_minus_mu_div_sigma) * inv_sigma; } return ops_partials.build(logp); } diff --git a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp index 44cc923c6ea..0c2258f6373 100644 --- a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp @@ -100,6 +100,20 @@ TEST(ProbDistributionsLogistic, opencl_matches_cpu_small) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogistic, opencl_matches_cpu_large_location) { + Eigen::VectorXd y(1); + y << 711; + Eigen::VectorXd mu(1); + mu << 710; + Eigen::VectorXd sigma(1); + sigma << 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lpdf_functor, y, mu, + sigma); + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lpdf_functor_propto, y, + mu, sigma); +} + TEST(ProbDistributionsLogistic, opencl_broadcast_y) { int N = 3; diff --git a/test/unit/math/rev/prob/logistic_lpdf_test.cpp b/test/unit/math/rev/prob/logistic_lpdf_test.cpp new file mode 100644 index 00000000000..67ad7ad8a1e --- /dev/null +++ b/test/unit/math/rev/prob/logistic_lpdf_test.cpp @@ -0,0 +1,14 @@ +#include +#include +#include +#include + +TEST_F(AgradRev, logistic_lpdf_location_gradient_large_values) { + stan::math::var mu = 710.0; + + stan::math::var logp = stan::math::logistic_lpdf(711.0, mu, 1.0); + logp.grad(); + + EXPECT_TRUE(std::isfinite(logp.val())); + EXPECT_DOUBLE_EQ(std::tanh(0.5), mu.adj()); +} From 777a274a6f491f1990df7ef989b95dcdf0bc02ce Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Sep 2026 08:15:13 -0400 Subject: [PATCH 02/10] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/opencl/prim/logistic_lpdf.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/opencl/prim/logistic_lpdf.hpp b/stan/math/opencl/prim/logistic_lpdf.hpp index f6a9b1dd23d..d976363bc69 100644 --- a/stan/math/opencl/prim/logistic_lpdf.hpp +++ b/stan/math/opencl/prim/logistic_lpdf.hpp @@ -83,8 +83,7 @@ inline return_type_t logistic_lpdf( auto y_deriv = elt_multiply( elt_divide(2.0, 1.0 + exp(y_minus_mu_div_sigma)) - 1.0, inv_sigma); - auto mu_deriv = elt_multiply( - tanh(0.5 * y_minus_mu_div_sigma), inv_sigma); + auto mu_deriv = elt_multiply(tanh(0.5 * y_minus_mu_div_sigma), inv_sigma); auto sigma_deriv = elt_multiply(-elt_multiply(y_deriv, y_minus_mu) - 1.0, inv_sigma); From 6c71ece7594fc6ded7025b80cd4192ff86491530 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 13 Sep 2026 17:25:00 +0300 Subject: [PATCH 03/10] Stabilize logistic log-CDF tails --- stan/math/opencl/prim/logistic_lccdf.hpp | 14 +++---- stan/math/opencl/prim/logistic_lcdf.hpp | 14 +++---- stan/math/prim/prob/logistic_lccdf.hpp | 39 +++++++------------ stan/math/prim/prob/logistic_lcdf.hpp | 37 +++++++----------- .../math/opencl/rev/logistic_lccdf_test.cpp | 10 +++++ .../math/opencl/rev/logistic_lcdf_test.cpp | 10 +++++ .../math/rev/prob/logistic_lccdf_test.cpp | 21 ++++++++++ .../unit/math/rev/prob/logistic_lcdf_test.cpp | 21 ++++++++++ 8 files changed, 102 insertions(+), 64 deletions(-) create mode 100644 test/unit/math/rev/prob/logistic_lccdf_test.cpp create mode 100644 test/unit/math/rev/prob/logistic_lcdf_test.cpp diff --git a/stan/math/opencl/prim/logistic_lccdf.hpp b/stan/math/opencl/prim/logistic_lccdf.hpp index 2c9155969cf..4552d58bbb7 100644 --- a/stan/math/opencl/prim/logistic_lccdf.hpp +++ b/stan/math/opencl/prim/logistic_lccdf.hpp @@ -66,16 +66,12 @@ inline return_type_t logistic_lccdf( auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto any_y_pos_inf = colwise_max(cast(y_val == INFTY)); auto inv_sigma = elt_divide(1.0, sigma_val); - auto mu_minus_y_div_sigma = elt_multiply(mu_val - y_val, inv_sigma); - auto exp_scaled_diff = exp(mu_minus_y_div_sigma); - auto Pn = 1.0 - elt_divide(1.0, 1.0 + exp_scaled_diff); - auto P_expr = colwise_sum(log(Pn)); - - auto mu_deriv = elt_divide( - exp(mu_minus_y_div_sigma - log(sigma_val) - 2.0 * log1p(exp_scaled_diff)), - Pn); + auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); + auto P_expr = colwise_sum(log1m_inv_logit(scaled_diff)); + + auto mu_deriv = elt_multiply(inv_logit(scaled_diff), inv_sigma); auto y_deriv = -mu_deriv; - auto sigma_deriv = elt_multiply(-mu_deriv, mu_minus_y_div_sigma); + auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff); matrix_cl any_y_neg_inf_cl; matrix_cl any_y_pos_inf_cl; diff --git a/stan/math/opencl/prim/logistic_lcdf.hpp b/stan/math/opencl/prim/logistic_lcdf.hpp index f16c9518bd1..06ae8dbec99 100644 --- a/stan/math/opencl/prim/logistic_lcdf.hpp +++ b/stan/math/opencl/prim/logistic_lcdf.hpp @@ -66,16 +66,12 @@ inline return_type_t logistic_lcdf( auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto cond = y_val == INFTY; auto inv_sigma = elt_divide(1.0, sigma_val); - auto mu_minus_y_div_sigma = elt_multiply(mu_val - y_val, inv_sigma); - auto exp_scaled_diff = exp(mu_minus_y_div_sigma); - auto Pn = elt_divide(1.0, 1.0 + exp_scaled_diff); - auto P_expr = colwise_sum(log(Pn)); - - auto y_deriv = elt_divide( - exp(mu_minus_y_div_sigma - log(sigma_val) - 2.0 * log1p(exp_scaled_diff)), - Pn); + auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); + auto P_expr = colwise_sum(log_inv_logit(scaled_diff)); + + auto y_deriv = elt_multiply(inv_logit(-scaled_diff), inv_sigma); auto mu_deriv = -y_deriv; - auto sigma_deriv = elt_multiply(y_deriv, mu_minus_y_div_sigma); + auto sigma_deriv = elt_multiply(-y_deriv, scaled_diff); matrix_cl any_y_neg_inf_cl; matrix_cl P_cl; diff --git a/stan/math/prim/prob/logistic_lccdf.hpp b/stan/math/prim/prob/logistic_lccdf.hpp index 016da2da39d..95e3bea4341 100644 --- a/stan/math/prim/prob/logistic_lccdf.hpp +++ b/stan/math/prim/prob/logistic_lccdf.hpp @@ -4,17 +4,14 @@ #include #include #include -#include -#include #include +#include #include #include #include #include #include -#include #include -#include namespace stan { namespace math { @@ -26,8 +23,6 @@ inline return_type_t logistic_lccdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; using T_sigma_ref = ref_type_t; @@ -70,26 +65,22 @@ inline return_type_t logistic_lccdf(const T_y& y, const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); - const T_partials_return sigma_dbl = sigma_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); + const T_partials_return scaled_diff + = (y_dbl - mu_dbl) * sigma_inv_vec; + P += log1m_inv_logit(scaled_diff); - // TODO(Andrew) Further simplify derivatives and log-scale below - const T_partials_return Pn - = 1.0 - inv_logit((y_dbl - mu_dbl) * sigma_inv_vec); - P += log(Pn); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - -= exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - -= -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - -= -(y_dbl - mu_dbl) * sigma_inv_vec - * exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; + if constexpr (is_any_autodiff_v) { + const T_partials_return deriv = inv_logit(scaled_diff) * sigma_inv_vec; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] -= deriv; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] += deriv; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] += scaled_diff * deriv; + } } } return ops_partials.build(P); diff --git a/stan/math/prim/prob/logistic_lcdf.hpp b/stan/math/prim/prob/logistic_lcdf.hpp index 10b9be9e78c..b0faa88fdc4 100644 --- a/stan/math/prim/prob/logistic_lcdf.hpp +++ b/stan/math/prim/prob/logistic_lcdf.hpp @@ -4,17 +4,14 @@ #include #include #include -#include -#include #include +#include #include #include #include #include #include -#include #include -#include namespace stan { namespace math { @@ -26,8 +23,6 @@ inline return_type_t logistic_lcdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; using T_sigma_ref = ref_type_t; @@ -70,25 +65,23 @@ inline return_type_t logistic_lcdf(const T_y& y, const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); - const T_partials_return sigma_dbl = sigma_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - // TODO(Andrew) Further simplify derivatives and log-scale below - const T_partials_return Pn = inv_logit((y_dbl - mu_dbl) * sigma_inv_vec); - P += log(Pn); + const T_partials_return scaled_diff + = (y_dbl - mu_dbl) * sigma_inv_vec; + P += log_inv_logit(scaled_diff); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - += exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += -(y_dbl - mu_dbl) * sigma_inv_vec - * exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; + if constexpr (is_any_autodiff_v) { + const T_partials_return deriv = inv_logit(-scaled_diff) * sigma_inv_vec; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] += deriv; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] -= deriv; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] -= scaled_diff * deriv; + } } } return ops_partials.build(P); diff --git a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp index 50796ef91b6..e4642438a83 100644 --- a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp @@ -79,6 +79,16 @@ TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_small) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_upper_tail) { + Eigen::VectorXd y(2); + y << 36.5, 40; + Eigen::VectorXd mu = Eigen::VectorXd::Zero(2); + Eigen::VectorXd sigma = Eigen::VectorXd::Ones(2); + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_small_y_neg_inf) { int N = 3; int M = 2; diff --git a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp index c3cb1cc9d4e..ed812f70d39 100644 --- a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp @@ -79,6 +79,16 @@ TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_small) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_lower_tail) { + Eigen::VectorXd y(2); + y << -745, -746; + Eigen::VectorXd mu = Eigen::VectorXd::Zero(2); + Eigen::VectorXd sigma = Eigen::VectorXd::Ones(2); + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_small_y_neg_inf) { int N = 3; int M = 2; diff --git a/test/unit/math/rev/prob/logistic_lccdf_test.cpp b/test/unit/math/rev/prob/logistic_lccdf_test.cpp new file mode 100644 index 00000000000..e08c2afff17 --- /dev/null +++ b/test/unit/math/rev/prob/logistic_lccdf_test.cpp @@ -0,0 +1,21 @@ +#include +#include +#include + +TEST_F(AgradRev, logistic_lccdf_upper_tail) { + for (double y_value : {36.5, 40.0}) { + stan::math::var y = y_value; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + const double deriv = stan::math::inv_logit(y_value); + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(y_value), log_ccdf.val()); + EXPECT_DOUBLE_EQ(-deriv, y.adj()); + EXPECT_DOUBLE_EQ(deriv, mu.adj()); + EXPECT_DOUBLE_EQ(y_value * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} diff --git a/test/unit/math/rev/prob/logistic_lcdf_test.cpp b/test/unit/math/rev/prob/logistic_lcdf_test.cpp new file mode 100644 index 00000000000..40146563763 --- /dev/null +++ b/test/unit/math/rev/prob/logistic_lcdf_test.cpp @@ -0,0 +1,21 @@ +#include +#include +#include + +TEST_F(AgradRev, logistic_lcdf_lower_tail) { + for (double y_value : {-745.0, -746.0}) { + stan::math::var y = y_value; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = stan::math::inv_logit(-y_value); + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(-y_value), log_cdf.val()); + EXPECT_DOUBLE_EQ(deriv, y.adj()); + EXPECT_DOUBLE_EQ(-deriv, mu.adj()); + EXPECT_DOUBLE_EQ(-y_value * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} From d8f37961ab37ad3acab809fecce05749e8482f8f Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Sep 2026 14:43:08 -0400 Subject: [PATCH 04/10] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/logistic_lccdf.hpp | 3 +-- stan/math/prim/prob/logistic_lcdf.hpp | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/prob/logistic_lccdf.hpp b/stan/math/prim/prob/logistic_lccdf.hpp index 95e3bea4341..13b4e05fa6c 100644 --- a/stan/math/prim/prob/logistic_lccdf.hpp +++ b/stan/math/prim/prob/logistic_lccdf.hpp @@ -66,8 +66,7 @@ inline return_type_t logistic_lccdf(const T_y& y, const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - const T_partials_return scaled_diff - = (y_dbl - mu_dbl) * sigma_inv_vec; + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; P += log1m_inv_logit(scaled_diff); if constexpr (is_any_autodiff_v) { diff --git a/stan/math/prim/prob/logistic_lcdf.hpp b/stan/math/prim/prob/logistic_lcdf.hpp index b0faa88fdc4..e0ebfd0b8e6 100644 --- a/stan/math/prim/prob/logistic_lcdf.hpp +++ b/stan/math/prim/prob/logistic_lcdf.hpp @@ -67,8 +67,7 @@ inline return_type_t logistic_lcdf(const T_y& y, const T_partials_return mu_dbl = mu_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - const T_partials_return scaled_diff - = (y_dbl - mu_dbl) * sigma_inv_vec; + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; P += log_inv_logit(scaled_diff); if constexpr (is_any_autodiff_v) { From fac79ffd35d50d28c4336747d23ee74cf24071f6 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 23:41:03 +0300 Subject: [PATCH 05/10] Stabilize logistic_cdf tails and derive all logistic_lpdf gradients from tanh --- stan/math/opencl/prim/logistic_cdf.hpp | 16 ++--- stan/math/opencl/prim/logistic_lcdf.hpp | 8 ++- stan/math/opencl/prim/logistic_lpdf.hpp | 8 ++- stan/math/prim/prob/logistic_cdf.hpp | 36 +++++------ stan/math/prim/prob/logistic_lccdf.hpp | 14 ++--- stan/math/prim/prob/logistic_lcdf.hpp | 3 +- stan/math/prim/prob/logistic_lpdf.hpp | 25 ++++---- .../math/opencl/rev/logistic_cdf_test.cpp | 26 ++++++++ .../math/opencl/rev/logistic_lccdf_test.cpp | 24 ++++++-- .../math/opencl/rev/logistic_lcdf_test.cpp | 20 ++++++- .../math/opencl/rev/logistic_lpdf_test.cpp | 12 ++-- test/unit/math/rev/prob/logistic_cdf_test.cpp | 58 ++++++++++++++++++ .../math/rev/prob/logistic_lccdf_test.cpp | 60 +++++++++++++++++++ .../unit/math/rev/prob/logistic_lcdf_test.cpp | 43 +++++++++++++ .../unit/math/rev/prob/logistic_lpdf_test.cpp | 51 +++++++++++++++- 15 files changed, 335 insertions(+), 69 deletions(-) create mode 100644 test/unit/math/rev/prob/logistic_cdf_test.cpp diff --git a/stan/math/opencl/prim/logistic_cdf.hpp b/stan/math/opencl/prim/logistic_cdf.hpp index 7d02536baa3..307f17460fd 100644 --- a/stan/math/opencl/prim/logistic_cdf.hpp +++ b/stan/math/opencl/prim/logistic_cdf.hpp @@ -66,16 +66,16 @@ inline return_type_t logistic_cdf( auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto cond = y_val == INFTY; auto inv_sigma = elt_divide(1.0, sigma_val); - auto mu_minus_y_div_sigma = elt_multiply(mu_val - y_val, inv_sigma); - auto exp_scaled_diff = exp(mu_minus_y_div_sigma); - auto Pn = elt_divide(1.0, 1.0 + exp_scaled_diff); + auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); + auto Pn = inv_logit(scaled_diff); auto P_expr = colwise_prod(select(cond, 1.0, Pn)); - auto y_deriv_tmp = select(cond, 0.0, - elt_divide(exp(mu_minus_y_div_sigma - log(sigma_val) - - 2.0 * log1p(exp_scaled_diff)), - Pn)); - auto sigma_deriv_tmp = elt_multiply(y_deriv_tmp, mu_minus_y_div_sigma); + // These are the log-scale derivatives; they are rescaled by the product P + // below. inv_logit(-scaled_diff) avoids the pdf / Pn quotient, which is + // 0 / 0 once Pn underflows. + auto deriv = elt_multiply(inv_logit(-scaled_diff), inv_sigma); + auto y_deriv_tmp = select(cond, 0.0, deriv); + auto sigma_deriv_tmp = select(cond, 0.0, elt_multiply(-deriv, scaled_diff)); matrix_cl any_y_neg_inf_cl; matrix_cl P_cl; diff --git a/stan/math/opencl/prim/logistic_lcdf.hpp b/stan/math/opencl/prim/logistic_lcdf.hpp index 06ae8dbec99..b9706fe715b 100644 --- a/stan/math/opencl/prim/logistic_lcdf.hpp +++ b/stan/math/opencl/prim/logistic_lcdf.hpp @@ -69,9 +69,13 @@ inline return_type_t logistic_lcdf( auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); auto P_expr = colwise_sum(log_inv_logit(scaled_diff)); - auto y_deriv = elt_multiply(inv_logit(-scaled_diff), inv_sigma); + // y == INFTY contributes log(1) = 0 to P and zero to every partial; without + // the select the scale partial would be 0 * INFTY = NaN, which prim (where + // the element is skipped outright) never produces. + auto deriv = elt_multiply(inv_logit(-scaled_diff), inv_sigma); + auto y_deriv = select(cond, 0.0, deriv); auto mu_deriv = -y_deriv; - auto sigma_deriv = elt_multiply(-y_deriv, scaled_diff); + auto sigma_deriv = select(cond, 0.0, elt_multiply(-deriv, scaled_diff)); matrix_cl any_y_neg_inf_cl; matrix_cl P_cl; diff --git a/stan/math/opencl/prim/logistic_lpdf.hpp b/stan/math/opencl/prim/logistic_lpdf.hpp index d976363bc69..3ef004d2df9 100644 --- a/stan/math/opencl/prim/logistic_lpdf.hpp +++ b/stan/math/opencl/prim/logistic_lpdf.hpp @@ -81,11 +81,13 @@ inline return_type_t logistic_lpdf( = colwise_sum(static_select::value>( logp1 - log(sigma_val), logp1)); - auto y_deriv = elt_multiply( - elt_divide(2.0, 1.0 + exp(y_minus_mu_div_sigma)) - 1.0, inv_sigma); + // d/dmu = tanh(z / 2) / sigma with z = (y - mu) / sigma. The y and sigma + // partials reuse it so that d/dy == -d/dmu exactly; the equivalent + // 2 / (1 + exp(z)) - 1 form loses all relative precision as z -> 0. auto mu_deriv = elt_multiply(tanh(0.5 * y_minus_mu_div_sigma), inv_sigma); + auto y_deriv = -mu_deriv; auto sigma_deriv - = elt_multiply(-elt_multiply(y_deriv, y_minus_mu) - 1.0, inv_sigma); + = elt_multiply(elt_multiply(mu_deriv, y_minus_mu) - 1.0, inv_sigma); matrix_cl logp_cl; matrix_cl y_deriv_cl; diff --git a/stan/math/prim/prob/logistic_cdf.hpp b/stan/math/prim/prob/logistic_cdf.hpp index 3465e44eae6..fe659f91333 100644 --- a/stan/math/prim/prob/logistic_cdf.hpp +++ b/stan/math/prim/prob/logistic_cdf.hpp @@ -4,16 +4,13 @@ #include #include #include -#include #include #include #include #include #include #include -#include #include -#include namespace stan { namespace math { @@ -26,7 +23,6 @@ inline return_type_t logistic_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; - using std::exp; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; using T_sigma_ref = ref_type_t; @@ -69,26 +65,26 @@ inline return_type_t logistic_cdf(const T_y& y, const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); - const T_partials_return sigma_dbl = sigma_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - - // TODO(Andrew) Further simplify derivatives and log scale below - const T_partials_return Pn = inv_logit((y_dbl - mu_dbl) * sigma_inv_vec); + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; + const T_partials_return Pn = inv_logit(scaled_diff); P *= Pn; - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - += exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += -(y_dbl - mu_dbl) * sigma_inv_vec - * exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; + // The partials accumulate d/d. log(Pn); they are rescaled by the product + // P below. Writing the log-scale derivative as inv_logit(-scaled_diff) + // avoids the pdf / Pn quotient, which is 0 / 0 once Pn underflows. + if constexpr (is_any_autodiff_v) { + const T_partials_return deriv = inv_logit(-scaled_diff) * sigma_inv_vec; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] += deriv; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] -= deriv; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] -= scaled_diff * deriv; + } } } diff --git a/stan/math/prim/prob/logistic_lccdf.hpp b/stan/math/prim/prob/logistic_lccdf.hpp index 95e3bea4341..a08dc16a892 100644 --- a/stan/math/prim/prob/logistic_lccdf.hpp +++ b/stan/math/prim/prob/logistic_lccdf.hpp @@ -49,25 +49,23 @@ inline return_type_t logistic_lccdf(const T_y& y, size_t N = max_size(y, mu, sigma); // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero + // The gradients are technically ill-defined, but treated as zero. for (size_t i = 0; i < stan::math::size(y); i++) { if (y_vec.val(i) == NEGATIVE_INFTY) { return ops_partials.build(0.0); } } - - for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + for (size_t i = 0; i < stan::math::size(y); i++) { + if (y_vec.val(i) == INFTY) { return ops_partials.build(negative_infinity()); } + } + for (size_t n = 0; n < N; n++) { const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - const T_partials_return scaled_diff - = (y_dbl - mu_dbl) * sigma_inv_vec; + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; P += log1m_inv_logit(scaled_diff); if constexpr (is_any_autodiff_v) { diff --git a/stan/math/prim/prob/logistic_lcdf.hpp b/stan/math/prim/prob/logistic_lcdf.hpp index b0faa88fdc4..e0ebfd0b8e6 100644 --- a/stan/math/prim/prob/logistic_lcdf.hpp +++ b/stan/math/prim/prob/logistic_lcdf.hpp @@ -67,8 +67,7 @@ inline return_type_t logistic_lcdf(const T_y& y, const T_partials_return mu_dbl = mu_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - const T_partials_return scaled_diff - = (y_dbl - mu_dbl) * sigma_inv_vec; + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; P += log_inv_logit(scaled_diff); if constexpr (is_any_autodiff_v) { diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp index ce4bb340016..a6891a8af17 100644 --- a/stan/math/prim/prob/logistic_lpdf.hpp +++ b/stan/math/prim/prob/logistic_lpdf.hpp @@ -70,22 +70,25 @@ inline return_type_t logistic_lpdf(const T_y& y, logp -= sum(log(sigma_val)) * N / math::size(sigma); } - if constexpr (is_any_autodiff_v) { - const auto& exp_y_minus_mu_div_sigma = exp(y_minus_mu_div_sigma); - const auto& y_deriv - = to_ref_if<(is_autodiff_v && is_autodiff_v)>( - (2 / (1 + exp_y_minus_mu_div_sigma) - 1) * inv_sigma); + // d/dmu = tanh(z / 2) / sigma with z = (y - mu) / sigma. The y and sigma + // partials are built from the same expression so that d/dy == -d/dmu + // exactly; the equivalent 2 / (1 + exp(z)) - 1 form loses all relative + // precision as z -> 0 (it returns 0 for z < eps). + if constexpr (is_any_autodiff_v) { + // to_ref, not to_ref_if: tanh() of an Eigen argument returns a holder that + // owns its operand, so the product has to be evaluated inside this + // full-expression rather than forwarded on as a lazy expression. + const auto& mu_deriv = to_ref(tanh(0.5 * y_minus_mu_div_sigma) * inv_sigma); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = y_deriv; + partials<0>(ops_partials) = -mu_deriv; + } + if constexpr (is_autodiff_v) { + edge<1>(ops_partials).partials_ = mu_deriv; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = (-y_deriv * y_minus_mu - 1) * inv_sigma; + partials<2>(ops_partials) = (mu_deriv * y_minus_mu - 1) * inv_sigma; } } - if constexpr (is_autodiff_v) { - edge<1>(ops_partials).partials_ - = tanh(0.5 * y_minus_mu_div_sigma) * inv_sigma; - } return ops_partials.build(logp); } diff --git a/test/unit/math/opencl/rev/logistic_cdf_test.cpp b/test/unit/math/opencl/rev/logistic_cdf_test.cpp index f279c45691e..2def176fb54 100644 --- a/test/unit/math/opencl/rev/logistic_cdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_cdf_test.cpp @@ -97,6 +97,32 @@ TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_small_y_neg_inf) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_lower_tail) { + Eigen::VectorXd y(3); + y << -700, -350, -1400; + Eigen::VectorXd mu(3); + mu << 0, 0, 0; + Eigen::VectorXd sigma(3); + sigma << 1, 0.5, 2; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_cdf_functor, y, mu, + sigma); +} + +// y == INFTY contributes a factor of 1 and zero partials; the scale partial +// used to be 0 * INFTY = NaN on the device while the CPU skipped the element. +TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_y_pos_inf) { + Eigen::VectorXd y(3); + y << 0.3, INFINITY, 1.0; + Eigen::VectorXd mu(3); + mu << 0.3, 0.8, 1.0; + Eigen::VectorXd sigma(3); + sigma << 0.3, 0.8, 1.0; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_cdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticCdf, opencl_broadcast_y) { int N = 3; diff --git a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp index e4642438a83..25807c3febf 100644 --- a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp @@ -80,10 +80,26 @@ TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_small) { } TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_upper_tail) { - Eigen::VectorXd y(2); - y << 36.5, 40; - Eigen::VectorXd mu = Eigen::VectorXd::Zero(2); - Eigen::VectorXd sigma = Eigen::VectorXd::Ones(2); + Eigen::VectorXd y(4); + y << 30, 36.5, 40, 85; + Eigen::VectorXd mu(4); + mu << 0, 0, 0, 5; + Eigen::VectorXd sigma(4); + sigma << 1, 1, 1, 2; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, + sigma); +} + +// A single y == INFTY makes the whole result -INFTY with zero partials; the +// partials of the finite elements preceding it must not survive. +TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_y_pos_inf) { + Eigen::VectorXd y(3); + y << 1.5, INFINITY, 1.0; + Eigen::VectorXd mu(3); + mu << 0.3, 0.8, 1.0; + Eigen::VectorXd sigma(3); + sigma << 0.3, 0.8, 1.0; stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, sigma); diff --git a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp index ed812f70d39..8c52f8c2dc4 100644 --- a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp @@ -80,10 +80,24 @@ TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_small) { } TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_lower_tail) { + Eigen::VectorXd y(3); + y << -745, -746, -1495; + Eigen::VectorXd mu(3); + mu << 0, 0, 5; + Eigen::VectorXd sigma(3); + sigma << 1, 1, 2; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, + sigma); +} + +TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_y_pos_inf) { Eigen::VectorXd y(2); - y << -745, -746; - Eigen::VectorXd mu = Eigen::VectorXd::Zero(2); - Eigen::VectorXd sigma = Eigen::VectorXd::Ones(2); + y << 1.5, INFINITY; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1, 2; stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, sigma); diff --git a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp index 0c2258f6373..678eeb2381a 100644 --- a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp @@ -101,12 +101,12 @@ TEST(ProbDistributionsLogistic, opencl_matches_cpu_small) { } TEST(ProbDistributionsLogistic, opencl_matches_cpu_large_location) { - Eigen::VectorXd y(1); - y << 711; - Eigen::VectorXd mu(1); - mu << 710; - Eigen::VectorXd sigma(1); - sigma << 1; + Eigen::VectorXd y(3); + y << 711, 712, 1e-14; + Eigen::VectorXd mu(3); + mu << 710, 710, 0; + Eigen::VectorXd sigma(3); + sigma << 1, 2, 1; stan::math::test::compare_cpu_opencl_prim_rev(logistic_lpdf_functor, y, mu, sigma); diff --git a/test/unit/math/rev/prob/logistic_cdf_test.cpp b/test/unit/math/rev/prob/logistic_cdf_test.cpp new file mode 100644 index 00000000000..44adb9c0a18 --- /dev/null +++ b/test/unit/math/rev/prob/logistic_cdf_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include +#include + +// The pdf / cdf quotient underflows to 0 / 0 in the lower tail. The partials +// are now built from inv_logit(-z) and only pick up the cdf itself as a +// factor, so they stay finite all the way down to the denormal range. +TEST_F(AgradRev, logistic_cdf_lower_tail) { + stan::math::var y = -745.0; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + // dF/dy = F * (1 - F) / sigma, and 1 - F == 1 to machine precision here + EXPECT_GT(cdf.val(), 0.0); + EXPECT_DOUBLE_EQ(cdf.val(), y.adj()); + EXPECT_DOUBLE_EQ(-cdf.val(), mu.adj()); + EXPECT_DOUBLE_EQ(745.0 * cdf.val(), sigma.adj()); +} + +// One ulp further out the cdf itself underflows to zero; the gradient must +// follow it to zero rather than become NaN. +TEST_F(AgradRev, logistic_cdf_lower_tail_underflow) { + stan::math::var y = -746.0; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + EXPECT_DOUBLE_EQ(0.0, cdf.val()); + EXPECT_DOUBLE_EQ(0.0, y.adj()); + EXPECT_DOUBLE_EQ(0.0, mu.adj()); + EXPECT_DOUBLE_EQ(0.0, sigma.adj()); +} + +TEST_F(AgradRev, logistic_cdf_location_scale) { + const double mu_value = 5.0; + const double sigma_value = 2.0; + const double scaled_diff = -1.5; + + stan::math::var y = mu_value + sigma_value * scaled_diff; + stan::math::var mu = mu_value; + stan::math::var sigma = sigma_value; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + const double cdf_value = stan::math::inv_logit(scaled_diff); + const double deriv = stan::math::inv_logit(-scaled_diff) / sigma_value; + EXPECT_DOUBLE_EQ(cdf_value, cdf.val()); + EXPECT_DOUBLE_EQ(deriv * cdf_value, y.adj()); + EXPECT_DOUBLE_EQ(-deriv * cdf_value, mu.adj()); + EXPECT_DOUBLE_EQ(-scaled_diff * deriv * cdf_value, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lccdf_test.cpp b/test/unit/math/rev/prob/logistic_lccdf_test.cpp index e08c2afff17..683c90f5f1f 100644 --- a/test/unit/math/rev/prob/logistic_lccdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lccdf_test.cpp @@ -1,6 +1,7 @@ #include #include #include +#include TEST_F(AgradRev, logistic_lccdf_upper_tail) { for (double y_value : {36.5, 40.0}) { @@ -19,3 +20,62 @@ TEST_F(AgradRev, logistic_lccdf_upper_tail) { stan::math::recover_memory(); } } + +// Same tail, but with mu != 0 and sigma != 1 so that the scaled difference is +// distinguishable from y and the 1 / sigma factors in the partials are +// exercised. +TEST_F(AgradRev, logistic_lccdf_upper_tail_location_scale) { + const double mu_value = 5.0; + const double sigma_value = 2.0; + + for (double scaled_diff : {36.5, 40.0}) { + stan::math::var y = mu_value + sigma_value * scaled_diff; + stan::math::var mu = mu_value; + stan::math::var sigma = sigma_value; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + const double deriv = stan::math::inv_logit(scaled_diff) / sigma_value; + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(scaled_diff), log_ccdf.val()); + EXPECT_DOUBLE_EQ(-deriv, y.adj()); + EXPECT_DOUBLE_EQ(deriv, mu.adj()); + EXPECT_DOUBLE_EQ(scaled_diff * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} + +// The 1 - inv_logit(z) cancellation degrades long before it returns -Inf: +// at z = 30 it gave -30.001021 for the value and -1.00102 for the y partial. +TEST_F(AgradRev, logistic_lccdf_upper_tail_moderate) { + stan::math::var y = 30.0; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + EXPECT_NEAR(-30.000000000000092, log_ccdf.val(), 1e-12); + EXPECT_NEAR(-0.99999999999990652, y.adj(), 1e-12); + EXPECT_NEAR(0.99999999999990652, mu.adj(), 1e-12); + EXPECT_NEAR(29.999999999997197, sigma.adj(), 1e-10); +} + +// An infinite element short-circuits the result; the partials of the finite +// elements that precede it must not leak into the returned gradient. +TEST_F(AgradRev, logistic_lccdf_pos_inf_zeroes_partials) { + Eigen::Matrix y(2); + y << 1.5, stan::math::INFTY; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + EXPECT_TRUE(std::isinf(log_ccdf.val())); + EXPECT_LT(log_ccdf.val(), 0.0); + EXPECT_DOUBLE_EQ(0.0, y(0).adj()); + EXPECT_DOUBLE_EQ(0.0, y(1).adj()); + EXPECT_DOUBLE_EQ(0.0, mu.adj()); + EXPECT_DOUBLE_EQ(0.0, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lcdf_test.cpp b/test/unit/math/rev/prob/logistic_lcdf_test.cpp index 40146563763..9a645e78c03 100644 --- a/test/unit/math/rev/prob/logistic_lcdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lcdf_test.cpp @@ -19,3 +19,46 @@ TEST_F(AgradRev, logistic_lcdf_lower_tail) { stan::math::recover_memory(); } } + +// Same tail, but with mu != 0 and sigma != 1 so that the scaled difference is +// distinguishable from y and the 1 / sigma factors in the partials are +// exercised. +TEST_F(AgradRev, logistic_lcdf_lower_tail_location_scale) { + const double mu_value = 5.0; + const double sigma_value = 2.0; + + for (double scaled_diff : {-745.0, -750.0}) { + stan::math::var y = mu_value + sigma_value * scaled_diff; + stan::math::var mu = mu_value; + stan::math::var sigma = sigma_value; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = stan::math::inv_logit(-scaled_diff) / sigma_value; + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(-scaled_diff), log_cdf.val()); + EXPECT_DOUBLE_EQ(deriv, y.adj()); + EXPECT_DOUBLE_EQ(-deriv, mu.adj()); + EXPECT_DOUBLE_EQ(-scaled_diff * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} + +// y == INFTY contributes log(1) = 0 and zero partials; it must not poison the +// gradient of the finite elements. +TEST_F(AgradRev, logistic_lcdf_pos_inf_element) { + Eigen::Matrix y(2); + y << 1.5, stan::math::INFTY; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = stan::math::inv_logit(-1.5); + EXPECT_DOUBLE_EQ(stan::math::log_inv_logit(1.5), log_cdf.val()); + EXPECT_DOUBLE_EQ(deriv, y(0).adj()); + EXPECT_DOUBLE_EQ(0.0, y(1).adj()); + EXPECT_DOUBLE_EQ(-deriv, mu.adj()); + EXPECT_DOUBLE_EQ(-1.5 * deriv, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lpdf_test.cpp b/test/unit/math/rev/prob/logistic_lpdf_test.cpp index 67ad7ad8a1e..28882ac4d87 100644 --- a/test/unit/math/rev/prob/logistic_lpdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lpdf_test.cpp @@ -4,11 +4,58 @@ #include TEST_F(AgradRev, logistic_lpdf_location_gradient_large_values) { + stan::math::var y = 711.0; stan::math::var mu = 710.0; + stan::math::var sigma = 1.0; - stan::math::var logp = stan::math::logistic_lpdf(711.0, mu, 1.0); + stan::math::var logp = stan::math::logistic_lpdf(y, mu, sigma); logp.grad(); + // z = 1: d/dmu = tanh(z / 2) / sigma, d/dy = -d/dmu, + // d/dsigma = (z * tanh(z / 2) - 1) / sigma + const double tanh_half = std::tanh(0.5); + EXPECT_DOUBLE_EQ(-1.0 - 2.0 * stan::math::log1p_exp(-1.0), logp.val()); + EXPECT_DOUBLE_EQ(tanh_half, mu.adj()); + EXPECT_DOUBLE_EQ(-tanh_half, y.adj()); + EXPECT_DOUBLE_EQ(tanh_half - 1.0, sigma.adj()); +} + +// The 2 / (1 + exp(z)) - 1 form of the y partial returns 0 for z below eps, +// while the location partial is computed from tanh(z / 2). The two must agree +// to the last bit: they are the same quantity with opposite sign. +TEST_F(AgradRev, logistic_lpdf_gradients_near_location) { + for (double scaled_diff : {1e-8, 1e-12, 1e-14, 1e-16}) { + stan::math::var y = scaled_diff; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var logp = stan::math::logistic_lpdf(y, mu, sigma); + logp.grad(); + + EXPECT_DOUBLE_EQ(std::tanh(0.5 * scaled_diff), mu.adj()); + EXPECT_DOUBLE_EQ(-mu.adj(), y.adj()); + stan::math::recover_memory(); + } +} + +// The location partial is assigned through the container path of the edge, so +// exercise it with a vector argument and a non-unit scale. +TEST_F(AgradRev, logistic_lpdf_gradients_vectorized) { + Eigen::Matrix y(2); + y << 711.0, 712.0; + stan::math::var mu = 710.0; + stan::math::var sigma = 2.0; + + stan::math::var logp = stan::math::logistic_lpdf(y, mu, sigma); + logp.grad(); + + const double d0 = std::tanh(0.25) / 2.0; // z = 0.5 + const double d1 = std::tanh(0.5) / 2.0; // z = 1.0 EXPECT_TRUE(std::isfinite(logp.val())); - EXPECT_DOUBLE_EQ(std::tanh(0.5), mu.adj()); + EXPECT_DOUBLE_EQ(d0 + d1, mu.adj()); + EXPECT_DOUBLE_EQ(-d0, y(0).adj()); + EXPECT_DOUBLE_EQ(-d1, y(1).adj()); + EXPECT_DOUBLE_EQ( + (0.5 * std::tanh(0.25) - 1.0) / 2.0 + (1.0 * std::tanh(0.5) - 1.0) / 2.0, + sigma.adj()); } From 4b0034b6e6f9c5fd9ff6f870d122ad3ab15e4520 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 16 Sep 2026 18:19:30 +0300 Subject: [PATCH 06/10] Cut comments Co-authored-by: Andrew Johnson --- stan/math/prim/prob/logistic_cdf.hpp | 3 --- stan/math/prim/prob/logistic_lpdf.hpp | 5 ----- 2 files changed, 8 deletions(-) diff --git a/stan/math/prim/prob/logistic_cdf.hpp b/stan/math/prim/prob/logistic_cdf.hpp index fe659f91333..6ae1a308eb2 100644 --- a/stan/math/prim/prob/logistic_cdf.hpp +++ b/stan/math/prim/prob/logistic_cdf.hpp @@ -71,9 +71,6 @@ inline return_type_t logistic_cdf(const T_y& y, P *= Pn; - // The partials accumulate d/d. log(Pn); they are rescaled by the product - // P below. Writing the log-scale derivative as inv_logit(-scaled_diff) - // avoids the pdf / Pn quotient, which is 0 / 0 once Pn underflows. if constexpr (is_any_autodiff_v) { const T_partials_return deriv = inv_logit(-scaled_diff) * sigma_inv_vec; if constexpr (is_autodiff_v) { diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp index a6891a8af17..027daa3b5af 100644 --- a/stan/math/prim/prob/logistic_lpdf.hpp +++ b/stan/math/prim/prob/logistic_lpdf.hpp @@ -70,14 +70,9 @@ inline return_type_t logistic_lpdf(const T_y& y, logp -= sum(log(sigma_val)) * N / math::size(sigma); } - // d/dmu = tanh(z / 2) / sigma with z = (y - mu) / sigma. The y and sigma - // partials are built from the same expression so that d/dy == -d/dmu - // exactly; the equivalent 2 / (1 + exp(z)) - 1 form loses all relative - // precision as z -> 0 (it returns 0 for z < eps). if constexpr (is_any_autodiff_v) { // to_ref, not to_ref_if: tanh() of an Eigen argument returns a holder that // owns its operand, so the product has to be evaluated inside this - // full-expression rather than forwarded on as a lazy expression. const auto& mu_deriv = to_ref(tanh(0.5 * y_minus_mu_div_sigma) * inv_sigma); if constexpr (is_autodiff_v) { partials<0>(ops_partials) = -mu_deriv; From 6f54a4c4d2452af91bc2dd01977b0e74c21b446e Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 16 Sep 2026 18:32:12 +0300 Subject: [PATCH 07/10] Keep the logistic tail derivative on the log scale where inv_logit underflows --- stan/math/prim/prob/logistic_cdf.hpp | 9 ++++++++- stan/math/prim/prob/logistic_lccdf.hpp | 9 ++++++++- stan/math/prim/prob/logistic_lcdf.hpp | 9 ++++++++- test/unit/math/rev/prob/logistic_cdf_test.cpp | 18 ++++++++++++++++++ .../unit/math/rev/prob/logistic_lccdf_test.cpp | 18 ++++++++++++++++++ test/unit/math/rev/prob/logistic_lcdf_test.cpp | 18 ++++++++++++++++++ 6 files changed, 78 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/prob/logistic_cdf.hpp b/stan/math/prim/prob/logistic_cdf.hpp index 6ae1a308eb2..3e162b24eaa 100644 --- a/stan/math/prim/prob/logistic_cdf.hpp +++ b/stan/math/prim/prob/logistic_cdf.hpp @@ -4,11 +4,15 @@ #include #include #include +#include +#include +#include #include #include #include #include #include +#include #include #include @@ -72,7 +76,10 @@ inline return_type_t logistic_cdf(const T_y& y, P *= Pn; if constexpr (is_any_autodiff_v) { - const T_partials_return deriv = inv_logit(-scaled_diff) * sigma_inv_vec; + const T_partials_return deriv + = value_of_rec(scaled_diff) > 700.0 + ? exp(log1m_inv_logit(scaled_diff) - log(sigma_vec.val(n))) + : inv_logit(-scaled_diff) * sigma_inv_vec; if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] += deriv; } diff --git a/stan/math/prim/prob/logistic_lccdf.hpp b/stan/math/prim/prob/logistic_lccdf.hpp index a08dc16a892..2bc889f27ac 100644 --- a/stan/math/prim/prob/logistic_lccdf.hpp +++ b/stan/math/prim/prob/logistic_lccdf.hpp @@ -4,13 +4,17 @@ #include #include #include +#include #include +#include #include +#include #include #include #include #include #include +#include #include namespace stan { @@ -69,7 +73,10 @@ inline return_type_t logistic_lccdf(const T_y& y, P += log1m_inv_logit(scaled_diff); if constexpr (is_any_autodiff_v) { - const T_partials_return deriv = inv_logit(scaled_diff) * sigma_inv_vec; + const T_partials_return deriv + = value_of_rec(scaled_diff) < -700.0 + ? exp(log_inv_logit(scaled_diff) - log(sigma_vec.val(n))) + : inv_logit(scaled_diff) * sigma_inv_vec; if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] -= deriv; } diff --git a/stan/math/prim/prob/logistic_lcdf.hpp b/stan/math/prim/prob/logistic_lcdf.hpp index e0ebfd0b8e6..ee64c6bd021 100644 --- a/stan/math/prim/prob/logistic_lcdf.hpp +++ b/stan/math/prim/prob/logistic_lcdf.hpp @@ -4,13 +4,17 @@ #include #include #include +#include #include +#include +#include #include #include #include #include #include #include +#include #include namespace stan { @@ -71,7 +75,10 @@ inline return_type_t logistic_lcdf(const T_y& y, P += log_inv_logit(scaled_diff); if constexpr (is_any_autodiff_v) { - const T_partials_return deriv = inv_logit(-scaled_diff) * sigma_inv_vec; + const T_partials_return deriv + = value_of_rec(scaled_diff) > 700.0 + ? exp(log1m_inv_logit(scaled_diff) - log(sigma_vec.val(n))) + : inv_logit(-scaled_diff) * sigma_inv_vec; if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] += deriv; } diff --git a/test/unit/math/rev/prob/logistic_cdf_test.cpp b/test/unit/math/rev/prob/logistic_cdf_test.cpp index 44adb9c0a18..ad1b424897f 100644 --- a/test/unit/math/rev/prob/logistic_cdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_cdf_test.cpp @@ -37,6 +37,24 @@ TEST_F(AgradRev, logistic_cdf_lower_tail_underflow) { EXPECT_DOUBLE_EQ(0.0, sigma.adj()); } +// inv_logit(-z) underflows to zero above z = 745; with the cdf itself equal +// to one there, the gradient is the only thing left to get right. +TEST_F(AgradRev, logistic_cdf_underflow_rescued_by_small_sigma) { + stan::math::var y = 8e-298; + stan::math::var mu = 0.0; + stan::math::var sigma = 1e-300; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + const double deriv = 3.6678745841780173e-48; + EXPECT_EQ(1.0, cdf.val()); + EXPECT_NEAR(deriv, y.adj(), 1e-10 * deriv); + EXPECT_NEAR(-deriv, mu.adj(), 1e-10 * deriv); + EXPECT_NEAR(-2.9342996673424137e-45, sigma.adj(), + 1e-10 * 2.9342996673424137e-45); +} + TEST_F(AgradRev, logistic_cdf_location_scale) { const double mu_value = 5.0; const double sigma_value = 2.0; diff --git a/test/unit/math/rev/prob/logistic_lccdf_test.cpp b/test/unit/math/rev/prob/logistic_lccdf_test.cpp index 683c90f5f1f..b8048359a56 100644 --- a/test/unit/math/rev/prob/logistic_lccdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lccdf_test.cpp @@ -61,6 +61,24 @@ TEST_F(AgradRev, logistic_lccdf_upper_tail_moderate) { EXPECT_NEAR(29.999999999997197, sigma.adj(), 1e-10); } +// inv_logit(z) underflows to zero below z = -745, but dividing by a small +// enough sigma brings the quotient back into range. +TEST_F(AgradRev, logistic_lccdf_underflow_rescued_by_small_sigma) { + stan::math::var y = -8e-298; + stan::math::var mu = 0.0; + stan::math::var sigma = 1e-300; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + const double deriv = 3.6678745841780173e-48; + EXPECT_EQ(0.0, stan::math::inv_logit(-8e-298 / 1e-300)); + EXPECT_NEAR(-deriv, y.adj(), 1e-10 * deriv); + EXPECT_NEAR(deriv, mu.adj(), 1e-10 * deriv); + EXPECT_NEAR(-2.9342996673424137e-45, sigma.adj(), + 1e-10 * 2.9342996673424137e-45); +} + // An infinite element short-circuits the result; the partials of the finite // elements that precede it must not leak into the returned gradient. TEST_F(AgradRev, logistic_lccdf_pos_inf_zeroes_partials) { diff --git a/test/unit/math/rev/prob/logistic_lcdf_test.cpp b/test/unit/math/rev/prob/logistic_lcdf_test.cpp index 9a645e78c03..f0936bcc2d5 100644 --- a/test/unit/math/rev/prob/logistic_lcdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lcdf_test.cpp @@ -44,6 +44,24 @@ TEST_F(AgradRev, logistic_lcdf_lower_tail_location_scale) { } } +// inv_logit(-z) underflows to zero above z = 745, but dividing by a small +// enough sigma brings the quotient back into range. +TEST_F(AgradRev, logistic_lcdf_underflow_rescued_by_small_sigma) { + stan::math::var y = 8e-298; + stan::math::var mu = 0.0; + stan::math::var sigma = 1e-300; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = 3.6678745841780173e-48; + EXPECT_EQ(0.0, stan::math::inv_logit(-8e-298 / 1e-300)); + EXPECT_NEAR(deriv, y.adj(), 1e-10 * deriv); + EXPECT_NEAR(-deriv, mu.adj(), 1e-10 * deriv); + EXPECT_NEAR(-2.9342996673424137e-45, sigma.adj(), + 1e-10 * 2.9342996673424137e-45); +} + // y == INFTY contributes log(1) = 0 and zero partials; it must not poison the // gradient of the finite elements. TEST_F(AgradRev, logistic_lcdf_pos_inf_element) { From 448dcc178c12d563d8284892df446afbc1008e02 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 16 Sep 2026 18:56:35 +0300 Subject: [PATCH 08/10] Guard the OpenCL logistic tail derivative the same way as the CPU --- .../kernel_generator/elt_function_cl.hpp | 6 +++ .../device_functions/logistic_tail_deriv.hpp | 43 +++++++++++++++++++ stan/math/opencl/prim/logistic_cdf.hpp | 2 +- stan/math/opencl/prim/logistic_lccdf.hpp | 2 +- stan/math/opencl/prim/logistic_lcdf.hpp | 2 +- .../math/opencl/rev/logistic_cdf_test.cpp | 12 ++++++ .../math/opencl/rev/logistic_lccdf_test.cpp | 12 ++++++ .../math/opencl/rev/logistic_lcdf_test.cpp | 12 ++++++ 8 files changed, 88 insertions(+), 3 deletions(-) create mode 100644 stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index 96bbabbaf65..cccc7401638 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -342,6 +343,11 @@ ADD_CLASSIFICATION_FUNCTION(isinf, ADD_CLASSIFICATION_FUNCTION(isnan, this->template get_arg<0>().extreme_diagonals()) +ADD_BINARY_FUNCTION_WITH_INCLUDES( + logistic_tail_deriv, opencl_kernels::log1p_exp_device_function, + opencl_kernels::log1m_inv_logit_device_function, + opencl_kernels::inv_logit_device_function, + opencl_kernels::logistic_tail_deriv_device_function) ADD_BINARY_FUNCTION_WITH_INCLUDES(fdim) ADD_BINARY_FUNCTION_WITH_INCLUDES(fmax) ADD_BINARY_FUNCTION_WITH_INCLUDES(fmin) diff --git a/stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp b/stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp new file mode 100644 index 00000000000..87f5ca52d90 --- /dev/null +++ b/stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp @@ -0,0 +1,43 @@ +#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_DERIV_HPP +#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_DERIV_HPP +#ifdef STAN_OPENCL + +#include +#include + +namespace stan { +namespace math { +namespace opencl_kernels { + +// \cond +static constexpr const char* logistic_tail_deriv_device_function + = "\n" + "#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_DERIV\n" + "#define " + "STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_" + "DERIV\n" STRINGIFY( + // \endcond + /** \ingroup opencl_kernels + * + * Return inv_logit(-x) / sigma. + * + * @param x scaled difference + * @param sigma scale + * @return inv_logit(-x) / sigma + */ + double logistic_tail_deriv(double x, double sigma) { + if (x > 700.0) { + return exp(log1m_inv_logit(x) - log(sigma)); + } + return inv_logit(-x) / sigma; + } + // \cond + ) "\n#endif\n"; // NOLINT +// \endcond + +} // namespace opencl_kernels +} // namespace math +} // namespace stan + +#endif +#endif diff --git a/stan/math/opencl/prim/logistic_cdf.hpp b/stan/math/opencl/prim/logistic_cdf.hpp index 307f17460fd..336191c827d 100644 --- a/stan/math/opencl/prim/logistic_cdf.hpp +++ b/stan/math/opencl/prim/logistic_cdf.hpp @@ -73,7 +73,7 @@ inline return_type_t logistic_cdf( // These are the log-scale derivatives; they are rescaled by the product P // below. inv_logit(-scaled_diff) avoids the pdf / Pn quotient, which is // 0 / 0 once Pn underflows. - auto deriv = elt_multiply(inv_logit(-scaled_diff), inv_sigma); + auto deriv = logistic_tail_deriv(scaled_diff, sigma_val); auto y_deriv_tmp = select(cond, 0.0, deriv); auto sigma_deriv_tmp = select(cond, 0.0, elt_multiply(-deriv, scaled_diff)); diff --git a/stan/math/opencl/prim/logistic_lccdf.hpp b/stan/math/opencl/prim/logistic_lccdf.hpp index 4552d58bbb7..cfe8b1c871b 100644 --- a/stan/math/opencl/prim/logistic_lccdf.hpp +++ b/stan/math/opencl/prim/logistic_lccdf.hpp @@ -69,7 +69,7 @@ inline return_type_t logistic_lccdf( auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); auto P_expr = colwise_sum(log1m_inv_logit(scaled_diff)); - auto mu_deriv = elt_multiply(inv_logit(scaled_diff), inv_sigma); + auto mu_deriv = logistic_tail_deriv(-scaled_diff, sigma_val); auto y_deriv = -mu_deriv; auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff); diff --git a/stan/math/opencl/prim/logistic_lcdf.hpp b/stan/math/opencl/prim/logistic_lcdf.hpp index b9706fe715b..841da7ae211 100644 --- a/stan/math/opencl/prim/logistic_lcdf.hpp +++ b/stan/math/opencl/prim/logistic_lcdf.hpp @@ -72,7 +72,7 @@ inline return_type_t logistic_lcdf( // y == INFTY contributes log(1) = 0 to P and zero to every partial; without // the select the scale partial would be 0 * INFTY = NaN, which prim (where // the element is skipped outright) never produces. - auto deriv = elt_multiply(inv_logit(-scaled_diff), inv_sigma); + auto deriv = logistic_tail_deriv(scaled_diff, sigma_val); auto y_deriv = select(cond, 0.0, deriv); auto mu_deriv = -y_deriv; auto sigma_deriv = select(cond, 0.0, elt_multiply(-deriv, scaled_diff)); diff --git a/test/unit/math/opencl/rev/logistic_cdf_test.cpp b/test/unit/math/opencl/rev/logistic_cdf_test.cpp index 2def176fb54..27d7c64912a 100644 --- a/test/unit/math/opencl/rev/logistic_cdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_cdf_test.cpp @@ -111,6 +111,18 @@ TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_lower_tail) { // y == INFTY contributes a factor of 1 and zero partials; the scale partial // used to be 0 * INFTY = NaN on the device while the CPU skipped the element. +TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_underflow_small_sigma) { + Eigen::VectorXd y(2); + y << 8e-298, 1.0; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_cdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_y_pos_inf) { Eigen::VectorXd y(3); y << 0.3, INFINITY, 1.0; diff --git a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp index 25807c3febf..a425e65e843 100644 --- a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp @@ -93,6 +93,18 @@ TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_upper_tail) { // A single y == INFTY makes the whole result -INFTY with zero partials; the // partials of the finite elements preceding it must not survive. +TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_underflow_small_sigma) { + Eigen::VectorXd y(2); + y << -8e-298, -1.0; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_y_pos_inf) { Eigen::VectorXd y(3); y << 1.5, INFINITY, 1.0; diff --git a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp index 8c52f8c2dc4..67d25a4f78c 100644 --- a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp @@ -91,6 +91,18 @@ TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_lower_tail) { sigma); } +TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_underflow_small_sigma) { + Eigen::VectorXd y(2); + y << 8e-298, 1.0; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_y_pos_inf) { Eigen::VectorXd y(2); y << 1.5, INFINITY; From b63ac78b5942ba62430edd67f4df08d7facc59dc Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 12 Oct 2026 23:31:41 +0300 Subject: [PATCH 09/10] Fix logistic edge cases found in review --- stan/math/opencl/prim/logistic_cdf.hpp | 5 +- stan/math/opencl/prim/logistic_lccdf.hpp | 23 +++---- stan/math/opencl/prim/logistic_lcdf.hpp | 8 +-- stan/math/opencl/prim/logistic_lpdf.hpp | 4 +- stan/math/prim/prob/logistic_cdf.hpp | 6 +- stan/math/prim/prob/logistic_lccdf.hpp | 16 +++-- stan/math/prim/prob/logistic_lcdf.hpp | 5 +- stan/math/prim/prob/logistic_lpdf.hpp | 8 ++- .../math/opencl/rev/logistic_cdf_test.cpp | 51 +++++++++++++++ .../math/opencl/rev/logistic_lccdf_test.cpp | 62 +++++++++++++++++- .../math/opencl/rev/logistic_lcdf_test.cpp | 40 ++++++++++++ .../math/opencl/rev/logistic_lpdf_test.cpp | 28 +++++++++ test/unit/math/rev/prob/logistic_cdf_test.cpp | 47 ++++++++++++++ .../math/rev/prob/logistic_lccdf_test.cpp | 63 +++++++++++++++++++ .../unit/math/rev/prob/logistic_lcdf_test.cpp | 42 +++++++++++++ .../unit/math/rev/prob/logistic_lpdf_test.cpp | 25 ++++++++ 16 files changed, 404 insertions(+), 29 deletions(-) diff --git a/stan/math/opencl/prim/logistic_cdf.hpp b/stan/math/opencl/prim/logistic_cdf.hpp index 336191c827d..4798f97e24d 100644 --- a/stan/math/opencl/prim/logistic_cdf.hpp +++ b/stan/math/opencl/prim/logistic_cdf.hpp @@ -73,9 +73,12 @@ inline return_type_t logistic_cdf( // These are the log-scale derivatives; they are rescaled by the product P // below. inv_logit(-scaled_diff) avoids the pdf / Pn quotient, which is // 0 / 0 once Pn underflows. + // At an infinite scaled_diff, either deriv or P is zero; the scale partial + // would be 0 * INFTY = NaN without the select. auto deriv = logistic_tail_deriv(scaled_diff, sigma_val); auto y_deriv_tmp = select(cond, 0.0, deriv); - auto sigma_deriv_tmp = select(cond, 0.0, elt_multiply(-deriv, scaled_diff)); + auto sigma_deriv_tmp + = select(isinf(scaled_diff), 0.0, elt_multiply(-deriv, scaled_diff)); matrix_cl any_y_neg_inf_cl; matrix_cl P_cl; diff --git a/stan/math/opencl/prim/logistic_lccdf.hpp b/stan/math/opencl/prim/logistic_lccdf.hpp index cfe8b1c871b..b6fa69bdb3e 100644 --- a/stan/math/opencl/prim/logistic_lccdf.hpp +++ b/stan/math/opencl/prim/logistic_lccdf.hpp @@ -63,17 +63,22 @@ inline return_type_t logistic_lccdf( = check_cl(function, "Scale parameter", sigma_val, "positive finite"); auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); - auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto any_y_pos_inf = colwise_max(cast(y_val == INFTY)); + auto y_neg_inf = y_val == NEGATIVE_INFTY; auto inv_sigma = elt_divide(1.0, sigma_val); auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); - auto P_expr = colwise_sum(log1m_inv_logit(scaled_diff)); + // y == NEGATIVE_INFTY contributes log(1) = 0 to P and zero to every partial + auto P_expr + = colwise_sum(select(y_neg_inf, 0.0, log1m_inv_logit(scaled_diff))); - auto mu_deriv = logistic_tail_deriv(-scaled_diff, sigma_val); + auto mu_deriv + = select(y_neg_inf, 0.0, logistic_tail_deriv(-scaled_diff, sigma_val)); auto y_deriv = -mu_deriv; - auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff); + // scaled_diff == NEGATIVE_INFTY (also from an overflow with finite y) has a + // zero mu partial, so the product would be 0 * INFTY = NaN + auto sigma_deriv = select(scaled_diff == NEGATIVE_INFTY, 0.0, + elt_multiply(mu_deriv, scaled_diff)); - matrix_cl any_y_neg_inf_cl; matrix_cl any_y_pos_inf_cl; matrix_cl P_cl; matrix_cl mu_deriv_cl; @@ -81,17 +86,13 @@ inline return_type_t logistic_lccdf( matrix_cl sigma_deriv_cl; results(check_y_not_nan, check_mu_finite, check_sigma_positive_finite, - any_y_neg_inf_cl, any_y_pos_inf_cl, P_cl, y_deriv_cl, mu_deriv_cl, - sigma_deriv_cl) + any_y_pos_inf_cl, P_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl) = expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_finite_expr, - any_y_neg_inf, any_y_pos_inf, P_expr, + any_y_pos_inf, P_expr, calc_if>(y_deriv), calc_if>(mu_deriv), calc_if>(sigma_deriv)); - if (from_matrix_cl(any_y_neg_inf_cl).maxCoeff()) { - return 0.0; - } if (from_matrix_cl(any_y_pos_inf_cl).maxCoeff()) { return NEGATIVE_INFTY; } diff --git a/stan/math/opencl/prim/logistic_lcdf.hpp b/stan/math/opencl/prim/logistic_lcdf.hpp index 841da7ae211..31d949630fd 100644 --- a/stan/math/opencl/prim/logistic_lcdf.hpp +++ b/stan/math/opencl/prim/logistic_lcdf.hpp @@ -64,14 +64,14 @@ inline return_type_t logistic_lcdf( auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val); auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); - auto cond = y_val == INFTY; auto inv_sigma = elt_divide(1.0, sigma_val); auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); + auto cond = scaled_diff == INFTY; auto P_expr = colwise_sum(log_inv_logit(scaled_diff)); - // y == INFTY contributes log(1) = 0 to P and zero to every partial; without - // the select the scale partial would be 0 * INFTY = NaN, which prim (where - // the element is skipped outright) never produces. + // scaled_diff == INFTY (y == INFTY, or an overflow with finite y) + // contributes log(1) = 0 to P and zero to every partial; without the select + // the scale partial would be 0 * INFTY = NaN. auto deriv = logistic_tail_deriv(scaled_diff, sigma_val); auto y_deriv = select(cond, 0.0, deriv); auto mu_deriv = -y_deriv; diff --git a/stan/math/opencl/prim/logistic_lpdf.hpp b/stan/math/opencl/prim/logistic_lpdf.hpp index 3ef004d2df9..2f333f528c7 100644 --- a/stan/math/opencl/prim/logistic_lpdf.hpp +++ b/stan/math/opencl/prim/logistic_lpdf.hpp @@ -76,7 +76,9 @@ inline return_type_t logistic_lpdf( auto y_minus_mu = y_val - mu_val; auto y_minus_mu_div_sigma = elt_multiply(y_minus_mu, inv_sigma); - auto logp1 = -y_minus_mu_div_sigma - 2.0 * log1p_exp(-y_minus_mu_div_sigma); + // The density is symmetric in z; -|z| avoids inf - inf at z = -inf. + auto abs_scaled_diff = fabs(y_minus_mu_div_sigma); + auto logp1 = -abs_scaled_diff - 2.0 * log1p_exp(-abs_scaled_diff); auto logp_expr = colwise_sum(static_select::value>( logp1 - log(sigma_val), logp1)); diff --git a/stan/math/prim/prob/logistic_cdf.hpp b/stan/math/prim/prob/logistic_cdf.hpp index 3e162b24eaa..0a0dbddf2fb 100644 --- a/stan/math/prim/prob/logistic_cdf.hpp +++ b/stan/math/prim/prob/logistic_cdf.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -87,7 +88,10 @@ inline return_type_t logistic_cdf(const T_y& y, partials<1>(ops_partials)[n] -= deriv; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] -= scaled_diff * deriv; + // At an overflowed scaled_diff, either deriv or P is zero + if (!is_inf(scaled_diff)) { + partials<2>(ops_partials)[n] -= scaled_diff * deriv; + } } } } diff --git a/stan/math/prim/prob/logistic_lccdf.hpp b/stan/math/prim/prob/logistic_lccdf.hpp index 2bc889f27ac..8b01cbbefd8 100644 --- a/stan/math/prim/prob/logistic_lccdf.hpp +++ b/stan/math/prim/prob/logistic_lccdf.hpp @@ -54,11 +54,6 @@ inline return_type_t logistic_lccdf(const T_y& y, // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero. - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == NEGATIVE_INFTY) { - return ops_partials.build(0.0); - } - } for (size_t i = 0; i < stan::math::size(y); i++) { if (y_vec.val(i) == INFTY) { return ops_partials.build(negative_infinity()); @@ -66,6 +61,12 @@ inline return_type_t logistic_lccdf(const T_y& y, } for (size_t n = 0; n < N; n++) { + // Explicit results for extreme values + // The gradients are technically ill-defined, but treated as zero + if (y_vec.val(n) == NEGATIVE_INFTY) { + continue; + } + const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); @@ -84,7 +85,10 @@ inline return_type_t logistic_lccdf(const T_y& y, partials<1>(ops_partials)[n] += deriv; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] += scaled_diff * deriv; + // deriv is zero at an overflowed scaled_diff == NEGATIVE_INFTY + if (scaled_diff != NEGATIVE_INFTY) { + partials<2>(ops_partials)[n] += scaled_diff * deriv; + } } } } diff --git a/stan/math/prim/prob/logistic_lcdf.hpp b/stan/math/prim/prob/logistic_lcdf.hpp index ee64c6bd021..d12dabca17e 100644 --- a/stan/math/prim/prob/logistic_lcdf.hpp +++ b/stan/math/prim/prob/logistic_lcdf.hpp @@ -86,7 +86,10 @@ inline return_type_t logistic_lcdf(const T_y& y, partials<1>(ops_partials)[n] -= deriv; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] -= scaled_diff * deriv; + // deriv is zero at an overflowed scaled_diff == INFTY + if (scaled_diff != INFTY) { + partials<2>(ops_partials)[n] -= scaled_diff * deriv; + } } } } diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp index 027daa3b5af..e49ed38f2a3 100644 --- a/stan/math/prim/prob/logistic_lpdf.hpp +++ b/stan/math/prim/prob/logistic_lpdf.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -64,8 +65,11 @@ inline return_type_t logistic_lpdf(const T_y& y, const auto& y_minus_mu_div_sigma = to_ref(y_minus_mu * inv_sigma); size_t N = max_size(y, mu, sigma); - T_partials_return logp = -sum(y_minus_mu_div_sigma) - - 2.0 * sum(log1p_exp(-y_minus_mu_div_sigma)); + // The density is symmetric in the scaled difference z. The -|z| form is + // -inf at an overflowed z = -inf, where -z - 2 log1p_exp(-z) is inf - inf. + const auto& abs_scaled_diff = to_ref(fabs(y_minus_mu_div_sigma)); + T_partials_return logp + = -sum(abs_scaled_diff) - 2.0 * sum(log1p_exp(-abs_scaled_diff)); if constexpr (include_summand::value) { logp -= sum(log(sigma_val)) * N / math::size(sigma); } diff --git a/test/unit/math/opencl/rev/logistic_cdf_test.cpp b/test/unit/math/opencl/rev/logistic_cdf_test.cpp index 27d7c64912a..14beab5a014 100644 --- a/test/unit/math/opencl/rev/logistic_cdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_cdf_test.cpp @@ -135,6 +135,57 @@ TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_y_pos_inf) { sigma); } +// A scaled difference that overflows to INFTY with finite arguments +// (y = 1e10, sigma = 1e-300) contributes a factor 1 and zero partials; one +// that overflows to NEGATIVE_INFTY makes the cdf and every partial 0. +// Absolute reference values from mpmath at z = 0.5, sigma = 2, because a +// comparison with the CPU passes when both are wrong in the same way. +TEST(ProbDistributionsLogisticCdf, opencl_overflow_reference) { + Eigen::VectorXd mu(2); + mu << 0.0, 0.5; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 2.0; + + for (double y0 : {1e10, -1e10}) { + Eigen::VectorXd y(2); + y << y0, 1.5; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_cdf_functor, y, mu, + sigma); + + stan::math::var_value> y_cl + = stan::math::to_matrix_cl(y); + stan::math::var_value> mu_cl + = stan::math::to_matrix_cl(mu); + stan::math::var_value> sigma_cl + = stan::math::to_matrix_cl(sigma); + stan::math::var cdf = stan::math::logistic_cdf(y_cl, mu_cl, sigma_cl); + cdf.grad(); + Eigen::VectorXd y_adj + = stan::math::from_matrix_cl(y_cl.adj()); + Eigen::VectorXd mu_adj + = stan::math::from_matrix_cl(mu_cl.adj()); + Eigen::VectorXd sigma_adj + = stan::math::from_matrix_cl(sigma_cl.adj()); + + EXPECT_EQ(0.0, y_adj(0)); + EXPECT_EQ(0.0, mu_adj(0)); + EXPECT_EQ(0.0, sigma_adj(0)); + if (y0 > 0) { + EXPECT_NEAR(0.62245933120185456, cdf.val(), 1e-15); + EXPECT_NEAR(0.11750185610079724, y_adj(1), 1e-15); + EXPECT_NEAR(-0.11750185610079724, mu_adj(1), 1e-15); + EXPECT_NEAR(-0.058750928050398622, sigma_adj(1), 1e-15); + } else { + EXPECT_EQ(0.0, cdf.val()); + EXPECT_EQ(0.0, y_adj(1)); + EXPECT_EQ(0.0, mu_adj(1)); + EXPECT_EQ(0.0, sigma_adj(1)); + } + stan::math::recover_memory(); + } +} + TEST(ProbDistributionsLogisticCdf, opencl_broadcast_y) { int N = 3; diff --git a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp index a425e65e843..d71c84a6a4b 100644 --- a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp @@ -91,8 +91,6 @@ TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_upper_tail) { sigma); } -// A single y == INFTY makes the whole result -INFTY with zero partials; the -// partials of the finite elements preceding it must not survive. TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_underflow_small_sigma) { Eigen::VectorXd y(2); y << -8e-298, -1.0; @@ -105,6 +103,8 @@ TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_underflow_small_sigma) { sigma); } +// A single y == INFTY makes the whole result -INFTY with zero partials; the +// partials of the finite elements preceding it must not survive. TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_y_pos_inf) { Eigen::VectorXd y(3); y << 1.5, INFINITY, 1.0; @@ -135,6 +135,64 @@ TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_small_y_neg_inf) { sigma.transpose().eval()); } +// y == NEGATIVE_INFTY and a scaled difference that overflows to +// NEGATIVE_INFTY (y = -1e10, sigma = 1e-300) contribute log(1) = 0 to the +// value and zero to every partial. Absolute reference values from mpmath at +// z = 0.5, sigma = 2, because a comparison with the CPU passes when both are +// wrong in the same way. +TEST(ProbDistributionsLogisticLccdf, opencl_neg_inf_elements_reference) { + Eigen::VectorXd y(3); + y << -INFINITY, -1e10, 1.5; + Eigen::VectorXd mu(3); + mu << 0.5, 0.0, 0.5; + Eigen::VectorXd sigma(3); + sigma << 2.0, 1e-300, 2.0; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, + sigma); + + stan::math::var_value> y_cl + = stan::math::to_matrix_cl(y); + stan::math::var_value> mu_cl + = stan::math::to_matrix_cl(mu); + stan::math::var_value> sigma_cl + = stan::math::to_matrix_cl(sigma); + stan::math::var log_ccdf = stan::math::logistic_lccdf(y_cl, mu_cl, sigma_cl); + log_ccdf.grad(); + Eigen::VectorXd y_adj + = stan::math::from_matrix_cl(y_cl.adj()); + Eigen::VectorXd mu_adj + = stan::math::from_matrix_cl(mu_cl.adj()); + Eigen::VectorXd sigma_adj + = stan::math::from_matrix_cl(sigma_cl.adj()); + + EXPECT_NEAR(-0.97407698418010668, log_ccdf.val(), 1e-15); + for (int i = 0; i < 2; ++i) { + EXPECT_EQ(0.0, y_adj(i)); + EXPECT_EQ(0.0, mu_adj(i)); + EXPECT_EQ(0.0, sigma_adj(i)); + } + EXPECT_NEAR(-0.31122966560092728, y_adj(2), 1e-15); + EXPECT_NEAR(0.31122966560092728, mu_adj(2), 1e-15); + EXPECT_NEAR(0.15561483280046364, sigma_adj(2), 1e-15); + stan::math::recover_memory(); +} + +TEST(ProbDistributionsLogisticLccdf, opencl_neg_inf_and_pos_inf_reference) { + Eigen::VectorXd y(2); + y << -INFINITY, INFINITY; + Eigen::VectorXd mu(2); + mu << 0.0, 0.0; + Eigen::VectorXd sigma(2); + sigma << 1.0, 1.0; + + stan::math::matrix_cl y_cl(y); + stan::math::matrix_cl mu_cl(mu); + stan::math::matrix_cl sigma_cl(sigma); + EXPECT_EQ(stan::math::NEGATIVE_INFTY, + stan::math::logistic_lccdf(y_cl, mu_cl, sigma_cl)); +} + TEST(ProbDistributionsLogisticLccdf, opencl_broadcast_y) { int N = 3; diff --git a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp index 67d25a4f78c..f148ded9461 100644 --- a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp @@ -115,6 +115,46 @@ TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_y_pos_inf) { sigma); } +// A scaled difference that overflows to INFTY with finite arguments +// (y = 1e10, sigma = 1e-300) contributes log(1) = 0 and zero partials. +// Absolute reference values from mpmath at z = 0.5, sigma = 2, because a +// comparison with the CPU passes when both are wrong in the same way. +TEST(ProbDistributionsLogisticLcdf, opencl_pos_overflow_reference) { + Eigen::VectorXd y(2); + y << 1e10, 1.5; + Eigen::VectorXd mu(2); + mu << 0.0, 0.5; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 2.0; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, + sigma); + + stan::math::var_value> y_cl + = stan::math::to_matrix_cl(y); + stan::math::var_value> mu_cl + = stan::math::to_matrix_cl(mu); + stan::math::var_value> sigma_cl + = stan::math::to_matrix_cl(sigma); + stan::math::var log_cdf = stan::math::logistic_lcdf(y_cl, mu_cl, sigma_cl); + log_cdf.grad(); + Eigen::VectorXd y_adj + = stan::math::from_matrix_cl(y_cl.adj()); + Eigen::VectorXd mu_adj + = stan::math::from_matrix_cl(mu_cl.adj()); + Eigen::VectorXd sigma_adj + = stan::math::from_matrix_cl(sigma_cl.adj()); + + EXPECT_NEAR(-0.47407698418010668, log_cdf.val(), 1e-15); + EXPECT_EQ(0.0, y_adj(0)); + EXPECT_EQ(0.0, mu_adj(0)); + EXPECT_EQ(0.0, sigma_adj(0)); + EXPECT_NEAR(0.18877033439907272, y_adj(1), 1e-15); + EXPECT_NEAR(-0.18877033439907272, mu_adj(1), 1e-15); + EXPECT_NEAR(-0.094385167199536359, sigma_adj(1), 1e-15); + stan::math::recover_memory(); +} + TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_small_y_neg_inf) { int N = 3; int M = 2; diff --git a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp index 678eeb2381a..bbc4861a148 100644 --- a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp @@ -208,4 +208,32 @@ TEST(ProbDistributionsLogistic, opencl_sigma_mu_scalar) { mu, sigma); } +// y - mu overflows to NEGATIVE_INFTY with finite arguments; the value is +// NEGATIVE_INFTY, not inf - inf = NaN. The second row checks a finite value +// on both sides of the location against mpmath (z = -20.25 and z = 1.25, +// sigma = 2), because a comparison with the CPU passes when both are wrong in +// the same way. +TEST(ProbDistributionsLogistic, opencl_value_reference) { + Eigen::VectorXd y(2); + y << -1e308, 3.0; + Eigen::VectorXd mu(2); + mu << 1e308, 0.5; + Eigen::VectorXd sigma(2); + sigma << 1.0, 2.0; + stan::math::matrix_cl y_cl(y); + stan::math::matrix_cl mu_cl(mu); + stan::math::matrix_cl sigma_cl(sigma); + EXPECT_EQ(stan::math::NEGATIVE_INFTY, + stan::math::logistic_lpdf(y_cl, mu_cl, sigma_cl)); + + y << -40.0, 3.0; + mu << 0.5, 0.5; + sigma << 2.0, 2.0; + stan::math::matrix_cl y2_cl(y); + stan::math::matrix_cl mu2_cl(mu); + stan::math::matrix_cl sigma2_cl(sigma); + EXPECT_NEAR(-23.390152527021093, + stan::math::logistic_lpdf(y2_cl, mu2_cl, sigma2_cl), 1e-14); +} + #endif diff --git a/test/unit/math/rev/prob/logistic_cdf_test.cpp b/test/unit/math/rev/prob/logistic_cdf_test.cpp index ad1b424897f..b27ff8e1156 100644 --- a/test/unit/math/rev/prob/logistic_cdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_cdf_test.cpp @@ -74,3 +74,50 @@ TEST_F(AgradRev, logistic_cdf_location_scale) { EXPECT_DOUBLE_EQ(-deriv * cdf_value, mu.adj()); EXPECT_DOUBLE_EQ(-scaled_diff * deriv * cdf_value, sigma.adj()); } + +// A scaled difference that overflows to INFTY with finite arguments +// (y = 1e10, sigma = 1e-300) contributes a factor 1 and zero partials. The +// second element is unaffected. Reference values from mpmath at z = 0.5, +// sigma = 2. +TEST_F(AgradRev, logistic_cdf_pos_overflow_scale_partial) { + using stan::math::var; + Eigen::Matrix y(2); + y << 1e10, 1.5; + Eigen::Matrix mu(2); + mu << 0.0, 0.5; + Eigen::Matrix sigma(2); + sigma << 1e-300, 2.0; + + var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + EXPECT_NEAR(0.62245933120185456, cdf.val(), 1e-15); + EXPECT_EQ(0.0, y(0).adj()); + EXPECT_EQ(0.0, mu(0).adj()); + EXPECT_EQ(0.0, sigma(0).adj()); + EXPECT_NEAR(0.11750185610079724, y(1).adj(), 1e-15); + EXPECT_NEAR(-0.11750185610079724, mu(1).adj(), 1e-15); + EXPECT_NEAR(-0.058750928050398622, sigma(1).adj(), 1e-15); +} + +// A scaled difference that overflows to NEGATIVE_INFTY makes the cdf 0; every +// partial is then 0, also the scale partial of the overflowed element. +TEST_F(AgradRev, logistic_cdf_neg_overflow_scale_partial) { + using stan::math::var; + Eigen::Matrix y(2); + y << -1e10, 1.5; + Eigen::Matrix mu(2); + mu << 0.0, 0.5; + Eigen::Matrix sigma(2); + sigma << 1e-300, 2.0; + + var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + EXPECT_EQ(0.0, cdf.val()); + for (int i = 0; i < 2; ++i) { + EXPECT_EQ(0.0, y(i).adj()); + EXPECT_EQ(0.0, mu(i).adj()); + EXPECT_EQ(0.0, sigma(i).adj()); + } +} diff --git a/test/unit/math/rev/prob/logistic_lccdf_test.cpp b/test/unit/math/rev/prob/logistic_lccdf_test.cpp index b8048359a56..2e0f9f7102d 100644 --- a/test/unit/math/rev/prob/logistic_lccdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lccdf_test.cpp @@ -97,3 +97,66 @@ TEST_F(AgradRev, logistic_lccdf_pos_inf_zeroes_partials) { EXPECT_DOUBLE_EQ(0.0, mu.adj()); EXPECT_DOUBLE_EQ(0.0, sigma.adj()); } + +// y == NEGATIVE_INFTY and a scaled difference that overflows to +// NEGATIVE_INFTY (y = -1e10, sigma = 1e-300) contribute log(1) = 0 to the +// value and zero to every partial. The third element is unaffected. +// Reference values from mpmath at z = 0.5, sigma = 2. +TEST_F(AgradRev, logistic_lccdf_neg_inf_elements) { + using stan::math::var; + Eigen::Matrix y(3); + y << stan::math::NEGATIVE_INFTY, -1e10, 1.5; + Eigen::Matrix mu(3); + mu << 0.5, 0.0, 0.5; + Eigen::Matrix sigma(3); + sigma << 2.0, 1e-300, 2.0; + + var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + EXPECT_NEAR(-0.97407698418010668, log_ccdf.val(), 1e-15); + for (int i = 0; i < 2; ++i) { + EXPECT_EQ(0.0, y(i).adj()); + EXPECT_EQ(0.0, mu(i).adj()); + EXPECT_EQ(0.0, sigma(i).adj()); + } + EXPECT_NEAR(-0.31122966560092728, y(2).adj(), 1e-15); + EXPECT_NEAR(0.31122966560092728, mu(2).adj(), 1e-15); + EXPECT_NEAR(0.15561483280046364, sigma(2).adj(), 1e-15); +} + +// y == INFTY makes the result NEGATIVE_INFTY also when another element is +// NEGATIVE_INFTY. +TEST_F(AgradRev, logistic_lccdf_neg_inf_and_pos_inf) { + using stan::math::var; + Eigen::Matrix y(2); + y << stan::math::NEGATIVE_INFTY, stan::math::INFTY; + var mu = 0.0; + var sigma = 1.0; + + var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + EXPECT_EQ(stan::math::NEGATIVE_INFTY, log_ccdf.val()); + EXPECT_EQ(0.0, y(0).adj()); + EXPECT_EQ(0.0, y(1).adj()); + EXPECT_EQ(0.0, mu.adj()); + EXPECT_EQ(0.0, sigma.adj()); +} + +// The scaled difference overflows to INFTY with finite arguments. The exact +// scale partial z * inv_logit(z) / sigma overflows too. +TEST_F(AgradRev, logistic_lccdf_pos_overflow_scale_partial) { + using stan::math::var; + var y = 1e10; + var mu = 0.0; + var sigma = 1e-300; + + var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + EXPECT_EQ(stan::math::NEGATIVE_INFTY, log_ccdf.val()); + EXPECT_NEAR(-1e300, y.adj(), 1e285); + EXPECT_NEAR(1e300, mu.adj(), 1e285); + EXPECT_EQ(stan::math::INFTY, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lcdf_test.cpp b/test/unit/math/rev/prob/logistic_lcdf_test.cpp index f0936bcc2d5..bff94a0cf55 100644 --- a/test/unit/math/rev/prob/logistic_lcdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lcdf_test.cpp @@ -80,3 +80,45 @@ TEST_F(AgradRev, logistic_lcdf_pos_inf_element) { EXPECT_DOUBLE_EQ(-deriv, mu.adj()); EXPECT_DOUBLE_EQ(-1.5 * deriv, sigma.adj()); } + +// A scaled difference that overflows to INFTY with finite arguments +// (y = 1e10, sigma = 1e-300) contributes log(1) = 0 and zero partials, as +// y == INFTY does. The second element is unaffected. Reference values from +// mpmath at z = 0.5, sigma = 2. +TEST_F(AgradRev, logistic_lcdf_pos_overflow_scale_partial) { + using stan::math::var; + Eigen::Matrix y(2); + y << 1e10, 1.5; + Eigen::Matrix mu(2); + mu << 0.0, 0.5; + Eigen::Matrix sigma(2); + sigma << 1e-300, 2.0; + + var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + EXPECT_NEAR(-0.47407698418010668, log_cdf.val(), 1e-15); + EXPECT_EQ(0.0, y(0).adj()); + EXPECT_EQ(0.0, mu(0).adj()); + EXPECT_EQ(0.0, sigma(0).adj()); + EXPECT_NEAR(0.18877033439907272, y(1).adj(), 1e-15); + EXPECT_NEAR(-0.18877033439907272, mu(1).adj(), 1e-15); + EXPECT_NEAR(-0.094385167199536359, sigma(1).adj(), 1e-15); +} + +// The scaled difference overflows to NEGATIVE_INFTY with finite arguments. +// The exact scale partial -z * inv_logit(-z) / sigma overflows too. +TEST_F(AgradRev, logistic_lcdf_neg_overflow_scale_partial) { + using stan::math::var; + var y = -1e10; + var mu = 0.0; + var sigma = 1e-300; + + var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + EXPECT_EQ(stan::math::NEGATIVE_INFTY, log_cdf.val()); + EXPECT_NEAR(1e300, y.adj(), 1e285); + EXPECT_NEAR(-1e300, mu.adj(), 1e285); + EXPECT_EQ(stan::math::INFTY, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lpdf_test.cpp b/test/unit/math/rev/prob/logistic_lpdf_test.cpp index 28882ac4d87..d162f0e139c 100644 --- a/test/unit/math/rev/prob/logistic_lpdf_test.cpp +++ b/test/unit/math/rev/prob/logistic_lpdf_test.cpp @@ -59,3 +59,28 @@ TEST_F(AgradRev, logistic_lpdf_gradients_vectorized) { (0.5 * std::tanh(0.25) - 1.0) / 2.0 + (1.0 * std::tanh(0.5) - 1.0) / 2.0, sigma.adj()); } + +// y - mu overflows to NEGATIVE_INFTY with finite arguments. The value is +// NEGATIVE_INFTY, not inf - inf = NaN. +TEST_F(AgradRev, logistic_lpdf_value_neg_overflow) { + using stan::math::var; + var y = -1e308; + var mu = 1e308; + var sigma = 1.0; + + var logp = stan::math::logistic_lpdf(y, mu, sigma); + logp.grad(); + + EXPECT_EQ(stan::math::NEGATIVE_INFTY, logp.val()); + EXPECT_EQ(1.0, y.adj()); + EXPECT_EQ(-1.0, mu.adj()); +} + +// Values on both sides of the location. Reference value from mpmath at +// z = -20.25 and z = 1.25, sigma = 2. +TEST_F(AgradRev, logistic_lpdf_value_both_signs) { + Eigen::VectorXd y(2); + y << -40.0, 3.0; + EXPECT_NEAR(-23.390152527021093, stan::math::logistic_lpdf(y, 0.5, 2.0), + 1e-14); +} From a61a2779f79e01fbbc0135857e91f1a63f72771e Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 12 Oct 2026 23:31:58 +0300 Subject: [PATCH 10/10] Test the second derivative of logistic_lpdf in nested forward mode --- .../unit/math/mix/prob/logistic_lpdf_test.cpp | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 test/unit/math/mix/prob/logistic_lpdf_test.cpp diff --git a/test/unit/math/mix/prob/logistic_lpdf_test.cpp b/test/unit/math/mix/prob/logistic_lpdf_test.cpp new file mode 100644 index 00000000000..769ae04d131 --- /dev/null +++ b/test/unit/math/mix/prob/logistic_lpdf_test.cpp @@ -0,0 +1,54 @@ +#include +#include +#include +#include +#include +#include + +// First and second derivatives of logistic_lpdf with respect to y in nested +// forward mode and in fvar, against mpmath: +// d/dy = -tanh(z / 2) / sigma and d2/dy2 = -sech(z / 2)^2 / (2 sigma^2) +// with z = (y - mu) / sigma. The second derivative is the derivative of +// tanh(z / 2), which must not round to 0 for large |z|. +TEST_F(AgradRev, logistic_lpdf_second_derivative_y) { + using stan::math::fvar; + using stan::math::var; + // y, mu, sigma, d/dy, d2/dy2 + std::vector> cases + = {{1.0, 0.0, 1.0, -0.46211715726000976, -0.39322386648296371}, + {-20.0, 0.0, 1.0, 0.99999999587769276, -4.1223072278836987e-9}, + {40.0, 0.0, 1.0, -0.99999999999999999, -8.4967085105831779e-18}, + {100.0, 0.0, 1.0, -1.0, -7.4401519520416719e-44}, + {80.5, 0.5, 2.0, -0.5, -2.1241771276457945e-18}}; + for (const auto& c : cases) { + const double y = c[0]; + const double mu = c[1]; + const double sigma = c[2]; + + fvar> y_ffd(fvar(y, 1.0), fvar(1.0, 0.0)); + fvar> lp_ffd = stan::math::logistic_lpdf(y_ffd, mu, sigma); + EXPECT_NEAR(c[3], lp_ffd.d_.val_, 1e-14 * std::fabs(c[3])) << "y = " << y; + EXPECT_NEAR(c[4], lp_ffd.d_.d_, 1e-14 * std::fabs(c[4])) << "y = " << y; + + var y_v = y; + fvar y_fv(y_v, 1.0); + fvar lp_fv = stan::math::logistic_lpdf(y_fv, mu, sigma); + lp_fv.d_.grad(); + EXPECT_NEAR(c[3], lp_fv.d_.val(), 1e-14 * std::fabs(c[3])) << "y = " << y; + EXPECT_NEAR(c[4], y_v.adj(), 1e-14 * std::fabs(c[4])) << "y = " << y; + stan::math::recover_memory(); + } +} + +// The same through the vectorized tanh of an Eigen argument. With the +// direction (1, 1) the second directional derivative is the sum of the two +// second derivatives. +TEST_F(AgradRev, logistic_lpdf_second_derivative_y_vector) { + using stan::math::fvar; + Eigen::Matrix>, Eigen::Dynamic, 1> y(2); + y << fvar>(fvar(40.0, 1.0), fvar(1.0, 0.0)), + fvar>(fvar(100.0, 1.0), fvar(1.0, 0.0)); + fvar> lp = stan::math::logistic_lpdf(y, 0.0, 1.0); + const double d2 = -8.4967085105831779e-18 + -7.4401519520416719e-44; + EXPECT_NEAR(d2, lp.d_.d_, 1e-14 * std::fabs(d2)); +}