Stan  2.10.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/fun/Eigen.hpp>
5 #include <stan/math/prim/mat/fun/LDLT_factor.hpp>
6 #include <stan/math/prim/scal/fun/constants.hpp>
7 #include <stan/math/prim/mat/meta/get.hpp>
8 #include <stan/math/prim/arr/meta/get.hpp>
9 #include <stan/math/prim/mat/meta/length.hpp>
10 #include <stan/math/prim/mat/meta/is_vector.hpp>
11 #include <stan/math/prim/mat/meta/is_vector_like.hpp>
12 #include <stan/math/prim/mat/fun/value_of_rec.hpp>
13 #include <stan/math/prim/scal/prob/normal_rng.hpp>
14 #include <stan/math/prim/mat/err/check_square.hpp>
15 #include <stan/math/prim/mat/err/check_lower_triangular.hpp>
16 #include <stan/math/prim/scal/err/check_finite.hpp>
17 #include <stan/math/prim/scal/err/check_not_nan.hpp>
18 #include <stan/math/prim/scal/err/check_positive.hpp>
19 #include <stan/math/prim/scal/err/check_size_match.hpp>
20 #include <stan/math/prim/scal/err/domain_error.hpp>
21 #include <stan/model/util.hpp>
23 #include <algorithm>
24 #include <ostream>
25 #include <vector>
26 
27 namespace stan {
28 
29  namespace variational {
30 
35  class normal_fullrank : public base_family {
36  private:
40  Eigen::VectorXd mu_;
41 
46  Eigen::MatrixXd L_chol_;
47 
51  const int dimension_;
52 
61  void validate_mean(const char* function,
62  const Eigen::VectorXd& mu) {
63  stan::math::check_not_nan(function, "Mean vector", mu);
64  stan::math::check_size_match(function,
65  "Dimension of input vector", mu.size(),
66  "Dimension of current vector", dimension_);
67  }
68 
83  void validate_cholesky_factor(const char* function,
84  const Eigen::MatrixXd& L_chol) {
85  stan::math::check_square(function, "Cholesky factor", L_chol);
86  stan::math::check_lower_triangular(function,
87  "Cholesky factor", L_chol);
88  stan::math::check_size_match(function,
89  "Dimension of mean vector", dimension_,
90  "Dimension of Cholesky factor", L_chol.rows());
91  stan::math::check_not_nan(function, "Cholesky factor", L_chol);
92  }
93 
94 
95  public:
103  explicit normal_fullrank(size_t dimension)
104  : mu_(Eigen::VectorXd::Zero(dimension)),
105  L_chol_(Eigen::MatrixXd::Zero(dimension, dimension)),
106  dimension_(dimension) {
107  }
108 
109 
116  explicit normal_fullrank(const Eigen::VectorXd& cont_params)
117  : mu_(cont_params),
118  L_chol_(Eigen::MatrixXd::Identity(cont_params.size(),
119  cont_params.size())),
120  dimension_(cont_params.size()) {
121  }
122 
137  normal_fullrank(const Eigen::VectorXd& mu,
138  const Eigen::MatrixXd& L_chol)
139  : mu_(mu), L_chol_(L_chol), dimension_(mu.size()) {
140  static const char* function = "stan::variational::normal_fullrank";
141  validate_mean(function, mu);
142  validate_cholesky_factor(function, L_chol);
143  }
144 
148  int dimension() const { return dimension_; }
149 
153  const Eigen::VectorXd& mu() const { return mu_; }
154 
158  const Eigen::MatrixXd& L_chol() const { return L_chol_; }
159 
167  void set_mu(const Eigen::VectorXd& mu) {
168  static const char* function = "stan::variational::set_mu";
169  validate_mean(function, mu);
170  mu_ = mu;
171  }
172 
182  void set_L_chol(const Eigen::MatrixXd& L_chol) {
183  static const char* function = "stan::variational::set_L_chol";
184  validate_cholesky_factor(function, L_chol);
185  L_chol_ = L_chol;
186  }
187 
192  void set_to_zero() {
193  mu_ = Eigen::VectorXd::Zero(dimension_);
194  L_chol_ = Eigen::MatrixXd::Zero(dimension_, dimension_);
195  }
196 
204  return normal_fullrank(Eigen::VectorXd(mu_.array().square()),
205  Eigen::MatrixXd(L_chol_.array().square()));
206  }
207 
219  return normal_fullrank(Eigen::VectorXd(mu_.array().sqrt()),
220  Eigen::MatrixXd(L_chol_.array().sqrt()));
221  }
222 
235  static const char* function =
236  "stan::variational::normal_fullrank::operator=";
237  stan::math::check_size_match(function,
238  "Dimension of lhs", dimension_,
239  "Dimension of rhs", rhs.dimension());
240  mu_ = rhs.mu();
241  L_chol_ = rhs.L_chol();
242  return *this;
243  }
244 
257  static const char* function =
258  "stan::variational::normal_fullrank::operator+=";
259  stan::math::check_size_match(function,
260  "Dimension of lhs", dimension_,
261  "Dimension of rhs", rhs.dimension());
262  mu_ += rhs.mu();
263  L_chol_ += rhs.L_chol();
264  return *this;
265  }
266 
280  static const char* function =
281  "stan::variational::normal_fullrank::operator/=";
282 
283  stan::math::check_size_match(function,
284  "Dimension of lhs", dimension_,
285  "Dimension of rhs", rhs.dimension());
286 
287  mu_.array() /= rhs.mu().array();
288  L_chol_.array() /= rhs.L_chol().array();
289  return *this;
290  }
291 
303  normal_fullrank& operator+=(double scalar) {
304  mu_.array() += scalar;
305  L_chol_.array() += scalar;
306  return *this;
307  }
308 
321  normal_fullrank& operator*=(double scalar) {
322  mu_ *= scalar;
323  L_chol_ *= scalar;
324  return *this;
325  }
326 
334  const Eigen::VectorXd& mean() const {
335  return mu();
336  }
337 
347  double entropy() const {
348  static double mult = 0.5 * (1.0 + stan::math::LOG_TWO_PI);
349  double result = mult * dimension_;
350  for (int d = 0; d < dimension_; ++d) {
351  double tmp = fabs(L_chol_(d, d));
352  if (tmp != 0.0) result += log(tmp);
353  }
354  return result;
355  }
356 
369  Eigen::VectorXd transform(const Eigen::VectorXd& eta) const {
370  static const char* function =
371  "stan::variational::normal_fullrank::transform";
372  stan::math::check_size_match(function,
373  "Dimension of input vector", eta.size(),
374  "Dimension of mean vector", dimension_);
375  stan::math::check_not_nan(function, "Input vector", eta);
376 
377  return (L_chol_ * eta) + mu_;
378  }
379 
389  template <class BaseRNG>
390  void sample(BaseRNG& rng, Eigen::VectorXd& eta) const {
391  for (int d = 0; d < dimension_; ++d)
392  eta(d) = stan::math::normal_rng(0, 1, rng);
393  eta = transform(eta);
394  }
395 
413  template <class M, class BaseRNG>
414  void calc_grad(normal_fullrank& elbo_grad,
415  M& m,
416  Eigen::VectorXd& cont_params,
417  int n_monte_carlo_grad,
418  BaseRNG& rng,
420  const {
421  static const char* function =
422  "stan::variational::normal_fullrank::calc_grad";
423  stan::math::check_size_match(function,
424  "Dimension of elbo_grad", elbo_grad.dimension(),
425  "Dimension of variational q", dimension_);
426  stan::math::check_size_match(function,
427  "Dimension of variational q", dimension_,
428  "Dimension of variables in model", cont_params.size());
429 
430  Eigen::VectorXd mu_grad = Eigen::VectorXd::Zero(dimension_);
431  Eigen::MatrixXd L_grad = Eigen::MatrixXd::Zero(dimension_, dimension_);
432  double tmp_lp = 0.0;
433  Eigen::VectorXd tmp_mu_grad = Eigen::VectorXd::Zero(dimension_);
434  Eigen::VectorXd eta = Eigen::VectorXd::Zero(dimension_);
435  Eigen::VectorXd zeta = Eigen::VectorXd::Zero(dimension_);
436 
437  // Naive Monte Carlo integration
438  static const int n_retries = 10;
439  for (int i = 0, n_monte_carlo_drop = 0; i < n_monte_carlo_grad; ) {
440  // Draw from standard normal and transform to real-coordinate space
441  for (int d = 0; d < dimension_; ++d) {
442  eta(d) = stan::math::normal_rng(0, 1, rng);
443  }
444  zeta = transform(eta);
445  try {
446  std::stringstream ss;
447  stan::model::gradient(m, zeta, tmp_lp, tmp_mu_grad, &ss);
448  if (ss.str().length() > 0)
449  message_writer(ss.str());
450  stan::math::check_finite(function, "Gradient of mu", tmp_mu_grad);
451 
452  mu_grad += tmp_mu_grad;
453  for (int ii = 0; ii < dimension_; ++ii) {
454  for (int jj = 0; jj <= ii; ++jj) {
455  L_grad(ii, jj) += tmp_mu_grad(ii) * eta(jj);
456  }
457  }
458  ++i;
459  } catch (const std::exception& e) {
460  ++n_monte_carlo_drop;
461  if (n_monte_carlo_drop >= n_retries * n_monte_carlo_grad) {
462  const char* name = "The number of dropped evaluations";
463  const char* msg1 = "has reached its maximum amount (";
464  int y = n_retries * n_monte_carlo_grad;
465  const char* msg2 = "). Your model may be either severely "
466  "ill-conditioned or misspecified.";
467  stan::math::domain_error(function, name, y, msg1, msg2);
468  }
469  }
470  }
471  mu_grad /= static_cast<double>(n_monte_carlo_grad);
472  L_grad /= static_cast<double>(n_monte_carlo_grad);
473 
474  // Add gradient of entropy term
475  L_grad.diagonal().array() += L_chol_.diagonal().array().inverse();
476 
477  elbo_grad.set_mu(mu_grad);
478  elbo_grad.set_L_chol(L_grad);
479  }
480  };
481 
493  return lhs += rhs;
494  }
495 
506  return lhs /= rhs;
507  }
508 
519  return rhs += scalar;
520  }
521 
532  return rhs *= scalar;
533  }
534 
535  }
536 }
537 #endif
const Eigen::VectorXd & mean() const
Returns the mean vector for this approximation.
normal_fullrank(size_t dimension)
Construct a variational distribution of the specified dimensionality with a zero mean and Cholesky fa...
Eigen::VectorXd transform(const Eigen::VectorXd &eta) const
Return the transform of the sepcified vector using the Cholesky factor and mean vector.
normal_fullrank & operator/=(const normal_fullrank &rhs)
Return this approximation after elementwise division by the specified approximation'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...
base_family operator/(base_family lhs, const base_family &rhs)
base_family operator*(double scalar, base_family rhs)
Probability, optimization and sampling library.
double entropy() const
Return the entropy of this approximation.
normal_fullrank square() const
Return a new full rank approximation resulting from squaring the entries in the mean and Cholesky fac...
Variational family approximation with full-rank multivariate normal distribution. ...
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: util.hpp:422
const Eigen::VectorXd & mu() const
Return the mean vector.
int dimension() const
Return the dimensionality of the 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
void sample(BaseRNG &rng, Eigen::VectorXd &eta) const
Set the specified vector to a draw from this variational approximation using the specified random num...
normal_fullrank & operator=(const normal_fullrank &rhs)
Return this approximation after setting its mean vector and Cholesky factor for covariance to the val...
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...
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(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.
normal_fullrank & operator+=(double scalar)
Return this approximation after adding the specified scalar to each entry in the mean and cholesky fa...
const Eigen::MatrixXd & L_chol() const
Return the Cholesky factor of the covariance matrix.
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.