Stan  2.14.0
probability, sampling & optimization
newton.hpp
Go to the documentation of this file.
1 #ifndef STAN_OPTIMIZATION_NEWTON_HPP
2 #define STAN_OPTIMIZATION_NEWTON_HPP
3 
6 #include <Eigen/Dense>
7 #include <Eigen/Cholesky>
8 #include <Eigen/Eigenvalues>
9 #include <vector>
10 
11 namespace stan {
12  namespace optimization {
13 
14  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_d;
15  typedef Eigen::Matrix<double, Eigen::Dynamic, 1> vector_d;
16 
17  namespace {
18  // Negates any positive eigenvalues in H so that H is negative
19  // definite, and then solves Hu = g and stores the result into
20  // g. Avoids problems due to non-log-concave distributions.
21  void make_negative_definite_and_solve(matrix_d& H, vector_d& g) {
22  Eigen::SelfAdjointEigenSolver<matrix_d> solver(H);
23  matrix_d eigenvectors = solver.eigenvectors();
24  vector_d eigenvalues = solver.eigenvalues();
25  vector_d eigenprojections = eigenvectors.transpose() * g;
26  for (int i = 0; i < g.size(); i++) {
27  eigenprojections[i] = -eigenprojections[i] / fabs(eigenvalues[i]);
28  }
29  g = eigenvectors * eigenprojections;
30  }
31  }
32 
33  template <typename M>
34  double newton_step(M& model,
35  std::vector<double>& params_r,
36  std::vector<int>& params_i,
37  std::ostream* output_stream = 0) {
38  std::vector<double> gradient;
39  std::vector<double> hessian;
40 
41  double f0
42  = stan::model::grad_hess_log_prob<true, false>(model,
43  params_r, params_i,
44  gradient, hessian);
45  matrix_d H(params_r.size(), params_r.size());
46  for (size_t i = 0; i < hessian.size(); i++) {
47  H(i) = hessian[i];
48  }
49  vector_d g(params_r.size());
50  for (size_t i = 0; i < gradient.size(); i++)
51  g(i) = gradient[i];
52  make_negative_definite_and_solve(H, g);
53 // H.ldlt().solveInPlace(g);
54 
55  std::vector<double> new_params_r(params_r.size());
56  double step_size = 2;
57  double min_step_size = 1e-50;
58  double f1 = -1e100;
59 
60  while (f1 < f0) {
61  step_size *= 0.5;
62  if (step_size < min_step_size)
63  return f0;
64 
65  for (size_t i = 0; i < params_r.size(); i++)
66  new_params_r[i] = params_r[i] - step_size * g[i];
67  try {
68  f1 = stan::model::log_prob_grad<true, false>(model,
69  new_params_r,
70  params_i, gradient);
71  } catch (std::exception& e) {
72  // FIXME: this is not a good way to handle a general exception
73  f1 = -1e100;
74  }
75  }
76  for (size_t i = 0; i < params_r.size(); i++)
77  params_r[i] = new_params_r[i];
78 
79  return f1;
80  }
81 
82  }
83 }
84 #endif
Probability, optimization and sampling library.
void hessian(const M &model, const Eigen::Matrix< double, Eigen::Dynamic, 1 > &x, double &f, Eigen::Matrix< double, Eigen::Dynamic, 1 > &grad_f, Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > &hess_f, std::ostream *msgs=0)
Definition: hessian.hpp:12
double newton_step(M &model, std::vector< double > &params_r, std::vector< int > &params_i, std::ostream *output_stream=0)
Definition: newton.hpp:34
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_d
Definition: newton.hpp:14
void gradient(const M &model, const Eigen::Matrix< double, Eigen::Dynamic, 1 > &x, double &f, Eigen::Matrix< double, Eigen::Dynamic, 1 > &grad_f, std::ostream *msgs=0)
Definition: gradient.hpp:14
Eigen::Matrix< double, Eigen::Dynamic, 1 > vector_d
Definition: newton.hpp:15

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