1 #ifndef STAN_OPTIMIZATION_NEWTON_HPP
2 #define STAN_OPTIMIZATION_NEWTON_HPP
6 #include <Eigen/Cholesky>
7 #include <Eigen/Eigenvalues>
11 namespace optimization {
13 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
matrix_d;
14 typedef Eigen::Matrix<double, Eigen::Dynamic, 1>
vector_d;
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]);
28 g = eigenvectors * eigenprojections;
34 std::vector<double>& params_r,
35 std::vector<int>& params_i,
36 std::ostream* output_stream = 0) {
41 = stan::model::grad_hess_log_prob<true, false>(model,
44 matrix_d H(params_r.size(), params_r.size());
45 for (
size_t i = 0; i < hessian.size(); i++) {
48 vector_d g(params_r.size());
49 for (
size_t i = 0; i < gradient.size(); i++)
51 make_negative_definite_and_solve(H, g);
54 std::vector<double> new_params_r(params_r.size());
56 double min_step_size = 1e-50;
61 if (step_size < min_step_size)
64 for (
size_t i = 0; i < params_r.size(); i++)
65 new_params_r[i] = params_r[i] - step_size * g[i];
67 f1 = stan::model::log_prob_grad<true, false>(model,
70 }
catch (std::exception& e) {
75 for (
size_t i = 0; i < params_r.size(); i++)
76 params_r[i] = new_params_r[i];
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)
double newton_step(M &model, std::vector< double > ¶ms_r, std::vector< int > ¶ms_i, std::ostream *output_stream=0)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_d
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)
Eigen::Matrix< double, Eigen::Dynamic, 1 > vector_d