Efficient solution for a discrete assignment problem with pairwise costs

Question

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.

Answer 1

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[9], {2}]}, 
   Riffle[pairs, Map[Reverse, pairs]]];
clist = Tuples[{edges, klpairs}];
vterms = Apply[v, Tuples[Range[9], {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[#[[2]]] &, 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[#[[2]]] &, 
    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[1781]= {0.005997, Null}

Out[1782]= -53884. *)

How well did we do?

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

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

Out[1786]= 4494 *)

Double check with cost2.

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

(* Out[1787]= 4494 *)

Random sampling does slightly better in this case.

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

(* Out[1790]= 4493

Out[1791]= {{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.

Answer 2

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, 
      MapThread[Rule, {VertexList@g, items[[#2]]}]} & @@ 
    First@Nest[minimizationStep, {initialCandidate, initialCandidate},
       iter]];

Module[{g, weight, sol},
 g = GridGraph[{3, 3}];
 {weight, sol} = mhGraphPairwiseMinimize[
   g, {1, 4, 4, 9, 9, 16, 16, 32, 64}, Abs[#1 - #2] &, 500, 100];
 {weight, Graph[g, VertexLabels -> sol]}]