Saturday, September 26, 2026

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.

Data


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






Dataset 
----     43 PARAMETER k                    =           25  number of points to select

----     43 PARAMETER p  points

                   x           y

point1    4.97311132  9.63862756
point2    8.73922492  8.44273500
point3    1.91226253  4.94526579
point4    2.34553216  3.80983451
point5    6.00887463  1.51476321
point6    5.07716895  4.63323882
point7    3.59064347  1.52625715
point8    3.66418827  8.26045475
point9    0.23893508  0.59261597
point10   2.84759145  1.14012489
point11   8.90454986  3.74340300
point12   2.34069434  1.47696552
point13   0.19119162  0.78401936
point14   9.27971966  5.55039788
point15   6.16361826  3.60477909
point16   5.21416943  4.32016162
point17   7.49929678  8.17376279
point18   0.16037081  3.72949918
point19   6.97091674  7.66895295
point20   8.87457833  8.11673997
point21   7.89650246  8.03634493
point22   0.16576575  3.93244606
point23   9.60138480  7.21520825
point24   8.14734411  8.37501191
point25   7.74541956  5.28258128
point26   6.55663865  3.92384422
point27   4.85091115  4.02642582
point28   4.34623796  1.93472329
point29   1.28959803  4.87601902
point30   7.91163148  6.85615242
point31   3.02621818  0.00063517
point32   4.46278421  8.14820541
point33   2.08022534  3.35953314
point34   2.81600217  1.56499252
point35   6.20283636  8.59944445
point36   6.15258941  3.85719583
point37   7.63343999  7.31506452
point38   0.79832230  0.82942221
point39   0.19238302  2.45380498
point40   3.34088869  3.03838803
point41   8.34902953  3.86320590
point42   0.53110352  5.05611011
point43   1.91904761  2.34489699
point44   9.34705030  7.86066922
point45   9.43453787  1.71763581
point46   1.81723665  2.02565095
point47   5.32342079  5.06474382
point48   2.10786534  6.38173925
point49   2.95917598  5.53252693
point50   0.34180259  8.07512765
point51   3.55101998  9.88688936
point52   8.13607683  7.15280359
point53   2.43498983  5.41651376
point54   5.01929436  1.10750815
point55   1.68573148  2.86927244
point56   8.41569886  8.06049381
point57   5.63100458  1.57040011
point58   2.71516744  0.33909707
point59   7.18412193  7.05045720
point60   1.57713422  7.52170649
point61   9.55352881  7.85491794
point62   2.56555975  0.88005074
point63   7.99748484  2.19292948
point64   9.55568192  8.03379715
point65   4.59853976  2.02151223
point66   8.95713092  5.24080749
point67   8.71950503  9.06200769
point68   3.51902443  3.22910422
point69   9.24501444  6.34123137
point70   4.14098498  2.68712111
point71   2.81651324  7.41232701
point72   6.57107119  4.55851885
point73   3.80600433  4.96108423
point74   8.99943277  5.32278562
point75   0.60195961  3.95965993
point76   4.30151391  9.59964284
point77   2.71652476  8.77681007
point78   1.87229129  4.95157553
point79   3.73429995  4.27888100
point80   9.51936360  5.11992531
point81   4.61819146  4.18987215
point82   0.91512513  9.23192650
point83   9.52838375  6.38224088
point84   5.35277514  8.98266937
point85   0.46190997  4.24850010
point86   8.16439743  4.83398475
point87   8.05382905  3.06655684
point88   8.59419222  7.31792459
point89   5.60900115  3.90522405
point90   2.74138705  4.44598243
point91   9.81040377  2.77587450
point92   7.98768312  9.73693294
point93   5.81026191  7.96489037
point94   8.28671680  7.17813456
point95   4.33880861  2.03345798
point96   2.70771867  6.64833521
point97   1.10849014  8.46666734
point98   2.17469964  2.35013718
point99   8.99662815  8.51057738
point100  8.92751427  2.23560422


----     53 PARAMETER box  min/max values

          min         max

x       0.160       9.810
y       0.001       9.887


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}]  & \forall i,c\\ & \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.  

Many solvers support implications as indicator constraints. Typically, they want a linear constraint on the right-hand side. That would mean two indicator constraints in this case. However, we can also easily rewrite the indicator constraint as simple big-M constraints. 

