Thursday, February 25, 2021

A strange regression problem: multiple line fitting

Consider the following data:


----     14 PARAMETER x  data

i1  17.175,    i2  84.327,    i3  55.038,    i4  30.114,    i5  29.221,    i6  22.405,    i7  34.983,    i8  85.627
i9   6.711,    i10 50.021,    i11 99.812,    i12 57.873,    i13 99.113,    i14 76.225,    i15 13.069,    i16 63.972
i17 15.952,    i18 25.008,    i19 66.893,    i20 43.536,    i21 35.970,    i22 35.144,    i23 13.149,    i24 15.010
i25 58.911,    i26 83.089,    i27 23.082,    i28 66.573,    i29 77.586,    i30 30.366,    i31 11.049,    i32 50.238
i33 16.017,    i34 87.246,    i35 26.511,    i36 28.581,    i37 59.396,    i38 72.272,    i39 62.825,    i40 46.380
i41 41.331,    i42 11.770,    i43 31.421,    i44  4.655,    i45 33.855,    i46 18.210,    i47 64.573,    i48 56.075
i49 76.996,    i50 29.781


----     14 PARAMETER y  data

i1    32.320,    i2   238.299,    i3    87.081,    i4    17.825,    i5   -10.154,    i6   -35.813,    i7    14.506
i8   204.903,    i9  -103.014,    i10   79.521,    i11 -153.679,    i12  -47.971,    i13  277.483,    i14  160.084
i15   28.687,    i16  133.650,    i17  -63.820,    i18    4.028,    i19  145.434,    i20   48.326,    i21   49.295
i22   35.478,    i23  -58.602,    i24   34.463,    i25  -50.666,    i26 -115.522,    i27  -19.540,    i28  134.829
i29  198.074,    i30   -7.351,    i31   40.786,    i32   82.107,    i33   27.661,    i34  -82.669,    i35    5.081
i36   -1.278,    i37  118.345,    i38  172.905,    i39  -57.943,    i40   56.981,    i41  -45.248,    i42   34.517
i43   11.780,    i44  -98.506,    i45   -1.064,    i46   33.323,    i47  139.050,    i48  -35.595,    i49 -102.580
i50    3.912


When we plot this, we see:




Saturday, February 20, 2021

2d bin packing with Google OR-Tools CP-SAT, how to fix?

In [1] I looked at some MIP formulations for a 2d bin packing problem. OR-Tools has a nice model.AddNoOverlap2D constraint function. So, I was excited to try this out. (In these times I get quickly enthusiastic about these things). 

The idea is as follows. ortools has this construct Constraint.OnlyEnforce If(boolean). This is essentially an indicator constraint: \((\color{darkred}x\ge 5)\text{.OnlyEnforceIf($\color{darkred}\delta$)}\) means \[\color{darkred}\delta=1 \Rightarrow \color{darkred}x \ge 5\] This inspires me to develop a model along the lines of: \[\begin{align} \min& \sum_k \color{darkred}u_k\\  & \color{darkred}\beta_{i,j} = 0 \Rightarrow \color{darkred}b_i\ne\color{darkred}b_j \\ & \color{darkred}\beta_{i,j} = 1 \Rightarrow \mathbf{NoOverlap2D}(i,j) \\ & \color{darkred}u_k = 0 \Rightarrow \color{darkred}b_i \lt k \end{align}\] Here \(\color{darkred}b_i\) is the bin number where item \(i\) is placed and \(\color{darkred}\beta_{i,j}\) is true if  \(\color{darkred}b_i=\color{darkred}b_j\).


OK, so I started coding:

Thursday, February 18, 2021

An easy variant of 2d bin packing

 In [1] the following problem is posted:

Consider the 2d bin packing problem [2]. Now we have the following additional restrictions: items of a different type or category cannot be below or above each other. In other words, we arrange items of the same type in columns. 

We use again our small data set: 


----     37 PARAMETER itemSize  sizes of bin and items

               w           h

