help intsolver

Title

intsolver -- Mata Functions and Structures for finding all solutions to nonlinear systems of equations

Introduction

intsolver is a collection of Mata routines for completely solving - finding all zeros - of nonlinear systems. intsolver uses
interval methods for solving systems of equations. Interval methods for solving systems are fairly accessible to the novice, and
are also numerically stable. Descriptions appear in Moore et. al. (1999), and Jaulin et al. (2001).

The specific method intsolver uses to completely solve systems is a form of constraint
propogation combined with a so-called Taylor inclusion for the function
of interest. The method is a simplified version of some of the methods described in Hentenryck et. al. (1997), and also bears resemblance to
the equation-by-equation interval Gauss-Seidel method described in Moore et. al. (1999).

An extensive description of interval arithmetic is described in the package int_utils, which,
along with the package rowmat_utils, must be installed for intsolver
to work (along with moremata).

For a brief overview of the method, consider the problem of determining whether or not a
zero of a function f(x) exists in some range [a,b]. Interval arithmetic uses outward rounding
to create a range of possible values for f(x) on [a,b], which in trade jargon is referred to as
an interval extension of f(x): sometimes written as [f]([x]).

In most cases, the interval extension [f]([x]) has two important properties. First, the bounds it generates
are not sharp in that the true mininum (maximum) value of f(x) over the range [x] is greater (less than) the bounds
provided by [f]([x]). This phenomenon is referred to as pessimism. Second, it can be shown that as the size
of the range [x] shrinks, it is the case that the bounds [f]([x])
imposes on the domain of f(x) shrink more rapidly towards the true domain of f(x) on [a,b].
This convergence typically occurs at a rate >2 for a Taylor inclusion function, which is the type of inclusion
function used by intsolver.

This latter feature of the problem forms the basis for algorithms in which intervals are continuously narrowed
via bisection or some other means. One maintains
a list of potential intervals where a zero might
be located, and then tests the intervals to see if [f]([x]) brackets zero.
If a tested interval contains zero, it can be shrunk by some method (perhaps bisection) into smaller intervals, which can
then be retested. If a tested interval does not bracket zero, it is discarded. This method can
be extended to multiple variables, and can also use different
sorts of information about the problem. The process terminates when all intervals remaining are of sufficiently small size.
These remaining intervals can be resolved into point solutions
via some technique such as Newton Iteration.

intsolver uses simple interval tests which are performed in parallel over sets of intervals,
Consider f(x) to be a vector-valued function of the vector x. A first-order Taylor expansion about
the point x0 can be written as f(x)~f(x0)+J(x0)(x-x0). Suppose [x] is an interval vector, with the first row containing
lower bounds and the second upper bounds (see int_utils for examples).
In interval arithmetic, it can be shown that f([x]) is
bounded by f(x0)+[J]([x])([x]-x0), where [J]([x]) is the interval Jacobian.
This principal is simple and owes to the mean value theorem: the value of the function cannot be greater (smaller) than its greatest (smallest)
value multiplied by the greatest (smallest) value of the derivative. The algorithm used by intsolver is a version of interval Gauss-Seidel
iteration (see, e.g., Moore et. al. 1999, p. 96).

There are other methods, but in applications it has been found by the author
that using this interval function equation by equation and variable by variable
is a relatively efficient means of solving most problems.

Usage

Syntax

Syntax is discussed in three steps:

Step 1: Preliminaries
Step 2: Initialization
Step 3: Solving a problem
Step 4: Refining and obtaining results

Step 1: Preliminaries

Ideas are best motivated by example. Accordingly, consider the problem of finding all the
zeros of the vector function f={1-2x1x2-3x1^2, f2=1-x1^2-3x2^2}
where both x1 and y1 are restricted to the [-100,100] interval. While simple, this
function reflects one of the situations for which intsolver was intended: Smaller
nonlinear problems that are difficult to solve completely by hand, that have multiple real solutions.
This function in fact has four roots in [-100,100].

