This is an old version, view current version.

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, theta, from data values, x_r, for efficiency in computing the gradients of 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.

9.1.2.1 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
  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

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.

References

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.