Modeling permutations

Introduction

In [1] the following question was posted:

Given a (square) matrix \(A\), reorder the columns such that the sum of the diagonals is maximized. We are working in an R environment.

My answer involves permutation matrices and assignment constraints. This is a useful concept for quite a few models.

Permutation Matrix

When dealing with permutations it is often a good idea to look at this in terms of a Permutation Matrix [2]. A permutation matrix \(P\) is an identity matrix with rows or columns permuted.

Interestingly, it does not matter whether we swap columns or rows here. The results are the same:

> n <- 5
> # identity matrix
> I <- diag(n)  
> print(I)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1
> # swap columns 2 and 4
> P1 <- I
> P1[,c(2,4)] <- I[,c(4,2)]
> print(P1)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    0    0    1    0
[3,]    0    0    1    0    0
[4,]    0    1    0    0    0
[5,]    0    0    0    0    1
> # swap rows 2 and 4
> P2 <- I
> P2[c(2,4),] <- I[c(4,2),]
> print(P2)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    0    0    1    0
[3,]    0    0    1    0    0
[4,]    0    1    0    0    0
[5,]    0    0    0    0    1
> 

A permutation matrix has the following properties:

1. It is a square \(n \times n\) matrix with \(n\) entries equal to 1 and the other entries equal to 0,
2. exactly one element is 1 in each row,
3. and exactly one element is 1 in each column. 

In an optimization model, this can be modeled as a set of assignment constraints: \[\begin{align}&\sum_j p_{i,j} = 1 && \forall i \\ &\sum_i p_{i,j} =1 && \forall j \\ &p_{i,j} \in \{0,1\}\end{align}\]

Row and column permutations

From [2] we can see:

• Applying a row permutation to an \((m\times n)\)  matrix \(A\) can be viewed as a pre-multiplication with an \((m \times m)\) permutation matrix: \[\widetilde{A}=PA\]
• and applying a column permutation to an \((m\times n)\)  matrix \(A\) can be viewed as a post-multiplication with an \((n \times n)\) permutation matrix: \[\widetilde{A}=AP\]

As seeing is believing, below is a little experiment. We start with a \(2 \times 3\) matrix \(A\) and permute some rows and columns. Note: the original question was about a square matrix, while here we use a non-square matrix to emphasize different sizes of the permutation matrices \(P_1\) and \(P_2\).

> m <- 2
> n <- 3
> # random matrix A
> set.seed(123)
> A <- matrix(runif(m*n,min=-10,max=+10),nrow=m,ncol=n)
> print(A)
          [,1]      [,2]      [,3]
[1,] -4.248450 -1.820462  8.809346
[2,]  5.766103  7.660348 -9.088870
> # swap rows 1 and 2
> P1 <- diag(m)
> P1[,c(1,2)] <- diag(m)[,c(2,1)]
> P1 %*% A
          [,1]      [,2]      [,3]
[1,]  5.766103  7.660348 -9.088870
[2,] -4.248450 -1.820462  8.809346
> # swap columns 1 and 3
> P2 <- diag(n)
> P2[,c(1,3)] <- diag(n)[,c(3,1)]
> A %*% P2
          [,1]      [,2]      [,3]
[1,]  8.809346 -1.820462 -4.248450
[2,] -9.088870  7.660348  5.766103
> 

Note: If we want to apply both row and column permutations, we would write: \[\widetilde{A}=P A Q\] with \(P\) and \(Q\) permutation matrices.

Optimization model

The original problem can now be formulated as a Mixed Integer Programming model. The mathematical model can look like:

Mixed Integer Programming Model
\[\begin{align}\max&\sum_i \color{darkred}y_{i,i} \\ & \color{darkred}y_{i,j} = \sum_k \color{darkblue}a_{i,k} \cdot \color{darkred}p_{k,j} \\ & \sum_j \color{darkred}p_{i,j} = 1 &&\forall i \\ & \sum_i \color{darkred}p_{i,j} = 1 &&\forall j \\ & \color{darkred}p_{i,j} \in \{0,1\} \\ & \color{darkred}y_{i,j} \>\mathbf{free}\end{align}\]

Implementation 1: OMPR

The model above can be implemented in a straightforward manner in OMPR.

> library(ompr)
> library(ompr.roi)
> library(dplyr)
> library(ROI)
> library(ROI.plugin.symphony)
> library(tidyr)
> n <- 10
> set.seed(123)
> a <- matrix(runif(n^2,min=-1,max=1),nrow=n,ncol=n)
> print(a)
             [,1]        [,2]        [,3]        [,4]       [,5]       [,6]       [,7]        [,8]       [,9]       [,10]
 [1,] -0.42484496  0.91366669  0.77907863  0.92604847 -0.7144000 -0.9083377  0.3302304  0.50895032 -0.5127611 -0.73860862
 [2,]  0.57661027 -0.09333169  0.38560681  0.80459809 -0.1709073 -0.1155999 -0.8103187  0.25844226  0.3361112  0.30620385
 [3,] -0.18204616  0.35514127  0.28101363  0.38141056 -0.1725513  0.5978497 -0.2320607  0.42036480 -0.1647064 -0.31296706
 [4,]  0.76603481  0.14526680  0.98853955  0.59093484 -0.2623091 -0.7562015 -0.4512327 -0.99875045  0.5763917  0.31351626
 [5,]  0.88093457 -0.79415063  0.31141160 -0.95077263 -0.6951105  0.1218960  0.6292801 -0.04936685 -0.7942707 -0.35925352
 [6,] -0.90888700  0.79964994  0.41706094 -0.04440806 -0.7223879 -0.5869372 -0.1029673 -0.55976223 -0.1302145 -0.62461776
 [7,]  0.05621098 -0.50782453  0.08813205  0.51691908 -0.5339318 -0.7449367  0.6201287 -0.24036692  0.9699140  0.56458860
 [8,]  0.78483809 -0.91588093  0.18828404 -0.56718413 -0.0680751  0.5066157  0.6247790  0.22554201  0.7861022 -0.81281003
 [9,]  0.10287003 -0.34415856 -0.42168053 -0.36363798 -0.4680547  0.7900907  0.5886846 -0.29640418  0.7729381 -0.06644192
