Microsoft has ownership of R?

From an (ex-) intellectual property lawyer at Adobe (http://www.infoworld.com/article/2940508/big-data/hey-microsoft-a-rewrite-of-the-r-language-is-a-silly-idea.html): 

“Microsoft, unlike Tibco, can claim ownership of the underlying R code. By buying Revolution Analytics, Microsoft gained that ownership”

I was worried. But luckily it looks like this was misinformation. The new version of the article says:

“Correction: An earlier version of this article incorrectly assumed that Microsoft acquired the rights to R Project from Revolution Analytics. As Revolution Analytics chief community officer David Smith has since stressed, "R is owned by the R Foundation," not Microsoft.”

MPS pictures

http://blogs.sas.com/content/operations/2015/06/25/finding-the-beauty-in-optimization-models-visualizing-mps-data-sets/

Somehow I find these pictures not very helpful in my modeling work. One reason is that problems that don’t display this structure may be very well decomposable. The ordering of rows and columns is very important to detect this structure visually. E.g. I often use a time index t as last index in my variables and equations. As the last index runs the fastest, that will not give these nice pictures. Also modeling systems will likely export all variables x() before y() – i.e. ordered by variable (equation) name. In general you will need to reorder rows and columns to make these pictures meaningful.

Nonlinear system test function

The following paper:

has an interesting test function for a large system of nonlinear equations:

This function originates from:

Here it has a fixed starting point:

Solving a triangular system

This system of nonlinear equations actually forms a triangular system. This means a good preprocessor can solve this just by solving small 1x1 problems. To be precise: first solve for x₁, then for x₂ etc. Note that these tiny 1x1 problems are possibly nonlinear. In this case they are for sure nonlinear. We write the model as:

I.e. 10k equations and 10k variables. Note that x(i-1) for the first i is zero in GAMS.

CONOPT

Using GAMS/CONOPT we see

CONOPT 3         24.4.5 r52258 Released May 26, 2015 WEI x86 64bit/MS Windows     C O N O P T 3   version 3.16D     Copyright (C)   ARKI Consulting and Development A/S                     Bagsvaerdvej 246 A                     DK-2880 Bagsvaerd, Denmark    Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK       0   0        7.3122859754E+06 (Input point)                                 Pre-triangular equations:    10000                                 Post-triangular equations:       0       1   0        0.0000000000E+00 (After pre-processing)       2   0        0.0000000000E+00 (After scaling) ** Feasible solution to a square system.

I.e. no “real” work was needed for this problem: the presolver was able to solve the complete problem. CONOPT is really the smart here compared to some other solvers. Below some results.

PATH, KNITRO

Note that with the PATH solver we were not able to solve this problem:

** EXIT - other error.

I was not able to see if the preprocessor did anything here.

With KNITRO we see:

KNITRO presolve eliminated 0 variables and 0 constraints.

It probably does only a linear presolve. This results in a subsequent failure:

EXIT: Terminate at infeasible point because the relative change in solution       estimate < xtol.

Tight big-M values

It is well known large big-M values can cause big headaches for MIP solvers. There are many examples where something like M=9999999999 is used, leading to very wrong results. In some cases we can use a MIP model to calculate tight values for big-M constants. In this case I have a model where the MIP model to calculate the tightest bound is actually using the big-M itself:

This is of course not a useful approach. Luckily there are a few alternatives:

• Keep (1) and most of (2), and we end up with a case where it is possible to find a bound by just looking at the data and so some simple calculations (no MIP is needed). From the problem (a design problem) we know this is a really good bound.
• More sophisticated we can relax only equations (3), and solve a MIP to find a slightly tighter bound. We need to do this for several y_k's. But the models are small and solve fast (we can even solve them in parallel). That is the approach I am going to take.

In GAMS we can write something like:

eq3(k)$(not relax(k)).. y(k) =L= <expr>;

to turn on or off an equation based on set relax. This allows us to re-use this equation unaltered in both this preprocessing step, in the real model and in some reporting step.

There actually is a “Lazy Bureaucrat Problem”

 

http://link.springer.com/chapter/10.1007/978-3-319-18008-3_18

Adding symmetry-breaking constraints

On some models symmetry poses a real problem for solvers. Here is an example of a sports related scheduling problem. We achieve a speedup of 35x by adding constraints that break a number of symmetries in the model. The results here are with Gurobi.

Modern solvers have also built-in facilities to detect and do something about symmetry. E.g. Gurobi has the option Symmetry which can be set to 2 (aggressive). I never had much luck with this. Indeed the results are:

Cplex shows very similar behavior:

Note that Cplex provides a large number of settings for the symmetry option:

That is probably just an educational tool to teach the user a lesson about combinatorial explosion: finding the right settings for your model.

Often solver suppliers say this is a solved problem:

This is an example that shows the fixes are at least insufficient. The algorithms do not handle this problem always in a satisfactory way.

CGE modeling in The Economist

http://www.economist.com/cge15. They even mention GTAP.

Data tables and data frames

Many modern languages have now something called a data frame or data table that can store tabular data with a mixture of numeric, string columns and columns of other types. This is much more flexible than a typical array (with elements having all the same type) and more organized than just say a list of lists. Of course relational databases have tables that have a similar structure.

R has data frames, which form arguably the most important data structure in R.

Python has a DataFrame as part of pandas, the toolkit for doing data analysis.

Matlab has a new Table data type. Their statistical toolbox still has a DataSet.

.Net has a DataTable class which I use often.

I am missing similar "standard" facilities in C++, Java (ResultSet is not really the same) and in modeling systems like AMPL and GAMS.

R inside SQL Server

 

http://blog.revolutionanalytics.com/2015/05/r-in-sql-server.html

My data is usually not that big that I cannot move it from the RDBMS to the place where I do the “real work”.

Model generation in Julia

In http://www.optimization-online.org/DB_HTML/2015/04/4891.html we see some interesting numbers on generation times of large models in different modeling environments. AMPL and Julia/JuMP are doing very well on these benchmarks. E.g.

As https://plus.google.com/+VictorZverovich/posts/RtHfqgwh5QX indicates: AMPL is faster than Gurobi/C++. This does not come as a surprise: AMPL is around quite a while and is known to be quite fast. GAMS is slow compared to AMPL especially on the lqcp and acpower instances. Note that in general solution times are much larger than the modeling system needs to generate the problem. Only very easy LPs can be relatively expensive to generate.

Why is GAMS so slow on these instances? In https://github.com/mlubin/JuMPSupplement we can find much information so we can reproduce things (very praiseworthy to make this available).

lqcp-2000

Erwin@EKlap2 ~ $ gams lqcp.gms lo=3  | ts -s 00:00:00 --- Job lqcp.gms Start 05/05/15 15:13:48 24.4.1 r50296 WEX-WEI x86 64bit/MS Windows 00:00:00 GAMS 24.4.1   Copyright (C) 1987-2014 GAMS Development. All rights reserved 00:00:00 Licensee: Erwin Kalvelagen                               G140808/0001CV-GEN 00:00:00           Amsterdam Optimization Modeling Group                     DC10455 00:00:00 --- Starting compilation 00:00:00 --- lqcp.gms(65) 3 Mb 00:00:00 --- Starting execution: elapsed 0:00:00.003 00:00:00 --- lqcp.gms(63) 260 Mb 00:00:00 --- Generating MIQCP model cont5 --- lqcp.gms(65) 2406 Mb) 284 Mb 00:01:03 ---   4,004,002 rows  4,006,002 columns  23,994,008 non-zeroes 00:01:03 ---   18,010 nl-code  4,001 nl-non-zeroes --- lqcp.gms(65) 2404 Mb) 2406 Mb 00:01:03 --- Executing GUROBI (Solvelink=5): elapsed 0:01:02.947 00:01:03 00:01:03 Gurobi           24.4.1 r50296 Released Dec 20, 2014 WEI x86 64bit/MS Windows 00:01:03 00:01:03 Gurobi full license. 00:01:03 Gurobi library version 6.0.0 00:01:06 Space for names approximately 174.82 Mb 00:01:06 Use option 'names no' to turn use of names off 00:01:07 Starting Gurobi... 00:01:07 Optimize a model with 4004001 rows, 4006001 columns and 23990006 nonzeros

