diff --git a/stan/math/fwd/fun/inv_logit.hpp b/stan/math/fwd/fun/inv_logit.hpp index 325a74a133a..2af1eb48962 100644 --- a/stan/math/fwd/fun/inv_logit.hpp +++ b/stan/math/fwd/fun/inv_logit.hpp @@ -18,8 +18,9 @@ namespace math { */ template * = nullptr> inline auto inv_logit(T&& x) { - return std::decay_t(inv_logit(x.val_), - x.d_ * inv_logit(x.val_) * (1 - inv_logit(x.val_))); + const auto val = inv_logit(x.val_); + // inv_logit(-x) rather than 1 - inv_logit(x), which cancels for x >> 0 + return std::decay_t(val, x.d_ * val * inv_logit(-x.val_)); } } // namespace math diff --git a/stan/math/fwd/fun/log_inv_logit_diff.hpp b/stan/math/fwd/fun/log_inv_logit_diff.hpp index e31f166d040..e7e962c5a00 100644 --- a/stan/math/fwd/fun/log_inv_logit_diff.hpp +++ b/stan/math/fwd/fun/log_inv_logit_diff.hpp @@ -21,11 +21,11 @@ namespace math { \f] \f[ - \frac{\partial }{\partial x} = -\frac{e^x}{e^y-e^x}-\frac{e^x}{e^x+1} + \frac{\partial}{\partial x} = \mbox{logit}^{-1}(-x) + \frac{1}{e^{x-y}-1} \f] \f[ - \frac{\partial }{\partial x} = -\frac{e^y}{e^x-e^y}-\frac{e^y}{e^y+1} + \frac{\partial}{\partial y} = -\mbox{logit}^{-1}(y) - \frac{1}{e^{x-y}-1} \f] * * @tparam T inner type of the fvar @@ -36,16 +36,17 @@ namespace math { */ template inline fvar log_inv_logit_diff(const fvar& x, const fvar& y) { + // both partials as sums of positive terms, no cancellation for x >> 0 + const auto c = inv(expm1(x.val_ - y.val_)); return fvar( log_inv_logit_diff(x.val_, y.val_), - -x.d_ * (inv(expm1(y.val_ - x.val_)) + inv_logit(x.val_)) - - y.d_ * (inv(expm1(x.val_ - y.val_)) + inv_logit(y.val_))); + x.d_ * (c + inv_logit(-x.val_)) - y.d_ * (c + inv_logit(y.val_))); } template inline fvar log_inv_logit_diff(const fvar& x, double y) { return fvar(log_inv_logit_diff(x.val_, y), - -x.d_ * (inv(expm1(y - x.val_)) + inv_logit(x.val_))); + x.d_ * (inv(expm1(x.val_ - y)) + inv_logit(-x.val_))); } template diff --git a/stan/math/opencl/rev/inv_logit.hpp b/stan/math/opencl/rev/inv_logit.hpp index aef65ee3c4a..7602590bef6 100644 --- a/stan/math/opencl/rev/inv_logit.hpp +++ b/stan/math/opencl/rev/inv_logit.hpp @@ -21,7 +21,7 @@ inline var_value> inv_logit(const var_value& A) { return make_callback_var( inv_logit(A.val()), [A](vari_value>& res) mutable { A.adj() += elt_multiply(res.adj(), - elt_multiply(res.val(), 1.0 - res.val())); + elt_multiply(res.val(), inv_logit(-A.val()))); }); } diff --git a/stan/math/opencl/rev/log_inv_logit_diff.hpp b/stan/math/opencl/rev/log_inv_logit_diff.hpp index e8eb5381f7a..82e87a534b0 100644 --- a/stan/math/opencl/rev/log_inv_logit_diff.hpp +++ b/stan/math/opencl/rev/log_inv_logit_diff.hpp @@ -38,13 +38,11 @@ inline var_value> log_inv_logit_diff(T_x&& x, T_y&& y) { return make_callback_var( log_inv_logit_diff(value_of(x_arena), value_of(y_arena)), [x_arena, y_arena](const vari_value>& res) mutable { + // both partials as sums of positive terms, no cancellation for x >> 0 + auto c = inv(expm1(value_of(x_arena) - value_of(y_arena))); adjoint_results(x_arena, y_arena) += expressions( - -elt_multiply(res.adj(), - inv(expm1(value_of(y_arena) - value_of(x_arena))) - + inv_logit(value_of(x_arena))), - -elt_multiply(res.adj(), - inv(expm1(value_of(x_arena) - value_of(y_arena))) - + inv_logit(value_of(y_arena)))); + elt_multiply(res.adj(), c + inv_logit(-value_of(x_arena))), + -elt_multiply(res.adj(), c + inv_logit(value_of(y_arena)))); }); } diff --git a/stan/math/prim/fun/log_inv_logit_diff.hpp b/stan/math/prim/fun/log_inv_logit_diff.hpp index 96bce37601b..ff12a5bb34f 100644 --- a/stan/math/prim/fun/log_inv_logit_diff.hpp +++ b/stan/math/prim/fun/log_inv_logit_diff.hpp @@ -20,11 +20,11 @@ namespace math { \f] \f[ - \frac{\partial }{\partial x} = -\frac{e^x}{e^y-e^x}-\frac{e^x}{e^x+1} + \frac{\partial}{\partial x} = \mbox{logit}^{-1}(-x) + \frac{1}{e^{x-y}-1} \f] \f[ - \frac{\partial }{\partial x} = -\frac{e^y}{e^x-e^y}-\frac{e^y}{e^y+1} + \frac{\partial}{\partial y} = -\mbox{logit}^{-1}(y) - \frac{1}{e^{x-y}-1} \f] * * @tparam T1 type of x argument diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp index 13b0e051b65..ed2f744d427 100644 --- a/stan/math/rev/fun/inv_logit.hpp +++ b/stan/math/rev/fun/inv_logit.hpp @@ -3,7 +3,12 @@ #include #include +#include +#include +#include #include +#include +#include namespace stan { namespace math { @@ -15,20 +20,27 @@ namespace math { * * The derivative of inverse logit is * - * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \mbox{logit}^{-1}(x) (1 - - * \mbox{logit}^{-1}(x))\f$. + * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \frac{e}{(1 + e)^2}\f$ + * with \f$e = \exp(-|x|)\f$. * - * @tparam T Arithmetic or a type inheriting from `EigenBase`. + * @tparam T a `var` or a `var_value` of a matrix type * @param a Argument variable. - * @return Inverse logit of argument. + * @return Inverse logit of argument, elementwise for a matrix. */ -template * = nullptr> +template < + typename T, require_var_t* = nullptr, + require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> inline auto inv_logit(T&& a) { - return make_callback_var(inv_logit(a.val()), [a](auto& vi) mutable { - as_array_or_scalar(a).adj() += as_array_or_scalar(vi.adj()) - * as_array_or_scalar(vi.val()) - * (1.0 - as_array_or_scalar(vi.val())); - }); + auto e = to_arena(exp(-abs(a.val()))); + return make_callback_var( + apply_scalar_binary( + [](double x, double e) { return (x >= 0 ? 1.0 : e) / (1.0 + e); }, + a.val(), e), + [a, e](auto& vi) mutable { + const auto& e_a = as_array_or_scalar(e); + as_array_or_scalar(a.adj()) + += as_array_or_scalar(vi.adj()) * e_a / ((1.0 + e_a) * (1.0 + e_a)); + }); } /** @@ -38,20 +50,28 @@ inline auto inv_logit(T&& a) { * * The derivative of inverse logit is * - * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \mbox{logit}^{-1}(x) (1 - - * \mbox{logit}^{-1}(x))\f$. + * \f$\frac{d}{dx} \mbox{logit}^{-1}(x) = \frac{e}{(1 + e)^2}\f$ + * with \f$e = \exp(-|x|)\f$. It is computed in the same pass as the value + * and stored for the reverse pass. * * @tparam T type of Eigen expression * @param x Eigen expression - * @return Inverse logit of argument. + * @return Inverse logit of each value in x. */ template * = nullptr> inline auto inv_logit(T&& x) { auto x_arena = to_arena(std::forward(x)); - arena_t ret = inv_logit(x_arena.val()); - reverse_pass_callback([x_arena, ret]() mutable { - x_arena.adj().array() - += ret.adj().array() * ret.val().array() * (1.0 - ret.val().array()); + arena_t ret(x_arena.rows(), x_arena.cols()); + arena_t> d(x_arena.rows(), + x_arena.cols()); + for (Eigen::Index i = 0; i < x_arena.size(); ++i) { + const double xi = x_arena.coeff(i).val(); + const double e = std::exp(-std::fabs(xi)); + ret.coeffRef(i) = (xi >= 0 ? 1.0 : e) / (1.0 + e); + d.coeffRef(i) = e / ((1.0 + e) * (1.0 + e)); + } + reverse_pass_callback([x_arena, ret, d]() mutable { + x_arena.adj().array() += ret.adj().array() * d.array(); }); return ret; } diff --git a/stan/math/rev/fun/log_inv_logit_diff.hpp b/stan/math/rev/fun/log_inv_logit_diff.hpp index 132c5d2c85b..b502c3a0614 100644 --- a/stan/math/rev/fun/log_inv_logit_diff.hpp +++ b/stan/math/rev/fun/log_inv_logit_diff.hpp @@ -21,11 +21,11 @@ namespace math { \f] \f[ - \frac{\partial }{\partial x} = -\frac{e^x}{e^y-e^x}-\frac{e^x}{e^x+1} + \frac{\partial}{\partial x} = \mbox{logit}^{-1}(-x) + \frac{1}{e^{x-y}-1} \f] \f[ - \frac{\partial }{\partial x} = -\frac{e^y}{e^x-e^y}-\frac{e^y}{e^y+1} + \frac{\partial}{\partial y} = -\mbox{logit}^{-1}(y) - \frac{1}{e^{x-y}-1} \f] * * @tparam T1 Type of x argument @@ -41,11 +41,10 @@ class log_inv_logit_diff_vv_vari : public op_vv_vari { log_inv_logit_diff_vv_vari(vari* avi, vari* bvi) : op_vv_vari(log_inv_logit_diff(avi->val_, bvi->val_), avi, bvi) {} void chain() { - avi_->adj_ - -= adj_ * (inv(expm1(bvi_->val_ - avi_->val_)) + inv_logit(avi_->val_)); - - bvi_->adj_ - -= adj_ * (inv(expm1(avi_->val_ - bvi_->val_)) + inv_logit(bvi_->val_)); + // both partials as sums of positive terms, no cancellation for x >> 0 + double c = inv(expm1(avi_->val_ - bvi_->val_)); + avi_->adj_ += adj_ * (c + inv_logit(-avi_->val_)); + bvi_->adj_ -= adj_ * (c + inv_logit(bvi_->val_)); } }; @@ -54,7 +53,8 @@ class log_inv_logit_diff_vd_vari : public op_vd_vari { log_inv_logit_diff_vd_vari(vari* avi, double b) : op_vd_vari(log_inv_logit_diff(avi->val_, b), avi, b) {} void chain() { - avi_->adj_ -= adj_ * (inv(expm1(bd_ - avi_->val_)) + inv_logit(avi_->val_)); + avi_->adj_ + += adj_ * (inv(expm1(avi_->val_ - bd_)) + inv_logit(-avi_->val_)); } }; diff --git a/test/unit/math/mix/constraint/lub_constrain_1_test.cpp b/test/unit/math/mix/constraint/lub_constrain_1_test.cpp index d76175a441c..09ffa4b2e9c 100644 --- a/test/unit/math/mix/constraint/lub_constrain_1_test.cpp +++ b/test/unit/math/mix/constraint/lub_constrain_1_test.cpp @@ -1,5 +1,6 @@ #include #include +#include TEST(mathMixMatFun, lub_constrain_scalars) { double x1 = 0.7; @@ -24,3 +25,12 @@ TEST(mathMixMatFun, lub_constrain_scalars) { lub_constrain_tests::expect(x1, lb_inf, ub_inf); lub_constrain_tests::expect(x2, lb_inf, ub_inf); } + +TEST(mathMixMatFun, lub_constrain_fwd_derivative_tails) { + using inv_logit_tail_refs::ref; + // derivatives of -2 + 5 inv_logit(x) + inv_logit_tail_refs::expect_fwd_derivatives( + [](const auto& x) { return stan::math::lub_constrain(x, -2.0, 3.0); }, + [](const ref& r) { return 5 * r.d1; }, + [](const ref& r) { return 5 * r.d2; }); +} diff --git a/test/unit/math/mix/fun/inv_logit_tail_refs.hpp b/test/unit/math/mix/fun/inv_logit_tail_refs.hpp new file mode 100644 index 00000000000..25163d6b495 --- /dev/null +++ b/test/unit/math/mix/fun/inv_logit_tail_refs.hpp @@ -0,0 +1,86 @@ +#ifndef TEST_UNIT_MATH_MIX_FUN_INV_LOGIT_TAIL_REFS_HPP +#define TEST_UNIT_MATH_MIX_FUN_INV_LOGIT_TAIL_REFS_HPP + +#include +#include +#include +#include +#include + +namespace inv_logit_tail_refs { + +// inv_logit(x), inv_logit(-x), inv_logit'(x) and inv_logit''(x), computed +// with mpmath at 60 digits and rounded to double (0 where they underflow). +struct ref { + double x; + double s; + double s_neg; + double d1; + double d2; +}; + +inline std::vector refs() { + return {{-800, 0, 1.0, 0, 0}, + {-40, 4.248354255291589e-18, 1.0, 4.248354255291589e-18, + 4.248354255291589e-18}, + {-37, 8.533047625744065e-17, 0.9999999999999999, + 8.533047625744065e-17, 8.533047625744063e-17}, + {-30, 9.357622968839299e-14, 0.9999999999999064, + 9.357622968838423e-14, 9.357622968836672e-14}, + {-20, 2.0611536181902037e-09, 0.9999999979388464, + 2.061153613941849e-09, 2.061153605445141e-09}, + {-1, 0.2689414213699951, 0.7310585786300049, 0.19661193324148185, + 0.09085774767294841}, + {0, 0.5, 0.5, 0.25, 0}, + {1, 0.7310585786300049, 0.2689414213699951, 0.19661193324148185, + -0.09085774767294841}, + {20, 0.9999999979388464, 2.0611536181902037e-09, + 2.061153613941849e-09, -2.061153605445141e-09}, + {30, 0.9999999999999064, 9.357622968839299e-14, 9.357622968838423e-14, + -9.357622968836672e-14}, + {37, 0.9999999999999999, 8.533047625744065e-17, 8.533047625744065e-17, + -8.533047625744063e-17}, + {40, 1.0, 4.248354255291589e-18, 4.248354255291589e-18, + -4.248354255291589e-18}, + {100, 1.0, 3.720075976020836e-44, 3.720075976020836e-44, + -3.720075976020836e-44}, + {700, 1.0, 9.85967654375977e-305, 9.85967654375977e-305, + -9.85967654375977e-305}, + {800, 1.0, 0, 0, 0}}; +} + +// Relative check, absolute below the smallest normal double. +inline void expect_rel(double expected, double actual, double x, + double tol = 1e-13) { + EXPECT_NEAR( + expected, actual, + tol * std::fmax(std::fabs(expected), std::numeric_limits::min())) + << "x = " << x; +} + +// First and second derivatives of f at every reference point through +// fvar, fvar> and fvar, against d1(r) and d2(r). +template +void expect_fwd_derivatives(const F& f, const D1& d1, const D2& d2) { + using stan::math::fvar; + using stan::math::var; + for (const auto& r : refs()) { + expect_rel(d1(r), f(fvar(r.x, 1)).d_, r.x); + + fvar> ff(fvar(r.x, 1), fvar(1, 0)); + fvar> y_ff = f(ff); + expect_rel(d1(r), y_ff.d_.val_, r.x); + expect_rel(d2(r), y_ff.d_.d_, r.x); + + fvar fv(r.x, 1); + fvar y_fv = f(fv); + expect_rel(d1(r), y_fv.d_.val(), r.x); + y_fv.d_.grad(); + expect_rel(d2(r), fv.val_.adj(), r.x); + stan::math::recover_memory(); + } +} + +} // namespace inv_logit_tail_refs + +#endif diff --git a/test/unit/math/mix/fun/inv_logit_test.cpp b/test/unit/math/mix/fun/inv_logit_test.cpp index 8fee66f7eeb..8a6a058461c 100644 --- a/test/unit/math/mix/fun/inv_logit_test.cpp +++ b/test/unit/math/mix/fun/inv_logit_test.cpp @@ -1,4 +1,5 @@ #include +#include TEST(mathMixMatFun, invLogit) { auto f = [](const auto& x1) { return stan::math::inv_logit(x1); }; @@ -12,3 +13,87 @@ TEST(mathMixMatFun, invLogit) { stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args)); stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args)); } + +TEST(mathMixMatFun, invLogitDerivativesLargeArgs) { + using inv_logit_tail_refs::ref; + for (const auto& r : inv_logit_tail_refs::refs()) { + stan::math::var a = r.x; + stan::math::inv_logit(a).grad(); + inv_logit_tail_refs::expect_rel(r.d1, a.adj(), r.x); + stan::math::recover_memory(); + } + inv_logit_tail_refs::expect_fwd_derivatives( + [](const auto& x) { return stan::math::inv_logit(x); }, + [](const ref& r) { return r.d1; }, [](const ref& r) { return r.d2; }); +} + +TEST(mathMixMatFun, invLogitDerivativesLargeArgsContainers) { + using inv_logit_tail_refs::expect_rel; + using stan::math::fvar; + using stan::math::inv_logit; + using stan::math::sum; + using stan::math::var; + std::vector ps = inv_logit_tail_refs::refs(); + int n = ps.size(); + Eigen::VectorXd x(n); + Eigen::VectorXd d1(n); + Eigen::VectorXd d2(n); + for (int i = 0; i < n; ++i) { + x(i) = ps[i].x; + d1(i) = ps[i].d1; + d2(i) = ps[i].d2; + } + + stan::math::var_value xv(x); + sum(inv_logit(xv)).grad(); + for (int i = 0; i < n; ++i) { + expect_rel(d1(i), xv.adj()(i), x(i)); + } + stan::math::recover_memory(); + + stan::math::var_value xr(x.transpose()); + sum(inv_logit(xr)).grad(); + for (int i = 0; i < n; ++i) { + expect_rel(d1(i), xr.adj()(i), x(i)); + } + stan::math::recover_memory(); + + // inv_logit' is even + Eigen::MatrixXd x2(n, 2); + x2 << x, -x; + stan::math::var_value xmv(x2); + sum(inv_logit(xmv)).grad(); + for (int i = 0; i < n; ++i) { + expect_rel(d1(i), xmv.adj()(i, 0), x(i)); + expect_rel(d1(i), xmv.adj()(i, 1), -x(i)); + } + stan::math::recover_memory(); + + Eigen::Matrix xm = x; + sum(inv_logit(xm)).grad(); + for (int i = 0; i < n; ++i) { + expect_rel(d1(i), xm(i).adj(), x(i)); + } + stan::math::recover_memory(); + + std::vector xs(x.data(), x.data() + n); + sum(inv_logit(xs)).grad(); + for (int i = 0; i < n; ++i) { + expect_rel(d1(i), xs[i].adj(), x(i)); + } + stan::math::recover_memory(); + + Eigen::Matrix, Eigen::Dynamic, 1> xf(n); + for (int i = 0; i < n; ++i) { + xf(i) = fvar(x(i), 1); + } + Eigen::Matrix, Eigen::Dynamic, 1> yf = inv_logit(xf); + for (int i = 0; i < n; ++i) { + expect_rel(d1(i), yf(i).d_.val(), x(i)); + } + sum(yf).d_.grad(); + for (int i = 0; i < n; ++i) { + expect_rel(d2(i), xf(i).val_.adj(), x(i)); + } + stan::math::recover_memory(); +} diff --git a/test/unit/math/mix/fun/log1m_inv_logit_test.cpp b/test/unit/math/mix/fun/log1m_inv_logit_test.cpp index 175ead491ac..6e730bfb0ae 100644 --- a/test/unit/math/mix/fun/log1m_inv_logit_test.cpp +++ b/test/unit/math/mix/fun/log1m_inv_logit_test.cpp @@ -1,4 +1,5 @@ #include +#include TEST(mathMixMatFun, log1mInvLogit) { auto f = [](const auto& x1) { return stan::math::log1m_inv_logit(x1); }; @@ -25,3 +26,11 @@ TEST(mathMixMatFun, log1minvlogit_varmat) { } expect_ad_vector_matvar(f, A); } + +TEST(mathMixMatFun, log1mInvLogitFwdDerivativeTails) { + using inv_logit_tail_refs::ref; + // d/dx = -inv_logit(x), d^2/dx^2 = -inv_logit'(x) + inv_logit_tail_refs::expect_fwd_derivatives( + [](const auto& x) { return stan::math::log1m_inv_logit(x); }, + [](const ref& r) { return -r.s; }, [](const ref& r) { return -r.d1; }); +} diff --git a/test/unit/math/mix/fun/log1p_exp_test.cpp b/test/unit/math/mix/fun/log1p_exp_test.cpp index 9641c036c0e..13167f8a464 100644 --- a/test/unit/math/mix/fun/log1p_exp_test.cpp +++ b/test/unit/math/mix/fun/log1p_exp_test.cpp @@ -1,4 +1,7 @@ #include +#include +#include +#include TEST(mathMixMatFun, log1pExp) { auto f = [](const auto& x1) { return stan::math::log1p_exp(x1); }; @@ -13,3 +16,46 @@ TEST(mathMixMatFun, log1pExp) { stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args)); stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args)); } + +TEST(mathMixMatFun, log1pExpHigherDerivativesLargeArgs) { + using stan::math::fvar; + using stan::math::log1p_exp; + using stan::math::var; + // {x, log1p_exp''(x), log1p_exp'''(x)}: mpmath, 50 digits; 0 where the + // value underflows (e^-800) + std::vector> points{ + {-800, 0, 0}, + {-40, 4.2483542552915889592e-18, 4.2483542552915889231e-18}, + {-20, 2.0611536139418493437e-9, 2.0611536054451408856e-9}, + {-1, 1.9661193324148185254e-1, 9.0857747672948409442e-2}, + {0, 0.25, 0}, + {1, 1.9661193324148185254e-1, -9.0857747672948409442e-2}, + {20, 2.0611536139418493437e-9, -2.0611536054451408856e-9}, + {30, 9.3576229688384233028e-14, -9.3576229688366720006e-14}, + {37, 8.533047625744064338e-17, -8.5330476257440628818e-17}, + {40, 4.2483542552915889592e-18, -4.2483542552915889231e-18}, + {100, 3.720075976020835963e-44, -3.720075976020835963e-44}, + {700, 9.8596765437597708567e-305, -9.8596765437597708567e-305}, + {800, 0, 0}}; + for (const auto& p : points) { + double x = p[0]; + double d2 = p[1]; + double d3 = p[2]; + double tol2 = 1e-13 * std::fmax(d2, std::numeric_limits::min()); + double tol3 + = 1e-13 * std::fmax(std::fabs(d3), std::numeric_limits::min()); + + fvar> ff(fvar(x, 1), fvar(1, 0)); + EXPECT_NEAR(log1p_exp(ff).d_.d_, d2, tol2) << "x = " << x; + + fvar fv(x, 1); + log1p_exp(fv).d_.grad(); + EXPECT_NEAR(fv.val_.adj(), d2, tol2) << "x = " << x; + stan::math::recover_memory(); + + fvar> ffv(fvar(x, 1), fvar(1, 0)); + log1p_exp(ffv).d_.d_.grad(); + EXPECT_NEAR(ffv.val_.val_.adj(), d3, tol3) << "x = " << x; + stan::math::recover_memory(); + } +} diff --git a/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp b/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp index c5c4f31f474..ca7bde6a05b 100644 --- a/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp +++ b/test/unit/math/mix/fun/log_inv_logit_diff_test.cpp @@ -1,5 +1,6 @@ #include #include +#include TEST(mathMixScalFun, logInvLogitDiff) { auto f = [](const auto& x1, const auto& x2) { @@ -28,3 +29,67 @@ TEST(mathMixScalFun, logInvLogitDiff_vec) { in2 << 0.5, -3.4; stan::test::expect_ad_vectorized_binary(f, in1, in2); } + +TEST(mathMixScalFun, logInvLogitDiffGradientLargeArgs) { + using stan::math::fvar; + using stan::math::log_inv_logit_diff; + using stan::math::var; + // {x, y, d/dx, d/dy}: mpmath, 50 digits + std::vector> points{ + {40, 0, 8.4967085105831779907e-18, -5.0000000000000000425e-1}, + {100, 1, 1.3832290902125321262e-43, -7.3105857863000487925e-1}, + {30, -30, 9.3576229688401746049e-14, -9.3576229688401746049e-14}, + {37, 36, 5.8197670686932650972e-1, -1.5819767068693261924}, + {-30, -40, 1.0000454019909161115, -4.5401991009692016683e-5}, + {1, 0, 8.5091812823932154513e-1, -1.0819767068693264244}, + {700, 20, 4.7835719068902113962e-296, -9.9999999793884638181e-1}, + {0.5, -1.0, 6.647575855870136797e-1, -5.5615833815886336509e-1}}; + for (const auto& p : points) { + double x = p[0]; + double y = p[1]; + double dx = p[2]; + double dy = p[3]; + + var xv = x; + var yv = y; + log_inv_logit_diff(xv, yv).grad(); + EXPECT_NEAR(xv.adj(), dx, 1e-13 * dx) << x << ", " << y; + EXPECT_NEAR(yv.adj(), dy, -1e-13 * dy) << x << ", " << y; + stan::math::recover_memory(); + + var xv2 = x; + log_inv_logit_diff(xv2, y).grad(); + EXPECT_NEAR(xv2.adj(), dx, 1e-13 * dx) << x << ", " << y; + stan::math::recover_memory(); + + EXPECT_NEAR(log_inv_logit_diff(fvar(x, 1), y).d_, dx, 1e-13 * dx) + << x << ", " << y; + EXPECT_NEAR(log_inv_logit_diff(fvar(x, 1), fvar(y, 0)).d_, + dx, 1e-13 * dx) + << x << ", " << y; + EXPECT_NEAR(log_inv_logit_diff(fvar(x, 0), fvar(y, 1)).d_, + dy, -1e-13 * dy) + << x << ", " << y; + } +} + +TEST(mathMixScalFun, logInvLogitDiffGradientEqualArgs) { + using stan::math::fvar; + using stan::math::log_inv_logit_diff; + using stan::math::var; + // behaviour change: d/dx at x == y is +inf, the limit from x > y (was -inf) + double inf = std::numeric_limits::infinity(); + var xv = 2.0; + var yv = 2.0; + log_inv_logit_diff(xv, yv).grad(); + EXPECT_EQ(xv.adj(), inf); + EXPECT_EQ(yv.adj(), -inf); + stan::math::recover_memory(); + + var xv2 = 2.0; + log_inv_logit_diff(xv2, 2.0).grad(); + EXPECT_EQ(xv2.adj(), inf); + stan::math::recover_memory(); + + EXPECT_EQ(log_inv_logit_diff(fvar(2, 1), 2.0).d_, inf); +} diff --git a/test/unit/math/mix/fun/log_inv_logit_test.cpp b/test/unit/math/mix/fun/log_inv_logit_test.cpp index f94478e21c9..3bad8322319 100644 --- a/test/unit/math/mix/fun/log_inv_logit_test.cpp +++ b/test/unit/math/mix/fun/log_inv_logit_test.cpp @@ -1,4 +1,5 @@ #include +#include TEST(mathMixMatFun, logInvLogit) { auto f = [](const auto& x1) { return stan::math::log_inv_logit(x1); }; @@ -25,3 +26,11 @@ TEST(mathMixMatFun, logInvLogitVarMat) { } expect_ad_vector_matvar(f, A); } + +TEST(mathMixMatFun, logInvLogitFwdDerivativeTails) { + using inv_logit_tail_refs::ref; + // d/dx = inv_logit(-x), d^2/dx^2 = -inv_logit'(x) + inv_logit_tail_refs::expect_fwd_derivatives( + [](const auto& x) { return stan::math::log_inv_logit(x); }, + [](const ref& r) { return r.s_neg; }, [](const ref& r) { return -r.d1; }); +} diff --git a/test/unit/math/opencl/rev/inv_logit_test.cpp b/test/unit/math/opencl/rev/inv_logit_test.cpp index b8af60f0778..8549b411aef 100644 --- a/test/unit/math/opencl/rev/inv_logit_test.cpp +++ b/test/unit/math/opencl/rev/inv_logit_test.cpp @@ -27,4 +27,38 @@ TEST(OpenCLinv_logit, prim_rev_values_large) { stan::math::test::compare_cpu_opencl_prim_rev(inv_logit_functor, a); } +TEST(OpenCLinv_logit, rev_derivative_large_args) { + using stan::math::matrix_cl; + using stan::math::var_value; + Eigen::VectorXd x(8); + x << -40, -20, 0, 20, 30, 37, 40, 100; + // inv_logit'(x): mpmath, 50 digits + Eigen::VectorXd d(8); + d << 4.2483542552915889592e-18, 2.0611536139418493437e-9, 0.25, + 2.0611536139418493437e-9, 9.3576229688384233028e-14, + 8.533047625744064338e-17, 4.2483542552915889592e-18, + 3.720075976020835963e-44; + const var_value> x_cl = stan::math::to_matrix_cl(x); + stan::math::sum(stan::math::inv_logit(x_cl)).grad(); + Eigen::VectorXd adj = stan::math::from_matrix_cl(x_cl.adj()); + for (int i = 0; i < x.size(); ++i) { + EXPECT_NEAR(adj(i), d(i), 1e-13 * d(i)) << "x = " << x(i); + } + stan::math::recover_memory(); +} + +TEST(OpenCLinv_logit, rev_non_const_lvalue) { + using stan::math::matrix_cl; + using stan::math::var_value; + Eigen::VectorXd x(2); + x << 0, 40; + var_value> x_cl = stan::math::to_matrix_cl(x); + var_value> y_cl = stan::math::inv_logit(x_cl); + stan::math::sum(y_cl).grad(); + Eigen::VectorXd adj = stan::math::from_matrix_cl(x_cl.adj()); + EXPECT_NEAR(adj(0), 0.25, 1e-15); + EXPECT_NEAR(adj(1), 4.2483542552915889592e-18, 1e-30); + stan::math::recover_memory(); +} + #endif diff --git a/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp b/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp index b8df60796da..3827d77ad0f 100644 --- a/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp +++ b/test/unit/math/opencl/rev/log_inv_logit_diff_test.cpp @@ -57,4 +57,34 @@ TEST(OpenCLMatrix_log_inv_logit_diff, prim_rev_scalar_values_large) { a); } +TEST(OpenCLMatrix_log_inv_logit_diff, rev_gradient_large_args) { + using stan::math::matrix_cl; + using stan::math::var_value; + Eigen::VectorXd x(8); + Eigen::VectorXd y(8); + x << 40, 100, 30, 37, -30, 1, 700, 0.5; + y << 0, 1, -30, 36, -40, 0, 20, -1; + // d/dx and d/dy: mpmath, 50 digits + Eigen::VectorXd dx(8); + Eigen::VectorXd dy(8); + dx << 8.4967085105831779907e-18, 1.3832290902125321262e-43, + 9.3576229688401746049e-14, 5.8197670686932650972e-1, + 1.0000454019909161115, 8.5091812823932154513e-1, + 4.7835719068902113962e-296, 6.647575855870136797e-1; + dy << -5.0000000000000000425e-1, -7.3105857863000487925e-1, + -9.3576229688401746049e-14, -1.5819767068693261924, + -4.5401991009692016683e-5, -1.0819767068693264244, + -9.9999999793884638181e-1, -5.5615833815886336509e-1; + var_value> x_cl = stan::math::to_matrix_cl(x); + var_value> y_cl = stan::math::to_matrix_cl(y); + stan::math::sum(stan::math::log_inv_logit_diff(x_cl, y_cl)).grad(); + Eigen::VectorXd x_adj = stan::math::from_matrix_cl(x_cl.adj()); + Eigen::VectorXd y_adj = stan::math::from_matrix_cl(y_cl.adj()); + for (int i = 0; i < x.size(); ++i) { + EXPECT_NEAR(x_adj(i), dx(i), 1e-13 * dx(i)) << x(i) << ", " << y(i); + EXPECT_NEAR(y_adj(i), dy(i), -1e-13 * dy(i)) << x(i) << ", " << y(i); + } + stan::math::recover_memory(); +} + #endif