1 #ifndef STAN_VARIATIONAL_NORMAL_FULLRANK_HPP
2 #define STAN_VARIATIONAL_NORMAL_FULLRANK_HPP
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>
29 namespace variational {
46 Eigen::MatrixXd L_chol_;
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_);
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);
104 : mu_(Eigen::VectorXd::Zero(dimension)),
105 L_chol_(Eigen::MatrixXd::Zero(dimension, dimension)),
106 dimension_(dimension) {
118 L_chol_(Eigen::MatrixXd::Identity(cont_params.size(),
119 cont_params.size())),
120 dimension_(cont_params.size()) {
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);
153 const Eigen::VectorXd&
mu()
const {
return mu_; }
158 const Eigen::MatrixXd&
L_chol()
const {
return L_chol_; }
168 static const char*
function =
"stan::variational::set_mu";
169 validate_mean(
function, mu);
183 static const char*
function =
"stan::variational::set_L_chol";
184 validate_cholesky_factor(
function, L_chol);
193 mu_ = Eigen::VectorXd::Zero(dimension_);
194 L_chol_ = Eigen::MatrixXd::Zero(dimension_, dimension_);
205 Eigen::MatrixXd(L_chol_.array().square()));
220 Eigen::MatrixXd(L_chol_.array().sqrt()));
235 static const char*
function =
236 "stan::variational::normal_fullrank::operator=";
237 stan::math::check_size_match(
function,
238 "Dimension of lhs", dimension_,
257 static const char*
function =
258 "stan::variational::normal_fullrank::operator+=";
259 stan::math::check_size_match(
function,
260 "Dimension of lhs", dimension_,
280 static const char*
function =
281 "stan::variational::normal_fullrank::operator/=";
283 stan::math::check_size_match(
function,
284 "Dimension of lhs", dimension_,
287 mu_.array() /= rhs.
mu().array();
288 L_chol_.array() /= rhs.
L_chol().array();
304 mu_.array() += scalar;
305 L_chol_.array() += scalar;
334 const Eigen::VectorXd&
mean()
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);
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);
377 return (L_chol_ * eta) + mu_;
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);
413 template <
class M,
class BaseRNG>
416 Eigen::VectorXd& cont_params,
417 int n_monte_carlo_grad,
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());
430 Eigen::VectorXd mu_grad = Eigen::VectorXd::Zero(dimension_);
431 Eigen::MatrixXd L_grad = Eigen::MatrixXd::Zero(dimension_, dimension_);
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_);
438 static const int n_retries = 10;
439 for (
int i = 0, n_monte_carlo_drop = 0; i < n_monte_carlo_grad; ) {
441 for (
int d = 0; d < dimension_; ++d) {
442 eta(d) = stan::math::normal_rng(0, 1, rng);
446 std::stringstream ss;
448 if (ss.str().length() > 0)
449 message_writer(ss.str());
450 stan::math::check_finite(
function,
"Gradient of mu", tmp_mu_grad);
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);
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);
471 mu_grad /=
static_cast<double>(n_monte_carlo_grad);
472 L_grad /=
static_cast<double>(n_monte_carlo_grad);
475 L_grad.diagonal().array() += L_chol_.diagonal().array().inverse();
477 elbo_grad.
set_mu(mu_grad);
519 return rhs += scalar;
532 return rhs *= scalar;
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)
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.
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...