diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp index ffd56a9297e..16647f1d144 100644 --- a/stan/math/fwd/fun.hpp +++ b/stan/math/fwd/fun.hpp @@ -72,6 +72,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/fwd/fun/log_cosh.hpp b/stan/math/fwd/fun/log_cosh.hpp new file mode 100644 index 00000000000..580e6efd368 --- /dev/null +++ b/stan/math/fwd/fun/log_cosh.hpp @@ -0,0 +1,27 @@ +#ifndef STAN_MATH_FWD_FUN_LOG_COSH_HPP +#define STAN_MATH_FWD_FUN_LOG_COSH_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the natural logarithm of the hyperbolic cosine of the argument. + * + * @tparam T inner type of the forward-mode autodiff variable + * @param x argument + * @return `log(cosh(x))`, with derivative `tanh(x)` + */ +template +inline fvar log_cosh(const fvar& x) { + return fvar(log_cosh(x.val_), x.d_ * tanh(x.val_)); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/fwd/fun/tanh.hpp b/stan/math/fwd/fun/tanh.hpp index 9699163b226..7ab38dceb97 100644 --- a/stan/math/fwd/fun/tanh.hpp +++ b/stan/math/fwd/fun/tanh.hpp @@ -4,17 +4,41 @@ #include #include #include +#include +#include #include +#include #include #include namespace stan { namespace math { +/** + * Return the hyperbolic tangent of the argument. + * + * For |x| >= 0.5 the derivative 1 - tanh(x)^2 is computed as + * 4 e / (1 + e)^2 with e = exp(-2 |x|): the difference 1 - tanh(x)^2 loses + * all digits for large |x|, and 1 / cosh(x)^2 makes the next derivative NaN + * once cosh(x)^2 overflows. For |x| < 0.5 the difference is accurate. + * Neither branch uses fabs(x), whose derivative is 0 at x = 0 and would + * make the third derivative of tanh at 0 vanish in nested forward mode. + * + * @tparam T inner type of the forward-mode autodiff variable + * @param[in] x argument + * @return hyperbolic tangent of the argument + */ template inline fvar tanh(const fvar& x) { - T u = tanh(x.val_); - return fvar(u, x.d_ * (1 - u * u)); + const T u = tanh(x.val_); + const double x_val = value_of_rec(x.val_); + if (std::fabs(x_val) < 0.5) { + return fvar(u, x.d_ * (1 - u * u)); + } + // -2 |x|, with the sign of x fixed in this branch + const T e = exp(x_val > 0 ? -2 * x.val_ : 2 * x.val_); + const T one_p_e = 1 + e; + return fvar(u, x.d_ * (4 * e / (one_p_e * one_p_e))); } /** diff --git a/stan/math/prim/constraint/cholesky_corr_constrain.hpp b/stan/math/prim/constraint/cholesky_corr_constrain.hpp index 8ed2b1aa454..a964b716387 100644 --- a/stan/math/prim/constraint/cholesky_corr_constrain.hpp +++ b/stan/math/prim/constraint/cholesky_corr_constrain.hpp @@ -3,77 +3,187 @@ #include #include -#include -#include -#include -#include +#include +#include +#include +#include +#include #include +#include namespace stan { namespace math { -template * = nullptr> +namespace internal { +/** + * Return the coefficients `c` for which the log absolute Jacobian + * determinant of `cholesky_corr_constrain(y, K)` is + * `-sum_n c_n log(cosh(y_n))`. + * + * With `z = tanh(y)`, row `i` holds `z_m sqrt(prod_{l < m} (1 - z_l^2))` in + * column `m < i`. The Jacobian has a factor `1 - z_m^2` for each element + * and a factor `sqrt(prod_{l < j} (1 - z_l^2))` for each column + * `0 < j < i`. With `log(1 - z_m^2) = -2 log(cosh(y_m))`, the coefficient + * of element `m` of row `i` is `2 + (i - 1 - m) = i + 1 - m`. + * + * @param K number of rows of the Cholesky factor + * @return coefficients in the order of the unconstrained vector + */ +inline Eigen::VectorXd cholesky_corr_log_jacobian_coefficients(int K) { + Eigen::VectorXd c((K * (K - 1)) / 2); + int k = 0; + for (int i = 1; i < K; ++i) { + for (int m = 0; m < i; ++m) { + c.coeffRef(k++) = i + 1 - m; + } + } + return c; +} + +/** + * Return `tanh(y)` and `sech(y) = 1 / cosh(y)` from one call of `expm1` or + * `exp`, with a relative error of a few ulp each. + * + * For `|y| < 1`, with `u = expm1(|y|)`, `cosh = 1 + u^2 / (2 (1 + u))` and + * `sinh = u (u + 2) / (2 (1 + u))`. For `|y| >= 1`, with `e = exp(-|y|)`, + * `tanh = (1 - e^2) / (1 + e^2)` and `sech = 2 e / (1 + e^2)`. + * + * @param y argument + * @return pair of `tanh(y)` and `sech(y)` + */ +inline std::pair tanh_and_sech(double y) { + const double a = std::fabs(y); + double tanh_a; + double sech; + if (a < 1.0) { + const double u = std::expm1(a); + const double cosh = 1.0 + u * u / (2.0 * (1.0 + u)); + tanh_a = u * (u + 2.0) / (2.0 * (1.0 + u)) / cosh; + sech = 1.0 / cosh; + } else { + const double e = std::exp(-a); + const double e2 = e * e; + tanh_a = (1.0 - e2) / (1.0 + e2); + sech = 2.0 * e / (1.0 + e2); + } + return {std::copysign(tanh_a, y), sech}; +} + +/** + * Return the position in the unconstrained vector of the partial + * correlation for element `(r, c)`, `c < r`, of the Cholesky factor. + * + * @tparam ColumnOrder `false` for the row order of + * `cholesky_corr_constrain` (row 1, then row 2, ...), `true` for the column + * order of `corr_matrix_constrain` (column 0, then column 1, ...) + * @param r row + * @param c column + * @param K number of rows + * @return position in the unconstrained vector + */ +template +inline Eigen::Index cholesky_corr_position(Eigen::Index r, Eigen::Index c, + Eigen::Index K) { + if constexpr (ColumnOrder) { + return c * (K - 1) - (c * (c - 1)) / 2 + (r - c - 1); + } else { + return (r * (r - 1)) / 2 + c; + } +} + +/** + * Return the Cholesky factor of a correlation matrix from unconstrained + * values `y`, whose partial correlations are `tanh(y)`. + * + * In row `r`, `1 - sum` of the squares of the elements before column `c` + * is `prod_{m < c} (1 - tanh(y_m)^2) = exp(-2 a)` with + * `a = sum_{m < c} log(cosh(y_m))`, so + * `L(r, c) = tanh(y_c) exp(-a)` and `L(r, r) = exp(-a)` over the whole + * row. Unlike `sqrt(1 - sum)`, this has no subtraction: the elements keep + * a small relative error when partial correlations are close to ±1 or a + * row is long, and the diagonal stays positive unless `a` exceeds about + * 745. + * + * @tparam ColumnOrder order of the partial correlations in `y`, see + * `cholesky_corr_position` + * @tparam EigVec type of the vector + * @param y unconstrained values, `K (K - 1) / 2` of them + * @param K number of rows + * @return Cholesky factor of a correlation matrix + */ +template * = nullptr> inline Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> -cholesky_corr_constrain(const EigVec& y, int K) { - using Eigen::Dynamic; - using Eigen::Matrix; - using std::sqrt; +cholesky_corr_factor(const EigVec& y, Eigen::Index K) { using T_scalar = value_type_t; - int k_choose_2 = (K * (K - 1)) / 2; - check_size_match("cholesky_corr_constrain", "constrain size", y.size(), - "k_choose_2", k_choose_2); - Matrix z = corr_constrain(y); - Matrix x(K, K); + Eigen::Matrix L + = Eigen::Matrix::Zero(K, K); if (K == 0) { - return x; + return L; } - x.setZero(); - x.coeffRef(0, 0) = 1; - int k = 0; - for (int i = 1; i < K; ++i) { - x.coeffRef(i, 0) = z.coeff(k++); - T_scalar sum_sqs = square(x.coeff(i, 0)); - for (int j = 1; j < i; ++j) { - x.coeffRef(i, j) = z.coeff(k++) * sqrt(1.0 - sum_sqs); - sum_sqs += square(x.coeff(i, j)); + const auto& y_ref = to_ref(y); + L.coeffRef(0, 0) = 1; + for (Eigen::Index r = 1; r < K; ++r) { + T_scalar a = 0; + for (Eigen::Index c = 0; c < r; ++c) { + const auto& y_rc + = y_ref.coeff(cholesky_corr_position(r, c, K)); + L.coeffRef(r, c) = tanh(y_rc) * exp(-a); + a += log_cosh(y_rc); } - x.coeffRef(i, i) = sqrt(1.0 - sum_sqs); + L.coeffRef(r, r) = exp(-a); } - return x; + return L; } +} // namespace internal -// FIXME to match above after debugged +/** + * Return the Cholesky factor of a `K x K` correlation matrix from the + * `K (K - 1) / 2` unconstrained values `y`, whose partial correlations are + * `tanh(y)` in row order. + * + * @tparam EigVec type of the vector + * @param y unconstrained values + * @param K number of rows + * @return Cholesky factor of a correlation matrix + */ +template * = nullptr> +inline Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> +cholesky_corr_constrain(const EigVec& y, int K) { + int k_choose_2 = (K * (K - 1)) / 2; + check_size_match("cholesky_corr_constrain", "constrain size", y.size(), + "k_choose_2", k_choose_2); + return internal::cholesky_corr_factor(y, K); +} + +/** + * Return the Cholesky factor of a `K x K` correlation matrix from the + * `K (K - 1) / 2` unconstrained values `y`, and add the log absolute + * Jacobian determinant of the transform to `lp`. + * + * @tparam EigVec type of the vector + * @tparam Lp scalar type of the log density accumulator + * @param y unconstrained values + * @param K number of rows + * @param[in,out] lp log density accumulator + * @return Cholesky factor of a correlation matrix + */ template * = nullptr, require_convertible_t, Lp>* = nullptr> inline Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> cholesky_corr_constrain(const EigVec& y, int K, Lp& lp) { - using Eigen::Dynamic; - using Eigen::Matrix; - using std::sqrt; - using T_scalar = value_type_t; int k_choose_2 = (K * (K - 1)) / 2; check_size_match("cholesky_corr_constrain", "y.size()", y.size(), "k_choose_2", k_choose_2); - Matrix z = corr_constrain(y, lp); - Matrix x(K, K); if (K == 0) { - return x; - } - x.setZero(); - x.coeffRef(0, 0) = 1; - int k = 0; - for (int i = 1; i < K; ++i) { - x.coeffRef(i, 0) = z.coeff(k++); - T_scalar sum_sqs = square(x.coeff(i, 0)); - for (int j = 1; j < i; ++j) { - lp += 0.5 * log1m(sum_sqs); - x.coeffRef(i, j) = z.coeff(k++) * sqrt(1.0 - sum_sqs); - sum_sqs += square(x.coeff(i, j)); - } - x.coeffRef(i, i) = sqrt(1.0 - sum_sqs); + return Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic>( + 0, 0); } - return x; + const auto& y_ref = to_ref(y); + lp -= dot_product(internal::cholesky_corr_log_jacobian_coefficients(K), + log_cosh(y_ref)); + return internal::cholesky_corr_factor(y_ref, K); } /** diff --git a/stan/math/prim/constraint/corr_constrain.hpp b/stan/math/prim/constraint/corr_constrain.hpp index c5cd1862139..dfc424b0a3b 100644 --- a/stan/math/prim/constraint/corr_constrain.hpp +++ b/stan/math/prim/constraint/corr_constrain.hpp @@ -2,8 +2,7 @@ #define STAN_MATH_PRIM_CONSTRAINT_CORR_CONSTRAIN_HPP #include -#include -#include +#include #include #include #include @@ -36,7 +35,11 @@ inline plain_type_t corr_constrain(T&& x) { * corr_constrain(T). The log absolute Jacobian * determinant is * - *

