# Efficient solution for a discrete assignment problem with pairwise costs

Let's consider a simple graph with N vertices, and a corresponding set of N items. The goal of the problem is to assign every item to a vertex on the graph so that sum of a per-edge (that is item-pairwise) cost function over vertices is minimized.

This problem clearly has a discrete search space of N! candidates. Exhaustive search of the optimal solution becomes practically impossible with very low values of N, and I'd expect that there would be methods to put more efficient algorithms at work on this problem.

A brute-force toy attempt on the problem is presented below. Here g is the graph, i are the item values, and w is an extremely simple pairwise cost function:

With[{
g = GridGraph[{3, 3}],
i = {1, 4, 4, 9, 9, 16, 16, 32, 64},
w = Apply[Abs@*Subtract]},

First@TakeSmallestBy[
SetProperty[g,
VertexLabels -> MapThread[Rule, {Range@VertexCount@g, #}]] & /@
Permutations@i,
Function[g,
Total[w[PropertyValue[{g, #}, VertexLabels] & /@ #] & /@
EdgeList@g]], 1]] Please note that i and w are just examples; i might consist of, say, images, and w might be an earth mover's distance function, which makes neat assignment seen above impossible.

My more clever attempts this far have been based on an assumption this problem could be solved with integer linear programming. Sadly every attempt I've made to rephrase the problem statement in a suitable way for LinearProgramming has been either incomplete (resulting inconsistent edges and vertices), or ended up with an amount of constraints growing so big it's just moving the complexity to a new place.

Just to clarify: I'm not looking for methods to extract the last drop of exhaustive-search performance. Instead, I'm looking for algorithmic improvements in cases where search space consists of easily $10^{50}$ permutations, or more.

• The solution may depend on the shape of your graph. Is it consistently the same general shape? – rcollyer Jul 13 '15 at 15:34
• @rcollyer With high likelihood a grid graph - but support for more complicated cases, such as general planar graphs would also be good. – kirma Jul 13 '15 at 15:37
• I don't understand how the cost is computed. Are there weights involved? How do edges affect cost? – Daniel Lichtblau Jul 13 '15 at 16:00
• @DanielLichtblau Simple graph - no weights. To put it simply: take every edge on the graph, and apply a cost (or maybe more intuitively, distance) function to pair of items assigned on its' vertices. Then sum all these individual pairwise costs up. That's all. – kirma Jul 13 '15 at 16:04
• I would forgo the use of Graph and instead work with the AdjacencyMatrix and use LinearProgramming. – David G. Stork Jul 13 '15 at 16:52

Here is a linear programming surrogate that seems to do tolerably well. I'll show the code for your example but with a different cost function. I use 0-1 variables v[j,k] to indicate vertex j gets item k.

edges = {{1, 2}, {2, 3}, {1, 4}, {2, 5}, {4, 5}, {3, 6}, {5, 6}, {4,
7}, {7, 8}, {5, 8}, {6, 9}, {8, 9}};

klpairs =
With[{pairs = Subsets[Range, {2}]},
Riffle[pairs, Map[Reverse, pairs]]];
clist = Tuples[{edges, klpairs}];
vterms = Apply[v, Tuples[Range, {2}], {1}];

c1 = Map[0 <= # <= 1 &, vterms];
c2 = Table[Total[Table[v[j, k], {k, 9}]] == 1, {j, 9}];
c3 = Table[Total[Table[v[j, k], {j, 9}]] == 1, {k, 9}];
constraints = Flatten[Join[c1, c2, c3]];
cost = Total[
Map[(v[#[[1, 1]], #[[2, 1]]] + v[#[[1, 2]], #[[2, 2]]])*
f[#[]] &, clist]] -
2 Sum[If[l == k, 0, v[i, k]*f[{k, l}] + v[i, l]*f[{l, k}]], {i,
9}, {k, 9}, {l, 9}];
cost2 = Total[
Map[(v[#[[1, 1]], #[[2, 1]]]*v[#[[1, 2]], #[[2, 2]]])*f[#[]] &,
clist]];


In the above cost2 is the actual (nonlinear) thing we'd like to minimize, and cost is a linear surrogate that seems to do well at giving a low value for cost2.

Timing[{min, vals} =
NMinimize[{cost, constraints}, vterms];]
min

(* Out= {0.005997, Null}

Out= -53884. *)


How well did we do?

v2 = Chop[
Select[vals,
Head[#[]] == v && Chop[#[] /. #, 10^(-7)] == 0 &],
10^(-7)];
vterms2 = vterms /. v2 /. 0 :> Sequence[];
perm = Sort[vterms2][[All, 2]]
Total[Map[f[perm[[#]]] &, edges]]

(* Out= {6, 3, 9, 5, 1, 2, 7, 4, 8}

Out= 4494 *)


Double check with cost2.

cost2 /. Thread[vterms2 -> 1] /. v[__] :> 0

(* Out= 4494 *)


Random sampling does slightly better in this case.

perms = Table[RandomSample[Range], {10^5}];
fvals = Table[
Total[Map[f[perms[[j]][[#]]] &, edges]], {j, Length[perms]}];
Min[fvals]
Map[perms[[#]] &, Position[fvals, Min[fvals]][[All, 1]]]

(* Out= 4493

Out= {{8, 4, 6, 2, 1, 5, 9, 3, 7}, {8, 2, 9, 4, 1, 3, 6, 5, 7}} *)


But this tends to be slow.

I guess a good question is whether there might be better linear surrogates. Maybe someone else will find improvements for this.

I will remark that one can set this up as an ILP. Then we get a guaranteed minimizer. But some testing indicated that that can run a very long time without finishing. I've not tried coding an explicit branch and prune loop but that might also just hang.

Another option involves coding as a nonlinear permutation problem. This means no guarantee of optimality, but possibly generally decent results in reasonable time.

Yet another way is to use nonlinear continuous methods with penalties on nonintegrality. I show code for this in the section " Independent vertices in a graph" here, also with notebook in our library.

• I managed to write an possibly hugely unoptimal ILP implementation for this. When I feed that to LinearProgramming, getting results takes 1-2 orders of magnitude longer (even in the 9-vertex case) than in the case of the rather brainless exhaustive search. Thus, it really seems intelligent approximative methods would be better. I did try some pretty dumb methods on this regard couple months ago; advice on intelligent methods would be welcome! – kirma Jul 14 '15 at 18:04

I did some experimentation on Metropolis-Hastings algorithm for stochastic minimization of the cost function:

ClearAll@mhGraphPairwiseMinimize;

(* minimize sum of per-edge costs (computed as pairwise vertex item
distances) by assigning items to vertices in a graph, using a
Metropolis-Hastings algorithm. higher alpha makes random walk penalize
non-improving steps more. *)
mhGraphPairwiseMinimize[g_Graph, items_List, distFunc_Function,
alpha_, iter_Integer] :=
Module[{d, vertexToPosition, permutationCost, prependCost,
symmetricRandomPermute, proposedCandidate, acceptanceProbability,
newCandidate, initialCandidate, minimizationStep},
vertexToPosition[v_] := First@FirstPosition[VertexList@g, v];

(* cost function, constructed from graph, uses distance matrix d *)
permutationCost[perm_List] :=
Evaluate[
Total[Quiet@
d[[perm[[vertexToPosition@#1]],
perm[[vertexToPosition@#2]]]] & @@@ EdgeList@g]];

prependCost[perm_List] := {permutationCost@perm, perm};

(* just exchange two item-vertex mappings with each other,
with flat probability. this is both symmetric and transitively ergodic. *)
symmetricRandomPermute[{_, perm_List}] :=
Permute[perm, Cycles[{RandomSample[perm, 2]}]];

proposedCandidate[candidate_List] :=
prependCost@symmetricRandomPermute[candidate];

(* this is the acceptance probability for symmetric permutation distributions *)
acceptanceProbability[candnew_List, candold_List] :=
Min[1, Exp[-alpha (First@candnew - First@candold)/First@candold]];

newCandidate[candidate_List] :=
With[{newcand = proposedCandidate@candidate},
RandomChoice[{#, 1 - #} &@
acceptanceProbability[newcand, candidate] -> {newcand,
candidate}]];

(* distance matrix *)
d = Table[distFunc[a, b], {a, items}, {b, items}];

(* a random permutation *)
initialCandidate = prependCost@RandomSample@Range@VertexCount@g;

(* calculate next (possibly unchanged) random walk value,
and update minimum *)
minimizationStep[{lastmin_List, candidate_List}] :=
With[{newcand = newCandidate@candidate}, {First@
TakeSmallestBy[{lastmin, newcand}, First, 1], newcand}];

(* return the minimum found weight and the corresponding vertex -> item list *)
{#1, 