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

The output of this model can look like:

---- 40 SETii1, i2, i3 ---- 40 SETjj1, j2, j3, j4, j5, j6 ---- 40 PARAMETERai1 1.000, i2 2.000, i3 3.000 ---- 40 PARAMETERbj1 4.000, j2 9.000, j3 5.000, j4 3.000, j5 2.000, j6 10.000 ---- 40 VARIABLEz.L = 16.000objective---- 40 VARIABLEx.Lassignmentj1 j4 j5 i1 1.000 i2 1.000 i3 1.000 ---- 40 VARIABLEv.Lposition of a(i) in bi1 1.000, i2 4.000, i3 5.000

- We can tighten the bounds on variable \(\color{darkred}v_i\) to \(\color{darkred}v_i \in [i,\color{darkblue}n-\color{darkblue}m+i]\).
- The variable \(\color{darkred}v_i\) is automatically integer-valued. We can declare it as an integer or continuous variable.
- If declared as a continuous variable, Cplex may make it an integer variable. We can see this in the log: Reduced MIP has 17737 binaries, 99 generals, 0 SOSs, and 0 indicators.
- A larger problem with \(\color{darkblue}m=50, \color{darkblue}n=500\) solves in about 50 seconds.

Assignment problem as LP | Assignment problem as MIP | Full MIP model | |
---|---|---|---|

Rows/Columns | 1,100/100,000 | 1,100/100,000 | 1,300/100,100 |

Objective | 11,524.178 | 11,524.178 | 14,371.455 |

Time (seconds) | 1 | 7 | 1,853 |

Nodes | - | 0 | 3,419 |

Iterations | 2,672 | 2,727 | 30,7429 |

**assignment problem**was formed by dropping the constraints involving \(\color{darkred}v_i\). Further, the problem was solved as an LP. It is noted that Cplex behaves differently when we use continuous variables \(\color{darkred}x_{i,j}\) or declare them as binary variables, When they are binary variables, the presolver works differently, and as a result, the solution of the relaxed problem is no longer optimal: it has to do some work to find the optimal integer solution. In this case, this was done in MIP preprocessing: no branching was needed. When solved as a pure LP, the optimal solution is integer-valued automatically. There is a difference in solution time: the LP solves in 1 second, while the MIP version takes 7 seconds.

Assignment problem as LP | Assignment problem as MIP | Full MIP model | Using initial solution | |
---|---|---|---|---|

Rows/Columns | 1,100/100,000 | 1,100/100,000 | 1,300/100,100 | 1,300/100,100 |

Objective | 11,524.178 | 11,524.178 | 14,371.455 | initial: 17,739.702 optimal: 14,371.455 |

Time (seconds) | 1 | 7 | 1,853 | 1,458 |

Nodes | - | 0 | 3,419 | 3,546 |

Iterations | 2,672 | 2,727 | 307,429 | 343,969 |

- Select 100 smallest values of \(\color{darkblue}b_j\). Assign \(\color{darkblue}a_i\) in order. This gives an initial solution of 17,739.702.
- Pick 200 of of the smallest values of \(\color{darkblue}b_j\). Solve our MIP problem on this subset. This gives an optimal solution of 14,451.861 in 118 seconds.
- Use this as an initial solution for the whole problem. This now solves in 972 seconds. Total solution time is 118+972=1,090 seconds.

Assignment problem as LP | Assignment problem as MIP | Full MIP model | Using initial solution | Using two models | |
---|---|---|---|---|---|

Rows/Columns | 1,100/100,000 | 1,100/100,000 | 1,300/100,100 | 1,300/100,100 | subset: 500/20,100 full: 1,300/100,100 |

Objective | 11,524.178 | 11,524.178 | 14,371.455 | initial: 17,739.702 optimal: 14,371.455 | initial: 17,739.702 subset: 14,451.861 full: 14,371.455 |

Time (seconds) | 1 | 7 | 1,853 | 1,458 | subset: 118 full: 972 |

Nodes | - | 0 | 3,419 | 3,546 | subset: 3,205 full: 3,500 |

Iterations | 2,672 | 2,727 | 307,429 | 343,969 | subset: 276,796 full: 456,537 |

#### Conclusions

