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.

#### Implementation

---- 15 SETii1, i2, i3 ---- 15 SETjj1, j2, j3, j4, j5, j6 ---- 15 PARAMETERai1 1.000, i2 2.000, i3 3.000 ---- 15 PARAMETERbj1 4.000, j2 9.000, j3 5.000, j4 3.000, j5 2.000, j6 10.000

---- 31 SETnnodesj1 j2 j3 j4 j5 j6 src YES i1 YES YES YES YES i2 YES YES YES YES i3 YES YES YES YES snk YES ---- 32 PARAMETERnumnodes= 14.000 ---- 41 SETearcssrc. .i1 .j1, src. .i1 .j2, src. .i1 .j3, src. .i1 .j4, i1 .j1.i2 .j2, i1 .j1.i2 .j3 i1 .j1.i2 .j4, i1 .j1.i2 .j5, i1 .j2.i2 .j3, i1 .j2.i2 .j4, i1 .j2.i2 .j5, i1 .j3.i2 .j4 i1 .j3.i2 .j5, i1 .j4.i2 .j5, i2 .j2.i3 .j3, i2 .j2.i3 .j4, i2 .j2.i3 .j5, i2 .j2.i3 .j6 i2 .j3.i3 .j4, i2 .j3.i3 .j5, i2 .j3.i3 .j6, i2 .j4.i3 .j5, i2 .j4.i3 .j6, i2 .j5.i3 .j6 i3 .j3.snk. , i3 .j4.snk. , i3 .j5.snk. , i3 .j6.snk. ---- 42 PARAMETERnumarcs= 28.000 ---- 47 PARAMETERccost of arcssrc. .i1.j1 4.000, src. .i1.j2 9.000, src. .i1.j3 5.000, src. .i1.j4 3.000, i1 .j1.i2.j2 18.000 i1 .j1.i2.j3 10.000, i1 .j1.i2.j4 6.000, i1 .j1.i2.j5 4.000, i1 .j2.i2.j3 10.000, i1 .j2.i2.j4 6.000 i1 .j2.i2.j5 4.000, i1 .j3.i2.j4 6.000, i1 .j3.i2.j5 4.000, i1 .j4.i2.j5 4.000, i2 .j2.i3.j3 15.000 i2 .j2.i3.j4 9.000, i2 .j2.i3.j5 6.000, i2 .j2.i3.j6 30.000, i2 .j3.i3.j4 9.000, i2 .j3.i3.j5 6.000 i2 .j3.i3.j6 30.000, i2 .j4.i3.j5 6.000, i2 .j4.i3.j6 30.000, i2 .j5.i3.j6 30.000

**source**and the

**sink node**are denoted by \(\color{darkblue}n_{'src',''}\) and \(\color{darkblue}n_{'snk',''}\). Again, the costs of visiting a node are allocated to the incoming arcs. Note that zero costs are not printed (the parameter is stored as a sparse data structure, so zero and does not exist is the same).

**LP model**for this problem can look like:

Shortest Path LP Model |
---|

\[\begin{align}\min& \sum_{i,j,i',j'}\color{darkred}f_{i,j,i',j'}\cdot\color{darkblue}c_{i,j,i',j'}\\ &\sum_{i',j'|e(i',j',i,j)} \color{darkred}f_{i',j',i,j} + \color{darkblue}g_{i,j} = \sum_{i',j'|e(i,j,i',j')} \color{darkred}f_{i,j,i',j'} &&\forall \color{darkblue}n_{i,j}\\ & \color{darkred}f_{i,j,i',j'} \in [0,1] \end{align}\] |

where \(\color{darkred}f\) is our

**flow variable**and \(\color{darkblue}g\) is

**exogenous inflow**. In our case:\[\color{darkblue}g_{i,j} = \begin{cases} 1 & \text{for the source node}\\ -1 & \text{for the sink node} \\ 0 & \text{for all other nodes}\end{cases}\]

---- 75 VARIABLEcost.L = 16.000objective---- 75 VARIABLEf.Lflowsrc. .i1 .j1 1.000, i1 .j1.i2 .j4 1.000, i2 .j4.i3 .j5 1.000, i3 .j5.snk. 1.000 ---- 79 SETassignassignments recovered from flowsj1 j4 j5 i1 YES i2 YES i3 YES