\f$\log | \frac{d}{dx} \tanh x | = \log (1 - \tanh^2 x)\f$. + *

\f$\log | \frac{d}{dx} \tanh x | = \log (1 - \tanh^2 x) + * = -2 \log \cosh x\f$. + * + * The difference \f$1 - \tanh^2 x\f$ has no correct digits from + * \f$|x| = 19\f$; `log_cosh` has a small relative error for every `x`. * * @tparam T_x Type of scalar or container * @param[in] x value or container @@ -44,9 +47,8 @@ inline plain_type_t corr_constrain(T&& x) { */ template inline auto corr_constrain(T_x&& x, T_lp& lp) { - plain_type_t tanh_x = tanh(std::forward(x)); - lp += sum(log1m(square(tanh_x))); - return tanh_x; + lp -= 2.0 * sum(log_cosh(x)); + return corr_constrain(std::forward(x)); } /** diff --git a/stan/math/prim/constraint/corr_matrix_constrain.hpp b/stan/math/prim/constraint/corr_matrix_constrain.hpp index a2d54131f6f..2b9c3cae8bd 100644 --- a/stan/math/prim/constraint/corr_matrix_constrain.hpp +++ b/stan/math/prim/constraint/corr_matrix_constrain.hpp @@ -4,7 +4,11 @@ #include #include #include +#include #include +#include +#include +#include #include #include #include @@ -12,6 +16,33 @@ namespace stan { namespace math { +namespace internal { +/** + * Return the coefficients `c` for which the log absolute Jacobian + * determinant of `corr_matrix_constrain(x, K)` is + * `-sum_n c_n log(cosh(x_n))`. + * + * The transform is `tanh` to canonical partial correlations, with + * `log(1 - tanh(x)^2) = -2 log(cosh(x))`, then `read_corr_matrix`, whose + * log Jacobian is `(K - j - 2) / 2 log(1 - z^2)` for a partial correlation + * `z` of column `j` (equation 11 of Lewandowski, Kurowicka and Joe, 2009). + * Column `j` holds `K - 1 - j` values, each with the coefficient + * `2 + (K - j - 2) = K - j`. + * + * @param K number of rows of the correlation matrix + * @return coefficients in the order of the unconstrained vector + */ +inline Eigen::VectorXd corr_matrix_log_jacobian_coefficients(Eigen::Index K) { + Eigen::VectorXd c((K * (K - 1)) / 2); + Eigen::Index pos = 0; + for (Eigen::Index j = 0; j < K - 1; ++j) { + c.segment(pos, K - 1 - j).setConstant(K - j); + pos += K - 1 - j; + } + return c; +} +} // namespace internal + /** * Return the correlation matrix of the specified dimensionality * derived from the specified vector of unconstrained values. The @@ -44,9 +75,11 @@ corr_matrix_constrain(T&& x, Eigen::Index k) { decltype(auto) x_ref = to_ref(std::forward(x)); check_size_match("cov_matrix_constrain", "x.size()", x_ref.size(), "k_choose_2", k_choose_2); - decltype(auto) cpc_ref - = to_ref(corr_constrain(std::forward(x_ref))); - return read_corr_matrix(cpc_ref, k); + if (k == 0) { + return {}; + } + return multiply_lower_tri_self_transpose( + internal::cholesky_corr_factor(x_ref, k)); } /** @@ -79,9 +112,13 @@ corr_matrix_constrain(T&& x, Eigen::Index k, Lp& lp) { decltype(auto) x_ref = to_ref(std::forward(x)); check_size_match("cov_matrix_constrain", "x.size()", x_ref.size(), "k_choose_2", k_choose_2); - decltype(auto) cpc_ref - = to_ref(corr_constrain(std::forward(x_ref), lp)); - return read_corr_matrix(cpc_ref, k, lp); + if (k == 0) { + return {}; + } + lp -= dot_product(internal::corr_matrix_log_jacobian_coefficients(k), + log_cosh(x_ref)); + return multiply_lower_tri_self_transpose( + internal::cholesky_corr_factor(x_ref, k)); } /** diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 8ef4204c81b..5fb8708e5e9 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -176,6 +176,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/log_cosh.hpp b/stan/math/prim/fun/log_cosh.hpp new file mode 100644 index 00000000000..cbf2f0be941 --- /dev/null +++ b/stan/math/prim/fun/log_cosh.hpp @@ -0,0 +1,92 @@ +#ifndef STAN_MATH_PRIM_FUN_LOG_COSH_HPP +#define STAN_MATH_PRIM_FUN_LOG_COSH_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +namespace internal { +/** + * Return `log(cosh(x))` and its derivative `tanh(x)`. + * + * For `|x| < 1`, with `u = expm1(|x|)`, `cosh(x) - 1 = u^2 / (2 (1 + u))` + * and `tanh(|x|) = u (u + 2) / (u (u + 2) + 2)`. For `|x| >= 1`, with + * `t = expm1(-2 |x|)`, `log(cosh(x)) = |x| + log1p(t / 2)` and + * `tanh(|x|) = -t / (2 + t)`. No branch subtracts two nearly equal numbers + * or overflows, so both results have a small relative error for every `x`. + * + * @param x argument + * @return pair of `log(cosh(x))` and `tanh(x)` + */ +inline std::pair log_cosh_and_tanh(double x) { + const double a = std::fabs(x); + double value; + double tanh_a; + if (a < 1.0) { + const double u = std::expm1(a); + const double w = u * (u + 2.0); + value = std::log1p(u * u / (2.0 * (1.0 + u))); + tanh_a = w / (w + 2.0); + } else { + const double t = std::expm1(-2.0 * a); + value = a + std::log1p(0.5 * t); + tanh_a = -t / (2.0 + t); + } + return {value, std::copysign(tanh_a, x)}; +} +} // namespace internal + +/** + * Return the natural logarithm of the hyperbolic cosine of the argument. + * + * The result does not overflow for large `|x|` and has a small relative + * error for small `|x|`. The derivative is `tanh(x)`. + * + * @tparam T arithmetic type + * @param x argument + * @return `log(cosh(x))` + */ +template * = nullptr> +inline double log_cosh(T x) { + return internal::log_cosh_and_tanh(static_cast(x)).first; +} + +/** + * Structure to wrap `log_cosh()` so it can be vectorized. + */ +struct log_cosh_fun { + /** + * Return the natural logarithm of the hyperbolic cosine of the argument. + * + * @tparam T type of argument + * @param x argument + * @return `log(cosh(x))` + */ + template + static inline auto fun(const T& x) { + return log_cosh(x); + } +}; + +/** + * Return the elementwise `log_cosh()` of the specified container. + * + * @tparam T type of container + * @param x container + * @return elementwise `log(cosh(x))` + */ +template * = nullptr, + require_not_var_matrix_t* = nullptr, + require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> +inline auto log_cosh(const T& x) { + return apply_scalar_unary::apply(x); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/constraint/cholesky_corr_constrain.hpp b/stan/math/rev/constraint/cholesky_corr_constrain.hpp index 576a1bb962e..dab6933210b 100644 --- a/stan/math/rev/constraint/cholesky_corr_constrain.hpp +++ b/stan/math/rev/constraint/cholesky_corr_constrain.hpp @@ -2,17 +2,79 @@ #define STAN_MATH_REV_CONSTRAINT_CHOLESKY_CORR_CONSTRAIN_HPP #include +#include +#include #include #include -#include -#include -#include -#include +#include +#include #include +#include namespace stan { namespace math { +namespace internal { +/** + * Return the Cholesky factor of a correlation matrix from unconstrained + * values `y`, whose partial correlations are `tanh(y)`. See the `prim` + * version for the formulas. + * + * Reverse pass: in row `r`, element `c` depends on its own `y` through + * `tanh(y) exp(-a)`, with derivative `exp(-2 log(cosh(y)) - a)`, and every + * later element of the row, the diagonal included, depends on it through + * `exp(-a)`, with derivative `-tanh(y)` times the element. + * + * @tparam ColumnOrder order of the partial correlations in `y`, see + * `cholesky_corr_position` + * @tparam T type of the `var_value` vector + * @param y unconstrained values, `K (K - 1) / 2` of them + * @param K number of rows + * @return Cholesky factor of a correlation matrix + */ +template * = nullptr> +inline var_value cholesky_corr_factor(const T& y, + Eigen::Index K) { + using Eigen::Index; + if (K == 0) { + return Eigen::MatrixXd(); + } + const Index n = y.size(); + arena_t tanh_y(n); + // derivative of each element with respect to its own unconstrained value + arena_t d_own(n); + arena_t L_val = Eigen::MatrixXd::Zero(K, K); + L_val.coeffRef(0, 0) = 1; + for (Index r = 1; r < K; ++r) { + // exp(-a), as the product of sech(y) over the earlier elements of the row + double exp_minus_a = 1.0; + for (Index c = 0; c < r; ++c) { + const Index k = cholesky_corr_position(r, c, K); + const std::pair t_s = tanh_and_sech(y.val().coeff(k)); + tanh_y.coeffRef(k) = t_s.first; + d_own.coeffRef(k) = t_s.second * t_s.second * exp_minus_a; + L_val.coeffRef(r, c) = t_s.first * exp_minus_a; + exp_minus_a *= t_s.second; + } + L_val.coeffRef(r, r) = exp_minus_a; + } + var_value L = L_val; + reverse_pass_callback([y, L, tanh_y, d_own, K]() mutable { + for (Index r = 1; r < K; ++r) { + double later = L.adj().coeff(r, r) * L.val().coeff(r, r); + for (Index c = r - 1; c >= 0; --c) { + const Index k = cholesky_corr_position(r, c, K); + const double adj_rc = L.adj().coeff(r, c); + y.adj().coeffRef(k) + += adj_rc * d_own.coeff(k) - tanh_y.coeff(k) * later; + later += adj_rc * L.val().coeff(r, c); + } + } + }); + return L; +} +} // namespace internal + /** * Return the Cholesky factor of the correlation matrix of the specified * size read from the unconstrained vector `y`. A total of K choose 2 @@ -26,54 +88,10 @@ namespace math { */ template * = nullptr> inline var_value cholesky_corr_constrain(const T& y, int K) { - using std::sqrt; int k_choose_2 = (K * (K - 1)) / 2; check_size_match("cholesky_corr_constrain", "y.size()", y.size(), "k_choose_2", k_choose_2); - - if (K == 0) { - return Eigen::MatrixXd(); - } - - var_value z = corr_constrain(y); - arena_t x_val = Eigen::MatrixXd::Zero(K, K); - - x_val.coeffRef(0, 0) = 1; - int k = 0; - for (int i = 1; i < K; ++i) { - x_val.coeffRef(i, 0) = z.val().coeff(k); - k++; - double sum_sqs = square(x_val.coeff(i, 0)); - for (int j = 1; j < i; ++j) { - x_val.coeffRef(i, j) = z.val().coeff(k) * sqrt(1.0 - sum_sqs); - k++; - sum_sqs += square(x_val.coeff(i, j)); - } - x_val.coeffRef(i, i) = sqrt(1.0 - sum_sqs); - } - - var_value x = x_val; - - reverse_pass_callback([z, K, x]() mutable { - size_t k = z.size(); - for (int i = K - 1; i > 0; --i) { - double sum_sqs_val = 1.0 - square(x.val().coeffRef(i, i)); - double sum_sqs_adj = -0.5 * x.adj().coeff(i, i) / x.val().coeff(i, i); - for (int j = i - 1; j > 0; --j) { - x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); - sum_sqs_val -= square(x.val().coeff(i, j)); - k--; - sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(k) - / sqrt(1.0 - sum_sqs_val); - z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); - } - x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); - k--; - z.adj().coeffRef(k) += x.adj().coeffRef(i, 0); - } - }); - - return x; + return internal::cholesky_corr_factor(y, K); } /** @@ -91,61 +109,15 @@ inline var_value cholesky_corr_constrain(const T& y, int K) { template * = nullptr> inline var_value cholesky_corr_constrain( const T& y, int K, scalar_type_t& lp) { - using Eigen::Dynamic; - using Eigen::Matrix; - using std::sqrt; int k_choose_2 = (K * (K - 1)) / 2; check_size_match("cholesky_corr_constrain", "y.size()", y.size(), "k_choose_2", k_choose_2); - if (K == 0) { return Eigen::MatrixXd(); } - - var_value z = corr_constrain(y, lp); - arena_t x_val = Eigen::MatrixXd::Zero(K, K); - - x_val.coeffRef(0, 0) = 1; - int k = 0; - double lp_val = 0.0; - for (int i = 1; i < K; ++i) { - x_val.coeffRef(i, 0) = z.val().coeff(k); - k++; - double sum_sqs = square(x_val.coeff(i, 0)); - for (int j = 1; j < i; ++j) { - lp_val += 0.5 * log1m(sum_sqs); - x_val.coeffRef(i, j) = z.val().coeff(k) * sqrt(1.0 - sum_sqs); - k++; - sum_sqs += square(x_val.coeff(i, j)); - } - x_val.coeffRef(i, i) = sqrt(1.0 - sum_sqs); - } - - lp += lp_val; - var_value x = x_val; - - reverse_pass_callback([z, K, x, lp]() mutable { - size_t k = z.size(); - for (int i = K - 1; i > 0; --i) { - double sum_sqs_val = 1.0 - square(x.val().coeffRef(i, i)); - double sum_sqs_adj = -0.5 * x.adj().coeff(i, i) / x.val().coeff(i, i); - for (int j = i - 1; j > 0; --j) { - x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); - sum_sqs_val -= square(x.val().coeff(i, j)); - - k--; - sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(k) - / sqrt(1.0 - sum_sqs_val); - z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); - sum_sqs_adj -= 0.5 * lp.adj() / (1 - sum_sqs_val); - } - x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); - k--; - z.adj().coeffRef(k) += x.adj().coeffRef(i, 0); - } - }); - - return x; + lp -= dot_product(internal::cholesky_corr_log_jacobian_coefficients(K), + log_cosh(y)); + return internal::cholesky_corr_factor(y, K); } } // namespace math diff --git a/stan/math/rev/constraint/corr_matrix_constrain.hpp b/stan/math/rev/constraint/corr_matrix_constrain.hpp index 8b8a028c9a5..d29898cabe3 100644 --- a/stan/math/rev/constraint/corr_matrix_constrain.hpp +++ b/stan/math/rev/constraint/corr_matrix_constrain.hpp @@ -2,10 +2,15 @@ #define STAN_MATH_REV_CONSTRAINT_CORR_MATRIX_CONSTRAIN_HPP #include +#include +#include +#include +#include #include #include #include #include +#include #include #include @@ -43,7 +48,11 @@ inline var_value corr_matrix_constrain(const T& x, Eigen::Index k_choose_2 = (k * (k - 1)) / 2; check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2", k_choose_2); - return read_corr_matrix(corr_constrain(x), k); + if (k == 0) { + return Eigen::MatrixXd(); + } + return multiply_lower_tri_self_transpose( + internal::cholesky_corr_factor(x, k)); } /** @@ -73,7 +82,13 @@ inline var_value corr_matrix_constrain(const T& x, Eigen::Index k_choose_2 = (k * (k - 1)) / 2; check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2", k_choose_2); - return read_corr_matrix(corr_constrain(x, lp), k, lp); + if (k == 0) { + return Eigen::MatrixXd(); + } + lp -= dot_product(internal::corr_matrix_log_jacobian_coefficients(k), + log_cosh(x)); + return multiply_lower_tri_self_transpose( + internal::cholesky_corr_factor(x, k)); } } // namespace math diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index e5ef938a302..ccb7f813778 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -105,6 +105,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/log_cosh.hpp b/stan/math/rev/fun/log_cosh.hpp new file mode 100644 index 00000000000..d3b820656c1 --- /dev/null +++ b/stan/math/rev/fun/log_cosh.hpp @@ -0,0 +1,55 @@ +#ifndef STAN_MATH_REV_FUN_LOG_COSH_HPP +#define STAN_MATH_REV_FUN_LOG_COSH_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the natural logarithm of the hyperbolic cosine of the argument. + * + * @tparam T arithmetic type + * @param x `var_value` with inner arithmetic type + * @return `log(cosh(x))`, with derivative `tanh(x)` + */ +template * = nullptr> +inline var log_cosh(const var_value& x) { + const std::pair value_and_tanh + = internal::log_cosh_and_tanh(x.val()); + const double tanh_x = value_and_tanh.second; + return make_callback_var(value_and_tanh.first, [x, tanh_x](auto& vi) mutable { + x.adj() += vi.adj() * tanh_x; + }); +} + +/** + * Return the elementwise natural logarithm of the hyperbolic cosine of the + * argument. + * + * @tparam T a type derived from `Eigen::EigenBase` + * @param x `var_value` with inner Eigen type + * @return elementwise `log(cosh(x))`, with derivative `tanh(x)` + */ +template * = nullptr> +inline auto log_cosh(const var_value& x) { + using val_t = plain_type_t; + val_t value(x.rows(), x.cols()); + arena_t tanh_x(x.rows(), x.cols()); + for (Eigen::Index i = 0; i < x.size(); ++i) { + const std::pair value_and_tanh + = internal::log_cosh_and_tanh(x.val().coeff(i)); + value.coeffRef(i) = value_and_tanh.first; + tanh_x.coeffRef(i) = value_and_tanh.second; + } + return make_callback_var(std::move(value), [x, tanh_x](auto& vi) mutable { + x.adj().array() += vi.adj().array() * tanh_x.array(); + }); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/mix/constraint/cholesky_corr_constrain_test.cpp b/test/unit/math/mix/constraint/cholesky_corr_constrain_test.cpp index 89e1ebb7e39..7618098e847 100644 --- a/test/unit/math/mix/constraint/cholesky_corr_constrain_test.cpp +++ b/test/unit/math/mix/constraint/cholesky_corr_constrain_test.cpp @@ -1,4 +1,5 @@ #include +#include namespace cholesky_corr_constrain_test { // easier than fiddling the quadratic equation @@ -114,3 +115,182 @@ TEST(mathMixMatFun, cholesky_corr_constrain_lp) { stan::test::expect_ad_matvar(f2(4), x1); stan::test::expect_ad_matvar(f2(3), x2); } + +namespace cholesky_corr_constrain_test { +// the strictly lower triangular elements of the Cholesky factor +struct strict_lower_of_constrained { + int K; + template + Eigen::Matrix operator()( + const Eigen::Matrix& y) const { + auto L = stan::math::cholesky_corr_constrain(y, K); + Eigen::Matrix v((K * (K - 1)) / 2); + int k = 0; + for (int j = 0; j < K; ++j) { + for (int i = j + 1; i < K; ++i) { + v(k++) = L(i, j); + } + } + return v; + } +}; +} // namespace cholesky_corr_constrain_test + +// lp is log |det J| of the map from y to the strictly lower triangular +// elements of the Cholesky factor. +TEST(mathMixMatFun, cholesky_corr_constrain_lp_is_log_det_jacobian) { + Eigen::VectorXd y(6); + y << 0.3, -1.2, 0.8, 2.1, -0.4, 1.5; + Eigen::VectorXd fx; + Eigen::MatrixXd J; + stan::math::jacobian( + cholesky_corr_constrain_test::strict_lower_of_constrained{4}, y, fx, J); + double lp = 0; + stan::math::cholesky_corr_constrain(y, 4, lp); + EXPECT_NEAR(std::log(std::fabs(J.determinant())), lp, 1e-12); +} + +// lp = sum_n w_n log(1 - tanh(y_n)^2) with weights (1, 1.5, 1) for K = 3, +// and its gradient, at large y against mpmath. Forming 1 - tanh(y)^2 and +// 1 - sum of squares loses digits from |y| = 8. +TEST(mathMixMatFun, cholesky_corr_constrain_lp_large_argument) { + using stan::math::fvar; + using stan::math::var; + Eigen::VectorXd y(3); + y << 15.0, 12.0, -0.5; + const double lp_ref = -62.77449311123027; + Eigen::VectorXd grad_ref(3); + grad_ref << -1.9999999999996257, -2.9999999997734919, 0.92423431452001952; + + double lp_d = 0; + stan::math::cholesky_corr_constrain(y, 3, lp_d); + EXPECT_NEAR(lp_ref, lp_d, 1e-13); + + Eigen::Matrix y_v = stan::math::to_var(y); + var lp_v = 0; + stan::math::cholesky_corr_constrain(y_v, 3, lp_v); + lp_v.grad(); + EXPECT_NEAR(lp_ref, lp_v.val(), 1e-13); + for (int i = 0; i < 3; ++i) { + EXPECT_NEAR(grad_ref(i), y_v(i).adj(), 1e-14); + } + stan::math::recover_memory(); + + stan::math::var_value y_m(y); + var lp_m = 0; + stan::math::cholesky_corr_constrain(y_m, 3, lp_m); + lp_m.grad(); + EXPECT_NEAR(lp_ref, lp_m.val(), 1e-13); + for (int i = 0; i < 3; ++i) { + EXPECT_NEAR(grad_ref(i), y_m.adj()(i), 1e-14); + } + stan::math::recover_memory(); + + for (int n = 0; n < 3; ++n) { + Eigen::Matrix, Eigen::Dynamic, 1> y_f(3); + for (int i = 0; i < 3; ++i) { + y_f(i) = fvar(y(i), i == n ? 1.0 : 0.0); + } + fvar lp_f = 0; + stan::math::cholesky_corr_constrain(y_f, 3, lp_f); + EXPECT_NEAR(lp_ref, lp_f.val_, 1e-13); + EXPECT_NEAR(grad_ref(n), lp_f.d_, 1e-14); + } +} + +namespace cholesky_corr_constrain_test { +// lower triangle (diagonal included), row by row +template +Eigen::VectorXd lower_triangle(const M& L) { + const int K = L.rows(); + Eigen::VectorXd v((K * (K + 1)) / 2); + int k = 0; + for (int i = 0; i < K; ++i) { + for (int j = 0; j <= i; ++j) { + v(k++) = stan::math::value_of_rec(L(i, j)); + } + } + return v; +} +} // namespace cholesky_corr_constrain_test + +// Values of L at large |y| against mpmath. sqrt(1 - sum of squares) had a +// relative error of 3.8e-4 in L(1, 1) at y = 15. +TEST(mathMixMatFun, cholesky_corr_constrain_values_large_argument) { + using cholesky_corr_constrain_test::lower_triangle; + using stan::math::fvar; + using stan::math::var; + Eigen::VectorXd y(6); + y << 15.0, 12.0, -0.5, 9.0, 2.0, -14.0; + Eigen::VectorXd ref(10); + ref << 1.0, 0.99999999999981285, 6.1180464100359433e-7, 0.99999999992449731, + -5.6786918924293556e-6, 1.0897607083695929e-5, 0.99999996954004097, + 0.00023794090595878951, -6.5605200973176475e-5, 1.091052174635906e-10; + auto expect_ref = [&ref](const Eigen::VectorXd& v, const char* mode) { + for (int i = 0; i < ref.size(); ++i) { + EXPECT_NEAR(ref(i), v(i), 4e-15 * std::fabs(ref(i))) + << mode << ", element " << i; + } + }; + expect_ref(lower_triangle(stan::math::cholesky_corr_constrain(y, 4)), + "double"); + Eigen::Matrix y_v = stan::math::to_var(y); + expect_ref(lower_triangle(stan::math::cholesky_corr_constrain(y_v, 4)), + "var"); + stan::math::var_value y_m(y); + expect_ref( + lower_triangle(stan::math::cholesky_corr_constrain(y_m, 4).val().eval()), + "var_value"); + Eigen::Matrix, Eigen::Dynamic, 1> y_f(6); + for (int i = 0; i < 6; ++i) { + y_f(i) = fvar(y(i), 1.0); + } + expect_ref(lower_triangle(stan::math::cholesky_corr_constrain(y_f, 4)), + "fvar"); + stan::math::recover_memory(); +} + +// K = 40 with unconstrained values in [-1.9, 1.9], as default inits: a row +// of 39 moderate values makes 1 - sum of squares about 1e-20, which +// sqrt(1 - sum) rounded to 0. The diagonal must stay positive, the rows +// must have unit norm, lp and the LKJ density must be finite, and the +// var_value gradient must equal the Eigen gradient. +TEST(mathMixMatFun, cholesky_corr_constrain_long_rows) { + using stan::math::var; + const int K = 40; + Eigen::VectorXd y((K * (K - 1)) / 2); + for (int i = 0; i < y.size(); ++i) { + y(i) = 1.9 * std::sin(1.3 * i + 0.7); + } + double lp = 0; + Eigen::MatrixXd L = stan::math::cholesky_corr_constrain(y, K, lp); + EXPECT_GT(L.diagonal().minCoeff(), 0.0); + EXPECT_LT((L.rowwise().squaredNorm().array() - 1).abs().maxCoeff(), 1e-14); + EXPECT_TRUE(std::isfinite(lp)); + EXPECT_TRUE(std::isfinite(stan::math::lkj_corr_cholesky_lpdf(L, 1.0))); + + // gradient of 2 sum(log(diag(L))) + lp: Eigen through the prim code, + // var_value through the hand-written reverse pass + Eigen::Matrix y_v = stan::math::to_var(y); + var lp_v = 0; + Eigen::Matrix L_v + = stan::math::cholesky_corr_constrain(y_v, K, lp_v); + var f_v = 2.0 * stan::math::sum(stan::math::log(L_v.diagonal())) + lp_v; + f_v.grad(); + Eigen::VectorXd g_v = stan::math::value_of(y_v.adj()); + stan::math::recover_memory(); + + stan::math::var_value y_m(y); + var lp_m = 0; + stan::math::var_value L_m + = stan::math::cholesky_corr_constrain(y_m, K, lp_m); + var f_m = 2.0 * stan::math::sum(stan::math::log(L_m.diagonal())) + lp_m; + f_m.grad(); + Eigen::VectorXd g_m = y_m.adj(); + stan::math::recover_memory(); + + EXPECT_TRUE(g_v.allFinite()); + EXPECT_TRUE(g_m.allFinite()); + EXPECT_LT((g_v - g_m).cwiseAbs().maxCoeff(), + 1e-12 * g_v.cwiseAbs().maxCoeff()); +} diff --git a/test/unit/math/mix/constraint/corr_constrain_test.cpp b/test/unit/math/mix/constraint/corr_constrain_test.cpp index 609062f9f73..73f19c0bbd1 100644 --- a/test/unit/math/mix/constraint/corr_constrain_test.cpp +++ b/test/unit/math/mix/constraint/corr_constrain_test.cpp @@ -1,5 +1,52 @@ #include +// log(1 - tanh(x)^2) and its derivative -2 tanh(x) at large |x|, against +// mpmath. Forming 1 - tanh(x)^2 loses all digits from |x| = 19.06. +TEST(mathMixMatFun, corr_constrain_lp_large_argument) { + using stan::math::fvar; + using stan::math::var; + Eigen::VectorXd x(3); + x << 15.0, -19.5, 30.0; + const double lp_ref = -124.84111691664052; + Eigen::VectorXd grad_ref(3); + grad_ref << -1.9999999999996257, 2.0, -2.0; + + double lp_d = 0; + stan::math::corr_constrain(x, lp_d); + EXPECT_NEAR(lp_ref, lp_d, 1e-13); + + Eigen::Matrix x_v = stan::math::to_var(x); + var lp_v = 0; + stan::math::corr_constrain(x_v, lp_v); + lp_v.grad(); + EXPECT_NEAR(lp_ref, lp_v.val(), 1e-13); + for (int i = 0; i < 3; ++i) { + EXPECT_NEAR(grad_ref(i), x_v(i).adj(), 1e-14); + } + stan::math::recover_memory(); + + stan::math::var_value x_m(x); + var lp_m = 0; + stan::math::corr_constrain(x_m, lp_m); + lp_m.grad(); + EXPECT_NEAR(lp_ref, lp_m.val(), 1e-13); + for (int i = 0; i < 3; ++i) { + EXPECT_NEAR(grad_ref(i), x_m.adj()(i), 1e-14); + } + stan::math::recover_memory(); + + for (int n = 0; n < 3; ++n) { + Eigen::Matrix, Eigen::Dynamic, 1> x_f(3); + for (int i = 0; i < 3; ++i) { + x_f(i) = fvar(x(i), i == n ? 1.0 : 0.0); + } + fvar lp_f = 0; + stan::math::corr_constrain(x_f, lp_f); + EXPECT_NEAR(lp_ref, lp_f.val_, 1e-13); + EXPECT_NEAR(grad_ref(n), lp_f.d_, 1e-14); + } +} + TEST(mathMixMatFun, corr_constrain) { auto f = [](const auto& x1) { stan::scalar_type_t lp = 0.0; diff --git a/test/unit/math/mix/constraint/corr_matrix_constrain_test.cpp b/test/unit/math/mix/constraint/corr_matrix_constrain_test.cpp index 2b766dd7660..b0539788946 100644 --- a/test/unit/math/mix/constraint/corr_matrix_constrain_test.cpp +++ b/test/unit/math/mix/constraint/corr_matrix_constrain_test.cpp @@ -1,4 +1,5 @@ #include +#include namespace corr_matrix_constrain_test { // easier than fiddling the quadratic equation @@ -117,3 +118,176 @@ TEST(mathMixMatFun, corr_matrix_constrain_lp) { stan::test::expect_ad_matvar(f2(4), x1); stan::test::expect_ad_matvar(f2(3), x2); } + +namespace corr_matrix_constrain_test { +// the strictly lower triangular elements of the correlation matrix +struct strict_lower_of_constrained { + int K; + template + Eigen::Matrix operator()( + const Eigen::Matrix& y) const { + auto Omega = stan::math::corr_matrix_constrain(y, K); + Eigen::Matrix v((K * (K - 1)) / 2); + int k = 0; + for (int j = 0; j < K; ++j) { + for (int i = j + 1; i < K; ++i) { + v(k++) = Omega(i, j); + } + } + return v; + } +}; +} // namespace corr_matrix_constrain_test + +// lp is log |det J| of the map from y to the strictly lower triangular +// elements of the correlation matrix. +TEST(mathMixMatFun, corr_matrix_constrain_lp_is_log_det_jacobian) { + Eigen::VectorXd y(6); + y << 0.3, -1.2, 0.8, 2.1, -0.4, 1.5; + Eigen::VectorXd fx; + Eigen::MatrixXd J; + stan::math::jacobian( + corr_matrix_constrain_test::strict_lower_of_constrained{4}, y, fx, J); + double lp = 0; + stan::math::corr_matrix_constrain(y, 4, lp); + EXPECT_NEAR(std::log(std::fabs(J.determinant())), lp, 1e-12); +} + +// lp = sum_n w_n log(1 - tanh(y_n)^2) with weights (1.5, 1.5, 1) for K = 3, +// and its gradient, at large y against mpmath. Forming 1 - tanh(y)^2 loses +// digits from |y| = 8. +TEST(mathMixMatFun, corr_matrix_constrain_lp_large_argument) { + using stan::math::fvar; + using stan::math::var; + Eigen::VectorXd y(3); + y << 15.0, 12.0, -0.5; + const double lp_ref = -77.081345930670418; + Eigen::VectorXd grad_ref(3); + grad_ref << -2.9999999999994385, -2.9999999997734919, 0.92423431452001952; + + double lp_d = 0; + stan::math::corr_matrix_constrain(y, 3, lp_d); + EXPECT_NEAR(lp_ref, lp_d, 1e-13); + + Eigen::Matrix y_v = stan::math::to_var(y); + var lp_v = 0; + stan::math::corr_matrix_constrain(y_v, 3, lp_v); + lp_v.grad(); + EXPECT_NEAR(lp_ref, lp_v.val(), 1e-13); + for (int i = 0; i < 3; ++i) { + EXPECT_NEAR(grad_ref(i), y_v(i).adj(), 1e-14); + } + stan::math::recover_memory(); + + stan::math::var_value y_m(y); + var lp_m = 0; + stan::math::corr_matrix_constrain(y_m, 3, lp_m); + lp_m.grad(); + EXPECT_NEAR(lp_ref, lp_m.val(), 1e-13); + for (int i = 0; i < 3; ++i) { + EXPECT_NEAR(grad_ref(i), y_m.adj()(i), 1e-14); + } + stan::math::recover_memory(); + + for (int n = 0; n < 3; ++n) { + Eigen::Matrix, Eigen::Dynamic, 1> y_f(3); + for (int i = 0; i < 3; ++i) { + y_f(i) = fvar(y(i), i == n ? 1.0 : 0.0); + } + fvar lp_f = 0; + stan::math::corr_matrix_constrain(y_f, 3, lp_f); + EXPECT_NEAR(lp_ref, lp_f.val_, 1e-13); + EXPECT_NEAR(grad_ref(n), lp_f.d_, 1e-14); + } +} + +namespace corr_matrix_constrain_test { +// strictly lower triangle, row by row +template +Eigen::VectorXd strict_lower(const M& Omega) { + const int K = Omega.rows(); + Eigen::VectorXd v((K * (K - 1)) / 2); + int k = 0; + for (int i = 0; i < K; ++i) { + for (int j = 0; j < i; ++j) { + v(k++) = stan::math::value_of_rec(Omega(i, j)); + } + } + return v; +} +} // namespace corr_matrix_constrain_test + +// Correlations at large |y| against mpmath. With sqrt(1 - z^2) from a +// rounded z, the small elements of the Cholesky factor had relative errors +// up to 3.8e-4, which reached the correlations close to 1. +TEST(mathMixMatFun, corr_matrix_constrain_values_large_argument) { + using corr_matrix_constrain_test::strict_lower; + using stan::math::fvar; + using stan::math::var; + Eigen::VectorXd y(6); + y << 15.0, 12.0, -0.5, 9.0, 2.0, -14.0; + Eigen::VectorXd ref(6); + ref << 0.99999999999981285, 0.99999999992449731, 0.99999999993182827, + -0.46211715726000976, -0.46211663421720721, -0.46210665234659286; + auto expect_ref = [&ref](const Eigen::VectorXd& v, const char* mode) { + for (int i = 0; i < ref.size(); ++i) { + EXPECT_NEAR(ref(i), v(i), 4e-16) << mode << ", element " << i; + } + }; + expect_ref(strict_lower(stan::math::corr_matrix_constrain(y, 4)), "double"); + Eigen::Matrix y_v = stan::math::to_var(y); + expect_ref(strict_lower(stan::math::corr_matrix_constrain(y_v, 4)), "var"); + stan::math::var_value y_m(y); + expect_ref( + strict_lower(stan::math::corr_matrix_constrain(y_m, 4).val().eval()), + "var_value"); + Eigen::Matrix, Eigen::Dynamic, 1> y_f(6); + for (int i = 0; i < 6; ++i) { + y_f(i) = fvar(y(i), 1.0); + } + expect_ref(strict_lower(stan::math::corr_matrix_constrain(y_f, 4)), "fvar"); + stan::math::recover_memory(); +} + +// K = 40 with unconstrained values in [-1.9, 1.9], as default inits: the +// diagonal must be 1, lp finite, and the var_value gradient must equal the +// Eigen gradient and be finite. +TEST(mathMixMatFun, corr_matrix_constrain_long_rows) { + using stan::math::var; + const int K = 40; + Eigen::VectorXd y((K * (K - 1)) / 2); + for (int i = 0; i < y.size(); ++i) { + y(i) = 1.9 * std::sin(1.3 * i + 0.7); + } + Eigen::MatrixXd W(K, K); + for (int i = 0; i < W.size(); ++i) { + W(i) = std::cos(0.7 * i); + } + double lp = 0; + Eigen::MatrixXd Omega = stan::math::corr_matrix_constrain(y, K, lp); + EXPECT_LT((Omega.diagonal().array() - 1).abs().maxCoeff(), 1e-14); + EXPECT_TRUE(std::isfinite(lp)); + + Eigen::Matrix y_v = stan::math::to_var(y); + var lp_v = 0; + Eigen::Matrix Omega_v + = stan::math::corr_matrix_constrain(y_v, K, lp_v); + var f_v = (W.array() * Omega_v.array()).sum() + lp_v; + f_v.grad(); + Eigen::VectorXd g_v = stan::math::value_of(y_v.adj()); + stan::math::recover_memory(); + + stan::math::var_value y_m(y); + var lp_m = 0; + stan::math::var_value Omega_m + = stan::math::corr_matrix_constrain(y_m, K, lp_m); + var f_m = stan::math::sum(stan::math::elt_multiply(W, Omega_m)) + lp_m; + f_m.grad(); + Eigen::VectorXd g_m = y_m.adj(); + stan::math::recover_memory(); + + EXPECT_TRUE(g_v.allFinite()); + EXPECT_TRUE(g_m.allFinite()); + EXPECT_LT((g_v - g_m).cwiseAbs().maxCoeff(), + 1e-12 * g_v.cwiseAbs().maxCoeff()); +} diff --git a/test/unit/math/mix/fun/log_cosh_test.cpp b/test/unit/math/mix/fun/log_cosh_test.cpp new file mode 100644 index 00000000000..0eb052845fe --- /dev/null +++ b/test/unit/math/mix/fun/log_cosh_test.cpp @@ -0,0 +1,63 @@ +#include +#include +#include +#include +#include + +TEST(mathMixMatFun, logCosh) { + auto f = [](const auto& x1) { return stan::math::log_cosh(x1); }; + stan::test::expect_common_unary_vectorized( + f); + stan::test::expect_unary_vectorized( + f, -30.0, -2.6, -1.0, -0.5, 0.0, 0.5, 0.999, 1.0, 1.001, 2.6, 30.0); +} + +TEST(mathMixMatFun, logCoshVarMat) { + using stan::math::vec_concat; + using stan::test::expect_ad_vector_matvar; + using stan::test::internal::common_nonzero_args; + auto f = [](const auto& x1) { return stan::math::log_cosh(x1); }; + std::vector com_args = common_nonzero_args(); + std::vector args{-30.0, -2.6, -0.5, 0.999, 1.001, 2.6, 30.0}; + auto all_args = vec_concat(com_args, args); + Eigen::VectorXd A(all_args.size()); + for (int i = 0; i < all_args.size(); ++i) { + A(i) = all_args[i]; + } + expect_ad_vector_matvar(f, A); +} + +// Values and derivatives against mpmath, on both sides of the branch at +// |x| = 1, at a tiny argument (where log(cosh(x)) = x^2 / 2), and where +// cosh(x) overflows. +TEST(mathMixMatFun, logCoshReference) { + using stan::math::fvar; + using stan::math::var; + // x, log(cosh(x)), tanh(x) + std::vector> cases + = {{1e-10, 5.0e-21, 1.0e-10}, + {0.5, 0.12011450695827752, 0.46211715726000976}, + {0.999, 0.43301944642088479, 0.76117386166056995}, + {1.0, 0.43378083048302719, 0.76159415595576489}, + {-3.0, 2.3093285045777851, -0.99505475368673045}, + {30.0, 29.306852819440055, 1.0}, + {800.0, 799.30685281944005, 1.0}}; + for (const auto& c : cases) { + const double x = c[0]; + const double tol_v = 4 * std::numeric_limits::epsilon() * c[1]; + const double tol_d + = 4 * std::numeric_limits::epsilon() * std::fabs(c[2]); + EXPECT_NEAR(c[1], stan::math::log_cosh(x), tol_v) << "x = " << x; + + var x_v = x; + var y_v = stan::math::log_cosh(x_v); + y_v.grad(); + EXPECT_NEAR(c[1], y_v.val(), tol_v) << "x = " << x; + EXPECT_NEAR(c[2], x_v.adj(), tol_d) << "x = " << x; + stan::math::recover_memory(); + + fvar y_f = stan::math::log_cosh(fvar(x, 1.0)); + EXPECT_NEAR(c[1], y_f.val_, tol_v) << "x = " << x; + EXPECT_NEAR(c[2], y_f.d_, tol_d) << "x = " << x; + } +} diff --git a/test/unit/math/mix/fun/tanh_test.cpp b/test/unit/math/mix/fun/tanh_test.cpp index d9e4ad5d567..9a894406fb9 100644 --- a/test/unit/math/mix/fun/tanh_test.cpp +++ b/test/unit/math/mix/fun/tanh_test.cpp @@ -1,4 +1,7 @@ #include +#include +#include +#include TEST(mathMixMatFun, tanh) { auto f = [](const auto& x1) { @@ -13,6 +16,71 @@ TEST(mathMixMatFun, tanh) { stan::test::expect_complex_common(f); } +// First and second derivatives at large |x| against mpmath. 1 - tanh(x)^2 +// is 0 in double from |x| = 19.06, and 1 / cosh(x)^2 makes the second +// derivative NaN once cosh(x)^2 overflows (|x| > 355). +TEST(mathMixMatFun, tanh_derivatives_large_argument) { + using stan::math::fvar; + using stan::math::var; + // x, d/dx tanh(x), d2/dx2 tanh(x) + std::vector> cases + = {{3.0, 0.0098660371654401913, -0.019634494363042436}, + {-20.0, 1.6993417021166356e-17, 3.3986834042332711e-17}, + {20.0, 1.6993417021166356e-17, -3.3986834042332711e-17}, + {300.0, 1.0601586212017243e-260, -2.1203172424034487e-260}, + {1000.0, 0.0, 0.0}}; + for (const auto& c : cases) { + const double x = c[0]; + const double d1 = c[1]; + const double d2 = c[2]; + + fvar y1 = stan::math::tanh(fvar(x, 1.0)); + EXPECT_NEAR(d1, y1.d_, 1e-14 * std::fabs(d1)) << "x = " << x; + + fvar> y2 = stan::math::tanh( + fvar>(fvar(x, 1.0), fvar(1.0, 0.0))); + EXPECT_NEAR(d1, y2.d_.val_, 1e-14 * std::fabs(d1)) << "x = " << x; + EXPECT_NEAR(d2, y2.d_.d_, 1e-14 * std::fabs(d2)) << "x = " << x; + + fvar x3(var(x), 1.0); + fvar y3 = stan::math::tanh(x3); + y3.d_.grad(); + EXPECT_NEAR(d1, y3.d_.val(), 1e-14 * std::fabs(d1)) << "x = " << x; + EXPECT_NEAR(d2, x3.val_.adj(), 1e-14 * std::fabs(d2)) << "x = " << x; + stan::math::recover_memory(); + } +} + +// First, second and third derivatives in fvar> against mpmath, +// at 0 (where an fabs(x) in the derivative would make the third derivative +// 0), on both sides of 0.5, and in the tail. +TEST(mathMixMatFun, tanh_third_derivative) { + using stan::math::fvar; + using stan::math::var; + // x, tanh'(x), tanh''(x), tanh'''(x) + std::vector> cases = { + {0.0, 1.0, 0.0, -2.0}, + {0.25, 0.94001484880637796, -0.46045435881856592, -1.5417081006333538}, + {0.4999, 0.78652041633736061, -0.72680544069311918, -0.56560452654352799}, + {0.5, 0.78644773296592741, -0.72686198138358728, -0.56520928825977036}, + {-0.5, 0.78644773296592741, 0.72686198138358728, -0.56520928825977036}, + {3.0, 0.0098660371654401913, -0.019634494363042436, 0.038880116525661682}, + {-20.0, 1.6993417021166356e-17, 3.3986834042332711e-17, + 6.7973668084665422e-17}}; + for (const auto& c : cases) { + const double x = c[0]; + var x_v = x; + fvar> x_ffv(fvar(x_v, 1.0), fvar(1.0, 0.0)); + fvar> y = stan::math::tanh(x_ffv); + y.d_.d_.grad(); + EXPECT_NEAR(c[1], y.d_.val_.val(), 1e-14 * std::fabs(c[1])) << "x = " << x; + EXPECT_NEAR(c[2], y.d_.d_.val(), 1e-14 * std::fabs(c[2]) + 1e-300) + << "x = " << x; + EXPECT_NEAR(c[3], x_v.adj(), 1e-14 * std::fabs(c[3])) << "x = " << x; + stan::math::recover_memory(); + } +} + TEST(mathMixMatFun, tanh_varmat) { using stan::math::vec_concat; using stan::test::expect_ad_vector_matvar; diff --git a/test/unit/math/rev/prob/lkj_corr_test.cpp b/test/unit/math/rev/prob/lkj_corr_test.cpp index 6c912bd949a..c56c7972bf0 100644 --- a/test/unit/math/rev/prob/lkj_corr_test.cpp +++ b/test/unit/math/rev/prob/lkj_corr_test.cpp @@ -65,7 +65,10 @@ TEST_F(AgradRev, ProbDistributionsLkjCorrCholesky_gradients1) { double fx_ad_1; finite_diff_gradient(test_func_1, x1, fx_1, grad_1); gradient(test_func_1, x1, fx_ad_1, grad_ad_1); - test_grad_eq(grad_1, grad_ad_1); + // x1(0) = 0, where the exact derivative is 0; the finite difference has + // rounding noise of about 1e-13 there, so compare with an absolute + // tolerance instead of float equality + test_grad_near(grad_1, grad_ad_1, 1e-10); EXPECT_FLOAT_EQ(fx_1, fx_ad_1); }