- Sometimes it makes sense to replace one big solve with a number of solves. Each one will start from and improve a previously found solution.
- In production environments, a setup like this has another advantage: if a single model fails, you still get good solutions.
- In the above table, the statistics for the number of nodes and Simplex iterations are not very indicative of the total effort. As MIP solvers become more complex, the best performance measure is solution time.
- An assignment problem (or other network problems) should be solved as an LP, not as a MIP.

#### References

- Minimize sum of product of two uneven consecutive arrays, https://stackoverflow.com/questions/65884107/minimize-sum-of-product-of-two-uneven-consecutive-arrays
- Assignment problem with a wrinkle formulated as a network problem, https://yetanothermathprogrammingconsultant.blogspot.com/2021/01/assignment-problem-with-wrinkle.html. Write-up on how to formulate this problem as a shortest path problem.

#### Appendix: GAMS model

$ontext 1. Find
a trivial initial solution (100 smallest values in b)2. Find
a better solution using 200 smallest values3.
Solve the full problem$offtext *-------------------------------------------------------* data*-------------------------------------------------------setsi /i1*i100/ j /j1*j1000/ ; parameterA(i) B(j) ; a(i) = uniform(0,100); b(j) = uniform(0,100); display a,b;*-------------------------------------------------------* optimization model*-------------------------------------------------------set subset(j) 'subset of b';binary variablex(i,j) 'assign'; variablev(i) 'position of a(i) in b'z 'objective'; * tight bounds on vv.lo(i) = ord(i);v.up(i) = card(j) - (card(i)-ord(i));equationsassign1(i) 'assignment constraint'assign2(j) 'assignment constraint'calcv 'caculate v'order 'ordering constraint'obj 'objective'; obj.. z =e= sum((i,subset(j)),x(i,j)*a(i)*b(j));assign1(i).. sum(subset(j), x(i,j)) =e= 1;assign2(subset(j)).. sum(i, x(i,j)) =l= 1;calcv(i).. v(i) =e= sum(subset(j), ord(j)*x(i,j));order(i-1).. v(i) =g= v(i-1)+1; model m /all/;option optcr=0,
threads=8;*-------------------------------------------------------* set initial values using trivial solution of* 100 smallest values in b.*-------------------------------------------------------* find 100 smallest values in bsubset(j) = no;scalarsbmin 'minimum value in remaining b'k 'loop variable'; for(k=1 to 100,bmin = smin(j$(not subset(j)), b(j));loop(j$(b(j)=bmin),subset(j) = yes;break;); ); * let x reflect this initial subsetloop(i,loop(j$subset(j),x.l(i,j) = 1; subset(j) = 0; break;); ); z.l = sum((i,j),x.l(i,j)*a(i)*b(j));v.l(i) = sum(j, ord(j)*x.l(i,j));option v:0;display "*********** initial",z.l,v.l;*-------------------------------------------------------* limit search to 200 values*-------------------------------------------------------for(k=1 to 200,bmin = smin(j$(not subset(j)), b(j));loop(j$(b(j)=bmin),subset(j) = yes;break;); ); m.optfile=1; solve m minimizing
z using mip;display "*********** solve
1",z.l,v.l;*-------------------------------------------------------* solve full model*-------------------------------------------------------subset(j) = yes;solve m minimizing
z using mip;display "*********** solve
2",z.l,v.l;$onecho > cplex.opt mipstart 1$offecho |

My colleague (and regular commenter on this blog) Rob Pratt suggested a purely graph model:

ReplyDeleteNode(i,j,k), where (j < k) indicates that a_i is assigned to b_j and a_{i+1} is assigned to b_k. The main arcs are (i,j,k)->(i+1,k,\ell) with k < \ell and cost a_{i+1} b_k. The arc from a dummy source node to (1,j,k) picks up the initial cost a_1 b_j + a_2 b_k, and the cost from the last real node (m-1,j,k) to a dummy sink node is 0. The the problem is to find the shortest path from the dummy source to the sink in a directed acyclic graph. For the large model that is 50 million nodes and 2.5 billion arcs. It would take a lot more memory than the milp, but I think it would be faster.

That is truly a large network. I am sure I have never solved an LP of this size, but may be a shortest path algorithm can handle this.

