Stan  2.14.0
probability, sampling & optimization
softabs_metric.hpp
Go to the documentation of this file.
1 #ifndef STAN_MCMC_HMC_HAMILTONIANS_SOFTABS_METRIC_HPP
2 #define STAN_MCMC_HMC_HAMILTONIANS_SOFTABS_METRIC_HPP
3 
4 #include <stan/math/mix/mat.hpp>
7 #include <boost/random/variate_generator.hpp>
8 #include <boost/random/normal_distribution.hpp>
9 
10 namespace stan {
11  namespace mcmc {
12 
13  template <typename Model>
14  struct softabs_fun {
15  const Model& model_;
16  std::ostream* o_;
17 
18  softabs_fun(const Model& m, std::ostream* out): model_(m), o_(out) {}
19 
20  template <typename T>
21  T operator()(Eigen::Matrix<T, Eigen::Dynamic, 1>& x) const {
22  return model_.template log_prob<true, true, T>(x, o_);
23  }
24  };
25 
26  // Riemannian manifold with SoftAbs metric
27  template <class Model, class BaseRNG>
29  : public base_hamiltonian<Model, softabs_point, BaseRNG> {
30  private:
31  typedef typename stan::math::index_type<Eigen::VectorXd>::type idx_t;
32  public:
33  explicit softabs_metric(const Model& model)
34  : base_hamiltonian<Model, softabs_point, BaseRNG>(model) {}
35 
36  double T(softabs_point& z) {
37  return this->tau(z) + 0.5 * z.log_det_metric;
38  }
39 
40  double tau(softabs_point& z) {
41  Eigen::VectorXd Qp = z.eigen_deco.eigenvectors().transpose() * z.p;
42  return 0.5 * Qp.transpose() * z.softabs_lambda_inv.cwiseProduct(Qp);
43  }
44 
45  double phi(softabs_point& z) {
46  return this->V(z) + 0.5 * z.log_det_metric;
47  }
48 
49  double dG_dt(softabs_point& z,
52  return 2 * T(z)
53  - z.q.dot(dtau_dq(z, info_writer, error_writer)
54  + dphi_dq(z, info_writer, error_writer));
55  }
56 
57  Eigen::VectorXd dtau_dq(
58  softabs_point& z,
61  Eigen::VectorXd a = z.softabs_lambda_inv.cwiseProduct(
62  z.eigen_deco.eigenvectors().transpose() * z.p);
63  Eigen::MatrixXd A = a.asDiagonal()
64  * z.eigen_deco.eigenvectors().transpose();
65  Eigen::MatrixXd B = z.pseudo_j.selfadjointView<Eigen::Lower>() * A;
66  Eigen::MatrixXd C = A.transpose() * B;
67 
68  Eigen::VectorXd b(z.q.size());
70  softabs_fun<Model>(this->model_, 0), z.q, C, b);
71 
72  return 0.5 * b;
73  }
74 
75  Eigen::VectorXd dtau_dp(softabs_point& z) {
76  return z.eigen_deco.eigenvectors()
77  * z.softabs_lambda_inv.cwiseProduct(
78  z.eigen_deco.eigenvectors().transpose() * z.p);
79  }
80 
81  Eigen::VectorXd dphi_dq(
82  softabs_point& z,
85  Eigen::VectorXd a
86  = z.softabs_lambda_inv.cwiseProduct(z.pseudo_j.diagonal());
87  Eigen::MatrixXd A = a.asDiagonal()
88  * z.eigen_deco.eigenvectors().transpose();
89  Eigen::MatrixXd B = z.eigen_deco.eigenvectors() * A;
90 
92  softabs_fun<Model>(this->model_, 0), z.q, B, a);
93 
94  return - 0.5 * a + z.g;
95  }
96 
97  void sample_p(softabs_point& z, BaseRNG& rng) {
98  boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
99  rand_unit_gaus(rng, boost::normal_distribution<>());
100 
101  Eigen::VectorXd a(z.p.size());
102 
103  for (idx_t n = 0; n < z.p.size(); ++n)
104  a(n) = sqrt(z.softabs_lambda(n)) * rand_unit_gaus();
105 
106  z.p = z.eigen_deco.eigenvectors() * a;
107  }
108 
109  void init(
110  softabs_point& z,
113  update_metric(z, info_writer, error_writer);
114  update_metric_gradient(z, info_writer, error_writer);
115  }
116 
118  softabs_point& z,
121  // Compute the Hessian
122  stan::math::hessian<double>(
123  softabs_fun<Model>(this->model_, 0), z.q, z.V, z.g, z.hessian);
124 
125  z.V = -z.V;
126  z.g = -z.g;
127  z.hessian = -z.hessian;
128 
129  // Compute the eigen decomposition of the Hessian,
130  // then perform the SoftAbs transformation
131  z.eigen_deco.compute(z.hessian);
132 
133  for (idx_t i = 0; i < z.q.size(); ++i) {
134  double lambda = z.eigen_deco.eigenvalues()(i);
135  double alpha_lambda = z.alpha * lambda;
136 
137  double softabs_lambda = 0;
138 
139  // Thresholds defined such that the approximation
140  // error is on the same order of double precision
141  if (std::fabs(alpha_lambda) < lower_softabs_thresh) {
142  softabs_lambda = (1.0
143  + (1.0 / 3.0) * alpha_lambda * alpha_lambda)
144  / z.alpha;
145  } else if (std::fabs(alpha_lambda) > upper_softabs_thresh) {
146  softabs_lambda = std::fabs(lambda);
147  } else {
148  softabs_lambda = lambda / std::tanh(alpha_lambda);
149  }
150 
151  z.softabs_lambda(i) = softabs_lambda;
152  z.softabs_lambda_inv(i) = 1.0 / softabs_lambda;
153  }
154 
155  // Compute the log determinant of the metric
156  z.log_det_metric = 0;
157  for (idx_t i = 0; i < z.q.size(); ++i)
158  z.log_det_metric += std::log(z.softabs_lambda(i));
159  }
160 
162  softabs_point& z,
165  // Compute the pseudo-Jacobian of the SoftAbs transform
166  for (idx_t i = 0; i < z.q.size(); ++i) {
167  for (idx_t j = 0; j <= i; ++j) {
168  double delta = z.eigen_deco.eigenvalues()(i)
169  - z.eigen_deco.eigenvalues()(j);
170 
171  if (std::fabs(delta) < jacobian_thresh) {
172  double lambda = z.eigen_deco.eigenvalues()(i);
173  double alpha_lambda = z.alpha * lambda;
174 
175  // Thresholds defined such that the approximation
176  // error is on the same order of double precision
177  if (std::fabs(alpha_lambda) < lower_softabs_thresh) {
178  z.pseudo_j(i, j) = (2.0 / 3.0) * alpha_lambda
179  * (1.0 - (2.0 / 15.0)
180  * alpha_lambda * alpha_lambda);
181  } else if (std::fabs(alpha_lambda) > upper_softabs_thresh) {
182  z.pseudo_j(i, j) = lambda > 0 ? 1 : -1;
183  } else {
184  double sdx = std::sinh(alpha_lambda) / lambda;
185  z.pseudo_j(i, j) = (z.softabs_lambda(i)
186  - z.alpha / (sdx * sdx) ) / lambda;
187  }
188  } else {
189  z.pseudo_j(i, j) = (z.softabs_lambda(i)
190  - z.softabs_lambda(j) ) / delta;
191  }
192  }
193  }
194  }
195 
197  softabs_point& z,
200  update_metric_gradient(z, info_writer, error_writer);
201  }
202 
203  // Threshold below which a power series
204  // approximation of the softabs function is used
205  static double lower_softabs_thresh;
206 
207  // Threshold above which an asymptotic
208  // approximation of the softabs function is used
209  static double upper_softabs_thresh;
210 
211  // Threshold below which an exact derivative is
212  // used in the Jacobian calculation instead of
213  // finite differencing
214  static double jacobian_thresh;
215  };
216 
217  template <class Model, class BaseRNG>
219 
220  template <class Model, class BaseRNG>
222 
223  template <class Model, class BaseRNG>
225  } // mcmc
226 } // stan
227 #endif
Eigen::VectorXd dphi_dq(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Probability, optimization and sampling library.
Eigen::VectorXd dtau_dq(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Eigen::VectorXd softabs_lambda_inv
Eigen::VectorXd q
Definition: ps_point.hpp:45
Eigen::VectorXd softabs_lambda
void init(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
double T(softabs_point &z)
Eigen::VectorXd p
Definition: ps_point.hpp:46
double dG_dt(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
void update_metric(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Eigen::SelfAdjointEigenSolver< Eigen::MatrixXd > eigen_deco
double tau(softabs_point &z)
Point in a phase space with a base Riemannian manifold with SoftAbs metric.
base_writer is an abstract base class defining the interface for Stan writer callbacks.
Definition: base_writer.hpp:20
T operator()(Eigen::Matrix< T, Eigen::Dynamic, 1 > &x) const
void update_gradients(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
void update_metric_gradient(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Eigen::VectorXd g
Definition: ps_point.hpp:49
void grad_tr_mat_times_hessian(const M &model, const Eigen::Matrix< double, Eigen::Dynamic, 1 > &x, const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > &X, Eigen::Matrix< double, Eigen::Dynamic, 1 > &grad_tr_X_hess_f, std::ostream *msgs=0)
Eigen::VectorXd dtau_dp(softabs_point &z)
softabs_fun(const Model &m, std::ostream *out)
double phi(softabs_point &z)
softabs_metric(const Model &model)
void sample_p(softabs_point &z, BaseRNG &rng)

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