A note on variable splitting: complementarity

This note is about a well-known fact in linear programming: variable splitting and minimizing absolute values lead to complementarity pairs \[\begin{align}& r_i^{\bf +} \cdot r_i^{\bf -} = 0\\ & r_i^{\bf +}, r_i^{\bf -} \ge 0\end{align}\] in the solution. I.e., only one of  \(r_i^{\bf +}, r_i^{\bf -}\) can be nonzero. A simple example of this absolute-value modeling technique is implementing an LAD (Least Absolute Deviation) regression model through linear programming.

Let's have a look.

Problem: LAD regression

I created an artificial dataset with a few possible outliers:

----     25 PARAMETER data  

                        x           y

observation1        0.442       6.448
observation2        3.823      10.020
observation3        6.231       7.789
observation4        5.727       9.704
observation5        4.607      12.430
observation6        8.779      11.618
observation7        9.277      14.677
observation8        5.934      28.955
observation9        2.727       6.060
observation10       8.390      15.843
observation11       7.665      11.440
observation12       5.798       8.444
observation13       9.096       7.256
observation14       7.773      11.288
observation15       4.459       7.596
observation16       1.593       6.750
observation17       3.951       7.491
observation18       2.772      11.990
observation19       1.752      11.374
observation20       4.911      11.558
observation21       5.613      12.916
observation22       6.168      11.604
observation23       7.551      11.726
observation24       5.057       6.750
observation25       2.547       6.266
observation26       8.157      11.052
observation27       2.398       8.139
observation28       5.162      10.734
observation29       8.787      12.720
observation30       2.949      11.679

The outliers will heavily impact standard least-squares regression, but using absolute deviations reduces its impact.

Solver hangs

I don't see this that often: a solver that completely hangs. 

It keeps using CPU cycles.

Without knowing the solver's internals, my (brave) guess is an infinite loop. Note to myself: in production code, loops should always be protected (simply by an iteration count, or maybe a timer). I guess some exceptions are the operating system and a GUI event loop.

Rectilinear Clustering: area or perimeter?

In the previous post, I discussed covering \(k\) of \(n\) 2d points with the smallest possible rectangle. Minimizing the area yields a nonconvex quadratic objective. Solving these models requires global solvers. If we use the perimeter as a measure of size, the problem becomes linear. An example suggested that the solutions of these two approaches are similar. It certainly looked that way:

The picture shows the smallest rectangle containing \(k=25\) of a hundred random points. Minimizing either the area or the perimeter yields very similar rectangles. This comes from [1], where I noted that in this case the quadratic model (min area) and the linear model (min perimeter) yield similar solutions, but this is not guaranteed in general.

Here I discuss a different, but related problem: cover all points with \(k\) rectangles of minimum size. This problem is a well-known clustering problem that uses Manhattan or \(L1\) distances instead of the more traditional Euclidean distance. Again, we can choose area- or perimeter-based objectives. It is interesting to see the difference.

 

minimum enclosing rectangle for k out of n points

The problem of finding the smallest enclosing axis-aligned rectangle around a given set of points is trivial. We don't need an optimization model for that.

Things are more interesting if we need to cover just k out of n points. The model needs to select the k points and then form a rectangle around them such that the size of the rectangle (the area) is minimized. 

In this example we have n=100 and k=25. 

There are algorithms for this (see e.g. [1,2]). Here, of course, we try to solve this using an optimization model that finds proven, globally optimal solutions.

Mixture models as math programming problem

We can formulate a linear least squares regression model as an optimization problem. This is not how these problems are solved in statistical packages. Often, they use a QR decomposition. A real optimization formulation can be useful when we need to add unusual constraints that the statistical package does not support directly, or when we need to optimize a special maximum-likelihood function. 

Here is another example of a least-squares regression problem where we can benefit from mathematical programming techniques.

Data

Statistical problems typically start with a data set:

----     79 PARAMETER data  

                 x           y