I implemented the model as:

 
Implemented model
\[\begin{align} \min\> &{\color{darkred}{\mathit{area}}} =  {\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} & \forall c \\ & {\color{darkred}{\mathit{rect}}}_{c,lo} \le {\color{darkred}\delta}_i {\color{darkblue}p}_{i,c} +  (1-{\color{darkred}\delta}_i) {\color{darkblue}{\mathit{box}}}_{c,up} & \forall i,c \\ & {\color{darkred}{\mathit{rect}}}_{c,up} \ge  {\color{darkred}\delta}_i {\color{darkblue}p}_{i,c} + (1-{\color{darkred}\delta}_i){\color{darkblue}{\mathit{box}}}_{c,lo}& \forall i,c\\ &  \sum_i{\color{darkred}\delta}_i = {\color{darkblue}k} \\& {\color{darkred}{\mathit{rect}}}_{c,lo},{\color{darkred}{\mathit{rect}}}_{c,up} \in [{\color{darkblue}{\mathit{box}}}_{c,lo},{\color{darkblue}{\mathit{box}}}_{c,up}] \\ & {\color{darkred}d}_c \in [0,{\color{darkblue}{\mathit{box}}}_{c,up}-{\color{darkblue}{\mathit{box}}}_{c,lo}] \\ & {\color{darkred}\delta}_i \in \{0,1\}\\ &{\color{darkblue}{\mathit{box}}}_{c,lo} = \min_i{\color{darkblue}p}_{i,c} \\ &{\color{darkblue}{\mathit{box}}}_{c,up} = \max_i{\color{darkblue}p}_{i,c}\\ & c \in \{x,y\} \end{align}\]

Notes:
  • If we want to avoid large big-M coefficients (and we should), scale the data points. In our case, we generated points inside a \(([0,10]\times[0,10])\) square.
  • All constraints are linear.
  • To help global solvers, it is a good idea to use tight bounds on all variables.
  • 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 rectangle's location seems to make sense: there is a bit of a cluster of points in the upper-right corner.

Notes:
  • Solutions will have at least one selected point located on each side of the optimal rectangle. This property is not used in this model, but is key in algorithms for this problem.
  • It is possible that non-selected (blue) points are inside the optimal rectangle. This is not the case in this particular solution. It is easy to see that this can happen if multiple points share the same coordinates.  

This model is not very difficult 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.476e+01           --
       4             96               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

Algorithms


Algorithms for solving this problem are typically organized around clever enumeration. Rectangles are not placed in a continuous space: they always have a point exactly on the left, right, top, and bottom. This property makes the problem discrete and lets us search possible rectangles quickly. Interestingly, our optimization model actually uses a continuous framework for the rectangles (only selecting the points is discrete).

A simple algorithm could look like:

  1. Sort the points \(p_i\) on \(y\) coordinate.
  2. Look at windows \(p_{\ell},p_u\) with \(u\ge \ell + k - 1\) (at least \(k\) points in between)
  3. For each \(y\) window, sort on \(x\), and look at sliding \(x\) windows of exactly \(k\) points
  4. Calculate the rectangle and area, and record the best. 

When I used some Python code suggested by an LLM (ChatGPT), indeed we found the optimal solution quickly:

n=100,k=25
x bounds: [7.18412193, 9.6013848]
y bounds: [4.833984750000001, 8.51057738]
width: 2.41726287
height: 3.67659263
area: 8.88729085261465
selected points: [2, 14, 17, 20, 21, 23, 24, 25, 30, 37, 44, 52, 56, 59, 61, 64, 66, 69, 74, 80, 83, 86, 88, 94, 99]


Variant 1: Next best cluster


I often find it instructive to produce a different, next-best solution. The difference should be real: I am not that interested in a solution that is just an epsilon perturbation. One way to define a different solution is to forbid overlap. We could enforce no overlap in the rectangle (that would require extra binary variables), but it's much simpler to require no overlap among the selected points. We can achieve this by fixing all previously selected points: \({\color{darkred}\delta}_i=0\) when the optimal values are \({\color{darkred}\delta}^*_i=1\), and then resolve. 

In practice, I always use the following trick for this:

   solve m minimizing area using miqcp;
   delta.fx(i)$(delta.l(i)>0.5) = 0;
   solve m minimizing area using miqcp;


This is because tolerances are involved: binary variables can assume values like 0.9999 or 0.0001. Instead of introducing my own tolerance, I just use 0.5 as threshold.

This shows the following two solutions:


The area of the next-best solution is significantly larger than the original, best solution. This may explain why the models are fairly easy to solve.  


GAMS Model

$onText

 

   Find smallest axis-aligned rectangle containing k out of n points.

    

   References:

      Timothy M. Chan, Sariel Har-Peled

      Smallest k-Enclosing Rectangle Revisited

      https://arxiv.org/abs/1903.06785

 

$offtext

 

*---------------------------------------------------------------

* options

*---------------------------------------------------------------

 

* 0: no plot

* 1: produce HTML plot)

$set htmlplot 0

 

option miqcp=antigone;

 

option seed=12345;

 

*---------------------------------------------------------------

* data

*---------------------------------------------------------------

 

