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  // 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  inline 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  template <typename M>
32  double newton_step(M& model,
33  std::vector<double>& params_r,
34  std::vector<int>& params_i,
35  std::ostream* output_stream = 0) {
36  std::vector<double> gradient;
37  std::vector<double> hessian;
38 
39  double f0
40  = stan::model::grad_hess_log_prob<true, false>(model,
41  params_r, params_i,
42  gradient, hessian);
43  matrix_d H(params_r.size(), params_r.size());
44  for (size_t i = 0; i < hessian.size(); i++) {
45  H(i) = hessian[i];
46  }
47  vector_d g(params_r.size());
48  for (size_t i = 0; i < gradient.size(); i++)
49  g(i) = gradient[i];
51 // H.ldlt().solveInPlace(g);
52 
53  std::vector<double> new_params_r(params_r.size());
54  double step_size = 2;
55  double min_step_size = 1e-50;
56  double f1 = -1e100;
57 
58  while (f1 < f0) {
59  step_size *= 0.5;
60  if (step_size < min_step_size)
61  return f0;
62 
63  for (size_t i = 0; i < params_r.size(); i++)
64  new_params_r[i] = params_r[i] - step_size * g[i];
65  try {
66  f1 = stan::model::log_prob_grad<true, false>(model,
67  new_params_r,
68  params_i, gradient);
69  } catch (std::exception& e) {
70  // FIXME: this is not a good way to handle a general exception
71  f1 = -1e100;
72  }
73  }
74  for (size_t i = 0; i < params_r.size(); i++)
75  params_r[i] = new_params_r[i];
76 
77  return f1;
78  }
79 
80  }
81 }
82 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void make_negative_definite_and_solve(matrix_d &H, vector_d &g)
Definition: newton.hpp:20
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:32
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_d
Definition: newton.hpp:14
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
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:15
Float_t f1
Eigen::Matrix< double, Eigen::Dynamic, 1 > vector_d
Definition: newton.hpp:15
Float_t e
Definition: plot.C:35
const XML_Char XML_Content * model
Definition: expat.h:151