1 #ifndef STAN_OPTIMIZATION_NEWTON_HPP 2 #define STAN_OPTIMIZATION_NEWTON_HPP 7 #include <Eigen/Cholesky> 8 #include <Eigen/Eigenvalues> 12 namespace optimization {
14 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
matrix_d;
15 typedef Eigen::Matrix<double, Eigen::Dynamic, 1>
vector_d;
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]);
29 g = eigenvectors * eigenprojections;
35 std::vector<double>& params_r,
36 std::vector<int>& params_i,
37 std::ostream* output_stream = 0) {
42 = stan::model::grad_hess_log_prob<true, false>(model,
45 matrix_d H(params_r.size(), params_r.size());
46 for (
size_t i = 0; i < hessian.size(); i++) {
49 vector_d g(params_r.size());
50 for (
size_t i = 0; i < gradient.size(); i++)
52 make_negative_definite_and_solve(H, g);
55 std::vector<double> new_params_r(params_r.size());
57 double min_step_size = 1e-50;
62 if (step_size < min_step_size)
65 for (
size_t i = 0; i < params_r.size(); i++)
66 new_params_r[i] = params_r[i] - step_size * g[i];
68 f1 = stan::model::log_prob_grad<true, false>(model,
71 }
catch (std::exception& e) {
76 for (
size_t i = 0; i < params_r.size(); i++)
77 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