scalars

   maxsize 'size of the outer box' / 10 /

   k       'number of points to select' / 25 /

;

 

sets

   i   'points' /point1*point100/

   c   'coordinates' /x,y/

   lu  'lower/upper values' /min,max/

;

 

parameter p(i,c) 'points';

p(i,c) = uniform(0,maxsize);

 

display k,p;

 

*---------------------------------------------------------------

* derived data

*---------------------------------------------------------------

 

parameter box(c,lu) 'min/max values';

box(c,'min') = smin(i, p(i,c));

box(c,'max') = smax(i, p(i,c));

 

display box;

 

*---------------------------------------------------------------

* model: enclose k points

*---------------------------------------------------------------

 

binary variable delta(i) 'point is selected';

 

variable

   rect(c,lu) 'rectangle to contain k points'

   d(c)       'length of sides'

   area       'size of rectangle'   

;

 

* bounds

rect.lo(c,lu) = box(c,'min');

rect.up(c,lu) = box(c,'max');

 

d.lo(c) = 0;

d.up(c) = box(c,'max') - box(c,'min');

 

 

equations

    obj      'objective: area of enclosing rectangle'

    count    'number of points to select'

    inside_min(i,c)  'points are inside'

    inside_max(i,c)  'points are inside'

    calcd(c)         'calculate d'

;

 

count..  sum(i, delta(i)) =e= k;

 

* inside: delta(i)=1  ==> rect(c,'min') <= p(i,c) <= rect(c,'max')

inside_min(i,c)..  rect(c,'min') =l= delta(i)*p(i,c) + (1-delta(i))*box(c,'max');

inside_max(i,c)..  rect(c,'max') =g= delta(i)*p(i,c) + (1-delta(i))*box(c,'min');

 

calcd(c).. d(c) =e= rect(c,'max')-rect(c,'min');

 

obj.. area =e= d('x')*d('y');

 

 

model m /all/;

solve m minimizing area using miqcp;

 

display rect.l,delta.l,area.l;

 

 

*---------------------------------------------------------------------

* algorithm

*---------------------------------------------------------------------

 

scalar n;

n = card(i);

 

embeddedCode Python:

 

n = int(list(gams.get("n"))[0]);

k = int(list(gams.get("k"))[0]);

print(f"n={n},k={k}")

 

# get points and store as list of tuples (id,x,y)

points=[]

pgen = gams.get("p",keyFormat=KeyFormat.FLAT)

for i in range(n):

    px = next(pgen)

    py = next(pgen)

    points.append((i+1,px[2],py[2]))

 

# the rest was largely generated by ChatGPT

 

y_values = sorted({y for _, _, y in points})

 

best = None

 

# Try every possible bottom and top edge.

for i, bottom in enumerate(y_values):

    for top in y_values[i:]:

        # Points inside this horizontal band, sorted by x.

        in_band = sorted(

            (x, number)

            for number, x, y in points

            if bottom <= y <= top

        )

 

        if len(in_band) < k:

            continue

 

        # Any minimal-width interval containing k points can be

        # represented by k consecutive points in this sorted list.

        for start in range(len(in_band) - k + 1):

            left = in_band[start][0]

            right = in_band[start + k - 1][0]

            width = right - left

            height = top - bottom

            area = width * height

 

            if best is None or area < best["area"]:

                best = {

                    "area": area,

                    "left": left,

                    "right": right,

                    "bottom": bottom,

                    "top": top,

                    "point_numbers": [

                        number

                        for _, number in in_band[start:start + k]

                    ],

                }

 

print(f"x bounds: [{best['left']}, {best['right']}]")

print(f"y bounds: [{best['bottom']}, {best['top']}]")

print(f"width: {best['right'] - best['left']}")

print(f"height: {best['top'] - best['bottom']}")

print(f"area: {best['area']}")

print(f"selected points: {sorted(best['point_numbers'])}")

 

endEmbeddedCode

 

 

*---------------------------------------------------------------------

* visualization

*---------------------------------------------------------------------

 

$set html  plot.html

$set data  data.js

 

$if %htmlplot%==0 $goto skipplot

 

 

file fdata /%data%/; put fdata;

 

* points

put "p=["/;

loop(i,

  put "  {i:'",i.tl:0,"',x:",p(i,'x'):0:4,",y:",p(i,'y'):0:4,",delta:";

  put$(delta.l(i)>0.5) "'red'";

  put$(delta.l(i)<0.5) "'darkblue'";

  put "},"/;

);

put "]"/;

