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