1 #ifndef STAN_VARIATIONAL_NORMAL_FULLRANK_HPP 2 #define STAN_VARIATIONAL_NORMAL_FULLRANK_HPP 4 #include <stan/math/prim/mat.hpp> 13 namespace variational {
30 Eigen::MatrixXd L_chol_;
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_);
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);
88 : mu_(Eigen::VectorXd::Zero(dimension)),
89 L_chol_(Eigen::MatrixXd::Zero(dimension, dimension)),
90 dimension_(dimension) {
102 L_chol_(Eigen::MatrixXd::Identity(cont_params.size(),
103 cont_params.size())),
104 dimension_(cont_params.size()) {
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);
137 const Eigen::VectorXd&
mu()
const {
return mu_; }
142 const Eigen::MatrixXd&
L_chol()
const {
return L_chol_; }
152 static const char*
function =
"stan::variational::set_mu";
153 validate_mean(
function, mu);
167 static const char*
function =
"stan::variational::set_L_chol";
168 validate_cholesky_factor(
function, L_chol);
177 mu_ = Eigen::VectorXd::Zero(dimension_);
178 L_chol_ = Eigen::MatrixXd::Zero(dimension_, dimension_);
189 Eigen::MatrixXd(L_chol_.array().square()));
204 Eigen::MatrixXd(L_chol_.array().sqrt()));
219 static const char*
function =
220 "stan::variational::normal_fullrank::operator=";
221 stan::math::check_size_match(
function,
222 "Dimension of lhs", dimension_,
241 static const char*
function =
242 "stan::variational::normal_fullrank::operator+=";
243 stan::math::check_size_match(
function,
244 "Dimension of lhs", dimension_,
264 static const char*
function =
265 "stan::variational::normal_fullrank::operator/=";
267 stan::math::check_size_match(
function,
268 "Dimension of lhs", dimension_,
271 mu_.array() /= rhs.
mu().array();
272 L_chol_.array() /= rhs.
L_chol().array();
288 mu_.array() += scalar;
289 L_chol_.array() += scalar;
318 const Eigen::VectorXd&
mean()
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);
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);
361 return (L_chol_ * eta) + mu_;
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);
397 template <
class M,
class BaseRNG>
400 Eigen::VectorXd& cont_params,
401 int n_monte_carlo_grad,
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());
414 Eigen::VectorXd mu_grad = Eigen::VectorXd::Zero(dimension_);
415 Eigen::MatrixXd L_grad = Eigen::MatrixXd::Zero(dimension_, dimension_);
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_);
422 static const int n_retries = 10;
423 for (
int i = 0, n_monte_carlo_drop = 0; i < n_monte_carlo_grad; ) {
425 for (
int d = 0; d < dimension_; ++d) {
426 eta(d) = stan::math::normal_rng(0, 1, rng);
430 std::stringstream ss;
432 if (ss.str().length() > 0)
433 message_writer(ss.str());
434 stan::math::check_finite(
function,
"Gradient of mu", tmp_mu_grad);
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);
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);
455 mu_grad /=
static_cast<double>(n_monte_carlo_grad);
456 L_grad /=
static_cast<double>(n_monte_carlo_grad);
459 L_grad.diagonal().array() += L_chol_.diagonal().array().inverse();
461 elbo_grad.
set_mu(mu_grad);
503 return rhs += scalar;
516 return rhs *= scalar;
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'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)
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.
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...