put "pnts={n:",(card(i)):0:0,

       ",map:'[0,",maxsize:0:0,"]&times;[0,",maxsize:0:0,"]'",

       ",minx:",box('x','min'):0:3,

       ",miny:",box('y','min'):0:3,

       ",maxx:",box('x','max'):0:3,

       ",maxy:",box('y','max'):0:3,

       ",avgx:",(sum(i,p(i,'x'))/card(i)):0:3,

       ",avgy:",(sum(i,p(i,'y'))/card(i)):0:3,

       ",k:",k:0:0,

       ",rxmin:",(rect.l('x','min')):0:3,

       ",rxmax:",(rect.l('x','max')):0:3,

       ",rymin:",(rect.l('y','min')):0:3,

       ",rymax:",(rect.l('y','max')):0:3,

       ",area:",(area.l):0:3,

       "}"/;

 

putclose;

 

$onecho > %html%

<html>

<script src="https://cdn.plot.ly/plotly-3.4.0.min.js" charset="utf-8"></script>

<script src="%data%" charset="utf-8"></script>

<style>

.brdr {

  border: 1px solid black;

  border-collapse: collapse;

}

th, td {

  padding-left: 4px;

  padding-right: 4px;

}

</style>

 

<h1>Data</h1>

<table><tr><td>

<table class="brdr">

<tr><th colspan=2 class="brdr">points</th></tr>

<tr><td class="brdr">&#x1D45B;</td><td class="brdr" id="d_n"></td></tr>

<tr><td class="brdr">map</td><td class="brdr" id="d_map"></td></tr>

<tr><td class="brdr">min &#x1D465;</td><td class="brdr" id="d_minx"></td></tr>

<tr><td class="brdr">max &#x1D465;</td><td class="brdr" id="d_maxx"></td></tr>

<tr><td class="brdr">avg &#x1D465;</td><td class="brdr" id="d_avgx"></td></tr>

<tr><td class="brdr">min &#x1D466;</td><td class="brdr" id="d_miny"></td></tr>

<tr><td class="brdr">max &#x1D466;</td><td class="brdr" id="d_maxy"></td></tr>

<tr><td class="brdr">avg &#x1D466;</td><td class="brdr" id="d_avgy"></td></tr>

</table>

</td><td>

<div id="plotDiv0" style="width: 800px; height: 600px;"></div>

</td></tr></table>

 

 

<h1>Smallest Rectangle</h1>

<table><tr><td>

<table class="brdr">

<tr><th colspan=2 class="brdr">solution</th></tr>

<tr><td class="brdr">&#x1D45B;</td><td class="brdr" id="d_n2"></td></tr>

<tr><td class="brdr">&#x1D458;</td><td class="brdr" id="d_k"></td></tr>

<tr><td class="brdr">&#x1D465;min</td><td class="brdr" id="xmin">0</td></tr>

<tr><td class="brdr">&#x1D465;max</td><td class="brdr" id="xmax">0</td></tr>

<tr><td class="brdr">&#x1D466;min</td><td class="brdr" id="ymin">0</td></tr>

<tr><td class="brdr">&#x1D466;max</td><td class="brdr" id="ymax">0</td></tr>

<tr><td class="brdr">area</td><td class="brdr" id="area">0</td></tr>

</table>

</td><td>

<div id="plotDiv1" style="width: 800px; height: 600px;"></div>

</td></tr></table>

 

 

<script>

 

document.getElementById("d_n").innerHTML = pnts['n'];

document.getElementById("d_map").innerHTML = pnts['map'];

document.getElementById("d_minx").innerHTML = pnts['minx'];

document.getElementById("d_maxx").innerHTML = pnts['maxx'];

document.getElementById("d_miny").innerHTML = pnts['miny'];

document.getElementById("d_maxy").innerHTML = pnts['maxy'];

document.getElementById("d_avgx").innerHTML = pnts['avgx'];

document.getElementById("d_avgy").innerHTML = pnts['avgy'];

 

document.getElementById("d_n2").innerHTML = pnts['n'];

document.getElementById("d_k").innerHTML = pnts['k'];

document.getElementById("xmin").innerHTML = pnts['rxmin'];

document.getElementById("xmax").innerHTML = pnts['rxmax'];

document.getElementById("ymin").innerHTML = pnts['rymin'];

document.getElementById("ymax").innerHTML = pnts['rymax'];

document.getElementById("area").innerHTML = pnts['area'];

 

 

 

px = p.map((x) => x['x'])

py = p.map((x) => x['y'])

delta = p.map((x) => x['delta'])

 

var data1 = {

  x: px,

  y: py,

  mode: 'markers',

  type: 'scatter',

  name: 'data points',

  color: 'darkblue'

 

};

 

var layout1 = {

  autosize: false,

  width: 650,

  height: 600,

  showlegend: false,

  title: {text:"Data Points"},

}

 

 