DeleteCan we not just have node(i,j) -> n(i+1,j+k) for k=1,2,... Each incoming (or each outgoing) arc has a cost a(i)/b(j).

DeleteYes, that smaller network seems correct: arc from source to (1,j) with cost a(1)b(j), arc from (i,j) to (i+1,k), where, j < k, with cost a(i+1)b(k), and arc from (m,j) to sink with cost 0.

DeleteThis comment has been removed by the author.

ReplyDeleteI tried this. Seems to work for small instances. The big one is still big. Too large for me to solve on my laptop (90k nodes, 4e7 arcs).

DeleteSolves quickly with Dijkstra on my laptop.

DeleteI tried to pass it on to Cplex as an LP (and then solve it with the network solver). I should feed it directly into a shortest path algorithm.

DeleteYes. Using a network solver would be pretty fast too (for the solve). But, there will be a lot of overhead to first find the network embedding from the matrix, translate the matrix into a graph and then, of course, network simplex is going to be much slower than a pure shortest path algorithm. Network simplex (a flow formulation) is definitely overkill here.

DeleteThat worked quite nicely.

ReplyDeleteNOTE: -----------------------------------------------------------------------------------------

NOTE: -----------------------------------------------------------------------------------------

NOTE: Running OPTNETWORK.

NOTE: -----------------------------------------------------------------------------------------

NOTE: -----------------------------------------------------------------------------------------

NOTE: Reading the links data.

NOTE: Data input used 1.32 (cpu: 5.87) seconds.

NOTE: Building the input graph storage used 1.59 (cpu: 14.16) seconds.

NOTE: The number of nodes in the input graph is 94953.

NOTE: The number of links in the input graph is 44601302.

NOTE: Processing shortest paths problem using 16 threads across 1 machines.

NOTE: Processing the shortest paths problem between 1 source nodes and 1 sink nodes.

Real

Algorithm Sources Complete Time

shortestPath 1 100% 0.31

NOTE: Processing the shortest paths problem used 0.31 (cpu: 0.31) seconds.

NOTE: The Cloud Analytic Services server processed the request in 3.509987 seconds.

NOTE: The data set SASCAS1.OUT has 101 observations and 6 variables.

NOTE: PROCEDURE OPTNETWORK used (Total process time):

real time 3.56 seconds

cpu time 0.04 seconds

Interesting. We seem to match pretty closely on real time. I'm using an ordinary HP desktop, and I think the graph library I'm using is single threaded -- definitely not 16 threads, I don't have that many cores -- so I'm surprised my time is competitive. My graph is also a bit smaller (but not much).

DeleteHi Paul.

DeleteIn this case, the threading is only for reading the input data (from a table of links) and building the graph data structures. That is ~2.9 seconds.

Dijkstra's runs single-threaded here (since there is only one source) and takes 0.31 seconds.

OPTNETWORK: https://go.documentation.sas.com/?cdcId=pgmsascdc&cdcVersion=v_008&docsetId=casnopt&docsetTarget=casnopt_optnet_overview.htm&locale=en

I think since the graph is directed acyclic, you can do better than Dijkstra's (utilizing topological sort). But, I have not implemented that yet. It's on my ToDo list.

DeleteIn any case, the bulk of the time will be building the graph in memory.

Interesting. I have not used these type of algorithms a lot, so this is let's say educational for me.

DeleteI just wrote some Java code to use a purely graph model, but it's a bit more compact -- basically what Erwin has in his comment. For a 100 x 1000 test case, the graph has roughly 90 thousand nodes and 40.2 million arcs. I get the optimal solution in about 3.5 seconds (including both graph construction and finding the shortest path). I ran some smaller examples against the MIP model from this post and confirmed that the objective values are correct. I'm hoping to post details in the next day or so on my blog.

ReplyDeleteAssuming you encoded the graph object directly from the data (rather than creating a list of edges and then reading into a graph library), then you would just want to look at the time for the shortest path solve. I got around 0.3 seconds. I don't expect decent implementations of single s-t Dijkstra's in C/C++/Java to vary that much. So, I would expect your time to be in the ballpark. Unless it exploits acyclic digraph, then maybe it could be a little bit faster.

Delete