bfgs.hpp
Go to the documentation of this file.
1 #ifndef STAN_SERVICES_OPTIMIZE_BFGS_HPP
2 #define STAN_SERVICES_OPTIMIZE_BFGS_HPP
3 
14 #include <fstream>
15 #include <iostream>
16 #include <iomanip>
17 #include <string>
18 #include <vector>
19 
20 namespace stan {
21  namespace services {
22  namespace optimize {
23 
24  /**
25  * Runs the BFGS algorithm for a model.
26  *
27  * @tparam Model A model implementation
28  * @param[in] model Input model to test (with data already instantiated)
29  * @param[in] init var context for initialization
30  * @param[in] random_seed random seed for the random number generator
31  * @param[in] chain chain id to advance the pseudo random number generator
32  * @param[in] init_radius radius to initialize
33  * @param[in] init_alpha line search step size for first iteration
34  * @param[in] tol_obj convergence tolerance on absolute changes in
35  * objective function value
36  * @param[in] tol_rel_obj convergence tolerance on relative changes
37  * in objective function value
38  * @param[in] tol_grad convergence tolerance on the norm of the gradient
39  * @param[in] tol_rel_grad convergence tolerance on the relative norm of
40  * the gradient
41  * @param[in] tol_param convergence tolerance on changes in parameter
42  * value
43  * @param[in] num_iterations maximum number of iterations
44  * @param[in] save_iterations indicates whether all the interations should
45  * be saved to the parameter_writer
46  * @param[in] refresh how often to write output to logger
47  * @param[in,out] interrupt callback to be called every iteration
48  * @param[in,out] logger Logger for messages
49  * @param[in,out] init_writer Writer callback for unconstrained inits
50  * @param[in,out] parameter_writer output for parameter values
51  * @return error_codes::OK if successful
52  */
53  template <class Model>
55  unsigned int random_seed, unsigned int chain,
56  double init_radius, double init_alpha, double tol_obj,
57  double tol_rel_obj, double tol_grad, double tol_rel_grad,
58  double tol_param, int num_iterations, bool save_iterations,
59  int refresh,
60  callbacks::interrupt& interrupt,
61  callbacks::logger& logger,
62  callbacks::writer& init_writer,
63  callbacks::writer& parameter_writer) {
64  boost::ecuyer1988 rng = util::create_rng(random_seed, chain);
65 
66  std::vector<int> disc_vector;
67  std::vector<double> cont_vector
68  = util::initialize<false>(model, init, rng, init_radius, false,
69  logger, init_writer);
70 
71  std::stringstream bfgs_ss;
74  Optimizer bfgs(model, cont_vector, disc_vector, &bfgs_ss);
75  bfgs._ls_opts.alpha0 = init_alpha;
76  bfgs._conv_opts.tolAbsF = tol_obj;
77  bfgs._conv_opts.tolRelF = tol_rel_obj;
78  bfgs._conv_opts.tolAbsGrad = tol_grad;
79  bfgs._conv_opts.tolRelGrad = tol_rel_grad;
80  bfgs._conv_opts.tolAbsX = tol_param;
81  bfgs._conv_opts.maxIts = num_iterations;
82 
83  double lp = bfgs.logp();
84 
85  std::stringstream initial_msg;
86  initial_msg << "Initial log joint probability = " << lp;
87  logger.info(initial_msg);
88 
89  std::vector<std::string> names;
90  names.push_back("lp__");
91  model.constrained_param_names(names, true, true);
92  parameter_writer(names);
93 
94  if (save_iterations) {
95  std::vector<double> values;
96  std::stringstream msg;
97  model.write_array(rng, cont_vector, disc_vector, values,
98  true, true, &msg);
99  if (msg.str().length() > 0)
100  logger.info(msg);
101 
102  values.insert(values.begin(), lp);
103  parameter_writer(values);
104  }
105  int ret = 0;
106 
107  while (ret == 0) {
108  interrupt();
109  if (refresh > 0
110  && (bfgs.iter_num() == 0
111  || ((bfgs.iter_num() + 1) % refresh == 0)))
112  logger.info(" Iter"
113  " log prob"
114  " ||dx||"
115  " ||grad||"
116  " alpha"
117  " alpha0"
118  " # evals"
119  " Notes ");
120 
121  ret = bfgs.step();
122  lp = bfgs.logp();
123  bfgs.params_r(cont_vector);
124 
125  if (refresh > 0
126  && (ret != 0
127  || !bfgs.note().empty()
128  || bfgs.iter_num() == 0
129  || ((bfgs.iter_num() + 1) % refresh == 0))) {
130  std::stringstream msg;
131  msg << " " << std::setw(7) << bfgs.iter_num() << " ";
132  msg << " " << std::setw(12) << std::setprecision(6)
133  << lp << " ";
134  msg << " " << std::setw(12) << std::setprecision(6)
135  << bfgs.prev_step_size() << " ";
136  msg << " " << std::setw(12) << std::setprecision(6)
137  << bfgs.curr_g().norm() << " ";
138  msg << " " << std::setw(10) << std::setprecision(4)
139  << bfgs.alpha() << " ";
140  msg << " " << std::setw(10) << std::setprecision(4)
141  << bfgs.alpha0() << " ";
142  msg << " " << std::setw(7)
143  << bfgs.grad_evals() << " ";
144  msg << " " << bfgs.note() << " ";
145  logger.info(msg);
146  }
147 
148  if (bfgs_ss.str().length() > 0) {
149  logger.info(bfgs_ss);
150  bfgs_ss.str("");
151  }
152 
153  if (save_iterations) {
154  std::vector<double> values;
155  std::stringstream msg;
156  model.write_array(rng, cont_vector, disc_vector, values,
157  true, true, &msg);
158  // This if is here to match the pre-refactor behavior
159  if (msg.str().length() > 0)
160  logger.info(msg);
161 
162  values.insert(values.begin(), lp);
163  parameter_writer(values);
164  }
165  }
166 
167  if (!save_iterations) {
168  std::vector<double> values;
169  std::stringstream msg;
170  model.write_array(rng, cont_vector, disc_vector, values,
171  true, true, &msg);
172  if (msg.str().length() > 0)
173  logger.info(msg);
174  values.insert(values.begin(), lp);
175  parameter_writer(values);
176  }
177 
178  int return_code;
179  if (ret >= 0) {
180  logger.info("Optimization terminated normally: ");
181  return_code = error_codes::OK;
182  } else {
183  logger.info("Optimization terminated with error: ");
184  return_code = error_codes::SOFTWARE;
185  }
186  logger.info(" " + bfgs.get_code_string(ret));
187 
188  return return_code;
189  }
190 
191  }
192  }
193 }
194 #endif
Scalar prev_step_size() const
Definition: bfgs.hpp:100
const VectorT & curr_g() const
Definition: bfgs.hpp:93
stan::optimization::BFGSMinimizer< stan::optimization::ModelAdaptor< Model >, stan::optimization::BFGSUpdate_HInv<> > Optimizer
const Scalar & alpha() const
Definition: bfgs.hpp:112
const Scalar & alpha0() const
Definition: bfgs.hpp:111
rosenbrock_model_namespace::rosenbrock_model Model
std::string get_code_string(int retCode)
Definition: bfgs.hpp:117
chain
Check that an output directory exists.
LSOptions< Scalar > _ls_opts
Definition: bfgs.hpp:85
ConvergenceOptions< Scalar > _conv_opts
Definition: bfgs.hpp:86
const std::string & note() const
Definition: bfgs.hpp:115
const size_t iter_num() const
Definition: bfgs.hpp:113
int bfgs(Model &model, stan::io::var_context &init, unsigned int random_seed, unsigned int chain, double init_radius, double init_alpha, double tol_obj, double tol_rel_obj, double tol_grad, double tol_rel_grad, double tol_param, int num_iterations, bool save_iterations, int refresh, callbacks::interrupt &interrupt, callbacks::logger &logger, callbacks::writer &init_writer, callbacks::writer &parameter_writer)
Definition: bfgs.hpp:54
const XML_Char XML_Content * model
Definition: expat.h:151
boost::ecuyer1988 create_rng(unsigned int seed, unsigned int chain)
Definition: create_rng.hpp:25
void refresh()
Definition: show_event.C:21