I am a full-time consultant and provide services related to the design, implementation and deployment of mathematical programming, optimization and data-science applications. I also teach courses and workshops. Usually I cannot blog about projects I am doing, but there are many technical notes I'd like to share. Not in the least so I have an easy way to search and find them again myself. You can reach me at erwin@amsterdamoptimization.com.

Friday, December 24, 2010

I use gmail almost exclusively and I am very happy with how it works. Recently there surfaced an interesting bug. Attachments can be downloaded individually or combined as a single zip file. It turns out gmail adds a few bytes to the zipped files. After unzipping I see:

Tuesday, December 21, 2010

For a scheduling that I could not solve in one swoop I tried to develop a small algorithm. Basically: schedule the m<n jobs that come first, then fix the m jobs and schedule the next m jobs. This seems to work with Cplex. The algorithm looks like:

Sunday, December 5, 2010

Wed Nov 24 12:59 UTC -- On the evening before Thanksgiving, an IP located in Tbilisi, Georgia started an attack targeting the savannah.nongnu.org website. The perpetrators used SQL injection attacks to download the entire database of usernames and hashed passwords, and we should assume anything else in the Savannah MySQL database.

One way to prevent this is to solely use prepared SQL statements when accessing the DB. This is the strategy I have used in http://yetanothermathprogrammingconsultant.blogspot.com/2010/10/running-large-models-using-web.html. (Other approaches often mentioned are using an appropriate quoting mechanism, e.g. through mysql_real_escape_string and/or disabling multi statements) I browsed through a few recent books on PHP for this project, and I was amazed how much text is devoted to security; it looked like a third of the material seems to be about security. It truly sends shivers down your spine. Programming web applications becomes more and more a question of being paranoid all the time; this is not really a good development as it takes the fun out of things and also reduces productivity (extra work for code that does not implement any new and possibly exciting features for regular, non-hostile users). From a different perspective one could also say that the development tools are still too primitive such that developers are not isolated from these security issues and can concentrate on the real task at hand.

Wednesday, December 1, 2010

For different projects I have been using some of the following document sharing tools:

DropBox. This allows you to have a directory that is shared between different users. Advantage: on Windows this just looks like a network drive and you can easily copy files from and to the drop box directory. Disadvantages: just files. Also it looks for windows as a local disk so by default dragging means “move” instead of “copy”. Site: http://www.dropbox.com/.