cat1       7.000      12.000
cat2       9.000       3.000
cat3       5.000      14.000
cat4      13.000       9.000
cat5       6.000       8.000
cat6      20.000       5.000
bin       40.000      60.000


----     37 SET itemMap  mapping between items and categories

item1 .cat1,    item2 .cat1,    item3 .cat1,    item4 .cat1,    item5 .cat1,    item6 .cat1,    item7 .cat1
item8 .cat1,    item9 .cat1,    item10.cat1,    item11.cat2,    item12.cat2,    item13.cat2,    item14.cat2
item15.cat2,    item16.cat2,    item17.cat2,    item18.cat2,    item19.cat2,    item20.cat2,    item21.cat3
item22.cat3,    item23.cat3,    item24.cat3,    item25.cat3,    item26.cat3,    item27.cat3,    item28.cat3
item29.cat3,    item30.cat3,    item31.cat4,    item32.cat4,    item33.cat4,    item34.cat4,    item35.cat4
item36.cat4,    item37.cat4,    item38.cat4,    item39.cat4,    item40.cat4,    item41.cat5,    item42.cat5
item43.cat5,    item44.cat5,    item45.cat5,    item46.cat6,    item47.cat6,    item48.cat6,    item49.cat6
item50.cat6


I'll try two different models:

  • Assign items to bins, arrange them in columns
  • Assign columns to bins

Tuesday, February 16, 2021

2d Bin Packing

 In [1], a question was posed about a particular type of 2d bin-packing. This is a good reason to play with some models.

First I discuss three formulations with continuous no-overlap constraints and one for discrete data. In the next post, I'll discuss the variant described by the post [1]. This actually leads to a very different, but easy problem.


Problem statement and example data


We consider \(n=50\) rectangular items of different sizes that have to be placed into bins. There are six different categories, each with a different size (height and width). All the bins have the same size. The goal is to use as few bins as possible.


----     37 PARAMETER itemSize  sizes of bin and items

               w           h

cat1       7.000      12.000
cat2       9.000       3.000
cat3       5.000      14.000
cat4      13.000       9.000
cat5       6.000       8.000
cat6      20.000       5.000
bin       40.000      60.000


----     37 SET itemMap  mapping between items and categories

item1 .cat1,    item2 .cat1,    item3 .cat1,    item4 .cat1,    item5 .cat1,    item6 .cat1,    item7 .cat1
item8 .cat1,    item9 .cat1,    item10.cat1,    item11.cat2,    item12.cat2,    item13.cat2,    item14.cat2
item15.cat2,    item16.cat2,    item17.cat2,    item18.cat2,    item19.cat2,    item20.cat2,    item21.cat3
item22.cat3,    item23.cat3,    item24.cat3,    item25.cat3,    item26.cat3,    item27.cat3,    item28.cat3
item29.cat3,    item30.cat3,    item31.cat4,    item32.cat4,    item33.cat4,    item34.cat4,    item35.cat4
item36.cat4,    item37.cat4,    item38.cat4,    item39.cat4,    item40.cat4,    item41.cat5,    item42.cat5
item43.cat5,    item44.cat5,    item45.cat5,    item46.cat6,    item47.cat6,    item48.cat6,    item49.cat6
item50.cat6


Furthermore, we assume items cannot be rotated.

Friday, February 5, 2021

Non-differentiable NLP vs LP

 In [1] the following simple problem is posed:


Given data 

p = [50, 50.5, 52, 53, 55, 55.5, 56, 57,60, 61]

try to track this as close as possible by minimizing the largest deviation while obeying the restrictions:


\[\begin{align}&\color{darkred}x_1 = \color{darkblue}p_1\\ &\color{darkred}x_i - \color{darkred}x_{i-1}\le 1.5\end{align}\]


Non-differentiable NLP model


An obvious way to formulate this is: \[\begin{align} \min&\max_i |\color{darkred}x_i-\color{darkblue}p_i| \\ &\color{darkred}x_1 = \color{darkblue}p_1\\ &\color{darkred}x_i - \color{darkred}x_{i-1}\le 1.5\end{align}\]