To use intsolver, one must program four Mata functions: one that returns the value of a particular equation at a
set of points, one that returns an entry of the jacobian at a set of points, an interval representation of values, and
an interval representation of jacobian entries. The actual function arguments should be passed to the function
in a single matrix. Each row of the matrix is a sequence of points
at which the value of one of the component functions can be calculated. The first
step is to program the two functions themselves. First, the function itself:

. mata:

mata (type end to exit)

: real matrix fun(real matrix x,real scalar i)
>{
>if (i==1) return(1:-2*x[,1]:*x[,2]:-3*x[,1]:^2)
>if (i==2) return(1:-x[,1]:^2:-3*x[,2]:^2)
>}
>: end

The function fun() takes in two arguments. The first, a matrix, holds (x1,x2) combinations
listed as its rows. The second argument is simply the equation number in the system.

The next step is programming the function providing values of the entries of the Jacobian.
For the example, f11=-2x2-6x1, f12=-2x1, f21=-2x1, f22=-6x2. So:

. mata:

mata (type end to exit)

: real matrix jac(real matrix x,real scalar i, real scalar j)
>{
>if (i==1 & j==1) return(-2*x[,2]:-6:*x[,1])
>if (i==1 & j==2) return(-2*x[,1])
>if (i==2 & j==1) return(-2*x[,1]])
>if (i==2 & j==2) return(-6*x[,2])
>}
: end

Next, one must program an interval version of fun(). This is more challenging,
but can be accomplished through application of int_utils. In
the author's experience, it is useful to break interval operations into steps. Some operations,
like interval-scalar multiplication, can be
treated on the fly, although one should exert care to avoid rounding errors, which are
crucial to interval arithmetic. Treating a constant like an interval with equal upper and
lower limits is helpful in avoiding these sorts of errors.
The arguments of the function once again are rows of an input matrix, but now the first
two columns are the lower and upper bound of the first variable, columns 3 and 4 are the bounds of the
second variable, etc. Accordingly, the interval functions should be set up to take in a
matrix as the first argument, a scalar indicating the equation number as the next argument, and a
final argument (d in the functions below). This last argument is a tolerance for
rounding that is used in interval computations. It is described further below and also the help entry
for int_utils.

. mata:

mata (type end to exit)

: real matrix fun_I(real matrix x,real scalar i, real scalar d)
>{
>X1=x[,1::2]
>X2=x[,3::4]
>if (i==1) {
>	A=2*int_mult(X1,X2,d)
>	I=J(rows(A),2,1)
>	B=int_sub(I,A,d)
>	C=3*int_mult(X1,X1,d)
>	D=int_sub(B,C,d)
>	return(D)
>}
>if (i==2) {
>	A=int_mult(X1,X1,d)
>	B=3*int_mult(X2,X2)
>	I=J(rows(A),2,1)
>	C=int_sub(I,A,d)
>	D=int_sub(C,B,d)
>	return(D)
>}
>}
: end

Some tricks of the trade that allow programming interval versions of functions are employed
in the above Mata program. The line I=J(rows(A),2,1) produces a two-column matrix of ones
that conforms with the arguments X1 and X2. The next line, int_sub(I,A,d) is
thus a way of performing interval subtraction, where one wishes to
subtract an interval from a constant.

The interval jacobian is programmed similarly:

. mata:

mata (type end to exit)

