## Wednesday, January 20, 2021

### Continuous max sum rectangle: MIQP vs GA

This is a follow-up on the post on the max sum submatrix problem [1].  There we looked at the problem of finding a contiguous submatrix that maximizes the sum of the values in this submatrix.

A little bit more complicated is the following continuous version of this problem.

Assume we have $$n$$ points with an $$x$$- and $$y$$-coordinate and a value. Find the rectangle such that the sum of the values of all points inside is maximized.

#### Data set

I generated a random data set with $$n=100$$ points:

----     10 PARAMETER p  points

x           y       value

i1        17.175       5.141       8.655
i2        84.327       0.601      -3.025
i3        55.038      40.123      -9.834
i4        30.114      51.988       8.977
i5        29.221      62.888       1.438
i6        22.405      22.575      -3.327
i7        34.983      39.612       9.675
i8        85.627      27.601       5.329
i9         6.711      15.237      -7.798
i10       50.021      93.632       9.896
i11       99.812      42.266       1.606
i12       57.873      13.466      -6.672
i13       99.113      38.606       2.867
i14       76.225      37.463      -3.114
i15       13.069      26.848       8.247
i16       63.972      94.837       8.001
i17       15.952      18.894      -9.675
i18       25.008      29.751      -2.627
i19       66.893       7.455       3.288
i20       43.536      40.135       1.868
i21       35.970      10.169      -9.309
i22       35.144      38.389       6.836
i23       13.149      32.409       8.642
i24       15.010      19.213       0.159
i25       58.911      11.237      -4.008
i26       83.089      59.656      -0.068
i27       23.082      51.145      -9.101
i28       66.573       4.507       5.474
i29       77.586      78.310       0.659
i30       30.366      94.575       4.935
i31       11.049      59.646       4.401
i32       50.238      60.734       2.632
i33       16.017      36.251      -7.702
i34       87.246      59.407       9.423
i35       26.511      67.985       4.135
i36       28.581      50.659       9.725
i37       59.396      15.925       7.096
i38       72.272      65.689       2.429
i39       62.825      52.388       4.026
i40       46.380      12.440       4.018
i41       41.331      98.672       5.814
i42       11.770      22.812       2.204
i43       31.421      67.565      -8.914
i44        4.655      77.678      -0.296
i45       33.855      93.245      -8.949
i46       18.210      20.124       3.972
i47       64.573      29.714      -6.104
i48       56.075      19.723      -5.479
i49       76.996      24.635       6.273
i50       29.781      64.648       9.835
i51       66.111      73.497       5.013
i52       75.582       8.544       4.367
i53       62.745      15.035      -9.988
i54       28.386      43.419      -4.723
i55        8.642      18.694       6.476
i56       10.251      69.269       6.391
i57       64.125      76.297       7.208
i58       54.531      15.481      -5.746
i59        3.152      38.938      -0.864
i60       79.236      69.543      -9.233
i61        7.277      84.581      -3.540
i62       17.566      61.272      -1.202
i63       52.563      97.597      -3.693
i64       75.021       2.689      -7.305
i65       17.812      18.745       6.219
i66        3.414       8.712      -1.664
i67       58.513      54.040      -7.164
i68       62.123      12.686      -0.689
i69       38.936      73.400      -4.340
i70       35.871      11.323       7.914
i71       24.303      48.835      -8.712
i72       24.642      79.560      -1.708
i73       13.050      49.205      -3.168
i74       93.345      53.356      -0.634
i75       37.994       1.062       2.853
i76       78.340      54.387       2.872
i77       30.003      45.113      -3.248
i78       12.548      97.533      -7.984
i79       74.887      18.385       8.117
i80        6.923      16.353      -5.653
i81       20.202       2.463       8.377
i82        0.507      17.782      -0.965
i83       26.961       6.132      -8.201
i84       49.985       1.664      -2.516
i85       15.129      83.565      -1.700
i86       17.417      60.166      -1.916
i87       33.064       2.702      -7.767
i88       31.691      19.609       5.023
i89       32.209      95.071       6.068
i90       96.398      33.554      -9.527
i91       99.360      59.426      -0.382
i92       36.990      25.919      -4.428
i93       37.289      64.063       8.032
i94       77.198      15.525      -9.648
i95       39.668      46.002       3.621
i96       91.310      39.334       9.018
i97       11.958      80.546       8.004
i98       73.548      54.099       7.976
i99        5.542      39.072       7.489
i100      57.630      55.782      -2.180


#### High-level model

Let's denote our points by: $$\color{darkblue}p_{i,a}$$, where $$a=\{x,y,{\mathit{value}}\}$$. Our decision variables are the corner points of the rectangle: $$\color{darkred}r_{c,q}$$ where $$c=\{x,y\}$$ and $$q=\{min,max\}$$. These variables are continuous between $$\color{darkblue}L=0$$ and  $$\color{darkblue}U=100$$.