This is a non-differentiable formulation, and we are asking for trouble when trying to solve this using a standard nonlinear NLP solver. For instance, when I solve this with CONOPT, we may see:


               S O L V E      S U M M A R Y

     MODEL   m1                  OBJECTIVE  z
     TYPE    DNLP                DIRECTION  MINIMIZE
     SOLVER  CONOPT              FROM LINE  25

**** SOLVER STATUS     4 Terminated By Solver      
**** MODEL STATUS      7 Feasible Solution         
**** OBJECTIVE VALUE                1.2500

 RESOURCE USAGE, LIMIT          0.062 10000000000.000
 ITERATION COUNT, LIMIT        31    2147483647
 EVALUATION ERRORS              0             0
CONOPT 3         33.2.0 r4f23b21 Released Dec 01, 2020 WEI x86 64bit/MS Window
 
 
    C O N O P T 3   version 3.17L
    Copyright (C)   ARKI Consulting and Development A/S
                    Bagsvaerdvej 246 A
                    DK-2880 Bagsvaerd, Denmark
 
 
    The model has 11 variables and 10 constraints
    with 29 Jacobian elements, 10 of which are nonlinear.
    The Hessian of the Lagrangian has 10 elements on the diagonal,
    45 elements below the diagonal, and 10 nonlinear variables.
 
                   Pre-triangular equations:   0
                   Post-triangular equations:  10
 
 
 ** Feasible solution. Convergence too slow. The change in objective
    has been less than 3.7500E-12 for 20 consecutive iterations


Conopt realizes it has some troubles and declares the solution as just feasible instead of optimal. Indeed, we can do better than an objective of 1.25. Other NLP solvers may not even know they are in trouble and just report a solution and declare it (locally) optimal. 

There is a big problem with this model. Most NLP solvers assume the objective and non-linear constraints are smooth (i.e. differentiable). In this case, we don't have this, and this lack of smoothness is often a recipe for disaster. This is a nice example of this. 

I often complain about non-differentiability in models I see on forums. Often the reaction is: "Who cares?". Indeed, often NLP solvers will deliver some solution without warning. This post is about a small problem, where we really can do much better.   

Note that the initial point also plays a role here. With a different starting point, you may find different solutions.

A global NLP solver may help for this model. Some of the ones I tried, don't have support for the max() function, so we need to reformulate the model. When doing that, we can actually do a bit better and form a linear programming problem.

Wednesday, February 3, 2021

Dijkstra confusion

 Dijkstra is a famous algorithm for the shortest path problem. The original paper is from 1959 [1]:


(The format of this paper is interesting: almost no mathematical notation).  Dijkstra's algorithm is not suited for all shortest path problems. Negative weights or lengths can cause problems. There are basically three different reasons I see mentioned:

  1. Dijkstra just does not work correctly if there are negative weights.
  2. Dijkstra does not work if there are negative cycles in the graph. This would imply that a problem with negative weights would work as long as there are no negative cycles. 
  3. Keep things vague and mumble a bit about negative weights and negative cycles.

I notice that option 3 seems very popular. As a result, looking at questions about this on stackexchange, I see there is much confusion about this. Also, I see vague, misleading messages in documentation and error/warning messages. For instance, let's have a look at [2]:

Also, this routine does not work for graphs with negative distances. Negative distances can lead to infinite cycles that must be handled by specialized algorithms such as Bellman-Ford’s algorithm or Johnson’s algorithm.

Well, this is indeed sufficiently vague. We deal with option 3: we don't know if this is referring to reason 1 or 2. The same function that this documentation is about can produce the following warning:

UserWarning: Graph has negative weights: dijkstra will give inaccurate results if the graph contains negative cycles. Consider johnson or bellman_ford. distmat,pred = scipy.sparse.csgraph.dijkstra(spmat,indices=source,return_predecessors=True)

