This is an old version, view current version.

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,
                              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, 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, 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, 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.

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, 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, real[] x_r, int[] x_i, real rel_tol, 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, type vector,

  • theta: parameters only, type vector,

  • x_r: real data only, type real[], and

  • x_i: integer data only, type int[].

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.

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, and others 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.

References

Guennebaud, Gaël, Benoît Jacob, and others. 2010. “Eigen V3.” http://eigen.tuxfamily.org.

Hindmarsh, Alan C, Peter N Brown, Keith E Grant, Steven L Lee, Radu Serban, Dan E Shumaker, and Carol S Woodward. 2005. “SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers.” ACM Transactions on Mathematical Software (TOMS) 31 (3): 363–96.

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.

Powell, Michael J. D. 1970. “A Hybrid Method for Nonlinear Equations.” In Numerical Methods for Nonlinear Algebraic Equations, edited by P. Rabinowitz. Gordon; Breach.