Loading [MathJax]/extensions/TeX/AMSmath.js
Automatic Differentiation
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
cov_matrix_constrain.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_CONSTRAINT_COV_MATRIX_CONSTRAIN_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_COV_MATRIX_CONSTRAIN_HPP
3
11#include <cmath>
12
13namespace stan {
14namespace math {
15
30template <typename T, require_eigen_col_vector_t<T>* = nullptr>
31inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
32cov_matrix_constrain(const T& x, Eigen::Index K) {
33 using Eigen::Dynamic;
34 using Eigen::Matrix;
35 using std::exp;
36
37 Matrix<value_type_t<T>, Dynamic, Dynamic> L(K, K);
38 check_size_match("cov_matrix_constrain", "x.size()", x.size(),
39 "K + (K choose 2)", (K * (K + 1)) / 2);
40 const auto& x_ref = to_ref(x);
41 int i = 0;
42 for (Eigen::Index m = 0; m < K; ++m) {
43 L.row(m).head(m) = x_ref.segment(i, m);
44 i += m;
45 L.coeffRef(m, m) = exp(x_ref.coeff(i++));
46 L.row(m).tail(K - m - 1).setZero();
47 }
49}
50
65template <typename T, typename Lp, require_eigen_col_vector_t<T>* = nullptr,
66 require_convertible_t<return_type_t<T>, Lp>* = nullptr>
67inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
68cov_matrix_constrain(const T& x, Eigen::Index K, Lp& lp) {
69 using Eigen::Dynamic;
70 using Eigen::Matrix;
71 using std::exp;
72 using std::log;
73 check_size_match("cov_matrix_constrain", "x.size()", x.size(),
74 "K + (K choose 2)", (K * (K + 1)) / 2);
75 Matrix<value_type_t<T>, Dynamic, Dynamic> L(K, K);
76 const auto& x_ref = to_ref(x);
77 int i = 0;
78 for (Eigen::Index m = 0; m < K; ++m) {
79 L.row(m).head(m) = x_ref.segment(i, m);
80 i += m;
81 L.coeffRef(m, m) = exp(x_ref.coeff(i++));
82 L.row(m).tail(K - m - 1).setZero();
83 }
84 // Jacobian for complete transform, including exp() above
85 lp += (K * LOG_TWO); // needless constant; want propto
86 for (Eigen::Index k = 0; k < K; ++k) {
87 lp += (K - k + 1) * log(L.coeff(k, k)); // only +1 because index from 0
88 }
90}
91
104template <typename T, require_std_vector_t<T>* = nullptr>
105inline auto cov_matrix_constrain(const T& x, Eigen::Index K) {
107 x, [K](auto&& v) { return cov_matrix_constrain(v, K); });
108}
109
125template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
126 require_convertible_t<return_type_t<T>, Lp>* = nullptr>
127inline auto cov_matrix_constrain(const T& x, Eigen::Index K, Lp& lp) {
129 x, [&lp, K](auto&& v) { return cov_matrix_constrain(v, K, lp); });
130}
131
152template <bool Jacobian, typename T, typename Lp,
154inline auto cov_matrix_constrain(const T& x, Eigen::Index K, Lp& lp) {
155 if constexpr (Jacobian) {
156 return cov_matrix_constrain(x, K, lp);
157 } else {
158 return cov_matrix_constrain(x, K);
159 }
160}
161
162} // namespace math
163} // namespace stan
164
165#endif
require_t< std::is_convertible< std::decay_t< T >, std::decay_t< S > > > require_convertible_t
Require types T and S satisfies std::is_convertible.
Eigen::Matrix< value_type_t< EigMat >, EigMat::RowsAtCompileTime, EigMat::RowsAtCompileTime > multiply_lower_tri_self_transpose(const EigMat &m)
Eigen::Matrix< value_type_t< T >, Eigen::Dynamic, Eigen::Dynamic > cov_matrix_constrain(const T &x, Eigen::Index K)
Return the symmetric, positive-definite matrix of dimensions K by K resulting from transforming the s...
fvar< T > log(const fvar< T > &x)
Definition log.hpp:18
static constexpr double LOG_TWO
The natural logarithm of 2, .
Definition constants.hpp:80
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:15
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...