- Cplex default LP solver (dual simplex)
- Cplex network solver
- Sparse version of Dijkstra's algorithm from scipy.sparse.csgraph. This requires some programming to get the data organized in a form that is accepted by this function.

MIP Model | Network LP Model | Dijkstra | ||
---|---|---|---|---|

solver | MIP | LP default | Network | Sparse |

a/b length | 50/500 | 50/500 | 50/500 | 50/500 |

nodes/arcs | 22,552/4,995,276 | 22,552/4,995,276 | 22,552/4,995,276 | |

rows/columns/nz | 650/25,025/100,147 | 22,553/4,995,277/14,985,378 | 22,553/4,995,277/14,985,378 | |

objective | 6466.673 | 6466.673 | 6466.673 | 6466.673 |

time | 59 | 721 | total:38 presolve:32 extract network:2 solve network:3 | 0.1 |

b&b nodes | 5,023 | |||

iterations | 60,199 | 11,535 | 100,446 |

This looks better. My setup is a bit slow using unoptimized Python code, but the raw solve is very fast. When trying the \(m=100,n=1000\) data set, I got:

MIP Model | Dijkstra | |
---|---|---|

a/b length | 100/1000 | 100/1000 |

variables/constraints | 1,300/100,100 | |

nodes/arcs | 90,102/40,230,551 | |

objective | 14,371.455 | 14,371.455 |

time | 1,853 | 0.9 |

Note that the shortest path version would result in an LP with 90k rows and 40 million columns. That is very large. So here a specialized shortest path algorithm is far superior to more general tools.

#### Negative lengths

- When using Dijkstra with negative lengths, I see the message:
<ipython-input-121-4c8fd24b72ef>:67: 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 message seems to indicate that negative lengths are only a problem if there are negative cycles. As there are no cycles in our problem, does this imply we can use Dijkstra? I suspect this is not the case.
- Q: Does Dijkstra require non-negative lengths or only non-negative cycles? I think the answer is: non-negative lengths.
- Even with Johnson we may see:
`NegativeCycleError: Negative cycle detected on node 0`

- This can't be true, as there are no cycles in our network. It is a DAG (Direct Acyclic Graph). I think that the function complains as soon as a path has a negative length. That is slightly different.
- This is a bit more delicate than I expected. Not in the least because the error and warning messages are not as precise as they could be. Formulating good messages requires some thought, and I think that often not much thought is going into this.

#### References

- An assignment problem with a wrinkle, https://yetanothermathprogrammingconsultant.blogspot.com/2021/01/an-assignment-problem-with-wrinkle.html
- Dijkstra's algorithm, https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
- A monotonic assignment problem, https://orinanobworld.blogspot.com/2021/01/a-monotonic-assignment-problem.html has a Java implementation of this problem.

#### Appendix 1: GAMS code to solve network problem