Wiki-type environments. I have experience with PmWiki (http://www.pmwiki.org) and MediaWiki (http://www.mediawiki.org). Advantage: good support for HTML presentation of rich text. Disadvantage: somewhat complicated to setup and to secure; file upload facilities may be limited depending on the hosting environment.

My current favorite: sites.google.com. Clean and well-thought through. Has most of the facilities I need. Easily shared with just the people you want to share it with. Limits: 20 GB files, and 100 GB per site; I had to park some of the larger data files somewhere else.

The original Cube parameter has 326,755 elements, and this operation basically reshuffles things into a lower-dimensional parameter with some mappings applied. Upon completion the new parameter Download_Data has the same number of elements: 326,755. Unfortunately this operation was very expensive, so much so that I never saw it finish.

GAMS has some built-in facilities to automatically rearrange assignments like these to achieve better performance, but in this case it was not very effective. So human intervention was needed to optimize this assignment.

Tuesday, November 16, 2010

We are excited to announce the immediate availability of Cluster GPU Instances for Amazon EC2, a new instance type designed to deliver the power of GPU processing in the cloud. GPUs are increasingly being used to accelerate the performance of many general purpose computing problems. However, for many organizations, GPU processing has been out of reach due to the unique infrastructural challenges and high cost of the technology. Amazon Cluster GPU Instances remove this barrier by providing developers and businesses immediate access to the highly tuned compute performance of GPUs with no upfront investment or long-term commitment.

Amazon Cluster GPU Instances provide 22 GB of memory, 33.5 EC2 Compute Units, and utilize the Amazon EC2 Cluster network, which provides high throughput and low latency for High Performance Computing (HPC) and data intensive applications. Each GPU instance features two NVIDIA Tesla® M2050 GPUs, delivering peak performance of more than one trillion double-precision FLOPS. Many workloads can be greatly accelerated by taking advantage of the parallel processing power of hundreds of cores in the new GPU instances. Many industries including oil and gas exploration, graphics rendering and engineering design are using GPU processors to improve the performance of their critical applications.

Amazon Cluster GPU Instances extend the options for running HPC workloads in the AWS cloud. Cluster Compute Instances, launched earlier this year, provide the ability to create clusters of instances connected by a low latency, high throughput network. Cluster GPU Instances give customers with HPC workloads an additional option to further customize their high performance clusters in the cloud. For those customers who have applications that can benefit from the parallel computing power of GPUs, Amazon Cluster GPU Instances can often lead to even further efficiency gains over what can be achieved with traditional processors. By leveraging both instance types, HPC customers can tailor their compute cluster to best meet the performance needs of their workloads. For more information on HPC capabilities provided by Amazon EC2, visit aws.amazon.com/ec2/hpc-applications.

Sunday, October 31, 2010

For a large LP model a customer of mine wanted to have instead of a V shaped penalty around a goal, something that approximates a U shaped penalty. With standard slacks:

we essentially have a V shaped penalty:

We can approximate a more U shaped form by introducing extra slacks with a smaller penalty:

with unit penalty q smaller than p, e.q. q=0.1*p. This will give:

The idea is that the algorithm will select the cheap slacks t first, before spending on the expensive slacks s. This will keep the LP linear and we don’t need extra binary variables. For q=0 we essentially allow small deviations for free, but start paying when we go beyond the threshold T.

This approach can easily be extended to more piecewise linear segments (as long as we keep things convex). Also the range T for the slacks t may be generalized to:

Tuesday, October 26, 2010

For a project I am considering a fairly large GAMS models where we want end-users to run certain scenarios. Of course there must be HTML GUI to allow the user to specify the input for each scenario to run. That is fairly standard stuff – once you know what type of parameters you want to expose to the user to change. The bigger issue is that the model is too large to run immediately if the user presses the Submit button. Rather we need to create a small queueing mechanism where jobs are waiting to run. After the model has been solved the user can inspect the results.
Running background jobs from a web-server environment is a little bit more exotic, although other examples exist (see e.g. NEOS and this effort). For this project we use a Windows IIS web server with PHP.

Architecture

A background solve process can look like:

If background processes are already running, stop

Start a background process as follows:

While there are waiting jobs

select a job where status=waiting

set status=running

solve model

store results

set status=solved

We run this thing after every INSERT in the job table.
Actually we can have several of those processes if the server is beefy enough to allow for parallel solution of different models. This does not only depends on the number of cores available but also memory considerations play a role. Essentially line 1 becomes:

Wednesday, October 20, 2010

Savepoint/loadpoint is an easy way to use an initial point for an NLP or advanced basis for an LP. The commands are:

The first time you run this the solution is saved in the GDX file. The next time around, this GDX is used to set as many levels and marginals that can be matched by name. For this very large MCP model it turned out to be useful to bring down solution times of related experiments.

Sorry for the ugly syntax. I would have preferred something more symmetric for reading and writing.

Monday, October 11, 2010

After spending time and effort on a model, it may well be worth to pay attention to presentation and dissemination of the model results. Currently I am working on a project where we can export results for a very large and complex GAMS model to an Adobe Flash application called StatPlanet. Obviously a web-based tool to inspect results has lots of advantages (nothing to install, available from everywhere).

Thursday, October 7, 2010

GAMS has a binary data file format for which an API is available. In the last few weeks I encountered a number of similar issues with this in different projects. The interface does require more knowledge of the internals of GAMS than most people have. An example is here. When we want to read a GAMS parameter A(I,J) into a C/Fortran/Delphi matrix A[1..m,1..n] we see this is not completely trivial. The basic obstacle is that GAMS has strings as indices, while most programming languages have simple integers as indices. I.e. we need to map J = /new-york, chicago, topeka/ into j=1..3. For large problems that means some kind of hashing or other way to quickly look up strings. There is also a “raw” mode interface that does not give strings but rather internal GAMS numbers. So we need to map J = /1001,1002,1005/ to j=1..3. This again is not a completely trivial task. If you are familiar with the internals of GAMS one can know that GAMS internal numbers (so-called UEL numbers) are by definition increasing. So one could use a binary search to quickly find numbers. A fortran routine for this can look like: Another possibility would by to have an array of the size of MAXUEL that has mapping information from UEL-> j=1..3. That requires more space but has quick look up. If the GDX file contains only one set, one could use socalled “mapped” mode, but that cannot be used in general.

It is also worth noting that GAMS only stores and exports nonzero elements. So we cannot assume all A(i,j)’s are in the GDX file.

So to recap, my suggested approach to read a matrix A(i,j) from a GDX file is:

Read 1D set I (in raw mode) and store UEL’s in an integer array IA

Read 1D set J (in raw mode) and store UEL’s in an integer array JA

Allocate array A(i,j) and set all elements to zero

Read 2D parameter A (in raw mode). For each record do:

Retrieve (ueli, uelj, aij)

i := bsearch(IA, length(IA), ueli)

j := bsearch(JA, length(JA), uelj)

A[i,j] := aij

If you don’t want to use binary search an alternative algorithm would be:

Retrieve MAXUEL (gdxSystemInfo)

Allocate integer arrays IA and JA of length MAXUEL

Read 1D set I (in raw mode). For each record do:

Retrieve uel

IA[uel] := recordnumber

Read 1D set J into array JA as described in step 3

Read 2D parameter A (in raw mode). For each record do:

Retrieve (ueli, uelj, aij)

i := IA[ueli]

j := JA[uelj]

A[i,j] := aij

Notes:

The generated interface files are quite ugly, have strange design and implementation and have no comments on usage. The code in the provided files just looks very awkward. Below are a few more specific issues illustrating this.

The filenames are incomprehensible: gdxdocpdef.pas, gdxdcpdef.pas ???

The functions should return boolean/logical values instead of integers.

The previous point is the more crucial as some routines have 0=success and others have a return code <>0 meaning success. That throws me off balance quite often.

Fortran support for many compilers is missing.

The examples are too simplistic. The above is not explained at all.

The generated interface contains some code that I would refuse from my students when I was teaching programming at undergraduate level. An example is: n: ShortStringA symbol n should always be an integer, never a string (or a floating point number).

The Delphi interface does not compile under the latest compiler.

Delphi 2010

GAMS users routinely generate very large GDX files with millions of records. To get data in and out of these files efficiently still requires substantial skills and effort.

Sunday, October 3, 2010

For a very large large MCP model (a system of nonlinear equations with complementarity conditions) we observed it became unsolvable, even with a very good starting point. The reason was that the model was kept as small as possible in terms of equations and variables. As a result the model was much more dense and nonlinear (in terms of nonlinear nonzero elements) than needed. Basically, modeling systems and solvers that exploit sparsity should be fed large, sparse models instead of corresponding smaller, denser models. After substituting out some common expressions into separate equations (hence increasing the size in terms of equations and variables) we could solve the larger models (≈400K variables and equations) fairly easily.

An example of a nonlinear construct that should raise a red flag is something like:

In many cases it is better to write this as:

This will keep the number of nonlinear nonzero elements smaller at the expense of extra variables and equations.

In the model I was looking at a combination of these two situations was implemented: a large linear expression inside a nonlinear expression, and this nonlinear expression was repeated several times in the model. This caused both extra nonzero elements, and they appeared also in very nonlinear fashion.

For a smaller version of the model, the differences in statistics is striking (I was unable to get the large version to solve with the original formulation):

Monday, September 27, 2010

I get many mails from students asking for help. I usually try to give a small hint and direct them to get more help from their supervisor. Here are some examples.

I work for a supply chain management company, and we need some help on the attached exercise we have been given as an internal training exercise. Can you help? Attached is the project and all of the data and mod files. The deadline is Friday. If you can help, please provide a price quotation.

Answer: I would suggest to takes this up with your instructor. I believe it can be conceived as inappropriate if I would step into this.

After my standard reply of "Talk to your supervisor. He/she is actually getting paid to help you." one sometimes get long stories of poor supervision by professors:

I am really, really sorry to bother you, but my supervisor does not know how to use CPlex, only one lecturer in Math department knows how to use CPlex, but he is very busy at the moment (his wife just got cancer), you are the only expert that I can count on, it is a critical part of my PhD thesis, I really, really need your help for my problem, would you please have a look at my question and give me some idea? Thank you billion times! Your help will be greatly appreciated!

Or sometimes shorter:

Me and my wife are in very difficult situation. She runs out of the scholarship and I have to go back to Kuwait with her and the family.

Of course if I send back a file with some results, I really should hide all details about who actually did the work. Secrecy is important here:

Also can you change the following to STUDENTSNAME

Licensee: Erwin Kalvelagen G080731/0001CJ-WIN

Nobody know you are helping me please. My future is between your hands

Some people are quite brazen in their request. This guy wants me to do his complete PhD work and he will find someone else to help me if I am too expensive:

What is your lump thump price tell I get my PhD to get two publication accepted and modify my GAMS RESULTS tell the committee is satisfied and teaching me what you did and prepare me also for defence?

What is the budget for the first publication is $? while how much you should be able to write the second publication for $?. you will make sure that those publications are of such a high standard that you will be able to find a journal which will accept both papers.

In case that this is more than my budget allows, then we can just do the GAMS part for me. Once this is done we might as well find someone else who can write you the two publications.

Does not look PhD supervisors are doing a good job if students think they can get away with such requests!

It does not really work very well with the GAMS IDE. You cannot click on the red error messages to bring you to the input line with the error.

It does not handle combined use in equations and parameter assignments in post-solve reporting very well.

The documentation of the new $macro facility in GAMS (http://www.gams.com/docs/release/release.htm, section “The GAMS Macro Facility” ) seems to indicate that it can be used to replace the preprocessor tool: just use $macro and use set.local where needed. Unfortunately that is not completely true. Especially the sentence “The .local feature has been added to ensure local use and eliminates the need for new alias definitions.” is plain wrong.

Clearly the GAMS-F preprocessor handles this case correctly by introducing an alias. The $macro/.local combination does not do this and gives the wrong results. Instead of this .local hack it is needed to manually introduce an alias and use the aliased set where needed.

It is noted that overuse of these tools can lead to models that appear smaller (fewer equations and variables). However those models may be much denser and much more nonlinear (say in terms of nonlinear nonzero elements). Often larger, sparser and less non-linear models solve easier.

Monday, September 20, 2010

The first OR-tool from Google is a preliminary constraint programming solver, available from here: https://code.google.com/p/or-tools/. The name “or-tools” seems to indicate we can expect more to come.

Of course if λ=0 or λ=1 is part of the λ’s we need to trace, then this actually does not require extra models to solve, just a different ordering: first do λ=0 and λ=1 to get f0 and g0 and then do the λ’s in between. As the λ’s are user specified, we needed a little bit of logic to handle this.

Intuitively this makes some sense: MINOS and SNOPT scale by default only the linear variables. This strategy is actually not without merit: the Jacobian elements for these nonlinear variables can change quite a lot during the optimization, so scaling based on the initial Jacobian elements may not be very wise. However, in some cases making this distinction between linear and nonlinear variables can make things worse. In this case it is better not to scale at all.

I was first tinkering with the optimality tolerance, as I have seen with other forestry models (LP and NLP) that for these ridiculously long time horizons, tightening the optimality tolerance can help getting better solutions. For this model that was not the case.

Thursday, August 19, 2010

In some cases a new version of a solver can perform worse on some models. Here is a particular unfortunate example, where Cplex 12.1 finds several very good solutions in the initial heuristics, but Cplex 12.2 just fails to find any integer solutions within the allotted time:

Such behavior is always difficult to explain to a customer. The argument that on average a new version tends to do better is of little consolation.

PS. The time limit was very tight on this model: 400 seconds (extremely tight considering the LP takes 100 seconds). This was no problem with previous versions of Cplex as the initial heuristic always found solutions very quickly. These solutions were so good they already obeyed the requested gap of 0.8%. Cplex 12.2 finds the solution relatively quickly in the B&B after about 500 seconds on my machine with 1 thread (the client on his machine obtained no solution in 600 seconds with 2 threads and some other solver options; I assume a little bit more time may help here).

Saturday, August 14, 2010

Memory estimate for Hessian is too small. Increase the <modelname>.workfactor option from the current value of 1

For 2000 stocks with WORKFACTOR=10 I got the same error.

Even if the model is formulated as a QCP model, the GAMS link is not smart enough in setting up the hessian efficiently. Instead it looks like the QP is handled as a general NLP.

There are two easy workarounds. First of all you can solve this with IPOPT. The GAMS link for IPOPT is smarter in dealing with second derivatives than the link for KNITRO (as of GAMS 23.5). Note that it is better to solve as a QCP model than as an NLP model: GAMS is not smart enough to recognize this as a QP and exploit that:

model type

GAMS generation time

IPOPT time

function evaluation time

Total solver time

NLP

38

21.8

39.1

61.0

QCP

38

21.6

9.8

31.4

A second possibility is to reformulate model (1):

into model (2):

This formulation can be solved with KNITRO without the problems with the hessian storage.

As an aside we can improve formulation 1 by replacing the variance-covariance matrix Q by

We exploit here symmetry in the quadratic form x’Qx. This exercise substantially reduces the GAMS generation time, and also further shaves off some time from the solver:

model type

GAMS generation time

IPOPT time

function evaluation time

Total solver time

QCP –triangular Q

18.8

21.3

4.6

26

We can also apply the triangular storage scheme to model (2). With a solver like SNOPT this leads to a very fast total solution time (i.e. generation time + solver time).

model type

GAMS generation time

SNOPT time

NLP –triangular Q in linear constraints

1

5

If needed we can even speed this up by not using the covariances but rather the mean adjusted returns directly:

I used T=300 time periods. This formulation yields:

model type

GAMS generation time

SNOPT time

NLP – Mean Adj Returns

0.5

2.3

Even for the simple QP portfolio model, careful modeling can help to improve the performance in a significant way.