We do this in 67 seconds: a little bit better than reported. Not much I can do to speed this up: the structure of the offending equation pde has leads and lags which is not something GAMS is always doing extremely fast. Usually my client models are large and sparse and do not have this structure.

I actually think the equation pde is not 100% correctly formulated. With these discretizations it is important to pay extra attention to the border. This is especially the case as GAMS will set all terms with leads and lags that exceed the boundary are implicitly set to zero. Always inspect the equation listing in GAMS when there are leads and lags involved. However these issues are probably not important for these measurements.

fac-100

Erwin@EKlap2 ~ $ gams facility.gms lo=3 | ts -s 00:00:00 --- Job facility.gms Start 05/05/15 16:01:12 24.4.1 r50296 WEX-WEI x86 64bit/MS Windows 00:00:00 GAMS 24.4.1   Copyright (C) 1987-2014 GAMS Development. All rights reserved 00:00:00 Licensee: Erwin Kalvelagen                               G140808/0001CV-GEN 00:00:00           Amsterdam Optimization Modeling Group                     DC10455 00:00:00 --- Starting compilation 00:00:00 --- facility.gms(41) 3 Mb 00:00:00 --- Starting execution: elapsed 0:00:00.008 00:00:00 --- facility.gms(39) 4 Mb 00:00:01 --- Generating MIQCP model facility --- facility.gms(41) 1324 Mb) 235 Mb 00:00:11 ---   4,090,601 rows  4,080,601 columns  11,221,100 non-zeroes 00:00:11 ---   10,201,001 nl-code  3,060,300 nl-non-zeroes 00:00:12 ---   1,020,100 discrete-columns --- facility.gms(41) 1271 Mb) 1323 Mb 00:00:13 --- Executing GUROBI: elapsed 0:00:11.390 00:00:16 00:00:16 Gurobi           24.4.1 r50296 Released Dec 20, 2014 WEI x86 64bit/MS Windows 00:00:16 00:00:16 Gurobi full license. 00:00:16 Gurobi library version 6.0.0 00:00:20 Space for names approximately 165.81 Mb 00:00:20 Use option 'names no' to turn use of names off 00:01:06 Starting Gurobi... 00:01:06 Optimize a model with 3070501 rows, 4080601 columns and 8160800 nonzeros