$ontext $offtext * macros for sizing the problem$set m 10 $set n 100 $set irange i1*i%m% $set jrange j1*j%n% setsdummy 'ordering just for better
looking output' /src,%irange%,%jrange%,snk/i0 'includes
src/snk' /%irange%,src,snk/j0 'includes
blank for use with src/snk' /%jrange%,''/i(i0) 'core network' /%irange%/j(j0) 'core network' /%jrange%/ifirst(i) /i1/ ilast(i) /i%m%/ ; alias(i,ii),(j,jj),(i0,ii0),(j0,jj0);* random dataparametera(i) b(j) ; a(i) = uniform(0,100); b(j) = uniform(0,100); display a,b;scalarslena 'length of a' /%m%/lenb 'length of b' /%n%/; *--------------------------------------------------------------* network topology*--------------------------------------------------------------scalarsnumnodes numarcs ; set n(i0,j0) 'nodes (correspond to
assignments)';n(i,j) = ord(j)>=ord(i) and ord(j)<=lenb-lena+ord(i);n('src','') = yes;n('snk','') = yes;numnodes = card(n);display$(numnodes<1000)
n;display numnodes;alias(n,nn);set e(i0,j0,ii0,jj0) 'arcs';e(n(i,j),nn(i+1,jj)) = ord(jj)>ord(j);e('src','',n(ifirst,j)) = yes;e(n(ilast,j),'snk','') = yes;option e:0:0:7;numarcs = card(e);display$(numarcs<10000)
e;display numarcs;parameter
c(i0,j0,ii0,jj0) 'cost
of arcs';c(e(n,i,j)) = a(i)*b(j); option c:3:0:7;display$(numarcs<1000)
c;*--------------------------------------------------------------* shortest path model*--------------------------------------------------------------parameter
exinflow(i0,j0) 'exogenous
inflow' /src.'' 1 snk.'' -1 /; binary variable f(i0,j0,i0,j0) 'flow';variable cost 'objective';equationsobj 'objective'flowbal(i0,j0) 'flow balance at node n'; obj.. cost =e= sum(e, f(e)*c(e));flowbal(n).. sum(e(nn,n),f(e)) + exinflow(n) =e= sum(e(n,nn),f(e));model m /all/;option
limrow=0,limcol=0,solprint=off;m.optfile=1; solve m minimizing
cost using rmip;option f:3:0:7;display cost.l,f.l;set assign(i0,j0) 'assignments recovered from
flows';assign(n) = sum(e(nn,n),c(e)*f.l(e));option
assign:0:0:7;display assign;$onecho > cplex.opt * save some memory/timenames no* select network solverlpmethod 3$offecho |

#### Appendix B: Python code

import sys import pickle import scipy.sparse import scipy.sparse.csgraph import time print(sys.version) # # load data # data = pickle.load( open( r"d:\projects\blog\R\save.p", "rb" ) ) a = data['a'] b = data['b'] m = len(a) n = len(b) print("m={},n={}".format(m,n)) # # core network # nodes = {(i,j):((n-m+1)*i+j-i) for i in range(m) for j in range(i,n-m+i+1)} arcs = {(i,j,i+1,k):a[i+1]*b[k] for i in range(m-1) for j in range(i,n-m+i+1) for k in range(j+1,n-m+i+2)} # value for nodes is the node sequence number # value for arcs is the cost coefficient # # add source and sink node # src = ("src","src") snk = ("snk","snk") nodes[src]=len(nodes) nodes[snk]=len(nodes) # add arcs related to src,snk nodes for j in range(0,n-m+1): arcs[('src','src',0,j)] = a[0]*b[j] for j in range(m-1,n): arcs[(m-1,j,'snk','snk')] = 0 numnodes,numarcs = len(nodes),len(arcs) print("numnodes={},numarcs={}".format(numnodes,numarcs)) # # setup sparse matrix # row = [nodes[(i,j)] for (i,j,k,l) in arcs ] col = [nodes[(k,l)] for (i,j,k,l) in arcs ] weight = list(arcs.values()) spmat = scipy.sparse.csr_matrix((weight,(row,col)),shape=(numnodes,numnodes)) # # solve # source = [nodes[src]] start = time.time() sol = scipy.sparse.csgraph.dijkstra(spmat,indices=source) t = time.time()-start print("shortest path length: {}".format(sol[0,nodes[snk]])) print("time: {0:.1f} seconds".format(t))

3.9.1 (tags/v3.9.1:1e5d33e, Dec 7 2020, 17:08:21) [MSC v.1927 64 bit (AMD64)] m=100,n=1000 numnodes=90102,numarcs=40230551 shortest path length: 14371.455440918733 time: 0.9 seconds

I can reduce the time a bit by adding a special variant of ShortestPath which takes advantage of the fact that we have a DAG. In this case, the user would need to tell the solver that the graph is acyclic. It will throw an error if it is not.

ReplyDeleteBefore --

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

After --

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

That performance is mind-boggling. Not sure how to calculate the number of possible paths. Should be a little bit less than 1000^100.

DeleteWell, a lot less.

DeleteAre you asking how to enumerate all possible paths?

DeleteOr, are you asking how the shortest path algorithm is so fast - which obviously does not do full enumeration?

The standard approach uses Dijkstra's algorithm which is O(m + n log n). The specialized algorithm for a DAG uses topological sort and then a labeling pass similar to Dijkstra and is O(m + n).

Don't worry. I was just thinking out loud.

Delete