Stan  2.14.0
probability, sampling & optimization
normal_fullrank.hpp
Go to the documentation of this file.
1 #ifndef STAN_VARIATIONAL_NORMAL_FULLRANK_HPP
2 #define STAN_VARIATIONAL_NORMAL_FULLRANK_HPP
3 
4 #include <stan/math/prim/mat.hpp>
7 #include <algorithm>
8 #include <ostream>
9 #include <vector>
10 
11 namespace stan {
12 
13  namespace variational {
14 
19  class normal_fullrank : public base_family {
20  private:
24  Eigen::VectorXd mu_;
25 
30  Eigen::MatrixXd L_chol_;
31 
35  const int dimension_;
36 
45  void validate_mean(const char* function,
46  const Eigen::VectorXd& mu) {
47  stan::math::check_not_nan(function, "Mean vector", mu);
48  stan::math::check_size_match(function,
49  "Dimension of input vector", mu.size(),
50  "Dimension of current vector", dimension_);
51  }
52 
67  void validate_cholesky_factor(const char* function,
68  const Eigen::MatrixXd& L_chol) {
69  stan::math::check_square(function, "Cholesky factor", L_chol);
70  stan::math::check_lower_triangular(function,
71  "Cholesky factor", L_chol);
72  stan::math::check_size_match(function,
73  "Dimension of mean vector", dimension_,
74  "Dimension of Cholesky factor", L_chol.rows());
75  stan::math::check_not_nan(function, "Cholesky factor", L_chol);
76  }
77 
78 
79  public:
87  explicit normal_fullrank(size_t dimension)
88  : mu_(Eigen::VectorXd::Zero(dimension)),
89  L_chol_(Eigen::MatrixXd::Zero(dimension, dimension)),
90  dimension_(dimension) {
91  }
92 
93 
100  explicit normal_fullrank(const Eigen::VectorXd& cont_params)
101  : mu_(cont_params),
102  L_chol_(Eigen::MatrixXd::Identity(cont_params.size(),
103  cont_params.size())),
104  dimension_(cont_params.size()) {
105  }
106 
121  normal_fullrank(const Eigen::VectorXd& mu,
122  const Eigen::MatrixXd& L_chol)
123  : mu_(mu), L_chol_(L_chol), dimension_(mu.size()) {
124  static const char* function = "stan::variational::normal_fullrank";
125  validate_mean(function, mu);
126  validate_cholesky_factor(function, L_chol);
127  }
128 
132  int dimension() const { return dimension_; }
133 
137  const Eigen::VectorXd& mu() const { return mu_; }
138 
142  const Eigen::MatrixXd& L_chol() const { return L_chol_; }
143 
151  void set_mu(const Eigen::VectorXd& mu) {
152  static const char* function = "stan::variational::set_mu";
153  validate_mean(function, mu);
154  mu_ = mu;
155  }
156 
166  void set_L_chol(const Eigen::MatrixXd& L_chol) {
167  static const char* function = "stan::variational::set_L_chol";
168  validate_cholesky_factor(function, L_chol);
169  L_chol_ = L_chol;
170  }
171 
176  void set_to_zero() {
177  mu_ = Eigen::VectorXd::Zero(dimension_);
178  L_chol_ = Eigen::MatrixXd::Zero(dimension_, dimension_);
179  }
180 
188  return normal_fullrank(Eigen::VectorXd(mu_.array().square()),
189  Eigen::MatrixXd(L_chol_.array().square()));
190  }
191 
203  return normal_fullrank(Eigen::VectorXd(mu_.array().sqrt()),
204  Eigen::MatrixXd(L_chol_.array().sqrt()));
205  }
206 
219  static const char* function =
220  "stan::variational::normal_fullrank::operator=";
221  stan::math::check_size_match(function,
222  "Dimension of lhs", dimension_,
223  "Dimension of rhs", rhs.dimension());
224  mu_ = rhs.mu();
225  L_chol_ = rhs.L_chol();
226  return *this;
227  }
228 
241  static const char* function =
242  "stan::variational::normal_fullrank::operator+=";
243  stan::math::check_size_match(function,
244  "Dimension of lhs", dimension_,
245  "Dimension of rhs", rhs.dimension());
246  mu_ += rhs.mu();
247  L_chol_ += rhs.L_chol();
248  return *this;
249  }
250 
264  static const char* function =
265  "stan::variational::normal_fullrank::operator/=";
266 
267  stan::math::check_size_match(function,
268  "Dimension of lhs", dimension_,
269  "Dimension of rhs", rhs.dimension());
270 
271  mu_.array() /= rhs.mu().array();
272  L_chol_.array() /= rhs.L_chol().array();
273  return *this;
274  }
275 
287  normal_fullrank& operator+=(double scalar) {
288  mu_.array() += scalar;
289  L_chol_.array() += scalar;
290  return *this;
291  }
292 
305  normal_fullrank& operator*=(double scalar) {
306  mu_ *= scalar;
307  L_chol_ *= scalar;
308  return *this;
309  }
310 
318  const Eigen::VectorXd& mean() const {
319  return mu();
320  }
321 
331  double entropy() const {
332  static double mult = 0.5 * (1.0 + stan::math::LOG_TWO_PI);
333  double result = mult * dimension_;
334  for (int d = 0; d < dimension_; ++d) {
335  double tmp = fabs(L_chol_(d, d));
336  if (tmp != 0.0) result += log(tmp);
337  }
338  return result;
339  }
340 
353  Eigen::VectorXd transform(const Eigen::VectorXd& eta) const {
354  static const char* function =
355  "stan::variational::normal_fullrank::transform";
356  stan::math::check_size_match(function,
357  "Dimension of input vector", eta.size(),
358  "Dimension of mean vector", dimension_);
359  stan::math::check_not_nan(function, "Input vector", eta);
360 
361  return (L_chol_ * eta) + mu_;
362  }
363 
373  template <class BaseRNG>
374  void sample(BaseRNG& rng, Eigen::VectorXd& eta) const {
375  for (int d = 0; d < dimension_; ++d)
376  eta(d) = stan::math::normal_rng(0, 1, rng);
377  eta = transform(eta);
378  }
379 
397  template <class M, class BaseRNG>
398  void calc_grad(normal_fullrank& elbo_grad,
399  M& m,
400  Eigen::VectorXd& cont_params,
401  int n_monte_carlo_grad,
402  BaseRNG& rng,
404  const {
405  static const char* function =
406  "stan::variational::normal_fullrank::calc_grad";
407  stan::math::check_size_match(function,
408  "Dimension of elbo_grad", elbo_grad.dimension(),
409  "Dimension of variational q", dimension_);
410  stan::math::check_size_match(function,
411  "Dimension of variational q", dimension_,
412  "Dimension of variables in model", cont_params.size());
413 
414  Eigen::VectorXd mu_grad = Eigen::VectorXd::Zero(dimension_);
415  Eigen::MatrixXd L_grad = Eigen::MatrixXd::Zero(dimension_, dimension_);
416  double tmp_lp = 0.0;
417  Eigen::VectorXd tmp_mu_grad = Eigen::VectorXd::Zero(dimension_);
418  Eigen::VectorXd eta = Eigen::VectorXd::Zero(dimension_);
419  Eigen::VectorXd zeta = Eigen::VectorXd::Zero(dimension_);
420 
421  // Naive Monte Carlo integration
422  static const int n_retries = 10;
423  for (int i = 0, n_monte_carlo_drop = 0; i < n_monte_carlo_grad; ) {
424  // Draw from standard normal and transform to real-coordinate space
425  for (int d = 0; d < dimension_; ++d) {
426  eta(d) = stan::math::normal_rng(0, 1, rng);
427  }
428  zeta = transform(eta);
429  try {
430  std::stringstream ss;
431  stan::model::gradient(m, zeta, tmp_lp, tmp_mu_grad, &ss);
432  if (ss.str().length() > 0)
433  message_writer(ss.str());
434  stan::math::check_finite(function, "Gradient of mu", tmp_mu_grad);
435 
436  mu_grad += tmp_mu_grad;
437  for (int ii = 0; ii < dimension_; ++ii) {
438  for (int jj = 0; jj <= ii; ++jj) {
439  L_grad(ii, jj) += tmp_mu_grad(ii) * eta(jj);
440  }
441  }
442  ++i;
443  } catch (const std::exception& e) {
444  ++n_monte_carlo_drop;
445  if (n_monte_carlo_drop >= n_retries * n_monte_carlo_grad) {
446  const char* name = "The number of dropped evaluations";
447  const char* msg1 = "has reached its maximum amount (";
448  int y = n_retries * n_monte_carlo_grad;
449  const char* msg2 = "). Your model may be either severely "
450  "ill-conditioned or misspecified.";
451  stan::math::domain_error(function, name, y, msg1, msg2);
452  }
453  }
454  }
455  mu_grad /= static_cast<double>(n_monte_carlo_grad);
456  L_grad /= static_cast<double>(n_monte_carlo_grad);
457 
458  // Add gradient of entropy term
459  L_grad.diagonal().array() += L_chol_.diagonal().array().inverse();
460 
461  elbo_grad.set_mu(mu_grad);
462  elbo_grad.set_L_chol(L_grad);
463  }
464  };
465 
477  return lhs += rhs;
478  }
479 
490  return lhs /= rhs;
491  }
492 
503  return rhs += scalar;
504  }
505 
516  return rhs *= scalar;
517  }
518 
519  }
520 }
521 #endif
normal_fullrank(size_t dimension)
Construct a variational distribution of the specified dimensionality with a zero mean and Cholesky fa...
void calc_grad(normal_fullrank &elbo_grad, M &m, Eigen::VectorXd &cont_params, int n_monte_carlo_grad, BaseRNG &rng, interface_callbacks::writer::base_writer &message_writer) const
Calculates the "blackbox" gradient with respect to BOTH the location vector (mu) and the cholesky fac...
normal_fullrank & operator/=(const normal_fullrank &rhs)
Return this approximation after elementwise division by the specified approximation&#39;s mean and Choles...
normal_fullrank & operator+=(const normal_fullrank &rhs)
Add the mean and Cholesky factor of the covariance matrix of the specified approximation to this appr...
void sample(BaseRNG &rng, Eigen::VectorXd &eta) const
Set the specified vector to a draw from this variational approximation using the specified random num...
base_family operator/(base_family lhs, const base_family &rhs)
base_family operator*(double scalar, base_family rhs)
Probability, optimization and sampling library.
Eigen::VectorXd transform(const Eigen::VectorXd &eta) const
Return the transform of the sepcified vector using the Cholesky factor and mean vector.
Variational family approximation with full-rank multivariate normal distribution. ...
double entropy() const
Return the entropy of this approximation.
normal_fullrank(const Eigen::VectorXd &cont_params)
Construct a variational distribution with specified mean vector and Cholesky factor for identity cova...
void gradient(const M &model, const Eigen::Matrix< double, Eigen::Dynamic, 1 > &x, double &f, Eigen::Matrix< double, Eigen::Dynamic, 1 > &grad_f, std::ostream *msgs=0)
Definition: gradient.hpp:14
const Eigen::VectorXd & mean() const
Returns the mean vector for this approximation.
base_family operator+(base_family lhs, const base_family &rhs)
base_writer is an abstract base class defining the interface for Stan writer callbacks.
Definition: base_writer.hpp:20
normal_fullrank & operator=(const normal_fullrank &rhs)
Return this approximation after setting its mean vector and Cholesky factor for covariance to the val...
void set_mu(const Eigen::VectorXd &mu)
Set the mean vector to the specified value.
void set_L_chol(const Eigen::MatrixXd &L_chol)
Set the Cholesky factor to the specified value.
normal_fullrank square() const
Return a new full rank approximation resulting from squaring the entries in the mean and Cholesky fac...
const Eigen::MatrixXd & L_chol() const
Return the Cholesky factor of the covariance matrix.
normal_fullrank(const Eigen::VectorXd &mu, const Eigen::MatrixXd &L_chol)
Construct a variational distribution with specified mean and Cholesky factor for covariance.
void set_to_zero()
Set the mean vector and Cholesky factor for the covariance matrix to zero.
int dimension() const
Return the dimensionality of the approximation.
normal_fullrank & operator+=(double scalar)
Return this approximation after adding the specified scalar to each entry in the mean and cholesky fa...
const Eigen::VectorXd & mu() const
Return the mean vector.
normal_fullrank sqrt() const
Return a new full rank approximation resulting from taking the square root of the entries in the mean...
normal_fullrank & operator*=(double scalar)
Return this approximation after multiplying by the specified scalar to each entry in the mean and cho...

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