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. I would argue that in production code, loops should always be protected (simply by an iteration count, or maybe a timer).

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. 

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 where instead of the more traditional Euclidean distances we use Manhattan of \(L1\) distances. Again, we can choose area- or perimeter-based objectives. It is interesting to see the difference.

A high-level area focused model where we want to cover \(n\) points by \(k\) rectangles can look like:

 

High-level model

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

Here \[{\color{darkred}\delta}_{i,k} = \begin{cases}1&\text{if point $i$ is assigned to rectangle $k$} \\ 0 & \text{otherwise}\end{cases}\] Each point is assigned to exactly one rectangle. The \(\color{darkred}{\mathit{rect}}\) variables indicate the location, height and width of each rectangle. We implicitly require here that \({\color{darkred}{\mathit{rect}}}_{k,c,up}\ge {\color{darkred}{\mathit{rect}}}_{k,c,lo}\).

If we want to optimize the perimeters of the rectangles, we need to replace the quadratic objective by a linear one: \[ \min\> {\color{darkred}{\mathit{perimeter}}} = 2 \sum_{k,c} \left  ({\color{darkred}{\mathit{rect}}}_{k,c,up}-{\color{darkred}{\mathit{rect}}}_{k,c,lo}\right )\] 

The actual implemented model is:

 

Implemented model

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

The perimeter objective would be: \[\min\>\sum_k 2\left({\color{darkred}d}_{k,x} +{\color{darkred}d}_{k,y} \right)\]We can add to the model a simple symmetry breaker: \[ {\color{darkred}{\mathit{rect}}}_{k,x,lo} \le {\color{darkred}{\mathit{rect}}}_{k+1,x,lo}\] This constraint orders the rectangles by their \(x\)-coordinate. It makes the solution more readable, but more importantly, it reduces symmetry in the model and can increase solver performance, sometimes by a large amount. 

Data

I generated a small dataset for this model. The points are drawn from a uniform distribution, so the clustering shouldn't be so obvious.

Dataset

----     37 SET k  rectangles

rectangle1,    rectangle2,    rectangle3,    rectangle4,    rectangle5


----     37 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

 

Let's see how confused the models are about this dataset.

Solutions

Run 1: minimize area

The first run is to "minimize area". I used the non-convex MIQCP model for this. If we look at the \({\color{darkred}\delta}_{i,k}\) variables we see:

----    113 PARAMETER assign  assignment=round(delta)

         rectangle1  rectangle2  rectangle3  rectangle4  rectangle5

point1                                                1
point2                                                            1
point3            1
point4            1
point5                                    1
point6                                                1
point7                                    1
point8                        1
point9            1
point10           1
point11                                                           1
point12           1
point13           1
point14                                                           1
point15                                                           1
point16                                               1
point17                                                           1
point18           1
point19                                                           1
point20                                                           1
point21                                                           1
point22           1
point23                                                           1
point24                                                           1
point25                                                           1
point26                                                           1
point27                                               1
point28                                   1
point29           1
point30                                                           1
point31           1
point32                       1
point33           1
point34           1
point35                                                           1
point36                                                           1
point37                                                           1
point38           1
point39           1
point40           1
point41                                                           1
point42           1
point43           1
point44                                                           1
point45                                   1
point46           1
point47                                               1
point48           1
point49           1
point50                       1
count            20           3           4           5          18

We have two large clusters (18 and 20 points), and three small ones (3,4 and 5 points). I rounded the solution, because MIP solvers often deliver solutions with very small values, such as 1e-10. Here we filtered those out. In addition we added the counts.

If we look at the rectangles, we see:

----    141 PARAMETER results  by cluster

INDEX 1 = min area obj

                     min         max         len      points        area   perimeter

rectangle1.x       0.160       3.341       3.181
rectangle1.y       0.001       6.382       6.381
rectangle1.                                           20.000      20.295      19.123
rectangle2.x       0.342       4.463       4.121
rectangle2.y       8.075       8.260       0.185
rectangle2.                                            3.000       0.764       8.613
rectangle3.x       3.591       9.435       5.844
rectangle3.y       1.515       1.935       0.420
rectangle3.                                            4.000       2.454      12.528
rectangle4.x       4.851       5.323       0.473
rectangle4.y       4.026       9.639       5.612
rectangle4.                                            5.000       2.652      12.169
rectangle5.x       6.153       9.601       3.449
rectangle5.y       3.605       8.599       4.995
rectangle5.                                           18.000      17.226      16.887
total     .                                           50.000      43.391      69.320

Obviously, the clusters containing many points, are also the largest ones in size (either area or perimeter). 

The highlighted number is the one the model tries to minimize.

As a side note: I spent a little bit of time on this table. In general, for more complex models, with more complex results, I think it is always a good idea to think a bit about reporting. Reports, in essence, are the model's user interface (GUI). So they deserve some attention.  

A picture of the results is:

This looks quite reasonable.

Run 2: minimize perimeter

The "minimize perimeter" is an easier, linear MIP model. Let's see what happens when we run this version. I expect we will observe a solution close to the one in the previous section.

The optimal assignment is as follows:

----    133 PARAMETER assign  assignment=round(delta)

         rectangle1  rectangle2  rectangle3  rectangle4  rectangle5

point1                                                1
point2            1
point3            1
point4            1
point5            1
point6            1
point7            1
point8            1
point9                        1
point10           1
point11           1
point12           1
point13           1
point14           1
point15           1
point16           1
point17           1
point18           1
point19           1
point20           1
point21           1
point22           1
point23                                                           1
point24           1
point25           1
point26           1
point27           1
point28           1
point29           1
point30           1
point31                                   1
point32           1
point33           1
point34           1
point35           1
point36           1
point37           1
point38           1
point39           1
point40           1
point41           1
point42           1
point43           1
point44           1
point45           1
point46           1
point47           1
point48           1
point49           1
point50           1
count            46           1           1           1           1

The details about the clusters are:

----    148 PARAMETER results  by cluster

INDEX 1 = min perim obj

                     min         max         len      points        area   perimeter

rectangle1.x       0.160       9.435       9.274
rectangle1.y       0.784       8.599       7.815
rectangle1.                                           46.000      72.482      34.179
rectangle2.x       0.239       0.239
rectangle2.y       0.593       0.593
rectangle2.                                            1.000
rectangle3.x       3.026       3.026 
rectangle3.y       0.001       0.001 
rectangle3.                                            1.000 
rectangle4.x       4.973       4.973 
rectangle4.y       9.639       9.639 
rectangle4.                                            1.000 
rectangle5.x       9.601       9.601
rectangle5.y       7.215       7.215
rectangle5.                                            1.000
total     .                                           50.000      72.482      34.179

This looks very different: we have one big cluster, and then some singletons. 

Why is this happening? In the "minimize area" picture, we see three long, skinny clusters. These clusters have a relatively small area but a large perimeter. So a "minimize perimeter" approach prefers short, stocky rectangles closer to squares.

Conclusion

This is a nice example where a linear approximation of a quadratic problem gives very different results. The quadratic model, which minimizes the sum of the areas, identifies the clusters in the data nicely. In contrast, the linear model uses the perimeter and essentially says there are no clusters.

In the literature, we see the linear approach mentioned often, while minimizing areas is rarely discussed. Maybe that is not warranted.