: real matrix jac_I(real matrix x,real scalar i,real scalar j, real scalar d)
>{
>X1=x[,1::2]
>X2=x[,3::4]
>r=rows(X1)
>if (i==1 & j==1) {
>	A=int_mult(J(r,2,-2),X2,d)
>	B=int_mult(J(r,2,-6),X1,d)
>	C=int_add(A,B,d)
>	return(C)
>}
>if (i==1 & j==2) {
>	A=int_mult(J(r,2,-2),X1,d)
>	return(A)
>}
>if (i==2 & j==1) {
>	A=int_mult(J(r,2,-2),X1,d)
>	return(A)
>if (i==2 & j==2) {
>	A=J(rows(X1),2,-6)
>	B=int_mult(A,X2,d)
>	return(B)
>}
>}
: end

Step 2: Initiation

Now that the basic problem has been defined, the user must describe several basic features of the problem.
These include the interval within which intsolver is to search
for solutions and various sorts of tolerances. One starts by defining a problem using the mata
function int_prob_init(), then filling in the number of arguments using int_prob_args()
and then passing the four required functions to the solver
using int_prob_f(), int_prob_jac(), int_prob_f_Iform(), and int_prob_jac_Iform(). Continuing the example:

. mata:

mata (type end to exit)

: Prob=int_prob_init()
: int_prob_args(Prob,2)
: int_prob_f_Iform(Prob,&fun_I())
: int_prob_jac_Iform(Prob,&jac_I())
: int_prob_f(Prob,&fun())
: int_prob_jac(Prob,&jac())
: end

Step 3: Solving a problem

The next step is define the intervals within which to search for solutions and apply the solver:

. mata:

mata (type end to exit)

: Ival=(-100,100) \ (-100,100)
: int_prob_ival(Prob,Ival)
: int_solve(Prob)
: end

The problem has been solved in interval form. One can view the solutions through use of
the utility int_prob_ints_vals(Prob). Entering the command produces the following:

. mata:

mata (type end to exit)

: int_prob_ints_vals(Prob)
1 2 3 4
+
+
1 | -.7249408 -.7249155 .3976849 .3976926 |
2 | -.7249408 -.7249155 .3976849 .3976926 |
3 | -.4291268 -.4291181 -.5214909 -.5214856 |
4 | .4291181 .4291268 .5214856 .5214909 |
5 | .7249155 .7249408 -.3977003 -.3976926 |
6 | .7249155 .7249408 -.3976926 -.3976849 |
+
+

Step 4: Refining and obtaining results

In this particular case, intsolver produces 6 intervals which potentially contain solutions. The user
might wish to refine these solutions, or get a more accurate count of the number of solutions to the system.
To this end, one can refine the intervals into points using int_newton_iter(),
which operates on the principle that interval solutions
define a neighborhood of convergence and hence newton iteration will converge rapidly, a principle for which
Pandian (1986) provides rigorous foundation. The
exact method used for newton iteration is discussed in the help entry rowmat_utils. Suffice it to say
here that intsolver uses rowmat_utils to perform Newton iteration in parallel. To illustrate:

: int_newton_iter(Prob)
: int_prob_pts_vals(Prob)
1 2
+
+
1 | -.7249360956 .3976881766 |
2 | -.4291212844 -.5214898284 |
3 | .4291212844 .5214898284 |
4 | .7249360956 -.3976881766 |
+
+

Thus, applying int_newton_iter() pares down the solution set and renders
the set of intervals into a set of points, so the equations have four solutions in the given range. In the accompanying .do file
IntsolverExamples.do, an example using Powell's badly-scaled function is presented. This example produces approximately 115 separate
intervals where solutions might be. Newton iteration pares this set down to two actual solutions.

Additional function arguments

In the event that one wishes to use a function that has additional arguments or parameters,
one should specify these in terms of a single argument (by, for example, collecting the additional arguments
in a pointer, matrix, vector, etc. and then unpacking the additional arguments
inside the body of the function). If it were the case, for example, that fun() had additional arguments,
instead of programming

. mata:

mata (type end to exit)

: real scalar fun(real matrix x, real scalar i)
>{

<body of function>

>}
>: end

One would instead collect arguments in some vector, matrix, pointer, or some other object -
Z, say - and code a function of the form:

. mata:

mata (type end to exit)

: real scalar fun(real matrix x, real scalar i, transmorphic Z)
>{

<body of function>

>}
>: end

Similarly, fun_I(), jac(), and jac_I() should be programmed with an additional argument Z:

. mata:

mata (type end to exit)

: real scalar fun_I(real matrix x, real scalar i, real scalar d, transmorphic Z)
>{

<body of function>

>}
>: end

. mata:

mata (type end to exit)

: real scalar jac(real matrix x, real scalar i, real scalar j, transmorphic Z)
>{

<body of function>

>}
>: end

. mata:

mata (type end to exit)

: real scalar jac_I(real matrix x, real scalar i, real scalar j, real scalar d, transmorphic Z)
>{

<body of function>

>}
>: end

The final thing that needs to be done if additional arguments are to be
passed to the function is intsolver must be informed of the additional information
needed to calculate the functions. This is done by using the function int_prob_addinfo().

As an example, consider Diamond's (1982) model of search and unemployment.
Diamond shows that the equations describing the model have
multiple solutions, and therefore in analysis of the model, one might like to
be able to find all of the solutions.

In the model, agents search for production opportunities. When an opportunity arises,
the costs of executing it is drawn from a distribution G(c), and the agent must
decide whether or not to do the project. Agents engaged in projects are freed from
projects at some rate. The measure of agents is normalized to one, and opportunities
arrive at the rate a. Moreover, opportunities end at the rate b(e), where e
is the employment rate and b'>0. Steady-state employment occurs when the flow out of
opportunities equals the flow into opportunities:

a(1-e*)G(c*)-e*b(e*)=0

In the above, c* is the threshold cost below which production opportunities are accepted,
so G(c*) is the probability a given opportunity is accepted. Diamond shows
that the cutoff cost c* is described by an equation which equates the flow value of the
current opportunity with cost c* with the flow value of searching for another opportunity
given projects with costs greater than c* will not be undertaken. This equation can be written:

c*(r+b(e)+aG(c*))=b(e)y+a*E[c|c<c*]

In the previous equation r is the discount factor, and y is the flow output generated
by a project. To parameterize the model, we will suppose that b(e)=be, where b is a constant, and
that G(c) is the pareto distribution with parameter alpha. Hence, G(c)=1-c^alpha.
Therefore, the model of interest has five parameters: a,alpha,b,r, and y.

The first step is to define the equilibrium conditions and the jacobian. The equilibrium conditions are:

. mata:

mata (type end to exit)

: real matrix eqconds(real matrix X, real scalar i, transmorphic Z)
{
>real scalar a,alpha,b,r
>real matrix eq
>a=Z[1]
>alpha=Z[2]
>b=Z[3]
>r=Z[4]
>y=Z[5]
>
>e=X[,1]
>c=X[,2]
>
>Gc=1:-c:^alpha
>if (i==1) {
>	eq=a*(1:-e):*Gc:-e:*b*e /* Equilibrium employment */
>}
>if (i==2) {
>	CE=(-alpha/(alpha+1)):*c:^(alpha+1)
>	eq=c*r:+b:*e*r:+a*c:*Gc-b:*e*y:-a:*CE	/* Equilibrium project adoption */
>}
>return(eq)
>}
: end

While the Jacobian is:

. mata:

mata (type end to exit)

:
: real matrix eqjac(real matrix X, real scalar i, real scalar j, transmorphic Z)
>{
>real scalar a,alpha,b,r
>real matrix eq
>a=Z[1]
>alpha=Z[2]
>b=Z[3]
>r=Z[4]
>y=Z[5]
>
>e=X[,1]
>c=X[,2]
>
>if (i==1 & j==1) {
>	Gc=1:-c:^alpha
>	eq=-a*Gc:-2*e:*b
>}
>if (i==1 & j==2) {
>	dGc=-alpha*c:^(alpha-1)
>	eq=a*(1:-e):*dGc
>}
>if (i==2 & j==1) {
>	eq=J(rows(X),1,b*r:-b*y)
>}
>if (i==2 & j==2) {
>	Gc=1:-c:^alpha
>	dGc=-alpha*c:^(alpha-1)
>	dCE=-alpha:*c:^alpha
>	eq=r:+a*Gc:+a*c:*dGc:-a:*dCE
>}
>return(eq)
>}
: end

Note how the functions now take an additional argument, Z, which contains the
five parameters of the model. Next, the interval version of the equilibrium conditions should
be programmed:

. mata:

mata (type end to exit)

: real matrix eqconds_I(real matrix X, real scalar i, real scalar digits, transmorphic Z)
>{
>real scalar a,alpha,b,r
>real matrix eq,Gc,A,B,C,D,CE,Alphas
>
>a=Z[1]
>alpha=Z[2]
>b=Z[3]
>r=Z[4]
>y=Z[5]
>
>e=X[,1::2]
>c=X[,3::4]
>
>Ones=J(rows(X),2,1)
>A=int_pow(c,alpha,digits)
>Gc=int_sub(Ones,A,digits)
>if (i==1) {
>	A=a*int_sub(Ones,e,digits)
>	B=b*int_mult(e,e,digits)
>	C=int_mult(A,Gc,digits)
>	eq=int_sub(C,B)
>}
>else if (i==2) {
>	Alphas=J(rows(X),2,-alpha/(alpha+1))
>	CE=int_mult(Alphas,int_pow(c,alpha+1,digits),digits)
>	A=r*int_add(c,b*e,digits)
>	B=a*int_mult(c,Gc,digits)
>	C=int_add(b*y*e,a*CE,digits)
>	D=int_sub(B,C,digits)
>	eq=int_add(A,D,digits)
>}
>return(eq)
>}
: end

The interval version of the jacobian must also be programmed:

. mata:

mata (type end to exit)

: real matrix eqjac_I(real matrix X, real scalar i, real scalar j, real scalar digits, transmorphic Z)
>{
>real scalar a,alpha,b,r
>real matrix eq
>
>a=Z[1]
>alpha=Z[2]
>b=Z[3]
>r=Z[4]
>y=Z[5]
>
>e=X[,1::2]
>c=X[,3::4]
>
>Ones=J(rows(X),2,1)
>A=int_pow(c,alpha,digits)
>Gc=int_sub(Ones,A,digits)
>if (i==1 & j==1) {
>	eq=int_sub(-a*Gc,2*e:*b,digits)
>}
>if (i==1 & j==2) {
>	dGc=-alpha*int_pow(c,alpha-1,digits)
>	eq=a*int_mult(int_sub(Ones,e),dGc,digits)
>}
>if (i==2 & j==1) {
>	eq=J(rows(X),2,b*r:-b*y)
>}
>if (i==2 & j==2) {
>	dGc=-alpha*int_pow(c,alpha-1,digits)
>	dCE=-alpha*int_pow(c,alpha,digits)
>	eq=r:+a*int_sub(int_add(Gc,int_mult(c,dGc,digits),digits),dCE,digits)
>}
>return(eq)
>}
: end

Next, the problem is initialized and solved over a reasonable range,
bounded away from zero:

. mata:

mata (type end to exit)

: Diamond=int_prob_init()
:int_prob_args(Diamond,2)
:int_prob_f(Diamond,&eqconds())
:int_prob_jac(Diamond,&eqjac())
:int_prob_f_Iform(Diamond,&eqconds_I())
:int_prob_jac_Iform(Diamond,&eqjac_I())
:bounds=0.000001,1 \ 0.000001,10
:int_prob_ival(Diamond,bounds)
:Z=2,-.12,3,.4,2 /* Parameter values */
:int_prob_addinfo(Diamond,Z)
:int_solve(Diamond)
:int_prob_ints_vals(Diamond)
1 2 3 4
+
+
1 | .027320859 .027321329 1.009640757 1.009643518 |
2 | .027320859 .027321329 1.009640757 1.009643518 |
3 | .24685418 .24685852 2.939442678 2.939458896 |
4 | .24685418 .24685852 2.939458895 2.939475113 |
+
+

:int_newton_iter(Diamond)
:int_prob_pts_vals(Diamond)
1 2
+
+
1 | .0273211682 1.009644396 |
2 | .2468549652 2.939448548 |
+
+

:end

One can see that the model, as parameterized by Z, has two equilibria.
One is a low-employment equililbrium in which only very low-cost opportunities are undertaken, while
the other is a high-employment equilibrium, in which higher-cost opportunities are undertaken.

Additional Notes

If one has a problem that does not necessarily require that all solutions be found,
one can just program the function and its jacobian.
Then, if a set (however large) of starting points, along with an initial range, are
provided to intsolver newton iteration can be performed in parallel on
the given points in the hopes that there will be some that are sufficiently close
to the solutions that they will be found. In the author's experience,
a little experimentation usually gives a good indication about how many starting
points are needed, and most if not all solutions can be found this way quite rapidly.

For example, suppose in the first example given above that only the functions
and its jacobian had been programmed.
Then, a "scatter shot" approach might be used, where random points around a solution might be passed to intsolver using int_prob_init_pts():

. mata:

mata (type end to exit)

: Prob=int_prob_init()
: int_prob_args(Prob,2)
: int_prob_f(Prob,&fun())
: int_prob_jac(Prob,&jac())
: randpoints=-1:+2*runiform(10,2)
: int_prob_init_pts(Prob,randpoints)
: int_prob_method(Prob2,"newton")
: int_newton_iter(Prob)
: int_prob_pts_vals(Prob)

This produces the same results as the complete solution method described above. int_newton_iter() does not
return as solutions as many points as it is given to start;
it automatically detects and lumps together answers that it deems to be the same.
How different points must be to be considered numerically different
can be controlled by the user. In applications, the user can calibrate the number of random points so that
seemingly all solutions are found.

Options

There are several things that the user may specify to control the solution process
that may greatly influence the speed at which
solutions are found. In terms of interval arithmetic, two important tolerances
are the tolerance used for outward rounding
in interval arithmetic, and the minimum width of a box. In terms of newton
iterations, the important numerical ideas are the convergence
tolerance and the maximum number of iterations to undergo. These four things may be set with the options:

int_prob_mbwidth( real scalar w)

Can be used to set the minimum box width. The default is 1e-4. Note that the box width
cannot be smaller than the number of digits used for outward rounding; otherwise
intsolver rounding down and up will always create
intervals bigger than the smallest box. Once all intervals are are reduced to
the specified size or smaller, intsolver stops bisection.
In practice, one can often get a sequence of intervals that are wide rather
quickly and then use the midpoints of these intervals in Newton iteration.

int_prob_digits( real scalar d)

Can be used to set the digits for outward rounding. The default is 1e-8. As an example,
when presented with an interval [.11223344556677,.11223344557788], if {bf:d=1e-8, the interval is rounded to
[.1122334,.1122335]. Interval arithmetic uses this "outward rounding" to avoid dropping potential solutions via rounding.

int_prob_maxit(real scalar its)

Can be used to set the maximum number of iterations performed by int_newton_iter(). The default
is 50.

int_prob_tol(real scalar tol

Can be used to set convergence criterion for the Newton iteration process. The default is 1e-12. This can
in practice be tuned much lower.

int_prob_tolsols(real scalar tolsols)

Can be used to lump solutions together. When newton iteration
produces multiple solutions, one would like to know if results
are in fact "the same" or not. If solutions lie within tolsols of one another,
then are deemed the same and one is dropped. The default
is 1e-2.

Test Problems and Examples

The supporting -.do file intsolver_examples.do contains an assortment of examples and test problems.

Further requirements

intsolver requires that the user install Ben Jann's help moremata set of commands, and also that the user install the
mata packages rowmat_utils and int_utils.

Author

Matthew J. Baker, Hunter College and the Graduate Center, CUNY, matthew.baker@hunter.cuny.edu.

Note from the author: Comments and suggestions are greatly appreciated, as usual. But I will also add that
if you are interested in using intsolver and need help setting
up your problem, please do not hesitate to contact me. I am hoping to build a larger library of examples, and also would like to know what sorts
of applications users have for the package.

References

Diamond, Peter. 1982. Aggregate demand management in search equilibrium. Journal of Political Economy 90(5): 881-94.

van Hentenryck, Pascal, David McAllester, and Deepak Kapur. 1997. Solving polynomial systems using a branch and prune approach. SIAM Journal on	Numerical Analysis 34(2): 797-827.

Jaulin, Luc, Michel Kieffer, Olivier Didrit, and Eric Walter. 2001. Applied Interval Analysis. London, Berlin, and Heidelberg: Spring-Verlag.

Moore, Ramon E., R. Baker Kearfott, and Michael J. Cloud. 2009. Introduction to Interval Analysis. Philadelpha: SIAM.

Pandian, Maharaja C. 1985. Convergence test and componentwise error estimates for Newton type methods. SIAM Journal on Numerical Analysis 22(4): 779-91.

Also see

rowmat_utils, int_utils, moremata (all are required for intsolver)