9.1 Algebraic Equation Solver
Stan provides a built-in algebraic equation solver. Although it looks like other function applications, the algebraic solver is special in two ways.
First, the algebraic solver is a higher-order function, i.e. it takes another function as one of its arguments. The only other functions in Stan which share this feature are the ordinary differential equation solvers (see section Ordinary Differential Equation 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, 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 for
theta, parameter values used to evaluate the algebraic system
x_r, data values used to evaluate the algebraic system
x_i, integer data used to evaluate the algebraic system
The algebraic system function separates parameter values,
from data values,
x_r, for efficiency in computing the gradients
of the algebraic system.
9.1.2 Call to the Algebraic Solver
(function algebra_system, vector y_guess, vector theta, real x_r, int x_i)
Solves the algebraic system, given an initial guess, using the Powell hybrid algorithm.
(function algebra_system, vector y_guess, vector theta, real x_r, int x_i, real rel_tol, 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.
22.214.171.124 Arguments to the Algebraic Solver
The arguments to the algebraic solver 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
- 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, type
theta: parameters only, type
x_r: real data only, type
x_i: integer data only, type
For more fine-grained control of the algebraic solver, these parameters can also be provided:
rel_tol: relative tolerance for the algebraic solver, type
real, data only,
function_tol: function tolerance for the algebraic solver, type
real, data only,
max_num_steps: maximum number of steps to take in the algebraic solver, type
int, data only.
126.96.36.199 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.
188.8.131.52 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.
184.108.40.206 Algorithmic details
The algebraic solver uses 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, and others 2010). This solver is in turn based on the algorithm developed for the package MINPACK-1 (Jorge J. More 1980).
The Jacobian of the solution with respect to auxiliary parameters is computed using the implicit function theorem. Intermediate Jacobians (of the 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.
Powell, Michael J. D. 1970. “A Hybrid Method for Nonlinear Equations.” In Numerical Methods for Nonlinear Algebraic Equations, edited by P. Rabinowitz. Gordon; Breach.
Guennebaud, Gaël, Benoît Jacob, and others. 2010. “Eigen V3.” http://eigen.tuxfamily.org.
Jorge J. More, Kenneth E. Hillstrom, Burton S. Garbow. 1980. User Guide for Minpack-1. 9700 South Cass Avenue, Argonne, Illinois 60439: Argonne National Laboratory.