Stan Math Library  2.20.0
reverse mode automatic differentiation
algebra_solver.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_ALGEBRA_SOLVER_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_ALGEBRA_SOLVER_HPP
3 
4 #include <stan/math/rev/meta.hpp>
7 #include <stan/math/rev/core.hpp>
13 #include <unsupported/Eigen/NonLinearOptimization>
14 #include <iostream>
15 #include <string>
16 #include <vector>
17 
18 namespace stan {
19 namespace math {
20 
27 template <typename Fs, typename F, typename T, typename Fx>
28 struct algebra_solver_vari : public vari {
30  vari** y_;
32  int y_size_;
34  int x_size_;
38  double* Jx_y_;
39 
40  algebra_solver_vari(const Fs& fs, const F& f, const Eigen::VectorXd& x,
41  const Eigen::Matrix<T, Eigen::Dynamic, 1>& y,
42  const std::vector<double>& dat,
43  const std::vector<int>& dat_int,
44  const Eigen::VectorXd& theta_dbl, Fx& fx,
45  std::ostream* msgs)
46  : vari(theta_dbl(0)),
47  y_(ChainableStack::instance_->memalloc_.alloc_array<vari*>(y.size())),
48  y_size_(y.size()),
49  x_size_(x.size()),
50  theta_(
51  ChainableStack::instance_->memalloc_.alloc_array<vari*>(x_size_)),
52  Jx_y_(ChainableStack::instance_->memalloc_.alloc_array<double>(
53  x_size_ * y_size_)) {
54  using Eigen::Map;
55  using Eigen::MatrixXd;
56  for (int i = 0; i < y.size(); ++i)
57  y_[i] = y(i).vi_;
58 
59  theta_[0] = this;
60  for (int i = 1; i < x.size(); ++i)
61  theta_[i] = new vari(theta_dbl(i), false);
62 
63  // Compute the Jacobian and store in array, using the
64  // implicit function theorem, i.e. Jx_y = Jf_y / Jf_x
66  Map<MatrixXd>(&Jx_y_[0], x_size_, y_size_)
67  = -mdivide_left(fx.get_jacobian(theta_dbl),
68  f_y(fs, f, theta_dbl, value_of(y), dat, dat_int, msgs)
69  .get_jacobian(value_of(y)));
70  }
71 
72  void chain() {
73  for (int j = 0; j < y_size_; j++)
74  for (int i = 0; i < x_size_; i++)
75  y_[j]->adj_ += theta_[i]->adj_ * Jx_y_[j * x_size_ + i];
76  }
77 };
78 
123 template <typename F, typename T>
124 Eigen::VectorXd algebra_solver(
125  const F& f, const Eigen::Matrix<T, Eigen::Dynamic, 1>& x,
126  const Eigen::VectorXd& y, const std::vector<double>& dat,
127  const std::vector<int>& dat_int, std::ostream* msgs = nullptr,
128  double relative_tolerance = 1e-10, double function_tolerance = 1e-6,
129  long int max_num_steps = 1e+3) { // NOLINT(runtime/int)
130  check_nonzero_size("algebra_solver", "initial guess", x);
131  for (int i = 0; i < x.size(); i++)
132  check_finite("algebra_solver", "initial guess", x(i));
133  for (int i = 0; i < y.size(); i++)
134  check_finite("algebra_solver", "parameter vector", y(i));
135  for (double i : dat)
136  check_finite("algebra_solver", "continuous data", i);
137  for (int x : dat_int)
138  check_finite("algebra_solver", "integer data", x);
139 
140  if (relative_tolerance < 0)
141  invalid_argument("algebra_solver", "relative_tolerance,",
142  relative_tolerance, "",
143  ", must be greater than or equal to 0");
144  if (function_tolerance < 0)
145  invalid_argument("algebra_solver", "function_tolerance,",
146  function_tolerance, "",
147  ", must be greater than or equal to 0");
148  if (max_num_steps <= 0)
149  invalid_argument("algebra_solver", "max_num_steps,", max_num_steps, "",
150  ", must be greater than 0");
151 
152  // Create functor for algebraic system
155  Fx fx(Fs(), f, value_of(x), y, dat, dat_int, msgs);
156  Eigen::HybridNonLinearSolver<Fx> solver(fx);
157 
158  // Check dimension unknowns equals dimension of system output
159  check_matching_sizes("algebra_solver", "the algebraic system's output",
160  fx.get_value(value_of(x)), "the vector of unknowns, x,",
161  x);
162 
163  // Compute theta_dbl
164  Eigen::VectorXd theta_dbl = value_of(x);
165  solver.parameters.xtol = relative_tolerance;
166  solver.parameters.maxfev = max_num_steps;
167  solver.solve(theta_dbl);
168 
169  // Check if the max number of steps has been exceeded
170  if (solver.nfev >= max_num_steps) {
171  std::ostringstream message;
172  message << "algebra_solver: max number of iterations: " << max_num_steps
173  << " exceeded.";
174  throw boost::math::evaluation_error(message.str());
175  }
176 
177  // Check solution is a root
178  double system_norm = fx.get_value(theta_dbl).stableNorm();
179  if (system_norm > function_tolerance) {
180  std::ostringstream message2;
181  message2 << "algebra_solver: the norm of the algebraic function is: "
182  << system_norm << " but should be lower than the function "
183  << "tolerance: " << function_tolerance << ". Consider "
184  << "decreasing the relative tolerance and increasing the "
185  << "max_num_steps.";
186  throw boost::math::evaluation_error(message2.str());
187  }
188 
189  return theta_dbl;
190 }
191 
235 template <typename F, typename T1, typename T2>
236 Eigen::Matrix<T2, Eigen::Dynamic, 1> algebra_solver(
237  const F& f, const Eigen::Matrix<T1, Eigen::Dynamic, 1>& x,
238  const Eigen::Matrix<T2, Eigen::Dynamic, 1>& y,
239  const std::vector<double>& dat, const std::vector<int>& dat_int,
240  std::ostream* msgs = nullptr, double relative_tolerance = 1e-10,
241  double function_tolerance = 1e-6,
242  long int max_num_steps = 1e+3) { // NOLINT(runtime/int)
243  Eigen::VectorXd theta_dbl
244  = algebra_solver(f, x, value_of(y), dat, dat_int, 0, relative_tolerance,
245  function_tolerance, max_num_steps);
246 
248 
249  // TODO(charlesm93): a similar object gets constructed inside
250  // the call to algebra_solver. Cache the previous result
251  // and use it here (if possible).
254  Fx fx(Fs(), f, value_of(x), value_of(y), dat, dat_int, msgs);
255 
256  // Construct vari
258  = new algebra_solver_vari<Fy, F, T2, Fx>(Fy(), f, value_of(x), y, dat,
259  dat_int, theta_dbl, fx, msgs);
260  Eigen::Matrix<var, Eigen::Dynamic, 1> theta(x.size());
261  theta(0) = var(vi0->theta_[0]);
262  for (int i = 1; i < x.size(); ++i)
263  theta(i) = var(vi0->theta_[i]);
264 
265  return theta;
266 }
267 } // namespace math
268 } // namespace stan
269 
270 #endif
Eigen::Matrix< fvar< T >, R1, C2 > mdivide_left(const Eigen::Matrix< fvar< T >, R1, C1 > &A, const Eigen::Matrix< fvar< T >, R2, C2 > &b)
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
void check_nonzero_size(const char *function, const char *name, const T_y &y)
Check if the specified matrix/vector is of non-zero size.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
A functor that allows us to treat either x or y as the independent variable.
Eigen::VectorXd algebra_solver(const F &f, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &x, const Eigen::VectorXd &y, const std::vector< double > &dat, const std::vector< int > &dat_int, std::ostream *msgs=nullptr, double relative_tolerance=1e-10, double function_tolerance=1e-6, long int max_num_steps=1e+3)
Return the solution to the specified system of algebraic equations given an initial guess...
double * Jx_y_
Jacobian of the solution w.r.t parameters.
The variable implementation base class.
Definition: vari.hpp:30
algebra_solver_vari(const Fs &fs, const F &f, const Eigen::VectorXd &x, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &y, const std::vector< double > &dat, const std::vector< int > &dat_int, const Eigen::VectorXd &theta_dbl, Fx &fx, std::ostream *msgs)
int x_size_
number of unknowns
friend class var
Definition: vari.hpp:32
vari ** y_
vector of parameters
A functor with the required operators to call Eigen&#39;s algebraic solver.
The vari class for the algebraic solver.
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw an invalid_argument exception with a consistently formatted message.
void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
vari ** theta_
vector of solution
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
vari(double x)
Construct a variable implementation from a value.
Definition: vari.hpp:58
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
void check_matching_sizes(const char *function, const char *name1, const T_y1 &y1, const char *name2, const T_y2 &y2)
Check if two structures at the same size.
int y_size_
number of parameters
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
This struct always provides access to the autodiff stack using the singleton pattern.

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