var data2 = {

  x: px,

  y: py,

  mode: 'markers',

  type: 'scatter',

  name: 'solution',

  marker: { color: delta },

  showlegend: false,

 

};

 

var layout2 = {

  autosize: false,

  width: 650,

  height: 600,

  showlegend: false,

  title: {text:"Solution"},

  shapes: [{type:'rect',

            xref:'x',

            yref:'y',

            x0:pnts['rxmin'],

            y0:pnts['rymin'],

            x1:pnts['rxmax'],

            y1:pnts['rymax'],

            line:{color:'darkred',width:2,opacity:0.5},

            fillcolor:'rgba(0, 255, 255, 0.5)',

            },]

}

 

 

var trc1 = [data1];

var options1 = {staticPlot: true, displayModeBar: false, responsive: false};

Plotly.newPlot('plotDiv0', trc1, layout1, options1);

 

var trc2 = [data2];

var options2 = {staticPlot: true, displayModeBar: false, responsive: false};

Plotly.newPlot('plotDiv1', trc2, layout2, options2);

 

 

</script>

</html>

$offecho

 

executetool 'win32.ShellExecute "%html%"';

 

$label skipplot

  

 



Variant 2: 3d weight constrained model


This section presents an extension of this model. We make the points three-dimensional, so we have \((x,y,z)\) coordinates. Also, we attach a weight \(w_i\) to each point. We want to select points such that the smallest box containing them has the smallest possible volume, and the total weight inside the box is at least a given value.

Dataset 
----     42 PARAMETER minWeight            =       25.000  required total weight

----     42 PARAMETER p  points

                   x           y           z           w