[10,] -0.08677053  0.90900730 -0.70577271 -0.53674843  0.7156554 -0.2510744 -0.1203366 -0.77772915 -0.6498947  0.02301092
> m <- MIPModel() %>%
+   add_variable(p[i,j], i=1:n, j=1:n, type="binary") %>%
+   add_variable(y[i,j], i=1:n, j=1:n) %>%
+   add_constraint(y[i,j] == sum_expr(a[i,k]*p[k,j],k=1:n),i=1:n,j=1:n) %>%
+   add_constraint(sum_expr(p[i,j],j=1:n) == 1,i=1:n) %>%
+   add_constraint(sum_expr(p[i,j],i=1:n) == 1,j=1:n) %>%
+   set_objective(sum_expr(y[i,i], i=1:n),"max") %>%
+   solve_model(with_ROI(solver = "symphony",verbosity=1)) 
Starting Preprocessing...
Preprocessing finished...
   with no modifications...
Problem has 
  120 constraints 
  200 variables 
  1300 nonzero coefficients

Total Presolve Time: 0.000000...

Solving...

granularity set at 0.000000
solving root lp relaxation
The LP value is: -7.422 [0,20]


****** Found Better Feasible Solution !
****** Cost: -7.422180


****************************************************
* Optimal Solution Found                           *
* Now displaying stats and best solution found...  *
****************************************************

====================== Misc Timing =========================
  Problem IO        0.000
======================= CP Timing ===========================
  Cut Pool                  0.000
====================== LP/CG Timing =========================
  LP Solution Time          0.000
  LP Setup Time             0.000
  Variable Fixing           0.000
  Pricing                   0.000
  Strong Branching          0.000
  Separation                0.000
  Primal Heuristics         0.000
  Communication             0.000
  Total User Time              0.000
  Total Wallclock Time         0.000

====================== Statistics =========================
Number of created nodes :       1
Number of analyzed nodes:       1
Depth of tree:                  0
Size of the tree:               1
Number of solutions found:      1
Number of solutions in pool:    1
Number of Chains:               1
Number of Diving Halts:         0
Number of cuts in cut pool:     0
Lower Bound in Root:            -7.422

======================= LP Solver =========================
Number of times LP solver called:                 1
Number of calls from feasibility pump:            0
Number of calls from strong branching:            0
Number of solutions found by LP solve:            1
Number of bounds changed by strong branching:     0
Number of nodes pruned by strong branching:       0
Number of bounds changed by branching presolver:  0
Number of nodes pruned by branching presolver:    0

==================== Primal Heuristics ====================
                             Time      #Called   #Solutions
Rounding I                   0.00                           
Rounding II                  0.00                           
Diving                       0.00 
Feasibility Pump             0.00 
Local Search                 0.00            1            0 
Restricted Search            0.00 
Rins Search                  0.00 
Local Branching              0.00 

=========================== Cuts ==========================
Accepted:                         0
Added to LPs:                     0
Deleted from LPs:                 0
Removed because of bad coeffs:    0
Removed because of duplicacy:     0
Insufficiently violated:          0
In root:                          0

Time in cut generation:              0.00
Time in checking quality and adding: 0.00

                   Time     #Called     In Root       Total
Gomory             0.00 
Knapsack           0.00 
Clique             0.00 
Probing            0.00 
Flowcover          0.00 
Twomir             0.00 
Oddhole            0.00 
Mir                0.00 
Rounding           0.00 
LandP-I            0.00 
LandP-II           0.00 
Redsplit           0.00 

===========================================================
Solution Found: Node 0, Level 0
Solution Cost: -7.4221803085
> cat("Status:",solver_status(m),"\n")
Status: optimal 
> cat("Objective:",objective_value(m),"\n")
Objective: 7.42218 
> df <- get_solution(m,y[i, j])
> spread(df,key=j,value=value)[,-c(1,2)]
             1           2           3           4           5           6          7          8          9         10
