Wednesday, May 20, 2015

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.

Friday, May 15, 2015

R inside SQL Server

 

image

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”.

Tuesday, May 5, 2015

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.

image

image

image

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.