\documentclass[11pt,twoside]{article} \usepackage{amsmath} \usepackage{amssymb} \usepackage{subeqn,url} \textheight 8.75in \textwidth 5.5in \topmargin 0.0in \oddsidemargin 0.5in \evensidemargin 0.5in \pagestyle{myheadings} \markboth{MS\&E 318 (CME 338) \quad Large-Scale Numerical Optimization} {Spring 2009 Project} \title {$ $ \\[-1.25in] {\large Stanford University, Dept of Management Science and Engineering} \\[3pt] {\Large MS\&E 318 (CME 338) \quad Large-Scale Numerical Optimization} } \author{Instructor: Michael Saunders \qquad Spring 2009} \date {{\bf Project} \quad Due Thursday June 11} \input{macros} \begin{document} \maketitle \begin{center} \Large\bf GAMS and Basis Pursuit \end{center} \thispagestyle{empty} GAMS \cite{GAMS} is one of the widely used algebraic modeling languages for formulating optimization problems of many types. Our class on large-scale numerical optimization would be incomplete without giving you some experience with running a GAMS model. Luckily GAMS is available on the campus Linux servers, and we want to keep it that way by actually using it(!). We could say the same about AMPL \cite{AMPL}. The project will also help you learn about Basis Pursuit (BP) and Basis Pursuit Denoising (BPDN) \cite{CDS01}, which helped popularize the concept of finding sparse solutions to under-determined systems (also known as \emph{$\ell^1$-minimization} and \emph{sparse recovery}).% \footnote{In turn, BP has led to \emph{compressed sensing} \cite{CanRT06, Don06} and such concepts as the \emph{exact reconstruction principle} (ERP), the \emph{compressible reconstruction principle} (CRP), and the \emph{uniform uncertainty principle} (UUP) \cite{tao-sparse.html}.} Given a signal $b \in R^m$, signal analysis seeks a linear combination of certain waveforms that reconstruct $b$. The waveforms may be sinusoids (the traditional Fourier representation) or other ``dictionaries'' such as wavelets, chirplets, warplets, curvelets, etc. Many of the dictionaries are already ``overcomplete'' in the sense of containing more than $m$ waveforms. Combining one or more dictionaries gives a matrix $A \in R^{m\times n}$ with $m < n$. Basis Pursuit solves the $\ell^1$ problem \begin{align} \label{eqn:BP} \min_x \ \norm{x}_1 \text{ \ subject to \ } Ax = b \tag*{(BP)} \end{align} to obtain a solution $x$ that is likely to be \emph{sparse}.% \footnote{If the objective were $\min \, \norm{x}_p$ with $p \rightarrow 0^+$, the solution would be as sparse as possible. However, the problem would no longer be convex and its solution would be \emph{much} harder to find; in fact NP-hard. The choice $p=1$ gives the closest convex approximation to the $\ell^0$ problem.} Problem (BP) is equivalent to the linear program \begin{align} \label{eqn:BPLP} \min_{v,\,w} \ e^T(v+w) \text{ \ subject to \ } Av - Aw = b, \ \; v,w \ge 0, \tag*{(BPLP)} \end{align} where $e$ is a vector of 1s and $x=v-w$. From linear programming theory we know that a solution exists with at most $m$ nonzero $v_j$ or $w_j$, and that at most one of each pair $(v_j, w_j)$ will be nonzero for all $j=1$ to $n$. If the columns of $A$ contain ``good'' basis functions for $b$, almost all elements of $v$ and $w$ will be zero. Since signals are likely to be noisy, a more realistic model is the BP denoising problem \begin{align} \label{eqn:BPDN} \min_{x,\,r} \ \lambda\norm{x}_1 + \half r\T r \text{ \ subject to \ } Ax + r = b, \tag*{(BPDN)} \end{align} in which the 1-norm of $x$ is balanced against the 2-norm (squared) of a residual vector $r$. The size of the scalar parameter $\lambda$ greatly affects the sparsity of $x$. \subsection*{The GAMS model ajax.gms} The GAMS Model Library \url{http://www.gams.com/modlib/modlib.htm} contains many types of optimization problem (LP, NLP, MIP, MINLP, \dots) arising from many application areas (agricultural economics, chemical engineering, finance, \dots). The Management Science and OR section contains a small linear programming model named ``ajax''. It models a paper manufacturer wishing to find a production plan that maximizes monthly profits as the company produces four grades of paper on three machines, given a fixed machine availability (hours per month) and a fixed demand (tons per month). Also known are the machine production rates (tons per hour) and production costs and selling price (dollars per ton) for each kind of paper. The variables are a 2D array \verb|outp(g,m)| of nonnegative outputs (tons per month of each grade of paper by each machine). \subsection*{File ajaxBP.gms} We have taken the ajax model as a simple but realistic example of a linear program \begin{align} \label{eqn:LP} \min_x \ c\T x \text{ \ subject to \ } Ax = b, \ x \ge 0. \tag*{(LP)} \end{align} Our aim is to learn some parts of the GAMS language while experimenting with a Basis Pursuit version of the LP model. Rather than writing new GAMS statements (nontrivial without a manual!),\ you will simply study and run the modified GAMS model and conduct some ``good guess'' reverse engineering. Download file {\tt ajaxBP.gms} from \url{http://www.stanford.edu/class/msande318/} by clicking on {\tt Homework} and then on the filename. \subsection*{Project tasks} \begin{enumerate} \item Write down the BP problem associated with problem \ref{eqn:LP} above. In other words, which LP problem is likely to give a sparse feasible solution for \ref{eqn:LP}? Note that we ignore the $c\T x$ objective, but the bounds $x \ge 0$ are part of problem \ref{eqn:LP}. The required BP problem doesn't need to be as big as problem \ref{eqn:BPLP}. \item The {\tt Sets} statement defines two sets {\tt m} and {\tt g}, giving names to each machine and each grade of paper. What are the names and dimensions of four associated data matrices? (They are defined with varying syntax.) \item Note the {\tt Variables} and {\tt Positive Variable} statements. In your judgement, which upper and lower bounds do the variables have? \item Note the {\tt Equations} statement. The original linear objective in the ajax model is defined by a certain linear equation. What is the name of that equation, and what is the equation itself? \item Suppose we know that {\tt profit} is going to be nonnegative. We might try to help the optimizer by giving {\tt profit} a lower bound of zero. Can you give a reason why this would not be a good idea? \item Note the {\tt Models} statement. Write down algebraically the linear constraints and bounds in the original ajax LP model. \item We have added two new equations {\tt BP} and {\tt BPDN} to define the objective function for each of two new models. Write these objective functions algebraically. \item Note that {\tt sqr(r)} is a built-in function for specifying $\hbox{\tt r}^2$. We could have written {\tt r**2}, but GAMS would evaluate this as {\tt exp[2*log(r)]}. Why would this be a danger? \item GAMS allows us to write {\tt power(r,n)} if {\tt n} is an integer (positive or negative). Can you suggest why {\tt sqr(r)} is better if we know that {\tt n} = 2? (Remember that GAMS is generating gradients for us behind the scenes.) \item Transfer {\tt ajaxBP.gms} to your account on the campus file system. (If you wish, you can copy it from \url{/afs/ir.stanford.edu/class/msande318/WWW/homework/ajaxBP.gms}.) Then logon to one of the {\tt bramble} machines (64-bit 4-core Intel Xeon 3.00GHz Dell PowerEdge 1850 servers running Ubuntu 9.04) and run our three models by typing {\tt gams ajaxBP} at the command line. The output file will be called {\tt ajaxBP.lst}. Note that each {\tt Solve} statement specifies a model, a problem type, and an objective function. Which solvers did GAMS use by default for each model? How many iterations did each solve require? \item Typically the solvers output logfiles as they iterate. You can tell GAMS to include the logfile in its {\tt .lst} file by setting {\tt option sysout = on} before the {\tt Solve} statement. It is also normal to specify certain runtime options for the solver within a solver-specific file named \verb|.opt|. Download file {\tt ajaxBPDN.gms}. This has the first two solves (and reporting) omitted, but includes the following lines: \begin{footnotesize} \begin{verbatim} option sysout = on; ajaxBPDN.optfile = 1; Solve ajaxBPDN using nlp minimizing BPDNobj; \end{verbatim} \end{footnotesize} Prepare a file called {\tt minos.opt} as follows: \begin{footnotesize} \begin{verbatim} LU factor tolerance 2.0 LU update tolerance 2.0 Log frequency 1 Solution Yes \end{verbatim} \end{footnotesize} Now type {\tt gams ajaxBPDN nlp=minos} to run the single nonlinear model using MINOS as solver. The aim is to gain experience with controlling GAMS and the solver as well as observing the solver progress. (The MINOS log lines are a bit longer than the default GAMS line-length, but you can edit them a little to give just one line per iteration.) \item Finally, run GAMS on the BPDN model with SNOPT as solver. The same solver options can be used (they must be in another file called {\tt snopt.opt}). Study the MINOS and SNOPT output files and write some comments about how they compare. \end{enumerate} \subsection*{Au revoir} Thank you for being a gentle class. I hope we will meet again every so often. I also hope you encounter lots of matrix problems and optimization problems in the coming years, and that you'll immediately know which software you can put to work on them! \begin{footnotesize} % \bibliographystyle{alpha} \frenchspacing \bibliographystyle{plain} \bibliography{../notes/refs} \end{footnotesize} \end{document} %%% Local Variables: %%% mode: latex %%% TeX-master: t %%% End: