9.1 Algebraic equation solver
Stan provides two built-in algebraic equation solvers, respectively based on Powell’s and Newton’s methods. The Newton method constitutes a more recent addition to Stan; its use is recommended for most problems. Although they look like other function applications, algebraic solvers are special in two ways.
First, an algebraic solver is a higher-order function, i.e. it takes another function as one of its arguments. Other functions in Stan which share this feature are the ordinary differential equation solvers (see section Ordinary Differential Equation (ODE) Solvers). Ordinary Stan functions do not allow functions as arguments.
Second, some of the arguments of the algebraic solvers are restricted to data only expressions. These expressions must not contain variables other than those declared in the data or transformed data blocks. Ordinary Stan functions place no restriction on the origin of variables in their argument expressions.
9.1.1 Specifying an algebraic equation as a function
An algebraic system is specified as an ordinary function in Stan within the function block. The algebraic system function must have this signature:
vector algebra_system(vector y, vector theta,
data real[] x_r, int[] x_i)The algebraic system function should return the value of the algebraic function which goes to 0, when we plug in the solution to the algebraic system.
The argument of this function are:
y, the unknowns we wish to solve fortheta, parameter values used to evaluate the algebraic systemx_r, data values used to evaluate the algebraic systemx_i, integer data used to evaluate the algebraic system
The algebraic system function separates parameter values, theta,
from data values, x_r, for efficiency in propagating the derivatives
through the algebraic system.
9.1.2 Call to the algebraic solver
vector algebra_solver(function algebra_system, vector y_guess, vector theta, data real[] x_r, int[] x_i)
Solves the algebraic system, given an initial guess, using the Powell
hybrid algorithm.
vector algebra_solver(function algebra_system, vector y_guess, vector theta, data real[] x_r, int[] x_i, data real rel_tol, data real f_tol, int max_steps)
Solves the algebraic system, given an initial guess, using the Powell
hybrid algorithm with additional control parameters for the solver.
Note: In future releases, the function algebra_solver will be deprecated
and replaced with algebra_solver_powell.
vector algebra_solver_newton(function algebra_system, vector y_guess, vector theta, data real[] x_r, int[] x_i)
Solves the algebraic system, given an initial guess, using Newton’s method.
vector algebra_solver_newton(function algebra_system, vector y_guess, vector theta, data real[] x_r, int[] x_i, data real rel_tol, data real f_tol, int max_steps)
Solves the algebraic system, given an initial guess, using Newton’s method
with additional control parameters for the solver.
9.1.2.1 Arguments to the algebraic solver
The arguments to the algebraic solvers are as follows:
algebra_system: function literal referring to a function specifying the system of algebraic equations with signature(vector, vector, real[], int[]):vector. The arguments represent (1) unknowns, (2) parameters, (3) real data, and (4) integer data, and the return value contains the value of the algebraic function, which goes to 0 when we plug in the solution to the algebraic system,y_guess: initial guess for the solution, typevector,theta: parameters only, typevector,x_r: real data only, typereal[], andx_i: integer data only, typeint[].
For more fine-grained control of the algebraic solver, these parameters can also be provided:
rel_tol: relative tolerance for the algebraic solver, typereal, data only,function_tol: function tolerance for the algebraic solver, typereal, data only,max_num_steps: maximum number of steps to take in the algebraic solver, typeint, data only.
9.1.2.2 Return value
The return value for the algebraic solver is an object of type
vector, with values which, when plugged in as y make the algebraic
function go to 0.
9.1.2.3 Sizes and parallel arrays
Certain sizes have to be consistent. The initial guess, return value of the solver, and return value of the algebraic function must all be the same size.
The parameters, real data, and integer data will be passed from the solver directly to the system function.
9.1.2.4 Algorithmic details
Stan offers two algebraic solvers: algebra_solver and algebra_solver_newton.
algebra_solver is baed on the Powell hybrid method (Powell 1970),
which in turn uses first-order derivatives. The Stan code builds on
the implementation of the hybrid solver in the unsupported module for
nonlinear optimization problems of the Eigen library (Guennebaud, Jacob et al. 2010).
This solver is in turn based on the algorithm developed for the
package MINPACK-1 (Jorge J. More 1980).
algebra_solver_newton, uses Newton’s method,
also a first-order derivative based numerical solver.
The Stan code builds on the implementation in KINSOL
from the SUNDIALS suite (Hindmarsh et al. 2005).
For many problems, we find that algebra_solver_newton is faster
than Powell’s method.
If however Newton’s method performs poorly, either failing to or requiring an excessively
long time to converge, the user should be prepared to switch
to algebra_solver.
For both solvers, the Jacobian of the solution with respect to auxiliary parameters is computed using the implicit function theorem. Intermediate Jacobians (of the algebraic function’s output with respect to the unknowns y and with respect to the auxiliary parameters theta) are computed using Stan’s automatic differentiation.