point1    4.97311132  9.63862756  8.73922492  0.48694222
point2    8.44273500  1.91226253  4.94526579  1.97392086
point3    2.34553216  3.80983451  6.00887463  1.88604759
point4    1.51476321  5.07716895  4.63323882  2.84102226
point5    3.59064347  1.52625715  3.66418827  0.28598810
point6    8.26045475  0.23893508  0.59261597  0.74594589
point7    2.84759145  1.14012489  8.90454986  1.43031631
point8    3.74340300  2.34069434  1.47696552  2.91848795
point9    0.19119162  0.78401936  9.27971966  0.24722803
point10   5.55039788  6.16361826  3.60477909  1.54572707
point11   5.21416943  4.32016162  7.49929678  1.44839493
point12   8.17376279  0.16037081  3.72949918  0.35158607
point13   6.97091674  7.66895295  8.87457833  1.40797918
point14   8.11673997  7.89650246  8.03634493  0.01746947
point15   0.16576575  3.93244606  9.60138480  2.82216428
point16   7.21520825  8.14734411  8.37501191  2.86088036
point17   7.74541956  5.28258128  6.55663865  2.15518773
point18   3.92384422  4.85091115  4.02642582  0.49230636
point19   4.34623796  1.93472329  1.28959803  2.66902215
point20   4.87601902  7.91163148  6.85615242  2.64424362
point21   3.02621818  0.00063517  4.46278421  1.06793976
point22   8.14820541  2.08022534  3.35953314  0.37095954
point23   2.81600217  1.56499252  6.20283636  2.20711958
point24   8.59944445  6.15258941  3.85719583  2.47151690
point25   7.63343999  7.31506452  0.79832230  0.09003760
point26   0.82942221  0.19238302  2.45380498  0.98465450
point27   3.34088869  3.03838803  8.34902953  2.26614420
point28   3.86320590  0.53110352  5.05611011  0.48269769
point29   1.91904761  2.34489699  9.34705030  2.36613572
point30   7.86066922  9.43453787  1.71763581  1.12082503
point31   1.81723665  2.02565095  5.32342079  2.22454539
point32   5.06474382  2.10786534  6.38173925  0.49201408
point33   2.95917598  5.53252693  0.34180259  1.39078968
point34   8.07512765  3.55101998  9.88688936  0.21046279
point35   8.13607683  7.15280359  2.43498983  1.08584045
point36   5.41651376  5.01929436  1.10750815  0.32816410
point37   1.68573148  2.86927244  8.41569886  0.70203212
point38   8.06049381  5.63100458  1.57040011  1.03604428
point39   2.71516744  0.33909707  7.18412193  0.67756099
point40   7.05045720  1.57713422  7.52170649  2.86265273
point41   9.55352881  7.85491794  2.56555975  1.67214859
point42   0.88005074  7.99748484  2.19292948  0.49649330
point43   9.55568192  8.03379715  4.59853976  2.83082090
point44   2.02151223  8.95713092  5.24080749  1.85253269
point45   8.71950503  9.06200769  3.51902443  2.85107000
point46   3.22910422  9.24501444  6.34123137  0.55814034
point47   4.14098498  2.68712111  2.81651324  1.00289323
point48   7.41232701  6.57107119  4.55851885  0.38501468
point49   3.80600433  4.96108423  8.99943277  2.99576587
point50   5.32278562  0.60195961  3.95965993  2.38442606
point51   4.30151391  9.59964284  2.71652476  2.04875708
point52   8.77681007  1.87229129  4.95157553  2.26472630
point53   3.73429995  4.27888100  9.51936360  1.21839747
point54   5.11992531  4.61819146  4.18987215  0.81940347
point55   0.91512513  9.23192650  9.52838375  0.04029898
point56   6.38224088  5.35277514  8.98266937  0.04989643
point57   0.46190997  4.24850010  8.16439743  2.37151490
point58   4.83398475  8.05382905  3.06655684  2.88987083
point59   8.59419222  7.31792459  5.60900115  0.02396737
point60   3.90522405  2.74138705  4.44598243  2.18304321
point61   9.81040377  2.77587450  7.98768312  1.13544649
point62   9.73693294  5.81026191  7.96489037  2.16547194
point63   8.28671680  7.17813456  4.33880861  2.37657744
point64   2.03345798  2.70771867  6.64833521  2.48492345
point65   1.10849014  8.46666734  2.17469964  1.44539345
point66   2.35013718  8.99662815  8.51057738  0.33858711
point67   8.92751427  2.23560422  2.24591679  0.84037876
point68   8.66033570  8.13866369  3.55150241  2.57050056
point69   2.38252571  6.77980742  0.04611403  2.81595325
point70   4.86160682  6.21125032  0.99224055  2.76802585
point71   6.95409069  1.54137047  4.72319908  0.51330419
point72   1.48556645  0.52138555  2.59667923  2.03764570
point73   5.23274123  5.67431764  0.51242147  2.65722595
point74   0.46637606  7.84666623  3.42950908  1.99782248
point75   9.12075148  6.62127069  3.16943360  1.22121036
point76   9.32398470  8.61141502  5.43944178  1.40091651
point77   9.90875460  8.64934732  5.54423034  0.04479335
point78   7.61158446  3.40928578  9.43670811  0.69092216
point79   8.94830907  0.49443490  4.59025930  0.04885652
point80   5.25949368  0.48432727  3.50423564  0.43704578
point81   9.71403687  5.77860854  5.44716786  2.60240596
point82   2.39148403  8.53557465  1.92326606  1.99617675
point83   9.86191789  3.84063928  1.07906650  2.81705489
point84   8.84041273  6.22256973  2.61438727  1.10294490
point85   3.15780224  9.94484490  1.77369425  2.61049940
point86   7.99741669  0.49593827  3.30467894  1.84808825
point87   5.49767336  9.31596887  8.18561520  0.88339756
point88   2.72818944  9.25260567  2.30602627  2.23725467
point89   8.25085350  2.50810826  1.73096670  0.26801162
point90   1.94732458  3.49788333  0.78445261  1.61420166
point91   6.78641732  7.26785086  7.38760084  1.02966832
point92   3.33428691  6.12756745  7.54948123  1.35768156
point93   1.80636836  2.54101527  0.33444540  1.09511773
point94   3.69935940  8.37069102  6.12492871  1.79299763
point95   0.67172074  0.18046144  7.82525390  1.76851896
point96   0.39180678  2.78360137  2.23733221  0.43721792
point97   4.74344419  8.89007749  1.08722063  1.24374664
point98   5.05399142  8.07553458  8.84330606  2.31860623
point99   4.60749515  0.89989148  2.21880935  2.98406425
point100  4.21063606  5.03441875  0.29718500  1.07116102

The model only needs minor changes.

 
3d weight constrained model
\[\begin{align} \min\> &{\color{darkred}{\mathit{volume}}} =  {\color{darkred}d}_x \cdot{\color{darkred}d}_y\cdot{\color{darkred}d}_z  \\ & {\color{darkred}d}_c = {\color{darkred}{\mathit{sbox}}}_{c,up}- {\color{darkred}{\mathit{sbox}}}_{c,lo} & \forall c \\ & {\color{darkred}{\mathit{sbox}}}_{c,lo} \le {\color{darkred}\delta}_i {\color{darkblue}p}_{i,c} +  (1-{\color{darkred}\delta}_i) {\color{darkblue}{\mathit{box}}}_{c,up} & \forall i,c \\ & {\color{darkred}{\mathit{sbox}}}_{c,up} \ge  {\color{darkred}\delta}_i {\color{darkblue}p}_{i,c} + (1-{\color{darkred}\delta}_i){\color{darkblue}{\mathit{box}}}_{c,lo}& \forall i,c\\ &  \sum_i{\color{darkblue}p}_{i,w} {\color{darkred}\delta}_i \ge {\color{darkblue}{\mathit{minWeight}}} \\& {\color{darkred}{\mathit{sbox}}}_{c,lo},{\color{darkred}{\mathit{sbox}}}_{c,up} \in [{\color{darkblue}{\mathit{box}}}_{c,lo},{\color{darkblue}{\mathit{box}}}_{c,up}] \\ & {\color{darkred}d}_c \in [0,{\color{darkblue}{\mathit{box}}}_{c,up}-{\color{darkblue}{\mathit{box}}}_{c,lo}] \\ & {\color{darkred}\delta}_i \in \{0,1\}\\ &{\color{darkblue}{\mathit{box}}}_{c,lo} = \min_i{\color{darkblue}p}_{i,c} \\ &{\color{darkblue}{\mathit{box}}}_{c,up} = \max_i{\color{darkblue}p}_{i,c}\\ & c \in \{x,y,z\} \end{align}\]