This seems to say we only need to worry about graphs with negative cycles. So that would be option 2. 


A small counter example


In [3] a small acyclic graph (a DAG: Directed Acyclic Graph) is presented:


When we try to find the shortest path from n0 to n2, using [2], we can see:

dijkstra
shortest path length: 1.0
johnson
shortest path length: -8.0
bellman_ford
shortest path length: -8.0
<ipython-input-25-7a1de3894d98>:13: UserWarning: Graph has negative weights: dijkstra will give inaccurate results if the graph contains negative cycles. Consider johnson or bellman_ford.
  distmat,pred = scipy.sparse.csgraph.dijkstra(spmat,indices=source,return_predecessors=True)


This example disproves that we just have to watch out for negative cycles.  That means the documentation and warning messages should state clearly that the function scipy.sparse.csgraph.dijkstra does not work correctly with negative weights, irrespective of whether there are negative cycles.

Thursday, January 28, 2021

Assignment problem with a wrinkle formulated as a network problem

In a previous post [1] I discussed a simple problem. But it turned out not so easy to solve for some of the larger data sets. Basically, it was an assignment problem with an extra condition. The problem was a follows:

Consider two arrays \(\color{darkblue}a_i\) (length \(\color{darkblue}m\)) and \(\color{darkblue}b_j\) (length \(\color{darkblue}n\)) with \(\color{darkblue}m \lt \color{darkblue}n\). Assign all values \(\color{darkblue}a_i\) to a \(\color{darkblue}b_j\) such that:

  • Each \(\color{darkblue}b_j\) can have 0 or 1 \(\color{darkblue}a_i\) assigned to it.
  • The assignments need to maintain the original order of \(\color{darkblue}a_i\). I.e. if \(\color{darkblue}a_i \rightarrow \color{darkblue}b_j\) then \(\color{darkblue}a_{i+1}\) must be assigned to a slot in \(\color{darkblue}b\) that is beyond slot \(j\). In the picture below that means that arrows cannot cross.
  • Do this while minimizing the sum of the products.


In [1], I attacked this as a mixed-integer programming problem. In this post, I want to see if we can solve this as a network problem. This was largely inspired by the comments in [1].


Shortest path algorithms and models


The shortest path problem can be solved in different ways. Dijkstra's algorithm is a popular choice due to its simplicity and speed. In [2] we can read a quote by Edsger Dijkstra:

What is the shortest way to travel from Rotterdam to Groningen, in general: from given city to given city. It is the algorithm for the shortest path, which I designed in about twenty minutes. One morning I was shopping in Amsterdam with my young fiancée, and tired, we sat down on the café terrace to drink a cup of coffee and I was just thinking about whether I could do this, and I then designed the algorithm for the shortest path. As I said, it was a twenty-minute invention. In fact, it was published in '59, three years later. The publication is still readable, it is, in fact, quite nice. One of the reasons that it is so nice was that I designed it without pencil and paper. I learned later that one of the advantages of designing without pencil and paper is that you are almost forced to avoid all avoidable complexities. Eventually, that algorithm became to my great amazement, one of the cornerstones of my fame.
For an early implementation in a computer program, Dijkstra used a set of 64 nodes corresponding to Dutch cities. 64 was chosen so the node number could be encoded in 6 bits [2]. 

The shortest path can also be formulated as an LP problem. Standard LP solvers can handle quite large problems (each node becomes a constraint, and each arc a variable). In addition, a solver like Cplex contains a specialized network solver for such LPs.  

Graph 


We start with the nodes. We denote the nodes by \(\color{darkblue}n_{i,j}\) representing: \(\color{darkblue}a_i\) is assigned to \(\color{darkblue}b_j\). Not all assignments are possible. For instance, we cannot assign \(\color{darkblue}a_2\) to \(\color{darkblue}b_1\). That means: node  \(\color{darkblue}n_{2,1}\)  does not exist.

We also need a source node and a sink node.

