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

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