This is no longer a quadratic model: the objective is cubic. The solution looks like:

----    102 VARIABLE sbox.L  box to contain selected points

          min         max

x       8.287       9.714
y       1.872       9.062
z       2.566       5.447


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

point2  1.000,    point24 1.000,    point41 1.000,    point43 1.000,    point45 1.000,    point52 1.000,    point63 1.000
point68 1.000,    point75 1.000,    point76 1.000,    point81 1.000,    point84 1.000


----    102 VARIABLE volume.L              =       29.571  size of box
            PARAMETER count                =       12.000  number of selected points
            PARAMETER totw                 =       25.339  total weight of selected points


The solver selected 12 points, with a total weight a little bit over 25.

The global solver has no problems with this model:


-------------------------------------------------------------------------------
Time (s) Nodes explored Nodes remaining Best possible   Best found Relative Gap
-------------------------------------------------------------------------------
 
       0              1               1    +0.000e+00   +1.127e+02           --
       4             41               0    +2.957e+01   +2.957e+01   +1.000e-04
 
-------------------------------------------------------------------------------
Termination Status : Global minimum
Best Feasible Point: +2.957114e+01
Best Possible Point: +2.956818e+01
       Relative Gap: +1.000000e-04


The optimization model for this problem is quite simple and straightforward, and a direct extension of the previous model. In contrast, it is possible, but not as easy, to create an enumeration scheme for this problem. 


GAMS Model

$onText

 

 

   Find smallest axis-aligned 3d containing box containing

   a subset of given 3d points such that the total weight

   of the points inside the box is larger than a given number.

    

 

$offText

 

 

*---------------------------------------------------------------

* options

*---------------------------------------------------------------

 

option minlp=antigone;

 

option seed=12345;

 

*---------------------------------------------------------------

* data

*---------------------------------------------------------------

 

scalars

   maxsize 'size of the outer box' / 10 /

   minWeight  'required total weight' / 25 /

;

 

sets

   i    'points' /point1*point100/

   a    'attributes (coordinates,weight)' /x,y,z,w/

   c(a) 'coordinates' /x,y,z/

   lu   'lower/upper values' /min,max/

;

 

parameters p(i,a) 'points';

     

p(i,c) = uniform(0,maxsize);

p(i,'w') = uniform(0,3);

display minWeight,p;

 

 

*---------------------------------------------------------------

* derived data

*---------------------------------------------------------------

 

parameter box(c,lu) 'outer box: min/max values';

box(c,'min') = smin(i, p(i,c));

box(c,'max') = smax(i, p(i,c));

display box;

 

*---------------------------------------------------------------

* model: enclose k points

*---------------------------------------------------------------

 

binary variable delta(i) 'point is selected';

 

variable

   sbox(c,lu)  'box to contain selected points'

   d(c)       'length of sides'

   volume     'size of box'   

;

 

* bounds

sbox.lo(c,lu) = box(c,'min');

sbox.up(c,lu) = box(c,'max');

 

d.lo(c) = 0;

d.up(c) = box(c,'max') - box(c,'min');

 

 

equations

    obj      'objective: area of enclosing rectangle'

    sumweight   'weight restriction'

    inside_min(i,c)  'points are inside'

    inside_max(i,c)  'points are inside'

    calcd(c)         'calculate d'

;

 

sumweight..  sum(i, p(i,'w')*delta(i)) =g= minWeight;

 

* inside: delta(i)=1  ==> sbox(c,'min') <= p(i,c) <= sbox(c,'max')

inside_min(i,c)..  sbox(c,'min') =l= delta(i)*p(i,c) + (1-delta(i))*box(c,'max');

inside_max(i,c)..  sbox(c,'max') =g= delta(i)*p(i,c) + (1-delta(i))*box(c,'min');

 

calcd(c).. d(c) =e= sbox(c,'max')-sbox(c,'min');

 