1   0.92604847 -0.73860862  0.50895032  0.77907863 -0.42484496  0.91366669 -0.5127611  0.3302304 -0.9083377 -0.7144000
2   0.80459809  0.30620385  0.25844226  0.38560681  0.57661027 -0.09333169  0.3361112 -0.8103187 -0.1155999 -0.1709073
3   0.38141056 -0.31296706  0.42036480  0.28101363 -0.18204616  0.35514127 -0.1647064 -0.2320607  0.5978497 -0.1725513
4   0.59093484  0.31351626 -0.99875045  0.98853955  0.76603481  0.14526680  0.5763917 -0.4512327 -0.7562015 -0.2623091
5  -0.95077263 -0.35925352 -0.04936685  0.31141160  0.88093457 -0.79415063 -0.7942707  0.6292801  0.1218960 -0.6951105
6  -0.04440806 -0.62461776 -0.55976223  0.41706094 -0.90888700  0.79964994 -0.1302145 -0.1029673 -0.5869372 -0.7223879
7   0.51691908  0.56458860 -0.24036692  0.08813205  0.05621098 -0.50782453  0.9699140  0.6201287 -0.7449367 -0.5339318
8  -0.56718413 -0.81281003  0.22554201  0.18828404  0.78483809 -0.91588093  0.7861022  0.6247790  0.5066157 -0.0680751
9  -0.36363798 -0.06644192 -0.29640418 -0.42168053  0.10287003 -0.34415856  0.7729381  0.5886846  0.7900907 -0.4680547
10 -0.53674843  0.02301092 -0.77772915 -0.70577271 -0.08677053  0.90900730 -0.6498947 -0.1203366 -0.2510744  0.7156554
> 

Implementation 2: CVXR

The model can also be expressed conveniently in CVXR:

> library(CVXR)
> set.seed(123)
> n <- 10
> A <- matrix(runif(n^2,min=-1,max=1),nrow=n,ncol=n)
> print(A)
             [,1]        [,2]        [,3]        [,4]       [,5]       [,6]       [,7]        [,8]       [,9]       [,10]
 [1,] -0.42484496  0.91366669  0.77907863  0.92604847 -0.7144000 -0.9083377  0.3302304  0.50895032 -0.5127611 -0.73860862
 [2,]  0.57661027 -0.09333169  0.38560681  0.80459809 -0.1709073 -0.1155999 -0.8103187  0.25844226  0.3361112  0.30620385
 [3,] -0.18204616  0.35514127  0.28101363  0.38141056 -0.1725513  0.5978497 -0.2320607  0.42036480 -0.1647064 -0.31296706
 [4,]  0.76603481  0.14526680  0.98853955  0.59093484 -0.2623091 -0.7562015 -0.4512327 -0.99875045  0.5763917  0.31351626
 [5,]  0.88093457 -0.79415063  0.31141160 -0.95077263 -0.6951105  0.1218960  0.6292801 -0.04936685 -0.7942707 -0.35925352
 [6,] -0.90888700  0.79964994  0.41706094 -0.04440806 -0.7223879 -0.5869372 -0.1029673 -0.55976223 -0.1302145 -0.62461776
 [7,]  0.05621098 -0.50782453  0.08813205  0.51691908 -0.5339318 -0.7449367  0.6201287 -0.24036692  0.9699140  0.56458860
 [8,]  0.78483809 -0.91588093  0.18828404 -0.56718413 -0.0680751  0.5066157  0.6247790  0.22554201  0.7861022 -0.81281003
 [9,]  0.10287003 -0.34415856 -0.42168053 -0.36363798 -0.4680547  0.7900907  0.5886846 -0.29640418  0.7729381 -0.06644192
[10,] -0.08677053  0.90900730 -0.70577271 -0.53674843  0.7156554 -0.2510744 -0.1203366 -0.77772915 -0.6498947  0.02301092
> cat("sum diag of A:",sum(diag(A)))
sum diag of A: 0.7133438
> P <- Variable(n,n,boolean=T)
> Y <- Variable(n,n)
> problem <- Problem(Maximize(matrix_trace(Y)),
+                    list(Y==A %*% P,
+                         sum_entries(P,axis=1) == 1,
+                         sum_entries(P,axis=2) == 1))
> result <- solve(problem,verbose=T)
GLPK Simplex Optimizer, v4.47
120 rows, 200 columns, 1300 non-zeros
      0: obj =  0.000000000e+000  infeas = 2.000e+001 (120)
