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.

The output of this model can look like:

----     40 SET i

i1,    i2,    i3

----     40 SET j

j1,    j2,    j3,    j4,    j5,    j6

----     40 PARAMETER a

i1 1.000,    i2 2.000,    i3 3.000

----     40 PARAMETER b

j1  4.000,    j2  9.000,    j3  5.000,    j4  3.000,    j5  2.000,    j6 10.000

----     40 VARIABLE z.L                   =       16.000  objective

----     40 VARIABLE x.L  assignment

j1          j4          j5

i1       1.000
i2                   1.000
i3                               1.000

----     40 VARIABLE v.L  position of a(i) in b

i1 1.000,    i2 4.000,    i3 5.000


Notes:
• 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.

It is interesting to see what happens when we feed this model an even larger data set. I generated some random data: \begin{align} & \color{darkblue}m=100\\ &\color{darkblue}n=1000 \\ & \color{darkblue}a_i \sim U(0,100) \\ & \color{darkblue}b_j \sim U(0,100)\end{align} This is not a small problem: we will have 100,000 binary variables. This model has no longer a network structure, so compared to a pure assignment problem we see much poorer performance:

Assignment
problem
as LP
Assignment
problem
as MIP
Full MIP
model
Rows/Columns1,100/100,0001,100/100,0001,300/100,100
Objective11,524.17811,524.17814,371.455
Time (seconds)171,853
Nodes-03,419
Iterations2,6722,72730,7429

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

Let's see if we can shave something off this 1,800 seconds solution time. The first experiment I did was to provide an initial solution. I took the 100 smallest values of $$\color{darkblue}b_j$$ and assigned the 100 $$\color{darkblue}a_i$$'s to these while maintaining the order. This is very simple, so this took no time. The MIP was solved with a MIPSTART option. The results are:

Assignment
problem
as LP
Assignment
problem
as MIP
Full MIP
model
Using initial
solution
Rows/Columns1,100/100,0001,100/100,0001,300/100,1001,300/100,100
Objective11,524.17811,524.17814,371.455initial: 17,739.702
optimal: 14,371.455
Time (seconds)171,8531,458
Nodes-03,4193,546
Iterations2,6722,727307,429343,969

This model solves faster. We don't see this in the node and iteration count, mainly because Cplex spends much time in preprocessing and restarts. Our initial objective is not that great. But this leads to the question, can we select the say 200 smallest values of $$\color{darkblue}b_j$$ and solve this smaller problem first. Then do a second solve with use all values of $$\color{darkblue}b_j$$, again using a MIPSTART. So the algorithm becomes:

1. 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.
2. 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.
3. 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.

Adding this to our table gives:

Assignment
problem
as LP
Assignment
problem
as MIP
Full MIP
model
Using initial
solution
Using two
models
Rows/Columns1,100/100,0001,100/100,0001,300/100,1001,300/100,100subset: 500/20,100
full: 1,300/100,100
Objective11,524.17811,524.17814,371.455initial: 17,739.702
optimal: 14,371.455
initial: 17,739.702
subset: 14,451.861
full: 14,371.455
Time (seconds)171,8531,458subset: 118
full: 972
Nodes-03,4193,546subset: 3,205
full: 3,500
Iterations2,6722,727307,429343,969subset: 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

1. Minimize sum of product of two uneven consecutive arrays, https://stackoverflow.com/questions/65884107/minimize-sum-of-product-of-two-uneven-consecutive-arrays
2. 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

This is the model that I used to try out the three-tier approach.

 $ontext Solve ordered assignment problem in three stages: 1. Find a trivial initial solution (100 smallest values in b) 2. Find a better solution using 200 smallest values 3. Solve the full problem$offtext *------------------------------------------------------- * data *------------------------------------------------------- sets   i /i1*i100/   j /j1*j1000/ ; parameter   A(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 variable    x(i,j) 'assign' ; variable    v(i) 'position of a(i) in b'    z    'objective' ; * tight bounds on v v.lo(i) = ord(i); v.up(i) = card(j) - (card(i)-ord(i)); equations     assign1(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 b subset(j) = no; scalars    bmin   '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 subset loop(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

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

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

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

2. Can 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).

3. Yes, 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.

2. This comment has been removed by the author.

1. I 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).

2. Solves quickly with Dijkstra on my laptop.

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

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

3. That worked quite nicely.

NOTE: -----------------------------------------------------------------------------------------
NOTE: -----------------------------------------------------------------------------------------
NOTE: Running OPTNETWORK.
NOTE: -----------------------------------------------------------------------------------------
NOTE: -----------------------------------------------------------------------------------------
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

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

2. Hi Paul.

In 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

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

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

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

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

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