minimum enclosing rectangle for k out of n points

Finding the smallest axis-aligned rectangle that contains all \(n\) given points is very easy. We can simply do:

\[\begin{align}& {\color{darkred}x}_{lo} := \min_i {\color{darkblue}p}_{i,x} \\ & {\color{darkred}x}_{up} := \max_i {\color{darkblue}p}_{i,x} \\ &  {\color{darkred}y}_{lo} := \min_i {\color{darkblue}p}_{i,y} \\ &{\color{darkred}y}_{up} := \max_i {\color{darkblue}p}_{i,y}  \\ & {\color{darkred}{\mathit{area}}} := ({\color{darkred}x}_{up}-{\color{darkred}x}_{lo}) \cdot ({\color{darkred}y}_{up}-{\color{darkred}y}_{lo})  \end{align}\] 

We don't need an optimization model for this. A much more interesting problem is: find the smallest enclosing axis-aligned rectangle for \(k\lt n\) points, where the algorithm or model can choose the subset of \(k\) points. An algorithm for this is [1]. Here, of course, we try to solve this using an optimization model.

Data

I use a small dataset here with \(n=100\) data points. 

Dataset

----     37 PARAMETER p  points

                   x           y

point1         4.973       9.639
point2         8.739       8.443
point3         1.912       4.945
point4         2.346       3.810
point5         6.009       1.515
point6         5.077       4.633
point7         3.591       1.526
point8         3.664       8.260
point9         0.239       0.593
point10        2.848       1.140
point11        8.905       3.743
point12        2.341       1.477
point13        0.191       0.784
point14        9.280       5.550
point15        6.164       3.605
point16        5.214       4.320
point17        7.499       8.174
point18        0.160       3.729
point19        6.971       7.669
point20        8.875       8.117
point21        7.897       8.036
point22        0.166       3.932
point23        9.601       7.215
point24        8.147       8.375
point25        7.745       5.283
point26        6.557       3.924
point27        4.851       4.026
point28        4.346       1.935
point29        1.290       4.876
point30        7.912       6.856
point31        3.026       0.001
point32        4.463       8.148
point33        2.080       3.360
point34        2.816       1.565
point35        6.203       8.599
point36        6.153       3.857
point37        7.633       7.315
point38        0.798       0.829
point39        0.192       2.454
point40        3.341       3.038
point41        8.349       3.863
point42        0.531       5.056
point43        1.919       2.345
point44        9.347       7.861
point45        9.435       1.718
point46        1.817       2.026
point47        5.323       5.065
point48        2.108       6.382
point49        2.959       5.533
point50        0.342       8.075
point51        3.551       9.887
point52        8.136       7.153
point53        2.435       5.417
point54        5.019       1.108
point55        1.686       2.869
point56        8.416       8.060
point57        5.631       1.570
point58        2.715       0.339
point59        7.184       7.050
point60        1.577       7.522
point61        9.554       7.855
point62        2.566       0.880
point63        7.997       2.193
point64        9.556       8.034
point65        4.599       2.022
point66        8.957       5.241
point67        8.720       9.062
point68        3.519       3.229
point69        9.245       6.341
point70        4.141       2.687
point71        2.817       7.412
point72        6.571       4.559
point73        3.806       4.961
point74        8.999       5.323
point75        0.602       3.960
point76        4.302       9.600
point77        2.717       8.777
point78        1.872       4.952
point79        3.734       4.279
point80        9.519       5.120
point81        4.618       4.190
point82        0.915       9.232
point83        9.528       6.382
point84        5.353       8.983
point85        0.462       4.249
point86        8.164       4.834
point87        8.054       3.067
point88        8.594       7.318
point89        5.609       3.905
point90        2.741       4.446
point91        9.810       2.776
point92        7.988       9.737
point93        5.810       7.965
point94        8.287       7.178
point95        4.339       2.033
point96        2.708       6.648
point97        1.108       8.467
point98        2.175       2.350
point99        8.997       8.511
point100       8.928       2.236

Nonconvex optimization model

We can design an optimization model around two concepts:

• Select \(k\) out of \(n\) data points
• Draw a rectangle around these points

Then minimize the size (i.e. area) of the rectangle.

A complete optimization model for the \(k\) out of \(n\) minimum rectangle problem can look like:

 

High-level model

\[\begin{align} \min\> &{\color{darkred}{\mathit{area}}} = \prod_c ({\color{darkred}{\mathit{rect}}}_{c,up}-{\color{darkred}{\mathit{rect}}}_{c,lo}) \\ & {\color{darkred}\delta}_i = 1 \implies  {\color{darkblue}p}_{i,c} \in [{\color{darkred}{\mathit{rect}}}_{c,lo},{\color{darkred}{\mathit{rect}}}_{c,up}] \\ & \sum_i{\color{darkred}\delta}_i =  {\color{darkblue}k}  \\ & {\color{darkred}\delta}_i \in \{0,1\} \\ & c \in \{x,y\} \end{align}\]

The objective is quadratic and non-convex, so we need to use a global solver. 

The indicator constraint can be rewritten as big-M constraints: \[\begin{align} &{\color{darkred}{\mathit{rect}}}_{c,lo} \le {\color{darkblue}p}_{i,c} {\color{darkred}\delta}_i + {\color{darkblue}M}(1-{\color{darkred}\delta}_i) \\ & {\color{darkred}{\mathit{rect}}}_{c,up} \ge {\color{darkblue}p}_{i,c} {\color{darkred}\delta}_i \end{align}\] Note that we know that \[{\color{darkblue}p}_{i,c} \in [0,{\color{darkblue}M}], {\color{darkblue}M}=10\] Obviously, \[{\color{darkred}{\mathit{rect}}}_{c,lo},{\color{darkred}{\mathit{rect}}}_{c,up} \in [0,{\color{darkblue}M}]\] If we want we can tighten this slightly: \[\begin{align}&{\color{darkred}{\mathit{rect}}}_{c,lo},{\color{darkred}{\mathit{rect}}}_{c,up} \in [\min_i {\color{darkblue}p}_{i,c},\max_i {\color{darkblue}p}_{i,c}]\end{align}\] It is always a very good idea to bound all decision variables carefully in non-convex models. 

