From 1952caf94101c363d7d7bcd530ffa59652d6d7f5 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Thu, 18 Apr 2019 11:22:36 -0700 Subject: [PATCH 01/20] Added in switching adaptation --- lib/stan_math | 2 +- .../mcmc/hmc/nuts/adapt_switching_e_nuts.hpp | 59 +++++ src/stan/mcmc/stepsize_switching_adapter.hpp | 48 +++++ src/stan/mcmc/switching_adaptation.hpp | 202 ++++++++++++++++++ .../sample/hmc_nuts_switching_e_adapt.hpp | 181 ++++++++++++++++ 5 files changed, 491 insertions(+), 1 deletion(-) create mode 100644 src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp create mode 100644 src/stan/mcmc/stepsize_switching_adapter.hpp create mode 100644 src/stan/mcmc/switching_adaptation.hpp create mode 100644 src/stan/services/sample/hmc_nuts_switching_e_adapt.hpp diff --git a/lib/stan_math b/lib/stan_math index 48d34c56fe..5fe4e2996c 160000 --- a/lib/stan_math +++ b/lib/stan_math @@ -1 +1 @@ -Subproject commit 48d34c56fe3750e30122ef35bc92350c2f4a0775 +Subproject commit 5fe4e2996cb93220ccff9b426ae837073dba3338 diff --git a/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp new file mode 100644 index 0000000000..c9f660b6e0 --- /dev/null +++ b/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp @@ -0,0 +1,59 @@ +#ifndef STAN_MCMC_HMC_NUTS_ADAPT_SWITCHING_E_NUTS_HPP +#define STAN_MCMC_HMC_NUTS_ADAPT_SWITCHING_E_NUTS_HPP + +#include +#include +#include + +namespace stan { + namespace mcmc { + /** + * The No-U-Turn sampler (NUTS) with multinomial sampling + * with a Gaussian-Euclidean disintegration and adaptive + * dense metric and adaptive step size + */ + template + class adapt_switching_e_nuts : public dense_e_nuts, + public stepsize_switching_adapter { + protected: + const Model& model_; + public: + adapt_switching_e_nuts(const Model& model, BaseRNG& rng) + : model_(model), dense_e_nuts(model, rng), + stepsize_switching_adapter(model.num_params_r()) {} + + ~adapt_switching_e_nuts() {} + + sample + transition(sample& init_sample, callbacks::logger& logger) { + sample s = dense_e_nuts::transition(init_sample, + logger); + + if (this->adapt_flag_) { + this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_, + s.accept_stat()); + + bool update = this->switching_adaptation_.learn_covariance( + model_, + this->z_.inv_e_metric_, + this->z_.q); + + if (update) { + this->init_stepsize(logger); + + this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_)); + this->stepsize_adaptation_.restart(); + } + } + return s; + } + + void disengage_adaptation() { + base_adapter::disengage_adaptation(); + this->stepsize_adaptation_.complete_adaptation(this->nom_epsilon_); + } + }; + + } // mcmc +} // stan +#endif diff --git a/src/stan/mcmc/stepsize_switching_adapter.hpp b/src/stan/mcmc/stepsize_switching_adapter.hpp new file mode 100644 index 0000000000..aaa35ad623 --- /dev/null +++ b/src/stan/mcmc/stepsize_switching_adapter.hpp @@ -0,0 +1,48 @@ +#ifndef STAN_MCMC_STEPSIZE_SWITCHING_ADAPTER_HPP +#define STAN_MCMC_STEPSIZE_SWITCHING_ADAPTER_HPP + +#include +#include +#include +#include + +namespace stan { + + namespace mcmc { + + class stepsize_switching_adapter: public base_adapter { + public: + explicit stepsize_switching_adapter(int n) + : switching_adaptation_(n) { + } + + stepsize_adaptation& get_stepsize_adaptation() { + return stepsize_adaptation_; + } + + switching_adaptation& get_switching_adaptation() { + return switching_adaptation_; + } + + void set_window_params(unsigned int num_warmup, + unsigned int init_buffer, + unsigned int term_buffer, + unsigned int base_window, + callbacks::logger& logger) { + switching_adaptation_.set_window_params(num_warmup, + init_buffer, + term_buffer, + base_window, + logger); + } + + protected: + stepsize_adaptation stepsize_adaptation_; + switching_adaptation switching_adaptation_; + }; + + } // mcmc + +} // stan + +#endif diff --git a/src/stan/mcmc/switching_adaptation.hpp b/src/stan/mcmc/switching_adaptation.hpp new file mode 100644 index 0000000000..7630cc7a62 --- /dev/null +++ b/src/stan/mcmc/switching_adaptation.hpp @@ -0,0 +1,202 @@ +#ifndef STAN_MCMC_SWITCHING_ADAPTATION_HPP +#define STAN_MCMC_SWITCHING_ADAPTATION_HPP + +#include +#include +#include + +namespace stan { + + namespace mcmc { + template + struct log_prob_wrapper_covar { + const Model& model_; + log_prob_wrapper_covar(const Model& model) : model_(model) {} + + template + T operator()(const Eigen::Matrix& q) const { + return model_.template log_prob(const_cast& >(q), &std::cout); + } + }; + + template + class scaled_hessian_vector { + private: + const Model& model_; + const Eigen::MatrixXd& L_; + const Eigen::VectorXd& q_; + public: + scaled_hessian_vector(const Model& model, + const Eigen::MatrixXd& L, + const Eigen::VectorXd& q) : model_(model), + L_(L), + q_(q) {} + + int rows() { return q_.size(); } + int cols() { return q_.size(); } + + void perform_op(const double* x_in, double* y_out) { + Eigen::Map x(x_in, cols()); + Eigen::Map y(y_out, rows()); + + double lp; + Eigen::VectorXd grad1; + Eigen::VectorXd grad2; + //stan::math::hessian_times_vector(log_prob_wrapper_covar(model), q, x, lp, Ax); + double dx = 1e-5; + Eigen::VectorXd dr = L_ * x * dx; + stan::math::gradient(log_prob_wrapper_covar(model_), q_ + dr / 2.0, lp, grad1); + stan::math::gradient(log_prob_wrapper_covar(model_), q_ - dr / 2.0, lp, grad2); + y = L_.transpose() * (grad1 - grad2) / dx; + } + }; + + class switching_adaptation: public windowed_adaptation { + public: + explicit switching_adaptation(int n) + : windowed_adaptation("covariance") {} + + Eigen::MatrixXd covariance(const Eigen::MatrixXd& Y) { + Eigen::MatrixXd centered = Y.colwise() - Y.rowwise().mean(); + return centered * centered.transpose() / std::max(centered.rows() - 1.0, 1.0); + } + + template + double top_eigenvalue(const Model& model, const Eigen::MatrixXd& L, const Eigen::VectorXd& q) { + Eigen::VectorXd eigenvalues; + Eigen::MatrixXd eigenvectors; + + scaled_hessian_vector op(model, L, q); + + Spectra::SymEigsSolver eigs(&op, 1, 2); + eigs.init(); + eigs.compute(); + + if(eigs.info() != Spectra::SUCCESSFUL) { + throw std::domain_error("Failed to compute eigenvalue of Hessian of log density. The switching metric requires these"); + } + + return eigs.eigenvalues()(0); + } + + double bottom_eigenvalue_estimate(const Eigen::MatrixXd& L, const Eigen::MatrixXd& covar) { + Eigen::MatrixXd S = L.template triangularView(). + solve(L.template triangularView().solve(covar).transpose()).transpose(); + + Spectra::DenseSymMatProd op(S); + Spectra::SymEigsSolver eigs(&op, 1, 2); + eigs.init(); + eigs.compute(); + + if(eigs.info() != Spectra::SUCCESSFUL) { + throw std::domain_error("Failed to compute eigenvalue of covariance of log density. The switching metric requires these"); + } + + return -1.0 / eigs.eigenvalues()(0); + } + + template + bool learn_covariance(const Model& model, Eigen::MatrixXd& covar, const Eigen::VectorXd& q) { + if (adaptation_window()) + qs_.push_back(q); + + if (end_adaptation_window()) { + compute_next_window(); + + int N = q.size(); + int M = qs_.size(); + + Eigen::MatrixXd Y = Eigen::MatrixXd::Zero(N, M); + std::vector idxs(M); + for(int i = 0; i < qs_.size(); i++) + idxs[i] = i; + + std::random_shuffle(idxs.begin(), idxs.end()); + for(int i = 0; i < qs_.size(); i++) + Y.block(0, i, N, 1) = qs_[idxs[i]]; + + bool use_dense = false; + for(auto state : { "selection", "refinement" }) { + Eigen::MatrixXd Ytrain; + Eigen::MatrixXd Ytest; + + if(state == "selection") { + int Ntest; + Ntest = int(0.2 * Y.cols()); + if(Ntest < 5) { + Ntest = 5; + } + + if(Y.cols() < 10) { + throw std::runtime_error("Each warmup stage must have at least 10 samples"); + } + + std::cout << "train: " << Y.cols() - Ntest << ", test: " << Ntest << std::endl; + Ytrain = Y.block(0, 0, N, Y.cols() - Ntest); + Ytest = Y.block(0, Ytrain.cols(), N, Ntest); + } else { + Ytrain = Y; + } + + Eigen::MatrixXd cov_train = covariance(Ytrain); + Eigen::MatrixXd cov_test = covariance(Ytest); + + Eigen::MatrixXd dense = (N / (N + 5.0)) * cov_train + + 1e-3 * (5.0 / (N + 5.0)) * Eigen::MatrixXd::Identity(cov_train.rows(), cov_train.cols()); + Eigen::MatrixXd diag = dense.diagonal().asDiagonal(); + + covar = dense; + + if(state == "selection") { + Eigen::MatrixXd L_dense = dense.llt().matrixL(); + Eigen::MatrixXd L_diag = diag.diagonal().array().sqrt().matrix().asDiagonal(); + + double low_eigenvalue_dense = bottom_eigenvalue_estimate(L_dense, cov_test); + double low_eigenvalue_diag = bottom_eigenvalue_estimate(L_diag, cov_test); + + double c_dense = 0.0; + double c_diag = 0.0; + for(int i = 0; i < 5; i++) { + double high_eigenvalue_dense = top_eigenvalue(model, L_dense, Ytest.block(0, i, N, 1)); + double high_eigenvalue_diag = top_eigenvalue(model, L_diag, Ytest.block(0, i, N, 1)); + + c_dense = std::max(c_dense, std::sqrt(high_eigenvalue_dense / low_eigenvalue_dense)); + c_diag = std::max(c_diag, std::sqrt(high_eigenvalue_diag / low_eigenvalue_diag)); + } + + std::cout << "adapt: " << adapt_window_counter_ << ", which: dense, max: " << c_dense << std::endl; + std::cout << "adapt: " << adapt_window_counter_ << ", which: diag, max: " << c_diag << std::endl; + + if(c_dense < c_diag) { + use_dense = true; + } else { + use_dense = false; + } + } else { + if(use_dense) { + covar = dense; + } else { + covar = diag; + } + } + } + + ++adapt_window_counter_; + qs_.clear(); + + return true; + } + + ++adapt_window_counter_; + return false; + } + + protected: + std::vector< Eigen::VectorXd > qs_; + }; + + } // mcmc + +} // stan + +#endif diff --git a/src/stan/services/sample/hmc_nuts_switching_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_switching_e_adapt.hpp new file mode 100644 index 0000000000..dea87cc6c3 --- /dev/null +++ b/src/stan/services/sample/hmc_nuts_switching_e_adapt.hpp @@ -0,0 +1,181 @@ +#ifndef STAN_SERVICES_SAMPLE_HMC_NUTS_SWITCHING_E_ADAPT_HPP +#define STAN_SERVICES_SAMPLE_HMC_NUTS_SWITCHING_E_ADAPT_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { + namespace services { + namespace sample { + + /** + * Runs HMC with NUTS with adaptation using dense Euclidean metric + * with a pre-specified Euclidean metric. + * + * @tparam Model Model class + * @param[in] model Input model to test (with data already instantiated) + * @param[in] init var context for initialization + * @param[in] init_inv_metric var context exposing an initial dense + inverse Euclidean metric (must be positive definite) + * @param[in] random_seed random seed for the random number generator + * @param[in] chain chain id to advance the pseudo random number generator + * @param[in] init_radius radius to initialize + * @param[in] num_warmup Number of warmup samples + * @param[in] num_samples Number of samples + * @param[in] num_thin Number to thin the samples + * @param[in] save_warmup Indicates whether to save the warmup iterations + * @param[in] refresh Controls the output + * @param[in] stepsize initial stepsize for discrete evolution + * @param[in] stepsize_jitter uniform random jitter of stepsize + * @param[in] max_depth Maximum tree depth + * @param[in] delta adaptation target acceptance statistic + * @param[in] gamma adaptation regularization scale + * @param[in] kappa adaptation relaxation exponent + * @param[in] t0 adaptation iteration offset + * @param[in] init_buffer width of initial fast adaptation interval + * @param[in] term_buffer width of final fast adaptation interval + * @param[in] window initial width of slow adaptation interval + * @param[in,out] interrupt Callback for interrupts + * @param[in,out] logger Logger for messages + * @param[in,out] init_writer Writer callback for unconstrained inits + * @param[in,out] sample_writer Writer for draws + * @param[in,out] diagnostic_writer Writer for diagnostic information + * @return error_codes::OK if successful + */ + template + int hmc_nuts_switching_e_adapt(Model& model, stan::io::var_context& init, + stan::io::var_context& init_inv_metric, + unsigned int random_seed, unsigned int chain, + double init_radius, int num_warmup, + int num_samples, int num_thin, + bool save_warmup, int refresh, double stepsize, + double stepsize_jitter, int max_depth, + double delta, double gamma, double kappa, + double t0, unsigned int init_buffer, + unsigned int term_buffer, unsigned int window, + callbacks::interrupt& interrupt, + callbacks::logger& logger, + callbacks::writer& init_writer, + callbacks::writer& sample_writer, + callbacks::writer& diagnostic_writer) { + boost::ecuyer1988 rng = util::create_rng(random_seed, chain); + + std::vector disc_vector; + std::vector cont_vector + = util::initialize(model, init, rng, init_radius, true, + logger, init_writer); + + Eigen::MatrixXd inv_metric; + try { + inv_metric = + util::read_dense_inv_metric(init_inv_metric, model.num_params_r(), + logger); + util::validate_dense_inv_metric(inv_metric, logger); + } catch (const std::domain_error& e) { + return error_codes::CONFIG; + } + + stan::mcmc::adapt_switching_e_nuts + sampler(model, rng); + + sampler.set_metric(inv_metric); + + sampler.set_nominal_stepsize(stepsize); + sampler.set_stepsize_jitter(stepsize_jitter); + sampler.set_max_depth(max_depth); + + sampler.get_stepsize_adaptation().set_mu(log(10 * stepsize)); + sampler.get_stepsize_adaptation().set_delta(delta); + sampler.get_stepsize_adaptation().set_gamma(gamma); + sampler.get_stepsize_adaptation().set_kappa(kappa); + sampler.get_stepsize_adaptation().set_t0(t0); + + sampler.set_window_params(num_warmup, init_buffer, term_buffer, + window, logger); + + util::run_adaptive_sampler(sampler, model, cont_vector, num_warmup, + num_samples, num_thin, refresh, save_warmup, + rng, interrupt, logger, + sample_writer, diagnostic_writer); + + return error_codes::OK; + } + + /** + * Runs HMC with NUTS with adaptation using dense Euclidean metric, + * with identity matrix as initial inv_metric. + * + * @tparam Model Model class + * @param[in] model Input model to test (with data already instantiated) + * @param[in] init var context for initialization + * @param[in] random_seed random seed for the random number generator + * @param[in] chain chain id to advance the pseudo random number generator + * @param[in] init_radius radius to initialize + * @param[in] num_warmup Number of warmup samples + * @param[in] num_samples Number of samples + * @param[in] num_thin Number to thin the samples + * @param[in] save_warmup Indicates whether to save the warmup iterations + * @param[in] refresh Controls the output + * @param[in] stepsize initial stepsize for discrete evolution + * @param[in] stepsize_jitter uniform random jitter of stepsize + * @param[in] max_depth Maximum tree depth + * @param[in] delta adaptation target acceptance statistic + * @param[in] gamma adaptation regularization scale + * @param[in] kappa adaptation relaxation exponent + * @param[in] t0 adaptation iteration offset + * @param[in] init_buffer width of initial fast adaptation interval + * @param[in] term_buffer width of final fast adaptation interval + * @param[in] window initial width of slow adaptation interval + * @param[in,out] interrupt Callback for interrupts + * @param[in,out] logger Logger for messages + * @param[in,out] init_writer Writer callback for unconstrained inits + * @param[in,out] sample_writer Writer for draws + * @param[in,out] diagnostic_writer Writer for diagnostic information + * @return error_codes::OK if successful + */ + template + int hmc_nuts_switching_e_adapt(Model& model, stan::io::var_context& init, + unsigned int random_seed, unsigned int chain, + double init_radius, int num_warmup, + int num_samples, int num_thin, + bool save_warmup, int refresh, double stepsize, + double stepsize_jitter, int max_depth, + double delta, double gamma, double kappa, + double t0, unsigned int init_buffer, + unsigned int term_buffer, unsigned int window, + callbacks::interrupt& interrupt, + callbacks::logger& logger, + callbacks::writer& init_writer, + callbacks::writer& sample_writer, + callbacks::writer& diagnostic_writer) { + stan::io::dump dmp = + util::create_unit_e_dense_inv_metric(model.num_params_r()); + stan::io::var_context& unit_e_metric = dmp; + + return hmc_nuts_switching_e_adapt(model, init, unit_e_metric, + random_seed, chain, init_radius, + num_warmup, num_samples, num_thin, + save_warmup, refresh, + stepsize, stepsize_jitter, max_depth, + delta, gamma, kappa, t0, + init_buffer, term_buffer, window, + interrupt, logger, + init_writer, sample_writer, + diagnostic_writer); + } + + } + } +} +#endif From 4dc0fff4991d9759f9e0fb98d16f402d4f9606a8 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sat, 20 Apr 2019 14:27:10 -0700 Subject: [PATCH 02/20] Simplified code and added some comments --- lib/stan_math | 2 +- .../hmc/hamiltonians/switching_e_metric.hpp | 76 ++++++++ .../hmc/hamiltonians/switching_e_point.hpp | 81 +++++++++ .../mcmc/hmc/nuts/adapt_switching_e_nuts.hpp | 11 +- src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp | 26 +++ src/stan/mcmc/switching_adaptation.hpp | 164 +++++++++++------- 6 files changed, 291 insertions(+), 69 deletions(-) create mode 100644 src/stan/mcmc/hmc/hamiltonians/switching_e_metric.hpp create mode 100644 src/stan/mcmc/hmc/hamiltonians/switching_e_point.hpp create mode 100644 src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp diff --git a/lib/stan_math b/lib/stan_math index 5fe4e2996c..c2fb25c3f5 160000 --- a/lib/stan_math +++ b/lib/stan_math @@ -1 +1 @@ -Subproject commit 5fe4e2996cb93220ccff9b426ae837073dba3338 +Subproject commit c2fb25c3f53d3d33157a27ce995911b0bc92d55d diff --git a/src/stan/mcmc/hmc/hamiltonians/switching_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/switching_e_metric.hpp new file mode 100644 index 0000000000..b86cf3a913 --- /dev/null +++ b/src/stan/mcmc/hmc/hamiltonians/switching_e_metric.hpp @@ -0,0 +1,76 @@ +#ifndef STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_METRIC_HPP +#define STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_METRIC_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace stan { + namespace mcmc { + + // Euclidean manifold with dense metric + template + class switching_e_metric + : public base_hamiltonian { + public: + explicit switching_e_metric(const Model& model) + : base_hamiltonian(model) {} + + double T(switching_e_point& z) { + return 0.5 * z.p.transpose() * z.inv_e_metric_ * z.p; + } + + double tau(switching_e_point& z) { + return T(z); + } + + double phi(switching_e_point& z) { + return this->V(z); + } + + double dG_dt(switching_e_point& z, callbacks::logger& logger) { + return 2 * T(z) - z.q.dot(z.g); + } + + Eigen::VectorXd dtau_dq(switching_e_point& z, callbacks::logger& logger) { + return Eigen::VectorXd::Zero(this->model_.num_params_r()); + } + + Eigen::VectorXd dtau_dp(switching_e_point& z) { + if(z.is_diagonal_) { + return z.inv_e_metric_.diagonal().cwiseProduct(z.p); + } else { + return z.inv_e_metric_ * z.p; + } + } + + Eigen::VectorXd dphi_dq(switching_e_point& z, callbacks::logger& logger) { + return z.g; + } + + void sample_p(switching_e_point& z, BaseRNG& rng) { + typedef typename stan::math::index_type::type idx_t; + boost::variate_generator > + rand_gaus(rng, boost::normal_distribution<>()); + + if(z.is_diagonal_) { + for (int i = 0; i < z.p.size(); ++i) + z.p(i) = rand_gaus() / sqrt(z.inv_e_metric_(i, i)); + } else { + Eigen::VectorXd u(z.p.size()); + + for (idx_t i = 0; i < u.size(); ++i) + u(i) = rand_gaus(); + + z.p = z.inv_e_metric_.llt().matrixU().solve(u); + } + } + }; + + } // mcmc +} // stan +#endif diff --git a/src/stan/mcmc/hmc/hamiltonians/switching_e_point.hpp b/src/stan/mcmc/hmc/hamiltonians/switching_e_point.hpp new file mode 100644 index 0000000000..3d72230fab --- /dev/null +++ b/src/stan/mcmc/hmc/hamiltonians/switching_e_point.hpp @@ -0,0 +1,81 @@ +#ifndef STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_POINT_HPP +#define STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_POINT_HPP + +#include +#include + +namespace stan { + namespace mcmc { + /** + * Point in a phase space with a base + * Euclidean manifold with switching metric + */ + class switching_e_point: public ps_point { + public: + /** + * Inverse mass matrix. + */ + Eigen::MatrixXd inv_e_metric_; + + /** + * Is inv_e_metric_ diagonal or not + */ + bool is_diagonal_; + + /** + * Construct a switching point in n-dimensional phase space + * with identity matrix as inverse mass matrix. + * + * @param n number of dimensions + */ + explicit switching_e_point(int n) + : ps_point(n), inv_e_metric_(n, n), is_diagonal_(true) { + inv_e_metric_.setIdentity(); + } + + /** + * Copy constructor which does fast copy of inverse mass matrix. + * + * @param z point to copy + */ + switching_e_point(const switching_e_point& z) + : ps_point(z), inv_e_metric_(z.inv_e_metric_.rows(), + z.inv_e_metric_.cols()) { + fast_matrix_copy_(inv_e_metric_, z.inv_e_metric_); + is_diagonal_ = z.is_diagonal_; + } + + /** + * Set elements of mass matrix + * + * @param inv_e_metric initial mass matrix + */ + void + set_metric(const Eigen::MatrixXd& inv_e_metric) { + inv_e_metric_ = inv_e_metric; + is_diagonal_ = false; + } + + /** + * Write elements of mass matrix to string and handoff to writer. + * + * @param writer Stan writer callback + */ + inline + void + write_metric(stan::callbacks::writer& writer) { + writer("Elements of inverse mass matrix:"); + for (int i = 0; i < inv_e_metric_.rows(); ++i) { + std::stringstream inv_e_metric_ss; + inv_e_metric_ss << inv_e_metric_(i, 0); + for (int j = 1; j < inv_e_metric_.cols(); ++j) + inv_e_metric_ss << ", " << inv_e_metric_(i, j); + writer(inv_e_metric_ss.str()); + } + } + }; + + } // mcmc +} // stan + +#endif diff --git a/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp index c9f660b6e0..b0a21c9e02 100644 --- a/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp +++ b/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp @@ -3,7 +3,7 @@ #include #include -#include +#include namespace stan { namespace mcmc { @@ -13,21 +13,21 @@ namespace stan { * dense metric and adaptive step size */ template - class adapt_switching_e_nuts : public dense_e_nuts, + class adapt_switching_e_nuts : public switching_e_nuts, public stepsize_switching_adapter { protected: const Model& model_; public: adapt_switching_e_nuts(const Model& model, BaseRNG& rng) - : model_(model), dense_e_nuts(model, rng), + : model_(model), switching_e_nuts(model, rng), stepsize_switching_adapter(model.num_params_r()) {} ~adapt_switching_e_nuts() {} sample transition(sample& init_sample, callbacks::logger& logger) { - sample s = dense_e_nuts::transition(init_sample, - logger); + sample s = switching_e_nuts::transition(init_sample, + logger); if (this->adapt_flag_) { this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_, @@ -36,6 +36,7 @@ namespace stan { bool update = this->switching_adaptation_.learn_covariance( model_, this->z_.inv_e_metric_, + this->z_.is_diagonal_, this->z_.q); if (update) { diff --git a/src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp new file mode 100644 index 0000000000..cd91df6fb8 --- /dev/null +++ b/src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp @@ -0,0 +1,26 @@ +#ifndef STAN_MCMC_HMC_NUTS_SWITCHING_E_NUTS_HPP +#define STAN_MCMC_HMC_NUTS_SWITCHING_E_NUTS_HPP + +#include +#include +#include +#include + +namespace stan { + namespace mcmc { + /** + * The No-U-Turn sampler (NUTS) with multinomial sampling + * with a Gaussian-Euclidean disintegration and dense metric + */ + template + class switching_e_nuts : public base_nuts { + public: + switching_e_nuts(const Model& model, BaseRNG& rng) + : base_nuts(model, rng) { } + }; + + } // mcmc +} // stan +#endif diff --git a/src/stan/mcmc/switching_adaptation.hpp b/src/stan/mcmc/switching_adaptation.hpp index 7630cc7a62..088543c27b 100644 --- a/src/stan/mcmc/switching_adaptation.hpp +++ b/src/stan/mcmc/switching_adaptation.hpp @@ -19,84 +19,121 @@ namespace stan { } }; - template - class scaled_hessian_vector { - private: - const Model& model_; - const Eigen::MatrixXd& L_; - const Eigen::VectorXd& q_; - public: - scaled_hessian_vector(const Model& model, - const Eigen::MatrixXd& L, - const Eigen::VectorXd& q) : model_(model), - L_(L), - q_(q) {} - - int rows() { return q_.size(); } - int cols() { return q_.size(); } - - void perform_op(const double* x_in, double* y_out) { - Eigen::Map x(x_in, cols()); - Eigen::Map y(y_out, rows()); - - double lp; - Eigen::VectorXd grad1; - Eigen::VectorXd grad2; - //stan::math::hessian_times_vector(log_prob_wrapper_covar(model), q, x, lp, Ax); - double dx = 1e-5; - Eigen::VectorXd dr = L_ * x * dx; - stan::math::gradient(log_prob_wrapper_covar(model_), q_ + dr / 2.0, lp, grad1); - stan::math::gradient(log_prob_wrapper_covar(model_), q_ - dr / 2.0, lp, grad2); - y = L_.transpose() * (grad1 - grad2) / dx; - } - }; - class switching_adaptation: public windowed_adaptation { public: explicit switching_adaptation(int n) : windowed_adaptation("covariance") {} + /** + * Compute the covariance of data in Y. Rows of Y are different data. Columns of Y are different variables. + * + * @param Y Data + * @return Covariance of Y + */ Eigen::MatrixXd covariance(const Eigen::MatrixXd& Y) { Eigen::MatrixXd centered = Y.colwise() - Y.rowwise().mean(); return centered * centered.transpose() / std::max(centered.rows() - 1.0, 1.0); } - - template - double top_eigenvalue(const Model& model, const Eigen::MatrixXd& L, const Eigen::VectorXd& q) { - Eigen::VectorXd eigenvalues; - Eigen::MatrixXd eigenvectors; - - scaled_hessian_vector op(model, L, q); - Spectra::SymEigsSolver eigs(&op, 1, 2); - eigs.init(); - eigs.compute(); - - if(eigs.info() != Spectra::SUCCESSFUL) { - throw std::domain_error("Failed to compute eigenvalue of Hessian of log density. The switching metric requires these"); + /** + * Compute the largest magnitude eigenvalue of a symmetric matrix using the power method. The function f + * should return the product of that matrix with an abitrary vector. + * + * f should take one Eigen::VectorXd argument, x, and return the product of a matrix with x as + * an Eigen::VectorXd argument of the same size. + * + * The eigenvalue is estimated iteratively. If the kth estimate is e_k, then the function returns when + * either abs(e_{k + 1} - e_k) < tol * abs(e_k) or the maximum number of iterations have been performed + * + * This means the returned eigenvalue might not be computed to full precision + * + * @param initial_guess Initial guess of the eigenvector of the largest eigenvalue + * @param max_iterations Maximum number of power iterations + * @param tol Relative tolerance + * @return Largest magnitude eigenvalue of operator f + */ + template + double power_method(F& f, const Eigen::VectorXd& initial_guess, int max_iterations, double tol) { + Eigen::VectorXd v = initial_guess; + double eval = 0.0; + + for(int i = 0; i < max_iterations; i++) { + Eigen::VectorXd Av = f(v); + double v_norm = v.norm(); + double new_eval = v.dot(Av) / (v_norm * v_norm); + if(std::abs(new_eval - eval) <= tol * std::abs(eval)) { + std::cout << "Converged at i = " << i << std::endl; + eval = new_eval; + break; + } + eval = new_eval; + v = Av / Av.norm(); } - return eigs.eigenvalues()(0); + return eval; } - double bottom_eigenvalue_estimate(const Eigen::MatrixXd& L, const Eigen::MatrixXd& covar) { - Eigen::MatrixXd S = L.template triangularView(). - solve(L.template triangularView().solve(covar).transpose()).transpose(); + /** + * Compute the largest eigenvalue of the Hessian of the log density rescaled by a metric, + * that is, the largest eigenvalue of L^T \nabla^2_{qq} H(q) L + * + * @tparam Model Type of model + * @param model Defines the log density + * @param q Point around which to compute the Hessian + * @param L Cholesky decomposition of Metric + * @return Largest eigenvalue + */ + template + double eigenvalue_scaled_hessian(const Model& model, const Eigen::MatrixXd& L, const Eigen::VectorXd& q) { + Eigen::VectorXd eigenvalues; + Eigen::MatrixXd eigenvectors; - Spectra::DenseSymMatProd op(S); - Spectra::SymEigsSolver eigs(&op, 1, 2); - eigs.init(); - eigs.compute(); + auto hessian_vector = [&](const Eigen::VectorXd& x) -> Eigen::VectorXd { + double lp; + Eigen::VectorXd grad1; + Eigen::VectorXd grad2; + //stan::math::hessian_times_vector(log_prob_wrapper_covar(model), q, x, lp, Ax); + double dx = 1e-5; + Eigen::VectorXd dr = L * x * dx; + stan::math::gradient(log_prob_wrapper_covar(model), q + dr / 2.0, lp, grad1); + stan::math::gradient(log_prob_wrapper_covar(model), q - dr / 2.0, lp, grad2); + return L.transpose() * (grad1 - grad2) / dx; + }; + + return power_method(hessian_vector, Eigen::VectorXd::Random(q.size()), 100, 1e-3); + } - if(eigs.info() != Spectra::SUCCESSFUL) { - throw std::domain_error("Failed to compute eigenvalue of covariance of log density. The switching metric requires these"); - } + /** + * Compute the largest eigenvalue of the sample covariance rescaled by a metric, + * that is, the largest eigenvalue of L^{-T} \Sigma L^{-1} + * + * @param L Cholesky decomposition of Metric + * @param Sigma Sample covariance + * @return Largest eigenvalue + */ + double eigenvalue_scaled_covariance(const Eigen::MatrixXd& L, const Eigen::MatrixXd& Sigma) { + Eigen::MatrixXd S = L.template triangularView(). + solve(L.template triangularView().solve(Sigma).transpose()).transpose(); - return -1.0 / eigs.eigenvalues()(0); + auto Sx = [&](Eigen::VectorXd x) -> Eigen::VectorXd { + return S * x; + }; + + return power_method(Sx, Eigen::VectorXd::Random(Sigma.cols()), 100, 1e-3); } + /** + * Update the metric if at the end of an adaptation window. + * + * @tparam Model Type of model + * @param model Defines the log density + * @param covar[out] New metric + * @param covar_is_diagonal[out] Set to true if metric is diagonal, false otherwise + * @param q New MCMC draw + * @return True if this was the end of an adaptation window, false otherwise + */ template - bool learn_covariance(const Model& model, Eigen::MatrixXd& covar, const Eigen::VectorXd& q) { + bool learn_covariance(const Model& model, Eigen::MatrixXd& covar, bool& covar_is_diagonal, const Eigen::VectorXd& q) { if (adaptation_window()) qs_.push_back(q); @@ -131,7 +168,6 @@ namespace stan { throw std::runtime_error("Each warmup stage must have at least 10 samples"); } - std::cout << "train: " << Y.cols() - Ntest << ", test: " << Ntest << std::endl; Ytrain = Y.block(0, 0, N, Y.cols() - Ntest); Ytest = Y.block(0, Ytrain.cols(), N, Ntest); } else { @@ -151,14 +187,14 @@ namespace stan { Eigen::MatrixXd L_dense = dense.llt().matrixL(); Eigen::MatrixXd L_diag = diag.diagonal().array().sqrt().matrix().asDiagonal(); - double low_eigenvalue_dense = bottom_eigenvalue_estimate(L_dense, cov_test); - double low_eigenvalue_diag = bottom_eigenvalue_estimate(L_diag, cov_test); + double low_eigenvalue_dense = -1.0 / eigenvalue_scaled_covariance(L_dense, cov_test); + double low_eigenvalue_diag = -1.0 / eigenvalue_scaled_covariance(L_diag, cov_test); double c_dense = 0.0; double c_diag = 0.0; for(int i = 0; i < 5; i++) { - double high_eigenvalue_dense = top_eigenvalue(model, L_dense, Ytest.block(0, i, N, 1)); - double high_eigenvalue_diag = top_eigenvalue(model, L_diag, Ytest.block(0, i, N, 1)); + double high_eigenvalue_dense = eigenvalue_scaled_hessian(model, L_dense, Ytest.block(0, i, N, 1)); + double high_eigenvalue_diag = eigenvalue_scaled_hessian(model, L_diag, Ytest.block(0, i, N, 1)); c_dense = std::max(c_dense, std::sqrt(high_eigenvalue_dense / low_eigenvalue_dense)); c_diag = std::max(c_diag, std::sqrt(high_eigenvalue_diag / low_eigenvalue_diag)); @@ -175,8 +211,10 @@ namespace stan { } else { if(use_dense) { covar = dense; + covar_is_diagonal = false; } else { covar = diag; + covar_is_diagonal = true; } } } From 800069148021c77f173248da14af273ecd3a11c9 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Sat, 20 Apr 2019 15:28:13 -0700 Subject: [PATCH 03/20] Changed 'switching' to 'auto' --- ...ing_adaptation.hpp => auto_adaptation.hpp} | 10 +-- ...itching_e_metric.hpp => auto_e_metric.hpp} | 30 +++---- ...switching_e_point.hpp => auto_e_point.hpp} | 14 +-- ...ching_e_nuts.hpp => adapt_auto_e_nuts.hpp} | 24 +++--- .../{switching_e_nuts.hpp => auto_e_nuts.hpp} | 14 +-- ..._adapter.hpp => stepsize_auto_adapter.hpp} | 20 ++--- ..._e_adapt.hpp => hmc_nuts_auto_e_adapt.hpp} | 86 +++++++++---------- 7 files changed, 99 insertions(+), 99 deletions(-) rename src/stan/mcmc/{switching_adaptation.hpp => auto_adaptation.hpp} (96%) rename src/stan/mcmc/hmc/hamiltonians/{switching_e_metric.hpp => auto_e_metric.hpp} (61%) rename src/stan/mcmc/hmc/hamiltonians/{switching_e_point.hpp => auto_e_point.hpp} (83%) rename src/stan/mcmc/hmc/nuts/{adapt_switching_e_nuts.hpp => adapt_auto_e_nuts.hpp} (63%) rename src/stan/mcmc/hmc/nuts/{switching_e_nuts.hpp => auto_e_nuts.hpp} (55%) rename src/stan/mcmc/{stepsize_switching_adapter.hpp => stepsize_auto_adapter.hpp} (60%) rename src/stan/services/sample/{hmc_nuts_switching_e_adapt.hpp => hmc_nuts_auto_e_adapt.hpp} (74%) diff --git a/src/stan/mcmc/switching_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp similarity index 96% rename from src/stan/mcmc/switching_adaptation.hpp rename to src/stan/mcmc/auto_adaptation.hpp index 088543c27b..9ff87bed57 100644 --- a/src/stan/mcmc/switching_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -1,5 +1,5 @@ -#ifndef STAN_MCMC_SWITCHING_ADAPTATION_HPP -#define STAN_MCMC_SWITCHING_ADAPTATION_HPP +#ifndef STAN_MCMC_AUTO_ADAPTATION_HPP +#define STAN_MCMC_AUTO_ADAPTATION_HPP #include #include @@ -19,9 +19,9 @@ namespace stan { } }; - class switching_adaptation: public windowed_adaptation { + class auto_adaptation: public windowed_adaptation { public: - explicit switching_adaptation(int n) + explicit auto_adaptation(int n) : windowed_adaptation("covariance") {} /** @@ -62,7 +62,7 @@ namespace stan { double v_norm = v.norm(); double new_eval = v.dot(Av) / (v_norm * v_norm); if(std::abs(new_eval - eval) <= tol * std::abs(eval)) { - std::cout << "Converged at i = " << i << std::endl; + //std::cout << "Converged at i = " << i << std::endl; eval = new_eval; break; } diff --git a/src/stan/mcmc/hmc/hamiltonians/switching_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp similarity index 61% rename from src/stan/mcmc/hmc/hamiltonians/switching_e_metric.hpp rename to src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp index b86cf3a913..bdacb513e7 100644 --- a/src/stan/mcmc/hmc/hamiltonians/switching_e_metric.hpp +++ b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp @@ -1,10 +1,10 @@ -#ifndef STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_METRIC_HPP -#define STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_METRIC_HPP +#ifndef STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_METRIC_HPP +#define STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_METRIC_HPP #include #include #include -#include +#include #include #include #include @@ -14,33 +14,33 @@ namespace stan { // Euclidean manifold with dense metric template - class switching_e_metric - : public base_hamiltonian { + class auto_e_metric + : public base_hamiltonian { public: - explicit switching_e_metric(const Model& model) - : base_hamiltonian(model) {} + explicit auto_e_metric(const Model& model) + : base_hamiltonian(model) {} - double T(switching_e_point& z) { + double T(auto_e_point& z) { return 0.5 * z.p.transpose() * z.inv_e_metric_ * z.p; } - double tau(switching_e_point& z) { + double tau(auto_e_point& z) { return T(z); } - double phi(switching_e_point& z) { + double phi(auto_e_point& z) { return this->V(z); } - double dG_dt(switching_e_point& z, callbacks::logger& logger) { + double dG_dt(auto_e_point& z, callbacks::logger& logger) { return 2 * T(z) - z.q.dot(z.g); } - Eigen::VectorXd dtau_dq(switching_e_point& z, callbacks::logger& logger) { + Eigen::VectorXd dtau_dq(auto_e_point& z, callbacks::logger& logger) { return Eigen::VectorXd::Zero(this->model_.num_params_r()); } - Eigen::VectorXd dtau_dp(switching_e_point& z) { + Eigen::VectorXd dtau_dp(auto_e_point& z) { if(z.is_diagonal_) { return z.inv_e_metric_.diagonal().cwiseProduct(z.p); } else { @@ -48,11 +48,11 @@ namespace stan { } } - Eigen::VectorXd dphi_dq(switching_e_point& z, callbacks::logger& logger) { + Eigen::VectorXd dphi_dq(auto_e_point& z, callbacks::logger& logger) { return z.g; } - void sample_p(switching_e_point& z, BaseRNG& rng) { + void sample_p(auto_e_point& z, BaseRNG& rng) { typedef typename stan::math::index_type::type idx_t; boost::variate_generator > rand_gaus(rng, boost::normal_distribution<>()); diff --git a/src/stan/mcmc/hmc/hamiltonians/switching_e_point.hpp b/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp similarity index 83% rename from src/stan/mcmc/hmc/hamiltonians/switching_e_point.hpp rename to src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp index 3d72230fab..69ba07fb35 100644 --- a/src/stan/mcmc/hmc/hamiltonians/switching_e_point.hpp +++ b/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp @@ -1,5 +1,5 @@ -#ifndef STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_POINT_HPP -#define STAN_MCMC_HMC_HAMILTONIANS_SWITCHING_E_POINT_HPP +#ifndef STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_POINT_HPP +#define STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_POINT_HPP #include #include @@ -8,9 +8,9 @@ namespace stan { namespace mcmc { /** * Point in a phase space with a base - * Euclidean manifold with switching metric + * Euclidean manifold with auto metric */ - class switching_e_point: public ps_point { + class auto_e_point: public ps_point { public: /** * Inverse mass matrix. @@ -23,12 +23,12 @@ namespace stan { bool is_diagonal_; /** - * Construct a switching point in n-dimensional phase space + * Construct a auto point in n-dimensional phase space * with identity matrix as inverse mass matrix. * * @param n number of dimensions */ - explicit switching_e_point(int n) + explicit auto_e_point(int n) : ps_point(n), inv_e_metric_(n, n), is_diagonal_(true) { inv_e_metric_.setIdentity(); } @@ -38,7 +38,7 @@ namespace stan { * * @param z point to copy */ - switching_e_point(const switching_e_point& z) + auto_e_point(const auto_e_point& z) : ps_point(z), inv_e_metric_(z.inv_e_metric_.rows(), z.inv_e_metric_.cols()) { fast_matrix_copy_(inv_e_metric_, z.inv_e_metric_); diff --git a/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp similarity index 63% rename from src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp rename to src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp index b0a21c9e02..4335f61df9 100644 --- a/src/stan/mcmc/hmc/nuts/adapt_switching_e_nuts.hpp +++ b/src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp @@ -1,9 +1,9 @@ -#ifndef STAN_MCMC_HMC_NUTS_ADAPT_SWITCHING_E_NUTS_HPP -#define STAN_MCMC_HMC_NUTS_ADAPT_SWITCHING_E_NUTS_HPP +#ifndef STAN_MCMC_HMC_NUTS_ADAPT_AUTO_E_NUTS_HPP +#define STAN_MCMC_HMC_NUTS_ADAPT_AUTO_E_NUTS_HPP #include -#include -#include +#include +#include namespace stan { namespace mcmc { @@ -13,27 +13,27 @@ namespace stan { * dense metric and adaptive step size */ template - class adapt_switching_e_nuts : public switching_e_nuts, - public stepsize_switching_adapter { + class adapt_auto_e_nuts : public auto_e_nuts, + public stepsize_auto_adapter { protected: const Model& model_; public: - adapt_switching_e_nuts(const Model& model, BaseRNG& rng) - : model_(model), switching_e_nuts(model, rng), - stepsize_switching_adapter(model.num_params_r()) {} + adapt_auto_e_nuts(const Model& model, BaseRNG& rng) + : model_(model), auto_e_nuts(model, rng), + stepsize_auto_adapter(model.num_params_r()) {} - ~adapt_switching_e_nuts() {} + ~adapt_auto_e_nuts() {} sample transition(sample& init_sample, callbacks::logger& logger) { - sample s = switching_e_nuts::transition(init_sample, + sample s = auto_e_nuts::transition(init_sample, logger); if (this->adapt_flag_) { this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_, s.accept_stat()); - bool update = this->switching_adaptation_.learn_covariance( + bool update = this->auto_adaptation_.learn_covariance( model_, this->z_.inv_e_metric_, this->z_.is_diagonal_, diff --git a/src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp similarity index 55% rename from src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp rename to src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp index cd91df6fb8..4ed11dff61 100644 --- a/src/stan/mcmc/hmc/nuts/switching_e_nuts.hpp +++ b/src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp @@ -1,9 +1,9 @@ -#ifndef STAN_MCMC_HMC_NUTS_SWITCHING_E_NUTS_HPP -#define STAN_MCMC_HMC_NUTS_SWITCHING_E_NUTS_HPP +#ifndef STAN_MCMC_HMC_NUTS_AUTO_E_NUTS_HPP +#define STAN_MCMC_HMC_NUTS_AUTO_E_NUTS_HPP #include #include -#include +#include #include namespace stan { @@ -13,11 +13,11 @@ namespace stan { * with a Gaussian-Euclidean disintegration and dense metric */ template - class switching_e_nuts : public base_nuts { + class auto_e_nuts : public base_nuts { public: - switching_e_nuts(const Model& model, BaseRNG& rng) - : base_nuts(model, rng) { } }; diff --git a/src/stan/mcmc/stepsize_switching_adapter.hpp b/src/stan/mcmc/stepsize_auto_adapter.hpp similarity index 60% rename from src/stan/mcmc/stepsize_switching_adapter.hpp rename to src/stan/mcmc/stepsize_auto_adapter.hpp index aaa35ad623..0db72870d3 100644 --- a/src/stan/mcmc/stepsize_switching_adapter.hpp +++ b/src/stan/mcmc/stepsize_auto_adapter.hpp @@ -1,27 +1,27 @@ -#ifndef STAN_MCMC_STEPSIZE_SWITCHING_ADAPTER_HPP -#define STAN_MCMC_STEPSIZE_SWITCHING_ADAPTER_HPP +#ifndef STAN_MCMC_STEPSIZE_AUTO_ADAPTER_HPP +#define STAN_MCMC_STEPSIZE_AUTO_ADAPTER_HPP #include #include #include -#include +#include namespace stan { namespace mcmc { - class stepsize_switching_adapter: public base_adapter { + class stepsize_auto_adapter: public base_adapter { public: - explicit stepsize_switching_adapter(int n) - : switching_adaptation_(n) { + explicit stepsize_auto_adapter(int n) + : auto_adaptation_(n) { } stepsize_adaptation& get_stepsize_adaptation() { return stepsize_adaptation_; } - switching_adaptation& get_switching_adaptation() { - return switching_adaptation_; + auto_adaptation& get_auto_adaptation() { + return auto_adaptation_; } void set_window_params(unsigned int num_warmup, @@ -29,7 +29,7 @@ namespace stan { unsigned int term_buffer, unsigned int base_window, callbacks::logger& logger) { - switching_adaptation_.set_window_params(num_warmup, + auto_adaptation_.set_window_params(num_warmup, init_buffer, term_buffer, base_window, @@ -38,7 +38,7 @@ namespace stan { protected: stepsize_adaptation stepsize_adaptation_; - switching_adaptation switching_adaptation_; + auto_adaptation auto_adaptation_; }; } // mcmc diff --git a/src/stan/services/sample/hmc_nuts_switching_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp similarity index 74% rename from src/stan/services/sample/hmc_nuts_switching_e_adapt.hpp rename to src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp index dea87cc6c3..464880d316 100644 --- a/src/stan/services/sample/hmc_nuts_switching_e_adapt.hpp +++ b/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp @@ -1,5 +1,5 @@ -#ifndef STAN_SERVICES_SAMPLE_HMC_NUTS_SWITCHING_E_ADAPT_HPP -#define STAN_SERVICES_SAMPLE_HMC_NUTS_SWITCHING_E_ADAPT_HPP +#ifndef STAN_SERVICES_SAMPLE_HMC_NUTS_AUTO_E_ADAPT_HPP +#define STAN_SERVICES_SAMPLE_HMC_NUTS_AUTO_E_ADAPT_HPP #include #include @@ -8,7 +8,7 @@ #include #include #include -#include +#include #include #include #include @@ -54,21 +54,21 @@ namespace stan { * @return error_codes::OK if successful */ template - int hmc_nuts_switching_e_adapt(Model& model, stan::io::var_context& init, - stan::io::var_context& init_inv_metric, - unsigned int random_seed, unsigned int chain, - double init_radius, int num_warmup, - int num_samples, int num_thin, - bool save_warmup, int refresh, double stepsize, - double stepsize_jitter, int max_depth, - double delta, double gamma, double kappa, - double t0, unsigned int init_buffer, - unsigned int term_buffer, unsigned int window, - callbacks::interrupt& interrupt, - callbacks::logger& logger, - callbacks::writer& init_writer, - callbacks::writer& sample_writer, - callbacks::writer& diagnostic_writer) { + int hmc_nuts_auto_e_adapt(Model& model, stan::io::var_context& init, + stan::io::var_context& init_inv_metric, + unsigned int random_seed, unsigned int chain, + double init_radius, int num_warmup, + int num_samples, int num_thin, + bool save_warmup, int refresh, double stepsize, + double stepsize_jitter, int max_depth, + double delta, double gamma, double kappa, + double t0, unsigned int init_buffer, + unsigned int term_buffer, unsigned int window, + callbacks::interrupt& interrupt, + callbacks::logger& logger, + callbacks::writer& init_writer, + callbacks::writer& sample_writer, + callbacks::writer& diagnostic_writer) { boost::ecuyer1988 rng = util::create_rng(random_seed, chain); std::vector disc_vector; @@ -86,7 +86,7 @@ namespace stan { return error_codes::CONFIG; } - stan::mcmc::adapt_switching_e_nuts + stan::mcmc::adapt_auto_e_nuts sampler(model, rng); sampler.set_metric(inv_metric); @@ -145,34 +145,34 @@ namespace stan { * @return error_codes::OK if successful */ template - int hmc_nuts_switching_e_adapt(Model& model, stan::io::var_context& init, - unsigned int random_seed, unsigned int chain, - double init_radius, int num_warmup, - int num_samples, int num_thin, - bool save_warmup, int refresh, double stepsize, - double stepsize_jitter, int max_depth, - double delta, double gamma, double kappa, - double t0, unsigned int init_buffer, - unsigned int term_buffer, unsigned int window, - callbacks::interrupt& interrupt, - callbacks::logger& logger, - callbacks::writer& init_writer, - callbacks::writer& sample_writer, - callbacks::writer& diagnostic_writer) { + int hmc_nuts_auto_e_adapt(Model& model, stan::io::var_context& init, + unsigned int random_seed, unsigned int chain, + double init_radius, int num_warmup, + int num_samples, int num_thin, + bool save_warmup, int refresh, double stepsize, + double stepsize_jitter, int max_depth, + double delta, double gamma, double kappa, + double t0, unsigned int init_buffer, + unsigned int term_buffer, unsigned int window, + callbacks::interrupt& interrupt, + callbacks::logger& logger, + callbacks::writer& init_writer, + callbacks::writer& sample_writer, + callbacks::writer& diagnostic_writer) { stan::io::dump dmp = util::create_unit_e_dense_inv_metric(model.num_params_r()); stan::io::var_context& unit_e_metric = dmp; - return hmc_nuts_switching_e_adapt(model, init, unit_e_metric, - random_seed, chain, init_radius, - num_warmup, num_samples, num_thin, - save_warmup, refresh, - stepsize, stepsize_jitter, max_depth, - delta, gamma, kappa, t0, - init_buffer, term_buffer, window, - interrupt, logger, - init_writer, sample_writer, - diagnostic_writer); + return hmc_nuts_auto_e_adapt(model, init, unit_e_metric, + random_seed, chain, init_radius, + num_warmup, num_samples, num_thin, + save_warmup, refresh, + stepsize, stepsize_jitter, max_depth, + delta, gamma, kappa, t0, + init_buffer, term_buffer, window, + interrupt, logger, + init_writer, sample_writer, + diagnostic_writer); } } From b1205f3422e425c5bf5e56345869cc1c70f377ca Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Mon, 22 Apr 2019 15:51:47 -0700 Subject: [PATCH 04/20] Added try/catch so that if auto adaptation fails it falls back to diagonal --- src/stan/mcmc/auto_adaptation.hpp | 113 ++++++++++++++++-------------- 1 file changed, 61 insertions(+), 52 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index 9ff87bed57..daae43b0d0 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -152,71 +152,80 @@ namespace stan { for(int i = 0; i < qs_.size(); i++) Y.block(0, i, N, 1) = qs_[idxs[i]]; - bool use_dense = false; - for(auto state : { "selection", "refinement" }) { - Eigen::MatrixXd Ytrain; - Eigen::MatrixXd Ytest; - - if(state == "selection") { - int Ntest; - Ntest = int(0.2 * Y.cols()); - if(Ntest < 5) { - Ntest = 5; - } - - if(Y.cols() < 10) { - throw std::runtime_error("Each warmup stage must have at least 10 samples"); - } + try { + bool use_dense = false; + for(auto state : { "selection", "refinement" }) { + Eigen::MatrixXd Ytrain; + Eigen::MatrixXd Ytest; + + if(state == "selection") { + int Ntest; + Ntest = int(0.2 * Y.cols()); + if(Ntest < 5) { + Ntest = 5; + } + + if(Y.cols() < 10) { + throw std::runtime_error("Each warmup stage must have at least 10 samples"); + } - Ytrain = Y.block(0, 0, N, Y.cols() - Ntest); - Ytest = Y.block(0, Ytrain.cols(), N, Ntest); - } else { - Ytrain = Y; - } + Ytrain = Y.block(0, 0, N, Y.cols() - Ntest); + Ytest = Y.block(0, Ytrain.cols(), N, Ntest); + } else { + Ytrain = Y; + } - Eigen::MatrixXd cov_train = covariance(Ytrain); - Eigen::MatrixXd cov_test = covariance(Ytest); + Eigen::MatrixXd cov_train = covariance(Ytrain); + Eigen::MatrixXd cov_test = covariance(Ytest); - Eigen::MatrixXd dense = (N / (N + 5.0)) * cov_train + - 1e-3 * (5.0 / (N + 5.0)) * Eigen::MatrixXd::Identity(cov_train.rows(), cov_train.cols()); - Eigen::MatrixXd diag = dense.diagonal().asDiagonal(); + Eigen::MatrixXd dense = (N / (N + 5.0)) * cov_train + + 1e-3 * (5.0 / (N + 5.0)) * Eigen::MatrixXd::Identity(cov_train.rows(), cov_train.cols()); + Eigen::MatrixXd diag = dense.diagonal().asDiagonal(); - covar = dense; + covar = dense; - if(state == "selection") { - Eigen::MatrixXd L_dense = dense.llt().matrixL(); - Eigen::MatrixXd L_diag = diag.diagonal().array().sqrt().matrix().asDiagonal(); + if(state == "selection") { + Eigen::MatrixXd L_dense = dense.llt().matrixL(); + Eigen::MatrixXd L_diag = diag.diagonal().array().sqrt().matrix().asDiagonal(); - double low_eigenvalue_dense = -1.0 / eigenvalue_scaled_covariance(L_dense, cov_test); - double low_eigenvalue_diag = -1.0 / eigenvalue_scaled_covariance(L_diag, cov_test); + double low_eigenvalue_dense = -1.0 / eigenvalue_scaled_covariance(L_dense, cov_test); + double low_eigenvalue_diag = -1.0 / eigenvalue_scaled_covariance(L_diag, cov_test); - double c_dense = 0.0; - double c_diag = 0.0; - for(int i = 0; i < 5; i++) { - double high_eigenvalue_dense = eigenvalue_scaled_hessian(model, L_dense, Ytest.block(0, i, N, 1)); - double high_eigenvalue_diag = eigenvalue_scaled_hessian(model, L_diag, Ytest.block(0, i, N, 1)); + double c_dense = 0.0; + double c_diag = 0.0; + for(int i = 0; i < 5; i++) { + double high_eigenvalue_dense = eigenvalue_scaled_hessian(model, L_dense, Ytest.block(0, i, N, 1)); + double high_eigenvalue_diag = eigenvalue_scaled_hessian(model, L_diag, Ytest.block(0, i, N, 1)); - c_dense = std::max(c_dense, std::sqrt(high_eigenvalue_dense / low_eigenvalue_dense)); - c_diag = std::max(c_diag, std::sqrt(high_eigenvalue_diag / low_eigenvalue_diag)); - } + c_dense = std::max(c_dense, std::sqrt(high_eigenvalue_dense / low_eigenvalue_dense)); + c_diag = std::max(c_diag, std::sqrt(high_eigenvalue_diag / low_eigenvalue_diag)); + } - std::cout << "adapt: " << adapt_window_counter_ << ", which: dense, max: " << c_dense << std::endl; - std::cout << "adapt: " << adapt_window_counter_ << ", which: diag, max: " << c_diag << std::endl; + std::cout << "adapt: " << adapt_window_counter_ << ", which: dense, max: " << c_dense << std::endl; + std::cout << "adapt: " << adapt_window_counter_ << ", which: diag, max: " << c_diag << std::endl; - if(c_dense < c_diag) { - use_dense = true; - } else { - use_dense = false; - } - } else { - if(use_dense) { - covar = dense; - covar_is_diagonal = false; + if(c_dense < c_diag) { + use_dense = true; + } else { + use_dense = false; + } } else { - covar = diag; - covar_is_diagonal = true; + if(use_dense) { + covar = dense; + covar_is_diagonal = false; + } else { + covar = diag; + covar_is_diagonal = true; + } } } + } catch(const std::exception& e) { + std::cout << e.what() << std::endl; + std::cout << "Exception while using auto adaptation, falling back to diagonal" << std::endl; + Eigen::MatrixXd cov = covariance(Y); + covar = ((M / (M + 5.0)) * cov.diagonal() + + 1e-3 * (5.0 / (M + 5.0)) * Eigen::VectorXd::Ones(cov.cols())).asDiagonal(); + covar_is_diagonal = true; } ++adapt_window_counter_; From baacc34304202facdabe5593f531b996728ccc67 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 11 Sep 2019 15:11:35 -0400 Subject: [PATCH 05/20] Fixed bug with regularization of automatically picked metric. Added tests for automatically picked metric --- src/stan/mcmc/auto_adaptation.hpp | 135 ++++++++------ .../good/model/correlated_gaussian.stan | 7 + .../good/model/independent_gaussian.stan | 7 + .../test-models/good/model/known_hessian.stan | 7 + ...ation_learn_covariance_pick_dense_test.cpp | 64 +++++++ ...tation_learn_covariance_pick_diag_test.cpp | 63 +++++++ src/test/unit/mcmc/auto_adaptation_test.cpp | 170 ++++++++++++++++++ 7 files changed, 396 insertions(+), 57 deletions(-) create mode 100644 src/test/test-models/good/model/correlated_gaussian.stan create mode 100644 src/test/test-models/good/model/independent_gaussian.stan create mode 100644 src/test/test-models/good/model/known_hessian.stan create mode 100644 src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp create mode 100644 src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp create mode 100644 src/test/unit/mcmc/auto_adaptation_test.cpp diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index daae43b0d0..f7f85164ec 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -1,7 +1,7 @@ #ifndef STAN_MCMC_AUTO_ADAPTATION_HPP #define STAN_MCMC_AUTO_ADAPTATION_HPP -#include +#include #include #include @@ -19,20 +19,23 @@ namespace stan { } }; - class auto_adaptation: public windowed_adaptation { - public: - explicit auto_adaptation(int n) - : windowed_adaptation("covariance") {} - + namespace internal { /** - * Compute the covariance of data in Y. Rows of Y are different data. Columns of Y are different variables. + * Compute the covariance of data in Y. + * + * Columns of Y are different variables. Rows are different samples. + * + * When there is only one row in Y, return a covariance matrix of the expected + * size filled with zeros. * * @param Y Data * @return Covariance of Y */ Eigen::MatrixXd covariance(const Eigen::MatrixXd& Y) { - Eigen::MatrixXd centered = Y.colwise() - Y.rowwise().mean(); - return centered * centered.transpose() / std::max(centered.rows() - 1.0, 1.0); + stan::math::check_nonzero_size("covariance", "Y", Y); + + Eigen::MatrixXd centered = Y.rowwise() - Y.colwise().mean(); + return centered.transpose() * centered / std::max(centered.rows() - 1.0, 1.0); } /** @@ -48,31 +51,59 @@ namespace stan { * This means the returned eigenvalue might not be computed to full precision * * @param initial_guess Initial guess of the eigenvector of the largest eigenvalue - * @param max_iterations Maximum number of power iterations - * @param tol Relative tolerance + * @param[in,out] max_iterations Maximum number of power iterations, on return number of iterations used + * @param[in,out] tol Relative tolerance, on return the relative error in the eigenvalue estimate * @return Largest magnitude eigenvalue of operator f */ template - double power_method(F& f, const Eigen::VectorXd& initial_guess, int max_iterations, double tol) { + double power_method(F& f, const Eigen::VectorXd& initial_guess, int& max_iterations, double& tol) { Eigen::VectorXd v = initial_guess; double eval = 0.0; + Eigen::VectorXd Av = f(v); + stan::math::check_matching_sizes("power_method", "matrix vector product", Av, "vector", v); - for(int i = 0; i < max_iterations; i++) { - Eigen::VectorXd Av = f(v); + int i = 0; + for(; i < max_iterations; ++i) { double v_norm = v.norm(); double new_eval = v.dot(Av) / (v_norm * v_norm); - if(std::abs(new_eval - eval) <= tol * std::abs(eval)) { - //std::cout << "Converged at i = " << i << std::endl; + if(i == max_iterations - 1 || std::abs(new_eval - eval) <= tol * std::abs(eval)) { + tol = std::abs(new_eval - eval) / std::abs(eval); eval = new_eval; + max_iterations = i + 1; break; } + eval = new_eval; v = Av / Av.norm(); + + Av = f(v); } return eval; } + /** + * Compute the largest eigenvalue of the sample covariance rescaled by a metric, + * that is, the largest eigenvalue of L^{-1} \Sigma L^{-T} + * + * @param L Cholesky decomposition of Metric + * @param Sigma Sample covariance + * @return Largest eigenvalue + */ + double eigenvalue_scaled_covariance(const Eigen::MatrixXd& L, const Eigen::MatrixXd& Sigma) { + Eigen::MatrixXd S = L.template triangularView(). + solve(L.template triangularView().solve(Sigma).transpose()).transpose(); + + auto Sx = [&](Eigen::VectorXd x) -> Eigen::VectorXd { + return S * x; + }; + + int max_iterations = 100; + double tol = 1e-3; + + return internal::power_method(Sx, Eigen::VectorXd::Random(Sigma.cols()), max_iterations, tol); + } + /** * Compute the largest eigenvalue of the Hessian of the log density rescaled by a metric, * that is, the largest eigenvalue of L^T \nabla^2_{qq} H(q) L @@ -99,29 +130,18 @@ namespace stan { stan::math::gradient(log_prob_wrapper_covar(model), q - dr / 2.0, lp, grad2); return L.transpose() * (grad1 - grad2) / dx; }; - - return power_method(hessian_vector, Eigen::VectorXd::Random(q.size()), 100, 1e-3); - } - - /** - * Compute the largest eigenvalue of the sample covariance rescaled by a metric, - * that is, the largest eigenvalue of L^{-T} \Sigma L^{-1} - * - * @param L Cholesky decomposition of Metric - * @param Sigma Sample covariance - * @return Largest eigenvalue - */ - double eigenvalue_scaled_covariance(const Eigen::MatrixXd& L, const Eigen::MatrixXd& Sigma) { - Eigen::MatrixXd S = L.template triangularView(). - solve(L.template triangularView().solve(Sigma).transpose()).transpose(); - - auto Sx = [&](Eigen::VectorXd x) -> Eigen::VectorXd { - return S * x; - }; + + int max_iterations = 100; + double tol = 1e-3; - return power_method(Sx, Eigen::VectorXd::Random(Sigma.cols()), 100, 1e-3); + return internal::power_method(hessian_vector, Eigen::VectorXd::Random(q.size()), max_iterations, tol); } + } + class auto_adaptation: public windowed_adaptation { + public: + explicit auto_adaptation(int n) + : windowed_adaptation("covariance") {} /** * Update the metric if at the end of an adaptation window. * @@ -140,17 +160,17 @@ namespace stan { if (end_adaptation_window()) { compute_next_window(); - int N = q.size(); - int M = qs_.size(); + int M = q.size(); + int N = qs_.size(); - Eigen::MatrixXd Y = Eigen::MatrixXd::Zero(N, M); - std::vector idxs(M); + Eigen::MatrixXd Y = Eigen::MatrixXd::Zero(M, N); + std::vector idxs(N); for(int i = 0; i < qs_.size(); i++) idxs[i] = i; std::random_shuffle(idxs.begin(), idxs.end()); for(int i = 0; i < qs_.size(); i++) - Y.block(0, i, N, 1) = qs_[idxs[i]]; + Y.block(0, i, M, 1) = qs_[idxs[i]]; try { bool use_dense = false; @@ -159,27 +179,28 @@ namespace stan { Eigen::MatrixXd Ytest; if(state == "selection") { - int Ntest; - Ntest = int(0.2 * Y.cols()); - if(Ntest < 5) { - Ntest = 5; + int Mtest; + Mtest = int(0.2 * Y.cols()); + if(Mtest < 5) { + Mtest = 5; } if(Y.cols() < 10) { throw std::runtime_error("Each warmup stage must have at least 10 samples"); } - Ytrain = Y.block(0, 0, N, Y.cols() - Ntest); - Ytest = Y.block(0, Ytrain.cols(), N, Ntest); + Ytrain = Y.block(0, 0, M, Y.cols() - Mtest); + Ytest = Y.block(0, Ytrain.cols(), M, Mtest); } else { Ytrain = Y; } - Eigen::MatrixXd cov_train = covariance(Ytrain); - Eigen::MatrixXd cov_test = covariance(Ytest); - + Eigen::MatrixXd cov_train = (Ytrain.cols() > 0) ? internal::covariance(Ytrain.transpose()) : Eigen::MatrixXd::Zero(M, M); + Eigen::MatrixXd cov_test = (Ytest.cols() > 0) ? internal::covariance(Ytest.transpose()) : Eigen::MatrixXd::Zero(M, M); + Eigen::MatrixXd dense = (N / (N + 5.0)) * cov_train + 1e-3 * (5.0 / (N + 5.0)) * Eigen::MatrixXd::Identity(cov_train.rows(), cov_train.cols()); + Eigen::MatrixXd diag = dense.diagonal().asDiagonal(); covar = dense; @@ -188,14 +209,14 @@ namespace stan { Eigen::MatrixXd L_dense = dense.llt().matrixL(); Eigen::MatrixXd L_diag = diag.diagonal().array().sqrt().matrix().asDiagonal(); - double low_eigenvalue_dense = -1.0 / eigenvalue_scaled_covariance(L_dense, cov_test); - double low_eigenvalue_diag = -1.0 / eigenvalue_scaled_covariance(L_diag, cov_test); + double low_eigenvalue_dense = -1.0 / internal::eigenvalue_scaled_covariance(L_dense, cov_test); + double low_eigenvalue_diag = -1.0 / internal::eigenvalue_scaled_covariance(L_diag, cov_test); double c_dense = 0.0; double c_diag = 0.0; for(int i = 0; i < 5; i++) { - double high_eigenvalue_dense = eigenvalue_scaled_hessian(model, L_dense, Ytest.block(0, i, N, 1)); - double high_eigenvalue_diag = eigenvalue_scaled_hessian(model, L_diag, Ytest.block(0, i, N, 1)); + double high_eigenvalue_dense = internal::eigenvalue_scaled_hessian(model, L_dense, Ytest.block(0, i, M, 1)); + double high_eigenvalue_diag = internal::eigenvalue_scaled_hessian(model, L_diag, Ytest.block(0, i, M, 1)); c_dense = std::max(c_dense, std::sqrt(high_eigenvalue_dense / low_eigenvalue_dense)); c_diag = std::max(c_diag, std::sqrt(high_eigenvalue_diag / low_eigenvalue_diag)); @@ -222,9 +243,9 @@ namespace stan { } catch(const std::exception& e) { std::cout << e.what() << std::endl; std::cout << "Exception while using auto adaptation, falling back to diagonal" << std::endl; - Eigen::MatrixXd cov = covariance(Y); - covar = ((M / (M + 5.0)) * cov.diagonal() - + 1e-3 * (5.0 / (M + 5.0)) * Eigen::VectorXd::Ones(cov.cols())).asDiagonal(); + Eigen::MatrixXd cov = internal::covariance(Y.transpose()); + covar = ((N / (N + 5.0)) * cov.diagonal() + + 1e-3 * (5.0 / (N + 5.0)) * Eigen::VectorXd::Ones(cov.cols())).asDiagonal(); covar_is_diagonal = true; } diff --git a/src/test/test-models/good/model/correlated_gaussian.stan b/src/test/test-models/good/model/correlated_gaussian.stan new file mode 100644 index 0000000000..fc55ebd83e --- /dev/null +++ b/src/test/test-models/good/model/correlated_gaussian.stan @@ -0,0 +1,7 @@ +parameters { + vector[2] x; +} + +model { + x ~ multi_normal([0.0, 0.0], [[1.0, 0.99], [0.99, 1.0]]); +} diff --git a/src/test/test-models/good/model/independent_gaussian.stan b/src/test/test-models/good/model/independent_gaussian.stan new file mode 100644 index 0000000000..451a1251c8 --- /dev/null +++ b/src/test/test-models/good/model/independent_gaussian.stan @@ -0,0 +1,7 @@ +parameters { + vector[2] x; +} + +model { + x ~ normal(0.0, 1.0); +} diff --git a/src/test/test-models/good/model/known_hessian.stan b/src/test/test-models/good/model/known_hessian.stan new file mode 100644 index 0000000000..da35f67443 --- /dev/null +++ b/src/test/test-models/good/model/known_hessian.stan @@ -0,0 +1,7 @@ +parameters { + real x[3]; +} + +model { + x ~ normal(0, 1); +} diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp new file mode 100644 index 0000000000..0a8bd02139 --- /dev/null +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp @@ -0,0 +1,64 @@ +#include +#include +#include +#include + +TEST(McmcVarAdaptation, learn_covariance_pick_dense) { + std::fstream data_stream(std::string("").c_str(), std::fstream::in); + stan::io::dump data_var_context(data_stream); + data_stream.close(); + + std::stringstream output; + correlated_gaussian_model_namespace::correlated_gaussian_model + correlated_gaussian_model(data_var_context, &output); + + stan::test::unit::instrumented_logger logger; + + const int M = 2; + const int N = 20; + Eigen::MatrixXd qs(N, M); + qs << 0.256173753306128, -0.0238087093098673, + -1.63218152810157, -1.5309929638363, + -0.812451331685826, -1.15062373620068, + -1.49814775191801, -1.51310110681996, + 0.738630631536685, 1.03588205799336, + 0.472288580035284, 0.250286770328584, + -1.63634486169493, -1.6222798835089, + -0.400790615207103, -0.337669147200631, + -0.568779612417544, -0.424833495378187, + 0.103690913176746, 0.272885200284842, + -0.453017424229528, -0.504634004215693, + 3.34484533887237, 3.29418872328382, + -1.3376507113241, -1.32724775403694, + -0.137543235057544, -0.0290938109919368, + -1.58194496352741, -1.39338740677379, + 0.312166136194586, 0.336989933768233, + -0.628941448228566, -0.850758612234264, + -0.766816808981044, -0.645020468024267, + -0.75078110234827, -0.502544092120385, + -0.00694807494461906, -0.186748159558166; + + Eigen::MatrixXd covar(M, M); + bool covar_is_diagonal; + + Eigen::MatrixXd target_covar(M, M); + + target_covar << 1.0311414783609130, 1.0100577463968425, + 1.0100577463968425, 1.0148380697138280; + + stan::mcmc::auto_adaptation adapter(M); + adapter.set_window_params(50, 0, 0, N, logger); + + for (int i = 0; i < N; ++i) { + Eigen::VectorXd q = qs.block(i, 0, 1, M).transpose(); + adapter.learn_covariance(correlated_gaussian_model, covar, covar_is_diagonal, q); + } + + for (int i = 0; i < covar.size(); ++i) { + EXPECT_FLOAT_EQ(target_covar(i), covar(i)); + } + + EXPECT_EQ(covar_is_diagonal, false); + + EXPECT_EQ(0, logger.call_count()); +} diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp new file mode 100644 index 0000000000..26f3da56ce --- /dev/null +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp @@ -0,0 +1,63 @@ +#include +#include +#include +#include + +TEST(McmcVarAdaptation, learn_covariance_pick_diagonal) { + std::fstream data_stream(std::string("").c_str(), std::fstream::in); + stan::io::dump data_var_context(data_stream); + data_stream.close(); + + std::stringstream output; + independent_gaussian_model_namespace::independent_gaussian_model + independent_gaussian_model(data_var_context, &output); + + stan::test::unit::instrumented_logger logger; + + const int M = 2; + const int N = 20; + Eigen::MatrixXd qs(N, M); + qs << 0.607446257145326, 0.338465765807058, + 1.47389672467345, -1.0577986841911, + 1.02886652895522, 0.364277500948572, + 0.492316893603469, 2.19693408641558, + -0.931854393410476, 1.62634580968769, + -0.443145375724188, 0.902790875582656, + 0.517782110245233, -1.56724331755861, + -1.7556390097031, 0.310274990315213, + 0.0394975482340945, 0.366999438969482, + 1.29372950054929, 0.361369734821582, + -0.258301497542829, 0.166994731172984, + 0.492639248874412, -0.659502589885556, + 0.913729457222598, 1.99580706461809, + 0.669655370469707, -0.509028392475839, + -0.626041244059129, -0.771981104624195, + -0.842385483586737, 0.337166271031201, + 0.548177804329155, -0.0462961925005498, + 0.955748803092952, 1.3141117316189, + 0.335670079140694, 1.09112083087171, + 0.759245358940033, -1.11318882201676; + + Eigen::MatrixXd covar(M, M); + bool covar_is_diagonal; + + Eigen::MatrixXd target_covar(M, M); + + target_covar << 0.55350038163333048, 0.0, 0.0, 0.86122545968912112; + + stan::mcmc::auto_adaptation adapter(M); + adapter.set_window_params(50, 0, 0, N, logger); + + for (int i = 0; i < N; ++i) { + Eigen::VectorXd q = qs.block(i, 0, 1, M).transpose(); + adapter.learn_covariance(independent_gaussian_model, covar, covar_is_diagonal, q); + } + + for (int i = 0; i < covar.size(); ++i) { + EXPECT_FLOAT_EQ(target_covar(i), covar(i)); + } + + EXPECT_EQ(covar_is_diagonal, true); + + EXPECT_EQ(0, logger.call_count()); +} diff --git a/src/test/unit/mcmc/auto_adaptation_test.cpp b/src/test/unit/mcmc/auto_adaptation_test.cpp new file mode 100644 index 0000000000..6989016325 --- /dev/null +++ b/src/test/unit/mcmc/auto_adaptation_test.cpp @@ -0,0 +1,170 @@ +#include +#include +#include +#include + +TEST(McmcAutoAdaptation, test_covariance_zero_rows_zero_cols) { + Eigen::MatrixXd X1(0, 5); + + EXPECT_THROW(stan::mcmc::internal::covariance(X1), std::invalid_argument); + + Eigen::MatrixXd X2(1, 0); + + EXPECT_THROW(stan::mcmc::internal::covariance(X2), std::invalid_argument); +} + +TEST(McmcAutoAdaptation, test_covariance_one_row_one_col) { + Eigen::MatrixXd X1(1, 2); + Eigen::MatrixXd X2(3, 1); + + X1 << 1.0, 2.0; + X2 << 1.0, 2.0, 3.0; + + Eigen::MatrixXd cov1 = stan::mcmc::internal::covariance(X1); + Eigen::MatrixXd cov2 = stan::mcmc::internal::covariance(X2); + + ASSERT_EQ(cov1.rows(), 2); + ASSERT_EQ(cov1.cols(), 2); + + ASSERT_EQ(cov2.rows(), 1); + ASSERT_EQ(cov2.cols(), 1); + + for(int i = 0; i < cov1.size(); ++i) { + ASSERT_FLOAT_EQ(cov1(i), 0.0); + } + + ASSERT_FLOAT_EQ(cov2(0), 1.0); +} + +TEST(McmcAutoAdaptation, test_covariance) { + Eigen::MatrixXd X1(3, 2); + Eigen::MatrixXd X2(2, 3); + + X1 << 0.0, -1.0, 0.5, -2.7, 3.0, 5.0; + X2 << 0.0, 3, -2.7, 0.5, -1, 5.0; + + Eigen::MatrixXd cov1 = stan::mcmc::internal::covariance(X1); + Eigen::MatrixXd cov2 = stan::mcmc::internal::covariance(X2); + + Eigen::MatrixXd cov1_ref(2, 2); + Eigen::MatrixXd cov2_ref(3, 3); + + cov1_ref << 2.5833333333333335, 6.0666666666666664, + 6.0666666666666664, 16.3633333333333333; + + cov2_ref << 0.125, -1.0, 1.925, + -1.000, 8.0, -15.4, + 1.925, -15.4, 29.645; + + ASSERT_EQ(cov1.rows(), cov1_ref.rows()); + ASSERT_EQ(cov1.cols(), cov1_ref.cols()); + + ASSERT_EQ(cov2.rows(), cov2_ref.rows()); + ASSERT_EQ(cov2.cols(), cov2_ref.cols()); + + for(int i = 0; i < cov1_ref.size(); ++i) { + ASSERT_FLOAT_EQ(cov1(i), cov1_ref(i)); + } + + for(int i = 0; i < cov2_ref.size(); ++i) { + ASSERT_FLOAT_EQ(cov2(i), cov2_ref(i)); + } +} + +TEST(McmcAutoAdaptation, power_method) { + Eigen::MatrixXd X(2, 2); + Eigen::VectorXd x0(2); + + X << 2.0, 0.5, 0.5, 1.0; + x0 << 1.0, 0.0; + + const int max_iterations = 10; + const double tol = 1e-10; + + auto Av = [&](const Eigen::VectorXd& v) { return X * v; }; + + int max_iterations_1 = max_iterations; + double tol_1 = tol; + + double eval = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); + + EXPECT_FLOAT_EQ(eval, 2.20710678118654746); +} + +TEST(McmcAutoAdaptation, power_method_tol_check) { + Eigen::MatrixXd X(2, 2); + Eigen::VectorXd x0(2); + + X << 2.0, 0.5, 0.5, 1.0; + x0 << 1.0, 0.0; + + const int max_iterations = 1000; + const double tol = 1e-12; + + auto Av = [&](const Eigen::VectorXd& v) { return X * v; }; + + int max_iterations_1 = max_iterations; + double tol_1 = tol; + double eval = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); + + EXPECT_LT(tol_1, tol); +} + +TEST(McmcAutoAdaptation, power_method_iter_check) { + Eigen::MatrixXd X(2, 2); + Eigen::VectorXd x0(2); + + X << 2.0, 0.5, 0.5, 1.0; + x0 << 1.0, 0.0; + + const int max_iterations = 10; + const double tol = 1e-50; + + auto Av = [&](const Eigen::VectorXd& v) { return X * v; }; + + int max_iterations_1 = max_iterations; + double tol_1 = tol; + double eval = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); + + EXPECT_GT(tol_1, tol); + EXPECT_EQ(max_iterations_1, max_iterations); +} + +// The checks in here are very coarse because eigenvalue_scaled_covariance +// only estimates things with low precision +TEST(McmcAutoAdaptation, eigenvalue_scaled_covariance) { + Eigen::MatrixXd L(2, 2), Sigma(2, 2); + + L << 1.0, 0.0, 0.5, 1.0; + Sigma << 2.0, 0.7, 0.7, 1.3; + + double eval = stan::mcmc::internal::eigenvalue_scaled_covariance(L, Sigma); + + EXPECT_LT(std::abs(eval - 2.0908326913195983) / eval, 1e-2); + + L << 2.0, 0.0, 0.7, 1.3; + + eval = stan::mcmc::internal::eigenvalue_scaled_covariance(L, Sigma); + + EXPECT_LT(std::abs(eval - 0.62426035502958577) / eval, 1e-2); +} + +// The checks in here are very coarse because eigenvalue_scaled_hessian +// only estimates things with low precision +TEST(McmcAutoAdaptation, eigenvalue_scaled_hessian) { + std::fstream data_stream(std::string("").c_str(), std::fstream::in); + stan::io::dump data_var_context(data_stream); + data_stream.close(); + + std::stringstream output; + known_hessian_model_namespace::known_hessian_model known_hessian_model(data_var_context, &output); + + Eigen::MatrixXd L(3, 3); + Eigen::VectorXd q(3); + L << 2.0, 0.0, 0.0, 0.7, 1.3, 0.0, -1.5, 2.0, 4.0; + q << 0.0, 0.0, 0.0; + + double eval = stan::mcmc::internal::eigenvalue_scaled_hessian(known_hessian_model, L, q); + + EXPECT_LT(std::abs(eval - 22.8141075806892850) / eval, 1e-2); +} From 989b396660f1e541ab85d94e71582b9381332a8f Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 29 Oct 2020 15:43:54 -0400 Subject: [PATCH 06/20] Updated test file with new model signature (Issue #2814) --- src/stan/mcmc/auto_adaptation.hpp | 2 +- .../mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index f7f85164ec..1863f17d6c 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -1,7 +1,7 @@ #ifndef STAN_MCMC_AUTO_ADAPTATION_HPP #define STAN_MCMC_AUTO_ADAPTATION_HPP -#include +#include #include #include diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp index 26f3da56ce..488a12b1f7 100644 --- a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp @@ -10,7 +10,7 @@ TEST(McmcVarAdaptation, learn_covariance_pick_diagonal) { std::stringstream output; independent_gaussian_model_namespace::independent_gaussian_model - independent_gaussian_model(data_var_context, &output); + independent_gaussian_model(data_var_context, 0, &output); stan::test::unit::instrumented_logger logger; From 0f36d293b7016449e7bd82ed24506daa56e8c203 Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 29 Oct 2020 16:06:57 -0400 Subject: [PATCH 07/20] Fix cpplint errors (Issue #2814) --- .../mcmc/hmc/hamiltonians/auto_e_metric.hpp | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp index bdacb513e7..1bca1961c8 100644 --- a/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp +++ b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp @@ -41,11 +41,11 @@ namespace stan { } Eigen::VectorXd dtau_dp(auto_e_point& z) { - if(z.is_diagonal_) { - return z.inv_e_metric_.diagonal().cwiseProduct(z.p); - } else { - return z.inv_e_metric_ * z.p; - } + if(z.is_diagonal_) { + return z.inv_e_metric_.diagonal().cwiseProduct(z.p); + } else { + return z.inv_e_metric_ * z.p; + } } Eigen::VectorXd dphi_dq(auto_e_point& z, callbacks::logger& logger) { @@ -57,17 +57,17 @@ namespace stan { boost::variate_generator > rand_gaus(rng, boost::normal_distribution<>()); - if(z.is_diagonal_) { - for (int i = 0; i < z.p.size(); ++i) - z.p(i) = rand_gaus() / sqrt(z.inv_e_metric_(i, i)); - } else { - Eigen::VectorXd u(z.p.size()); + if(z.is_diagonal_) { + for (int i = 0; i < z.p.size(); ++i) + z.p(i) = rand_gaus() / sqrt(z.inv_e_metric_(i, i)); + } else { + Eigen::VectorXd u(z.p.size()); - for (idx_t i = 0; i < u.size(); ++i) - u(i) = rand_gaus(); + for (idx_t i = 0; i < u.size(); ++i) + u(i) = rand_gaus(); - z.p = z.inv_e_metric_.llt().matrixU().solve(u); - } + z.p = z.inv_e_metric_.llt().matrixU().solve(u); + } } }; From 5d37ce41d0d0f160b2883655760aa8540ab5a60a Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 29 Oct 2020 16:13:38 -0400 Subject: [PATCH 08/20] Moved format in front of lint --- Jenkinsfile | 78 ++++++++++++++++++++++++++--------------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/Jenkinsfile b/Jenkinsfile index 99185b0468..f343e5545e 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -74,45 +74,6 @@ pipeline { } } } - stage('Linting & Doc checks') { - agent any - steps { - script { - sh "printenv" - retry(3) { checkout scm } - sh """ - make math-revert - make clean-all - git clean -xffd - """ - utils.checkout_pr("math", "lib/stan_math", params.math_pr) - stash 'StanSetup' - setupCXX(true, env.GCC) - parallel( - CppLint: { sh "make cpplint" }, - API_docs: { sh 'make doxygen' }, - ) - } - } - post { - always { - - recordIssues id: "lint_doc_checks", - name: "Linting & Doc checks", - enabledForFailure: true, - aggregatingResults : true, - tools: [ - cppLint(id: "cpplint", name: "Linting & Doc checks@CPPLINT") - ], - blameDisabled: false, - qualityGates: [[threshold: 1, type: 'TOTAL', unstable: true]], - healthy: 10, unhealthy: 100, minimumSeverity: 'HIGH', - referenceJobName: env.BRANCH_NAME - - deleteDir() - } - } - } stage("Clang-format") { agent any steps { @@ -161,6 +122,45 @@ pipeline { } } } + stage('Linting & Doc checks') { + agent any + steps { + script { + sh "printenv" + retry(3) { checkout scm } + sh """ + make math-revert + make clean-all + git clean -xffd + """ + utils.checkout_pr("math", "lib/stan_math", params.math_pr) + stash 'StanSetup' + setupCXX(true, env.GCC) + parallel( + CppLint: { sh "make cpplint" }, + API_docs: { sh 'make doxygen' }, + ) + } + } + post { + always { + + recordIssues id: "lint_doc_checks", + name: "Linting & Doc checks", + enabledForFailure: true, + aggregatingResults : true, + tools: [ + cppLint(id: "cpplint", name: "Linting & Doc checks@CPPLINT") + ], + blameDisabled: false, + qualityGates: [[threshold: 1, type: 'TOTAL', unstable: true]], + healthy: 10, unhealthy: 100, minimumSeverity: 'HIGH', + referenceJobName: env.BRANCH_NAME + + deleteDir() + } + } + } stage('Verify changes') { agent { label 'linux' } steps { From 70c5417959945ecb5230f02cf23eeaabff7b2fa8 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 29 Oct 2020 20:18:28 +0000 Subject: [PATCH 09/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- src/stan/mcmc/auto_adaptation.hpp | 541 ++++++++++-------- .../mcmc/hmc/hamiltonians/auto_e_metric.hpp | 95 ++- .../mcmc/hmc/hamiltonians/auto_e_point.hpp | 129 ++--- src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp | 83 ++- src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp | 29 +- src/stan/mcmc/stepsize_auto_adapter.hpp | 63 +- .../services/sample/hmc_nuts_auto_e_adapt.hpp | 279 +++++---- ...ation_learn_covariance_pick_dense_test.cpp | 44 +- ...tation_learn_covariance_pick_diag_test.cpp | 40 +- src/test/unit/mcmc/auto_adaptation_test.cpp | 45 +- 10 files changed, 670 insertions(+), 678 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index 1863f17d6c..95627263b3 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -7,264 +7,305 @@ namespace stan { - namespace mcmc { - template - struct log_prob_wrapper_covar { - const Model& model_; - log_prob_wrapper_covar(const Model& model) : model_(model) {} - - template - T operator()(const Eigen::Matrix& q) const { - return model_.template log_prob(const_cast& >(q), &std::cout); - } - }; - - namespace internal { - /** - * Compute the covariance of data in Y. - * - * Columns of Y are different variables. Rows are different samples. - * - * When there is only one row in Y, return a covariance matrix of the expected - * size filled with zeros. - * - * @param Y Data - * @return Covariance of Y - */ - Eigen::MatrixXd covariance(const Eigen::MatrixXd& Y) { - stan::math::check_nonzero_size("covariance", "Y", Y); - - Eigen::MatrixXd centered = Y.rowwise() - Y.colwise().mean(); - return centered.transpose() * centered / std::max(centered.rows() - 1.0, 1.0); - } +namespace mcmc { +template +struct log_prob_wrapper_covar { + const Model& model_; + log_prob_wrapper_covar(const Model& model) : model_(model) {} + + template + T operator()(const Eigen::Matrix& q) const { + return model_.template log_prob( + const_cast&>(q), &std::cout); + } +}; + +namespace internal { +/** + * Compute the covariance of data in Y. + * + * Columns of Y are different variables. Rows are different samples. + * + * When there is only one row in Y, return a covariance matrix of the expected + * size filled with zeros. + * + * @param Y Data + * @return Covariance of Y + */ +Eigen::MatrixXd covariance(const Eigen::MatrixXd& Y) { + stan::math::check_nonzero_size("covariance", "Y", Y); + + Eigen::MatrixXd centered = Y.rowwise() - Y.colwise().mean(); + return centered.transpose() * centered / std::max(centered.rows() - 1.0, 1.0); +} + +/** + * Compute the largest magnitude eigenvalue of a symmetric matrix using the + * power method. The function f should return the product of that matrix with an + * abitrary vector. + * + * f should take one Eigen::VectorXd argument, x, and return the product of a + * matrix with x as an Eigen::VectorXd argument of the same size. + * + * The eigenvalue is estimated iteratively. If the kth estimate is e_k, then the + * function returns when either abs(e_{k + 1} - e_k) < tol * abs(e_k) or the + * maximum number of iterations have been performed + * + * This means the returned eigenvalue might not be computed to full precision + * + * @param initial_guess Initial guess of the eigenvector of the largest + * eigenvalue + * @param[in,out] max_iterations Maximum number of power iterations, on return + * number of iterations used + * @param[in,out] tol Relative tolerance, on return the relative error in the + * eigenvalue estimate + * @return Largest magnitude eigenvalue of operator f + */ +template +double power_method(F& f, const Eigen::VectorXd& initial_guess, + int& max_iterations, double& tol) { + Eigen::VectorXd v = initial_guess; + double eval = 0.0; + Eigen::VectorXd Av = f(v); + stan::math::check_matching_sizes("power_method", "matrix vector product", Av, + "vector", v); + + int i = 0; + for (; i < max_iterations; ++i) { + double v_norm = v.norm(); + double new_eval = v.dot(Av) / (v_norm * v_norm); + if (i == max_iterations - 1 + || std::abs(new_eval - eval) <= tol * std::abs(eval)) { + tol = std::abs(new_eval - eval) / std::abs(eval); + eval = new_eval; + max_iterations = i + 1; + break; + } - /** - * Compute the largest magnitude eigenvalue of a symmetric matrix using the power method. The function f - * should return the product of that matrix with an abitrary vector. - * - * f should take one Eigen::VectorXd argument, x, and return the product of a matrix with x as - * an Eigen::VectorXd argument of the same size. - * - * The eigenvalue is estimated iteratively. If the kth estimate is e_k, then the function returns when - * either abs(e_{k + 1} - e_k) < tol * abs(e_k) or the maximum number of iterations have been performed - * - * This means the returned eigenvalue might not be computed to full precision - * - * @param initial_guess Initial guess of the eigenvector of the largest eigenvalue - * @param[in,out] max_iterations Maximum number of power iterations, on return number of iterations used - * @param[in,out] tol Relative tolerance, on return the relative error in the eigenvalue estimate - * @return Largest magnitude eigenvalue of operator f - */ - template - double power_method(F& f, const Eigen::VectorXd& initial_guess, int& max_iterations, double& tol) { - Eigen::VectorXd v = initial_guess; - double eval = 0.0; - Eigen::VectorXd Av = f(v); - stan::math::check_matching_sizes("power_method", "matrix vector product", Av, "vector", v); - - int i = 0; - for(; i < max_iterations; ++i) { - double v_norm = v.norm(); - double new_eval = v.dot(Av) / (v_norm * v_norm); - if(i == max_iterations - 1 || std::abs(new_eval - eval) <= tol * std::abs(eval)) { - tol = std::abs(new_eval - eval) / std::abs(eval); - eval = new_eval; - max_iterations = i + 1; - break; - } - - eval = new_eval; - v = Av / Av.norm(); - - Av = f(v); - } - - return eval; + eval = new_eval; + v = Av / Av.norm(); + + Av = f(v); + } + + return eval; +} + +/** + * Compute the largest eigenvalue of the sample covariance rescaled by a metric, + * that is, the largest eigenvalue of L^{-1} \Sigma L^{-T} + * + * @param L Cholesky decomposition of Metric + * @param Sigma Sample covariance + * @return Largest eigenvalue + */ +double eigenvalue_scaled_covariance(const Eigen::MatrixXd& L, + const Eigen::MatrixXd& Sigma) { + Eigen::MatrixXd S = L.template triangularView() + .solve(L.template triangularView() + .solve(Sigma) + .transpose()) + .transpose(); + + auto Sx = [&](Eigen::VectorXd x) -> Eigen::VectorXd { return S * x; }; + + int max_iterations = 100; + double tol = 1e-3; + + return internal::power_method(Sx, Eigen::VectorXd::Random(Sigma.cols()), + max_iterations, tol); +} + +/** + * Compute the largest eigenvalue of the Hessian of the log density rescaled by + * a metric, that is, the largest eigenvalue of L^T \nabla^2_{qq} H(q) L + * + * @tparam Model Type of model + * @param model Defines the log density + * @param q Point around which to compute the Hessian + * @param L Cholesky decomposition of Metric + * @return Largest eigenvalue + */ +template +double eigenvalue_scaled_hessian(const Model& model, const Eigen::MatrixXd& L, + const Eigen::VectorXd& q) { + Eigen::VectorXd eigenvalues; + Eigen::MatrixXd eigenvectors; + + auto hessian_vector = [&](const Eigen::VectorXd& x) -> Eigen::VectorXd { + double lp; + Eigen::VectorXd grad1; + Eigen::VectorXd grad2; + // stan::math::hessian_times_vector(log_prob_wrapper_covar(model), q, + // x, lp, Ax); + double dx = 1e-5; + Eigen::VectorXd dr = L * x * dx; + stan::math::gradient(log_prob_wrapper_covar(model), q + dr / 2.0, lp, + grad1); + stan::math::gradient(log_prob_wrapper_covar(model), q - dr / 2.0, lp, + grad2); + return L.transpose() * (grad1 - grad2) / dx; + }; + + int max_iterations = 100; + double tol = 1e-3; + + return internal::power_method( + hessian_vector, Eigen::VectorXd::Random(q.size()), max_iterations, tol); +} +} // namespace internal + +class auto_adaptation : public windowed_adaptation { + public: + explicit auto_adaptation(int n) : windowed_adaptation("covariance") {} + /** + * Update the metric if at the end of an adaptation window. + * + * @tparam Model Type of model + * @param model Defines the log density + * @param covar[out] New metric + * @param covar_is_diagonal[out] Set to true if metric is diagonal, false + * otherwise + * @param q New MCMC draw + * @return True if this was the end of an adaptation window, false otherwise + */ + template + bool learn_covariance(const Model& model, Eigen::MatrixXd& covar, + bool& covar_is_diagonal, const Eigen::VectorXd& q) { + if (adaptation_window()) + qs_.push_back(q); + + if (end_adaptation_window()) { + compute_next_window(); + + int M = q.size(); + int N = qs_.size(); + + Eigen::MatrixXd Y = Eigen::MatrixXd::Zero(M, N); + std::vector idxs(N); + for (int i = 0; i < qs_.size(); i++) + idxs[i] = i; + + std::random_shuffle(idxs.begin(), idxs.end()); + for (int i = 0; i < qs_.size(); i++) + Y.block(0, i, M, 1) = qs_[idxs[i]]; + + try { + bool use_dense = false; + for (auto state : {"selection", "refinement"}) { + Eigen::MatrixXd Ytrain; + Eigen::MatrixXd Ytest; + + if (state == "selection") { + int Mtest; + Mtest = int(0.2 * Y.cols()); + if (Mtest < 5) { + Mtest = 5; + } + + if (Y.cols() < 10) { + throw std::runtime_error( + "Each warmup stage must have at least 10 samples"); + } + + Ytrain = Y.block(0, 0, M, Y.cols() - Mtest); + Ytest = Y.block(0, Ytrain.cols(), M, Mtest); + } else { + Ytrain = Y; + } + + Eigen::MatrixXd cov_train + = (Ytrain.cols() > 0) ? internal::covariance(Ytrain.transpose()) + : Eigen::MatrixXd::Zero(M, M); + Eigen::MatrixXd cov_test + = (Ytest.cols() > 0) ? internal::covariance(Ytest.transpose()) + : Eigen::MatrixXd::Zero(M, M); + + Eigen::MatrixXd dense + = (N / (N + 5.0)) * cov_train + + 1e-3 * (5.0 / (N + 5.0)) + * Eigen::MatrixXd::Identity(cov_train.rows(), + cov_train.cols()); + + Eigen::MatrixXd diag = dense.diagonal().asDiagonal(); + + covar = dense; + + if (state == "selection") { + Eigen::MatrixXd L_dense = dense.llt().matrixL(); + Eigen::MatrixXd L_diag + = diag.diagonal().array().sqrt().matrix().asDiagonal(); + + double low_eigenvalue_dense + = -1.0 + / internal::eigenvalue_scaled_covariance(L_dense, cov_test); + double low_eigenvalue_diag + = -1.0 + / internal::eigenvalue_scaled_covariance(L_diag, cov_test); + + double c_dense = 0.0; + double c_diag = 0.0; + for (int i = 0; i < 5; i++) { + double high_eigenvalue_dense + = internal::eigenvalue_scaled_hessian( + model, L_dense, Ytest.block(0, i, M, 1)); + double high_eigenvalue_diag = internal::eigenvalue_scaled_hessian( + model, L_diag, Ytest.block(0, i, M, 1)); + + c_dense = std::max(c_dense, std::sqrt(high_eigenvalue_dense + / low_eigenvalue_dense)); + c_diag = std::max(c_diag, std::sqrt(high_eigenvalue_diag + / low_eigenvalue_diag)); + } + + std::cout << "adapt: " << adapt_window_counter_ + << ", which: dense, max: " << c_dense << std::endl; + std::cout << "adapt: " << adapt_window_counter_ + << ", which: diag, max: " << c_diag << std::endl; + + if (c_dense < c_diag) { + use_dense = true; + } else { + use_dense = false; + } + } else { + if (use_dense) { + covar = dense; + covar_is_diagonal = false; + } else { + covar = diag; + covar_is_diagonal = true; + } + } + } + } catch (const std::exception& e) { + std::cout << e.what() << std::endl; + std::cout + << "Exception while using auto adaptation, falling back to diagonal" + << std::endl; + Eigen::MatrixXd cov = internal::covariance(Y.transpose()); + covar = ((N / (N + 5.0)) * cov.diagonal() + + 1e-3 * (5.0 / (N + 5.0)) * Eigen::VectorXd::Ones(cov.cols())) + .asDiagonal(); + covar_is_diagonal = true; } - /** - * Compute the largest eigenvalue of the sample covariance rescaled by a metric, - * that is, the largest eigenvalue of L^{-1} \Sigma L^{-T} - * - * @param L Cholesky decomposition of Metric - * @param Sigma Sample covariance - * @return Largest eigenvalue - */ - double eigenvalue_scaled_covariance(const Eigen::MatrixXd& L, const Eigen::MatrixXd& Sigma) { - Eigen::MatrixXd S = L.template triangularView(). - solve(L.template triangularView().solve(Sigma).transpose()).transpose(); - - auto Sx = [&](Eigen::VectorXd x) -> Eigen::VectorXd { - return S * x; - }; - - int max_iterations = 100; - double tol = 1e-3; - - return internal::power_method(Sx, Eigen::VectorXd::Random(Sigma.cols()), max_iterations, tol); - } + ++adapt_window_counter_; + qs_.clear(); - /** - * Compute the largest eigenvalue of the Hessian of the log density rescaled by a metric, - * that is, the largest eigenvalue of L^T \nabla^2_{qq} H(q) L - * - * @tparam Model Type of model - * @param model Defines the log density - * @param q Point around which to compute the Hessian - * @param L Cholesky decomposition of Metric - * @return Largest eigenvalue - */ - template - double eigenvalue_scaled_hessian(const Model& model, const Eigen::MatrixXd& L, const Eigen::VectorXd& q) { - Eigen::VectorXd eigenvalues; - Eigen::MatrixXd eigenvectors; - - auto hessian_vector = [&](const Eigen::VectorXd& x) -> Eigen::VectorXd { - double lp; - Eigen::VectorXd grad1; - Eigen::VectorXd grad2; - //stan::math::hessian_times_vector(log_prob_wrapper_covar(model), q, x, lp, Ax); - double dx = 1e-5; - Eigen::VectorXd dr = L * x * dx; - stan::math::gradient(log_prob_wrapper_covar(model), q + dr / 2.0, lp, grad1); - stan::math::gradient(log_prob_wrapper_covar(model), q - dr / 2.0, lp, grad2); - return L.transpose() * (grad1 - grad2) / dx; - }; - - int max_iterations = 100; - double tol = 1e-3; - - return internal::power_method(hessian_vector, Eigen::VectorXd::Random(q.size()), max_iterations, tol); - } + return true; } - class auto_adaptation: public windowed_adaptation { - public: - explicit auto_adaptation(int n) - : windowed_adaptation("covariance") {} - /** - * Update the metric if at the end of an adaptation window. - * - * @tparam Model Type of model - * @param model Defines the log density - * @param covar[out] New metric - * @param covar_is_diagonal[out] Set to true if metric is diagonal, false otherwise - * @param q New MCMC draw - * @return True if this was the end of an adaptation window, false otherwise - */ - template - bool learn_covariance(const Model& model, Eigen::MatrixXd& covar, bool& covar_is_diagonal, const Eigen::VectorXd& q) { - if (adaptation_window()) - qs_.push_back(q); - - if (end_adaptation_window()) { - compute_next_window(); - - int M = q.size(); - int N = qs_.size(); - - Eigen::MatrixXd Y = Eigen::MatrixXd::Zero(M, N); - std::vector idxs(N); - for(int i = 0; i < qs_.size(); i++) - idxs[i] = i; - - std::random_shuffle(idxs.begin(), idxs.end()); - for(int i = 0; i < qs_.size(); i++) - Y.block(0, i, M, 1) = qs_[idxs[i]]; - - try { - bool use_dense = false; - for(auto state : { "selection", "refinement" }) { - Eigen::MatrixXd Ytrain; - Eigen::MatrixXd Ytest; - - if(state == "selection") { - int Mtest; - Mtest = int(0.2 * Y.cols()); - if(Mtest < 5) { - Mtest = 5; - } - - if(Y.cols() < 10) { - throw std::runtime_error("Each warmup stage must have at least 10 samples"); - } - - Ytrain = Y.block(0, 0, M, Y.cols() - Mtest); - Ytest = Y.block(0, Ytrain.cols(), M, Mtest); - } else { - Ytrain = Y; - } - - Eigen::MatrixXd cov_train = (Ytrain.cols() > 0) ? internal::covariance(Ytrain.transpose()) : Eigen::MatrixXd::Zero(M, M); - Eigen::MatrixXd cov_test = (Ytest.cols() > 0) ? internal::covariance(Ytest.transpose()) : Eigen::MatrixXd::Zero(M, M); - - Eigen::MatrixXd dense = (N / (N + 5.0)) * cov_train + - 1e-3 * (5.0 / (N + 5.0)) * Eigen::MatrixXd::Identity(cov_train.rows(), cov_train.cols()); - - Eigen::MatrixXd diag = dense.diagonal().asDiagonal(); - - covar = dense; - - if(state == "selection") { - Eigen::MatrixXd L_dense = dense.llt().matrixL(); - Eigen::MatrixXd L_diag = diag.diagonal().array().sqrt().matrix().asDiagonal(); - - double low_eigenvalue_dense = -1.0 / internal::eigenvalue_scaled_covariance(L_dense, cov_test); - double low_eigenvalue_diag = -1.0 / internal::eigenvalue_scaled_covariance(L_diag, cov_test); - - double c_dense = 0.0; - double c_diag = 0.0; - for(int i = 0; i < 5; i++) { - double high_eigenvalue_dense = internal::eigenvalue_scaled_hessian(model, L_dense, Ytest.block(0, i, M, 1)); - double high_eigenvalue_diag = internal::eigenvalue_scaled_hessian(model, L_diag, Ytest.block(0, i, M, 1)); - - c_dense = std::max(c_dense, std::sqrt(high_eigenvalue_dense / low_eigenvalue_dense)); - c_diag = std::max(c_diag, std::sqrt(high_eigenvalue_diag / low_eigenvalue_diag)); - } - - std::cout << "adapt: " << adapt_window_counter_ << ", which: dense, max: " << c_dense << std::endl; - std::cout << "adapt: " << adapt_window_counter_ << ", which: diag, max: " << c_diag << std::endl; - - if(c_dense < c_diag) { - use_dense = true; - } else { - use_dense = false; - } - } else { - if(use_dense) { - covar = dense; - covar_is_diagonal = false; - } else { - covar = diag; - covar_is_diagonal = true; - } - } - } - } catch(const std::exception& e) { - std::cout << e.what() << std::endl; - std::cout << "Exception while using auto adaptation, falling back to diagonal" << std::endl; - Eigen::MatrixXd cov = internal::covariance(Y.transpose()); - covar = ((N / (N + 5.0)) * cov.diagonal() - + 1e-3 * (5.0 / (N + 5.0)) * Eigen::VectorXd::Ones(cov.cols())).asDiagonal(); - covar_is_diagonal = true; - } - - ++adapt_window_counter_; - qs_.clear(); - - return true; - } - - ++adapt_window_counter_; - return false; - } + ++adapt_window_counter_; + return false; + } - protected: - std::vector< Eigen::VectorXd > qs_; - }; + protected: + std::vector qs_; +}; - } // mcmc +} // namespace mcmc -} // stan +} // namespace stan #endif diff --git a/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp index 1bca1961c8..cff85d1dca 100644 --- a/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp +++ b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp @@ -10,67 +10,62 @@ #include namespace stan { - namespace mcmc { +namespace mcmc { - // Euclidean manifold with dense metric - template - class auto_e_metric - : public base_hamiltonian { - public: - explicit auto_e_metric(const Model& model) - : base_hamiltonian(model) {} +// Euclidean manifold with dense metric +template +class auto_e_metric : public base_hamiltonian { + public: + explicit auto_e_metric(const Model& model) + : base_hamiltonian(model) {} - double T(auto_e_point& z) { - return 0.5 * z.p.transpose() * z.inv_e_metric_ * z.p; - } + double T(auto_e_point& z) { + return 0.5 * z.p.transpose() * z.inv_e_metric_ * z.p; + } - double tau(auto_e_point& z) { - return T(z); - } + double tau(auto_e_point& z) { return T(z); } - double phi(auto_e_point& z) { - return this->V(z); - } + double phi(auto_e_point& z) { return this->V(z); } - double dG_dt(auto_e_point& z, callbacks::logger& logger) { - return 2 * T(z) - z.q.dot(z.g); - } + double dG_dt(auto_e_point& z, callbacks::logger& logger) { + return 2 * T(z) - z.q.dot(z.g); + } - Eigen::VectorXd dtau_dq(auto_e_point& z, callbacks::logger& logger) { - return Eigen::VectorXd::Zero(this->model_.num_params_r()); - } + Eigen::VectorXd dtau_dq(auto_e_point& z, callbacks::logger& logger) { + return Eigen::VectorXd::Zero(this->model_.num_params_r()); + } - Eigen::VectorXd dtau_dp(auto_e_point& z) { - if(z.is_diagonal_) { - return z.inv_e_metric_.diagonal().cwiseProduct(z.p); - } else { - return z.inv_e_metric_ * z.p; - } - } + Eigen::VectorXd dtau_dp(auto_e_point& z) { + if (z.is_diagonal_) { + return z.inv_e_metric_.diagonal().cwiseProduct(z.p); + } else { + return z.inv_e_metric_ * z.p; + } + } - Eigen::VectorXd dphi_dq(auto_e_point& z, callbacks::logger& logger) { - return z.g; - } + Eigen::VectorXd dphi_dq(auto_e_point& z, callbacks::logger& logger) { + return z.g; + } - void sample_p(auto_e_point& z, BaseRNG& rng) { - typedef typename stan::math::index_type::type idx_t; - boost::variate_generator > - rand_gaus(rng, boost::normal_distribution<>()); + void sample_p(auto_e_point& z, BaseRNG& rng) { + typedef typename stan::math::index_type::type idx_t; + boost::variate_generator > rand_gaus( + rng, boost::normal_distribution<>()); - if(z.is_diagonal_) { - for (int i = 0; i < z.p.size(); ++i) - z.p(i) = rand_gaus() / sqrt(z.inv_e_metric_(i, i)); - } else { - Eigen::VectorXd u(z.p.size()); + if (z.is_diagonal_) { + for (int i = 0; i < z.p.size(); ++i) + z.p(i) = rand_gaus() / sqrt(z.inv_e_metric_(i, i)); + } else { + Eigen::VectorXd u(z.p.size()); - for (idx_t i = 0; i < u.size(); ++i) - u(i) = rand_gaus(); + for (idx_t i = 0; i < u.size(); ++i) + u(i) = rand_gaus(); - z.p = z.inv_e_metric_.llt().matrixU().solve(u); - } - } - }; + z.p = z.inv_e_metric_.llt().matrixU().solve(u); + } + } +}; - } // mcmc -} // stan +} // namespace mcmc +} // namespace stan #endif diff --git a/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp b/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp index 69ba07fb35..13e24cc44f 100644 --- a/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp +++ b/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp @@ -5,77 +5,74 @@ #include namespace stan { - namespace mcmc { - /** - * Point in a phase space with a base - * Euclidean manifold with auto metric - */ - class auto_e_point: public ps_point { - public: - /** - * Inverse mass matrix. - */ - Eigen::MatrixXd inv_e_metric_; +namespace mcmc { +/** + * Point in a phase space with a base + * Euclidean manifold with auto metric + */ +class auto_e_point : public ps_point { + public: + /** + * Inverse mass matrix. + */ + Eigen::MatrixXd inv_e_metric_; - /** - * Is inv_e_metric_ diagonal or not - */ - bool is_diagonal_; + /** + * Is inv_e_metric_ diagonal or not + */ + bool is_diagonal_; - /** - * Construct a auto point in n-dimensional phase space - * with identity matrix as inverse mass matrix. - * - * @param n number of dimensions - */ - explicit auto_e_point(int n) - : ps_point(n), inv_e_metric_(n, n), is_diagonal_(true) { - inv_e_metric_.setIdentity(); - } + /** + * Construct a auto point in n-dimensional phase space + * with identity matrix as inverse mass matrix. + * + * @param n number of dimensions + */ + explicit auto_e_point(int n) + : ps_point(n), inv_e_metric_(n, n), is_diagonal_(true) { + inv_e_metric_.setIdentity(); + } - /** - * Copy constructor which does fast copy of inverse mass matrix. - * - * @param z point to copy - */ - auto_e_point(const auto_e_point& z) - : ps_point(z), inv_e_metric_(z.inv_e_metric_.rows(), - z.inv_e_metric_.cols()) { - fast_matrix_copy_(inv_e_metric_, z.inv_e_metric_); - is_diagonal_ = z.is_diagonal_; - } + /** + * Copy constructor which does fast copy of inverse mass matrix. + * + * @param z point to copy + */ + auto_e_point(const auto_e_point& z) + : ps_point(z), + inv_e_metric_(z.inv_e_metric_.rows(), z.inv_e_metric_.cols()) { + fast_matrix_copy_(inv_e_metric_, z.inv_e_metric_); + is_diagonal_ = z.is_diagonal_; + } - /** - * Set elements of mass matrix - * - * @param inv_e_metric initial mass matrix - */ - void - set_metric(const Eigen::MatrixXd& inv_e_metric) { - inv_e_metric_ = inv_e_metric; - is_diagonal_ = false; - } + /** + * Set elements of mass matrix + * + * @param inv_e_metric initial mass matrix + */ + void set_metric(const Eigen::MatrixXd& inv_e_metric) { + inv_e_metric_ = inv_e_metric; + is_diagonal_ = false; + } - /** - * Write elements of mass matrix to string and handoff to writer. - * - * @param writer Stan writer callback - */ - inline - void - write_metric(stan::callbacks::writer& writer) { - writer("Elements of inverse mass matrix:"); - for (int i = 0; i < inv_e_metric_.rows(); ++i) { - std::stringstream inv_e_metric_ss; - inv_e_metric_ss << inv_e_metric_(i, 0); - for (int j = 1; j < inv_e_metric_.cols(); ++j) - inv_e_metric_ss << ", " << inv_e_metric_(i, j); - writer(inv_e_metric_ss.str()); - } - } - }; + /** + * Write elements of mass matrix to string and handoff to writer. + * + * @param writer Stan writer callback + */ + inline void write_metric(stan::callbacks::writer& writer) { + writer("Elements of inverse mass matrix:"); + for (int i = 0; i < inv_e_metric_.rows(); ++i) { + std::stringstream inv_e_metric_ss; + inv_e_metric_ss << inv_e_metric_(i, 0); + for (int j = 1; j < inv_e_metric_.cols(); ++j) + inv_e_metric_ss << ", " << inv_e_metric_(i, j); + writer(inv_e_metric_ss.str()); + } + } +}; - } // mcmc -} // stan +} // namespace mcmc +} // namespace stan #endif diff --git a/src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp index 4335f61df9..56ebe6c160 100644 --- a/src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp +++ b/src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp @@ -6,55 +6,52 @@ #include namespace stan { - namespace mcmc { - /** - * The No-U-Turn sampler (NUTS) with multinomial sampling - * with a Gaussian-Euclidean disintegration and adaptive - * dense metric and adaptive step size - */ - template - class adapt_auto_e_nuts : public auto_e_nuts, - public stepsize_auto_adapter { - protected: - const Model& model_; - public: - adapt_auto_e_nuts(const Model& model, BaseRNG& rng) - : model_(model), auto_e_nuts(model, rng), +namespace mcmc { +/** + * The No-U-Turn sampler (NUTS) with multinomial sampling + * with a Gaussian-Euclidean disintegration and adaptive + * dense metric and adaptive step size + */ +template +class adapt_auto_e_nuts : public auto_e_nuts, + public stepsize_auto_adapter { + protected: + const Model& model_; + + public: + adapt_auto_e_nuts(const Model& model, BaseRNG& rng) + : model_(model), + auto_e_nuts(model, rng), stepsize_auto_adapter(model.num_params_r()) {} - ~adapt_auto_e_nuts() {} + ~adapt_auto_e_nuts() {} - sample - transition(sample& init_sample, callbacks::logger& logger) { - sample s = auto_e_nuts::transition(init_sample, - logger); + sample transition(sample& init_sample, callbacks::logger& logger) { + sample s = auto_e_nuts::transition(init_sample, logger); - if (this->adapt_flag_) { - this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_, - s.accept_stat()); + if (this->adapt_flag_) { + this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_, + s.accept_stat()); - bool update = this->auto_adaptation_.learn_covariance( - model_, - this->z_.inv_e_metric_, - this->z_.is_diagonal_, - this->z_.q); + bool update = this->auto_adaptation_.learn_covariance( + model_, this->z_.inv_e_metric_, this->z_.is_diagonal_, this->z_.q); - if (update) { - this->init_stepsize(logger); + if (update) { + this->init_stepsize(logger); - this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_)); - this->stepsize_adaptation_.restart(); - } - } - return s; + this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_)); + this->stepsize_adaptation_.restart(); } - - void disengage_adaptation() { - base_adapter::disengage_adaptation(); - this->stepsize_adaptation_.complete_adaptation(this->nom_epsilon_); - } - }; - - } // mcmc -} // stan + } + return s; + } + + void disengage_adaptation() { + base_adapter::disengage_adaptation(); + this->stepsize_adaptation_.complete_adaptation(this->nom_epsilon_); + } +}; + +} // namespace mcmc +} // namespace stan #endif diff --git a/src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp index 4ed11dff61..2974c8df07 100644 --- a/src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp +++ b/src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp @@ -7,20 +7,19 @@ #include namespace stan { - namespace mcmc { - /** - * The No-U-Turn sampler (NUTS) with multinomial sampling - * with a Gaussian-Euclidean disintegration and dense metric - */ - template - class auto_e_nuts : public base_nuts { - public: - auto_e_nuts(const Model& model, BaseRNG& rng) - : base_nuts(model, rng) { } - }; +namespace mcmc { +/** + * The No-U-Turn sampler (NUTS) with multinomial sampling + * with a Gaussian-Euclidean disintegration and dense metric + */ +template +class auto_e_nuts + : public base_nuts { + public: + auto_e_nuts(const Model& model, BaseRNG& rng) + : base_nuts(model, rng) {} +}; - } // mcmc -} // stan +} // namespace mcmc +} // namespace stan #endif diff --git a/src/stan/mcmc/stepsize_auto_adapter.hpp b/src/stan/mcmc/stepsize_auto_adapter.hpp index 0db72870d3..7645ec460a 100644 --- a/src/stan/mcmc/stepsize_auto_adapter.hpp +++ b/src/stan/mcmc/stepsize_auto_adapter.hpp @@ -8,41 +8,32 @@ namespace stan { - namespace mcmc { - - class stepsize_auto_adapter: public base_adapter { - public: - explicit stepsize_auto_adapter(int n) - : auto_adaptation_(n) { - } - - stepsize_adaptation& get_stepsize_adaptation() { - return stepsize_adaptation_; - } - - auto_adaptation& get_auto_adaptation() { - return auto_adaptation_; - } - - void set_window_params(unsigned int num_warmup, - unsigned int init_buffer, - unsigned int term_buffer, - unsigned int base_window, - callbacks::logger& logger) { - auto_adaptation_.set_window_params(num_warmup, - init_buffer, - term_buffer, - base_window, - logger); - } - - protected: - stepsize_adaptation stepsize_adaptation_; - auto_adaptation auto_adaptation_; - }; - - } // mcmc - -} // stan +namespace mcmc { + +class stepsize_auto_adapter : public base_adapter { + public: + explicit stepsize_auto_adapter(int n) : auto_adaptation_(n) {} + + stepsize_adaptation& get_stepsize_adaptation() { + return stepsize_adaptation_; + } + + auto_adaptation& get_auto_adaptation() { return auto_adaptation_; } + + void set_window_params(unsigned int num_warmup, unsigned int init_buffer, + unsigned int term_buffer, unsigned int base_window, + callbacks::logger& logger) { + auto_adaptation_.set_window_params(num_warmup, init_buffer, term_buffer, + base_window, logger); + } + + protected: + stepsize_adaptation stepsize_adaptation_; + auto_adaptation auto_adaptation_; +}; + +} // namespace mcmc + +} // namespace stan #endif diff --git a/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp index 464880d316..addd110753 100644 --- a/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp +++ b/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp @@ -16,166 +16,147 @@ #include namespace stan { - namespace services { - namespace sample { +namespace services { +namespace sample { - /** - * Runs HMC with NUTS with adaptation using dense Euclidean metric - * with a pre-specified Euclidean metric. - * - * @tparam Model Model class - * @param[in] model Input model to test (with data already instantiated) - * @param[in] init var context for initialization - * @param[in] init_inv_metric var context exposing an initial dense - inverse Euclidean metric (must be positive definite) - * @param[in] random_seed random seed for the random number generator - * @param[in] chain chain id to advance the pseudo random number generator - * @param[in] init_radius radius to initialize - * @param[in] num_warmup Number of warmup samples - * @param[in] num_samples Number of samples - * @param[in] num_thin Number to thin the samples - * @param[in] save_warmup Indicates whether to save the warmup iterations - * @param[in] refresh Controls the output - * @param[in] stepsize initial stepsize for discrete evolution - * @param[in] stepsize_jitter uniform random jitter of stepsize - * @param[in] max_depth Maximum tree depth - * @param[in] delta adaptation target acceptance statistic - * @param[in] gamma adaptation regularization scale - * @param[in] kappa adaptation relaxation exponent - * @param[in] t0 adaptation iteration offset - * @param[in] init_buffer width of initial fast adaptation interval - * @param[in] term_buffer width of final fast adaptation interval - * @param[in] window initial width of slow adaptation interval - * @param[in,out] interrupt Callback for interrupts - * @param[in,out] logger Logger for messages - * @param[in,out] init_writer Writer callback for unconstrained inits - * @param[in,out] sample_writer Writer for draws - * @param[in,out] diagnostic_writer Writer for diagnostic information - * @return error_codes::OK if successful - */ - template - int hmc_nuts_auto_e_adapt(Model& model, stan::io::var_context& init, - stan::io::var_context& init_inv_metric, - unsigned int random_seed, unsigned int chain, - double init_radius, int num_warmup, - int num_samples, int num_thin, - bool save_warmup, int refresh, double stepsize, - double stepsize_jitter, int max_depth, - double delta, double gamma, double kappa, - double t0, unsigned int init_buffer, - unsigned int term_buffer, unsigned int window, - callbacks::interrupt& interrupt, - callbacks::logger& logger, - callbacks::writer& init_writer, - callbacks::writer& sample_writer, - callbacks::writer& diagnostic_writer) { - boost::ecuyer1988 rng = util::create_rng(random_seed, chain); +/** + * Runs HMC with NUTS with adaptation using dense Euclidean metric + * with a pre-specified Euclidean metric. + * + * @tparam Model Model class + * @param[in] model Input model to test (with data already instantiated) + * @param[in] init var context for initialization + * @param[in] init_inv_metric var context exposing an initial dense + inverse Euclidean metric (must be positive definite) + * @param[in] random_seed random seed for the random number generator + * @param[in] chain chain id to advance the pseudo random number generator + * @param[in] init_radius radius to initialize + * @param[in] num_warmup Number of warmup samples + * @param[in] num_samples Number of samples + * @param[in] num_thin Number to thin the samples + * @param[in] save_warmup Indicates whether to save the warmup iterations + * @param[in] refresh Controls the output + * @param[in] stepsize initial stepsize for discrete evolution + * @param[in] stepsize_jitter uniform random jitter of stepsize + * @param[in] max_depth Maximum tree depth + * @param[in] delta adaptation target acceptance statistic + * @param[in] gamma adaptation regularization scale + * @param[in] kappa adaptation relaxation exponent + * @param[in] t0 adaptation iteration offset + * @param[in] init_buffer width of initial fast adaptation interval + * @param[in] term_buffer width of final fast adaptation interval + * @param[in] window initial width of slow adaptation interval + * @param[in,out] interrupt Callback for interrupts + * @param[in,out] logger Logger for messages + * @param[in,out] init_writer Writer callback for unconstrained inits + * @param[in,out] sample_writer Writer for draws + * @param[in,out] diagnostic_writer Writer for diagnostic information + * @return error_codes::OK if successful + */ +template +int hmc_nuts_auto_e_adapt( + Model& model, stan::io::var_context& init, + stan::io::var_context& init_inv_metric, unsigned int random_seed, + unsigned int chain, double init_radius, int num_warmup, int num_samples, + int num_thin, bool save_warmup, int refresh, double stepsize, + double stepsize_jitter, int max_depth, double delta, double gamma, + double kappa, double t0, unsigned int init_buffer, unsigned int term_buffer, + unsigned int window, callbacks::interrupt& interrupt, + callbacks::logger& logger, callbacks::writer& init_writer, + callbacks::writer& sample_writer, callbacks::writer& diagnostic_writer) { + boost::ecuyer1988 rng = util::create_rng(random_seed, chain); - std::vector disc_vector; - std::vector cont_vector - = util::initialize(model, init, rng, init_radius, true, - logger, init_writer); + std::vector disc_vector; + std::vector cont_vector = util::initialize( + model, init, rng, init_radius, true, logger, init_writer); - Eigen::MatrixXd inv_metric; - try { - inv_metric = - util::read_dense_inv_metric(init_inv_metric, model.num_params_r(), - logger); - util::validate_dense_inv_metric(inv_metric, logger); - } catch (const std::domain_error& e) { - return error_codes::CONFIG; - } - - stan::mcmc::adapt_auto_e_nuts - sampler(model, rng); + Eigen::MatrixXd inv_metric; + try { + inv_metric = util::read_dense_inv_metric(init_inv_metric, + model.num_params_r(), logger); + util::validate_dense_inv_metric(inv_metric, logger); + } catch (const std::domain_error& e) { + return error_codes::CONFIG; + } - sampler.set_metric(inv_metric); + stan::mcmc::adapt_auto_e_nuts sampler(model, rng); - sampler.set_nominal_stepsize(stepsize); - sampler.set_stepsize_jitter(stepsize_jitter); - sampler.set_max_depth(max_depth); + sampler.set_metric(inv_metric); - sampler.get_stepsize_adaptation().set_mu(log(10 * stepsize)); - sampler.get_stepsize_adaptation().set_delta(delta); - sampler.get_stepsize_adaptation().set_gamma(gamma); - sampler.get_stepsize_adaptation().set_kappa(kappa); - sampler.get_stepsize_adaptation().set_t0(t0); + sampler.set_nominal_stepsize(stepsize); + sampler.set_stepsize_jitter(stepsize_jitter); + sampler.set_max_depth(max_depth); - sampler.set_window_params(num_warmup, init_buffer, term_buffer, - window, logger); + sampler.get_stepsize_adaptation().set_mu(log(10 * stepsize)); + sampler.get_stepsize_adaptation().set_delta(delta); + sampler.get_stepsize_adaptation().set_gamma(gamma); + sampler.get_stepsize_adaptation().set_kappa(kappa); + sampler.get_stepsize_adaptation().set_t0(t0); - util::run_adaptive_sampler(sampler, model, cont_vector, num_warmup, - num_samples, num_thin, refresh, save_warmup, - rng, interrupt, logger, - sample_writer, diagnostic_writer); + sampler.set_window_params(num_warmup, init_buffer, term_buffer, window, + logger); - return error_codes::OK; - } + util::run_adaptive_sampler( + sampler, model, cont_vector, num_warmup, num_samples, num_thin, refresh, + save_warmup, rng, interrupt, logger, sample_writer, diagnostic_writer); - /** - * Runs HMC with NUTS with adaptation using dense Euclidean metric, - * with identity matrix as initial inv_metric. - * - * @tparam Model Model class - * @param[in] model Input model to test (with data already instantiated) - * @param[in] init var context for initialization - * @param[in] random_seed random seed for the random number generator - * @param[in] chain chain id to advance the pseudo random number generator - * @param[in] init_radius radius to initialize - * @param[in] num_warmup Number of warmup samples - * @param[in] num_samples Number of samples - * @param[in] num_thin Number to thin the samples - * @param[in] save_warmup Indicates whether to save the warmup iterations - * @param[in] refresh Controls the output - * @param[in] stepsize initial stepsize for discrete evolution - * @param[in] stepsize_jitter uniform random jitter of stepsize - * @param[in] max_depth Maximum tree depth - * @param[in] delta adaptation target acceptance statistic - * @param[in] gamma adaptation regularization scale - * @param[in] kappa adaptation relaxation exponent - * @param[in] t0 adaptation iteration offset - * @param[in] init_buffer width of initial fast adaptation interval - * @param[in] term_buffer width of final fast adaptation interval - * @param[in] window initial width of slow adaptation interval - * @param[in,out] interrupt Callback for interrupts - * @param[in,out] logger Logger for messages - * @param[in,out] init_writer Writer callback for unconstrained inits - * @param[in,out] sample_writer Writer for draws - * @param[in,out] diagnostic_writer Writer for diagnostic information - * @return error_codes::OK if successful - */ - template - int hmc_nuts_auto_e_adapt(Model& model, stan::io::var_context& init, - unsigned int random_seed, unsigned int chain, - double init_radius, int num_warmup, - int num_samples, int num_thin, - bool save_warmup, int refresh, double stepsize, - double stepsize_jitter, int max_depth, - double delta, double gamma, double kappa, - double t0, unsigned int init_buffer, - unsigned int term_buffer, unsigned int window, - callbacks::interrupt& interrupt, - callbacks::logger& logger, - callbacks::writer& init_writer, - callbacks::writer& sample_writer, - callbacks::writer& diagnostic_writer) { - stan::io::dump dmp = - util::create_unit_e_dense_inv_metric(model.num_params_r()); - stan::io::var_context& unit_e_metric = dmp; + return error_codes::OK; +} - return hmc_nuts_auto_e_adapt(model, init, unit_e_metric, - random_seed, chain, init_radius, - num_warmup, num_samples, num_thin, - save_warmup, refresh, - stepsize, stepsize_jitter, max_depth, - delta, gamma, kappa, t0, - init_buffer, term_buffer, window, - interrupt, logger, - init_writer, sample_writer, - diagnostic_writer); - } +/** + * Runs HMC with NUTS with adaptation using dense Euclidean metric, + * with identity matrix as initial inv_metric. + * + * @tparam Model Model class + * @param[in] model Input model to test (with data already instantiated) + * @param[in] init var context for initialization + * @param[in] random_seed random seed for the random number generator + * @param[in] chain chain id to advance the pseudo random number generator + * @param[in] init_radius radius to initialize + * @param[in] num_warmup Number of warmup samples + * @param[in] num_samples Number of samples + * @param[in] num_thin Number to thin the samples + * @param[in] save_warmup Indicates whether to save the warmup iterations + * @param[in] refresh Controls the output + * @param[in] stepsize initial stepsize for discrete evolution + * @param[in] stepsize_jitter uniform random jitter of stepsize + * @param[in] max_depth Maximum tree depth + * @param[in] delta adaptation target acceptance statistic + * @param[in] gamma adaptation regularization scale + * @param[in] kappa adaptation relaxation exponent + * @param[in] t0 adaptation iteration offset + * @param[in] init_buffer width of initial fast adaptation interval + * @param[in] term_buffer width of final fast adaptation interval + * @param[in] window initial width of slow adaptation interval + * @param[in,out] interrupt Callback for interrupts + * @param[in,out] logger Logger for messages + * @param[in,out] init_writer Writer callback for unconstrained inits + * @param[in,out] sample_writer Writer for draws + * @param[in,out] diagnostic_writer Writer for diagnostic information + * @return error_codes::OK if successful + */ +template +int hmc_nuts_auto_e_adapt( + Model& model, stan::io::var_context& init, unsigned int random_seed, + unsigned int chain, double init_radius, int num_warmup, int num_samples, + int num_thin, bool save_warmup, int refresh, double stepsize, + double stepsize_jitter, int max_depth, double delta, double gamma, + double kappa, double t0, unsigned int init_buffer, unsigned int term_buffer, + unsigned int window, callbacks::interrupt& interrupt, + callbacks::logger& logger, callbacks::writer& init_writer, + callbacks::writer& sample_writer, callbacks::writer& diagnostic_writer) { + stan::io::dump dmp + = util::create_unit_e_dense_inv_metric(model.num_params_r()); + stan::io::var_context& unit_e_metric = dmp; - } - } + return hmc_nuts_auto_e_adapt( + model, init, unit_e_metric, random_seed, chain, init_radius, num_warmup, + num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter, + max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window, + interrupt, logger, init_writer, sample_writer, diagnostic_writer); } + +} // namespace sample +} // namespace services +} // namespace stan #endif diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp index 0a8bd02139..b645e17496 100644 --- a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp @@ -10,48 +10,42 @@ TEST(McmcVarAdaptation, learn_covariance_pick_dense) { std::stringstream output; correlated_gaussian_model_namespace::correlated_gaussian_model - correlated_gaussian_model(data_var_context, &output); + correlated_gaussian_model(data_var_context, &output); stan::test::unit::instrumented_logger logger; const int M = 2; const int N = 20; Eigen::MatrixXd qs(N, M); - qs << 0.256173753306128, -0.0238087093098673, - -1.63218152810157, -1.5309929638363, - -0.812451331685826, -1.15062373620068, - -1.49814775191801, -1.51310110681996, - 0.738630631536685, 1.03588205799336, - 0.472288580035284, 0.250286770328584, - -1.63634486169493, -1.6222798835089, - -0.400790615207103, -0.337669147200631, - -0.568779612417544, -0.424833495378187, - 0.103690913176746, 0.272885200284842, - -0.453017424229528, -0.504634004215693, - 3.34484533887237, 3.29418872328382, - -1.3376507113241, -1.32724775403694, - -0.137543235057544, -0.0290938109919368, - -1.58194496352741, -1.39338740677379, - 0.312166136194586, 0.336989933768233, - -0.628941448228566, -0.850758612234264, - -0.766816808981044, -0.645020468024267, - -0.75078110234827, -0.502544092120385, - -0.00694807494461906, -0.186748159558166; + qs << 0.256173753306128, -0.0238087093098673, -1.63218152810157, + -1.5309929638363, -0.812451331685826, -1.15062373620068, + -1.49814775191801, -1.51310110681996, 0.738630631536685, 1.03588205799336, + 0.472288580035284, 0.250286770328584, -1.63634486169493, -1.6222798835089, + -0.400790615207103, -0.337669147200631, -0.568779612417544, + -0.424833495378187, 0.103690913176746, 0.272885200284842, + -0.453017424229528, -0.504634004215693, 3.34484533887237, + 3.29418872328382, -1.3376507113241, -1.32724775403694, -0.137543235057544, + -0.0290938109919368, -1.58194496352741, -1.39338740677379, + 0.312166136194586, 0.336989933768233, -0.628941448228566, + -0.850758612234264, -0.766816808981044, -0.645020468024267, + -0.75078110234827, -0.502544092120385, -0.00694807494461906, + -0.186748159558166; Eigen::MatrixXd covar(M, M); bool covar_is_diagonal; Eigen::MatrixXd target_covar(M, M); - target_covar << 1.0311414783609130, 1.0100577463968425, - 1.0100577463968425, 1.0148380697138280; + target_covar << 1.0311414783609130, 1.0100577463968425, 1.0100577463968425, + 1.0148380697138280; stan::mcmc::auto_adaptation adapter(M); adapter.set_window_params(50, 0, 0, N, logger); for (int i = 0; i < N; ++i) { Eigen::VectorXd q = qs.block(i, 0, 1, M).transpose(); - adapter.learn_covariance(correlated_gaussian_model, covar, covar_is_diagonal, q); + adapter.learn_covariance(correlated_gaussian_model, covar, + covar_is_diagonal, q); } for (int i = 0; i < covar.size(); ++i) { @@ -59,6 +53,6 @@ TEST(McmcVarAdaptation, learn_covariance_pick_dense) { } EXPECT_EQ(covar_is_diagonal, false); - + EXPECT_EQ(0, logger.call_count()); } diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp index 488a12b1f7..52aab40076 100644 --- a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp @@ -10,33 +10,26 @@ TEST(McmcVarAdaptation, learn_covariance_pick_diagonal) { std::stringstream output; independent_gaussian_model_namespace::independent_gaussian_model - independent_gaussian_model(data_var_context, 0, &output); + independent_gaussian_model(data_var_context, 0, &output); stan::test::unit::instrumented_logger logger; const int M = 2; const int N = 20; Eigen::MatrixXd qs(N, M); - qs << 0.607446257145326, 0.338465765807058, - 1.47389672467345, -1.0577986841911, - 1.02886652895522, 0.364277500948572, - 0.492316893603469, 2.19693408641558, - -0.931854393410476, 1.62634580968769, - -0.443145375724188, 0.902790875582656, - 0.517782110245233, -1.56724331755861, - -1.7556390097031, 0.310274990315213, - 0.0394975482340945, 0.366999438969482, - 1.29372950054929, 0.361369734821582, - -0.258301497542829, 0.166994731172984, - 0.492639248874412, -0.659502589885556, - 0.913729457222598, 1.99580706461809, - 0.669655370469707, -0.509028392475839, - -0.626041244059129, -0.771981104624195, - -0.842385483586737, 0.337166271031201, - 0.548177804329155, -0.0462961925005498, - 0.955748803092952, 1.3141117316189, - 0.335670079140694, 1.09112083087171, - 0.759245358940033, -1.11318882201676; + qs << 0.607446257145326, 0.338465765807058, 1.47389672467345, + -1.0577986841911, 1.02886652895522, 0.364277500948572, 0.492316893603469, + 2.19693408641558, -0.931854393410476, 1.62634580968769, + -0.443145375724188, 0.902790875582656, 0.517782110245233, + -1.56724331755861, -1.7556390097031, 0.310274990315213, + 0.0394975482340945, 0.366999438969482, 1.29372950054929, + 0.361369734821582, -0.258301497542829, 0.166994731172984, + 0.492639248874412, -0.659502589885556, 0.913729457222598, + 1.99580706461809, 0.669655370469707, -0.509028392475839, + -0.626041244059129, -0.771981104624195, -0.842385483586737, + 0.337166271031201, 0.548177804329155, -0.0462961925005498, + 0.955748803092952, 1.3141117316189, 0.335670079140694, 1.09112083087171, + 0.759245358940033, -1.11318882201676; Eigen::MatrixXd covar(M, M); bool covar_is_diagonal; @@ -50,7 +43,8 @@ TEST(McmcVarAdaptation, learn_covariance_pick_diagonal) { for (int i = 0; i < N; ++i) { Eigen::VectorXd q = qs.block(i, 0, 1, M).transpose(); - adapter.learn_covariance(independent_gaussian_model, covar, covar_is_diagonal, q); + adapter.learn_covariance(independent_gaussian_model, covar, + covar_is_diagonal, q); } for (int i = 0; i < covar.size(); ++i) { @@ -58,6 +52,6 @@ TEST(McmcVarAdaptation, learn_covariance_pick_diagonal) { } EXPECT_EQ(covar_is_diagonal, true); - + EXPECT_EQ(0, logger.call_count()); } diff --git a/src/test/unit/mcmc/auto_adaptation_test.cpp b/src/test/unit/mcmc/auto_adaptation_test.cpp index 6989016325..7c7234596d 100644 --- a/src/test/unit/mcmc/auto_adaptation_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_test.cpp @@ -28,8 +28,8 @@ TEST(McmcAutoAdaptation, test_covariance_one_row_one_col) { ASSERT_EQ(cov2.rows(), 1); ASSERT_EQ(cov2.cols(), 1); - - for(int i = 0; i < cov1.size(); ++i) { + + for (int i = 0; i < cov1.size(); ++i) { ASSERT_FLOAT_EQ(cov1(i), 0.0); } @@ -39,7 +39,7 @@ TEST(McmcAutoAdaptation, test_covariance_one_row_one_col) { TEST(McmcAutoAdaptation, test_covariance) { Eigen::MatrixXd X1(3, 2); Eigen::MatrixXd X2(2, 3); - + X1 << 0.0, -1.0, 0.5, -2.7, 3.0, 5.0; X2 << 0.0, 3, -2.7, 0.5, -1, 5.0; @@ -49,24 +49,22 @@ TEST(McmcAutoAdaptation, test_covariance) { Eigen::MatrixXd cov1_ref(2, 2); Eigen::MatrixXd cov2_ref(3, 3); - cov1_ref << 2.5833333333333335, 6.0666666666666664, - 6.0666666666666664, 16.3633333333333333; + cov1_ref << 2.5833333333333335, 6.0666666666666664, 6.0666666666666664, + 16.3633333333333333; - cov2_ref << 0.125, -1.0, 1.925, - -1.000, 8.0, -15.4, - 1.925, -15.4, 29.645; + cov2_ref << 0.125, -1.0, 1.925, -1.000, 8.0, -15.4, 1.925, -15.4, 29.645; ASSERT_EQ(cov1.rows(), cov1_ref.rows()); ASSERT_EQ(cov1.cols(), cov1_ref.cols()); ASSERT_EQ(cov2.rows(), cov2_ref.rows()); ASSERT_EQ(cov2.cols(), cov2_ref.cols()); - - for(int i = 0; i < cov1_ref.size(); ++i) { + + for (int i = 0; i < cov1_ref.size(); ++i) { ASSERT_FLOAT_EQ(cov1(i), cov1_ref(i)); } - for(int i = 0; i < cov2_ref.size(); ++i) { + for (int i = 0; i < cov2_ref.size(); ++i) { ASSERT_FLOAT_EQ(cov2(i), cov2_ref(i)); } } @@ -77,7 +75,7 @@ TEST(McmcAutoAdaptation, power_method) { X << 2.0, 0.5, 0.5, 1.0; x0 << 1.0, 0.0; - + const int max_iterations = 10; const double tol = 1e-10; @@ -85,9 +83,10 @@ TEST(McmcAutoAdaptation, power_method) { int max_iterations_1 = max_iterations; double tol_1 = tol; - - double eval = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); - + + double eval + = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); + EXPECT_FLOAT_EQ(eval, 2.20710678118654746); } @@ -105,8 +104,9 @@ TEST(McmcAutoAdaptation, power_method_tol_check) { int max_iterations_1 = max_iterations; double tol_1 = tol; - double eval = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); - + double eval + = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); + EXPECT_LT(tol_1, tol); } @@ -124,8 +124,9 @@ TEST(McmcAutoAdaptation, power_method_iter_check) { int max_iterations_1 = max_iterations; double tol_1 = tol; - double eval = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); - + double eval + = stan::mcmc::internal::power_method(Av, x0, max_iterations_1, tol_1); + EXPECT_GT(tol_1, tol); EXPECT_EQ(max_iterations_1, max_iterations); } @@ -157,14 +158,16 @@ TEST(McmcAutoAdaptation, eigenvalue_scaled_hessian) { data_stream.close(); std::stringstream output; - known_hessian_model_namespace::known_hessian_model known_hessian_model(data_var_context, &output); + known_hessian_model_namespace::known_hessian_model known_hessian_model( + data_var_context, &output); Eigen::MatrixXd L(3, 3); Eigen::VectorXd q(3); L << 2.0, 0.0, 0.0, 0.7, 1.3, 0.0, -1.5, 2.0, 4.0; q << 0.0, 0.0, 0.0; - double eval = stan::mcmc::internal::eigenvalue_scaled_hessian(known_hessian_model, L, q); + double eval = stan::mcmc::internal::eigenvalue_scaled_hessian( + known_hessian_model, L, q); EXPECT_LT(std::abs(eval - 22.8141075806892850) / eval, 1e-2); } From 9213b08f3c6ce66a86c16dd5b754c3de50e89a6e Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 29 Oct 2020 16:27:51 -0400 Subject: [PATCH 10/20] cpplint fixes (Issue #2814) --- src/stan/mcmc/auto_adaptation.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index 1863f17d6c..aba982a222 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -11,7 +11,7 @@ namespace stan { template struct log_prob_wrapper_covar { const Model& model_; - log_prob_wrapper_covar(const Model& model) : model_(model) {} + explicit log_prob_wrapper_covar(const Model& model) : model_(model) {} template T operator()(const Eigen::Matrix& q) const { @@ -180,7 +180,7 @@ namespace stan { if(state == "selection") { int Mtest; - Mtest = int(0.2 * Y.cols()); + Mtest = static_cast(0.2 * Y.cols()); if(Mtest < 5) { Mtest = 5; } From fa20a9f526a3d13d377cd69859f333800a690d04 Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 29 Oct 2020 16:34:11 -0400 Subject: [PATCH 11/20] Cpplint fix (Issue #2814) --- src/stan/mcmc/auto_adaptation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index cb51a54a08..0265ea227d 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -11,7 +11,7 @@ namespace mcmc { template struct log_prob_wrapper_covar { const Model& model_; - log_prob_wrapper_covar(const Model& model) : model_(model) {} + explicit log_prob_wrapper_covar(const Model& model) : model_(model) {} template T operator()(const Eigen::Matrix& q) const { From 0a2f629fdc0971ce8424e01869b3b0326c423023 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 30 Oct 2020 13:06:19 -0400 Subject: [PATCH 12/20] Updated tests (Issue #2814) --- .../mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp | 2 +- src/test/unit/mcmc/auto_adaptation_test.cpp | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp index b645e17496..72f1c7b444 100644 --- a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp @@ -10,7 +10,7 @@ TEST(McmcVarAdaptation, learn_covariance_pick_dense) { std::stringstream output; correlated_gaussian_model_namespace::correlated_gaussian_model - correlated_gaussian_model(data_var_context, &output); + correlated_gaussian_model(data_var_context, 0, &output); stan::test::unit::instrumented_logger logger; diff --git a/src/test/unit/mcmc/auto_adaptation_test.cpp b/src/test/unit/mcmc/auto_adaptation_test.cpp index 7c7234596d..301a99fa3d 100644 --- a/src/test/unit/mcmc/auto_adaptation_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_test.cpp @@ -158,8 +158,7 @@ TEST(McmcAutoAdaptation, eigenvalue_scaled_hessian) { data_stream.close(); std::stringstream output; - known_hessian_model_namespace::known_hessian_model known_hessian_model( - data_var_context, &output); + known_hessian_model_namespace::known_hessian_model known_hessian_model(data_var_context, 0, &output); Eigen::MatrixXd L(3, 3); Eigen::VectorXd q(3); From 90216026eed5100711c5bd94acee92492cda2d34 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 30 Oct 2020 17:09:44 +0000 Subject: [PATCH 13/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- .../mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp | 2 +- src/test/unit/mcmc/auto_adaptation_test.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp index 72f1c7b444..8b278d5395 100644 --- a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_dense_test.cpp @@ -10,7 +10,7 @@ TEST(McmcVarAdaptation, learn_covariance_pick_dense) { std::stringstream output; correlated_gaussian_model_namespace::correlated_gaussian_model - correlated_gaussian_model(data_var_context, 0, &output); + correlated_gaussian_model(data_var_context, 0, &output); stan::test::unit::instrumented_logger logger; diff --git a/src/test/unit/mcmc/auto_adaptation_test.cpp b/src/test/unit/mcmc/auto_adaptation_test.cpp index 301a99fa3d..a281b40cc6 100644 --- a/src/test/unit/mcmc/auto_adaptation_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_test.cpp @@ -158,7 +158,8 @@ TEST(McmcAutoAdaptation, eigenvalue_scaled_hessian) { data_stream.close(); std::stringstream output; - known_hessian_model_namespace::known_hessian_model known_hessian_model(data_var_context, 0, &output); + known_hessian_model_namespace::known_hessian_model known_hessian_model( + data_var_context, 0, &output); Eigen::MatrixXd L(3, 3); Eigen::VectorXd q(3); From 8cf629f11b219e84107b645b8d627a656e40c0f8 Mon Sep 17 00:00:00 2001 From: Ben Date: Sat, 31 Oct 2020 12:51:08 -0400 Subject: [PATCH 14/20] Switched strings for characters in state machine (Issue #2814) --- src/stan/mcmc/auto_adaptation.hpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index 0265ea227d..cae6ce42e3 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -192,11 +192,14 @@ class auto_adaptation : public windowed_adaptation { try { bool use_dense = false; - for (auto state : {"selection", "refinement"}) { + // 's' stands for selection + // 'r' stands for refinement + for (char state : { 's', 'r' }) { Eigen::MatrixXd Ytrain; Eigen::MatrixXd Ytest; - if (state == "selection") { + // If in selection state + if (state == 's') { int Mtest; Mtest = static_cast(0.2 * Y.cols()); if (Mtest < 5) { @@ -231,7 +234,8 @@ class auto_adaptation : public windowed_adaptation { covar = dense; - if (state == "selection") { + // If in selection state + if (state == 's') { Eigen::MatrixXd L_dense = dense.llt().matrixL(); Eigen::MatrixXd L_diag = diag.diagonal().array().sqrt().matrix().asDiagonal(); From ba03ec66c417cdd8c034d0474580dccc6183f989 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 31 Oct 2020 16:51:56 +0000 Subject: [PATCH 15/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- src/stan/mcmc/auto_adaptation.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index cae6ce42e3..98ff193058 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -192,13 +192,13 @@ class auto_adaptation : public windowed_adaptation { try { bool use_dense = false; - // 's' stands for selection - // 'r' stands for refinement - for (char state : { 's', 'r' }) { + // 's' stands for selection + // 'r' stands for refinement + for (char state : {'s', 'r'}) { Eigen::MatrixXd Ytrain; Eigen::MatrixXd Ytest; - // If in selection state + // If in selection state if (state == 's') { int Mtest; Mtest = static_cast(0.2 * Y.cols()); @@ -234,7 +234,7 @@ class auto_adaptation : public windowed_adaptation { covar = dense; - // If in selection state + // If in selection state if (state == 's') { Eigen::MatrixXd L_dense = dense.llt().matrixL(); Eigen::MatrixXd L_diag From b8707022b58b4985a52b7422aa64638d4f967923 Mon Sep 17 00:00:00 2001 From: Ben Date: Sat, 31 Oct 2020 13:49:02 -0400 Subject: [PATCH 16/20] Updated tolerances on power method (Issue #2814) --- src/stan/mcmc/auto_adaptation.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index cae6ce42e3..7b5fb33f1a 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -109,8 +109,8 @@ double eigenvalue_scaled_covariance(const Eigen::MatrixXd& L, auto Sx = [&](Eigen::VectorXd x) -> Eigen::VectorXd { return S * x; }; - int max_iterations = 100; - double tol = 1e-3; + int max_iterations = 200; + double tol = 1e-5; return internal::power_method(Sx, Eigen::VectorXd::Random(Sigma.cols()), max_iterations, tol); @@ -147,8 +147,8 @@ double eigenvalue_scaled_hessian(const Model& model, const Eigen::MatrixXd& L, return L.transpose() * (grad1 - grad2) / dx; }; - int max_iterations = 100; - double tol = 1e-3; + int max_iterations = 200; + double tol = 1e-5; return internal::power_method( hessian_vector, Eigen::VectorXd::Random(q.size()), max_iterations, tol); From e3d12212f1a6d31c830058e36daef66c814eeee7 Mon Sep 17 00:00:00 2001 From: Ben Date: Sat, 31 Oct 2020 14:37:31 -0400 Subject: [PATCH 17/20] Updated test to make randomness matter less (Issue #2814) --- ...tation_learn_covariance_pick_diag_test.cpp | 35 +++++++++++-------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp index 52aab40076..1c5f8bf02a 100644 --- a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp @@ -17,26 +17,33 @@ TEST(McmcVarAdaptation, learn_covariance_pick_diagonal) { const int M = 2; const int N = 20; Eigen::MatrixXd qs(N, M); - qs << 0.607446257145326, 0.338465765807058, 1.47389672467345, - -1.0577986841911, 1.02886652895522, 0.364277500948572, 0.492316893603469, - 2.19693408641558, -0.931854393410476, 1.62634580968769, - -0.443145375724188, 0.902790875582656, 0.517782110245233, - -1.56724331755861, -1.7556390097031, 0.310274990315213, - 0.0394975482340945, 0.366999438969482, 1.29372950054929, - 0.361369734821582, -0.258301497542829, 0.166994731172984, - 0.492639248874412, -0.659502589885556, 0.913729457222598, - 1.99580706461809, 0.669655370469707, -0.509028392475839, - -0.626041244059129, -0.771981104624195, -0.842385483586737, - 0.337166271031201, 0.548177804329155, -0.0462961925005498, - 0.955748803092952, 1.3141117316189, 0.335670079140694, 1.09112083087171, - 0.759245358940033, -1.11318882201676; + qs << -1.12310858, -0.98044486, + -0.40288484, -0.31777537, + -0.46665535, -0.41468940, + 0.77996512, 0.79419535, + -0.08336907, 0.12997631, + 0.25331851, 0.17888355, + -0.02854676, -0.25660897, + -0.04287046, 0.06199044, + 1.36860228, 1.16082198, + -0.22577099, -0.27199475, + 1.51647060, 1.46738068, + -1.54875280, -1.42235482, + 0.58461375, 0.40408060, + 0.12385424, 0.12959917, + 0.21594157, 0.18045828, + 0.37963948, 0.34225195, + -0.50232345, -0.41356307, + -0.33320738, -0.33695265, + -1.01857538, -0.85228019, + -1.07179123, -0.98666076; Eigen::MatrixXd covar(M, M); bool covar_is_diagonal; Eigen::MatrixXd target_covar(M, M); - target_covar << 0.55350038163333048, 0.0, 0.0, 0.86122545968912112; + target_covar << 0.50172809066980716963, 0.0, 0.0, 0.41270751419364998247; stan::mcmc::auto_adaptation adapter(M); adapter.set_window_params(50, 0, 0, N, logger); From b07569ab968bcf29703332930e9f437976d16a49 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 31 Oct 2020 19:37:51 +0000 Subject: [PATCH 18/20] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- ...tation_learn_covariance_pick_diag_test.cpp | 27 +++++-------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp index 1c5f8bf02a..64e3a946be 100644 --- a/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp +++ b/src/test/unit/mcmc/auto_adaptation_learn_covariance_pick_diag_test.cpp @@ -17,26 +17,13 @@ TEST(McmcVarAdaptation, learn_covariance_pick_diagonal) { const int M = 2; const int N = 20; Eigen::MatrixXd qs(N, M); - qs << -1.12310858, -0.98044486, - -0.40288484, -0.31777537, - -0.46665535, -0.41468940, - 0.77996512, 0.79419535, - -0.08336907, 0.12997631, - 0.25331851, 0.17888355, - -0.02854676, -0.25660897, - -0.04287046, 0.06199044, - 1.36860228, 1.16082198, - -0.22577099, -0.27199475, - 1.51647060, 1.46738068, - -1.54875280, -1.42235482, - 0.58461375, 0.40408060, - 0.12385424, 0.12959917, - 0.21594157, 0.18045828, - 0.37963948, 0.34225195, - -0.50232345, -0.41356307, - -0.33320738, -0.33695265, - -1.01857538, -0.85228019, - -1.07179123, -0.98666076; + qs << -1.12310858, -0.98044486, -0.40288484, -0.31777537, -0.46665535, + -0.41468940, 0.77996512, 0.79419535, -0.08336907, 0.12997631, 0.25331851, + 0.17888355, -0.02854676, -0.25660897, -0.04287046, 0.06199044, 1.36860228, + 1.16082198, -0.22577099, -0.27199475, 1.51647060, 1.46738068, -1.54875280, + -1.42235482, 0.58461375, 0.40408060, 0.12385424, 0.12959917, 0.21594157, + 0.18045828, 0.37963948, 0.34225195, -0.50232345, -0.41356307, -0.33320738, + -0.33695265, -1.01857538, -0.85228019, -1.07179123, -0.98666076; Eigen::MatrixXd covar(M, M); bool covar_is_diagonal; From 9a3cb36f2530651eb56942737e21c132a6ceefc4 Mon Sep 17 00:00:00 2001 From: Ben Date: Sat, 31 Oct 2020 17:02:56 -0400 Subject: [PATCH 19/20] Fixed header problems (Issue #2814) --- src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp | 3 +-- src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp | 12 ------------ src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp | 4 ++-- 3 files changed, 3 insertions(+), 16 deletions(-) diff --git a/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp index cff85d1dca..3041252a65 100644 --- a/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp +++ b/src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp @@ -2,12 +2,11 @@ #define STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_METRIC_HPP #include -#include +#include #include #include #include #include -#include namespace stan { namespace mcmc { diff --git a/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp b/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp index 13e24cc44f..3a0c2fa8e1 100644 --- a/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp +++ b/src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp @@ -33,18 +33,6 @@ class auto_e_point : public ps_point { inv_e_metric_.setIdentity(); } - /** - * Copy constructor which does fast copy of inverse mass matrix. - * - * @param z point to copy - */ - auto_e_point(const auto_e_point& z) - : ps_point(z), - inv_e_metric_(z.inv_e_metric_.rows(), z.inv_e_metric_.cols()) { - fast_matrix_copy_(inv_e_metric_, z.inv_e_metric_); - is_diagonal_ = z.is_diagonal_; - } - /** * Set elements of mass matrix * diff --git a/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp index addd110753..d66693426c 100644 --- a/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp +++ b/src/stan/services/sample/hmc_nuts_auto_e_adapt.hpp @@ -1,11 +1,11 @@ #ifndef STAN_SERVICES_SAMPLE_HMC_NUTS_AUTO_E_ADAPT_HPP #define STAN_SERVICES_SAMPLE_HMC_NUTS_AUTO_E_ADAPT_HPP -#include -#include +#include #include #include #include +#include #include #include #include From c805e39a90cb485944ed52e342bf920ac64638b9 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Wed, 3 Aug 2022 13:16:47 -0400 Subject: [PATCH 20/20] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- src/stan/mcmc/auto_adaptation.hpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/stan/mcmc/auto_adaptation.hpp b/src/stan/mcmc/auto_adaptation.hpp index f174f70026..8ff6df92ca 100644 --- a/src/stan/mcmc/auto_adaptation.hpp +++ b/src/stan/mcmc/auto_adaptation.hpp @@ -224,11 +224,10 @@ class auto_adaptation : public windowed_adaptation { = (Ytest.cols() > 0) ? internal::covariance(Ytest.transpose()) : Eigen::MatrixXd::Zero(M, M); - Eigen::MatrixXd dense - = (N / (N + 5.0)) * cov_train - + 1e-3 * (5.0 / (N + 5.0)) - * Eigen::MatrixXd::Identity(cov_train.rows(), - cov_train.cols()); + Eigen::MatrixXd dense = (N / (N + 5.0)) * cov_train + + 1e-3 * (5.0 / (N + 5.0)) + * Eigen::MatrixXd::Identity( + cov_train.rows(), cov_train.cols()); Eigen::MatrixXd diag = dense.diagonal().asDiagonal();