This is a note about BRATIO, a very exotic GAMS option related to providing an advanced basis to a Simplex-based LP (or NLP) solver.

I am a full-time consultant and provide services related to the design, implementation and deployment of mathematical programming, optimization and data-science applications. I also teach courses and workshops. Usually I cannot blog about projects I am doing, but there are many technical notes I'd like to share. Not in the least so I have an easy way to search and find them again myself. You can reach me at erwin@amsterdamoptimization.com.

## Thursday, April 29, 2021

## Tuesday, April 20, 2021

### Inverting a matrix by LP

## Monday, April 19, 2021

### Parallel machine scheduling II, three more formulations

In [1] two formulations were discussed for a **scheduling problem with multiple machines**. Here we add a few more. Some of them are a bit strange. None of them really works better than the ones in [1].

So this is a collection of formulations that may sound reasonable at first, but are really not suited for this problem. If you want to read about good models, skip this post.

## Monday, April 5, 2021

### This is really a modeling issue

In [1], the following scenario is described:

- A MIP model is solved, delivering a
**feasible**solution. - All integer variables are fixed to their levels.
- The resulting LP model is solved in order to produce duals. This model turns out to be
**infeasible**.

## Monday, March 29, 2021

### Mathematical Notation

I often tell beginners in optimization modeling that before starting to code, they should get a piece of paper and write down the mathematical model. There are several reasons for this:

## Sunday, March 28, 2021

### Parallel machine scheduling I, two formulations

## Tuesday, March 23, 2021

### SOS1 sets: when to use them

Short answer: almost never.

Beginners in optimization are always very excited when they read about SOS sets in the documentation. This looks like a very efficient tool to use for constructs like pick (at most) 1 out of \(n\) (SOS1) or for interpolation (SOS2). However, experienced modelers use SOS sets very sparingly. Why? Well, solvers like binary variables much better than SOS sets. Binary variables yield great bounding and cuts (there has been enormous progress in this area), while SOS constructs are really just about branching.

## Wednesday, March 17, 2021

### MST + few Steiner points: convex quadratic model

In [1], I discussed a model for the **Euclidean Streiner Tree problem**. A non-convex integer programming model from the literature was reformulated into a MISOCP (not a version of a Japanese dish, but a **Mixed-Integer Second-Order Cone Program**). Together with a symmetry-breaking constraint and some high-performance solvers, this can lead to being able to solve larger models.

One of the disadvantages of the model in [1] is that it only works with the full number of Steiner points \(n-2\) where \(n\) is the number of given (original) points. Here I want to discuss how we can use Minimum Spanning Tree (MST) models [2], as a building block for a model where we can add just a few Steiner points. One reason this is interesting: the first few Steiner points have the most impact on the objective.

## Tuesday, March 16, 2021

### Excel Never Dies

In a post by (famous economist) Paul Krugman:

I’m sometimes embarrassed at how much I use Excel. But it turns out to be a work of genius.

The online article he refers to is:

Packy McCormick, Ben Rollert,

Excel Never Dies, https://www.notboring.co/p/excel-never-dies

It is a nice read.

In my work, I see Excel intensively used by (almost) all clients I work with or talk to. One sometimes wonders: what would happen if Excel would just suddenly stop working everywhere?

## Monday, March 15, 2021

### Euclidean Steiner tree problems as MISOCP

In [1] I looked at MIP models for the Minimum Spanning Tree problem. A related, but much more difficult problem is the **Euclidean Steiner Tree Problem**. Here we are allowed to add points to make the tree shorter. Here are two obvious examples [2]:

## Friday, March 12, 2021

## Monday, March 8, 2021

### Minimum Spanning Trees in Math Programming Models

Algorithms for the Minimum Spanning Tree (MST) problem are readily available. However, sometimes we want to solve this problem inside a Mathematical Programming model. Usually, this is for two reasons:

- We have some side constraints
- Or as part of a larger model

**Data**

---- 298 SETcitiesManchester, N.H. , Montpelier, Vt. , Detroit, Mich. , Cleveland, Ohio , Charleston, W.Va. Louisville, Ky. , Indianapolis, Ind. , Chicago, Ill. , Milwaukee, Wis. , Minneapolis, Minn. Pierre, S.D. , Bismarck, N.D. , Helena, Mont. , Seattle, Wash. , Portland, Ore. Boise, Idaho , Salt Lake City, Utah, Carson City, Nevada , Los Angeles, Calif. , Phoenix, Ariz. Santa Fe, N.M. , Denver, Colo. , Cheyenne, Wyo. , Omaha, Neb. , Des Moines, Iowa Kansas City, Mo. , Topeka, Kans. , Oklahoma City, Okla., Dallas, Tex. , Little Rock, Ark. Memphis, Tenn. , Jackson, Miss. , New Orleans, La. , Birmingham, Ala. , Atlanta, Ga. Jacksonville, Fla. , Columbia, S.C. , Raleigh, N.C. , Richmond, Va. , Washington, D.C. Boston, Mass. , Portland, Me.

**dantzig42**from TSPLIB [8].

## Monday, March 1, 2021

### A difficult multi-line regression data set

In [1], I discussed some models for a small 2 line regression problem. Here I want to dive into a more challenging problem.

#### Data generation process

## Thursday, February 25, 2021

### A strange regression problem: multiple-line fitting

Consider the following data:

---- 14 PARAMETERxdatai1 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 PARAMETERydatai1 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 PARAMETERitemSizesizes of bin and itemsw 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 SETitemMapmapping between items and categoriesitem1 .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 PARAMETERitemSizesizes of bin and itemsw 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 SETitemMapmapping between items and categoriesitem1 .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:

#### Non-differentiable NLP model

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

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

**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:

- Dijkstra just does not work correctly if there are negative weights.
- 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.
- 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:

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.

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

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.

#### Graph

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

**source node**and a

**sink node**.

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

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

## 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}\] |

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

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

---- 10 PARAMETERppointsx 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

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}\] |

## Tuesday, January 19, 2021

### Submatrix with Largest Sum

#### Problem

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

#### Example data

---- 8 PARAMETERamatrix (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

#### Non-contiguous submatrix

**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}\] |