At first reading it may look like I have a single big-M value \(M\). This is not the case. The coefficients of the \({\color{darkred}\delta}_i\) variables are actually \({\color{darkblue}p}_{i,c}-{\color{darkblue}M}\) and  \({\color{darkblue}p}_{i,c}\). These are usually smaller in size. 

Some solvers may benefit from an extra linear constraint: \[{\color{darkred}{\mathit{rect}}}_{c,lo} \le {\color{darkred}{\mathit{rect}}}_{c,up}\]

We can simplify the objective a bit: \[\begin{align}\min \>& {\color{darkred}d}_{x}\cdot {\color{darkred}d}_{y} \\ &{\color{darkred}d}_{c} = {\color{darkred}{\mathit{rect}}}_{c,up}-{\color{darkred}{\mathit{rect}}}_{c,lo} \\ & {\color{darkred}d}_{c} \in [0,{\color{darkblue}M}] \end{align}\] Now we no longer need that last extra constraint.

Note that all constraints are linear. 

ChatGPT suggested to explicitly allow more than \(k\) points in the rectangle by changing the equality constraint into an inequality. We already allow more points in the rectangle. In the picture below, we could have some blue points inside the rectangle. As long as this does not affect optimality, that is fine. I am not sure why that proposed change is a good idea.

 

Solution

----    101 VARIABLE rect2.L  contain k points

          min         max

x       7.184       9.601
y       4.834       8.511


----    101 VARIABLE delta.L  point is selected

point2  1.000,    point14 1.000,    point17 1.000,    point20 1.000,    point21 1.000,    point23 1.000,    point24 1.000
point25 1.000,    point30 1.000,    point37 1.000,    point44 1.000,    point52 1.000,    point56 1.000,    point59 1.000
point61 1.000,    point64 1.000,    point66 1.000,    point69 1.000,    point74 1.000,    point80 1.000,    point83 1.000
point86 1.000,    point88 1.000,    point94 1.000,    point99 1.000


----    101 VARIABLE area.L                =        8.887  size of rectangle

Looking at the distribution of the data points, the location of the rectangle seems to make sense: there is a bit of a cluster of points there in the right upper corner.

This is not a very difficult model to solve. The global MINLP solver Antigone shows:

-------------------------------------------------------------------------------
Time (s) Nodes explored Nodes remaining Best possible   Best found Relative Gap
-------------------------------------------------------------------------------
 
       0              1               1    +0.000e+00   +1.743e+01           --
       3             70               0    +8.886e+00   +8.887e+00   +1.000e-04
 
-------------------------------------------------------------------------------
Termination Status : Global minimum
Best Feasible Point: +8.887291e+00
Best Possible Point: +8.886402e+00
       Relative Gap: +1.000000e-04

Conclusion

This is a small, simple model that allows us to find global, proven solutions. 

References

1. Timothy M. Chan and Sariel Har-Peled, Smallest k-Enclosing Rectangle Revisited, 2019, https://arxiv.org/abs/1903.06785.

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.

Minimum enclosing circle/ellipse 2

In [1] where I discussed how to find the minimum enclosing circle and minimum enclosing ellipse around a set of points. This is a follow-up post where I extend this to sets of circles and ellipses. 

1. MINIMUM ENCLOSING CIRCLE

Here our data is a set of \(n\) circles (or disks) of different size. We want to find the smallest circle that contains all these circles. 

An example data set with random circles can look like:

----     30 PARAMETER circles  coordinates of center and radius

                   x           y           r

circle1        4.294      21.082       1.663
circle2       13.759       7.528       4.014
circle3        7.305       5.601       1.961
circle4        8.746      21.407       6.235
circle5        1.678      12.505       2.591
circle6       24.953      14.468       2.715
circle7       24.778      19.056       4.564
circle8        3.267      15.993       5.336
circle9        3.988       6.252       4.769
circle10      16.723      10.884       3.783
circle11       8.993       8.786       3.480
circle12       3.287       3.753       1.706
circle13      14.728      20.772       2.885
circle14       5.770      16.643       1.279
circle15      19.396       7.591       3.031

Minimum enclosing ellipse

Minimum encompassing circle and ellipse

This is again about finding the smallest geometric shape containing all our data points. Here I focus on easy, convex cases: a circle and an ellipse. After the last post, where I was struggling mightily with a non-convex version of this model, this should be a breeze. 

I'll discuss SDP (semidefinite programming), SOCP (second order cone programming), rotated SOCP and traditional NLP approaches using gradient based solvers. So there is a lot of ground to cover.

First let's do circles.

Revisiting a crazy global NLP problem

Minimum encompassing triangle

This looks like a simple problem. Given \(n\) 2d points, find the smallest encompassing triangle. I follow the formulation from [1]. This post is self-contained, so you don't have to go back to [1]. The main contribution here is how one could work around the performance issues of the original formulation. 

Summary

The problem can be formulated as a nonconvex quadratic problem. However, this leads to models that are very difficult for global NLP solvers: global optimality can not be proven in reasonable time, even for very small data sets. One possible way to attack the problem is to use a different measure for size of the triangle: circumference instead of the more obvious area of a triangle. Although there are explicit algorithms for this problem [3,4], here we are trying to use an optimization model.