obj.. volume =e= d('x')*d('y')*d('z');

 

model m /all/;

solve m minimizing volume using minlp;

 

scalars

  count 'number of selected points'

  totw  'total weight of selected points' 

;

count = sum(i$(delta.l(i)>0.5),1);

totw = sum(i$(delta.l(i)>0.5),p(i,'w'));

 

display sbox.l,delta.l,volume.l,count,totw;

 

 



Variant 3: max number of points


We can turn our model around and instead of:

    minimize area
    subject to number of points in rectangle is given

we can do:

    maximize number of points in rectangle
    subject to area of rectangle is given


 
Max Covered Points model
\[\begin{align} \max\> &{\color{darkred}{k}}=\sum_i{\color{darkred}\delta}_i \\ & {\color{darkblue}{\mathit{area}}} =  {\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} & \forall c \\ & {\color{darkred}{\mathit{rect}}}_{c,lo} \le {\color{darkred}\delta}_i {\color{darkblue}p}_{i,c} +  (1-{\color{darkred}\delta}_i) {\color{darkblue}{\mathit{box}}}_{c,up} & \forall i,c \\ & {\color{darkred}{\mathit{rect}}}_{c,up} \ge  {\color{darkred}\delta}_i {\color{darkblue}p}_{i,c} + (1-{\color{darkred}\delta}_i){\color{darkblue}{\mathit{box}}}_{c,lo}& \forall i,c\\& {\color{darkred}{\mathit{rect}}}_{c,lo},{\color{darkred}{\mathit{rect}}}_{c,up} \in [{\color{darkblue}{\mathit{box}}}_{c,lo},{\color{darkblue}{\mathit{box}}}_{c,up}] \\ & {\color{darkred}d}_c \in [0,{\color{darkblue}{\mathit{box}}}_{c,up}-{\color{darkblue}{\mathit{box}}}_{c,lo}] \\ & {\color{darkred}\delta}_i \in \{0,1\} \\ &{\color{darkblue}{\mathit{box}}}_{c,lo} = \min_i{\color{darkblue}p}_{i,c}\\ &{\color{darkblue}{\mathit{box}}}_{c,up} = \max_i{\color{darkblue}p}_{i,c}\\ & c \in \{x,y\} \end{align}\]


Global solvers have no problem with the non-convex quadratic objective becoming a constraint.

If we set the size of the problem again to \({\color{darkblue}n}=100\) and further \({\color{darkblue}{\mathit{area}}}=10\), we see:



I.e., we can find a rectangle with a given area of 10 that contains 26 points. In general, the solution is not unique (we may move the rectangle around a little bit).

Note that we no longer have the property that data points are on the border (all sides) of the rectangle. The enumeration scheme we used before is no longer directly applicable.  

Solver issues


  • Somehow Baron gives up on these models with the message: 

    A potentially catastrophic access violation just took place

    I have passed on the model to the GAMS people for further investigation. I don't have a workaround for this.
  • Xpress gives for the second-best model a slightly non-optimal solution:



    Other solvers find a second solution of 12.99. This may be a tolerance issue, although Xpress thinks this is globally optimal:

    ... Status: Normal (1) / Global Optimal (1)
    ... Objective Value:  1.305122276856e+01
    ... Objective Bound:  1.305122276856e+01
    ... Optimality Gap:   0.000000000000e+00 (0.00%)

    After changing the optimality tolerance from 1e-6 to 1e-7, Xpress delivers the solution of 12.99.

    ... Status: Normal (1) / Global Optimal (1)
    ... Objective Value:  1.299035980013e+01
    ... Objective Bound:  1.299035980013e+01
    ... Optimality Gap:   0.000000000000e+00 (0.00%)

    In general, I don't like changing solver tolerances, but this was a useful exercise. It confirms that we are dealing with a tolerance issue.

Conclusion


This small, simple model lets us find global, proven solutions. Global MINLP solvers do a really good job on this. Global MINLP solvers are readily available. Examples are: Antigone, Baron, Xpress, Gurobi and SCIP. An alternative is to use a clever enumeration scheme. Of course, we can easily make the problem (much) more difficult by increasing the size.

The enumeration algorithm works very well for these small data sets. However, it is not easy to adapt to model variants. The optimization model is more flexible: with just some minor updates, we can largely reuse the base model.

References


  1. Eppstein, D., Erickson, J. Iterated nearest neighbors and finding minimal polytopes. Discrete Comput Geom 11, 321–350 (1994). https://doi.org/10.1007/BF02574012
  2. Timothy M. Chan and Sariel Har-Peled, Smallest k-Enclosing Rectangle Revisited, 2019, https://arxiv.org/abs/1903.06785.
at

No comments:

Post a Comment