## 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 .  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 . 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 .  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 :

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,x)
+   xmax <- max(x,x)
+   ymin <- min(x,x)
+   ymax <- max(x,x)
+   ok <- (df$x <= xmax) & (df$x >= xmin) & (df$y <= ymax) & (df$y >= ymin)
+   sum(ok * df$value) + } > > # call the ga solver > system.time(result <- ga(type = "real-valued", + fitness = f, + lower = rep(0,4), upper = rep(100,4), + popSize = 100, maxiter = 500, monitor = T, + seed = 12345)) GA | iter = 1 | Mean = 8.058819 | Best = 48.175182 GA | iter = 2 | Mean = 6.413503 | Best = 48.175182 GA | iter = 3 | Mean = 8.234065 | Best = 48.175182 GA | iter = 4 | Mean = 5.209888 | Best = 48.175182 GA | iter = 5 | Mean = 7.171262 | Best = 48.175182 GA | iter = 6 | Mean = 9.767468 | Best = 48.175182 GA | iter = 7 | Mean = 11.41823 | Best = 48.17518 GA | iter = 8 | Mean = 9.89723 | Best = 48.20186 GA | iter = 9 | Mean = 14.19131 | Best = 53.02362 GA | iter = 10 | Mean = 20.83307 | Best = 53.61161 . . . GA | iter = 490 | Mean = 87.13611 | Best = 102.23405 GA | iter = 491 | Mean = 83.61628 | Best = 102.23405 GA | iter = 492 | Mean = 87.16728 | Best = 102.23405 GA | iter = 493 | Mean = 86.52578 | Best = 102.23405 GA | iter = 494 | Mean = 83.02602 | Best = 102.23405 GA | iter = 495 | Mean = 86.58533 | Best = 102.23405 GA | iter = 496 | Mean = 89.53308 | Best = 102.23405 GA | iter = 497 | Mean = 89.4938 | Best = 102.2341 GA | iter = 498 | Mean = 88.29249 | Best = 102.23405 GA | iter = 499 | Mean = 88.39676 | Best = 102.23405 GA | iter = 500 | Mean = 85.33931 | Best = 102.23405 user system elapsed 6.78 0.27 7.32  Note that I did not specify the constraints $$\color{darkred}{\mathit{xmin}} \le \color{darkred}{\mathit{xmax}}$$ and $$\color{darkred}{\mathit{ymin}} \le \color{darkred}{\mathit{ymax}}$$. Instead, when the fitness function receives the vector $$x$$ it just sorts the values. This is somewhat of a trick. But it helps in allowing us to use an unconstrained optimization problem. Additionally, we have good bounds on the four decision variables. These properties are usually a big win for heuristic solvers like this. Because we know the optimal solution, we can indeed verify that this heuristic also finds the optimal solution. A nice feature is to plot the results: It finds 5 solutions with the best objective: > result@solution x1 x2 x3 x4 [1,] 8.058577 91.76226 15.88715 95.49885 [2,] 8.072476 91.76382 15.87492 95.50330 [3,] 8.086532 91.76143 15.91187 95.52999 [4,] 8.071193 91.76518 15.88455 95.50384 [5,] 7.978493 92.05050 15.92354 95.49153  These are essentially the same. #### Larger data set When using a larger random data set with $$n=1,000$$ points, we end up with a large MIQP Model A. It has 10k rows and 6k columns. After Cplex reformulates this into a linear model, we have 11.5k rows and 7k columns. This solves to optimality in 2,000 seconds with a proven optimal objective of 209.426. ---- 63 VARIABLE z.L = 209.426 objective ---- 63 VARIABLE r.L rectangle min max x 7.045 97.298 y 2.812 58.167  MIQP Model B is smaller. Before linearization, we have 2k rows and 4k columns. This translates to 3,495 rows, 4,982 columns after (automatic) linearization and presolving. The model solves in about 100 seconds and finds and proves the global optimum of 209.426. ---- 88 VARIABLE z.L = 209.426 objective ---- 98 PARAMETER results model B results x y min 7.098 2.812 max 97.215 57.956  When I feed this into the GA solver, with a limit of 5,000 iterations, we see: -- Genetic Algorithm ------------------- GA settings: Type = real-valued Population size = 10 Number of generations = 5000 Elitism = 1 Crossover probability = 0.8 Mutation probability = 0.1 Search domain = x1 x2 x3 x4 lower 0 0 0 0 upper 100 100 100 100 GA results: Iterations = 5000 Fitness function value = 205.203 Solutions = x1 x2 x3 x4 [1,] 80.91125 7.089162 58.0342 11.90685 [2,] 80.91125 7.089162 58.0342 11.90685 [3,] 80.91125 7.089162 58.0342 11.90685 [4,] 80.91125 7.089162 58.0342 11.90685  Note that in this case, $$x_1,x_2,x_3,x_4$$ corresponds to $$xmax,xmin,ymax,ymin$$. We see that the $$x$$ part of the rectangle is close, but $$ymin=11.9$$ is a bit off compared to the optimal MIQP solution. Indeed the objective 205.203 is a little bit worse than 209.426 which we found earlier. On the other hand, GA just needed about 10 seconds to find this solution (depending on the population size). #### Conclusions • Formulating this problem as a non-convex MIQP model requires some thought (and some work). • But it can provide proven optimal solutions (or good solutions if optimality is too costly). • The sparse matrix formulation (model B) works better than the rectangle formulation (model A). • Using a Genetic Algorithm based heuristic makes the modeling quite easy. There is one trick involved so we have an unconstrained problem. After that, the GA solver finds the optimal solution fairly quickly. Of course, it does not prove optimality. We only know this because we also have a mathematical programming model that was solved before. • The GA model and the MIQP model have very little in common. It is almost always a bad idea to reuse a MIP model and feed it directly to a meta-heuristic. • When developing heuristics for a large and difficult problem, I also like to implement a mathematical programming model. This allows us to compare solutions. First, this is a good debugging aid. But also this can give us some feedback on the quality of the solutions found with the heuristic, even if only for smaller data sets. #### References 1. Submatrix with largest sum, http://yetanothermathprogrammingconsultant.blogspot.com/2021/01/submatrix-with-largest-sum.html #### Appendix A: GAMS formulation A $ontext    Given: n points with (x,y) coordinates and a value    (value can be positive or negative).    Find a rectangle such that the sum of the values of the    points inside is maximized.    This formulation works directly on the points. Each point    is in one of three regions for each coordinate:            L -------- min -------- max --------- U   region:       k1           k2           k3 $offtext *---------------------------------------------------------------- * data *---------------------------------------------------------------- set i 'number of points' /i1*i100/ a 'attributes' /x,y,value/ c(a) 'coordinates' /x,y/ k 'segment' /k1 'before' k2 'in between' k3 'after'/ q 'limit' /min,max/ ; * * select (x,y) coordinates from [L,U]x[L,U] * scalars L /0/ U /100/ ; parameter p(*,*) 'points'; p(i,'x') = uniform(L,U); p(i,'y') = uniform(L,U); p(i,'value') = uniform(-10,10); display p; *---------------------------------------------------------------- * model *---------------------------------------------------------------- variable z 'objective' r(c,q) 'rectangle' ; r.lo(c,q) = L; r.up(c,q) = U; binary variable delta(c,i,k) 'region of point i'; equations obj 'objective' e1(c,i) 'region k1' e2(c,i) 'region k2' e3(c,i) 'region k2' e4(c,i) 'region k3' one(c,i) 'point is in one region' minmax(c) 'min <= max' ; obj.. z=e= sum(i, p(i,'value')*delta('x',i,'k2')*delta('y',i,'k2')); e1(c,i).. r(c,'min') =g= p(i,c)*delta(c,i,'k1')+L*(1-delta(c,i,'k1')); e2(c,i).. r(c,'min') =l= p(i,c)*delta(c,i,'k2')+U*(1-delta(c,i,'k2')); e3(c,i).. r(c,'max') =g= p(i,c)*delta(c,i,'k2')+L*(1-delta(c,i,'k2')); e4(c,i).. r(c,'max') =l= p(i,c)*delta(c,i,'k3')+U*(1-delta(c,i,'k3')); one(c,i).. sum(k,delta(c,i,k)) =e= 1; minmax(c).. r(c,'min') =l= r(c,'max'); model mod /all/; option optcr=0,threads=8,miqcp=cplex; solve mod maximizing z using miqcp; display z.l,r.l; #### Appendix B: GAMS formulation B $ontext    Given: n points with (x,y) coordinates and a value    (value can be positive or negative).    Find a rectangle such that the sum of the values of the    points inside is maximized.    This formulation sorts the points by x and y coordinates    to form a (very) sparse matrix (one element per row and column). $offtext *------------------------------------------------------ * data *------------------------------------------------------ set i 'number of points' /i1*i100/ a 'attributes' /x,y,value/ c(a) 'coordinates' /x,y/ ; * * select (x,y) coordinates from [L,U]x[L,U] * scalars L /0/ U /100/ ; parameter pt(i,a) 'points'; pt(i,'x') = uniform(L,U); pt(i,'y') = uniform(L,U); pt(i,'value') = uniform(-10,10); display pt; *------------------------------------------------------ * create a sparse matrix representation *------------------------------------------------------ alias(i,j,k); parameter spMap(i,j,k,a) 'sparse matrix mapper'; EmbeddedCode Python: # copy pt to Python p = gams.get("pt") # index i indx = [t for t in p if t=="x"] n = len(indx) # get x,y,value x = [t for t in p if t=="x"] y = [t for t in p if t=="y"] v = [t for t in p if t=="value"] # get rank of x,y xrank = [sorted(x).index(i) for i in x] yrank = [sorted(y).index(i) for i in y] # create a parameter q that has x and y sorted q = [] for k in range(n): i0 = indx[k] i = indx[xrank[k]] j = indx[yrank[k]] q.append((i0,i,j,"x",x[k])) q.append((i0,i,j,"y",y[k])) q.append((i0,i,j,"value",v[k])) gams.set("spmap",q) endEmbeddedCode spmap option spmap:3:3:1; display$(card(spmap)<=500) spmap; parameter spA(i,j) 'sparse matrix representation'; spA(i,j) = sum(k,spmap(k,i,j,'value')); display$(card(spA)<=100) spA; *------------------------------------------------------ * model *------------------------------------------------------ binary variable x(i) 'selected rows' y(j) 'selected columns' p(i) 'start of submatrix' q(j) 'start of submatrix' ; variable z 'objective'; equation obj 'objective' pbound 'x(i-1)=0,x(i)=1 => p(i)=1' qbound 'y(i-1)=0,y(i)=1 => q(i)=1' plim 'sum(p) <= 1' qlim 'sum(q) <= 1' ; * objective is linearized automatically by Cplex obj.. z =e= sum((i,j),spA(i,j)*x(i)*y(j)); pbound(i).. p(i) =g= x(i)-x(i-1); qbound(j).. q(j) =g= y(j)-y(j-1); plim.. sum(i,p(i)) =l= 1; qlim.. sum(j,q(j)) =l= 1; model m /all/; option miqcp=cplex,optcr=0,threads=8; solve m maximizing z using miqcp; display z.l; parameter results(*,*) 'model B results'; results('min','x') = smin(i$(x.l(i)>0.5),sum((k,j),spmap(k,i,j,'x'))); results('max','x') = smax(i$(x.l(i)>0.5),sum((k,j),spmap(k,i,j,'x'))); results('min','y') = smin(j$(y.l(j)>0.5),sum((k,i),spmap(k,i,j,'y'))); results('max','y') = smax(j\$(y.l(j)>0.5),sum((k,i),spmap(k,i,j,'y'))); display results;

1. Can you please share the larger data?

1. https://amsterdamoptimization.com/data/p1000.csv

2. Thanks for sharing the data. My colleague Imre Polik suggested treating the points as elements of a sparse matrix and using your formulation from http://yetanothermathprogrammingconsultant.blogspot.com/2021/01/submatrix-with-largest-sum.html. This yielded a 20x speedup in solution time for us. For n points, the model has 4n variables and 2n+2 constraints. After linearizing the n products (not n^2), the model has 5n variables and 5n+2 constraints.

3. That is a great observation. I get a similar performance with Cplex.