High-level model
\begin{align}\max & \sum_{i|\color{darkred}r(c,min) \le \color{darkblue}p(i,c)\le \color{darkred}r(c,max)} \color{darkblue}p_{i,{\mathit{value}}}\\ & \color{darkred}r_{c,min} \le \color{darkred}r_{c,max} \\ &\color{darkred}r_{c,q}\in [L,U]\end{align}

This looks easy. Well, not so fast...

#### Development

The first thing to do is to introduce binary variables that indicate if a point $$i$$ is in between $$\color{darkred}r_{c,min}$$ and $$\color{darkred}r_{c,max}$$ for both $$c \in \{x,y\}$$. It is actually simpler to consider three cases:

The $$x$$-coordinate of a given point $$\color{darkblue}p_{i,x}$$ is in one of the three segments denoted by $$\color{darkred}\delta_{i,x,k}$$. Of course, similar for the $$y$$-coordinate. Note that this is a little bit different than usual. Normally we have a point which is a decision variable and the limits for the segments that are constants. Here we have a fixed data point, but the segment boundaries $$\color{darkred}r_{x,min}, \color{darkred}r_{x,max}$$ are variable.

So, looking at the picture, we need to implement: \begin{align} &\color{darkred}\delta_{i,x,1}=1 \Rightarrow \color{darkblue}L \le \color{darkblue} p_{i,x} \le \color{darkred}r_{x,min} \\ & \color{darkred}\delta_{i,x,2}=1 \Rightarrow \color{darkred}r_{x,min} \le \color{darkblue}p_{i,x} \le \color{darkred}r_{x,max} \\ &\color{darkred}\delta_{i,x,3}=1\Rightarrow \color{darkred}r_{x,max} \le \color{darkblue} p_{i,x} \le \color{darkblue}U \\ & \sum_k \color{darkred}\delta_{i,x,k}= 1\end{align} We can formulate the implications as the following inequalities: \begin{align} &\color{darkred}r_{x,min} \ge \color{darkblue}p_{i,x} \cdot \color{darkred}\delta_{i,x,1} + \color{darkblue}L \cdot (1-\color{darkred}\delta_{i,x,1}) \\ & \color{darkred}r_{x,min} \le \color{darkblue} p_{i,x} \cdot \color{darkred}\delta_{i,x,2} + \color{darkblue}U \cdot (1-\color{darkred}\delta_{i,x,2}) \\& \color{darkred}r_{x,max} \ge \color{darkblue}p_{i,x} \cdot \color{darkred}\delta_{i,x,2} + \color{darkblue}L \cdot (1-\color{darkred}\delta_{i,x,2}) \\ & \color{darkred}r_{x,max} \le \color{darkblue} p_{i,x} \cdot \color{darkred}\delta_{i,x,3} + \color{darkblue}U \cdot (1-\color{darkred}\delta_{i,x,3})\end{align}

Of course, we need to repeat this for the $$y$$ direction. Once we have all variables $$\color{darkred}\delta_{i,x,k}$$ and $$\color{darkred}\delta_{i,y,k}$$, we still need to combine them. A point $$i$$ is inside the rectangle $$\color{darkred}r$$ if and only if $$\color{darkred}\delta_{i,x,2}\cdot \color{darkred}\delta_{i,y,2} = 1$$.

With this, our model can look like:

MIQP Model A
\begin{align}\max & \sum_i \color{darkblue}p_{i,{\mathit{value}}} \cdot \color{darkred}\delta_{i,x,2} \cdot \color{darkred}\delta_{i,y,2}\\ & \color{darkred}r_{c,min} \ge \color{darkblue} p_{i,c} \cdot \color{darkred}\delta_{i,c,1} + \color{darkblue}L \cdot (1-\color{darkred}\delta_{i,c,1}) && \forall i,c\\ & \color{darkred}r_{c,{\mathit{min}}} \le \color{darkblue} p_{i,c} \cdot \color{darkred}\delta_{i,c,2} + \color{darkblue}U \cdot (1-\color{darkred}\delta_{i,c,2}) && \forall i,c\\& \color{darkred}r_{c,{\mathit{max}}} \ge \color{darkblue} p_{i,c} \cdot \color{darkred}\delta_{i,c,2} + \color{darkblue}L \cdot (1-\color{darkred}\delta_{i,c,2}) && \forall i,c \\ & \color{darkred}r_{c,{\mathit{max}}} \le \color{darkblue} p_{i,c} \cdot \color{darkred}\delta_{i,c,3} + \color{darkblue}U \cdot (1-\color{darkred}\delta_{i,c,3}) &&\forall i,c \\ & \sum_k \color{darkred}\delta_{i,c,k}=1 && \forall i,c \\ & \color{darkred}r_{c,min} \le \color{darkred}r_{c,max} && \forall c \\ & \color{darkred}\delta_{i,c,k} \in \{0,1\} \\ &\color{darkred}r_{c,q} \in [\color{darkblue}L,\color{darkblue}U] \end{align}