*   147: obj = -8.522275744e-002  infeas = 1.743e-014 (1)
*   213: obj = -7.422180308e+000  infeas = 4.854e-016 (1)
OPTIMAL SOLUTION FOUND
GLPK Integer Optimizer, v4.47
120 rows, 200 columns, 1300 non-zeros
100 integer variables, all of which are binary
Integer optimization begins...
+   213: mip =     not found yet >=              -inf        (1; 0)
+   213: >>>>> -7.422180308e+000 >= -7.422180308e+000   0.0% (1; 0)
+   213: mip = -7.422180308e+000 >=     tree is empty   0.0% (0; 1)
INTEGER OPTIMAL SOLUTION FOUND
> cat("status:",result$status)
status: optimal
> cat("objective:",result$value)
objective: 7.42218
> print(result$getValue(Y))
             [,1]        [,2]        [,3]        [,4]        [,5]        [,6]       [,7]       [,8]       [,9]      [,10]
 [1,]  0.92604847 -0.73860862  0.50895032  0.77907863 -0.42484496  0.91366669 -0.5127611  0.3302304 -0.9083377 -0.7144000
 [2,]  0.80459809  0.30620385  0.25844226  0.38560681  0.57661027 -0.09333169  0.3361112 -0.8103187 -0.1155999 -0.1709073
 [3,]  0.38141056 -0.31296706  0.42036480  0.28101363 -0.18204616  0.35514127 -0.1647064 -0.2320607  0.5978497 -0.1725513
 [4,]  0.59093484  0.31351626 -0.99875045  0.98853955  0.76603481  0.14526680  0.5763917 -0.4512327 -0.7562015 -0.2623091
 [5,] -0.95077263 -0.35925352 -0.04936685  0.31141160  0.88093457 -0.79415063 -0.7942707  0.6292801  0.1218960 -0.6951105
 [6,] -0.04440806 -0.62461776 -0.55976223  0.41706094 -0.90888700  0.79964994 -0.1302145 -0.1029673 -0.5869372 -0.7223879
 [7,]  0.51691908  0.56458860 -0.24036692  0.08813205  0.05621098 -0.50782453  0.9699140  0.6201287 -0.7449367 -0.5339318
 [8,] -0.56718413 -0.81281003  0.22554201  0.18828404  0.78483809 -0.91588093  0.7861022  0.6247790  0.5066157 -0.0680751
 [9,] -0.36363798 -0.06644192 -0.29640418 -0.42168053  0.10287003 -0.34415856  0.7729381  0.5886846  0.7900907 -0.4680547
[10,] -0.53674843  0.02301092 -0.77772915 -0.70577271 -0.08677053  0.90900730 -0.6498947 -0.1203366 -0.2510744  0.7156554
> 

It helps that we can use the built-in functions matrix_trace and sum_entries. A CVXR model in pure matrix algebra would be more difficult to read. A model in pure matrix algebra could look like:

Model in Matrix Algebra
\[\begin{align}\max\>&\mathbf{tr}(\color{darkred}Y) \\ & \color{darkred}Y = \color{darkblue}A \color{darkred}P \\ & \color{darkred}P^T \color{darkblue}e = \color{darkblue}e \\ & \color{darkred}P\color{darkblue}e = \color{darkblue}e \\ & \color{darkred}P \in \{0,1\}^{n \times n} \\ & \color{darkred}Y \>\mathbf{free}\end{align}\]

where \(e\) is a column of ones.

CVXR (and its Python cousin CVXPY) often works very nicely when working with 1d and 2d data structures: vectors and matrices. Unfortunately, it does not generalize to higher-dimensional problems. Many practical problems need more than 2 indices, but here we are lucky: this is a matrix problem where CVXR is at its best.

It is noted that this size (\(n=10\)) can still be handled by complete enumeration. There are \(n!=3,628,800\) column permutations.

Optimize the formulation

It is possible to optimize the formulation. First, we can observe that not all \(y\)'s are used. Only the diagonal elements are of interest. We can remove all off-diagonal variables. Second, we can substitute out the remaining \(y\) variables.

This would lead to:

Mixed Integer Programming Model v2
\[\begin{align}\max&\sum_i \sum_k \color{darkblue}a_{i,k} \cdot \color{darkred}p_{k,i} \\ & \sum_j \color{darkred}p_{i,j} = 1 &&\forall i \\ & \sum_i \color{darkred}p_{i,j} = 1 &&\forall j \\ & \color{darkred}p_{i,j} \in \{0,1\} \end{align}\]

In practical models, I am usually not so worried about a bunch of extra continuous variables. Often I leave them in, so I can inspect their values if I need to debug the model. In the implementations above I printed out the \(y\) variables to observe how columns were permuted.

This is indeed a standard assignment problem. We can relax the integer variables (they will be integer-valued automatically) and also use specialized algorithms. Note that the solver logs from the previous model also indicate the MIP solvers can solve that problem without doing any branching.

References

1. Maximise diagonal of matrix by permuting columns in R, https://stackoverflow.com/questions/61565176/maximise-diagonal-of-matrix-by-permuting-columns-in-r
2. Permutation matrix, https://en.wikipedia.org/wiki/Permutation_matrix

Random Sparse Arcs in GAMS

I was playing a bit with generating a large network (graph) with \(n=5000\) nodes and say \[\frac{n^2}{100}= 250,000\] arcs. The standard way to generate this in GAMS is:

* * nodes * set i 'nodes' /node1*node5000/; alias(i,j); * * random arcs (1%) * set A(i,j) 'arcs'; A(i,j) = uniform(0,1)<0.01; scalar numArcs 'number of arcs'; numArcs = card(A); display numArcs;

Basically, this approach generates \(n^2\) random numbers  \( r_{i,j} \sim U(0,1)\) (uniform distribution) and picks the arcs related to a value \(r_{i,j}\lt 0.01\). This is a randomized process, so we don't see exactly 250,000 arcs, but rather:

----     16 PARAMETER numArcs              =   249170.000  number of arcs

This operation is not exactly cheap. It takes 11.2 seconds on my laptop.

There are a few other approaches we can try.

Random locations

Instead of generating \(n^2\) random numbers, we can generate \(k=250,000\) pairs of \(i,j\) values (integers between 1 and \(n\)). A straigthforward implementation would be:

