what spamfunc is

In developing optimization algorithms, one of the most tedious parts is trying different examples, each of which might have its own starting points or upper or lower bounds or other information. The tedium really starts when your algorithm requires first or second order information, which might be tricky to calculate correctly. These bugs can be pernicious, because it might be difficult to differentiate between a bug in your algorithm and a bug in your objective or constraint evaluation. Handily, Northwestern Professor Robert Fourer wrote a language called AMPL, which takes a programming-language specification of objective and constraints and calculates derivatives as needed. The official amplfunc/spamfunc reference is contained in Hooking Your Solver to AMPL, but I’m shooting for a more low-key introduction.

You’ll need to compile spamfunc before using it, but I’m saving that for last–until after I’ve shown off it’s power a little bit. You’ll also need a .nl file for spamfunc; these are generated by AMPL’s command line program from an AMPL .mod source file. I’ll be using rosenb.nl as my example here, for the Rosenbrock function.

>> [x0,bl,bu,lam,cl,cu]=spamfunc('rosenb.nl');

This serves as the initialization. It’s worth noting the mathematical form AMPL bends your problem into. In particular, AMPL’s spamfunc returns your problem in the form ( \begin{split} \min \qquad & f(x)\ \text{s.t.}\qquad & c_L \leq c(x) \leq c_U\ & x_L \leq x \leq x_U. \end{split} )

So we’ve initialized 6 variables for our environment: the initial (x_0) we specified in the mod file, lower($x_L$) and upper bounds on those variables, the initial Lagrange multipler(lam, short for $\lambda$), and the lower($c_L$) and upper($c_U$) bounds on the constraints. To get zeroth order information(that is, $f(x)$ and $c(x)$) for a particular point $x$, just call spamfunc again with the the $x$ you’d like to evaluate at, as in

>> [f,c] = spamfunc(x,0);

To get first order information(the gradient of the objective and Jacobian of constraints), call spamfunc as

>> [g,A]=spamfunc(x,1);

so that the gradient $\nabla f(x)$ is stored in g, and the $i^{th}$ row of A is the transposed-gradient $(\nabla c_i(x))^T$. That is, the entry $A_{ij}$ of the matrix $A$ is the derivative $\frac{\partial c_i(x)}{\partial x_j}$.

Finally, to get the Hessian of the Lagrangian, call spamfunc with just the multipliers $\lambda$, as in

>> WL = spamfunc(lam);
>> WF = spamfunc(0*lam); % if you just want \nabla^2 f(x).

Here, spamfunc reuses the last $x$ at which spamfunc(x,1) was evaluated at.

how to use spamfunc in your own code

Fine and good, but how do you use this in your code? I’ll consider a very simple unconstrained optimization problem(and finally define the Rosenbrock function, if you’ve been waiting with bated breath). The Rosenbrock function–or, Rosenbrock’s banana function, as I’ll be calling it from now on, is defined as ( f(x_1, x_2) = 100(x_2 - x_1^2)^2 + (1-x_1)^2 ) If we wanted to use our spamfunc interface to solve this, we’d generate the .nl file as described below, then do a standard Newton method iteration. Since Newton steps are defined as $p_k = -\nabla^2 f(x_k) \nabla f(x_k)$ and then the iterate is updated as $x_{k+1}=x_k + p_k$, the most basic possible form of the algorithm is

[x,bl,bu,lam,cl,cu]=spamfunc('rosenb.nl');

grad = inf;
while norm(grad) > 1e-6
	    [grad,A] = spamfunc(x,1);
	    W = spamfunc(lam);
	    pk = -W\grad;
	    x = x + pk
end

Note: This is a terrible algorithm; there’re no bells and whistles like line search, Hessian modification, etc. But it demonstrates calling at least a couple of spamfunc options.

how to compile spamfunc

I recovered an outline for building spamfunc and amplfunc from a project I’ve worked on; it should get you a compiled spamfunc binary. You need to make sure to have mex (MATLAB’s compiler wrapper) in your path and properly configured.

wget --ftp-user anonymous --ftp-password traviscj@traviscj.com ftp://netlib.org/ampl.tar
# could do the .gz version, but only saves 5mb
tar xf ampl.tar
cd ampl/solvers
sh configurehere
mv makefile makefile.dist
# possibly fixed, needed for old osx, or linux?
sed "s/CFLAGS = -O/CFLAGS = -fPIC -O/g" makefile.dist > Makefile
make
cd examples
# following line should reflect where mex(1) is.
export PATH=/meas/bin:$PATH
mv makefile makefile.dist
# first sed cleans up evalchk, second one maybe needed for linux??
sed "s/- o/-o/g" makefile.dist | sed "s/mex -I/mex -ldl -I/g" > Makefile
make spamfunc.mex amplfunc.mex

It’s worth noting the difference between amplfunc and spamfunc; amplfunc returns dense matrices while spamfunc returns them in sparse format. MATLAB can do linear algebra on sparse format matrices much more efficiently than dense format matrices(as long as the matrices themselves are sparse enough!)

how to generate nl files

In the above examples, I started with an AMPL source file, rosen.mod.

// rosenb.mod
var x {1..2};

minimize obj: 100*(x[2] - x[1]^2)^2 + (1-x[1])^2;

let x[1] := 1.2;
let x[2] := 1.2;

This file tells AMPL what how many variables you need, how you’d like to name them, what objective function you’re minimizing, and what initial values they should take. To compile it to a .nl file, use

$ ampl -P -ogrosenb rosenb.mod

Congratulations! You now have a .nl file to feed into your super sweet solver!