This model can easily be linearized using the techniques demonstrated in [1]. Here we will assume the solver is linearizing this automatically. This model solves very quickly, and the results are:

----     62 VARIABLE z.L                   =      102.234  objective

----     62 VARIABLE r.L  rectangle

min         max

x       7.277      93.345
y      15.525      97.533


The solver (Cplex in this case), linearized the problem for us and solved it as a linear MIP.  It found and proved the global optimal solution in 0 nodes (i.e. all work was done during preprocessing) and 1,249 iterations. The solution time was 2 seconds.

#### Alternative approach

In the comments, Imre Polik suggested populating a (very) sparse matrix using the values and then apply a method from [1].  The sparse matrix has one element per row and column, and the points are sorted by their $$x$$- and $$y$$-coordinate. The following small example shows how this is done:

----     17 PARAMETER pt  points

x           y       value

i1       17.175      99.812      -2.806
i2       84.327      57.873      -2.971
i3       55.038      99.113      -7.370
i4       30.114      76.225      -6.998
i5       29.221      13.069       1.782
i6       22.405      63.972       6.618
i7       34.983      15.952      -5.384
i8       85.627      25.008       3.315
i9        6.711      66.893       5.517
i10      50.021      43.536      -3.927

----     52 PARAMETER spMap  sparse matrix mapper

x           y       value

i1 .i2 .i10      17.175      99.812      -2.806
i2 .i9 .i5       84.327      57.873      -2.971
i3 .i8 .i9       55.038      99.113      -7.370
i4 .i5 .i8       30.114      76.225      -6.998
i5 .i4 .i1       29.221      13.069       1.782
i6 .i3 .i6       22.405      63.972       6.618
i7 .i6 .i2       34.983      15.952      -5.384
i8 .i10.i3       85.627      25.008       3.315
i9 .i1 .i7        6.711      66.893       5.517
i10.i7 .i4       50.021      43.536      -3.927

----     56 PARAMETER spA  sparse matrix representation

i1          i2          i3          i4          i5          i6          i7          i8          i9         i10

i1                                                                                5.517
i2                                                                                                                   -2.806
i3                                                                    6.618
i4        1.782
i5                                                                                           -6.998
i6                   -5.384
i7                                           -3.927
i8                                                                                                       -7.370
i9                                                       -2.971
i10                               3.315


The first parameter holds our points. The second parameter shows the mapping between the original record number (column one) and the location in the matrix $$(i,j)$$: the second and third index. The final parameter shows the matrix in a more standard representation. (Note: this is not so easy to do in GAMS. You could use gdxrank. I used a little bit of Python code to do the sorting.) With this we can use the model from [1]:

MIQP model B
\begin{align} \max& \sum_{i,j} \color{darkblue}a_{i,j} \cdot \color{darkred}x_i \cdot \color{darkred}y_j\\ & \color{darkred}p_i \ge \color{darkred}x_{i} - \color{darkred}x_{i-1} \\ & \color{darkred}q_j \ge \color{darkred}y_{j} - \color{darkred}y_{i-1} \\ & \sum_i \color{darkred}p_i \le 1 \\ & \sum_j \color{darkred}q_j \le 1\\ &\color{darkred}x_i, \color{darkred}y_j,\color{darkred}p_i, \color{darkred}q_j \in \{0,1\} \end{align}

We assume $$\color{darkred}x_0 = \color{darkred}y_0= 0$$. Of course, we can apply a standard linearization or let the solver linearize. Note that we only have $$n$$ quadratic terms (only for the nonzero elements). When we find the optimal solution for this model, we can retrieve the $$x$$- and $$y$$-coordinate of the first and last selected row and column to form our rectangle.

With our 100 point data set, we find the optimal solution in 0 nodes and 759 iterations. This took less than a second.

#### Genetic algorithm

Here we try to solve the problem using R's GA package. The fitness function is basically the same as our objective in the high-level model. The code (and output) can look like:

> library(GA)
>
> # data is stored in data frame
> str(df)
'data.frame':	100 obs. of  4 variables:
$point: chr "i1" "i2" "i3" "i4" ...$ x    : num  17.2 84.3 55 30.1 29.2 ...
$y : num 5.141 0.601 40.123 51.988 62.888 ...$ value: num  8.66 -3.02 -9.83 8.98 1.44 ...
>
> # fitness function
> f <- function(x) {
+   xmin <- min(x[1],x[2])
+   xmax <- max(x[1],x[2])
+   ymin <- min(x[3],x[4])
+   ymax <- max(x[3],x[4])
+   ok <- (df$x <= xmax) & (df$x >= xmin) & (df$y <= ymax) & (df$y >= ymin)