scalars k,n,nn,ni,nj; n = card(i); nn = n*n/100; for (k = 1 to nn,    ni = uniformint(1,n);    nj = uniformint(1,n);    A(i,j)$(ord(i)=ni and ord(j)=nj) = yes; );

Unfortunately, this is extremely slow. I stopped it after 2,000 seconds: it was still not finished. GAMS is horribly slow when doing loops like this.

Crazy indexing

A variant of the above approach is to try to trick GAMS into a vectorized version of this. Well, this is not so easy. Here is a version that uses leads/lags when indexing:

sets    A(*,*) 'arcs'    k /node1*node5000,k5001*k250000/    ij /i,j/ ; parameter r(k,ij) 'random offset from k'; r(k,ij) = uniformint(1,card(i)) - ord(k); A(k+r(k,'i'),k+r(k,'j')) = yes;

This requires some explanation. 

• The set A is no longer domain-checked. We use funky indexing here, so we need loosen things up.
• The set k has 250,000 elements. The first 5000 are identical to \(i,j\).
• The parameter r contains \(2\times 250,000\) random integers between 1 and \(n\). We store them as offsets from the index \(k\). This sounds crazy but the next statement explains why this is done.
• The assignment to A is using leads. Think of it as a variant of A(k,k) = yes.  This would populate a diagonal. Similarly, A(k+1,k) = yes shifts the diagonal by one. Using A(k+r(k,'i'),k+r(k,'j')) we index exactly the correct location.
• We may generate duplicates, so the number will be slightly below 250k arcs.

This beast actually works, but the performance is not great. The fragment takes about 13 seconds, so no gain compared to our original approach.

Python loop

Let's try some Python. We can build up sets using Python within a GAMS model:

set A(i,j) 'arcs'; $onEmbeddedCode Python: from random import seed,randint seed(999) n = len(list(gams.get("i"))) nn = n*n//100 a = {} for k in range(nn):    ni = randint(1,n)    nj = randint(1,n)    elem = ("node{}".format(ni),"node{}".format(nj))    a[elem] = 1 gams.set("A",list(a)) $offEmbeddedCode A

This is just a k-loop. We use here a Python dictionary to store elements as this handles possible duplicates correctly. This approach also takes about 11 seconds, so no gain compared to our first approach.

R to the rescue

So let's see if we can use R to speed things up:

$set n    5000 $set inc  data.inc * * R script * $onecho > script.R library(data.table) n <- %n% nn <- n^2/100 df <- data.frame(          ni=sample(n,nn,replace=TRUE),          nj=sample(n,nn,replace=TRUE)) df <- df[order(df$ni,df$nj),] v <- unique(paste0("node",df$ni,".node",df$nj)) fwrite(list(v),"%inc%",col.names=F,quote=F) $offecho $call '"c:\program files\R\R-3.5.0\bin\Rscript.exe" script.R' * * nodes and arcs * set i 'nodes' /node1*node%n%/; alias(i,j); set A(i,j) 'arcs' / $offlisting $include %inc% $onlisting /;

We write from inside GAMS an R script that gets executed by rscript.exe. The steps are:

• sample the arcs,
• sort in the order that GAMS likes,
• create a string representation of the tuple, e.g. "node1.node2",
• remove duplicates (this means we end up with a number of arcs that is smaller than 250k),
• write to an include file (note: the function fwrite from the data.table package seems a bit faster than the function writeLines from the base package)

The include file is then processed by GAMS. To speed things up we remove echoing to the listing file. 

This method is the winner: it takes just 3.5 seconds. Sometimes using a plain old text file as a communication channel is not so bad.

GDX file is slower

The code:

$set n    5000 $set gdx  data.gdx * * R script * $onecho > script.R library(gdxrrw) library(dplyr) n <- 5000 nn <- n^2/100 df <- data.frame(   ni=paste0("node",sample(n,nn,replace=TRUE)),   nj=paste0("node",sample(n,nn,replace=TRUE)),   stringsAsFactors = F) df <- distinct(df,ni,nj) nr <- nrow(df) wgdx("%gdx%",list(name="A", type="set", dim=2, form="sparse", ts="arcs", val=cbind(1:nr,1:nr), uels=c(list(df$ni),list(df$nj)))) $offecho $call '"c:\program files\R\R-3.5.0\bin\Rscript.exe" script.R' * * nodes and arcs * set i 'nodes' /node1*node%n%/; alias(i,j); set A(i,j) 'arcs'; $gdxin %gdx% $loaddc A

is slower than using a text file. The reason is that wgdx is not as efficient.

Mondriaan Tilings 2

In [1], I showed a Mixed Integer Programming formulation for the Mondriaan Tiling Problem. To summarize, the problem can be stated as:

1. Given an \(n\times n\) square area,
2. Partition this area into differently sized rectangles,
3. Such that the difference in size between the largest and the smallest rectangle is minimized.

Instead of partitioning, we can look at the problem as a tiling problem: given a set of differently sized tiles, place a subset of them on the area such that the whole area is covered and the difference in the size of used tiles (largest minus smallest) is minimized.

The model in [1] looked like:

Original MIP Model

