In this post, I want to discuss an observation about the root node when solving a MIP model.
Problem Description
We have a large matrix \(\color{darkblue}A=\color{darkblue}a_{i,j}\) with values 0 and 1. In addition, we have a minimum on the row and column totals. These are called \(\color{darkblue}r_i\) and \(\color{darkblue}c_j\). The goal is to remove as many 1's in the matrix \(\color{darkblue}A\) subject to these minimum row and column totals.
Mathematical Model
Now, we can write:
Mathematical Model 

\[ \begin{align} \min& \sum_{i,j\color{darkblue}a_{i,j}=1} \color{darkred}x_{i,j} \\ & \sum_{j\color{darkblue}a_{i,j}=1} \color{darkred}x_{i,j} \ge \color{darkblue}r_i && \forall i \\ & \sum_{i\color{darkblue}a_{i,j}=1} \color{darkred}x_{i,j} \ge \color{darkblue}c_j && \forall j \\ & \color{darkred}x_{i,j} \in \{0,1\} \end{align}\] 
If you prefer, we can also write this as:
Mathematical Model 

\[ \begin{align} \min& \sum_{i,j} \color{darkblue}a_{i,j} \cdot \color{darkred}x_{i,j} \\ & \sum_j \color{darkblue}a_{i,j}\cdot \color{darkred}x_{i,j} \ge \color{darkblue}r_i && \forall i \\ & \sum_i \color{darkblue}a_{i,j}\cdot \color{darkred}x_{i,j} \ge \color{darkblue}c_j && \forall j \\ & \color{darkred}x_{i,j} \in \{0,1\} \end{align}\] 
The post in [1] mentions that the goal is to solve this model with a matrix size of 20,000 rows and 500 columns. As I used a fraction of 10% ones in the matrix \(\color{darkblue}A\), the number of (discrete) variables in the model is close to million: 999,469 (we could have predicted this from \(20,000 \times 500 \times 0.1\)). The number of constraints is 20,500.
Observation
 All matrix coefficients are \(1\) or \(+1\) (check),
 the righthandside should be integer valued (check),
 there are two nonzero coefficients in each column of the constraint matrix (check),
 and these two nonzero coefficients have the value \(1\) and \(+1\).
 88 PARAMETER results MIP LP NETWORK status Optimal Optimal Optimal obj 504711.000 504711.000 504711.000 time 5141.515 14.437 4.000 iter 571595.000 38276.000 nodes NA NA
Iteration log . . . Iteration: 1 Dual objective = 1079.000000 Perturbation started. Iteration: 606 Dual objective = 499871.000000 Iteration: 2167 Dual objective = 499871.008110 Iteration: 3647 Dual objective = 499871.016553 . . . Iteration: 33697 Dual objective = 504711.126635 Iteration: 33962 Dual objective = 504711.126636 Removing perturbation. Primal simplex solved model. Root relaxation solution time = 16.92 sec. (4357.36 ticks) Nodes Cuts/ Node Left Objective IInf Best Integer Best Bound ItCnt Gap * 0+ 0 999469.0000 0.0000 100.00% Found incumbent of value 999469.000000 after 20.62 sec. (6918.29 ticks) * 0+ 0 522140.0000 0.0000 100.00% Found incumbent of value 522140.000000 after 5140.89 sec. (6925.92 ticks) * 0 0 integral 0 504711.0000 504711.0000 571595 0.00% Elapsed time = 5140.97 sec. (7024.08 ticks, tree = 0.00 MB) Found incumbent of value 504711.000000 after 5140.97 sec. (7024.08 ticks) Root node processing (before b&c): Real time = 5140.98 sec. (7037.52 ticks) Parallel b&c, 16 threads: Real time = 0.00 sec. (0.00 ticks) Sync time (average) = 0.00 sec. Wait time (average) = 0.00 sec.  Total (root+branch&cut) = 5140.98 sec. (7037.52 ticks)
Conclusion
Update: some other solvers do not seem to have this extreme behavior. They solve the MIP and LP models in about the same time. It looks like the poor performance on the MIP model is a particularity of Cplex.
References
 Efficiently remove a maximum amount of binary elements while keeping row and column sums above a certain level, https://stackoverflow.com/questions/76105697/efficientlyremoveamaximumamountofbinaryelementswhilekeepingrowandcol