Here we are actually slower than reported. The pure generation time is not that bad: 11 seconds. The main problem we see here is that the GAMS/Gurobi link is inefficient in handling the setup of the quadratic terms for these large models. This turned out to be a problem with Gurobi: adding individual quadratic constraints is expensive. A fix is on its way. Update: fixed in Gurobi 6.0.4.

acpower-100

Erwin@EKlap2 ~ $ gams opf_66200bus.gms lo=3 | ts -s 00:00:00 --- Job opf_66200bus.gms Start 05/06/15 10:29:26 24.4.1 r50296 WEX-WEI x86 64bit/MS Windows 00:00:00 GAMS 24.4.1   Copyright (C) 1987-2014 GAMS Development. All rights reserved 00:00:00 Licensee: Erwin Kalvelagen                               G140808/0001CV-GEN 00:00:00           Amsterdam Optimization Modeling Group                     DC10455 00:00:00 --- Starting compilation 00:00:00 --- opf_66200bus.gms(17) 2 Mb --- .IEEE66200.bus.gms(1059233) 15 Mb3) 15 Mb 00:00:01 --- opf_66200bus.gms(46) 15 Mb --- .IEEE66200.branch.gms(1322139) 59 Mb2) 41 Mb 00:00:01 --- opf_66200bus.gms(239) 59 Mb 00:00:02 --- Starting execution: elapsed 0:00:01.581 00:00:02 --- opf_66200bus.gms(2381609) 93 Mb 00:00:02 --- Generating NLP model acpower --- opf_66200bus.gms(2381610) 313 Mb8) 107 Mb 00:04:54 ---   132,401 rows  416,601 columns  1,975,401 non-zeroes 00:04:54 ---   14,576,003 nl-code  1,960,800 nl-non-zeroes 00:04:54 --- opf_66200bus.gms(2381610) 239 Mb 00:04:54 --- Executing IPOPT (Solvelink=5): elapsed 0:04:54.060 00:04:55 00:04:55 COIN-OR Ipopt    24.4.1 r50296 Released Dec 20, 2014 WEI x86 64bit/MS Windows 00:04:55 00:04:55 COIN-OR Interior Point Optimizer (Ipopt Library 3.11) 00:04:55 written by A. Waechter. 00:04:55 00:04:55 Note: This is the free version IPOPT, but you could also use the commercially supported and potentially higher performance version IPOPTH. 00:05:51 00:05:51 ****************************************************************************** 00:05:51 This program contains Ipopt, a library for large-scale nonlinear optimization. 00:05:51  Ipopt is released as open source code under the Eclipse Public License (EPL). 00:05:51          For more information visit http://projects.coin-or.org/Ipopt 00:05:51 ****************************************************************************** 00:05:51 00:05:51 This is Ipopt version 3.11, running with linear solver mumps. 00:05:51 00:05:51 Space for names approximately 14.45 MB. 00:05:51 Use statement '<modelname>.dictfile=0;' to turn dictionary off. 00:05:51 Number of nonzeros in equality constraint Jacobian...:  1068500 00:05:51 Number of nonzeros in inequality constraint Jacobian.:        0 00:05:51 Number of nonzeros in Lagrangian Hessian.............:   805300 00:05:51 00:05:54 Total number of variables............................:   161800 00:05:54                      variables with only lower bounds:        0 00:05:54                 variables with lower and upper bounds:    95700 00:05:54                      variables with only upper bounds:        0 00:05:54 Total number of equality constraints.................:   132400 00:05:54 Total number of inequality constraints...............:        0 00:05:54         inequality constraints with only lower bounds:        0 00:05:54    inequality constraints with lower and upper bounds:        0 00:05:54         inequality constraints with only upper bounds:        0 00:05:54