\[\begin{align}\min\> & \color{darkred}u - \color{darkred}\ell\\ & \sum_{k,r,i,j|\color{darkblue}{\mathit{cover}}(k,r,i,j,i',j')}\color{darkred}x_{k,r,i',j'}=1 && \forall i',j'&&\text{cover cell $(i',j')$}\\ &\color{darkred}y_k = \sum_{r,i,j|\color{darkblue}{\mathit{ok}}(k,r,i,j)}\color{darkred}x_{k,r,i,j}&& \forall k && \text{select max one type of rectangle} \\ & \color{darkred}\ell \le \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k + \color{darkblue}M (1-\color{darkred}y_k) && \forall k && \text{bound smallest area}\\ & \color{darkred}u \ge \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k  && \forall k && \text{bound largest area}\\ & \sum_k \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k = \color{darkblue} n^2 && && \text{extra constraint}\\ & \color{darkred}x_{k,r,i,j} \in \{0,1\} \\&  \color{darkred}y_k \in \{0,1\} \\ & \color{darkred}u, \color{darkred}\ell \ge 0 \\& \color{darkblue}M = \max_k \color{darkblue}{\mathit{area}}_k\end{align}\]

The main (binary) variables are:

Decision Variables
\[\color{darkred} y_k = \begin{cases} 1 & \text{if tile $k$ is selected} \\ 0 & \text{otherwise} \end{cases}\]
\[\color{darkred} x_{k,r,i,j} = \begin{cases} 1 & \text{if tile $k$ (orientation $r$) is placed at location $(i,j)$} \\ 0 & \text{otherwise} \end{cases}\]

Here I want to discuss a different approach. Instead of solving the problem in one swoop for both \(y\) (select tiles) and \(x\) (place tiles), we can design a simple algorithm that iterates between two different models until we find a feasible integer solution. By ordering things from best to worse, as soon as we find a feasible solution, we have the optimal solution.

The first model is dedicated to select \(y\) variables (i.e. select tiles) such that \[ \sum_k \mathit{area}_k \cdot y_k = n^2\] and where the objective \[ \min \> u-\ell\] is optimized. The second model takes a solution for \(y\) and checks if it is feasible. If not feasible, we go back to model 1 and generate a new, next best solution for \(y\).

Model 1: selection of tiles, can look like:

MIP Model 1: Generate proposals for \(y\)

\[\begin{align}\min\> & \color{darkred}u - \color{darkred}\ell\\ & \color{darkred}\ell \le \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k + \color{darkblue}M (1-\color{darkred}y_k) && \forall k && \text{bound smallest area}\\ & \color{darkred}u \ge \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k  && \forall k && \text{bound largest area}\\ & \sum_k \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k = \color{darkblue} n^2 && && \text{area constraint}\\ & \sum_k (2 \color{darkblue}Y^*_{p,k} -1) \color{darkred}y_k \le \sum_k \color{darkblue}Y^*_{p,k}-1 && p \in P && \text{cut off previously found solutions} \\ &  \color{darkred}y_k \in \{0,1\} \\ & \color{darkred}u, \color{darkred}\ell \ge 0 \\& \color{darkblue}M = \max_k \color{darkblue}{\mathit{area}}_k\end{align}\]

Previously found solutions for \(y\) are stored in the data structure \(Y^*\). The cut will forbid any previously found solution to be considered again. In the first solve, we can use \(P=\emptyset\) (no previously found solutions).

The second model establishes whether the solution \(y^*\) in model 1 is feasible. This model, given a selection of tiles (\(y^*\)) will try to place the tiles such that all locations of the whole area are covered by exactly one tile.

MIP Model 2: Check feasibility of \(y^*\)

\[\begin{align}\min\> & 0 \\ & \sum_{k^*,r,i,j|\color{darkblue}{\mathit{cover}}(k^*,r,i,j,i',j')}\color{darkred}x_{k^*,r,i',j'}=1 && \forall i',j' \text{ and } k^* \in \color{darkblue} K^* &&\text{cover cell $(i',j')$}\\ &\sum_{r,i,j|\color{darkblue}{\mathit{ok}}(k^*,r,i,j)}\color{darkred}x_{k^*,r,i,j} = 1&& \text{where } k^* \in \color{darkblue} K^* && \text{select one type of rectangle}\\ & \color{darkred}x_{k^*,r,i,j} \in \{0,1\}  \\ & \color{darkblue} K^* = \{k|\color{darkblue}y_k^*=1\} && &&\text{set of selected tiles}\end{align}\]

The algorithm can look like:

1. Initialization: \(P = \emptyset\).
2. Solve model 1
3. Record the solution \(y^*_k\), augment the set \(P\) and add this solution to \(Y^*_{p,k}\).
4. Solve model 2
5. If feasible: STOP. We have found the optimal solution.
6. Go to step 2.

Small example

When we solve a \(10\times 10\) problem, we need 10 cycles:

----    227 PARAMETER ysolrep  solutions for m1 (ordered by obj)

              obj    smallest     largest      status

sol1            3          24          27  Infeasible
sol2            6          18          24  Infeasible
sol3            6          18          24  Infeasible
sol4            7           9          16  Infeasible
sol5            7           9          16  Infeasible
sol6            7          21          28  Infeasible
sol7            7          21          28  Infeasible
sol8            7          14          21  Infeasible
sol9            7          14          21  Infeasible
sol10           8           4          12    Feasible