case1       20.202      85.162
case2        0.507       2.103
case3       26.961      55.969
case4       49.985      44.690
case5       15.129      86.515
case6       17.417      79.866
case7       33.064      56.328
case8       31.691      29.422
case9       32.209      64.021
case10      96.398      85.191
case11      99.360      68.235
case12      36.990      57.516
case13      37.289      25.884
case14      77.198      56.157
case15      39.668      58.398
case16      91.310      66.205
case17      11.958      93.742
case18      73.548      28.178
case19       5.542       5.788
case20      57.630      60.830
case21       5.141      53.988
case22       0.601      42.559
case23      40.123      61.928
case24      51.988      42.984
case25      62.888      58.308
case26      22.575       4.414
case27      39.612      67.282
case28      27.601      56.445
case29      15.237       0.218
case30      93.632      11.896
case31      42.266      60.515
case32      13.466      51.721
case33      38.606      65.392
case34      37.463      16.978
case35      26.848      74.588
case36      94.837      -0.803
case37      18.894      60.060
case38      29.751      14.005
case39       7.455      60.066
case40      40.135      62.898

It looks like there is not a single regression line, but rather three of them hidden in this data set. So we should find something like:

 

The question is: how can we estimate the intercept and slope of these three lines?

MINLP instead of indicator constraints?

In this post, I want to discuss indicator constraints and how to replace them with simple multiplications. As we shall see, this is a somewhat harebrained but still interesting idea. 

Indicator constraints are implications of the form: \[\delta=0 \implies \text{linear constraint}\] or \[\delta=1 \implies \text{linear constraint}\] where \(\delta \in \{0,1\}\) is a binary decision variable.

There are two aspects of indicator constraints: 

• Indicator constraints help with MIP models where otherwise we would use big-M constraints. This will help address the numerical issues that result from using big-M constraints. These include small (and sometimes not-so-small) solution values where you expect zero, and poor solver performance. With big-M constraints, we need to pay much attention to the size of the big-M constants. This is the algorithmic aspect.
• Indicator constraints provide a convenient modeling construct. They form a useful abstraction that makes MIP modeling easier and more straightforward. This is the modeling aspect.

Another example where a solver feature has such a dual benefit (algorithmic and modeling) are SOS2 sets for piecewise linear functions. 

Experiments with Hostile Brothers nonconvex NLP model

In this post, let's do some experiments with the "hostile brothers problem." We have a plot of land. In my test models \([L,U]\times [L,U]\) with \(L=0\), \(U=100\)). We want to place \(n\) brothers on this plot. As they don't get along, we want to spread them out by maximizing the distance between neighbors. In modeling terms, we create a maximin model that maximizes the minimum distance between any two brothers. This yields some non-convex problems, so I'll make the problem small: we just have \(n=10\) brothers. This is a very small, but awe-inspiring problem.

Largest Empty Rotated Square: lots of trigonometry

Here, I delve into the problem of finding the largest empty rotated square. Given \(n\) data points, find the largest square rotated by an angle \(\theta\), such that the square does not contain any of the points.  In [1], I discussed two special cases: 

• an axis-aligned square (angle is 0°),
• a diamond shape, which is a square rotated by 45°.

Here we focus on arbitrary angles. A picture can help a bit: this is the largest square of all the squares with a 25° angle, which does not cover any of the given data points and is complete inside the outer box:

There are two cases I want to consider:

1. The angle is given. In modeling terms, the angle is exogenous.
2. Let the model find the best angle. I.e., the angle is endogenous.

Largest empty shapes

Finding the largest empty shapes

This post is about finding empty regions (square, rectangle, or circle) in a (large) collection of given points. These are interesting little optimization problems.

Data

I generated \(n=100\) data points \((x_i,y_i), i=1,\dots,n\) drawn from a uniform distribution \[\begin{align}&x_i \sim U(0,10) \\ & y_i \sim U(0,10)\end{align} \] 

It is always a good idea to have a careful look at the data. I.e., print and plot it (for very large datasets, print or plot a subset), and collect some elementary statistics (min, max, average, etc.). In practical modeling, it is surprising how many problems arise from simple data errors.  

Convex hull models

Convex hull as an optimization problem

In the previous posts, constructing a convex hull (or rather, the set of extreme points) played a significant role: it was an easy way to reduce the size of the data sets, often by a large amount. Somehow, it escaped me that we can try to formulate this as what turns out to be a conceptually rather simple optimization problem. This has probably little or no practical value. But it remains an interesting exercise.