The arcs indicate: after assigning \(\color{darkblue}a_i \rightarrow \color{darkblue}b_j\) we need to assign the next item \(\color{darkblue}a_{i+1}\rightarrow \color{darkblue}b_{j+k}\) for some \(k\ge 1\). In addition we need to connect the source node to all nodes with \(i=1\) and the sink node to all nodes with \(i=3\) (our last \(\color{darkblue}a_i\)).  So our network looks like:

Network representation


Note how any arc not connected to the sink or source node, goes to the right and downwards.

We have costs for visiting each node: \(\color{darkblue}a_i\cdot\color{darkblue}b_j\).  As we want to formulate this as a shortest path problem, we need to allocate these costs to arcs. I used the incoming arcs for this. So any arc \(e_{i,j,i',j'}\) gets a cost \(\color{darkblue}a_{i'}\cdot\color{darkblue}b_{j'}\).  The arcs to the sink node have a zero cost. 

Note: this graph is acyclic, so negative arc lengths are ok. We will never see negative cycles.

As you can see: because the nodes have two indices, the arcs have four! That will be fun.

Wednesday, January 27, 2021

An assignment problem with a wrinkle

In [1], the following problem is posed:

Consider two arrays \(\color{darkblue}a_i\) (length \(\color{darkblue}m\)) and \(\color{darkblue}b_j\) (length \(\color{darkblue}n\)) with \(\color{darkblue}m \lt \color{darkblue}n\). Assign all values \(\color{darkblue}a_i\) to a \(\color{darkblue}b_j\) such that:

  • Each \(\color{darkblue}b_j\) can have 0 or 1 \(\color{darkblue}a_i\) assigned to it.
  • The assignments need to maintain the original order of \(\color{darkblue}a_i\). I.e. if \(\color{darkblue}a_i \rightarrow \color{darkblue}b_j\) then \(\color{darkblue}a_{i+1}\) must be assigned to a slot in \(\color{darkblue}b\) that is beyond slot \(j\). In the picture below that means that arrows cannot cross.
  • Do this while minimizing the sum of the products.


Mixed-integer programming model

This problem can be viewed as an assignment problem with a side constraint:


MIP Model
\[\begin{align}\min& \sum_{i,j}\color{darkred}x_{i,j}\cdot\color{darkblue}a_i\cdot\color{darkblue}b_j \\ &\sum_j \color{darkred}x_{i,j}=1 &&\forall i\\ & \sum_i \color{darkred}x_{i,j}\le 1 &&\forall j\\ & \color{darkred}v_i = \sum_j j \cdot \color{darkred} x_{i,j}\\ & \color{darkred}v_i \ge \color{darkred}v_{i-1}+1\\ & \color{darkred}x_{i,j} \in \{0,1\} \\ & \color{darkred}v_i \ge 1\end{align}\]


Here \(\color{darkred}x_{i,j}\) indicates the assignment \(\color{darkblue}a_i \rightarrow \color{darkblue}b_j\). The variable \(\color{darkred}v_i\) represents the position in \(\color{darkblue}b\) to which \(\color{darkblue}a_i\) is assigned.

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

Tuesday, January 19, 2021

Submatrix with Largest Sum

Problem


Given a data matrix \(\color{darkblue}a_{i,j}\), find a submatrix such that the sum of the elements is maximized. The term "submatrix" can indicate a contiguous or a non-contiguous subset of rows and columns. 

This is a generalization of the 1d maximum subarray problem [1].

Looks easy, but let's make sure it is.

Example data


Of course, the problem is trivial if all elements are non-negative: just select the whole matrix. It is only interesting if there are negative elements. So here is my random matrix:


----      8 PARAMETER a  matrix (data)

             c1          c2          c3          c4          c5          c6          c7          c8          c9         c10