The status column indicates what model 2 thinks about this proposal. The first 9 generated proposals were infeasible. The first tiling with a difference of 8 was feasible, and we can declare that an optimal solution. We have proven that no better solution exists.

The results look like:

Solution for 10x10 problem

Large example

The results for a \(25\times 25\) problem are:

Optimal solution for 25x25 problem

This problem required 246 cycles (i.e. 496 models to be solved).

References

1. Mondriaan Tilings, https://yetanothermathprogrammingconsultant.blogspot.com/2020/04/mondriaan-tilings.html The problem and a mathematical model are detailed in this post. 

On the importance of well-written error messages

This is funny:

"I also tried to get Fortran programs running on Imperial's IBM 7010, since I sort of knew Fortran, certainly better than I knew Cobol, and Fortran would have been better suited for data analysis. It was only after weeks of fighting JCL, IBM's Job Control Language, that I deduced that there was no Fortran compiler on the 7010, but the JCL error messages were so inscrutable that no one else had figured that out either." [1]

Indeed, even today error messages are often poorly formulated and leaving users in the dark.

IBM 7010

References

1. Brian Kernighan, UNIX, A History and a Memoir, 2020. 

Mondriaan Tilings

Introduction

Piet Mondriaan (later Mondrian) (1872-1944) was a Dutch painter, famous for his abstract works [1]. 

A tiling problem is related to this [2]:

The problem can be stated as follows:

• Given is a square \((n \times n)\) area,
• Partition this square into multiple, integer-sided, non-congruent rectangles,
• Minimize the difference in size between the largest and the smallest rectangle.

Non-congruent means: only one \((p \times q)\) or \((q\times p)\) rectangle can be used.

It is interesting to make some predictions.

• We try not to use a combination of very small and very large rectangles. That is obvious.
• Only using small rectangles will always lead to a larger difference, as we need many small rectangles to cover the complete square. This will lead automatically to a larger range of sizes.
• Very large rectangles are probably no good either: high probability we need some very small rectangles to fill up any remaining holes.
• So I guess, we will use medium-sized rectangles of similar size. 

Example

Consider a small problem with \(n=5\). The possible sides for each rectangle are \(1,\dots,n\). From this, we can form a set of possible rectangles (I automated this step):

----     65 SET s  sizes of rectangles (height or width)

s1,    s2,    s3,    s4,    s5


----     65 SET k  rectangle id

k1 ,    k2 ,    k3 ,    k4 ,    k5 ,    k6 ,    k7 ,    k8 ,    k9 ,    k10,    k11,    k12,    k13,    k14


----     65 SET rec  rectangles

                s1          s2          s3          s4

k1 .s1         YES
k2 .s2         YES
k3 .s2                     YES
k4 .s3         YES
k5 .s3                     YES
k6 .s3                                 YES
k7 .s4         YES
k8 .s4                     YES
k9 .s4                                 YES
k10.s4                                             YES
k11.s5         YES
k12.s5                     YES
k13.s5                                 YES
k14.s5                                             YES

We see there are 14 different rectangles we can use. Note that the rectangle with size \((5 \times 5)\) is excluded. It would completely fill the area, while we need to place multiple rectangles. Also note: we only have \((p \times q)\) but not \((q \times p)\).  We cannot use both such rectangles as they are congruent. We can rotate a rectangle when placing it.

Optimal solution

An optimal solution for this problem is:

These are rectangles \(k5, k6, k12\) with \(k5\) being rotated. Obviously, this solution is not unique.

Optimization model

To model this, I use a binary variable: \[ x_{k,r,i,j} = \begin{cases} 1 & \text{if rectangle $k$ (rotated if $r=\mathrm{y}$) is placed at location $(i,j)$} \\ 0 & \text{otherwise}\end{cases}\] Here \(r = \{\mathrm{n},\mathrm{y}\}\) indicates if a rectangle is rotated. Note that not all rectangles can be rotated (it does not make sense to rotate a square, or there may not be space left to rotate for a given position).

For convenience, I also introduce \[y_{k} = \begin{cases} 1 & \text{if rectangle $k$ is used} \\  0 & \text{otherwise}\end{cases}\] This is essentially an aggration of \(x\).

In addition, I use a few more complex data structures: \[\mathit{ok}_{k,r,i,j} = \begin{cases} \mathit{Yes} & \text{if we can place rectangle $k$ at location $(i,j)$} \\ \mathit{No} & \text{otherwise}\end{cases}\] Again, \(r\) indicates if rectangle \(k\) is rotated. I implemented this as a GAMS set. You can consider it as a Boolean data structure. Furthermore: \[\mathit{cover}_{k,r,i,j,i',j'} = \begin{cases} \mathit{Yes} & \text{if cell $(i',j')$ is covered by placing rectangle $k$ at location $(i,j)$} \\ \mathit{No} & \text{otherwise}\end{cases}\]

These sets are large. The set ok looks like:

----    103 SET ok  usable position

INDEX 1 = k1

              j1          j2          j3          j4          j5

n.i1         YES         YES         YES         YES         YES
n.i2         YES         YES         YES         YES         YES
n.i3         YES         YES         YES         YES         YES
n.i4         YES         YES         YES         YES         YES
n.i5         YES         YES         YES         YES         YES