Appendix: GAMS model
$onText
A very large MIP. It takes some time to solve. However the LP is actually automatically integer valued. This means the relaxation is no good? Finally we also solve as a network model.
Reference: https://stackoverflow.com/questions/76105697/efficientlyremoveamaximumamountofbinaryelementswhilekeepingrowandcol
$offText
* * Data *
set i /i1*i20000/ j /j1*j500/ ;
set a(i,j) 'randomly filled with 10% ones'; a(i,j)$(uniform(0,1)<0.1) = YES;
* calculate rows and column sums parameter rowsum(i), colsum(j); rowsum(i) = sum(a(i,j),1); colsum(j) = sum(a(i,j),1);
parameter r(i),c(j); r(i) = ceil(rowsum(i)/2); c(j) = ceil(colsum(j)/2);
* * MIP Model *
binary variable x(i,j) 'only used when a(i,j)=1'; variable z 'objective';
equations obj erowsum(i) 'minimum number of ones in row' ecolsum(j) 'minimum number of ones in column' ;
obj.. z =e= sum(a,x(a)); erowsum(i).. sum(a(i,j),x(i,j)) =g= r(i); ecolsum(j).. sum(a(i,j),x(i,j)) =g= c(j);
model m /all/;
* * Reporting macro *
acronym Optimal;
parameter results; $macro report(m,label) \ results('status',label) = m.modelstat; \ results('status',label)$(m.modelstat=1) = Optimal; \ results('obj',label) = z.l; \ results('frac',label) = sum((i,j)$(x.l(i,j)>0.001 and x.l(i,j)<0.999),1); \ results('time',label) = m.resusd; \ results('iter',label) = m.iterusd; \ results('nodes',label) = m.nodusd;
* * solve as MIP, LP and Network model *
m.solprint = %solprint.Silent%; option threads = 0, bratio = 1;
solve m minimizing z using mip; report(m,'MIP')
solve m minimizing z using rmip; report(m,'LP')
m.optfile=1; solve m minimizing z using rmip; report(m,'NETWORK')
display results;
* * network option file, for use with third solve *
$onecho > cplex.opt lpmethod 3 $offecho

Appendix Highs open source solver
C:\tmp\HiGHS\build\RELEASE\bin>highs parallel on solver ipm
\tmp\test\Free.mps 
Just I guess, but I suspect that the CPLEX does a bunch of presolve operations before it even attempts to solve the root LP. That would account, at least in part, for the time differential, and I think it might also account for the zero lower bound. Note that the log lines with zero lower bound also report zero (blank) iteration count. I'm guessing the zero lower bound came from some very basic inference by the presolver (the sum of a bunch of binary variables cannot go below 0) and predated any attempt to solve the root LP.
ReplyDeleteThis could support that:
DeleteParallel mode: deterministic, using up to 16 threads.
Parallel mode: deterministic, using up to 3 threads for parallel tasks at root LP.
Parallel mode: deterministic, using up to 13 threads for concurrent optimization
So 3 threads are devoted for doing things other than solving the LP (via concurrent optimization). But still. the LP takes only 17 seconds. They could check integrality in the 18th second and just kill running these expensive (5000 seconds) heuristics. Could be that these heuristics are usually cheap, so normally not a problem to wait for them to finish.
This comment has been removed by the author.
ReplyDeleteI have seen that Cplex seems to be stuck for a very long time after the root LP has been solved many times. What often helps me in those cases is to disable "auxiliary tasks" at the root node. I am not sure what those tasks are exactly. It was once told to me by Cplex support the they include "extra heuristics running at root node only". One can disable those tasks via GAMS/CPLEX option auxrootthreads=1. The line
ReplyDeleteParallel mode: deterministic, using up to 3 threads for parallel tasks at root LP.
should then disappear and for me the MIP solves a bit faster than the LP (33 sec. vs 42 sec.).
Good to know. You may want to share this problem with the Cplex people as it is embarrasing enough to warrant some attention.
Delete