Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions stan/math/opencl/kernel_generator/elt_function_cl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include <stan/math/opencl/kernels/device_functions/log1m_exp.hpp>
#include <stan/math/opencl/kernels/device_functions/log1m_inv_logit.hpp>
#include <stan/math/opencl/kernels/device_functions/log1p_exp.hpp>
#include <stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp>
#include <stan/math/opencl/kernels/device_functions/logit.hpp>
#include <stan/math/opencl/kernels/device_functions/multiply_log.hpp>
#include <stan/math/opencl/kernels/device_functions/Phi.hpp>
Expand Down Expand Up @@ -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)
Expand Down
43 changes: 43 additions & 0 deletions stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp
Original file line number Diff line number Diff line change
@@ -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 <stan/math/opencl/stringify.hpp>
#include <string>

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
16 changes: 8 additions & 8 deletions stan/math/opencl/prim/logistic_cdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,16 +66,16 @@ inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> logistic_cdf(
auto any_y_neg_inf = colwise_max(cast<char>(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 = 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));

matrix_cl<char> any_y_neg_inf_cl;
matrix_cl<double> P_cl;
Expand Down
14 changes: 5 additions & 9 deletions stan/math/opencl/prim/logistic_lccdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,16 +66,12 @@ inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> logistic_lccdf(
auto any_y_neg_inf = colwise_max(cast<char>(y_val == NEGATIVE_INFTY));
auto any_y_pos_inf = colwise_max(cast<char>(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 = logistic_tail_deriv(-scaled_diff, sigma_val);
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<char> any_y_neg_inf_cl;
matrix_cl<char> any_y_pos_inf_cl;
Expand Down
18 changes: 9 additions & 9 deletions stan/math/opencl/prim/logistic_lcdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,16 +66,16 @@ inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> logistic_lcdf(
auto any_y_neg_inf = colwise_max(cast<char>(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));

// 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 = logistic_tail_deriv(scaled_diff, sigma_val);
auto y_deriv = select(cond, 0.0, deriv);
auto mu_deriv = -y_deriv;
auto sigma_deriv = elt_multiply(y_deriv, mu_minus_y_div_sigma);
auto sigma_deriv = select(cond, 0.0, elt_multiply(-deriv, scaled_diff));

matrix_cl<char> any_y_neg_inf_cl;
matrix_cl<double> P_cl;
Expand Down
17 changes: 6 additions & 11 deletions stan/math/opencl/prim/logistic_lpdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -81,18 +81,13 @@ inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> logistic_lpdf(
= colwise_sum(static_select<include_summand<propto, T_scale_cl>::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);
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);
// 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<double> logp_cl;
matrix_cl<double> y_deriv_cl;
Expand Down
38 changes: 19 additions & 19 deletions stan/math/prim/prob/logistic_cdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,16 @@
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/exp.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/log1m_inv_logit.hpp>
#include <stan/math/prim/fun/max_size.hpp>
#include <stan/math/prim/fun/scalar_seq_view.hpp>
#include <stan/math/prim/fun/size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <stan/math/prim/fun/value_of.hpp>
#include <stan/math/prim/fun/value_of_rec.hpp>
#include <stan/math/prim/fun/inv_logit.hpp>
#include <stan/math/prim/prob/logistic_lpdf.hpp>
#include <stan/math/prim/functor/partials_propagator.hpp>
#include <cmath>

namespace stan {
namespace math {
Expand All @@ -26,7 +27,6 @@ inline return_type_t<T_y, T_loc, T_scale> logistic_cdf(const T_y& y,
const T_loc& mu,
const T_scale& sigma) {
using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
using std::exp;
using T_y_ref = ref_type_t<T_y>;
using T_mu_ref = ref_type_t<T_loc>;
using T_sigma_ref = ref_type_t<T_scale>;
Expand Down Expand Up @@ -69,26 +69,26 @@ inline return_type_t<T_y, T_loc, T_scale> 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;
Comment on lines +74 to 76

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can use the same log-calculations as the lcdf (see below), accumulating on the log scale and exponentiating the result (for the derivatives as well)

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Claude assisted answer:

Measured against a long double reference, exp(Σ log inv_logit) is 749 ulp at N=100, z=−0.5 against 3.0 for the product, and 641 vs 1.0 at N=10, z=−70 — exp amplifies the summed log by |log P|·eps. It also doesn't extend the range: both forms underflow at the same points, including for a single observation at z = −745.2. inv_logit has no cancellation to avoid, unlike the difference in exp_mod_normal_cdf where the log scale genuinely buys something.

The rescaling partials *= cdf does drop the case where the cdf underflows but cdf × ∂log F/∂θ is representable. That's the same σ < 1e-16 window as the tail fix, and there the cdf returns exactly 0, so logistic_lcdf is the right function. If you want this to be changes, it needs sign tracking to exponentiate each partial separately.


if constexpr (is_autodiff_v<T_y>) {
partials<0>(ops_partials)[n]
+= exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn;
}
if constexpr (is_autodiff_v<T_loc>) {
partials<1>(ops_partials)[n]
+= -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn;
}
if constexpr (is_autodiff_v<T_scale>) {
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<T_y, T_loc, T_scale>) {
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<T_y>) {
partials<0>(ops_partials)[n] += deriv;
}
if constexpr (is_autodiff_v<T_loc>) {
partials<1>(ops_partials)[n] -= deriv;
}
if constexpr (is_autodiff_v<T_scale>) {
partials<2>(ops_partials)[n] -= scaled_diff * deriv;
}
}
}

Expand Down
54 changes: 25 additions & 29 deletions stan/math/prim/prob/logistic_lccdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,17 @@
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/exp.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/inv_logit.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/log1m_inv_logit.hpp>
#include <stan/math/prim/fun/log_inv_logit.hpp>
#include <stan/math/prim/fun/max_size.hpp>
#include <stan/math/prim/fun/scalar_seq_view.hpp>
#include <stan/math/prim/fun/size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <stan/math/prim/fun/value_of.hpp>
#include <stan/math/prim/prob/logistic_lpdf.hpp>
#include <stan/math/prim/fun/value_of_rec.hpp>
#include <stan/math/prim/functor/partials_propagator.hpp>
#include <cmath>

namespace stan {
namespace math {
Expand All @@ -26,8 +27,6 @@ inline return_type_t<T_y, T_loc, T_scale> logistic_lccdf(const T_y& y,
const T_loc& mu,
const T_scale& sigma) {
using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
using std::exp;
using std::log;
using T_y_ref = ref_type_t<T_y>;
using T_mu_ref = ref_type_t<T_loc>;
using T_sigma_ref = ref_type_t<T_scale>;
Expand All @@ -54,42 +53,39 @@ inline return_type_t<T_y, T_loc, T_scale> 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_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<T_y>) {
partials<0>(ops_partials)[n]
-= exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn;
}
if constexpr (is_autodiff_v<T_loc>) {
partials<1>(ops_partials)[n]
-= -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn;
}
if constexpr (is_autodiff_v<T_scale>) {
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<T_y, T_loc, T_scale>) {
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<T_y>) {
partials<0>(ops_partials)[n] -= deriv;
}
if constexpr (is_autodiff_v<T_loc>) {
partials<1>(ops_partials)[n] += deriv;
}
if constexpr (is_autodiff_v<T_scale>) {
partials<2>(ops_partials)[n] += scaled_diff * deriv;
}
}
}
return ops_partials.build(P);
Expand Down
41 changes: 20 additions & 21 deletions stan/math/prim/prob/logistic_lcdf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,17 @@
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/exp.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/inv_logit.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/log1m_inv_logit.hpp>
#include <stan/math/prim/fun/log_inv_logit.hpp>
#include <stan/math/prim/fun/scalar_seq_view.hpp>
#include <stan/math/prim/fun/max_size.hpp>
#include <stan/math/prim/fun/size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <stan/math/prim/fun/value_of.hpp>
#include <stan/math/prim/prob/logistic_lpdf.hpp>
#include <stan/math/prim/fun/value_of_rec.hpp>
#include <stan/math/prim/functor/partials_propagator.hpp>
#include <cmath>

namespace stan {
namespace math {
Expand All @@ -26,8 +27,6 @@ inline return_type_t<T_y, T_loc, T_scale> logistic_lcdf(const T_y& y,
const T_loc& mu,
const T_scale& sigma) {
using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
using std::exp;
using std::log;
using T_y_ref = ref_type_t<T_y>;
using T_mu_ref = ref_type_t<T_loc>;
using T_sigma_ref = ref_type_t<T_scale>;
Expand Down Expand Up @@ -70,25 +69,25 @@ inline return_type_t<T_y, T_loc, T_scale> 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<T_y>) {
partials<0>(ops_partials)[n]
+= exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn;
}
if constexpr (is_autodiff_v<T_loc>) {
partials<1>(ops_partials)[n]
+= -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn;
}
if constexpr (is_autodiff_v<T_scale>) {
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<T_y, T_loc, T_scale>) {
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<T_y>) {
partials<0>(ops_partials)[n] += deriv;
}
if constexpr (is_autodiff_v<T_loc>) {
partials<1>(ops_partials)[n] -= deriv;
}
if constexpr (is_autodiff_v<T_scale>) {
partials<2>(ops_partials)[n] -= scaled_diff * deriv;
}
}
}
return ops_partials.build(P);
Expand Down
Loading