This is pretty bad for a somewhat small and somewhat dense model. We can make the model easily larger and hopefully sparser and “less” non-linear. After adding some extra options to reduce the size of the listing file we see:

Erwin@EKlap2 ~ $ gams acpower_ek.gms lo=3 limrow=0 limcol=0 | ts -s 00:00:00 --- Job acpower_ek.gms Start 05/06/15 10:57:38 24.4.1 r50296 WEX-WEI x86 64bit/MS Windows 00:00:00 GAMS 24.4.1   Copyright (C) 1987-2014 GAMS Development. All rights reserved 00:00:00 Licensee: Erwin Kalvelagen                               G140808/0001CV-GEN 00:00:00           Amsterdam Optimization Modeling Group                     DC10455 00:00:01 --- Starting compilation 00:00:01 --- acpower_ek.gms(17) 2 Mb --- .IEEE66200.bus.gms(1059233) 15 Mb2) 15 Mb 00:00:01 --- acpower_ek.gms(46) 15 Mb --- .IEEE66200.branch.gms(1322139) 59 Mb7) 43 Mb 00:00:02 --- acpower_ek.gms(263) 59 Mb 00:00:04 --- Starting execution: elapsed 0:00:01.566 --- acpower_ek.gms(2381632) 108 Mb4) 100 Mb 00:00:07 --- Generating NLP model acpower --- acpower_ek.gms(2381634) 445 Mb4) 118 Mb 00:00:57 ---   671,601 rows  955,801 columns  3,683,001 non-zeroes 00:00:57 ---   11,999,503 nl-code  3,063,000 nl-non-zeroes 00:00:57 --- acpower_ek.gms(2381634) 384 Mb 00:00:57 --- Executing IPOPT (Solvelink=5): elapsed 0:00:56.518 00:00:58 00:00:58 COIN-OR Ipopt    24.4.1 r50296 Released Dec 20, 2014 WEI x86 64bit/MS Windows 00:00:58 00:00:58 COIN-OR Interior Point Optimizer (Ipopt Library 3.11) 00:00:58 written by A. Waechter. 00:00:58 00:00:58 Note: This is the free version IPOPT, but you could also use the commercially supported and potentially higher performance version IPOPTH. 00:03:56 00:03:56 ****************************************************************************** 00:03:56 This program contains Ipopt, a library for large-scale nonlinear optimization. 00:03:56  Ipopt is released as open source code under the Eclipse Public License (EPL). 00:03:56          For more information visit http://projects.coin-or.org/Ipopt 00:03:56 ****************************************************************************** 00:03:56 00:03:56 This is Ipopt version 3.11, running with linear solver mumps. 00:03:56 00:03:56 Number of nonzeros in equality constraint Jacobian...:  2586800 00:03:56 Number of nonzeros in inequality constraint Jacobian.:        0 00:03:56 Number of nonzeros in Lagrangian Hessian.............:  3461700 00:03:56 00:04:08 Total number of variables............................:   701000 00:04:08                      variables with only lower bounds:        0 00:04:08                 variables with lower and upper bounds:    95700 00:04:08                      variables with only upper bounds:        0 00:04:08 Total number of equality constraints.................:   671600 00:04:08 Total number of inequality constraints...............:        0 00:04:08         inequality constraints with only lower bounds:        0 00:04:08    inequality constraints with lower and upper bounds:        0 00:04:08         inequality constraints with only upper bounds:        0

We see again that GAMS generates the model quickly, but that inside the solver link we have some inefficiencies. I am not sure what is happening in the 3 lost minutes, but it could be related in setting up the second derivatives. Further improvements can be obtained by telling GAMS to consider fixed variables as parameters.

Some general points:

• Generation time is often a minor part of the total turn around time
• Solvers are in most cases much slower to solve the problem than the time the modeling system needs to generate the model
• In many models we spend more time in data manipulation than in pure model generation. Many practical models have most of the code devoted to data manipulation (data import, data conversion, checks, aggregation, consolidation etc.)
• I believe some models will behave better with upcoming version of GAMS and Gurobi. These very large examples are excellent to find weak spots in the modeling software and will lead to improvements.
• I am really impressed with the JuMP performance.
• The models are very large. Lots of meaningful models are way smaller than these (but may be more complex). Systems like Pyomo should not be dismissed just on these numbers.

The paper is an interesting read. It touches upon many concepts that are not widely known.