Stan Math Library  2.9.0
reverse mode automatic differentiation
cholesky_factor_constrain.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_CHOLESKY_FACTOR_CONSTRAIN_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_CHOLESKY_FACTOR_CONSTRAIN_HPP
3 
6 #include <cmath>
7 #include <stdexcept>
8 #include <vector>
9 
10 namespace stan {
11 
12  namespace math {
13 
14  // CHOLESKY FACTOR
15 
27  template <typename T>
28  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
29  cholesky_factor_constrain(const Eigen::Matrix<T, Eigen::Dynamic, 1>& x,
30  int M,
31  int N) {
32  using std::exp;
33  if (M < N)
34  throw std::domain_error("cholesky_factor_constrain: "
35  "num rows must be >= num cols");
36  if (x.size() != ((N * (N + 1)) / 2 + (M - N) * N))
37  throw std::domain_error("cholesky_factor_constrain: x.size() must"
38  " be (N * (N + 1)) / 2 + (M - N) * N");
39  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> y(M, N);
40  T zero(0);
41  int pos = 0;
42  // upper square
43  for (int m = 0; m < N; ++m) {
44  for (int n = 0; n < m; ++n)
45  y(m, n) = x(pos++);
46  y(m, m) = exp(x(pos++));
47  for (int n = m + 1; n < N; ++n)
48  y(m, n) = zero;
49  }
50  // lower rectangle
51  for (int m = N; m < M; ++m)
52  for (int n = 0; n < N; ++n)
53  y(m, n) = x(pos++);
54  return y;
55  }
56 
71  template <typename T>
72  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
73  cholesky_factor_constrain(const Eigen::Matrix<T, Eigen::Dynamic, 1>& x,
74  int M,
75  int N,
76  T& lp) {
77  // cut-and-paste from above, so checks twice
78 
79  using stan::math::sum;
80  if (x.size() != ((N * (N + 1)) / 2 + (M - N) * N))
81  throw std::domain_error("cholesky_factor_constrain: x.size() "
82  "must be (k choose 2) + k");
83  int pos = 0;
84  std::vector<T> log_jacobians(N);
85  for (int n = 0; n < N; ++n) {
86  pos += n;
87  log_jacobians[n] = x(pos++);
88  }
89  lp += sum(log_jacobians); // optimized for autodiff vs. direct lp +=
90  return cholesky_factor_constrain(x, M, N);
91  }
92 
93 
94  }
95 
96 }
97 
98 #endif
fvar< T > sum(const std::vector< fvar< T > > &m)
Return the sum of the entries of the specified standard vector.
Definition: sum.hpp:20
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > cholesky_factor_constrain(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &x, int M, int N)
Return the Cholesky factor of the specified size read from the specified vector.
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
void domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.

     [ Stan Home Page ] © 2011–2015, Stan Development Team.