r1       -6.565       6.865       1.008      -3.977      -4.156      -5.519      -3.003       7.125      -8.658       0.004
r2        9.962       1.575       9.823       5.245      -7.386       2.794      -6.810      -4.998       3.379      -1.293
r3       -2.806      -2.971      -7.370      -6.998       1.782       6.618      -5.384       3.315       5.517      -3.927
r4       -7.790       0.048      -6.797       7.449      -4.698      -4.284       1.879       4.454       2.565      -0.724
r5       -1.734      -7.646      -3.716      -9.069      -3.229      -6.358       2.915       1.215       5.399      -4.044
r6        3.222       5.116       2.549      -4.323      -8.272      -7.950       2.825       0.906      -9.370       5.847
r7       -8.545      -6.487       0.513       5.004      -6.438      -9.317       1.703       2.425      -2.213      -2.826
r8       -5.139      -5.072      -7.390       8.669      -2.401       5.668      -3.999      -7.490       4.977      -8.615
r9       -5.960      -9.899      -4.608      -0.003      -6.974      -6.517      -3.387      -3.662      -3.558       9.280
r10       9.872      -2.602      -2.542       5.440      -2.066       8.262      -7.608       4.710      -8.892       1.526
r11      -8.972      -9.880      -1.975       0.398       2.578      -5.485      -2.078      -4.480      -6.953       8.726
r12      -1.547      -7.307      -2.279      -2.507      -4.630       8.967      -6.221      -4.050      -8.509      -1.973
r13      -7.966      -2.322      -3.518      -6.157      -7.753       1.931       0.229      -9.099       5.662       8.915
r14       1.929       2.147      -2.750       1.881       3.597       0.132      -6.815       3.138       0.478      -7.512
r15       9.734      -5.438       3.513       5.536       8.649      -5.975      -4.057      -6.055      -5.073       2.930
r16       4.699      -8.291      -6.993      -1.316      -6.261       3.854       5.259      -6.904      -2.212       3.909
r17       6.916       2.254       9.519      -9.462      -6.251      -8.258       0.808      -7.463       4.680      -7.735
r18      -0.233       5.912      -0.159       0.671      -9.788       0.877      -0.977       9.507      -6.323      -6.729
r19      -9.507      -6.444      -8.774      -9.667       6.713       2.033      -9.460      -6.078       9.014      -3.289
r20       1.885      -4.816       2.813      -6.895      -0.800      -2.133       6.109       0.820      -2.186       1.156
r21       8.655      -3.025      -9.834       8.977       1.438      -3.327       9.675       5.329      -7.798       9.896
r22       1.606      -6.672       2.867      -3.114       8.247       8.001      -9.675      -2.627       3.288       1.868
r23      -9.309       6.836       8.642       0.159      -4.008      -0.068      -9.101       5.474       0.659       4.935
r24       4.401       2.632      -7.702       9.423       4.135       9.725       7.096       2.429       4.026       4.018
r25       5.814       2.204      -8.914      -0.296      -8.949       3.972      -6.104      -5.479       6.273       9.835
r26       5.013       4.367      -9.988      -4.723       6.476       6.391       7.208      -5.746      -0.864      -9.233
r27      -3.540      -1.202      -3.693      -7.305       6.219      -1.664      -7.164      -0.689      -4.340       7.914
r28      -8.712      -1.708      -3.168      -0.634       2.853       2.872      -3.248      -7.984       8.117      -5.653
r29       8.377      -0.965      -8.201      -2.516      -1.700      -1.916      -7.767       5.023       6.068      -9.527
r30      -0.382      -4.428       8.032      -9.648       3.621       9.018       8.004       7.976       7.489      -2.180


The values are drawn from a uniform distribution, \(\color{darkblue}a_{i,j} \sim U(-10,10)\).

Non-contiguous submatrix


The non-convex MIQP formulation is easy. Just introduce binary variables that indicate if a row or column is selected. A cell is selected if both its corresponding row and column are selected.

MIQP model A
\[\begin{align} \max& \sum_{i,j} \color{darkblue}a_{i,j} \cdot \color{darkred}x_i \cdot \color{darkred}y_j\\ &\color{darkred}x_i, \color{darkred}y_j \in \{0,1\} \end{align}\]