INDEX 1 = k2

              j1          j2          j3          j4          j5

n.i1         YES         YES         YES         YES         YES
n.i2         YES         YES         YES         YES         YES
n.i3         YES         YES         YES         YES         YES
n.i4         YES         YES         YES         YES         YES
y.i1         YES         YES         YES         YES
y.i2         YES         YES         YES         YES
y.i3         YES         YES         YES         YES
y.i4         YES         YES         YES         YES
y.i5         YES         YES         YES         YES

...

INDEX 1 = k13

              j1          j2          j3

n.i1         YES         YES         YES
y.i1         YES
y.i2         YES
y.i3         YES

INDEX 1 = k14

              j1          j2

n.i1         YES         YES
y.i1         YES
y.i2         YES

The set cover is even larger:

----    103 SET cover  coverage of cells

INDEX 1 = k1  INDEX 2 = n  INDEX 3 = i1  INDEX 4 = j1

            j1

i1         YES

INDEX 1 = k1  INDEX 2 = n  INDEX 3 = i1  INDEX 4 = j2

            j2

i1         YES

...

INDEX 1 = k14  INDEX 2 = y  INDEX 3 = i1  INDEX 4 = j1

            j1          j2          j3          j4          j5

i1         YES         YES         YES         YES         YES
i2         YES         YES         YES         YES         YES
i3         YES         YES         YES         YES         YES
i4         YES         YES         YES         YES         YES

INDEX 1 = k14  INDEX 2 = y  INDEX 3 = i2  INDEX 4 = j1

            j1          j2          j3          j4          j5

i2         YES         YES         YES         YES         YES
i3         YES         YES         YES         YES         YES
i4         YES         YES         YES         YES         YES
i5         YES         YES         YES         YES         YES

The reason to develop these sets is to make the equations simpler. In addition, we can develop, test and debug the calculation of these sets in isolation of the test of the model and before the equations are written.

With this, we can formulate:

MIP Model

\[\begin{align}\min\> & \color{darkred}u - \color{darkred}\ell\\ & \sum_{k,r,i,j|\color{darkblue}{\mathit{cover}}(k,r,i,j,i',j')}\color{darkred}x_{k,r,i',j'}=1 && \forall i',j'&&\text{cover cell $(i',j')$}\\ &\color{darkred}y_k = \sum_{r,i,j|\color{darkblue}{\mathit{ok}}(k,r,i,j)}\color{darkred}x_{k,r,i,j}&& \forall k && \text{select max one type of rectangle} \\ & \color{darkred}\ell \le \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k + \color{darkblue}M (1-\color{darkred}y_k) && \forall k && \text{bound smallest area}\\ & \color{darkred}u \ge \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k  && \forall k && \text{bound largest area}\\ & \sum_k \color{darkblue}{\mathit{area}}_k \cdot\color{darkred}y_k = \color{darkblue} n^2 && && \text{extra constraint}\\ & \color{darkred}x_{k,r,i,j} \in \{0,1\} \\&  \color{darkred}y_k \in \{0,1\} \\ & \color{darkred}u, \color{darkred}\ell \ge 0 \\& \color{darkblue}M = \max_k \color{darkblue}{\mathit{area}}_k\end{align}\]

When solving this model, we see:

----    146 VARIABLE x.L  place tile

                  j1          j4

k5 .y.i4           1
k6 .n.i1           1
k12.n.i1                       1


----    149 VARIABLE y.L  tile is used

k5  1,    k6  1,    k12 1


----    157 PARAMETER v  tile numbers 

            j1          j2          j3          j4          j5

i1           6           6           6          12          12
i2           6           6           6          12          12
i3           6           6           6          12          12
i4           5           5           5          12          12
i5           5           5           5          12          12


----    159 VARIABLE smallest.L            =        6.000  smallest used tile
            VARIABLE largest.L             =       10.000  largest used tile
            VARIABLE z.L                   =        4.000  objective

It is noted that these types of formulations lead to very large models [5].

One remaining question is whether we can use some symmetry-breaking constraint to speed up solution times.

A larger problem

It looks like \(n=19\) is interesting [3]. This is about the limit that we can solve comfortably (a few hours to global optimality) with a MIP solver. This model had about 36,000 binary variables.

An optimal 19x19 problem

The colors are not part of the model: they were added by hand to make the picture more Mondrianic.

A list of optimal objective values can be found in [4].

References

1. Piet Mondrian, https://en.wikipedia.org/wiki/Piet_Mondrian 
2. Mondrian Puzzle -- Numberphile, https://www.youtube.com/watch?v=49KvZrioFB0
3. Brady Haran, Mondrian Art Puzzle, https://www.bradyharanblog.com/blog/mondrian-art-puzzle
4. A276523 in The Online Encyclopedia of Sequences, https://oeis.org/A276523
5. Tiling, https://yetanothermathprogrammingconsultant.blogspot.com/2020/03/tiling.html This post shows a comparison of the above formulation with a "no-overlap" formulation for a related problem.
6. Mondriaan Tilings 2, http://yetanothermathprogrammingconsultant.blogspot.com/2020/04/mondriaan-tilings-2.html. This post demonstrates an alternative approach to this problem.